Support hierarchical triple systems with dual binary components (BINARY + BINARY2)

CHANGELOG-unreleased.md

@@ -17,6 +17,7 @@ the released changes.
 ### Added
 - Plot whitened DM residuals in pintk.
 - `ssb_to_psb_xyz_ECL` and `ssb_to_psb_xyz_ICRS` are now cached
+- Support for hierarchical triple systems: a second (outer) binary component can be added via a `BINARY2` line with `_2`-suffixed orbital parameters (e.g. `PB_2`, `A1_2`). Outer orbit delay is computed before, and propagated into, the inner binary. Outer wrappers `BinaryDD2` and `BinaryBT2` are provided.
 ### Fixed
 - `WidebandTOAFitter` raises a warning if the model has correlated errors (It used to give wrong results before).
 - Fixed bug where "include_bipm" flag was being ignored when loading Fermi TOAs with weights, now defaults to using EPHEM, CLOCK and PLANET_SHAPIRO from the timing model

src/pint/models/__init__.py

@@ -21,8 +21,8 @@
 
 # Import all standard model components here
 from pint.models.astrometry import AstrometryEcliptic, AstrometryEquatorial
-from pint.models.binary_bt import BinaryBT, BinaryBTPiecewise
-from pint.models.binary_dd import BinaryDD, BinaryDDGR, BinaryDDH, BinaryDDS
+from pint.models.binary_bt import BinaryBT, BinaryBT2, BinaryBTPiecewise
+from pint.models.binary_dd import BinaryDD, BinaryDD2, BinaryDDGR, BinaryDDH, BinaryDDS
 from pint.models.binary_ddk import BinaryDDK
 from pint.models.binary_ell1 import BinaryELL1, BinaryELL1H, BinaryELL1k
 from pint.models.chromatic_model import ChromaticCM, ChromaticCMX

src/pint/models/binary_bt.py

@@ -66,18 +66,46 @@ def validate(self):
         """Validate BT model parameters"""
         super().validate()
         for p in ("T0", "A1"):
-            if getattr(self, p).value is None:
+            if self._bp(p).value is None:
                 raise MissingParameter("BT", p, f"{p} is required for BT")
 
         # If any *DOT is set, we need T0
         for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT"):
-            if getattr(self, p).value is None:
-                getattr(self, p).value = "0"
-                getattr(self, p).frozen = True
+            if self._bp(p).value is None:
+                self._bp(p).value = "0"
+                self._bp(p).frozen = True
 
-        if self.GAMMA.value is None:
-            self.GAMMA.value = "0"
-            self.GAMMA.frozen = True
+        if self._bp("GAMMA").value is None:
+            self._bp("GAMMA").value = "0"
+            self._bp("GAMMA").frozen = True
+
+
+class BinaryBT2(BinaryBT):
+    """Outer-orbit Blandford and Teukolsky model for a hierarchical triple.
+
+    This is identical to :class:`pint.models.binary_bt.BinaryBT` except that all
+    of its parameters carry a ``_2`` suffix (``PB_2``, ``A1_2``, ``T0_2``, ...)
+    and it is selected with the ``BINARY2`` parfile parameter instead of
+    ``BINARY``. See :class:`pint.models.binary_dd.BinaryDD2` for a description of
+    how the outer orbit couples into the inner binary.
+
+    Orbital-frequency (``FBn``) and ``ORBWAVE`` parameterizations are not
+    supported for the outer orbit.
+
+    Parameters supported:
+
+    .. paramtable::
+        :class: pint.models.binary_bt.BinaryBT2
+    """
+
+    register = True
+    category = "pulsar_system_outer"
+    param_suffix = "_2"
+    binary_param_tag = "BINARY2"
+
+    def __init__(self):
+        super().__init__()
+        self._apply_param_suffix()
 
 
 class BinaryBTPiecewise(PulsarBinary):

src/pint/models/binary_dd.py

@@ -115,21 +115,59 @@ def validate(self):
         self.check_required_params(["T0", "A1"])
         # If any *DOT is set, we need T0
         for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT"):
-            if hasattr(self, p) and getattr(self, p).value is None:
-                getattr(self, p).value = 0.0
-                getattr(self, p).frozen = True
+            if self._hasbp(p) and self._bp(p).value is None:
+                self._bp(p).value = 0.0
+                self._bp(p).frozen = True
 
-        if hasattr(self, "GAMMA") and self.GAMMA.value is None:
-            self.GAMMA.value = 0.0
-            self.GAMMA.frozen = True
+        if self._hasbp("GAMMA") and self._bp("GAMMA").value is None:
+            self._bp("GAMMA").value = 0.0
+            self._bp("GAMMA").frozen = True
 
         # If eccentricity is zero, freeze some parameters to 0
         # OM = 0 -> T0 = TASC
-        if self.ECC.value == 0 or self.ECC.value is None:
+        if self._bp("ECC").value == 0 or self._bp("ECC").value is None:
             for p in ("ECC", "OM", "OMDOT", "EDOT"):
