diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md
index 37f3843..780ec2a 100644
--- a/.claude/CLAUDE.md
+++ b/.claude/CLAUDE.md
@@ -88,10 +88,12 @@ tests/testthat/gradient_golden.tsv
                                on that platform only, like tests/cpp/golden/ --
                                these are derivatives of argmax-evaluated outputs,
                                so they inherit its sqrt-amplified class and
-                               disagree cross-platform by up to 1.3e-3 -- ten
-                               times the solved outputs, because a finite
-                               difference divides the solver floor by the step
-                               and `R_d_25`'s step is the smallest here
+                               disagree cross-platform by up to 2.34e-3 (was
+                               1.3e-3 before #87 added `profit`) -- twenty times
+                               the solved outputs, because a finite difference
+                               divides the solver floor by the step, `R_d_25`'s
+                               step is the smallest here, and `profit` is the
+                               smallest magnitude carrying that same floor
 tests/validate/                R scripts comparing against plant (needs R)
 CMakeLists.txt                 the no-R build: C++ and Python consumers, and the
                                thing that makes "does not need R" runnable
diff --git a/NAMESPACE b/NAMESPACE
index 3a6f740..c616d2d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,6 +3,7 @@
 S3method(print,leaf_batch)
 export(Leaf)
 export(RootNetwork)
+export(gradient_output_names)
 export(gradient_par_names)
 export(leaf_batch)
 export(leaf_control)
diff --git a/NEWS.md b/NEWS.md
index 919a946..0ed7814 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -162,6 +162,66 @@ of the wrong shape — see the note in `update_temperature_dependent_params()`.
 `rd_to_vcmax_ratio` is still not a `leaf_traits()` member, so it cannot yet be
 fitted or differentiated; it is set as a field.
 
+## The gradient differentiates `profit`, which is what a demographic caller bills
+
+`leaf_gradient()` and `leaf_gradient_batch()` return a fifth column. The four
+that were there — `A`, `gc`, `psi_stem`, `collar` — are what a gas-exchange
+calibration observes, and they were chosen for the customer this feature was
+built for. They are **disjoint** from what `plant` reads off a solved leaf: its
+carbon budget is `leaf.profit_` (not `assim_colimited_`) and its water budget is
+`leaf.soil_consumption_`. So no trait gradient this package produced reached a
+demographic model at all, at the optimum or anywhere else
+([#87](https://github.com/traitecoevo/phylloptim/issues/87)).
+
+⚠️ **`profit` is the one output the envelope theorem reaches, and the only place
+this package uses it.** At an interior optimum `dprofit/dpsi = 0`, so the
+indirect term `(dprofit/dpsi)(dpsi*/dtheta)` vanishes identically and
+`dprofit/dtheta` is the direct partial at fixed ψ — no `dY/dpsi`, no `−M/H`. It
+is set from that term rather than computed through the composite, exactly as
+`collar` is set from `dpsi*/dtheta` rather than differenced.
+
+**That is a numerical decision, not a tidiness one, and the measurement is the
+reason.** The dropped term is *noise*, not an `h²` truncation: `profit` is the
+maximum, so it is flat, and a central difference of it divides the solve's ~1e-09
+floor by a ~1e-06 step. Over the golden grid's 136 interior rows —
+
+| | median | max |
+|---|---|---|
+| `\|dprofit/dpsi\|`, exact (forward AD) | 4.8e-15 | 5.4e-10 |
+| `\|dprofit/dpsi\|`, central difference | 7.8e-10 | 2.1e-04 |
+| relative move in `dprofit/dtheta` if kept | 2.7e-10 | 8.0e-05 |
+
+— eleven orders between the two instruments at the median, and the worst row sits
+in the band this repo calls a real difference rather than rounding. The identity
+is applied only where `status == "interior"`; at a pinned optimum `psi*` is a
+trait-dependent bound, `dprofit/dpsi` is not zero, and `profit` takes the same
+finite-difference fallback as the other four.
+
+The four existing columns are **bit-identical** — this is additive, like
+appending to `gradient_par_names()`. `tests/testthat/gradient_golden.tsv` gains a
+column and no existing cell moved, checked against master rather than against the
+branch point.
+
+⚠️ **The shut-down row's profit column is asserted against a closed form, not
+only recorded.** It is the one regime where `profit_` is written by a branch that
+leaves the other outputs alone (hazard 8), so a hex with nothing saying what it
+ought to be would pin a number rather than a fact. There `E = 0`, so `A = -R_d`
+exactly and the hydraulic cost does not depend on `R_d_25`: `dprofit/dR_d_25` is
+**−1**. And the shut-down collar is pinned at `psi_crit`, so
+`dcollar/dpsi_crit` is **1** — which is why `psi_crit` alone carries a non-zero
+profit gradient there.
+
+**`gradient_output_names()` is exported**, and R now *reads* the list rather than
+keeping a second copy. `gradient_par_names()` has to be duplicated-and-compared
+because R builds `theta` before any C++ call; the outputs have no such
+constraint, so adding one is a single edit.
+
+⚠️ **`uptake` was considered and is not here.** Every output must be a field R
+*copies* out of `operating_point_values()`; `uptake` is one R *computes*, by
+summing over the finite soil layers, so adding it means reproducing that
+summation and its order on the C++ side too. That is a separate decision from
+this one.
+
 ## An out-of-domain transport lookup says which spline, and which caller
 
 The stem curve is the only interpolator here built with extrapolation disabled, so
diff --git a/R/RcppExports.R b/R/RcppExports.R
index bffcb1a..db96ca4 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -774,6 +774,26 @@ gradient_par_names <- function() {
     .Call('_phylloptim_gradient_par_names', PACKAGE = 'phylloptim')
 }
 
+#' The differentiated outputs, in the order C++ indexes them
+#'
+#' `A`, `gc`, `psi_stem`, `collar` and `profit`. Unlike [gradient_par_names()],
+#' which R holds its own copy of and a test compares, R **reads** this one — so
+#' the list exists once and adding an output is one edit rather than two that
+#' can disagree.
+#'
+#' The first four are what a gas-exchange calibration observes; `profit` is what
+#' a demographic consumer bills, and it is the one output the envelope theorem
+#' reaches. See [leaf_gradient()].
+#'
+#' @return A character vector of five names.
+#' @seealso [leaf_gradient()], [leaf_gradient_batch()]
+#' @examples
+#' gradient_output_names()
+#' @export
+gradient_output_names <- function() {
+    .Call('_phylloptim_gradient_output_names', PACKAGE = 'phylloptim')
+}
+
 gradient_batch_prepare <- function(root_network, psi_soil, soil_depth, PPFD, atm_vpd, ca, leaf_temp, atm_o2_kpa, atm_kpa) {
     .Call('_phylloptim_gradient_batch_prepare', PACKAGE = 'phylloptim', root_network, psi_soil, soil_depth, PPFD, atm_vpd, ca, leaf_temp, atm_o2_kpa, atm_kpa)
 }
diff --git a/R/gradient-batch.R b/R/gradient-batch.R
index c5a9c56..5f641cc 100644
--- a/R/gradient-batch.R
+++ b/R/gradient-batch.R
@@ -188,7 +188,7 @@ print.leaf_batch <- function(x, ...) {
 ##' observation.
 ##'
 ##' @section What it computes, and where the maths is written down:
-##' The same four derivatives at the same solved operating point, by the same two
+##' The same five derivatives at the same solved operating point, by the same two
 ##' routes, with the same active-set test deciding between them. **Read
 ##' [leaf_gradient()]'s documentation for all of it** — the implicit function
 ##' theorem, why the second term is not a correction, why the premise is tested
@@ -239,9 +239,9 @@ print.leaf_batch <- function(x, ...) {
 ##'
 ##' @return A list with
 ##'   \describe{
-##'     \item{`gradient`}{an `n` × `length(pars)` × 4 array, the last dimension
-##'       being `A`, `gc`, `psi_stem` and `collar`}
-##'     \item{`value`}{an `n` × 4 matrix of the solved outputs}
+##'     \item{`gradient`}{an `n` × `length(pars)` × 5 array, the last dimension
+##'       being `A`, `gc`, `psi_stem`, `collar` and `profit`}
+##'     \item{`value`}{an `n` × 5 matrix of the solved outputs}
 ##'     \item{`method`}{`"ift"` or `"fd"` per observation}
 ##'     \item{`status`}{`"interior"`, `"pinned"`, `"no-gradient"` or `"error"`}
 ##'     \item{`H`, `stationarity`}{the curvature and the implied Newton step
@@ -314,8 +314,8 @@ leaf_gradient_batch <- function(batch,
                             match(pars, par_names) - 1L, step,
                             stationarity_tol, method, fast_stem_curve)
 
-  dimnames(res$gradient) <- list(NULL, pars, .gradient_output_names)
-  dimnames(res$value) <- list(NULL, .gradient_output_names)
+  dimnames(res$gradient) <- list(NULL, pars, .gradient_output_names())
+  dimnames(res$value) <- list(NULL, .gradient_output_names())
   res
 }
 
diff --git a/R/gradient.R b/R/gradient.R
index a3e6580..6f2a5c8 100644
--- a/R/gradient.R
+++ b/R/gradient.R
@@ -96,8 +96,9 @@ set_traits <- function(x, traits) {
 
 ##' Trait gradients of a solved operating point
 ##'
-##' The derivatives of the solved outputs with respect to the traits: `dA/dtheta`,
-##' `dgc/dtheta`, `dpsi_stem/dtheta` and `dcollar/dtheta`, at one operating point.
+##' The derivatives of the solved outputs with respect to the traits:
+##' `dA/dtheta`, `dgc/dtheta`, `dpsi_stem/dtheta`, `dcollar/dtheta` and
+##' `dprofit/dtheta`, at one operating point.
 ##'
 ##' @section The maths, and why it is not just a finite difference:
 ##' The reported outputs are evaluated at the profit-maximising collar potential
@@ -120,6 +121,44 @@ set_traits <- function(x, traits) {
 ##' The second term is not a correction. For `cost_scale_TF24`, `beta2`, `stem_b`
 ##' and `stem_c` it is 100% of the answer, and for `vcmax_25` 52%.
 ##'
+##' @section profit, which is the one output the envelope theorem reaches:
+##' The first four outputs are what a gas-exchange calibration observes. `profit`
+##' is here because it is what a **demographic** consumer bills: plant's carbon
+##' budget reads `leaf.profit_`, not `assim_colimited_`, so before it was added
+##' the two sets were disjoint and no gradient from this package reached a
+##' demographic model at all.
+##'
+##' It is also the cheapest of the five, because it is the objective. At an
+##' interior optimum `dprofit/dpsi = 0`, so the indirect term above vanishes
+##' identically and
+##'
+##' \deqn{dprofit/d\theta = \partial profit/\partial\theta|_{\psi^*}}
+##'
+##' -- the direct partial alone, with no `dY/dpsi` and no `-M/H`. That is the
+##' envelope theorem, and it is the only place this package uses it.
+##'
+##' ⚠️ **The term this drops is noise, not an `h^2` truncation, and that is why
+##' dropping it matters.** `profit` is the maximum, so it is flat, and a central
+##' difference of it at `psi*` divides the solve's ~1e-09 floor by a ~1e-06 step.
+##' Over the golden grid's 136 interior rows:
+##'
+##' | | median | max |
+##' |---|---|---|
+##' | `\|dprofit/dpsi\|`, exact (forward AD) | 4.8e-15 | 5.4e-10 |
+##' | `\|dprofit/dpsi\|`, central difference | 7.8e-10 | 2.1e-04 |
+##' | relative move in `dprofit/dtheta` if kept | 2.7e-10 | 8.0e-05 |
+##'
+##' Eleven orders between the two instruments at the median, and the worst row
+##' sits in the band this package calls a real difference. So `dprofit/dtheta` is
+##' set from the direct term rather than computed through the composite, exactly
+##' as `collar` is set from `dpsi*/dtheta` rather than differenced.
+##'
+##' ⚠️ **The envelope does not survive the active set.** At a pinned optimum
+##' `psi*` is a trait-dependent bound and `dprofit/dpsi` is not zero there, so the
+##' indirect term is real. This applies the identity only where `status` is
+##' `"interior"`; everywhere else `profit` goes through the same finite-difference
+##' fallback as the other four.
+##'
 ##' @section The active set, which is the reason this function is careful:
 ##' Stationarity is the premise of the whole derivation, and it fails when the
 ##' optimum is pinned to an end of the feasible collar interval. There `psi*` is a
@@ -288,7 +327,8 @@ set_traits <- function(x, traits) {
 ##' @return A list with
 ##'   \describe{
 ##'     \item{`gradient`}{a matrix, one row per trait in `pars`, with columns
-##'       `A` (umol C m^-2 s^-1 per trait unit), `gc`, `psi_stem` and `collar`}
+##'       `A` (umol C m^-2 s^-1 per trait unit), `gc`, `psi_stem`, `collar` and
+##'       `profit` (also umol C m^-2 s^-1 per trait unit)}
 ##'     \item{`value`}{the solved outputs the gradient is taken at}
 ##'     \item{`method`}{`"ift"` if the implicit-function composite was used,
 ##'       `"fd"` if the fallback was}
@@ -497,12 +537,40 @@ leaf_gradient <- function(psi_soil,
       status <- "pinned"
     } else {
       dY_dpsi <- (hi - lo) / (2 * h_psi)
+      # `dprofit_droot_collar_psi` is EXACT in psi -- forward AD plus the IFT at
+      # the ci root-find -- so for profit alone the package has something better
+      # than a difference of the same quantity, and it is already computed. The
+      # other four have no such route and must be differenced.
+      #
+      # ⚠️ NOTHING CONSUMES THIS TODAY, and it is here rather than deleted for #88.
+      # The only reader is `.gradient_ift()` with `envelope = FALSE`, i.e. a FORCED
+      # method = "ift" at a pinned point -- and the 288-point grid test records that
+      # forcing it there throws ("narrower than one step") at all 42 pinned rows
+      # before this value is reached. So it is unexercised, not load-bearing: do
+      # not read a green suite as evidence about it.
+      dY_dpsi[["profit"]] <- resid
     }
   }
 
+  # ⚠️ THE ENVELOPE THEOREM, and the only place this package uses it. At a
+  # STATIONARY point dprofit/dpsi is analytically zero, so profit's indirect term
+  # vanishes identically and dprofit/dtheta is the direct partial alone. The
+  # composite is told to ASSIGN that column rather than reach it by multiplying a
+  # near-zero dY/dpsi -- the same treatment `collar` gets, and for the same
+  # reason: an identity is stated, not arrived at.
+  #
+  # ⚠️ It is conditional on stationarity and not on `use_ift`. The identity comes
+  # from dprofit/dpsi == 0; at a pinned optimum psi* is a theta-dependent BOUND,
+  # dprofit/dpsi is not zero there, and the indirect term survives. Someone who
+  # forces method = "ift" at such a point already gets a confidently wrong number
+  # and should not get a differently wrong one for this column alone.
+  #
+  # The measured size of what this removes is in ?leaf_gradient.
+  envelope <- isTRUE(usable && stationarity <= stationarity_tol)
+
   grad <- if (use_ift) {
     .gradient_ift(l, reset, theta, pars, psi_star, H, dY_dpsi, step,
-                  fast_stem_curve)
+                  fast_stem_curve, envelope = envelope)
   } else {
     .gradient_fd(l, reset, theta, pars, step, fast_stem_curve)
   }
@@ -517,11 +585,34 @@ leaf_gradient <- function(psi_soil,
 
 # --- internals ---------------------------------------------------------------
 
-# The four differentiated outputs, in one place so the composite and the fallback
-# cannot disagree about what they are. `collar` is psi* itself, which makes
-# dcollar/dtheta equal to dpsi*/dtheta -- so the two routes compute the same
-# quantity by different means, and a test can compare them.
-.gradient_output_names <- c("A", "gc", "psi_stem", "collar")
+# The differentiated outputs, READ from C++ rather than restated here.
+#
+# ⚠️ This is deliberately NOT the pattern `.gradient_par_names` uses. That one is a
+# second copy kept honest by a test comparing it with `gradient_par_names()`, and
+# it has to be, because R builds `theta` from `leaf_traits()` and so needs the
+# order before any C++ call. The OUTPUT list has no such constraint, so it is a
+# single definition and there is no hazard comment to write: adding an output is
+# one edit in `gradient.hpp`, not two that can disagree.
+#
+# `collar` is psi* itself, which makes dcollar/dtheta equal to dpsi*/dtheta -- so
+# the two routes compute the same quantity by different means and a test can
+# compare them. `profit` is what plant CONSUMES (`leaf.profit_`, not
+# `assim_colimited_`) and is the one output the envelope theorem reaches; see
+# ?leaf_gradient.
+#
+# ⚠️ Cached at FIRST CALL, not at build time, for the reason
+# `.gradient_outputs_idx()` below records -- and here there is a second reason:
+# calling into the shared library while the package is still being sourced is not
+# something to rely on.
+.gradient_output_names <- local({
+  nms <- NULL
+  function() {
+    if (is.null(nms)) {
+      nms <<- gradient_output_names()
+    }
+    nms
+  }
+})
 
 # ⚠️ ONE call, not four field reads. Every `l$field` is an R6 ACTIVE BINDING -- a
 # closure call wrapping a `.Call` -- and this function runs once per perturbation, so
@@ -541,7 +632,7 @@ leaf_gradient <- function(psi_soil,
   idx <- NULL
   function() {
     if (is.null(idx)) {
-      idx <<- match(.gradient_output_names, .operating_point_names)
+      idx <<- match(.gradient_output_names(), .operating_point_names)
     }
     idx
   }
@@ -549,7 +640,7 @@ leaf_gradient <- function(psi_soil,
 
 .gradient_outputs <- function(l) {
   v <- l$operating_point_values()[.gradient_outputs_idx()]
-  names(v) <- .gradient_output_names
+  names(v) <- .gradient_output_names()
   v
 }
 
@@ -815,7 +906,7 @@ leaf_gradient <- function(psi_soil,
 # neither of which re-solves the model: `dprofit` at the UNPERTURBED psi* gives
 # the mixed partial, and the outputs at that same psi* give the direct term.
 .gradient_ift <- function(l, reset, theta, pars, psi_star, H, dY_dpsi, step,
-                          fast_stem_curve = TRUE) {
+                          fast_stem_curve = TRUE, envelope = FALSE) {
   seat <- .gradient_reseat_base(reset, theta, pars, fast_stem_curve)
   out <- t(vapply(seq_along(pars), function(k) {
     p <- pars[[k]]
@@ -846,8 +937,16 @@ leaf_gradient <- function(psi_soil,
     # difference of two identical numbers.
     g <- direct + dY_dpsi * dpsi_dtheta
     g[["collar"]] <- dpsi_dtheta
-    g[.gradient_output_names]
-  }, numeric(length(.gradient_output_names))))
+    # The envelope theorem, ASSIGNED for the same reason `collar` is: profit's
+    # indirect term is identically zero at a stationary point, so stating that is
+    # better than multiplying a measured near-zero by dpsi*/dtheta. It is also
+    # immune to a non-finite dpsi*/dtheta, where `0 * x` would be NaN in this one
+    # column while the other four carried +-Inf.
+    if (envelope) {
+      g[["profit"]] <- direct[["profit"]]
+    }
+    g[.gradient_output_names()]
+  }, numeric(length(.gradient_output_names()))))
   reset(theta)
   rownames(out) <- pars
   out
@@ -869,8 +968,8 @@ leaf_gradient <- function(psi_soil,
       l$find_root_collar_psi()
       .gradient_outputs(l)
     }
-    ((side(1) - side(-1)) / (2 * h))[.gradient_output_names]
-  }, numeric(length(.gradient_output_names))))
+    ((side(1) - side(-1)) / (2 * h))[.gradient_output_names()]
+  }, numeric(length(.gradient_output_names()))))
   rownames(out) <- pars
   reset(theta)
   out
diff --git a/README.md b/README.md
index 33947ed..e98d9f8 100644
--- a/README.md
+++ b/README.md
@@ -324,10 +324,15 @@ traits, which is what a gradient-based optimiser or a Hamiltonian sampler wants:
 ```r
 g <- leaf_gradient(psi_soil = 2.0, PPFD = 900,
                    pars = c("vcmax_25", "stem_b", "cost_scale_TF24"))
-g$gradient   # rows: parameters.  columns: A, gc, psi_stem, collar
+g$gradient   # rows: parameters.  columns: A, gc, psi_stem, collar, profit
 g$method     # "ift" or "fd" -- see below
 ```
 
+The first four columns are what a gas-exchange calibration observes. `profit` is
+the objective, and it is there for a demographic consumer: `plant` bills carbon
+from the leaf's profit rather than its assimilation, so without that column no
+gradient from this package reached a demographic model at all.
+
 `pars` is not restricted to traits: `leaf_specific_conductance_max` and, on the
 single-potential path, `resistance` are differentiable too, because a calibration
 fits them and nothing in the derivation cares whether a parameter is a trait.
@@ -345,6 +350,11 @@ moving that optimum — and for `cost_scale_TF24`, `beta2`, `stem_b` and `stem_c
 the second route is **100%** of the answer. Differentiating the optimality
 condition rather than the solved output gets both terms exactly.
 
+`profit` is the exception, and it is the cheapest column for the reason it is the
+exception: it *is* the objective, so at an interior optimum the second route
+contributes nothing and its gradient is the direct partial alone. That is the
+envelope theorem, and the only place this package uses it.
+
 That derivation assumes the optimum is interior, and at the dry end it often is
 not: with the optimum pinned to the edge of the feasible range the formula returns
 a confidently wrong number, off by up to seven orders of magnitude. So the
diff --git a/inst/include/phylloptim/gradient.hpp b/inst/include/phylloptim/gradient.hpp
index 224d33c..8abbdc0 100644
--- a/inst/include/phylloptim/gradient.hpp
+++ b/inst/include/phylloptim/gradient.hpp
@@ -105,29 +105,44 @@ inline const std::vector<std::string>& par_names() {
   return names;
 }
 
-// --- the four differentiated outputs -----------------------------------------
+// --- the five differentiated outputs ------------------------------------------
 //
-// A, gc, psi_stem and collar, in that order, which is R's
+// A, gc, psi_stem, collar and profit, in that order, which is R's
 // `.gradient_output_names`. `collar` is psi* itself, which is what makes
 // dcollar/dtheta equal dpsi*/dtheta and lets the two routes below compute the
 // same quantity by different means.
-inline constexpr int n_outputs = 4;
+//
+// ⚠️ APPENDING IS SAFE AND REORDERING IS NOT, exactly as for `par_names` above.
+//
+// The first four are what a gas-exchange calibration OBSERVES. `profit` is here
+// because it is what plant CONSUMES: `leaf.profit_`, not `assim_colimited_`, is
+// the carbon that reaches its mass budget, so until #87 the two sets were
+// disjoint and no gradient this package produced reached a demographic model.
+inline constexpr int n_outputs = 5;
 inline constexpr int out_collar = 3;
+inline constexpr int out_profit = 4;
 
 inline const std::vector<std::string>& output_names() {
-  static const std::vector<std::string> names{"A", "gc", "psi_stem", "collar"};
+  static const std::vector<std::string> names{"A", "gc", "psi_stem", "collar",
+                                              "profit"};
   return names;
 }
 
 // Read straight off the members rather than through `operating_point_values()`,
-// which is what R has to use. Bit-identical: that reader copies these same four
-// fields into positions 3, 5, 0 and 1 of its twelve, and the three columns it
+// which is what R has to use. Bit-identical: that reader copies these same five
+// fields into positions 3, 5, 0, 1 and 6 of its twelve, and the three columns it
 // computes rather than copies (uptake, lambda, g1_eff) are not among them.
+//
+// ⚠️ That is why `profit` was cheap to add here and `uptake` would not be. Every
+// output in this list has to be a field R COPIES; `uptake` is one R sums over
+// the finite soil layers, so adding it means reproducing that summation -- and
+// its order -- on this side too.
 inline void outputs(const Leaf& l, double* y) {
   y[0] = l.assim_colimited_;
   y[1] = l.stom_cond_CO2_;
   y[2] = l.opt_psi_stem_;
   y[3] = l.opt_root_psi_;
+  y[4] = l.profit_;
 }
 
 // The outputs with the collar held at `psi` rather than optimised. False when
@@ -228,9 +243,10 @@ struct Settings {
 };
 
 struct Result {
-  // The four outputs the gradient is taken at.
+  // The five outputs the gradient is taken at.
   double value[n_outputs];
-  // npars * n_outputs, parameter-major: d(output j)/d(pars[k]) at [k * 4 + j].
+  // npars * n_outputs, parameter-major: d(output j)/d(pars[k]) at
+  // [k * n_outputs + j].
   std::vector<double> grad;
   Status status = Status::Error;
   bool used_ift = false;
@@ -337,7 +353,7 @@ inline bool takes_shortcut(int par, const Settings& s) {
 inline void gradient_ift(Leaf& l, const double* theta, const Drivers& d,
                          bool single, const int* pars, std::size_t npars,
                          double psi_star, double H, const double* dY_dpsi,
-                         const Settings& s, double* out) {
+                         const Settings& s, bool envelope, double* out) {
   double th[n_pars];
   double up[1 + n_outputs];
   double dn[1 + n_outputs];
@@ -378,6 +394,15 @@ inline void gradient_ift(Leaf& l, const double* theta, const Drivers& d,
     // dpsi*/dtheta. Set explicitly rather than left as the difference of two
     // identical numbers.
     out[k * n_outputs + out_collar] = dpsi_dtheta;
+    // The envelope theorem, ASSIGNED for the same reason `collar` is: profit's
+    // indirect term is identically zero at a stationary point, so stating that
+    // beats multiplying a measured near-zero by dpsi*/dtheta. It is also immune
+    // to a non-finite dpsi*/dtheta, where `0 * x` would be NaN in this one column
+    // while the other four carried +-Inf.
+    if (envelope) {
+      out[k * n_outputs + out_profit] =
+          (up[1 + out_profit] - dn[1 + out_profit]) / (2.0 * h);
+    }
   }
   apply(l, theta, d, single, -1, s.fast_stem_curve);
 }
@@ -501,13 +526,39 @@ inline void at(Leaf& l, const double* theta, const Drivers& d, bool single,
       for (int j = 0; j < n_outputs; ++j) {
         dY_dpsi[j] = (hi[j] - lo[j]) / (2.0 * h_psi);
       }
+      // `dprofit_droot_collar_psi` is EXACT in psi -- forward AD plus the IFT at
+      // the ci root-find -- so for profit alone there is something better than a
+      // difference of the same quantity, and it is already computed. The other
+      // four have no such route and must be differenced.
+      //
+      // ⚠️ NOTHING CONSUMES THIS TODAY, and it is here rather than deleted for
+      // #88. The only reader is `gradient_ift` with `envelope` false, i.e. a
+      // FORCED Method::Ift at a pinned point -- and the 288-point grid test
+      // records that forcing it there throws ("narrower than one step") at all
+      // 42 pinned rows before this value is reached. So it is unexercised, not
+      // load-bearing: do not read a green suite as evidence about it.
+      dY_dpsi[out_profit] = resid;
     }
   }
 
+  // ⚠️ THE ENVELOPE THEOREM, and the only place this package uses it. At a
+  // STATIONARY point dprofit/dpsi is analytically zero, so profit's indirect
+  // term vanishes identically and dprofit/dtheta is the direct partial alone.
+  // `gradient_ift` is told to ASSIGN that column rather than reach it by
+  // multiplying a near-zero dY/dpsi -- the same treatment `collar` gets, for the
+  // same reason: an identity is stated, not arrived at.
+  //
+  // ⚠️ Conditional on stationarity, not on `use_ift`. The identity comes from
+  // dprofit/dpsi == 0; at a pinned optimum psi* is a theta-dependent BOUND,
+  // dprofit/dpsi is not zero there, and the indirect term survives. Someone
+  // forcing Method::Ift there already gets a confidently wrong number and should
+  // not get a differently wrong one for this column alone.
+  const bool envelope = usable && out.stationarity <= s.stationarity_tol;
+
   out.used_ift = use_ift;
   if (use_ift) {
     gradient_ift(l, theta, d, single, pars, npars, psi_star, H, dY_dpsi, s,
-                 out.grad.data());
+                 envelope, out.grad.data());
   } else {
     gradient_fd(l, theta, d, single, pars, npars, s, out.grad.data());
   }
diff --git a/man/gradient_output_names.Rd b/man/gradient_output_names.Rd
new file mode 100644
index 0000000..beb4906
--- /dev/null
+++ b/man/gradient_output_names.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{gradient_output_names}
+\alias{gradient_output_names}
+\title{The differentiated outputs, in the order C++ indexes them}
+\usage{
+gradient_output_names()
+}
+\value{
+A character vector of five names.
+}
+\description{
+`A`, `gc`, `psi_stem`, `collar` and `profit`. Unlike [gradient_par_names()],
+which R holds its own copy of and a test compares, R **reads** this one — so
+the list exists once and adding an output is one edit rather than two that
+can disagree.
+}
+\details{
+The first four are what a gas-exchange calibration observes; `profit` is what
+a demographic consumer bills, and it is the one output the envelope theorem
+reaches. See [leaf_gradient()].
+}
+\examples{
+gradient_output_names()
+}
+\seealso{
+[leaf_gradient()], [leaf_gradient_batch()]
+}
diff --git a/man/gradient_par_names.Rd b/man/gradient_par_names.Rd
index 058d15d..7a1d697 100644
--- a/man/gradient_par_names.Rd
+++ b/man/gradient_par_names.Rd
@@ -7,10 +7,10 @@
 gradient_par_names()
 }
 \value{
-A character vector of fifteen names.
+A character vector of sixteen names.
 }
 \description{
-The thirteen [leaf_traits()] in `set_traits()`'s argument order, then
+The fourteen [leaf_traits()] in `set_traits()`'s argument order, then
 `leaf_specific_conductance_max` and `resistance`. Exported so that R's own
 copy of this order can be compared against it in a test: R passes integer
 positions into this enumeration, so appending to it is safe and reordering it
diff --git a/man/leaf_gradient.Rd b/man/leaf_gradient.Rd
index 21bd86c..65f8542 100644
--- a/man/leaf_gradient.Rd
+++ b/man/leaf_gradient.Rd
@@ -116,7 +116,8 @@ the two routes to agree. `stem_c` has no such identity and always rebuilds.}
 A list with
   \describe{
     \item{`gradient`}{a matrix, one row per trait in `pars`, with columns
-      `A` (umol C m^-2 s^-1 per trait unit), `gc`, `psi_stem` and `collar`}
+      `A` (umol C m^-2 s^-1 per trait unit), `gc`, `psi_stem`, `collar` and
+      `profit` (also umol C m^-2 s^-1 per trait unit)}
     \item{`value`}{the solved outputs the gradient is taken at}
     \item{`method`}{`"ift"` if the implicit-function composite was used,
       `"fd"` if the fallback was}
@@ -127,8 +128,9 @@ A list with
   }
 }
 \description{
-The derivatives of the solved outputs with respect to the traits: `dA/dtheta`,
-`dgc/dtheta`, `dpsi_stem/dtheta` and `dcollar/dtheta`, at one operating point.
+The derivatives of the solved outputs with respect to the traits:
+`dA/dtheta`, `dgc/dtheta`, `dpsi_stem/dtheta`, `dcollar/dtheta` and
+`dprofit/dtheta`, at one operating point.
 }
 \section{The maths, and why it is not just a finite difference}{
 
@@ -153,6 +155,46 @@ The second term is not a correction. For `cost_scale_TF24`, `beta2`, `stem_b`
 and `stem_c` it is 100% of the answer, and for `vcmax_25` 52%.
 }
 
+\section{profit, which is the one output the envelope theorem reaches}{
+
+The first four outputs are what a gas-exchange calibration observes. `profit`
+is here because it is what a **demographic** consumer bills: plant's carbon
+budget reads `leaf.profit_`, not `assim_colimited_`, so before it was added
+the two sets were disjoint and no gradient from this package reached a
+demographic model at all.
+
+It is also the cheapest of the five, because it is the objective. At an
+interior optimum `dprofit/dpsi = 0`, so the indirect term above vanishes
+identically and
+
+\deqn{dprofit/d\theta = \partial profit/\partial\theta|_{\psi^*}}
+
+-- the direct partial alone, with no `dY/dpsi` and no `-M/H`. That is the
+envelope theorem, and it is the only place this package uses it.
+
+⚠️ **The term this drops is noise, not an `h^2` truncation, and that is why
+dropping it matters.** `profit` is the maximum, so it is flat, and a central
+difference of it at `psi*` divides the solve's ~1e-09 floor by a ~1e-06 step.
+Over the golden grid's 136 interior rows:
+
+| | median | max |
+|---|---|---|
+| `\|dprofit/dpsi\|`, exact (forward AD) | 4.8e-15 | 5.4e-10 |
+| `\|dprofit/dpsi\|`, central difference | 7.8e-10 | 2.1e-04 |
+| relative move in `dprofit/dtheta` if kept | 2.7e-10 | 8.0e-05 |
+
+Eleven orders between the two instruments at the median, and the worst row
+sits in the band this package calls a real difference. So `dprofit/dtheta` is
+set from the direct term rather than computed through the composite, exactly
+as `collar` is set from `dpsi*/dtheta` rather than differenced.
+
+⚠️ **The envelope does not survive the active set.** At a pinned optimum
+`psi*` is a trait-dependent bound and `dprofit/dpsi` is not zero there, so the
+indirect term is real. This applies the identity only where `status` is
+`"interior"`; everywhere else `profit` goes through the same finite-difference
+fallback as the other four.
+}
+
 \section{The active set, which is the reason this function is careful}{
 
 Stationarity is the premise of the whole derivation, and it fails when the
diff --git a/man/leaf_gradient_batch.Rd b/man/leaf_gradient_batch.Rd
index 156dc02..5b662a6 100644
--- a/man/leaf_gradient_batch.Rd
+++ b/man/leaf_gradient_batch.Rd
@@ -58,9 +58,9 @@ the two routes to agree. `stem_c` has no such identity and always rebuilds.}
 \value{
 A list with
   \describe{
-    \item{`gradient`}{an `n` × `length(pars)` × 4 array, the last dimension
-      being `A`, `gc`, `psi_stem` and `collar`}
-    \item{`value`}{an `n` × 4 matrix of the solved outputs}
+    \item{`gradient`}{an `n` × `length(pars)` × 5 array, the last dimension
+      being `A`, `gc`, `psi_stem`, `collar` and `profit`}
+    \item{`value`}{an `n` × 5 matrix of the solved outputs}
     \item{`method`}{`"ift"` or `"fd"` per observation}
     \item{`status`}{`"interior"`, `"pinned"`, `"no-gradient"` or `"error"`}
     \item{`H`, `stationarity`}{the curvature and the implied Newton step
@@ -75,7 +75,7 @@ observation.
 }
 \section{What it computes, and where the maths is written down}{
 
-The same four derivatives at the same solved operating point, by the same two
+The same five derivatives at the same solved operating point, by the same two
 routes, with the same active-set test deciding between them. **Read
 [leaf_gradient()]'s documentation for all of it** — the implicit function
 theorem, why the second term is not a correction, why the premise is tested
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index a086b7e..de15051 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -2175,6 +2175,16 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// gradient_output_names
+std::vector<std::string> gradient_output_names();
+RcppExport SEXP _phylloptim_gradient_output_names() {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    rcpp_result_gen = Rcpp::wrap(gradient_output_names());
+    return rcpp_result_gen;
+END_RCPP
+}
 // gradient_batch_prepare
 SEXP gradient_batch_prepare(Rcpp::List root_network, Rcpp::List psi_soil, Rcpp::List soil_depth, Rcpp::NumericVector PPFD, Rcpp::NumericVector atm_vpd, Rcpp::NumericVector ca, Rcpp::NumericVector leaf_temp, Rcpp::NumericVector atm_o2_kpa, Rcpp::NumericVector atm_kpa);
 RcppExport SEXP _phylloptim_gradient_batch_prepare(SEXP root_networkSEXP, SEXP psi_soilSEXP, SEXP soil_depthSEXP, SEXP PPFDSEXP, SEXP atm_vpdSEXP, SEXP caSEXP, SEXP leaf_tempSEXP, SEXP atm_o2_kpaSEXP, SEXP atm_kpaSEXP) {
@@ -2428,6 +2438,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_phylloptim_Leaf__lambda_multilayer__get", (DL_FUNC) &_phylloptim_Leaf__lambda_multilayer__get, 1},
     {"_phylloptim_Leaf__g1_eff__get", (DL_FUNC) &_phylloptim_Leaf__g1_eff__get, 1},
     {"_phylloptim_gradient_par_names", (DL_FUNC) &_phylloptim_gradient_par_names, 0},
+    {"_phylloptim_gradient_output_names", (DL_FUNC) &_phylloptim_gradient_output_names, 0},
     {"_phylloptim_gradient_batch_prepare", (DL_FUNC) &_phylloptim_gradient_batch_prepare, 9},
     {"_phylloptim_gradient_batch_check", (DL_FUNC) &_phylloptim_gradient_batch_check, 1},
     {"_phylloptim_gradient_batch_run", (DL_FUNC) &_phylloptim_gradient_batch_run, 8},
diff --git a/src/gradient.cpp b/src/gradient.cpp
index cbcedce..a52b926 100644
--- a/src/gradient.cpp
+++ b/src/gradient.cpp
@@ -71,6 +71,27 @@ std::vector<std::string> gradient_par_names() {
   return phylloptim::gradient::par_names();
 }
 
+//' The differentiated outputs, in the order C++ indexes them
+//'
+//' `A`, `gc`, `psi_stem`, `collar` and `profit`. Unlike [gradient_par_names()],
+//' which R holds its own copy of and a test compares, R **reads** this one — so
+//' the list exists once and adding an output is one edit rather than two that
+//' can disagree.
+//'
+//' The first four are what a gas-exchange calibration observes; `profit` is what
+//' a demographic consumer bills, and it is the one output the envelope theorem
+//' reaches. See [leaf_gradient()].
+//'
+//' @return A character vector of five names.
+//' @seealso [leaf_gradient()], [leaf_gradient_batch()]
+//' @examples
+//' gradient_output_names()
+//' @export
+// [[Rcpp::export]]
+std::vector<std::string> gradient_output_names() {
+  return phylloptim::gradient::output_names();
+}
+
 // Convert N observations' drivers to C++ once. Everything is already resolved
 // and recycled to length N on the R side by `leaf_batch()`, which calls
 // `.resolve_drivers()` per row -- so the defaults (1 m layers, the nominal
diff --git a/tests/testthat/gradient_golden.tsv b/tests/testthat/gradient_golden.tsv
index bc5db9d..c705884 100644
--- a/tests/testthat/gradient_golden.tsv
+++ b/tests/testthat/gradient_golden.tsv
@@ -1,21 +1,21 @@
-case	status	method	par	A	gc	psi_stem	collar
-interior-1layer	interior	ift	vcmax_25	0x1.56ff4d6eecf96p-6	0x1.6279997564ff5p-15	0x1.2fe9f734f59d8p-8	0x1.08c7ad051fb7dp-9
-interior-1layer	interior	ift	stem_b	0x1.6c5b48adf2fp+1	0x1.7ab19a189bb21p-7	0x1.d9c2457829579p-1	0x1.1add81b8d5d2fp-1
-interior-1layer	interior	ift	psi_crit	0x0p+0	0x0p+0	0x0p+0	0x0p+0
-interior-1layer	interior	ift	R_d_25	-0x1.fd6a09155bf1cp-3	-0x1.90ddd402a632ep-12	-0x1.57b035f654ff8p-5	-0x1.2b6ec9ddabc73p-6
-interior-5layer	interior	ift	vcmax_25	0x1.a789fd22240c2p-5	0x1.f668e2d6ae349p-16	0x1.4bbb591c8bbecp-10	0x1.9618142da1224p-11
-interior-5layer	interior	ift	stem_b	0x1.bf615f57bdaf1p+0	0x1.40cec78217423p-6	0x1.35865423734e6p-1	0x1.03506a09d13edp-1
-interior-5layer	interior	ift	psi_crit	0x0p+0	0x0p+0	0x0p+0	0x0p+0
-interior-5layer	interior	ift	R_d_25	-0x1.53043af91afd9p-1	-0x1.75c8e3e4fefbp-10	-0x1.ed9b3452956e1p-5	-0x1.2e209f7e6b335p-5
-pinned-dry-3layer	pinned	fd	vcmax_25	0x1.1f8090c0306p-6	0x0p+0	0x0p+0	0x0p+0
-pinned-dry-3layer	pinned	fd	stem_b	0x1.2ccb4ca42ff14p+2	0x1.64b64bdad4946p-6	0x1.33ff011ff4998p-14	0x1.c3ef65a37fee1p-2
-pinned-dry-3layer	pinned	fd	psi_crit	0x1.41dade66b0dcp+0	0x1.7db02a11ba0cbp-8	0x1.fff9995bb6188p-1	0x1.e3a34a2bcbb21p-4
-pinned-dry-3layer	pinned	fd	R_d_25	-0x1.b3f6e43dfed39p-3	0x0p+0	0x0p+0	0x0p+0
-shutdown-1layer	no-gradient	fd	vcmax_25	0x0p+0	0x0p+0	0x0p+0	0x0p+0
-shutdown-1layer	no-gradient	fd	stem_b	0x0p+0	0x0p+0	0x0p+0	0x0p+0
-shutdown-1layer	no-gradient	fd	psi_crit	0x0p+0	0x0p+0	0x1.000000001aa4dp+0	0x1.000000001aa4dp+0
-shutdown-1layer	no-gradient	fd	R_d_25	-0x1.ffffffff73c8bp-1	0x0p+0	0x0p+0	0x0p+0
-single-potential	interior	ift	vcmax_25	0x1.8dbb95a5a7daap-11	0x1.8665f5c42e12fp-20	0x1.0f2ea755b7941p-11	0x1.f6f5831d8118cp-12
-single-potential	interior	ift	leaf_specific_conductance_max	0x1.61e2e05d142cfp+11	0x1.116a882643229p+3	0x1.c7141229c1e05p+10	0x1.603f57174c2c4p+11
-single-potential	interior	ift	resistance	-0x1.f8de9025064bep-14	-0x1.861138fe7344ap-22	-0x1.34162dac5854ep-14	-0x1.0ed249f513811p-14
-single-potential	interior	ift	R_d_25	-0x1.090f44c09b126p-5	-0x1.d8fdb4470b455p-15	-0x1.488db9f00ccp-6	-0x1.30aebd98fd4a3p-6
+case	status	method	par	A	gc	psi_stem	collar	profit
+interior-1layer	interior	ift	vcmax_25	0x1.56ff4d6eecf96p-6	0x1.6279997564ff5p-15	0x1.2fe9f734f59d8p-8	0x1.08c7ad051fb7dp-9	0x1.58f0981adp-7
+interior-1layer	interior	ift	stem_b	0x1.6c5b48adf2fp+1	0x1.7ab19a189bb21p-7	0x1.d9c2457829579p-1	0x1.1add81b8d5d2fp-1	0x1.6b78ea6fb8b06p+1
+interior-1layer	interior	ift	psi_crit	0x0p+0	0x0p+0	0x0p+0	0x0p+0	0x0p+0
+interior-1layer	interior	ift	R_d_25	-0x1.fd6a09155bf1cp-3	-0x1.90ddd402a632ep-12	-0x1.57b035f654ff8p-5	-0x1.2b6ec9ddabc73p-6	-0x1.3c91baea59239p-3
+interior-5layer	interior	ift	vcmax_25	0x1.a789fd22240c2p-5	0x1.f668e2d6ae349p-16	0x1.4bbb591c8bbecp-10	0x1.9618142da1224p-11	0x1.91a4ee5117eaap-5
+interior-5layer	interior	ift	stem_b	0x1.bf615f57bdaf1p+0	0x1.40cec78217423p-6	0x1.35865423734e6p-1	0x1.03506a09d13edp-1	0x1.203ce9bef1614p+1
+interior-5layer	interior	ift	psi_crit	0x0p+0	0x0p+0	0x0p+0	0x0p+0	0x0p+0
+interior-5layer	interior	ift	R_d_25	-0x1.53043af91afd9p-1	-0x1.75c8e3e4fefbp-10	-0x1.ed9b3452956e1p-5	-0x1.2e209f7e6b335p-5	-0x1.11dbf0824c12bp-1
+pinned-dry-3layer	pinned	fd	vcmax_25	0x1.1f8090c0306p-6	0x0p+0	0x0p+0	0x0p+0	0x1.1f8090c0306p-6
+pinned-dry-3layer	pinned	fd	stem_b	0x1.2ccb4ca42ff14p+2	0x1.64b64bdad4946p-6	0x1.33ff011ff4998p-14	0x1.c3ef65a37fee1p-2	0x1.750f766b59533p+2
+pinned-dry-3layer	pinned	fd	psi_crit	0x1.41dade66b0dcp+0	0x1.7db02a11ba0cbp-8	0x1.fff9995bb6188p-1	0x1.e3a34a2bcbb21p-4	0x1.03cb5d0db7589p-1
+pinned-dry-3layer	pinned	fd	R_d_25	-0x1.b3f6e43dfed39p-3	0x0p+0	0x0p+0	0x0p+0	-0x1.b3f6e43dfed39p-3
+shutdown-1layer	no-gradient	fd	vcmax_25	0x0p+0	0x0p+0	0x0p+0	0x0p+0	0x0p+0
+shutdown-1layer	no-gradient	fd	stem_b	0x0p+0	0x0p+0	0x0p+0	0x0p+0	0x1.21143bf2a3122p+0
+shutdown-1layer	no-gradient	fd	psi_crit	0x0p+0	0x0p+0	0x1.000000001aa4dp+0	0x1.000000001aa4dp+0	-0x1.7fef216dfc87fp-1
+shutdown-1layer	no-gradient	fd	R_d_25	-0x1.ffffffff73c8bp-1	0x0p+0	0x0p+0	0x0p+0	-0x1.00000000636f2p+0
+single-potential	interior	ift	vcmax_25	0x1.8dbb95a5a7daap-11	0x1.8665f5c42e12fp-20	0x1.0f2ea755b7941p-11	0x1.f6f5831d8118cp-12	0x1.222aede6bp-12
+single-potential	interior	ift	leaf_specific_conductance_max	0x1.61e2e05d142cfp+11	0x1.116a882643229p+3	0x1.c7141229c1e05p+10	0x1.603f57174c2c4p+11	0x1.1bcb9d6cef5e3p+10
+single-potential	interior	ift	resistance	-0x1.f8de9025064bep-14	-0x1.861138fe7344ap-22	-0x1.34162dac5854ep-14	-0x1.0ed249f513811p-14	-0x1.b3ad2d1358p-15
+single-potential	interior	ift	R_d_25	-0x1.090f44c09b126p-5	-0x1.d8fdb4470b455p-15	-0x1.488db9f00ccp-6	-0x1.30aebd98fd4a3p-6	-0x1.c00a32cb74c72p-7
diff --git a/tests/testthat/helper-golden.R b/tests/testthat/helper-golden.R
index 204fb8b..bbba82b 100644
--- a/tests/testthat/helper-golden.R
+++ b/tests/testthat/helper-golden.R
@@ -56,10 +56,18 @@ golden_tolerance <- function(field) {
 # already moved 1.3% while `vcmax_25` has moved 0.2%. That column goes
 # noise-dominated first, and it is the one that sets this number.
 #
-# 5e-03 leaves ~4x headroom over the one observation there is. It cannot hide a real
-# change: #41's own reallocation moved these cells by 16% to 250%, two orders above
-# this. ⚠️ It IS one CI observation, so read the worst-difference line the test
-# prints on every run rather than assuming the headroom is still there.
+# ⚠️ THE HEADROOM IS 2.1x, NOT THE ~4x THIS COMMENT USED TO CLAIM, and what halved it
+# was #87 adding the `profit` column. Same mechanism one step further: at a given row
+# `profit` is SMALLER in magnitude than `A` while carrying the same absolute floor, so
+# its relative disagreement is larger. Measured on Linux CI, the whole-file worst:
+#
+#   master (four columns)  1.30e-03
+#   with `profit`          2.34e-03
+#
+# It still cannot hide a real change -- #41's own reallocation moved these cells by 16%
+# to 250%, two orders above this. But the margin is now one order, not two, so read
+# the worst-difference line the test prints on every run; do not assume the headroom
+# is still there, and do not add a sixth column without re-reading it.
 gradient_golden_tolerance <- function() 5.0e-3
 
 # expect_identical where it can hold, expect_equal with the measured per-field
diff --git a/tests/testthat/test-gradient-batch.R b/tests/testthat/test-gradient-batch.R
index c421d7a..c0c77f2 100644
--- a/tests/testthat/test-gradient-batch.R
+++ b/tests/testthat/test-gradient-batch.R
@@ -227,7 +227,7 @@ test_that("the recorded gradients have not moved", {
     expect_identical(g$status[[1]], rows$status[[1]], label = nm)
     expect_identical(g$method[[1]], rows$method[[1]], label = nm)
     for (i in seq_len(nrow(rows))) {
-      for (out in c("A", "gc", "psi_stem", "collar")) {
+      for (out in c("A", "gc", "psi_stem", "collar", "profit")) {
         got <- g$gradient[1, rows$par[[i]], out]
         expect_golden(got, rows[[out]][[i]], out,
                       paste(nm, rows$par[[i]]),
@@ -272,7 +272,7 @@ test_that("an unsolvable row costs that row and not the batch", {
   solo <- vapply(ok, function(p) {
     as.vector(leaf_gradient_batch(leaf_batch(psi_soil = p, PPFD = 900),
                                   pars = pars, method = "ift")$gradient[1, , ])
-  }, numeric(length(pars) * 4L))
+  }, numeric(length(pars) * dim(g$gradient)[[3]]))
   for (i in seq_along(ok)) {
     expect_identical(as.vector(g$gradient[2 * i - 1, , ]), solo[, i],
                      label = paste("row", 2 * i - 1))
@@ -391,12 +391,13 @@ test_that("the result is shaped and named for a caller applying a Jacobian", {
   pars <- c("vcmax_25", "stem_b", "cost_scale_TF24")
   b <- leaf_batch(psi_soil = c(1.0, 1.5, 2.0, 2.5), PPFD = 900)
   g <- leaf_gradient_batch(b, pars = pars)
-  expect_identical(dim(g$gradient), c(4L, 3L, 4L))
+  expect_identical(dim(g$gradient), c(4L, 3L, 5L))
   expect_identical(dimnames(g$gradient)[[2]], pars)
   expect_identical(dimnames(g$gradient)[[3]],
-                   c("A", "gc", "psi_stem", "collar"))
-  expect_identical(dim(g$value), c(4L, 4L))
-  expect_identical(colnames(g$value), c("A", "gc", "psi_stem", "collar"))
+                   c("A", "gc", "psi_stem", "collar", "profit"))
+  expect_identical(dim(g$value), c(4L, 5L))
+  expect_identical(colnames(g$value),
+                   c("A", "gc", "psi_stem", "collar", "profit"))
   for (f in c("status", "method", "message")) {
     expect_length(g[[f]], 4L)
   }
diff --git a/tests/testthat/test-gradient.R b/tests/testthat/test-gradient.R
index 440a01e..a3e982b 100644
--- a/tests/testthat/test-gradient.R
+++ b/tests/testthat/test-gradient.R
@@ -220,7 +220,7 @@ test_that("a shut-down operating point reports no gradient and still differences
   # exactly here, so `dA/dR_d_25` is -1 and `dA/dvcmax_25` is EXACTLY zero --
   # vcmax_25 does not reach A at all at a shut-down point, so both perturbed solves
   # return the same bits. The -1 is a central difference and lands within ~6e-11.
-  g <- grid_gradient(6.0, pars = c("vcmax_25", "stem_b", "R_d_25"))
+  g <- grid_gradient(6.0, pars = c("vcmax_25", "stem_b", "R_d_25", "psi_crit"))
   expect_identical(g$status, "no-gradient")
   expect_identical(g$method, "fd")
   expect_equal(g$value[["A"]], -leaf_traits()$R_d_25)
@@ -228,6 +228,34 @@ test_that("a shut-down operating point reports no gradient and still differences
   expect_identical(g$gradient["vcmax_25", "A"], 0)
   expect_equal(g$gradient["stem_b", "A"], 0)
 
+  # ⚠️ `profit` NEEDS ITS OWN CLOSED FORM HERE, AND HAS ONE. This is the single
+  # regime where `profit_` is written by a branch that leaves the other outputs
+  # alone (hazard 8), so a recorded hex with nothing saying what it ought to be
+  # would pin a number rather than a fact. Two identities close it:
+  #
+  #   E = 0, so A = -R_d exactly and R_d(25 C) = R_d_25 -- and the hydraulic cost
+  #   does not depend on R_d_25 at all. So dprofit/dR_d_25 = dA/dR_d_25 = -1.
+  l <- leaf_model()
+  do.call(set_drivers, c(list(l), grid_drivers(6.0)))
+  l$find_root_collar_psi()
+  expect_identical(l$transpiration_, 0)
+  expect_identical(l$assim_colimited_, -l$R_d_)
+  expect_identical(g$value[["profit"]],
+                   l$assim_colimited_ - l$hydraulic_cost_)
+  expect_equal(g$gradient["R_d_25", "profit"], -1, tolerance = 1e-8)
+  expect_equal(g$gradient["R_d_25", "profit"], g$gradient["R_d_25", "A"],
+               tolerance = 1e-8)
+
+  #   The shut-down collar is PINNED AT psi_crit, so dcollar/dpsi_crit is exactly
+  #   1 -- which is why psi_crit, alone among the four, carries a non-zero profit
+  #   gradient here: it moves the collar, and the collar sets the hydraulic cost.
+  #   That is the whole explanation of a column that would otherwise look like
+  #   noise, and it is asserted rather than described.
+  expect_equal(g$gradient["psi_crit", "collar"], 1, tolerance = 1e-8)
+  expect_equal(g$gradient["psi_crit", "psi_stem"], 1, tolerance = 1e-8)
+  expect_lt(g$gradient["psi_crit", "profit"], 0)
+  expect_identical(g$gradient["psi_crit", "A"], 0)
+
   # Forcing the composite here is an error rather than a wrong number: unlike a
   # pinned point, there is no curvature to divide by at all.
   expect_error(grid_gradient(6.0, pars = "vcmax_25", method = "ift"),
@@ -370,7 +398,8 @@ test_that("the two non-trait parameters agree with a resolved reference", {
       a[[par]] <- value
     }
     x <- do.call(leaf_solve, a)
-    c(A = x$A, gc = x$gc, psi_stem = x$psi_stem, collar = x$collar)
+    c(A = x$A, gc = x$gc, psi_stem = x$psi_stem, collar = x$collar,
+      profit = x$profit)
   }
 
   g <- do.call(leaf_gradient,
@@ -495,16 +524,108 @@ test_that("leaf_gradient() rejects bad arguments", {
 
 test_that("the gradient is reported for every output the fit needs", {
   # leaf-calibration fits three responses -- A, gs and psi_leaf -- so all three
-  # are differentiated, not just A. `collar` comes along because it is psi*.
+  # are differentiated, not just A. `collar` comes along because it is psi*, and
+  # `profit` because it is what a DEMOGRAPHIC consumer bills: plant's carbon is
+  # `leaf.profit_`, not `assim_colimited_`, so without it the four columns above
+  # -- the calibration set -- were disjoint from plant's (#87).
   g <- grid_gradient(2.0, pars = "vcmax_25")
-  expect_identical(colnames(g$gradient), c("A", "gc", "psi_stem", "collar"))
-  expect_identical(names(g$value), c("A", "gc", "psi_stem", "collar"))
+  expect_identical(colnames(g$gradient),
+                   c("A", "gc", "psi_stem", "collar", "profit"))
+  expect_identical(names(g$value),
+                   c("A", "gc", "psi_stem", "collar", "profit"))
   expect_true(all(is.finite(g$gradient)))
   # Raising vcmax_25 raises assimilation and opens the stomata, and the leaf pays
   # for it with a more negative water potential (a larger positive magnitude).
   expect_gt(g$gradient["vcmax_25", "A"], 0)
   expect_gt(g$gradient["vcmax_25", "gc"], 0)
   expect_gt(g$gradient["vcmax_25", "psi_stem"], 0)
+  # More photosynthetic capacity is worth having: profit is A minus the hydraulic
+  # cost, and the extra carbon exceeds the extra cost. It is BELOW dA/dvcmax_25,
+  # which is the statement that the cost is not zero -- the cheapest available
+  # check that this column is the objective and not a copy of A.
+  expect_gt(g$gradient["vcmax_25", "profit"], 0)
+  expect_lt(g$gradient["vcmax_25", "profit"], g$gradient["vcmax_25", "A"])
+})
+
+test_that("profit's gradient is the direct term alone at an interior optimum", {
+  # The envelope theorem, which is the ONE place this package uses it. At an
+  # interior optimum dprofit/dpsi == 0, so the indirect term
+  # (dprofit/dpsi)(dpsi*/dtheta) vanishes IDENTICALLY, and dprofit/dtheta is the
+  # direct partial at fixed psi. `leaf_gradient()` encodes that by ASSIGNING the
+  # profit column from the direct term -- it does not zero dY_dpsi["profit"],
+  # which carries the exact dprofit/dpsi for the pinned route -- so the check is
+  # that the reported column equals a direct central difference with the collar
+  # HELD at psi*.
+  #
+  # ⚠️ THE OPERATING POINT IS CHOSEN, NOT ARBITRARY, and choosing it is what makes
+  # this a test rather than a formality. The assignment only matters where the term
+  # it removes is big enough to see, and that term is NOISE rather than an h^2
+  # truncation: `profit` is the maximum, so it is flat, and a central difference
+  # of it divides the solve's ~1e-9 floor by a ~1e-6 step. `?leaf_gradient` has
+  # the distribution over the golden grid's 136 interior rows; what matters here
+  # is that `psi_soil = 0.5, vpd = 2, 3 layers` is the WORST of them at 8.0e-05,
+  # so a tolerance three orders inside that fails if the assignment is removed.
+  # At the suite's usual `grid_drivers(2.0)` the same term is 2.9e-10 and this
+  # test would pass either way, which is the version of it written first.
+  #
+  # ⚠️ Those figures are macOS/arm64's, and the second half of this test says why
+  # that matters. Read it before adding an assertion on a magnitude here.
+  d <- grid_drivers(0.5, vpd = 2.0, layers = 3L)
+  g <- do.call(leaf_gradient, c(d, list(pars = "vcmax_25")))
+  expect_identical(g$status, "interior")
+
+  l <- leaf_model(leaf_traits(), leaf_control(), leaf_supply_multilayer())
+  do.call(set_drivers, c(list(l), d))
+  l$find_root_collar_psi()
+  psi_star <- l$opt_root_psi_
+
+  # The two instruments, both read on a leaf still at base traits.
+  hp <- max(abs(psi_star), 1) * 1e-6
+  exact <- l$dprofit_droot_collar_psi(psi_star)
+  l$evaluate_root_collar_psi(psi_star + hp)
+  hi <- l$profit_
+  l$evaluate_root_collar_psi(psi_star - hp)
+  lo <- l$profit_
+  fd <- (hi - lo) / (2 * hp)
+
+  # The answer itself, against a direct difference at fixed collar. THIS IS THE
+  # PART THAT HOLDS EVERYWHERE -- it is the claim the feature makes, and it does
+  # not depend on how big the term that was dropped happens to be.
+  v <- leaf_traits()$vcmax_25
+  h <- max(abs(v), 1) * 1e-6
+  at <- function(x) {
+    set_traits(l, leaf_traits(vcmax_25 = x))
+    do.call(set_drivers, c(list(l), d))
+    l$evaluate_root_collar_psi(psi_star)
+    l$profit_
+  }
+  direct <- (at(v + h) - at(v - h)) / (2 * h)
+  expect_equal(g$gradient["vcmax_25", "profit"], direct, tolerance = 1e-7)
+
+  # ⚠️ THE MAGNITUDES BELOW ARE PLATFORM-SPECIFIC, AND THE FIRST VERSION OF THIS
+  # TEST ASSERTED THEM EVERYWHERE. It passed on macOS/arm64 and failed on Linux
+  # CI, where the same operating point gives `exact` = 6.6e-11 and `fd` = 8.0e-10
+  # rather than 2.4e-15 and 2.1e-04.
+  #
+  # That is not a different answer, it is a different NOISE FLOOR: which side of
+  # the collar solver's tolerance the root-find lands on is set by libm's exp/pow,
+  # and those are not reproducible between Apple's arm64 libm and glibc on
+  # x86-64. The package already has a name for that -- it is the same reason the
+  # golden files compare bit-exactly on one platform and with a tolerance
+  # elsewhere -- so the same predicate gates it here.
+  #
+  # ⚠️ So this test has TEETH only on the platform the measurement was made on.
+  # Off it, the check above still holds and the one below is skipped; a Linux-only
+  # run would not catch the assignment being removed. Said plainly rather than left
+  # for someone to discover from a green CI.
+  skip_if_not(golden_bit_exact_platform(),
+              "the noise floor these two numbers measure is macOS/arm64's")
+  # The exact instrument says stationary to solver precision; the differenced one
+  # says nothing of the kind. Asserting BOTH is the point -- it is the difference
+  # between "the dropped term is small" and "the dropped term is unmeasurable by
+  # the route that would have supplied it".
+  expect_lt(abs(exact), 1e-12)
+  expect_gt(abs(fd), 1e-5)
 })
 
 # ---------------------------------------------------------------------------
diff --git a/tools/gradient_golden.R b/tools/gradient_golden.R
index d2671b1..9be9f64 100644
--- a/tools/gradient_golden.R
+++ b/tools/gradient_golden.R
@@ -98,6 +98,7 @@ out <- do.call(rbind, lapply(cases, function(cs) {
                gc = sprintf("%a", g$gradient[1, p, "gc"]),
                psi_stem = sprintf("%a", g$gradient[1, p, "psi_stem"]),
                collar = sprintf("%a", g$gradient[1, p, "collar"]),
+               profit = sprintf("%a", g$gradient[1, p, "profit"]),
                stringsAsFactors = FALSE)
   }))
 }))
diff --git a/vignettes/fitting.Rmd b/vignettes/fitting.Rmd
index 88684dc..e4f063a 100644
--- a/vignettes/fitting.Rmd
+++ b/vignettes/fitting.Rmd
@@ -108,6 +108,12 @@ data.frame(
 `leaf_gradient()` gives exact derivatives of each output, so the objective's
 gradient follows by the chain rule — no differencing of the solve.
 
+This fit names three of the five columns, and that is the normal case here: `A`,
+`gc` and `psi_stem` are measured, `collar` and `profit` are not. `profit` is in
+the matrix for a *demographic* consumer rather than a calibration one — `plant`
+bills carbon from it — so a fit against gas exchange simply ignores that column.
+It costs nothing to leave there; the whole row is read off one solve.
+
 ```{r}
 ssq_grad <- function(theta, l) {
   tr <- do.call(leaf_traits, as.list(exp(theta)))
diff --git a/vignettes/phylloptim.Rmd b/vignettes/phylloptim.Rmd
index c2631b6..b1a9946 100644
--- a/vignettes/phylloptim.Rmd
+++ b/vignettes/phylloptim.Rmd
@@ -374,10 +374,17 @@ tibble::tibble(
 So in a calibration, name only the parameters you are fitting. Asking for all
 fourteen when you fit three is several times the work.
 
-⚠️ **The four outputs are not selectable** — every call returns the derivative of
-`A`, `gc`, `psi_stem` and `collar` with respect to each requested parameter,
-because they are read off the same solve and cost nothing extra. `pars` chooses
-the rows of that matrix, never the columns.
+⚠️ **The five outputs are not selectable** — every call returns the derivative of
+`A`, `gc`, `psi_stem`, `collar` and `profit` with respect to each requested
+parameter, because they are read off the same solve and cost nothing extra.
+`pars` chooses the rows of that matrix, never the columns.
+
+The first four are what a gas-exchange calibration observes. `profit` is the
+objective, and it is there for a demographic consumer: `plant` bills carbon from
+the leaf's `profit`, not its assimilation. It is also the one output the envelope
+theorem reaches — at an interior optimum its gradient is the direct partial
+alone — which `?leaf_gradient` explains along with the measurement that makes
+that worth doing rather than merely true.
 
 Naming something undifferentiable is an error rather than a silent omission, and
 it tells you what is available: