diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md
index fd6837924..ed3cb87c7 100644
--- a/CHANGELOG-unreleased.md
+++ b/CHANGELOG-unreleased.md
@@ -20,6 +20,7 @@ the released changes.
 - `ssb_to_psb_xyz_ECL` and `ssb_to_psb_xyz_ICRS` are now cached
 - ELL1H with H3+STIGMA: add opt-in ``ell1h_shapiro="absorbed"`` to select Freire & Wex Eq. (28) (Tempo2 ELL1H/T2 mode 1). Default remains Eq. (29) ``"full"`` (`get_model` / `get_model_and_toas` / `ModelBuilder`).
 ### Fixed
+- Remove spurious ``/ Tsun`` factor from analytic DDH ``∂delay/∂STIGMA`` (design matrix / GLS for free ``STIGMA`` was wrong by ``1/Tsun`` since the Maple rewrite in PINT ≥ 1.0).
 - Align ``d_delayS3p_H3_STIGMA_exact_d_STIGMA`` with Eq. (28): ``cos(2*Phi)``.
 - Prefer DD over BT when guessing the binary model for Tempo2 `T2` par files (`allow_T2`), matching Tempo2's `allTerms=1` behavior
 - `WidebandTOAFitter` raises a warning if the model has correlated errors (It used to give wrong results before).
diff --git a/src/pint/models/stand_alone_psr_binaries/DDH_model.py b/src/pint/models/stand_alone_psr_binaries/DDH_model.py
index 1fdc40de5..4405d31a0 100644
--- a/src/pint/models/stand_alone_psr_binaries/DDH_model.py
+++ b/src/pint/models/stand_alone_psr_binaries/DDH_model.py
@@ -99,37 +99,18 @@ def d_delayS_d_par(self, par):
         cOmega = np.cos(self.omega())
         TM2 = self.M2.value * Tsun
 
-        logNum = (
-            1
-            - e * cE
-            - self.SINI * (sOmega * (cE - e) + (1 - e**2) ** 0.5 * cOmega * sE)
-        )
+        # Geometric factor in the Shapiro log argument (same as DD).
+        geo = sOmega * (cE - e) + (1 - e**2) ** 0.5 * cOmega * sE
+        logNum = 1 - e * cE - self.SINI * geo
         with u.set_enabled_equivalencies(u.dimensionless_angles()):
             dH3_dpar = self.prtl_der("H3", par)
             dsDelay_dH3 = -2 * np.log(logNum) / self.STIGMA**3
             dSTIGMA_dpar = self.prtl_der("STIGMA", par)
-            dsDelay_dSTIGMA = 6 * self.H3 / self.STIGMA**4 / Tsun.value * np.log(
-                1
-                - e * cE
-                - 2
-                * self.STIGMA
-                / (self.STIGMA**2 + 1)
-                * (sOmega * (cE - e) + (-(e**2) + 1) ** 0.5e0 * cOmega * sE)
-            ) - 2 * self.H3 / self.STIGMA**3 / Tsun.value * (
-                -2
-                / (self.STIGMA**2 + 1)
-                * (sOmega * (cE - e) + (-(e**2) + 1) ** 0.5e0 * cOmega * sE)
-                + 4
-                * self.STIGMA**2
-                / (self.STIGMA**2 + 1) ** 2
-                * (sOmega * (cE - e) + (-(e**2) + 1) ** 0.5e0 * cOmega * sE)
-            ) / (
-                1
-                - e * cE
-                - 2
-                * self.STIGMA
-                / (self.STIGMA**2 + 1)
-                * (sOmega * (cE - e) + (-(e**2) + 1) ** 0.5e0 * cOmega * sE)
+            dSINI_dSTIGMA = 2 * (1 - self.STIGMA**2) / (1 + self.STIGMA**2) ** 2
+            d_logNum_d_STIGMA = -dSINI_dSTIGMA * geo
+            dsDelay_dSTIGMA = (
+                6 * self.H3 / self.STIGMA**4 * np.log(logNum)
+                - 2 * self.H3 / self.STIGMA**3 * d_logNum_d_STIGMA / logNum
             )
 
             decc_dpar = self.prtl_der("ecc", par)
diff --git a/tests/test_ddh.py b/tests/test_ddh.py
index 40a9a1929..88849da9b 100644
--- a/tests/test_ddh.py
+++ b/tests/test_ddh.py
@@ -94,3 +94,29 @@ def test_ddh_sim():
     assert np.isclose(f.resids.calc_chi2(), f2.resids.calc_chi2(), atol=0.5)
     assert np.isclose(f.model.M2.value, f2.model.M2.value, atol=0.01)
     assert np.isclose(f.model.SINI.value, f2.model.SINI.value, atol=0.01)
+
+
+def test_ddh_stigma_derivative_matches_finite_difference():
+    """Analytic ∂delay/∂STIGMA (and H3) must match PINT's numerical derivative.
+
+    Regression for a spurious ``/ Tsun.value`` factor formerly present in
+    ``DDHmodel.d_delayS_d_par`` for STIGMA (analytic column was ~1/Tsun too
+    large). H3 is a non-regression control. Compare RMS ratios so a few
+    superior-conjunction TOAs do not dominate the check.
+    """
+    m = get_model(
+        io.StringIO(
+            f"{parDD}\nBINARY DD\nSINI {np.sin(i).value}\nA1 {A1.value}\n"
+            f"PB {PB.value}\nM2 {Mc.value}\n"
+        )
+    )
+    t = pint.simulation.make_fake_toas_uniform(
+        50000, 51000, 200, m, add_noise=False, error=1 * u.us
+    )
+    m2 = pint.binaryconvert.convert_binary(m, "DDH")
+
+    for param in ("STIGMA", "H3"):
+        analytic = m2.d_delay_d_param(t, param)
+        numerical = m2.d_delay_d_param_num(t, param)
+        rms_ratio = np.sqrt(np.mean(analytic**2)) / np.sqrt(np.mean(numerical**2))
+        assert abs(rms_ratio - 1.0) < 1e-3