Better modelling of hierarchical triple systems

[coclar]

While writing this paper, we found that better modelling of multiple orbits was necessary in the T2 binary model.

This system consists of an 8-day MSP-WD binary system, which is also in a 70-year eccentric orbit with a third body. The outer orbit is large enough to cause a significant time-varying Doppler shift of the inner orbit.

To account for this, we first need to compute the Roemer delay of the outer orbit, $d_{\rm out}$, subtract this from the barycentric arrival time, $t_{\rm ssb}$, and compute the inner orbital phase at the retarded time ($t_{\rm ssb} - d_{\rm out}$). The T2 model was previously computing the inner orbital phase at $t_{\rm ssb}$ instead. For (approximately) constant $d_{\rm out}$, this effect is absorbed by a shift in $T_0$, but time-dependent changes in $d_{\rm out}$ were not accounted for.

For fitting, we also needed to update the derivatives. Since changes in the outer orbit parameters now also affect the inner-orbital Roemer delay, we need to apply a chain-rule correction. The inner orbital delay is now $d_{\rm in}(t_{\rm ssb} - T_0 - d_{\rm out}$). The derivative of the total delay, $d_{\rm in} + d_{\rm out}$, w.r.t. some outer-orbital parameter $p_{\rm out}$ is:

$$\frac{\partial d_{\rm tot}}{\partial p_{\rm out}} = \frac{\partial d_{\rm out}}{\partial p_{\rm out}} + \frac{\partial d_{\rm in}}{\partial p_{\rm out}}$$

Since a change in $d_{\rm out}$ is equivalent to a change in $T_0$ for the inner orbit, we have:

$$ \frac{\partial d_{\rm in}}{\partial p_{\rm out}} = \frac{\partial d_{\rm in}}{\partial d_{\rm out}} \frac{\partial d_{\rm out}}{\partial p_{\rm out}} = \frac{\partial d_{\rm in}}{\partial T_0} \frac{\partial d_{\rm out}}{\partial p_{\rm out}}$$

so the total derivative is:

$$\frac{\partial d_{\rm tot}}{\partial p_{\rm out}} = \frac{\partial d_{\rm out}}{\partial p_{\rm out}}\left(1 + \frac{\partial d_{\rm in}}{\partial T_0}\right)$$

The new fac term in T2model.C implements the term in brackets here.

For all this to work, it is necessary that the inner-binary parameters have no suffix (e.g., PB, A1, etc.), while outer-orbital parameters need suffixes that increase while moving from inner- to outer-orbits. (e.g. PB_2, PB_3, etc.).

Finally, for this pulsar we also had to add an X2DOT parameter to the T2model.C (second derivative of the semi-major axis). This parameter already existed for the MSS binary model, but now is also used in the T2 model.

With these changes, the T2 model is consistent with an equivalent model added to PINT, verified by @vhaasteren.