diff --git a/docs/joint-command-gain-prototype.md b/docs/joint-command-gain-prototype.md
new file mode 100644
index 0000000..fb9daf3
--- /dev/null
+++ b/docs/joint-command-gain-prototype.md
@@ -0,0 +1,29 @@
+# Joint phase and command-response estimation: development prototype
+
+This draft adds an offline estimator, not a runtime controller or a new claimed scientific method. The response parameter is a positive scalar multiplying commanded OPD. It is not the optical gain of a pyramid wavefront sensor.
+
+For nominal command u_i, the forward model is I(a - g u_i). The existing Poisson-Gaussian camera likelihood is reused while fitting six phase coefficients a and one scalar response g. The final nominal correction is clipped(a_hat / g_hat). True phase and physical response are never estimator inputs.
+
+`NominalObservation` copies and freezes camera arrays and nominal command values. Its legacy `applied_opd_nm` alias is explicitly nominal: it exists only to feed the unchanged baseline estimator without leaking physical response labels. This is not a claim that a requested command was physically applied.
+
+## Executed checks
+
+The added test file passed 16 tests in a local focused environment using byte-verified copies of the existing optical model, likelihood and nominal estimator at commit `9406f5b5b04f9692cd194599277109423a2dbee6`. Checks include all seven gradient coordinates, reduction to the nominal objective at unit gain, input immutability, refusal of unchanging commands, compensation algebra and a bright synthetic recovery fixture.
+
+```sh
+PYTHONPATH=src python -m pytest tests/test_gain_estimation.py -q
+```
+
+The entire repository suite, supported-version matrix, CI and hardware execution were not run in that local environment. Passing these tests does not qualify the prototype for deployment.
+
+## Development evidence and promotion boundary
+
+A separate retained local experiment compares unchanged nominal fixed and adaptive phase-diversity estimators with this joint estimator. It is not integrated into the repository's release or reserved evaluation cohorts. No existing controller, policy, result, frozen data or default behavior is changed by this draft.
+
+Before promoting scientific claims in the repository, port and review the complete acquisition harness and raw artifacts, independently reproduce the comparison, and admit a new disjoint cohort. The local study uses a static six-mode, matched pupil and detector model, known photon flux, constant scalar gain and offline fitting. It does not establish performance under unknown flux, unmodeled optical aberrations, temporal gain variation, misregistration, mode-dependent response or real instrument measurements.
+
+## Known limitations
+
+Command diversity and bounded gain do not guarantee global identifiability. Three local starts are not proof of global optimality. The finite best fit is returned even when its optimizer reports failure, matching the old estimator's convention; every caller must inspect diagnostics. No confidence-calibrated withholding rule or safe command-authority gate is provided. The nominal model may perform better when calibration is already correct because the joint model estimates an additional parameter.
+
+Phase diversity and deformable-mirror calibration have substantial prior art. This is a targeted engineering extension to ProbeBench, not a novelty or learned-control claim.
diff --git a/src/probebench/gain_estimation.py b/src/probebench/gain_estimation.py
new file mode 100644
index 0000000..5d2fe25
--- /dev/null
+++ b/src/probebench/gain_estimation.py
@@ -0,0 +1,138 @@
+"""Joint phase/command-response estimation; no evaluator truth is accepted.
+
+The gain here multiplies commanded OPD. It is NOT pyramid-WFS optical gain.
+This is an offline calibration prototype, not a qualified runtime controller.
+"""
+from __future__ import annotations
+
+from dataclasses import dataclass
+from time import perf_counter
+from typing import Sequence
+import numpy as np
+from scipy.optimize import minimize
+from .likelihood import poisson_gaussian_nll
+from .optics import coefficients, OpticalModel
+
+
+@dataclass(frozen=True)
+class NominalObservation:
+    """Camera data and nominal commands only; never physical response labels."""
+    nominal_command_nm: np.ndarray
+    image_electrons: np.ndarray
+    expected_photons: float
+
+    def __post_init__(self):
+        command = coefficients(self.nominal_command_nm)
+        image = np.array(self.image_electrons, dtype=float, copy=True)
+        if image.ndim != 2 or not image.size or not np.isfinite(image).all():
+            raise ValueError('A finite two-dimensional camera array is required')
+        if not np.isfinite(self.expected_photons) or self.expected_photons <= 0:
+            raise ValueError('Expected photons must be positive and finite')
+        command.setflags(write=False)
+        image.setflags(write=False)
+        object.__setattr__(self, 'nominal_command_nm', command)
+        object.__setattr__(self, 'image_electrons', image)
+
+    @property
+    def applied_opd_nm(self):
+        """Legacy estimator adapter: nominal value, not actual physical OPD.
+
+        The legacy estimator assumes calibration is correct. This property exists
+        only to compare it on the same camera data without leaking true gain.
+        """
+        return self.nominal_command_nm
+
+
+def joint_objective(model: OpticalModel, observations: Sequence[NominalObservation],
+                    read_noise_electrons: float, parameters: np.ndarray):
+    """Exact Poisson-Gaussian NLL; x[:6]=phase/100, x[6]=(gain-1)/0.1."""
+    x = np.asarray(parameters, dtype=float)
+    if x.shape != (7,) or not np.isfinite(x).all():
+        raise ValueError('Seven finite scaled parameters are required')
+    phase, gain = 100*x[:6], 1 + 0.1*x[6]
+    photons = sum(o.expected_photons for o in observations)
+    if photons <= 0 or gain <= 0:
+        raise ValueError('Positive photons and gain are required')
+    cost, gradient = 0.0, np.zeros(7)
+    for obs in observations:
+        intensity, jac = model.intensity_and_jacobian(phase-gain*obs.nominal_command_nm)
+        value, pixel_gradient = poisson_gaussian_nll(
+            obs.expected_photons*intensity, obs.image_electrons, read_noise_electrons)
+        derivative = obs.expected_photons*np.einsum('ij,mij->m', pixel_gradient, jac)
+        cost += value
+        gradient[:6] += 100*derivative
+        gradient[6] -= 0.1*np.dot(obs.nominal_command_nm, derivative)
+    return cost/photons, gradient/photons
+
+
+def estimate_phase_and_gain(model: OpticalModel,
+                            observations: Sequence[NominalObservation],
+                            read_noise_electrons: float,
+                            coefficient_limit_nm: float = 200.0,
+                            gain_bounds: tuple[float, float] = (0.7, 1.3),
+                            max_iterations: int = 200):
+    """Fit six phase coefficients and one positive scalar response from images.
+
+    A bounded gain prior and command diversity are required, but do not guarantee
+    global identifiability. All local-fit outcomes are returned. The finite best
+    fit is retained even when its optimizer reports failure, matching the existing
+    nominal estimator; callers must check optimizer_success before deployment.
+    """
+    if len(observations) < 2:
+        raise ValueError('At least two observations are required')
+    if (not np.isfinite(coefficient_limit_nm) or coefficient_limit_nm <= 0
+            or type(max_iterations) is not int or max_iterations <= 0):
+        raise ValueError('Positive coefficient and iteration limits are required')
+    lo, hi = gain_bounds
+    if not (np.isfinite([lo, hi]).all() and 0 < lo < 1 < hi):
+        raise ValueError('Gain bounds must be finite, positive, and contain one')
+    commands = np.array([o.nominal_command_nm for o in observations])
+    if np.all(commands == commands[0]):
+        raise ValueError('Unchanging commands cannot separate phase from command gain')
+    expected_shape = (model.config.pupil_pixels*model.config.oversampling,)*2
+    if any(o.image_electrons.shape != expected_shape for o in observations):
+        raise ValueError('Camera shape does not match the optical model')
+    if not np.isfinite(read_noise_electrons) or read_noise_electrons < 0:
+        raise ValueError('Read noise must be finite and nonnegative')
+    bounds = [(-coefficient_limit_nm/100, coefficient_limit_nm/100)]*6
+    bounds += [((lo-1)/0.1, (hi-1)/0.1)]
+    offset = min(0.5, coefficient_limit_nm/200)
+    fits, records = [], []
+    started = perf_counter()
+    for sign in (0, 1, -1):
+        initial = np.zeros(7)
+        initial[0] = sign*offset
+        fit = minimize(lambda x: joint_objective(model, observations, read_noise_electrons, x),
+                       initial, jac=True, method='L-BFGS-B', bounds=bounds,
+                       options={'maxiter': max_iterations, 'ftol': 1e-12,
+                                'gtol': 1e-8, 'maxls': 30})
+        finite = bool(np.isfinite(fit.fun) and np.isfinite(fit.x).all())
+        fits.append(fit)
+        records.append({'initial_opd_nm':(100*initial[:6]).tolist(),
+                        'initial_gain':1.0, 'success':bool(fit.success),
+                        'finite':finite, 'cost':float(fit.fun) if finite else None,
+                        'parameters':fit.x.tolist() if finite else None,
+                        'function_evaluations':int(fit.nfev),
+                        'iterations':int(fit.nit), 'message':str(fit.message)})
+    valid = [i for i, r in enumerate(records) if r['finite']]
+    if not valid:
+        raise FloatingPointError('No finite joint fit; retain this episode as a failure')
+    chosen = min(valid, key=lambda i:fits[i].fun)
+    best = fits[chosen]
+    gain = float(1+0.1*best.x[6])
+    diagnostics = {'optimizer_success':bool(best.success), 'chosen_start':chosen,
+                   'cost':float(best.fun), 'starts':records,
+                   'reported_function_evaluations':sum(r['function_evaluations'] for r in records),
+                   'fit_elapsed_seconds':perf_counter()-started,
+                   'gain_at_bound':bool(min(gain-lo, hi-gain) < 1e-5),
+                   'gain_bounds':list(gain_bounds),
+                   'objective':'poisson_gaussian_convolution_nll_per_expected_photon',
+                   'qualification':'development_only_not_global_identifiability_or_runtime_assurance'}
+    return 100*best.x[:6], gain, diagnostics
+
+
+def compensated_command(phase_nm: np.ndarray, gain: float, limit_nm: float = 200.0):
+    """Convert estimated physical phase to nominal actuator command, with bounds."""
+    if not np.isfinite(gain) or gain <= 0 or not np.isfinite(limit_nm) or limit_nm <= 0:
+        raise ValueError('Gain and command limit must be positive and finite')
+    return np.clip(coefficients(phase_nm)/gain, -limit_nm, limit_nm)
diff --git a/tests/test_gain_estimation.py b/tests/test_gain_estimation.py
new file mode 100644
index 0000000..4bda758
--- /dev/null
+++ b/tests/test_gain_estimation.py
@@ -0,0 +1,83 @@
+import numpy as np
+import pytest
+from probebench.optics import OpticalModel, OpticsConfig
+from probebench.gain_estimation import (NominalObservation, joint_objective,
+                                       estimate_phase_and_gain, compensated_command)
+from probebench.estimation import phase_objective
+
+
+def fixture(read_noise=1.0):
+    model = OpticalModel(OpticsConfig(pupil_pixels=12))
+    phase = np.array([15., 10., -3., 4., -5., 2.])
+    rng = np.random.default_rng(104)
+    obs = []
+    for d in (0., 80., -80.):
+        command = np.array([d, 0, 0, 0, 0, 0])
+        mean = 20000*model.intensity(phase-0.9*command)
+        image = rng.poisson(mean).astype(float) + rng.normal(0, read_noise, mean.shape)
+        obs.append(NominalObservation(command, image, 20000))
+    return model, obs
+
+
+@pytest.mark.parametrize('index', range(7))
+def test_joint_gradient_matches_central_difference(index):
+    model, obs = fixture()
+    x = np.array([.12,.08,-.02,.03,-.04,.01,-.5])
+    _, grad = joint_objective(model, obs, 1., x)
+    step = np.eye(7)[index]*1e-5
+    numerical = (joint_objective(model, obs, 1., x+step)[0]
+                 - joint_objective(model, obs, 1., x-step)[0])/2e-5
+    assert grad[index] == pytest.approx(numerical, rel=3e-5, abs=1e-7)
+
+
+def test_at_unit_gain_exactly_matches_existing_nominal_objective():
+    model, obs = fixture()
+    x = np.array([.12,.08,-.02,.03,-.04,.01,0.])
+    new, grad = joint_objective(model, obs, 1., x)
+    old, oldgrad = phase_objective(model, obs, 1., x[:6])
+    assert new == old
+    np.testing.assert_allclose(grad[:6], oldgrad, rtol=2e-14, atol=2e-15)
+
+
+def test_no_diversity_refused():
+    model, obs = fixture()
+    with pytest.raises(ValueError, match='Unchanging'):
+        estimate_phase_and_gain(model, [obs[0], obs[0]], 1.)
+
+
+def test_snapshot_does_not_mutate_or_accept_truth():
+    model, obs = fixture()
+    assert not hasattr(obs[0], 'true_gain')
+    assert not obs[0].image_electrons.flags.writeable
+    with pytest.raises(ValueError):
+        obs[0].nominal_command_nm[0] = 20
+    with pytest.raises(TypeError):
+        NominalObservation(np.zeros(6), np.ones((12,12)), 1., true_gain=.9)
+
+
+@pytest.mark.parametrize('gain', [0., -1., np.nan])
+def test_invalid_compensation_refused(gain):
+    with pytest.raises(ValueError):
+        compensated_command(np.ones(6), gain)
+
+
+def test_gain_compensation_cancels_ideal_plant_and_respects_limit():
+    phase = np.array([15., 10., -3., 4., -5., 2.])
+    np.testing.assert_allclose(phase-.9*compensated_command(phase,.9), 0., atol=1e-14)
+    assert np.max(np.abs(compensated_command(phase,.01,10))) <= 10
+
+
+def test_power_and_finite_pupil_wfe():
+    model, _ = fixture()
+    phase = np.array([15., 10., -3., 4., -5., 2.])
+    assert model.intensity(phase).sum() == pytest.approx(1., abs=2e-14)
+    wf = np.einsum('m,mij->ij', phase, model.basis)[model.pupil]
+    assert model.evaluate(phase)['residual_wfe_nm'] == pytest.approx(np.std(wf))
+
+
+def test_joint_fit_recovers_response_in_bright_development_fixture():
+    model, obs = fixture()
+    phase, gain, diag = estimate_phase_and_gain(model, obs, 1.)
+    assert diag['optimizer_success']
+    assert gain == pytest.approx(.9, abs=.025)
+    assert model.evaluate(np.array([15.,10.,-3.,4.,-5.,2.])-.9*compensated_command(phase,gain))['residual_wfe_nm'] < 2.