-                if hasattr(self, p):
-                    getattr(self, p).value = 0.0
-                    getattr(self, p).frozen = True
+                if self._hasbp(p):
+                    self._bp(p).value = 0.0
+                    self._bp(p).frozen = True
+
+
+class BinaryDD2(BinaryDD):
+    """Outer-orbit Damour and Deruelle model for a hierarchical triple system.
+
+    This is identical to :class:`pint.models.binary_dd.BinaryDD` except that all
+    of its parameters carry a ``_2`` suffix (``PB_2``, ``A1_2``, ``T0_2``, ...)
+    and it is selected with the ``BINARY2`` parfile parameter instead of
+    ``BINARY``. It is intended to model the *outer* orbit of a hierarchical
+    triple, alongside a normal inner binary component.
+
+    Because this component belongs to the ``pulsar_system_outer`` category,
+    which is ordered before ``pulsar_system`` in
+    :data:`pint.models.timing_model.DEFAULT_ORDER`, its delay is accumulated
+    before the inner binary's delay. PINT evaluates each binary at
+    ``barycentric time - accumulated delay``, so the outer orbit's light-travel
+    delay automatically shifts the epoch at which the inner orbit is evaluated.
+    This reproduces the physical coupling of a hierarchical triple (the wide
+    outer orbit Doppler-shifting the inner orbit), rather than naively adding two
+    independent binary delays.
+
+    Orbital-frequency (``FBn``) and ``ORBWAVE`` parameterizations are not
+    supported for the outer orbit.
+
+    Parameters supported:
+
+    .. paramtable::
+        :class: pint.models.binary_dd.BinaryDD2
+    """
+
+    register = True
+    category = "pulsar_system_outer"
+    param_suffix = "_2"
+    binary_param_tag = "BINARY2"
+
+    def __init__(self):
+        super().__init__()
+        self._apply_param_suffix()
 
 
 class BinaryDDS(BinaryDD):

src/pint/models/model_builder.py

@@ -271,8 +271,8 @@ def _validate_components(self):
                     f" class, please set register to 'False' in the class"
                     f" of component {k}."
                 )
-                if v.category == "pulsar_system":
-                    # The pulsar system will be selected by parameter BINARY
+                if v.category in ("pulsar_system", "pulsar_system_outer"):
+                    # The pulsar system is selected by parameter BINARY/BINARY2
                     continue
                 else:
                     raise ComponentConflict(m)
@@ -488,6 +488,12 @@ def choose_model(self, param_inpar, force_binary_model=None, allow_T2=False):
                 self.choose_binary_model(param_inpar, force_binary_model, allow_T2)
             )
 
+        # Outer-orbit binary for hierarchical triple systems.
+        binary2 = param_inpar.get("BINARY2", None)
+        if binary2:
+            binary2 = binary2[0]
+            selected_components.add(self.choose_outer_binary_model(param_inpar))
+
         # 2. Get the component list from the parameters in the parfile.
         # 2.1 Check the aliases of input parameters.
         # This does not include the repeating parameters, but it should not
@@ -531,6 +537,17 @@ def choose_model(self, param_inpar, force_binary_model=None, allow_T2=False):
                     )
                 else:
                     continue
+            # Outer-orbit binary component, controlled by the BINARY2 tag.
+            if self.all_components.components[cps[0]].category == "pulsar_system_outer":
+                if binary2 is None:
+                    raise MissingBinaryError(
+                        f"The outer-orbit binary model is decided by the"
+                        f" parameter 'BINARY2'. Please indicate the outer"
+                        f" binary model before using parameter {k}, which"
+                        f" is an outer binary model parameter."
+                    )
+                else:
+                    continue
 
             if len(cps) == 1:  # No conflict, parameter only shows in one component.
                 selected_components.add(cps[0])
@@ -645,6 +662,32 @@ def choose_binary_model(self, param_inpar, force_binary_model=None, allow_T2=Fal
 
         return binary_cp.__class__.__name__
 
+    def choose_outer_binary_model(self, param_inpar):
+        """Choose the outer-orbit BINARY2 model for a hierarchical triple.
+
+        Parameters
+        ----------
+        param_inpar: dict
+            Dictionary of the unique parameters in .par file with the key being
+            the parfile line. :func:`parse_parfile` returns this dictionary.
+
+        Returns
+        -------
+        str
+            Name of the outer binary component class.
+
+        Note
+        ----
+        The outer model is taken verbatim from the ``BINARY2`` parameter (no T2
+        guessing is performed) and must correspond to a registered component in
+        the ``pulsar_system_outer`` category (e.g. ``DD`` -> ``BinaryDD2``).
+        """
+        binary2 = param_inpar["BINARY2"][0]
+        binary_cp = self.all_components.search_binary_components(
+            binary2, category="pulsar_system_outer"
+        )
+        return binary_cp.__class__.__name__
+
     def _setup_model(
         self,
         timing_model,