Better modelling of hierarchical triple systems - #9
Conversation
…l to account for Doppler shift due to outer orbits
There was a problem hiding this comment.
Pull request overview
This PR updates the TEMPO2 T2 binary model to better support hierarchical multi-orbit (e.g., triple) systems by computing inner-orbit phases at a retarded time that accounts for outer-orbit Roemer delays, and by extending derivative handling (including adding X2DOT support) so fitting remains consistent under the new timing relationship.
Changes:
- Reorders/updates T2 orbit evaluation so inner orbits use barycentric times retarded by accumulated outer-orbit delays, and applies a chain-rule correction factor (
fac) to derivatives. - Adds X2DOT (second derivative of projected semi-major axis) support to T2model parameter propagation and derivative reporting.
- Extends parfile parsing and parameter initialization/allocation to support per-companion XDOT_n, X2DOT_n, and PBDOT_n.
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 3 comments.
| File | Description |
|---|---|
| T2model.C | Implements outer-to-inner delay evaluation, chain-rule derivative scaling (fac), and adds X2DOT propagation/derivatives. |
| readParfile.C | Adds parsing for per-companion PBDOT_n, XDOT_n, and X2DOT_n parameters. |
| initialise.C | Adds per-companion labels for XDOT_n/X2DOT_n/PBDOT_n and allocates companion-array storage for XDOT/X2DOT/PBDOT. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| else if (strcasecmp(str,"PBDOT")!=0 && (str[0]=='P' || str[0]=='p') && (str[1]=='B' || str[1] == 'b') && (str[2]=='D' || str[2]=='d')) | ||
| { | ||
| int val; | ||
| if (sscanf(str+6,"%d",&val)==1) | ||
| { | ||
| if (val-1<psr->param[param_pbdot].aSize) | ||
| readValue(psr,str,fin,&(psr->param[param_pbdot]),val-1); | ||
| } | ||
| } |
There was a problem hiding this comment.
@coclar - I think that we probably don't care about this, i.e. we can assume that the orbits are in the tempo2 units, rather than worrying about tempo backwards compatibility for a feature that I think isn't in tempo. Please can you confirm.
| else | ||
| { | ||
| com1 = arr; | ||
| com1 = 0; | ||
| com2 = arr+1; | ||
| } |
There was a problem hiding this comment.
@coclar - to be honest I suspect this is just a limitation of the code, but can you check this is the intended fitting behaviour before we merge.
There was a problem hiding this comment.
This was a mistake on my part - for a while I had the orbit order going in the other direction (e.g. PB_2 was for the inner orbit, rather than outer), but decided to switch them around. I forgot to fix these for-loop limits when switching the order. Copilot is correct! I have fixed this now.
| else if (strcasecmp(str,"XDOT")!=0 && (str[0]=='X' || str[0]=='x') && | ||
| (str[1]=='D' || str[1]=='d')) | ||
| { | ||
| int val; | ||
| if (sscanf(str+5,"%d",&val)==1) | ||
| { | ||
| if (val-1<psr->param[param_a1dot].aSize) | ||
| readValue(psr,str,fin,&(psr->param[param_a1dot]),val-1); | ||
| } | ||
| } |
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:
Since a change in$d_{\rm out}$ is equivalent to a change in $T_0$ for the inner orbit, we have:
so the total derivative is:
The new
facterm 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.