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

CHANGELOG-unreleased.md

@@ -15,9 +15,11 @@ the released changes.
 - Updated GMRT coordinates.
 - Replaced custom ``pint.ls`` with astropy ``u.lsec``
 - Updated code to remove deprecation warnings during CI
+- `d_delay_d_param` now applies the chain rule for delay components that respond to delays accumulated from earlier components (e.g. a binary delay's dependence on its evaluation epoch). Design-matrix entries for binary pulsars change at the ~1e-5 relative level. The parameter-independent ingredients are computed once per design matrix via the new `TimingModel.delay_deriv_chain`; the unused `acc_delay` argument of `d_delay_d_param` was replaced by the optional `chain` argument.
 ### 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`, `BinaryBT2`, and `BinaryELL12` are provided. Delay derivatives account for the outer→inner coupling (chain rule through the previous delay). The projected semi-major axis includes a second time derivative `A1DOT2` (alias `X2DOT`).
 ### 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
@@ -36,4 +38,5 @@ the released changes.
 - Make VLBI frame rotation work correctly when proper motion is present.
 - Changed some API to pass Mac CI
 - Log-separated frequency computation for red noise components.
+- Place ``solar_windx`` before the binary in ``DEFAULT_ORDER`` so SolarWindDispersionX delays and derivatives chain-rule through the binary the same way as ``solar_wind``.
 ### Removed

docs/explanation.rst

@@ -588,6 +588,9 @@ components.
 
    - If the ``BINARY`` line is present in the parameter file, its value
      determines which binary model to use; if not, no binary model is used.
+     An optional ``BINARY2`` line selects an outer-orbit component for a
+     hierarchical triple (parameters are suffixed with ``_2``, e.g. ``PB_2``;
+     outer wrappers include ``BinaryDD2``, ``BinaryBT2``, and ``BinaryELL12``).
    - Each model component has one or more "special parameters" or families of
      parameters identified by a common prefix. If a par file contains a special
      parameter, or a known alias of one, then the timing model uses the

docs/timingmodels.rst

@@ -13,8 +13,10 @@ file, these are selected based on the parameters present. Binary models are
 selected explicitly using the ``BINARY`` parameter, while each non-binary
 component is selected if some parameter unique to it is included (for example
 if ``ELAT`` is present, :class:`~pint.models.astrometry.AstrometryEcliptic` is
-selected). Ambiguous or contradictory parameter files are possible, and for
-these PINT raises an exception.
+selected). Hierarchical triples may also include an outer orbit selected with
+``BINARY2`` (parameters carry a ``_2`` suffix, e.g. ``PB_2``). Ambiguous or
+contradictory parameter files are possible, and for these PINT raises an
+exception.
 
 .. componentlist::
 

src/pint/binaryconvert.py

@@ -7,7 +7,7 @@
 """
 
 import copy
-from typing import List, Optional, Tuple, Union
+from typing import Tuple
 
 import numpy as np
 from astropy import units as u
@@ -598,9 +598,16 @@ def convert_binary(
     if not model.is_binary:
         raise AttributeError("Input model is not a binary")
 
-    binary_component_name = [
+    binary_component_names = [
         x for x in model.components.keys() if x.startswith("Binary")
-    ][0]
+    ]
+    if len(binary_component_names) > 1:
+        raise ValueError(
+            "convert_binary does not support hierarchical triple systems "
+            f"with multiple binary components ({binary_component_names}); "
+            "convert each orbit separately or remove BINARY2 first."
+        )
+    binary_component_name = binary_component_names[0]
     binary_component = model.components[binary_component_name]
     if binary_component.binary_model_name == output:
         log.debug(

src/pint/models/__init__.py

@@ -21,10 +21,10 @@
 
 # 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.binary_ell1 import BinaryELL1, BinaryELL1H, BinaryELL1k, BinaryELL12
 from pint.models.chromatic_model import ChromaticCM, ChromaticCMX
 from pint.models.cmwavex import CMWaveX
 from pint.models.dispersion_model import (

src/pint/models/binary_bt.py

@@ -66,18 +66,47 @@ def validate(self):
         """Validate BT model parameters"""
         super().validate()
         for p in ("T0", "A1"):
-            if getattr(self, p).value is None:
-                raise MissingParameter("BT", p, f"{p} is required for BT")
+            if self._bp(p).value is None:
+                pname = self._bp(p).name
+                raise MissingParameter("BT", pname, f"{pname} 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
+        for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT", "A1DOT2"):
+            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

@@ -114,22 +114,60 @@ def validate(self):
         super().validate()
         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
+        for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT", "A1DOT2"):
+            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/binary_ell1.py

@@ -218,12 +218,16 @@ def validate(self):
         """Validate parameters."""
         super().validate()
 
-        if self.TASC.value is None:
+        if self._bp("TASC").value is None:
             raise MissingParameter("ELL1", "TASC", "TASC is required for ELL1 model.")
-        for p in ["EPS1", "EPS2"]:
-            pm = getattr(self, p)
+        for p in ("EPS1", "EPS2"):
+            pm = self._bp(p)
             if pm.value is None:
                 pm.value = 0
+        for p in ("A1DOT", "A1DOT2"):
+            if self._hasbp(p) and self._bp(p).value is None:
+                self._bp(p).value = 0.0
+                self._bp(p).frozen = True
 
     def change_binary_epoch(self, new_epoch):
         """Change the epoch for this binary model.
@@ -307,6 +311,34 @@ def change_binary_epoch(self, new_epoch):
         return dt_integer_orbits
 
 
+class BinaryELL12(BinaryELL1):
+    """Outer-orbit ELL1 model for a hierarchical triple system.
+
+    This is identical to :class:`pint.models.binary_ell1.BinaryELL1` except that
+    all of its parameters carry a ``_2`` suffix (``PB_2``, ``A1_2``,
+    ``TASC_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_ell1.BinaryELL12
+    """
+
+    register = True
+    category = "pulsar_system_outer"
+    param_suffix = "_2"
+    binary_param_tag = "BINARY2"
+
+    def __init__(self):
+        super().__init__()
+        self._apply_param_suffix()
+
+
 class BinaryELL1H(BinaryELL1):
     """ELL1 modified to use H3 parameter for Shapiro delay.
 

src/pint/models/chromatic_model.py

@@ -7,9 +7,9 @@
 from loguru import logger as log
 
 from pint import DMconst
-from pint.exceptions import MissingParameter
+from pint.exceptions import MissingParameter, MissingTOAs
 from pint.models.parameter import MJDParameter, floatParameter, prefixParameter
-from pint.models.timing_model import DelayComponent, MissingParameter, MissingTOAs
+from pint.models.timing_model import DelayComponent
 from pint.toa_select import TOASelect
 from pint.utils import split_prefixed_name, taylor_horner, taylor_horner_deriv
 

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,