diff --git a/R/data.processing.R b/R/data.processing.R
index 2b429b9..c7d67a6 100644
--- a/R/data.processing.R
+++ b/R/data.processing.R
@@ -16,26 +16,9 @@
   ## Safe upper bound for arguments passed to exp() -- exp(700) ~= 1.01e304,
   ## comfortably under .Machine$double.xmax (~1.80e308).
   EXP_CLIP_UPPER =  700,
-  ## What one OBSERVATION contributes when a parameter is outside its domain.
-  ## Two jobs pull against each other: be decisively worse than any legitimate
-  ## optimum, and stay small enough that optim()'s finite-difference gradient
-  ## across the boundary is still usable. Swept on five tHN fits that all reach
-  ## their domain guard (see notes/code_history/sfm.md); fits keeping their
-  ## standard errors, out of 5, against the SUMMED penalty:
-  ##   xmax 0/5 | 1e14 0/5 | 1e12 0/5 | 1e10 3/5 | 1e8 4/5 | 1e6 5/5 | 1e4 5/5
-  ## The fitted log-likelihood was IDENTICAL at every one, so the choice costs
-  ## no accuracy -- it decides only whether the Hessian is computable.
-  ##
-  ## PER OBSERVATION rather than a fixed total, because a fixed total keeps its
-  ## margin only while n is small: |logLik| grows with n (2861 at n = 1000 on
-  ## these data), so a total of 1e6 is 350x worse than the optimum at n = 1000
-  ## but only 3.5x at n = 1e5. Per observation the ratio is constant. Swept
-  ## across n, fits keeping standard errors:
-  ##             n = 200   n = 500   n = 1000
-  ##   1e2/obs     5/5       5/5       5/5
-  ##   1e4/obs     5/5       5/5       5/5
-  ##   1e6/obs     4/5       5/5       5/5
-  ## 1e4 is the largest that is 5/5 everywhere tested.
+  ## Per-observation penalty when a parameter leaves its domain: decisively
+  ## worse than any optimum, small enough to keep optim()'s gradient usable.
+  ## Swept; notes/code_history/data.processing.md.
   DOMAIN_PENALTY_PER_OBS = 1e4
   # HALTON_PRIMES = c(2, 3)
 )
diff --git a/R/sfareg_methods.R b/R/sfareg_methods.R
index f7bc55c..f3b0b57 100644
--- a/R/sfareg_methods.R
+++ b/R/sfareg_methods.R
@@ -118,12 +118,8 @@ vcov.sfareg <- function(object, type = c("hessian", "bhhh", "sandwich", "cluster
         call. = FALSE
       )
     }
-    ## Both failure paths below used to advise type = "bhhh" unconditionally,
-    ## on the reasoning that it needs no Hessian. That is true and not
-    ## sufficient: BHHH still has to invert crossprod(G), and the fits that
-    ## reach these errors are frequently flat in a parameter, which kills the
-    ## outer product too. Advising it blind sends the user to a second error.
-    ## So ask first, and only recommend what is actually available.
+    ## BHHH needs no Hessian but still inverts crossprod(G), which these flat
+    ## fits often kill too. notes/code_history/sfareg_methods.md.
     bhhh_ok <- !is.null(tryCatch(solve(crossprod(G)), error = function(e) NULL))
     bhhh_hint <- if (bhhh_ok) {
       " type = \"bhhh\" needs no Hessian and is defined here."
diff --git a/R/sfm.R b/R/sfm.R
index 4d06fce..3b5a163 100644
--- a/R/sfm.R
+++ b/R/sfm.R
@@ -950,15 +950,8 @@ sfm <- function(formula,
         nu <- x[3]
         like <- if (!is.finite(sig_v) || !is.finite(sig_u) || !is.finite(nu) ||
           sig_v <= 0 || sig_u <= 0 || nu <= 2) {
-          ## NOT .Machine$double.xmax / n: summed over n that IS double.xmax,
-          ## and optim() differences the objective, so the gradient overflows to
-          ## a genuinely non-finite value. Measured at n = 300: 5.99e305 per
-          ## observation, 1.7977e308 summed, (f - 295.5)/1e-3 = Inf. This guard
-          ## is reached in practice -- optim()'s ndeps = 1e-3 step takes sig_u
-          ## from its 1e-7 bound to -0.0009999 -- and instrumenting it gave a
-          ## 9/9 correlation between reaching it and losing every standard
-          ## error. The sibling branches' literal 1e12 does NOT fix that; see
-          ## .SFA_CONSTANTS$DOMAIN_PENALTY for the sweep behind this value.
+          ## NOT double.xmax / n: optim() differences the objective and the
+          ## gradient overflows. notes/code_history/sfm.md (A56a).
           rep(-.SFA_CONSTANTS$DOMAIN_PENALTY_PER_OBS, length(eps))
         } else {
           .log_d_thn(eps, sig_v, sig_u, nu)