# Precision Predictions for Charged Higgs Boson Decays in the Real and Complex
NMSSM
Talk presented at the International Workshop on Future Linear Colliders
(LCWS2021), 15-18 March 2021, Online Conference, C21-03-15.1.

Thi Nhung Dao${}^{1\,}$111E-mail<EMAIL_ADDRESS>Margarete
Mühlleitner${}^{2\,}$222E-mail<EMAIL_ADDRESS>Shruti
Patel${}^{2,3\,}$333E-mail<EMAIL_ADDRESS>Kodai
Sakurai${}^{2\,}$444Speaker 555Address after April 2021, Department of
Physics, Tohoku University, Sendai, Miyagi 980-8578, Japan; E-mail:
<EMAIL_ADDRESS>
1Institute For Interdisciplinary Research in Science and Education, ICISE,
590000 Quy Nhon, Vietnam.
2Institute for Theoretical Physics, Karlsruhe Institute of Technology,
Wolfgang-Gaede-Str. 1, 76131 Karlsruhe, Germany.
3Institute for Nuclear Physics, Karlsruhe Institute of Technology, 76344
Karlsruhe, Germany.

###### Abstract

We present the full next-to-leading order (NLO) supersymmetric (SUSY)
electroweak and SUSY-QCD corrections to the decay widths of the charged Higgs
boson decays into on-shell final states in the framework of the CP-conserving
and CP-violating Next-to-Minimal Supersymmetric Model (NMSSM) of Ref.
Dao:2020dfb . The newly calculated corrections have been implemented in the
code NMSSMCALCEW. In these proceedings, we discuss the impact of the NLO
corrections on the charged Higgs boson branching ratios in a wide range of the
parameter space that is still compatible with the experimental constraints. We
also investigate the effect of CP violation in these corrections.

##  I Introduction

The Standard Model (SM) has been established as a low-energy effective theory
that describes electroweak symmetry breaking (EWSB). However, there exist
several unsolved phenomena, such as neutrino oscillations, dark matter, baryon
asymmetry of the Universe, and inflation. Since they are not explained within
the framework of the SM, new physics is needed to address these puzzles. One
among the models beyond the SM that solve the hierarchy problem is given by
supersymmetry (SUSY) Golfand:1971iw ; Volkov:1973ix ; Wess:1974tw ;
Fayet:1974pd ; Fayet:1977yc ; Fayet:1976cr ; Nilles:1983ge ; Haber:1984rc ;
Sohnius:1985qm ; Gunion:1984yn ; Gunion:1986nh . Its simplest version is
realized by the Minimal Supersymmetric SM (MSSM), where two Higgs doublet
fields are required in order to ensure supersymmetry and to keep the theory
anomaly-free. In this work, we focus on the Next-to minimal Supersymmetric SM
(NMSSM) Fayet:1974pd ; Barbieri:1982eh ; Dine:1981rt ; Nilles:1982dy ;
Frere:1983ag ; Derendinger:1983bz ; Ellis:1988er ; Drees:1988fc ;
Ellwanger:1993xa ; Ellwanger:1995ru ; Ellwanger:1996gw ; Elliott:1994ht ;
King:1995vk ; Franke:1995tc ; Maniatis:2009re ; Ellwanger:2009dp . Besides the
virtue of solving the so-called $\mu$ problem Kim:1983dt , the NMSSM disposes
of an interesting Higgs boson phenomenology. Its Higgs sector contains a
complex singlet field in addition to the two Higgs doublet fields which yields
three CP-even Higgs bosons, two CP-odd Higgs bosons and two charged Higgs
bosons after EWSB, implying interesting Higgs physics full of variety.
Furthermore, the Higgs sector can be CP-violating already at the tree level
unlike the MSSM.

So far the LHC experiments have not found any direct evidence of new particles
beyond the SM, leading to lower bounds on the masses of new particles. In
addition, measurements of the observables for the discovered Higgs boson (e.g.
production cross sections, decay branching ratios and couplings) are
consistent with the predictions of the SM within the current experimental
uncertainties, meaning that the discovered Higgs boson behaves very SM-like.
While the SUSY particles are required to be heavy enough to comply with the
constraints from the direct searches they may still give significant
contributions in the higher-order corrections to Higgs observables. In
addition, the properties of the discovered Higgs boson will be measured very
precisely in future collider experiments, such as the high-luminosity (HL-)
LHC ApollinariG.:2017ojx , the International Linear Collider (ILC)
Baer:2013cma , the Future Circular Collider (FCC-ee) Gomez-Ceballos:2013zzn ,
or the Circular Electron Positron Collider (CEPC) CEPC-SPPCStudyGroup:2015csa
. Hence, precise theoretical calculations are inevitable in order to be able
to identify indirect signs of new physics. In the context of the NMSSM,
higher-order corrections have been calculated to the NMSSM Higgs boson masses
up to two-loop accuracy both in the CP-conserving and CP-violating NMSSM. For
a recent overview, see Slavich:2020zjv . For the Higgs boson decays, full one-
loop corrections to neutral Higgs boson decays into two gauge bosons and two
fermions have been calculated in Ref. Domingo:2018uim  in the CP-violating
NMSSM. The effects of Sudakov logarithms on fermionic decays of heavy Higgs
bosons have been studied in Ref. Domingo:2019vit . For Higgs-to-Higgs decays,
not only the one-loop corrections Nhung:2013lpa ; Belanger:2017rgu  but also
the two loop corrections with $\mathcal{O}(\alpha_{S}\alpha_{t})$
Muhlleitner:2015dua  have been studied. Furthermore, in Ref. Baglio:2019nlc ,
all on-shell two-body decays of neutral Higgs bosons have been evaluated at
NLO electroweak (EW) and NLO SUSY-QCD order in the NMSSM with CP violation. In
Ref. Dao:2019nxi , the NLO EW corrections to charged Higgs boson decays into a
$W$ boson and a neutral Higgs boson have been calculated and the gauge
dependence, which arises from higher-order corrections to the mass and due to
the resummed $Z$ factor of the external Higgs boson has been studied. This
problem and possible strategies to preserve or restore gauge invariance has
been further investigated in Domingo:2020wiy . In Ref. Dao:2020dfb , we
completed the computation of the full NLO corrections to the on-shell two-body
decays of the charged Higgs bosons in the NMSSM with CP violation. We
evaluated the NLO SUSY-EW and SUSY-QCD corrections to the decays of the
charged Higgs bosons into quarks, leptons, electroweakinos, squarks and
sleptons for the first time 666In the previous implementation in the program
NMSSMCALC Baglio:2013iia , the charged Higgs decay widths included the state-
of-the art QCD corrections to the decays into quarks and squarks as well as
the resummed SUSY-QCD corrections and SUSY-EW corrections for the decays into
quarks and leptons.. These new NLO corrections have been implemented in the
new version of the program NMSSMCALCEW Baglio:2019nlc  777The program is
available at the URL: https://www.itp.kit.edu/~maggie/NMSSMCALCEW/. In these
proceedings, we summarize our results obtained in Ref. Dao:2020dfb , and
discuss the sizes of the pure NLO corrections for the charged Higgs boson
decays into quarks and electroweakinos. We furthermore show numerical results
for the effects of CP violation on the loop-corrections to the decays into
quarks.

##  II The Lagrangian of the NMSSM

We briefly describe the Higgs potential of the NMSSM and the electroweakino
sector in order to fix our conventions and input parameters. We work in the
framework of the scale-invariant NMSSM with $Z_{3}$ symmetry. The tree-level
Higgs potential obtained from the $F$\- and $D$-terms of the supersymmetric
Lagrangian and the soft SUSY-breaking Lagrangian is given by

$\displaystyle V_{H}$ $\displaystyle=|\lambda
S|^{2}\left(H_{u}^{\dagger}H_{u}+H_{d}^{\dagger}H_{d}\right)+\left|-\epsilon_{ij}\lambda\left(H_{d}^{i}H_{u}^{j}\right)+\kappa
S^{2}\right|^{2}$
$\displaystyle+\frac{1}{2}g_{2}^{2}\left|H_{u}^{\dagger}H_{d}\right|^{2}+\frac{1}{8}(g_{1}^{2}+g_{2}^{2})\left(H_{u}^{\dagger}H_{u}-H_{d}^{\dagger}H_{d}\right)^{2}$
$\displaystyle+m_{H_{u}}^{2}H_{u}^{\dagger}H_{u}+m_{H_{d}}^{2}H_{d}^{\dagger}H_{d}+m_{S}^{2}\left|S\right|^{2}+\left[-\epsilon_{ij}\lambda
A_{\lambda}\left(H_{d}^{i}H_{u}^{j}\right)S+\frac{1}{3}\kappa
A_{\kappa}S^{3}+{\rm h.c.}\right],$ (1)

where $\lambda$ and $\kappa$ are complex dimensionless parameters defined in
the superpotential, and the corresponding terms proportional to $A_{\lambda}$
and $A_{\kappa}$ arise from the soft SUSY-breaking terms. The gauge couplings
of the $U(1)_{Y}$ and $SU(2)_{L}$ symmetry are denoted by $g_{1}$ and $g_{2}$,
respectively. The Higgs doublet fields $H_{u}$ and $H_{d}$ and the complex
singlet field $S$ are expressed in terms of the component fields and vacuum
expectation values (VEVs) as

$\displaystyle H_{u}=e^{i\varphi_{u}}\begin{pmatrix}h_{u}^{+}\\\
\frac{1}{\sqrt{2}}(v_{u}+h_{u}+ia_{u})\end{pmatrix},\
H_{d}=\begin{pmatrix}\frac{1}{\sqrt{2}}(v_{d}+h_{d}+ia_{d})\\\
h_{d}^{-}\end{pmatrix},\
S=\frac{1}{\sqrt{2}}e^{i\varphi_{s}}(v_{s}+h_{s}+ia_{s})\;,$ (2)

where $v_{u},v_{d}$ and $v_{s}$ are the VEVs of $H_{u}$, $H_{d}$, and $S$,
respectively. The EW VEV is given by
$v^{2}=v_{u}^{2}+v_{d}^{2}\simeq(246~{}{\rm GeV})^{2}$. The two CP-violating
phases $\varphi_{u}$ and $\varphi_{s}$ denote the phase differences between
the VEVs. The neutral component fields are transformed into the mass
eigenstates through the orthogonal rotation matrix $\mathcal{R}$,

$\displaystyle(h_{1},h_{2},h_{3},h_{4},h_{5},G^{0})^{T}=\mathcal{R}(h_{d},h_{u},h_{s},a_{d},a_{u},a_{s})^{T},$
(3)

where $G^{0}$ is the neutral Nambu-Goldstone (NG) boson. Similarly, the mass
eigenstates for the charged Higgs bosons are obtained by,

$\displaystyle\begin{pmatrix}G^{\pm}\\\
H^{\pm}\end{pmatrix}=\begin{pmatrix}-c_{\beta}&s_{\beta}\\\
s_{\beta}&c_{\beta}\end{pmatrix}\begin{pmatrix}h_{d}^{\pm}\\\
h_{u}^{\pm}\end{pmatrix},$ (4)

where the mixing angle $\beta$ is defined by $\tan\beta=v_{u}/v_{d}$.

The tree-level Higgs sector of the CP-violating NMSSM is described by eighteen
independent input parameters which we choose as

$\displaystyle
m_{H_{d}}^{2},m_{H_{u}}^{2},m_{S}^{2},M_{W}^{2},M_{Z}^{2},e,\tan\beta,v_{s},\varphi_{s},\varphi_{u},\lvert\lambda\rvert,\varphi_{\lambda},\lvert\kappa\rvert,\varphi_{\kappa},{\rm
Re}A_{\lambda},{{\rm Im}A_{\lambda}},{\rm Re}A_{\kappa},{{\rm
Im}A_{\kappa}}~{}.$ (5)

Here the three Lagrangian parameters $g_{1}$, $g_{2}$ and $v$ have been
replaced by the three physical observables $M_{W}$, $M_{Z}$ and the electric
coupling $e$. The effective $\mu$ parameter can be expressed as

$\displaystyle\mu_{\rm eff}=\lvert\mu_{\rm eff}\rvert
e^{i\varphi_{\mu}}=\frac{|\lambda|v_{s}}{\sqrt{2}}e^{i(\varphi_{\lambda}+\varphi_{s})}\;.$
(6)

The real part of $A_{\lambda}$, ${\rm Re}A_{\lambda}$, can be replaced by the
charged Higgs boson mass $M^{2}_{H^{\pm}}$ through

$\displaystyle M^{2}_{H^{\pm}}=M_{W}^{2}+\frac{\lambda
v_{S}}{\sin(2\beta)}\left(\sqrt{2}{\rm Re}A_{\lambda}+\kappa
v_{S}\right)-\frac{\lambda^{2}v^{2}}{2}~{}.$ (7)

The soft SUSY-breaking parameters $m_{H_{d}}^{2},m_{H_{u}}^{2},m_{S}^{2},{\rm
Im}A_{\lambda},{\rm Im}A_{\kappa}$ are fixed by the tadpole conditions for the
five neutral Higgs bosons.

Adding the superchiral singlet field $\hat{S}$ introduces an additional degree
of freedom, the singlino $\tilde{S}$, in the electroweakino sector. This
results in five neutralino mass eigenstates, which are related to the gauge
eigenstates via the 5$\times$5 unitary matrix $N$,

$\displaystyle\chi^{0}=N\psi^{0}\;,$ (8)

with
$\psi^{0}=(\tilde{B},~{}\tilde{W}^{3},~{}\tilde{H}_{d},~{}\tilde{H}_{u},~{}\tilde{S})^{T}$,
where $\tilde{B},~{}\tilde{W}^{3}$ stand for neutral gaugino states,
$\tilde{S}$ for the singlino state and $~{}\tilde{H}_{d}$, and $\tilde{H}_{u}$
are the Higgsino states. The chargino mass eigenstates in the basis of the
Weyl spinors for the charginos,
$\tilde{\chi}^{+}_{L}=(\tilde{\chi}^{+}_{L_{1}},\tilde{\chi}^{+}_{L_{2}})^{T},\
\tilde{\chi}^{-}_{R}=(\tilde{\chi}^{-}_{R_{1}},\tilde{\chi}^{-}_{R_{2}})^{T}$,
are obtained by rotating the spinors in gauge basis,
$\psi^{-}_{R}=(\tilde{W}^{-},\tilde{H}^{-}_{d})^{T},\
\psi^{+}_{L}=(\tilde{W}^{+},\tilde{H}^{+}_{u})^{T}$, as

$\displaystyle\tilde{\chi}_{L}^{+}=V\psi_{L}^{+},\ \ \
\tilde{\chi}_{R}^{-}=U\psi_{R}^{-}\;,$ (9)

with the 2$\times$2 unitary matrices $U$, $V$. The specific expressions for
the mass matrices of the neutralinos and charginos in terms of the input
parameters can be found in Eqs. (17) and (21) of Ref. Dao:2020dfb . By
definition, they are diagonalized by the unitary matrices $N$, $U$ and $V$.

##  III The Charged Higgs Boson Decay Widths Including Higher-Order
Corrections

In Ref. Dao:2020dfb , we evaluated the NLO SUSY-EW corrections and SUSY-QCD
corrections to the two-body on-shell decay widths of the charged Higgs bosons,
i.e., $H^{+}\to t\bar{b}$, $H^{+}\to\nu\bar{\tau}$,
$H^{+}\to\chi_{i}^{+}\chi_{j}^{0}$ ($i=1,2,~{}j=1,...5$),
$H^{+}\to\tilde{t}\tilde{b}$ and $H^{+}\to\tilde{\tau}\tilde{\nu}$. In order
to get UV finite results, we used a mixed OS and $\overline{\rm DR}$
renormalization scheme for the Higgs sector and the OS scheme for the gauge
sector and the SM fermion sector. For the renormalization of the
electroweakino, the squark and the slepton sectors, both the OS scheme and the
$\overline{\rm DR}$ scheme were utilized. In this section, we give schematic
formulae for the partial decay widths into $t\bar{b}$ and electroweakinos,
including the higher-order corrections.

For the decay $H^{+}\to t\bar{b}$, the loop-corrected partial decay width can
schematically be written as

$\displaystyle\Gamma(H^{+}\to t\bar{b})^{\rm NLO}=\Gamma^{\rm
LO~{}imp.}_{H^{+}\to t\bar{b}}+\Gamma^{\rm SUSYQCD}_{H^{+}\to
t\bar{b}}+\Gamma^{\rm SUSYEW}_{H^{+}\to t\bar{b}}.$ (10)

The improved leading-order (LO) contribution $\Gamma^{\rm LO~{}imp.}_{H^{+}\to
t\bar{b}}$, which already includes the state-of-art QCD corrections and the
resummed SUSY-QCD and SUSY-EW corrections, is given in terms of
$\mu_{t}=m_{t}^{2}/M_{H^{\pm}}^{2}$ and $\mu_{b}=m_{b}^{2}/M_{H^{\pm}}^{2}$ by

$\displaystyle\Gamma^{\rm LO~{}imp.}_{H^{+}\rightarrow t\bar{b}}$
$\displaystyle=$
$\displaystyle\frac{3G_{F}M^{3}_{H^{\pm}}}{4\sqrt{2}\pi}|V_{tb}|^{2}\,\beta^{1/2}\left(\mu_{t},\mu_{b}\right)\,\left\\{(1-\mu_{t}-\mu_{b})\left[\mu_{t}\frac{1}{\tan^{2}\beta}\left(1+\frac{4}{3}\frac{\alpha_{s}}{\pi}\delta_{tb}^{+}\right)\right.\right.$
(11) $\displaystyle\left.\left.+\mu_{b}\tan^{2}\beta
R^{2}\left(1+\frac{4}{3}\frac{\alpha_{s}}{\pi}\delta_{bt}^{+}\right)\right]-4\mu_{t}\mu_{b}R\left(1+\frac{4}{3}\frac{\alpha_{s}}{\pi}\delta_{tb}^{-}\right)\right\\}\;,$

with

$\beta^{1/2}(x,y)=(1-x-y)^{2}-4xy,$ (12)

and where $G_{F}$ denotes the Fermi constant, $\alpha_{s}$ the strong coupling
constant, and $V_{tb}$ the top-bottom CKM matrix element. Explicit expressions
for the QCD correction factors $\delta_{tb}^{+},\delta_{bt}^{+}$ and
$\delta_{tb}^{-}$ can be found in Ref. Spira:2016ztx . The universal SUSY-QCD
and SUSY-EW corrections that are enhanced in the large $\tan\beta$ regime are
resummed into effective bottom Yukawa couplings. These $\Delta_{b}$
corrections Hall:1993gn ; Hempfling:1993kv ; Carena:1994bv ; Pierce:1996zz ;
Carena:1998gk ; Carena:1999py ; Carena:2002bb ; Guasch:2003cv ; Spira:2016ztx
are included in the decay width through the $R$ factor,
$R=\frac{1}{1+\Delta_{b}}\left[1-\frac{\Delta_{b}}{\tan^{2}\beta}\right]$. The
explicit formula for $\Delta_{b}$ is given in Ref. Baglio:2013iia . The
contributions $\Gamma^{\rm SUSYQCD}_{H^{+}\to t\bar{b}}$ and $\Gamma^{\rm
SUSYEW}_{H^{+}\to t\bar{b}}$ correspond to one-loop SUSY-QCD and SUSY-EW
corrections, respectively, where subtraction terms are added to avoid double
counting arising from the $\Delta_{b}$ corrections. The SUSY-EW corrections
contain, besides the pure vertex correction, the external leg corrections from
the mixing between the charged Higgs boson and the charged NG boson and the
$W$ boson, respectively, and the real photon emission $\Gamma_{H^{+}\to
t\bar{b}\gamma}$. Specific expressions for these NLO corrections are given in
Ref. Dao:2020dfb .

The partial decay widths for the decays into electroweakinos are written in
the same notation as

$\displaystyle\Gamma(H^{+}\to\chi_{i}^{+}\chi_{j}^{0})^{\rm NLO}=\Gamma^{\rm
LO}_{H^{+}\to\chi_{i}^{+}\chi_{j}^{0}}+\Gamma^{\rm
SUSYEW}_{H^{+}\to\chi_{i}^{+}\chi_{j}^{0}},$ (13)

$(i=1,2,~{}j=1,...5)$. The LO decay width is given by

$\displaystyle\Gamma^{\rm
LO}_{H^{+}\to\tilde{\chi}_{i}^{0}\tilde{\chi}^{+}_{j}}$ $\displaystyle=$
$\displaystyle\frac{\beta^{1/2}\left(\mu_{i},\mu_{j}\right)M_{H^{\pm}}}{16\pi}\Bigg{[}\left(1-\mu_{i}-\mu_{j}\right)\left(\lvert
g_{H^{+}\tilde{\chi}_{i}^{0}\tilde{\chi}^{-}_{j}}^{L,{\rm
tree}}\rvert^{2}+\lvert{g}_{H^{+}\tilde{\chi}_{i}^{0}\tilde{\chi}^{-}_{j}}^{R,{\rm
tree}}\rvert^{2}\right)$
$\displaystyle-2\sqrt{\mu_{i}\mu_{j}}\left({g}_{H^{+}\tilde{\chi}_{i}^{0}\tilde{\chi}^{-}_{j}}^{L,{\rm
tree}}(g_{H^{+}\tilde{\chi}_{i}^{0}\tilde{\chi}^{-}_{j}}^{R,{\rm
tree}})^{*}+{g}_{H^{+}\tilde{\chi}_{i}^{0}\tilde{\chi}^{-}_{j}}^{R,{\rm
tree}}(g_{H^{+}\tilde{\chi}_{i}^{0}\tilde{\chi}^{-}_{j}}^{L,{\rm
tree}})^{*}\right)\Bigg{]}\,,$

with $\mu_{i}=M_{\tilde{\chi}_{i}^{0}}^{2}/M_{H^{\pm}}^{2}$ and
$\mu_{j}=M_{\tilde{\chi}_{j}^{+}}^{2}/M_{H^{\pm}}^{2}$, where the tree-level
couplings are given by

$\displaystyle g_{H^{+}\tilde{\chi}_{i}^{0}\tilde{\chi}^{-}_{j}}^{R,{\rm
tree}}$
$\displaystyle=-\frac{1}{\sqrt{2}}{V}_{j2}\Big{(}(g_{1}N_{i1}+g_{2}N_{i2})c_{\beta}e^{i\varphi_{u}}+\sqrt{2}\lambda^{\ast}N_{i5}s_{\beta}\Big{)}-g_{2}e^{i\varphi_{u}}{V}_{j1}N_{i4}c_{\beta}\,,$
(15) $\displaystyle g_{H^{+}\tilde{\chi}_{i}^{0}\tilde{\chi}^{-}_{j}}^{L,{\rm
tree}}$
$\displaystyle=-\frac{1}{\sqrt{2}}{U}_{j2}^{\ast}\Big{(}{-}g_{1}N_{i1}^{\ast}s_{\beta}{-}g_{2}N_{i2}^{\ast}s_{\beta}+\sqrt{2}e^{i\varphi_{u}}\lambda
N_{i5}^{\ast}c_{\beta}\Big{)}-g_{2}{U}_{j1}^{\ast}N_{i3}^{\ast}s_{\beta}.$
(16)

There are no NLO SUSY-QCD corrections, and the NLO SUSY-EW corrections
comprise besides the vertex corrections the mixing between $H^{+}$ and the
charged Goldstone and $W$ bosons, respectively, and the real photon emission
$\Gamma_{H^{+}\to\chi_{i}^{+}\chi_{j}^{0}\gamma}$. The concrete formulae are
also given in Ref. Dao:2020dfb .

##  IV Numerical results

###  IV.1 Branching ratios for $H^{+}$ at LO

Figure 1: The branching ratios of the charged Higgs boson $H^{+}$ decays
$H^{+}\to\tau^{+}\nu_{\tau}$ (pink), $H^{+}\to t\bar{b}$ (gray), $H^{+}\to WH$
(orange), $H^{+}\to\tilde{\chi}^{+}_{1}\tilde{\chi}^{0}$ (cyan),
$H^{+}\to\tilde{\chi}^{+}_{2}\tilde{\chi}^{0}$ (violet),
$H^{+}\to\tilde{t}_{1}\tilde{b}_{1}^{\ast}$ (black),
$H^{+}\to\tilde{\tau}\tilde{\nu_{\tau}}^{\ast}$ (blue), at LO as a function of
the charged Higgs boson mass $M_{H^{\pm}}$. .

We show the branching ratios of the charged Higgs bosons at LO in Fig. 1 to
illustrate the behavior as a function of the charged Higgs boson mass. In the
plot, pure NLO corrections are switched off for all decay processes. The input
parameters for the Higgs sector and the electroweakino sector are fixed as

$\displaystyle\begin{array}[]{llllllllllll}\tan\beta&=&10.14\;,&\;|\lambda|&=&0.093\;,&\;|\kappa|&=&-0.0821\;,\\\
M_{H^{+}}&=&1500\ {\rm GeV},&\mu_{\rm eff}&=&-891\ {\rm GeV},&\;{\rm
Re}A_{\kappa}&=&-1.6\ {\rm TeV},\\\ |M_{1}|&=&752\ {\rm GeV},&\;|M_{2}|&=&806\
{\rm GeV}.\end{array}$ (20)

The soft SUSY-breaking parameters for the sfermions and the gluino, and the
complex phases are fixed as

$\displaystyle|M_{3}|=2334\ {\rm GeV},\;|A_{t}|=3.7\ {\rm TeV},\;|A_{b}|=2\
{\rm TeV},\;|A_{\tau}|=2\ {\rm TeV},$ $\displaystyle m_{\tilde{Q_{3}}}=1.43\
{\rm TeV},\;m_{\tilde{t}_{R}}=1.93\ {\rm TeV},\;m_{\tilde{b}_{R}}=2.16\ {\rm
TeV},\;m_{\tilde{L}_{3}}=1.22\ {\rm TeV},\;$ $\displaystyle
m_{\tilde{\tau}_{R}}=1.21\ {\rm
TeV},\;m_{\tilde{u}_{R},\tilde{c}_{R}}=m_{\tilde{d}_{R},\tilde{s}_{R}}=m_{\tilde{Q}_{1,2}}=m_{\tilde{L}_{1,2}}=m_{\tilde{e}_{R},\tilde{\mu}_{R}}=3\
{\rm TeV},$
$\displaystyle\varphi_{M_{1},M_{2},M_{3}}=\varphi_{A_{t},A_{b},A_{\tau}}=\varphi_{\mu}=\varphi_{\kappa}=0\;.$
(21)

In this plot, the decays into the neutral Higgs bosons plus the $W$ boson and
the electroweakinos, respectively, are summed with respect to the final states
as

$\displaystyle\mbox{BR}(H^{+}\to W^{+}H)$
$\displaystyle\equiv\sum_{i=1}^{3}\mbox{BR}(H^{+}\to
W^{+}H_{i})+\sum_{j=1}^{2}\mbox{BR}(H^{+}\to W^{+}A_{j})\;,$
$\displaystyle\mbox{BR}(H^{+}\to\tilde{\chi}^{+}_{1,2}\tilde{\chi}^{0})$
$\displaystyle\equiv\sum_{i=1}^{5}\mbox{BR}(H^{+}\to\tilde{\chi}^{+}_{1,2}\tilde{\chi}^{0}_{i})\;.$
(22)

Similarly, the decays into the slepton final states are summed up. The mass
spectrum for this bench mark point is given by

$\displaystyle M_{H_{1}}=123.33\mbox{ GeV},\;M_{H_{2}}=1.09\mbox{
TeV},\;M_{H_{3}}=M_{A_{1}}=1.50\mbox{ TeV},\;$ $\displaystyle
M_{A_{2}}=1.93\mbox{ TeV}\;,M_{\chi^{+}_{1}}=817.87~{}\mbox{
GeV},\;M_{\chi^{+}_{2}}=927.87~{}\mbox{ GeV},\;$ $\displaystyle
M_{\chi^{0}_{1}}=747.69~{}{\rm GeV},\;M_{\chi^{0}_{2}}=820.07~{}{\rm
GeV},\;M_{\chi^{0}_{3}}=894.20~{}{\rm GeV},\;$ $\displaystyle
M_{\chi^{0}_{4}}=929.24~{}{\rm GeV},\;M_{\chi^{0}_{5}}=1.56~{}{\rm
TeV},\;M_{\tilde{t}_{1}}=1.38~{}{\rm TeV},\;$
$\displaystyle\;M_{\tilde{b}_{1}}=1.43~{}{\rm
TeV},M_{\tilde{\tau}_{1}}=M_{\tilde{\tau}_{2}}=M_{\tilde{\nu}_{\tau}}=1.2~{}{\rm
TeV}\;.$ (23)

As can be inferred from the plot, for $M_{H^{+}}\lesssim 1.9\ {\rm TeV}$, the
decay into the top-bottom final state is the main decay mode because of the
$\tan\beta$ enhancement for the bottom Yukawa coupling. Most of the decay
channels into the electroweakinos are kinematically open for $1.9\lesssim
M_{H^{+}}\ {\rm TeV}$, so that the branching ratios into electoroweakinos
become more important than those into top-bottom.

###  IV.2 Impact of the NLO Corrections on the Charged Higgs Branching Ratios

Figure 2: Correlation between the relative change for the decay branching
ratio, $\Delta_{\rm BR}$, and the NLO branching ratio for $H^{+}\to t\bar{b}$
(left panel) and $H^{+}\to\tilde{W}^{+}\tilde{B}$ (right panel). The color bar
shows the size of the relative NLO corrections for the decay width
$\delta_{\Gamma}$.

In this section, we discuss the impact of the NLO corrections on the branching
ratios of the charged Higgs boson for the decays into the top and bottom quark
pair and the electroweakino final states. In order to describe the size of the
NLO corrections for the branching ratios as well as the partial widths of the
decays $H^{+}\to XY$, we introduce the following quantities,

$\displaystyle\Delta_{\rm BR}(H^{+}XY)$ $\displaystyle=\frac{{\rm BR}^{\rm
NLO}(H^{+}\to XY)-{\rm BR}^{\rm LO}(H^{+}\to XY)}{{\rm max}({\rm BR}^{\rm
NLO}(H^{+}\to XY),{\rm BR}^{\rm LO}(H^{+}\to XY))}\;,$ (24)
$\displaystyle\delta_{\Gamma}(H^{+}XY)$ $\displaystyle=\frac{\Gamma(H^{+}\to
XY)^{\rm NLO}}{\Gamma(H^{+}\to XY)^{{\rm LO}}}-1\;.$ (25)

The relative change $\Delta_{\rm BR}$ is normalized to the maximum of the
branching ratio at NLO and LO, respectively, to avoid an enhancement due to
the smallness of ${\rm BR}^{\rm LO}$. The notation ’LO’ basically means that
we refer to the old implementation in NMSSMCALC without the genuine SUSY-EW
and SUSY-QCD vertex corrections. This means that the LO quantities are
calculated with ’Higgs effective tree-level couplings’, meaning that the Higgs
tree-level rotation matrix elements have been replaced by the loop-corrected
ones. Furthermore, these ’LO’ quantities also include the QCD corrections and
resummed SUSY-EW and SUSY-QCD corrections in effective quark couplings as
already implemented in the first release of NMSSMCALC Baglio:2013iia . Hence,
the relative changes $\Delta_{\rm BR}$ and $\delta_{\Gamma}$ quantify the
effect of the pure NLO corrections that were newly evaluated in Ref.
Dao:2020dfb . For the following scatter plots, we scanned the NMSSM input
parameters in the ranges given in Table. 1 together with the parameter
settings in Eqs. (128) and (129) of Ref. Dao:2020dfb  to understand the
typical size of the NLO corrections in a wide range of the parameter space. We
only retained points that are compatible with the LHC Higgs data. For details,
see Ref. Dao:2020dfb . For the calculation of the NLO corrections, we used the
OS1 renormalization scheme for the eletroweakino sector, the $\overline{\rm
DR}$ scheme for the stop and sbottom sector, and the OS scheme for the stau
sector. The definition of each scheme can be found in Ref. Baglio:2019nlc

In the left panel of Fig. 2, we show the relative change $\Delta_{\rm BR}$ for
$H^{+}\to t\bar{b}$ as a function of the branching ratio at NLO. The color
code denotes the relative change of the decay width. As written in Eq. (10),
the total NLO corrections for the partial width are determined by the sum of
the SUSY-QCD and the SUSY-EW corrections. The latter gives negative
corrections ranging between $-15\%$ and $-2\%$, and the former remains between
$-15\%$ and $7\%$ so that the relative change $\delta_{\Gamma}$ is negative in
most of the parameter region and the maximum size
$\lvert\delta_{\Gamma}\rvert$ is $29\%$. For this decay channel, the relative
change of the branching ratio $\Delta_{\rm BR}$ is basically correlated with
that of the partial decay width.

In the right panel of Fig. 2, the relative change $\Delta_{\rm BR}$ for the
decay into the charged Wino and the Bino final state,
$H^{+}\to\tilde{W}^{+}\tilde{B}$, is shown as a function of the corresponding
branching ratio at NLO. We select the charged Wino-like chargino and Bino-like
neutralino states from the electroweakino mass eigenstates. They are
identified by the square of the mixing matrix elements, i.e. $|U_{i1}|^{2}$
and $|N_{j1}|^{2}$, which are required to exceed 0.5 for $\tilde{W}^{+}$ and
$\tilde{B}$, respectively. In the plot, we applied cuts for the left tree-
level coupling of Eq. (16) as well as for the mass difference among the
electroweakinos, because we found that artificially enhanced NLO corrections
can appear without such a cut, which are discussed in Sec. 5.5 of Ref.
Dao:2020dfb . As can be inferred from the plot, the relative change of the
branching ratio is $\lvert\Delta_{\rm BR}\rvert<50\%$ for the bulk of the
parameter points. On the other hand, if the size of the branching ratio at NLO
becomes smaller than about $5\times 10^{-2}$, outliers appear and the relative
change $\Delta_{\rm BR}$ can reach almost $-100\%$. These large corrections
arise from large contributions due to the wave function renormalization
constants (WFRCs) for electroweakinos, which are involved in $\Gamma^{\rm
SUSYEW}_{H^{+}\to\chi_{i}^{+}\chi_{j}^{0}}$ in Eq. (13). The WFRCs contains
mixing contributions between another electroweakino state and the one in
question, and they can be significant if the other electroweakino has a large
coupling with the charged Higgs boson.

Figure 3:  CP-violating phase dependence of the partial decay width for
$H^{+}\to t\bar{b}$ at LO (dashed lines) and NLO (solid lines) for the CP-
violating phases of $\mu_{\rm eff}$ (red), $A_{t}$ (blue) and $M_{2}$ (green).
The lower panel displays the dependence of the relative change of the partial
decay width $\delta_{\Gamma}$.

Before we close this section, we discuss how CP-violating phases can affect
the NLO corrections to the decays of the charged Higgs boson. We just examine
the theoretical behavior of the relative change of the partial width when the
phase parameters are varied, and do not take into account experimental
constraints by the measurement of the electric dipole moment (EDM) in this
analysis. In order to illustrate the effect of CP violation, we show the
partial decay width for $H^{+}\to t\bar{b}$ as a function of complex phases
for $\mu_{\rm eff}$, $A_{t}$ and $M_{2}$ in Fig. 3. In the plot, the benchmark
scenario for the CP-conserving case presented in Sec. 5.7 of Ref. Dao:2020dfb
is used, and the three phases are varied individually. Note that we set the
phase $\varphi_{\lambda}$ such that we do not encounter CP violation at tree
level. The nevertheless observed subtle phase dependence on $\varphi_{\mu}$
and $\varphi_{A_{t}}$, respectively, originates from SUSY loop contributions
that are resummed in the $\Delta_{b}$ corrections. Once the pure NLO
corrections are included, we see a stronger dependence on the CP-violating
phases. The relative change in the partial decay width can reach e.g. around
$-10\%$ at $\varphi_{\mu}\simeq-0.46$ while it is around -20% at
$\varphi_{\mu}=0$.

##  V Conclusion

We have summarised the results for the NLO SUSY-QCD and NLO SUSY-EW
corrections for various on-shell two body decays of the charged Higgs boson,
presented in Ref. Dao:2020dfb . The newly evaluated corrections are included
in the latest version of NMSSMCALCEW. The NLO corrections are evaluated by
using a mixed OS and $\overline{\rm DR}$ renormalization scheme in the Higgs
sector. For the results presented here, in the electroweakino sector the OS1
scheme is used while the $\overline{\rm DR}$ (OS) scheme is chosen for the
squarks (slepton) sector. We analysed the relative change of the branching
ratios and the partial decay widths, which are defined in Eqs. (24) and (25),
for $H^{+}\to t\bar{b}$ and $H^{+}\to\tilde{W}^{+}\tilde{B}$ in a wide range
of the parameter space. We found that the NLO corrections to ${\rm
BR}(H^{+}\to t\bar{b})$ are moderate and can reach almost $-30\%$ at ${\rm
BR}^{\rm NLO}\simeq 6\times 10^{-2}$. On the other hand, we found that for
$H^{+}\to\tilde{W}^{+}\tilde{B}$, negative large corrections can appear in a
corner of the parameter space where the size of the branching ratio at NLO
becomes small. We also presented the effects of CP-violating phases on the NLO
corrections to $H^{+}\to t\bar{b}$. We found that the changes of the CP-
violating phases can significantly affect the NLO corrections to the decays of
the charged Higgs boson.

### Acknowledgments

In this work, T.N.D has been funded by the Vietnam National Foundation for
Science and Technology Development (NAFOSTED) under grant number
103.01-2020.17. M.M. and K.S. acknowledge support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762
- TRR 257.

## References

  * (1) T. N. Dao, M. Muhlleitner, S. Patel and K. Sakurai, Eur. Phys. J. C 81 (2021) no.4, 340 doi:10.1140/epjc/s10052-021-09129-1 [arXiv:2012.14889 [hep-ph]].
  * (2) Yu. A. Golfand and E. P. Likhtman, _Extension of the Algebra of Poincare Group Generators and Violation of p Invariance_ , JETP Lett.  13 (1971) 323–326.
  * (3) D. Volkov and V. Akulov, _Is the Neutrino a Goldstone Particle?_ , Phys.Lett.  B46 (1973) 109–110.
  * (4) J. Wess and B. Zumino, _Supergauge Transformations in Four-Dimensions_ , Nucl.Phys.  B70 (1974) 39–50.
  * (5) P. Fayet, _Supergauge Invariant Extension of the Higgs Mechanism and a Model for the electron and Its Neutrino_ , Nucl. Phys.  B90 (1975) 104–124.
  * (6) P. Fayet, _Spontaneously Broken Supersymmetric Theories of Weak, Electromagnetic and Strong Interactions_ , Phys. Lett.  69B (1977) 489.
  * (7) P. Fayet and S. Ferrara, _Supersymmetry_ , Phys. Rept.  32 (1977) 249–334.
  * (8) H. P. Nilles, _Supersymmetry, Supergravity and Particle Physics_ , Phys.Rept.  110 (1984) 1–162.
  * (9) H. E. Haber and G. L. Kane, _The Search for Supersymmetry: Probing Physics Beyond the Standard Model_ , Phys.Rept.  117 (1985) 75–263.
  * (10) M. Sohnius, _Introducing Supersymmetry_ , Phys.Rept.  128 (1985) 39–204.
  * (11) J. Gunion and H. E. Haber, _Higgs Bosons in Supersymmetric Models. 1._ , Nucl.Phys.  B272 (1986) 1.
  * (12) J. Gunion and H. E. Haber, _Higgs Bosons in Supersymmetric Models. 2. Implications for Phenomenology_ , Nucl.Phys.  B278 (1986) 449.
  * (13) R. Barbieri, S. Ferrara and C. A. Savoy, Phys. Lett. B 119 (1982), 343 doi:10.1016/0370-2693(82)90685-2.
  * (14) M. Dine, W. Fischler and M. Srednicki, Phys. Lett. B 104 (1981), 199-202 doi:10.1016/0370-2693(81)90590-6.
  * (15) H. P. Nilles, M. Srednicki and D. Wyler, Phys. Lett. B 120 (1983), 346 doi:10.1016/0370-2693(83)90460-4.
  * (16) J. M. Frere, D. R. T. Jones and S. Raby, Nucl. Phys. B 222 (1983), 11-19 doi:10.1016/0550-3213(83)90606-5.
  * (17) J. P. Derendinger and C. A. Savoy, Nucl. Phys. B 237 (1984), 307-328 doi:10.1016/0550-3213(84)90162-7.
  * (18) J. R. Ellis, J. F. Gunion, H. E. Haber, L. Roszkowski and F. Zwirner, Phys. Rev. D 39 (1989), 844 doi:10.1103/PhysRevD.39.844.
  * (19) M. Drees, Int. J. Mod. Phys. A 4 (1989), 3635 doi:10.1142/S0217751X89001448.
  * (20) U. Ellwanger, M. Rausch de Traubenberg and C. A. Savoy, Phys. Lett. B 315 (1993), 331-337 doi:10.1016/0370-2693(93)91621-S [arXiv:hep-ph/9307322 [hep-ph]].
  * (21) U. Ellwanger, M. Rausch de Traubenberg and C. A. Savoy, Z. Phys. C 67 (1995), 665-670 doi:10.1007/BF01553993 [arXiv:hep-ph/9502206 [hep-ph]].
  * (22) U. Ellwanger, M. Rausch de Traubenberg and C. A. Savoy, Nucl. Phys. B 492 (1997), 21-50 doi:10.1016/S0550-3213(97)00128-4 [arXiv:hep-ph/9611251 [hep-ph]].
  * (23) T. Elliott, S. F. King and P. L. White, Phys. Lett. B 351 (1995), 213-219 doi:10.1016/0370-2693(95)00381-T [arXiv:hep-ph/9406303 [hep-ph]].
  * (24) S. F. King and P. L. White, Phys. Rev. D 52 (1995), 4183-4216 doi:10.1103/PhysRevD.52.4183 [arXiv:hep-ph/9505326 [hep-ph]].
  * (25) F. Franke and H. Fraas, Int. J. Mod. Phys. A 12 (1997), 479-534 doi:10.1142/S0217751X97000529 [arXiv:hep-ph/9512366 [hep-ph]].
  * (26) M. Maniatis, Int. J. Mod. Phys. A 25 (2010), 3505-3602 doi:10.1142/S0217751X10049827 [arXiv:0906.0777 [hep-ph]].
  * (27) U. Ellwanger, C. Hugonie and A. M. Teixeira, Phys. Rept. 496 (2010), 1-77 doi:10.1016/j.physrep.2010.07.001 [arXiv:0910.1785 [hep-ph]].
  * (28) J. E. Kim and H. P. Nilles, Phys. Lett. B 138 (1984), 150-154 doi:10.1016/0370-2693(84)91890-2.
  * (29) G. Apollinari, I. Béjar Alonso, O. Brüning, P. Fessia, M. Lamont, L. Rossi and L. Tavian, doi:10.23731/CYRM-2017-004.
  * (30) H. Baer, T. Barklow, K. Fujii, Y. Gao, A. Hoang, S. Kanemura, J. List, H. E. Logan, A. Nomerotski and M. Perelstein, et al. [arXiv:1306.6352 [hep-ph]].
  * (31) M. Bicer et al. [TLEP Design Study Working Group], JHEP 01 (2014), 164 doi:10.1007/JHEP01(2014)164 [arXiv:1308.6176 [hep-ex]].
  * (32) M. Ahmad, D. Alves, H. An, Q. An, A. Arhrib, N. Arkani-Hamed, I. Ahmed, Y. Bai, R. B. Ferroli and Y. Ban, et al. IHEP-CEPC-DR-2015-01.
  * (33) P. Slavich, S. Heinemeyer, E. Bagnaschi, H. Bahl, M. Goodsell, H. E. Haber, T. Hahn, R. Harlander, W. Hollik and G. Lee, et al. [arXiv:2012.15629 [hep-ph]].
  * (34) F. Domingo, S. Heinemeyer, S. Paßehr and G. Weiglein, Eur. Phys. J. C 78 (2018) no.11, 942 doi:10.1140/epjc/s10052-018-6400-1 [arXiv:1807.06322 [hep-ph]].
  * (35) F. Domingo, S. Heinemeyer, S. Paßehr and G. Weiglein, _Decays of the neutral Higgs bosons into SM fermions and gauge bosons in the $\mathcal{CP}$-violating NMSSM_, _Eur. Phys. J. C_ 78 (2018) 942, [1807.06322].
  * (36) F. Domingo and S. Paßehr, Eur. Phys. J. C 79, no.11, 905 (2019) doi:10.1140/epjc/s10052-019-7415-y [arXiv:1907.05468 [hep-ph]].
  * (37) D. T. Nhung, M. Muhlleitner, J. Streicher and K. Walz, JHEP 11, 181 (2013) doi:10.1007/JHEP11(2013)181 [arXiv:1306.3926 [hep-ph]].
  * (38) G. Bélanger, V. Bizouard, F. Boudjema and G. Chalons, Phys. Rev. D 96, no.1, 015040 (2017) doi:10.1103/PhysRevD.96.015040 [arXiv:1705.02209 [hep-ph]].
  * (39) M. Mühlleitner, D. T. Nhung and H. Ziesche, JHEP 12, 034 (2015) doi:10.1007/JHEP12(2015)034 [arXiv:1506.03321 [hep-ph]].
  * (40) J. Baglio, T. N. Dao and M. Mühlleitner, Eur. Phys. J. C 80 (2020) no.10, 960 doi:10.1140/epjc/s10052-020-08520-8 [arXiv:1907.12060 [hep-ph]].
  * (41) T. N. Dao, L. Fritz, M. Krause, M. Mühlleitner and S. Patel, Eur. Phys. J. C 80 (2020) no.3, 292 doi:10.1140/epjc/s10052-020-7837-6 [arXiv:1911.07197 [hep-ph]].
  * (42) F. Domingo and S. Paßehr, Eur. Phys. J. C 80 (2020) no.12, 1124 doi:10.1140/epjc/s10052-020-08655-8 [arXiv:2007.11010 [hep-ph]].
  * (43) J. Baglio, R. Gröber, M. Mühlleitner, D. T. Nhung, H. Rzehak, M. Spira, J. Streicher and K. Walz, Comput. Phys. Commun. 185 (2014) no.12, 3372-3391 doi:10.1016/j.cpc.2014.08.005 [arXiv:1312.4788 [hep-ph]].
  * (44) M. Spira, Prog. Part. Nucl. Phys. 95 (2017), 98-159 doi:10.1016/j.ppnp.2017.04.001 [arXiv:1612.07651 [hep-ph]].
  * (45) L. J. Hall, R. Rattazzi and U. Sarid, Phys. Rev. D 50 (1994), 7048-7065 doi:10.1103/PhysRevD.50.7048 [arXiv:hep-ph/9306309 [hep-ph]].
  * (46) R. Hempfling, Phys. Rev. D 49 (1994), 6168-6172 doi:10.1103/PhysRevD.49.6168.
  * (47) M. Carena, M. Olechowski, S. Pokorski and C. E. M. Wagner, Nucl. Phys. B 426 (1994), 269-300 doi:10.1016/0550-3213(94)90313-1 [arXiv:hep-ph/9402253 [hep-ph]].
  * (48) D. M. Pierce, J. A. Bagger, K. T. Matchev and R. j. Zhang, Nucl. Phys. B 491 (1997), 3-67 doi:10.1016/S0550-3213(96)00683-9 [arXiv:hep-ph/9606211 [hep-ph]].
  * (49) M. Carena, S. Mrenna and C. E. M. Wagner, Phys. Rev. D 60 (1999), 075010 doi:10.1103/PhysRevD.60.075010 [arXiv:hep-ph/9808312 [hep-ph]].
  * (50) M. Carena, D. Garcia, U. Nierste and C. E. M. Wagner, Nucl. Phys. B 577 (2000), 88-120 doi:10.1016/S0550-3213(00)00146-2 [arXiv:hep-ph/9912516 [hep-ph]].
  * (51) M. Carena, J. R. Ellis, S. Mrenna, A. Pilaftsis and C. E. M. Wagner, Nucl. Phys. B 659 (2003), 145-178 doi:10.1016/S0550-3213(03)00192-5 [arXiv:hep-ph/0211467 [hep-ph]].
  * (52) J. Guasch, P. Hafliger and M. Spira, Phys. Rev. D 68 (2003), 115001 doi:10.1103/PhysRevD.68.115001 [arXiv:hep-ph/0305101 [hep-ph]].


# The neutron star mass, distance, and inclination from precision timing of
the brilliant millisecond pulsar J0437$-$4715

Daniel J. Reardon Centre for Astrophysics and Supercomputing, Swinburne
University of Technology, P.O. Box 218, Hawthorn, Victoria 3122, Australia
OzGrav: The Australian Research Council Centre of Excellence for Gravitational
Wave Discovery, Hawthorn VIC 3122, Australia  Matthew Bailes Centre for
Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box
218, Hawthorn, Victoria 3122, Australia  OzGrav: The Australian Research
Council Centre of Excellence for Gravitational Wave Discovery, Hawthorn VIC
3122, Australia  Ryan M. Shannon Centre for Astrophysics and Supercomputing,
Swinburne University of Technology, P.O. Box 218, Hawthorn, Victoria 3122,
Australia  OzGrav: The Australian Research Council Centre of Excellence for
Gravitational Wave Discovery, Hawthorn VIC 3122, Australia  Chris Flynn Centre
for Astrophysics and Supercomputing, Swinburne University of Technology, P.O.
Box 218, Hawthorn, Victoria 3122, Australia  OzGrav: The Australian Research
Council Centre of Excellence for Gravitational Wave Discovery, Hawthorn VIC
3122, Australia  Jacob Askew Centre for Astrophysics and Supercomputing,
Swinburne University of Technology, P.O. Box 218, Hawthorn, Victoria 3122,
Australia  OzGrav: The Australian Research Council Centre of Excellence for
Gravitational Wave Discovery, Hawthorn VIC 3122, Australia  N. D. Ramesh Bhat
International Centre for Radio Astronomy Research, Curtin University, Bentley,
WA 6102, Australia  Zu-Cheng Chen Department of Physics and Synergetic
Innovation Center for Quantum Effects and Applications, Hunan Normal
University, Changsha, Hunan 410081, China  Institute of Interdisciplinary
Studies, Hunan Normal University, Changsha, Hunan 410081, China  Małgorzata
Curyło Astronomical Observatory, University of Warsaw, Aleje Ujazdowskie 4,
00-478 Warsaw, Poland  Yi Feng Research Center for Astronomical Computing,
Zhejiang Laboratory, Hangzhou 311100, China  George B. Hobbs Australia
Telescope National Facility, CSIRO, Space and Astronomy, P.O. Box 76, Epping,
NSW 1710, Australia  Agastya Kapur Australia Telescope National Facility,
CSIRO, Space and Astronomy, P.O. Box 76, Epping, NSW 1710, Australia  Matthew
Kerr Space Science Division, US Naval Research Laboratory, 4555 Overlook Ave
SW, Washington DC 20375, USA  Xiaojin Liu Faculty of Arts and Sciences,
Beijing Normal University, Zhuhai 519087, People‘s Republic of China  Richard
N. Manchester Australia Telescope National Facility, CSIRO, Space and
Astronomy, P.O. Box 76, Epping, NSW 1710, Australia  Rami Mandow Department of
Physics and Astronomy and MQ Research Centre in Astronomy, Astrophysics and
Astrophotonics, Macquarie University, NSW 2109, Australia  Australia Telescope
National Facility, CSIRO, Space and Astronomy, P.O. Box 76, Epping, NSW 1710,
Australia  Saurav Mishra Centre for Astrophysics and Supercomputing, Swinburne
University of Technology, P.O. Box 218, Hawthorn, Victoria 3122, Australia
Christopher J. Russell CSIRO Scientific Computing, Australian Technology Park,
Locked Bag 9013, Alexandria, NSW 1435, Australia  Mohsen Shamohammadi Centre
for Astrophysics and Supercomputing, Swinburne University of Technology, P.O.
Box 218, Hawthorn, Victoria 3122, Australia  OzGrav: The Australian Research
Council Centre of Excellence for Gravitational Wave Discovery, Hawthorn VIC
3122, Australia  Lei Zhang National Astronomical Observatories, Chinese
Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100101,
People‘s Republic of China  Centre for Astrophysics and Supercomputing,
Swinburne University of Technology, P.O. Box 218, Hawthorn, Victoria 3122,
Australia  Andrew Zic Australia Telescope National Facility, CSIRO, Space and
Astronomy, P.O. Box 76, Epping, NSW 1710, Australia

(Received 2024 April 3; Revised 2024 June 12; Accepted 2024 July 1)

###### Abstract

The observation of neutron stars enables the otherwise impossible study of
fundamental physical processes. The timing of binary radio pulsars is
particularly powerful, as it enables precise characterization of their (three-
dimensional) positions and orbits. PSR J0437$-$4715 is an important
millisecond pulsar for timing array experiments and is also a primary target
for the Neutron Star Interior Composition Explorer (NICER). The main aim of
the NICER mission is to constrain the neutron star equation of state by
inferring the compactness ($M_{p}/R$) of the star. Direct measurements of the
mass $M_{p}$ from pulsar timing therefore substantially improve constraints on
the radius $R$ and the equation of state. Here we use observations spanning 26
years from Murriyang, the 64-m Parkes radio telescope, to improve the timing
model for this pulsar. Among the new precise measurements are the pulsar mass
$M_{p}=1.418\pm 0.044$ $M_{\odot}$, distance $D=156.96\pm 0.11$ pc, and
orbital inclination angle $i=137.506\pm 0.016^{\circ}$, which can be used to
inform the X-ray pulse profile models inferred from NICER observations. We
demonstrate that these results are consistent between multiple data sets from
the Parkes Pulsar Timing Array (PPTA), each modeled with different noise
assumptions. Using the longest available PPTA data set, we measure an apparent
second derivative of the pulsar spin frequency and discuss how this can be
explained either by kinematic effects due to the proper motion and radial
velocity of the pulsar or excess low-frequency noise such as a gravitational-
wave background.

dense matter — equation of state — millisecond pulsars — pulsars: general —
pulsars: individual (PSR J0437$-$4715) — pulsar timing method

††journal: ApJ Letters

##  1 Introduction

The equation of state (EoS) of neutron stars is an important probe of the
physics at supranuclear densities and low (nonrelativistic) temperatures that
are inaccessible to laboratory experiments (Özel et al., 2016). Precision
radio-frequency timing of millisecond pulsars (MSPs) can place stringent
constraints on the EoS through neutron star mass measurements (e.g. Demorest
et al., 2010; Cromartie et al., 2020; Fonseca et al., 2021).

Pulsar timing arrays (PTAs) regularly monitor a set of MSPs over years,
primarily to act as a Galactic-scale detector of nanohertz-frequency
gravitational waves (e.g. Agazie et al., 2023a; Antoniadis et al., 2023a;
Reardon et al., 2023a; Xu et al., 2023). The long-duration PTA data sets are
also ideally suited for studying the astrometry and orbits of MSPs, leading to
tests of theories of gravity and the physics of neutron star interiors. The
Parkes Pulsar Timing Array (PPTA) was the first such long-term PTA observing
project (Manchester et al., 2013). One of the most important pulsars in the
array is the nearest and brightest MSP, J0437$-$4715 (Johnston et al., 1993).
The proximity, long timing baseline ($\gtrsim 20$ yr) and high-precision
arrival times for PSR J0437$-$4715 provide one of the best opportunities to
constrain the neutron star EoS, in combination with a radius measurement.

The Neutron Star Interior Composition Explorer (NICER; Gendreau et al., 2012,
2016) is a NASA telescope on the International Space Station that observes and
times MSPs in soft X-rays ($\sim 0.3$ \- $3$ keV). The spectral-timing pulse
profile data from NICER is modeled to infer X-ray emission region geometry and
relativistic effects that yield the neutron star mass and radius. NICER has
already constrained the mass and radius for the isolated PSR J0030+0451 (Riley
et al., 2019; Miller et al., 2019) and the more massive PSR J0740+6620 (Riley
et al., 2021; Miller et al., 2021) by leveraging a mass constraint from radio-
frequency pulsar timing (Cromartie et al., 2020; Fonseca et al., 2021).
However, PSR J0437$-$4715 remains a primary target for the mission because of
its proximity, brightness, and the precise independent mass measurement
enabled by the extremely accurate pulsar timing. The mass of PSR J0437$-$4715
is near the peak of the neutron star mass distribution (Özel & Freire, 2016),
meaning that a radius measurement will constrain the neutron star EoS at the
crucial intermediate densities that describe much of the population. A
significant challenge in determining the mass and radius of PSR J0437$-$4715
from NICER data alone is the presence of a bright active galactic nucleus in
the field that contributes additional background noise. The NICER models are
improved further by prior knowledge of the pulsar distance and orbital
inclination angle, which can also be inferred from PPTA observations (Reardon
et al., 2016).

The times of arrival (ToAs) of MSPs reveal not just the deterministic physics
captured by timing model parameters but also numerous stochastic processes
(Agazie et al., 2023b; Reardon et al., 2023b; Antoniadis et al., 2023b).
Examples of these processes include interstellar dispersion measure (DM) and
scattering variations (e.g. Keith et al., 2013; Donner et al., 2020; Turner et
al., 2021), pulsar rotation irregularities (e.g. Shannon & Cordes, 2010),
pulse jitter (e.g. Liu et al., 2012; Shannon et al., 2014; Lam et al., 2019;
Parthasarathy et al., 2021), radio-frequency interference (RFI), sudden and
secular profile shape variations (e.g. Shannon et al., 2016; Jennings et al.,
2024), calibration errors in specific observing systems (e.g. van Straten,
2013; Lentati et al., 2016), and nanohertz-frequency gravitational waves. PSR
J0437$-$4715, shows evidence for most of these noise processes. Determining
precise and accurate measurements of the neutron star mass, distance, and
inclination requires both accurate ToAs and noise models.

Throughout the NICER mission, the modeling of PSR J0437$-$4715 observations
has been conducted using informative prior probability distributions derived
from PPTA observations. In this Letter, we describe these PPTA data sets and
noise models (Section 2), the timing model parameters and inference technique
(Section 3), and our resulting measurements and derived quantities (Section
4). Our results are discussed in Section 5 and we conclude in Section 6.

##  2 Data sets

Pulsar timing model parameters can change significantly between different
analyses. For this work, we critically analyze different PPTA data sets to
demonstrate that our results are robust. The measured ToAs as a function of
time and frequency, for all data sets discussed in this section, are shown in
Figure 1.

Figure 1: Measured ToAs for PSR J0437$-$4715 across multiple PPTA data sets.
The colors of the points indicate the observing frequency on a linear scale.
The data set is labelled on each panel, and described in Section 2.

###  2.1 PPTA first data release

The PPTA commenced regular observations of a set of MSPs in 2004 and produced
a first data release of 20 pulsars after six years of observations (Manchester
et al., 2013). This first data release (PPTA-DR1) was combined with earlier
“legacy” observations of many pulsars taken under different projects, which
for PSR J0437$-$4715 dated back to 1996 (Verbiest et al., 2008). The
combination of these data with PPTA-DR1 is referred to as the “extended” first
data release (PPTA-DR1e). A timing and noise analysis of the 20 MSPs in PPTA-
DR1e, including PSR J0437$-$4715, was conducted by Reardon et al. (2016) and
included precise measurements of the pulsar mass, $M_{p}$, and distance, $D$.
For PPTA-DR1 and all subsequent data releases to date, the pulsar was timed
using the polarization invariant interval (Britton, 2000) to mitigate
polarization calibration errors (Manchester et al., 2013).

###  2.2 PPTA second data release: Noise model and legacy data extension

The PPTA second data release (PPTA-DR2) is described in Kerr et al. (2020) and
included observations of 26 pulsars, spanning up to 14 years. The timing
residuals were analyzed in detail for noise processes by Goncharov et al.
(2021). Despite using the most comprehensive noise model ever created for a
pulsar, the normalized and noise-subtracted timing residuals for PSR
J0437$-$4715 failed to pass any of the statistical tests used to assess if the
model fully characterizes the noise. These challenges have been noted for PSR
J0437$-$4715 during numerous earlier analyses (e.g. Verbiest et al., 2008; van
Straten, 2013; Lentati et al., 2014b, 2016; Reardon et al., 2016).

For a detailed description of the PPTA-DR2 noise models, see Goncharov et al.
(2021). In brief, the noise models include white-noise parameters (that
account for uncorrelated instrumental noise and pulse jitter) and Gaussian
random processes with a power-law power spectral density (PSD; e.g. van
Haasteren & Levin, 2013; Lentati et al., 2013). These Gaussian processes are
used to account for achromatic noise (for example spin noise and gravitational
waves), DM variations (scaling in amplitude $\propto\nu^{-2}$, for observing
frequency $\nu$), scattering variations (scaling with $\nu^{-4}$), and excess
achromatic noise in specific observing bands (band noise) and observing
systems (system noise). The parameters inferred for each such process are the
PSD power law amplitude, $A$, and spectral exponent, $\gamma$.

The PPTA-DR2 data set was then extended with the legacy observations in
Reardon et al. (2021), forming PPTA-DR2e, for all pulsars except PSR
J0437$-$4715\. The timing models for these pulsars were updated, with the
parameters inferred in the presence of the Goncharov et al. (2021) noise
models with least-squares regression in the pulsar timing package tempo2
(Hobbs et al., 2006; Edwards et al., 2006). The white-noise parameters, as
well as $A$ and $\gamma$ for each Gaussian process, were fixed at their
inferred maximum likelihood values. The Gaussian processes are implemented in
the timing model as a Fourier series with $N$ component frequencies, from
$f_{\rm min}=1/T_{\rm span}$ to $f_{\rm max}=N/T_{\rm span}$, for total data
set time span, $T_{\rm span}$. The amplitudes of the Fourier basis functions
are included in the timing model but constrained to follow the PSD described
by $A$ and $\gamma$. The noise models from PPTA-DR2 were extrapolated for use
with PPTA-DR2e by increasing the number of Fourier components, $N$, such that
the highest-frequency component, $f_{\rm max}$, was similar for both data sets
(Reardon et al., 2021).

For this work, we repeated this process for PSR J0437$-$4715, to produce a
22-year PPTA-DR2e data set and an extrapolated noise model. The timing model
in this data set assumes the DE436 solar system ephemeris model (an updated
version of DE430; Folkner et al., 2014; Park et al., 2021) since this was
assumed for the noise analysis (Goncharov et al., 2021). As with PPTA-DR2, the
noise model whitens the residuals on many timescales, but the noise-subtracted
residuals remain non-Gaussian and heteroscedastic. This is likely due to
changes in the noise statistics of the multiple observing systems over the two
decades of PPTA observations. While the effect of such noise model
deficiencies on the timing model parameters are expected to be small, based on
attempts to model the timing residuals with a non-Gaussian likelihood (Lentati
et al., 2014b), we endeavored to improve the residuals before these parameters
were adopted for the NICER mission.

###  2.3 Reprocessing the PPTA second data release: PPTA-DR2.5

To better understand the incomplete model of noise in PPTA-DR2, we produced an
ad hoc data set by reprocessing the PPTA-DR2 pulse profiles of PSR
J0437$-$4715, with the goal to first produce the highest-quality ToAs and
uncertainties. We refer to this data set as “PPTA-DR2.5” as it was produced
specifically to achieve new timing model measurements for the NICER mission,
prior to the publication of the PPTA third data release (PPTA-DR3).

The published PPTA-DR2 data set included calibrated (polarization and flux
density) and RFI-excised pulse profiles with 32 frequency channels and eight
subintegrations in time. We developed a pipeline extension to the PPTA-DR2
data reduction pipeline to reprocess these pulse profiles. This extension was
used in PPTA-DR3, and it is therefore also discussed in Zic et al. (2023,
Section 2.2). Below we summarize the key pipeline changes and two differences
between our ad hoc PPTA-DR2.5 (produced only for PSR J0437$-$4715) and the
published PPTA-DR3:

  * •

Both PPTA-DR2.5 and PPTA-DR3 only include observations from the more modern
signal processing systems (backends) as the earlier systems contribute
unidentified excess noise to the ToAs (Goncharov et al., 2021). The data sets
therefore start in July 2007 with observations recorded with the PDFB2
backend. PPTA-DR3 additionally includes $\sim 3.3$ years of observations from
the ultrawide-band low (UWL) receiver (Hobbs et al., 2020).

  * •

PPTA-DR2.5 retained all 32 frequency channels of the pulse profiles, while
PPTA-DR3 was reduced to four channels per observation (or per sub-band for UWL
observations). These differences are apparent in Figure 1.

  * •

For each combination of observing receiver and backend, we created a separate
template profile by summing the brightest 100 observations in time and
smoothing the resulting 32-channel profile with wavelet smoothing using the
psrsmooth of the psrchive software package (van Straten et al., 2012). For
PPTA-DR3, a single analytical wide-band profile (derived from UWL
observations; Curyło et al., 2023) was used.

  * •

Each channel in the pulse profile was timed with the corresponding channel
from the template, with the ToA and uncertainty measured with the Fourier-
domain Monte Carlo method in pat (van Straten et al., 2012).

The noise in PPTA-DR2.5 was modeled using temponest (Lentati et al., 2014a).
Power-law models were used to describe the DM variations, achromatic red
noise, and band noise in two bands: $\nu<1$ GHz and $1\,{\rm GHz}<\nu<2$ GHz.
White-noise parameters were inferred for each combination of receiver and
backend. We found that PPTA-DR2.5 did not require any additional timing model
parameters to account for profile evolution with frequency, while PPTA-DR3 did
(Zic et al., 2023).

Following this noise analysis of PPTA-DR2.5, the frequency-averaged, whitened,
and normalized residuals were found to be statistically white and Gaussian
distributed using the same tests as in Goncharov et al. (2021); Reardon et al.
(2021). Inclusion of ToAs from earlier observing systems caused these tests to
fail, which motivated the use of only these observing systems for PPTA-DR3.

###  2.4 PPTA third data release: PPTA-DR3

Zic et al. (2023) describes the third PPTA data release, which was analyzed
for noise processes in Reardon et al. (2023b), using the enterprise package
(Ellis et al., 2019). We do not further modify the noise model or ToAs from
those already published. The noise analysis for PPTA-DR3 assumed the DE440 JPL
solar system ephemeris (Park et al., 2021).

The timing residuals for PSR J0437$-$4715 from PPTA-DR3 are shown in Figure 2.
The frequency-averaged, noise-subtracted (whitened), and normalized residuals
in the bottom panel pass the same tests as PPTA-DR2.5 and so are also
statistically white and Gaussian distributed as expected from a complete noise
model.

Figure 2: Band-averaged timing residuals for PSR J0437$-$4715 from PPTA-DR3.
Residuals at observing frequencies near 800 MHz, 1400 MHz, and 3100 MHz are
shown in red, teal, and blue, respectively. Top panel shows post-fit residuals
with all modeled Gaussian random processes present, including the dominant DM
variations that increase in amplitude $\propto\nu^{-2}$. The middle panel
shows the residuals after removal of the maximum likelihood realizations of
these Gaussian processes. The bottom panel shows these noise-subtracted
residuals divided by the uncertainty. The top and middle panels have labels
that show the weighted rms of the residuals in each band.

##  3 Timing model

PSR J0437$-$4715 is in a $5.74$-day, near-circular orbit (eccentricity, $e\sim
1.9\times 10^{-5}$) with a helium white dwarf companion Bell et al. (1993).
The brightness and proximity of the pulsar, as well as the geometry of the
orbit, are favorable for measuring useful relativistic and kinematic effects
from precisely measured ToAs. The timing model implemented in tempo2 is
described in detail in Edwards et al. (2006), with the post-Newtonian orbital
mechanics derived in Damour & Deruelle (1986). We use the T2 binary model from
tempo2, which is an extension of the model described in Damour & Taylor (1992)
to include kinematic effects (e.g. Section 3.2). Here we provide a summary of
the most important timing model parameters relevant to this work.

###  3.1 Shapiro Delay

General relativity predicts that spacetime can be measurably compressed around
massive and compact objects. As radio waves from PSR J0437$-$4715 propagate
through the compressed spacetime surrounding the white dwarf companion, they
are delayed relative to propagation through uncompressed spacetime. This is
the well-known Shapiro delay (Shapiro, 1964), and it is observed most readily
in binary MSPs with orbits that are observed near edge-on (highly inclined).
It is parameterized in the timing model of PSR J0437$-$4715 by the companion
mass $M_{c}$ and shape $s=\sin{i}$, where $i$ is the orbital inclination
angle. A precise measure of $M_{c}$ and $\sin{i}$ therefore yields the pulsar
mass $M_{p}$ through the binary mass function. However, $i$ can also be
measured independently through kinematic effects that change the projected
shape of the orbit on the sky, thus eliminating the strong covariance between
$M_{c}$ and $i$ due to Shapiro delay.

###  3.2 Annual-orbital parallax

For the nearest and most precisely timed binary pulsars, it is possible to
measure annual variations in the projected shape of the orbit. This is the
annual-orbital parallax and leads to variations in the projected semi-major
axis, $x=a\sin{i}$ (where $a$ is the semi-major axis), and longitude of
periastron $\omega$ (Kopeikin, 1995). In addition, secular variations in $x$
and $\omega$ originate from the proper motion of the binary slowly altering
the projection of the orbit (Kopeikin, 1996). These effects are parametrized
in tempo2 by the orbital inclination angle $i$ (which simultaneously describes
the Shapiro delay shape) and the longitude of ascending node, $\Omega$. Thus,
these “Kopeikin terms” account for the annual-orbital parallax, the kinematic
contributions to the time derivatives of $x$ and $\omega$, and also tightly
constrain the Shapiro delay.

Additional contributions to the observed time derivatives of $x$ and $\omega$
can arise, for example, from relativistic effects, or spin-orbit coupling. For
PSR J0437$-$4715, we detect a significant relativistic periastron advance,
$\dot{\omega}$, and include this in the binary model. However, the only
measurable component to $\dot{x}$ is the kinematic term that is accounted for
with the Kopeikin terms; thus, we do not include an additional $\dot{x}$ in
the model.

###  3.3 Shklovskii and parallax distances

The proper motion of a pulsar not only changes the projection of the orbit but
is also associated with an apparent radial acceleration that causes an
increasing Doppler shift to the spin and orbital periods. This Shklovskii
effect (Shklovskii, 1970) can be used to derive the pulsar distance from a
measurement of the resulting orbital period derivative $\dot{P}_{b}^{\rm
shk}$, with (Bell & Bailes, 1996)

$D_{\rm shk}=\frac{c}{\mu^{2}}\frac{\dot{P}_{b}^{\rm shk}}{P_{b}},$ (1)

where $\mu$ is the magnitude of the proper motion, $P_{b}$ is the orbital
period, and $c$ is the speed of light. However, the observed orbital period
derivative, $\dot{P}_{b}^{\rm obs}$, has additional contributions, including
from the radiation of gravitational waves by the system, $\dot{P}_{b}^{\rm
GW}$, and from differential acceleration in the Galaxy, $\dot{P}_{b}^{\rm
Gal}$. The Shklovskii component can be computed by correcting for these
effects

$\dot{P}_{b}^{\rm shk}=\dot{P}_{b}^{\rm obs}-\dot{P}_{b}^{\rm
GW}-\dot{P}_{b}^{\rm Gal}.$ (2)

We use galpy (Bovy, 2015) to measure and correct for the Galactic acceleration
term, $\dot{P}_{b}^{\rm Gal}$. We propagate uncertainties in this correction
by assuming a 3% uncertainty for the circular differential acceleration (Bovy
et al., 2012) and 10% for the vertical acceleration in the Galactic potential
(Holmberg & Flynn, 2004), finding $\dot{P}_{b}^{\rm Gal}=-(2.24\pm 0.20)\times
10^{-14}$. The expected $\dot{P}_{b}$ from gravitational-wave emission is two
orders of magnitude smaller, $\dot{P}_{b}^{\rm GW}\sim-3.2\times 10^{-16}$.

An independent distance can be found from the inverse of parallax, $\pi$,
which in pulsar timing is a detection of the time delays of curved wave fronts
(with respect to plane waves) arriving from a finite distance. However, the
sensitivity of a pulsar timing data set to parallax $D_{\pi}\propto T_{\rm
span}^{1/2}$ does not grow as rapidly with time as $D_{\rm shk}$, which
improves $\propto T_{\rm span}^{5/2}$.

###  3.4 Radial velocity

The radial motion of the pulsar causes an apparent second derivative of the
pulsar spin period, $\ddot{f}$ (van Straten, 2003; Liu et al., 2018). It also
causes a proper motion derivative (van Straten, 2003), which will have a
smaller effect on the timing residuals. However, the observed $\ddot{f}$
includes contributions from multiple other sources,

$\ddot{f}_{\rm obs}=\left(1-\frac{v_{r}}{c}\right)\ddot{f}_{\rm
int}+\ddot{f}_{\rm shk}+\ddot{f}_{\parallel}+\ddot{f}_{\perp}+\ddot{f}_{\rm
acc}+\ddot{f}_{\rm jerk}+\ddot{f}_{\rm red},$ (3)

where $v_{r}$ is the radial velocity; $\ddot{f}_{\rm int}$ comes from the
intrinsic spin-down of the pulsar; $\ddot{f}_{\rm shk}$ originates from the
Shklovskii effect acting on the intrinsic $\dot{f}$; $\ddot{f}_{\parallel}$
and $\ddot{f}_{\perp}$ are additional Shklovskii terms that depend on the
radial velocity and transverse acceleration, respectively; and $\ddot{f}_{\rm
acc}$ and $\ddot{f}_{\rm jerk}$ depend, respectively, on differential
acceleration and jerk (for details on these terms, see Liu et al., 2018).
$\ddot{f}_{\rm red}$ is an additional noise component due to unmodeled low-
frequency (red) noise in the timing residuals.

For PSR J0437$-$4715, the dominant term (assuming an accurate noise model,
$\ddot{f}_{\rm red}=0$) in Equation 3 is $\ddot{f}_{\parallel}$, which is
proportional to the radial velocity. Under typical spin-down scenarios, the
expected intrinsic $\ddot{f}_{\rm int}$ for low-magnetic-field MSPs is small.
This and the other terms in Equation 3 are estimated to contribute $\lesssim
1\%$ to $\ddot{f}_{\rm obs}$ (Liu et al., 2018, Table 1). For this estimate,
the jerk and acceleration terms were calculated for the Galactic potential
only. Interactions with field stars or a second companion in a wide orbit
could potentially increase these components, but this is unlikely to be
significant for PSR J0437$-$4715\. One Galactic (field) MSP was found to be in
a wide binary orbit that manifested as a significant $\ddot{f}$ due to
gravitational jerk (Bassa et al., 2016; Kaplan et al., 2016). However, unlike
PSR J0437$-$4715, this pulsar also had an apparent $\dot{f}$ that was
inconsistent with spin-down alone because of a significant contribution from
gravitational acceleration, and the wide orbit companion was identified
optically.

The radial velocity is then calculated assuming
$\ddot{f}_{\parallel}\approx\ddot{f}_{\rm obs}$ with

$v_{r}=\frac{\ddot{f}_{\parallel}}{f_{0}}\frac{c}{3\mu^{2}}.$ (4)

However, covariance between $\ddot{f}$ and red noise processes makes the
measurement challenging, particularly if the noise has a steep spectrum (Liu
et al., 2019; Keith & Niţu, 2023). Indeed leakage of red noise into the
$\ddot{f}$ measurement due to model misspecification, $\ddot{f}_{\rm red}$,
may be one of the dominant terms in Equation 3. The precision in $\ddot{f}$,
and hence $v_{r}$, is estimated to improve $\propto T_{\rm span}^{\sim 2}$ for
PSR J0437$-$4715, after accounting for covariance with the achromatic red
noise (Liu et al., 2019), which has spectral exponent $\gamma\sim-3$
(Goncharov et al., 2021; Reardon et al., 2023b). If the red noise in PSR
J0437$-$4715 becomes dominated by the steeper “common-spectrum” process
observed using PPTA-DR3 (which may be a gravitational-wave background; Reardon
et al., 2023a), then the precision in $\ddot{f}$ would begin to improve more
slowly, $\propto T_{\rm span}^{\sim 1.5}$, for $\gamma\sim-3.9$. With the
present sensitivity of PPTA data sets, the magnitude of $v_{r}$ can only be
estimated under an assumed noise model. Confirmation of $v_{r}$ from pulsar
timing will be possible upon measurement of the proper motion derivative
(requiring $\gtrsim 30$ years of timing baseline), which has a purely
geometrical origin (van Straten, 2003).

###  3.5 Parameter inference

Bayesian inference is commonly used for noise modeling of MSPs but can be
computationally expensive for large data sets like those of the PSR
J0437$-$4715 ToAs. This is especially true when greatly increasing the
dimensionality of the parameter space by choosing to sample over timing model
parameters. Fortunately, for the case of timing model parameters that are
measured with high significance, generalized least-squares fitting of the
linearized timing model in codes such as tempo2 is an efficient alternative
and provides accurate measurements and uncertainties, provided the noise model
is well specified and the fixed noise parameters (e.g., the amplitude $A$ and
spectral exponent $\gamma$ of a power-law process) are uncorrelated with the
timing model parameters. We verified this using temponest to sample $M_{c}$,
$\dot{P}_{b}$, and $i$ simultaneously with the noise model for PPTA-DR2.5,
finding negligible correlations between the parameters. We further verified
that the results of the least-squares fit in tempo2 (including the noise
model) were consistent with the Bayesian posteriors (see also Lentati et al.,
2014a).

For the analysis of all data sets, we therefore used this generalized least-
squares method after validation of the noise models. The Gaussian process
components of the noise models are realized simultaneously with the timing
model fit, as described in Section 2.2 for the PPTA-DR2e analysis of other
MSPs in the PPTA (Reardon et al., 2021). To derive further quantities from our
measured timing model parameters (e.g. $M_{p}$ and $D$), we perform a Monte
Carlo simulation of our measurements by drawing 1 million samples from
Gaussian distributions described by the timing model parameter values and
their uncertainties. For each set of parameter samples drawn, we derived the
relevant quantity. The quoted derived quantities in Section 4 represent the
mean and standard deviation of the distribution of these derived quantities.

##  4 Results

We present timing model parameter measurements derived from the two most
sensitive data sets: PPTA-DR2e and PPTA-DR3. As the data set with the longest
timing baseline, PPTA-DR2e is more sensitive to parameters that grow rapidly
with $T_{\rm span}$. PPTA-DR3 has higher-quality ToAs from modern backends and
the UWL receiver, as well as a complete noise model suitable for accurate
measurements of effects at the smaller timescales, like the Shapiro delay and
annual-orbital parallax. The parameters inferred from these data sets are
presented in Table 1.

The DM is taken from PPTA-DR2.5, which did not require the use of parameters
that describe frequency-dependent (FD) pulse profile evolution (Arzoumanian et
al., 2018). DM derivative terms, magnetospheric exponential dip parameters
(Goncharov et al., 2021), observing system jumps, and FD parameters are not
presented here because we consider these part of the noise models. The ToAs
from each data set and the complete pulsar ephemerides, which includes system
jumps, noise model parameters, and other ancillary tempo2 parameters, are
available from the CSIRO Data Access Portal (Hobbs et al., 2011) at
https://doi.org/10.25919/20rx-5f63.

Table 1: Measured timing model parameters and derived quantities for PSR
J0437$-$4715 from PPTA-DR2e and PPTA-DR3. Values in parentheses are the
1-$\sigma$ standard errors on the last quoted decimal place. Values in bold
are the most precise of the two data sets. †The measured DM comes from PPTA-
DR2.5, which did not require profile frequency evolution parameters that are
highly covariant with DM. ‡The radial velocity is derived assuming no excess
noise contributes to the measured $\ddot{f}$, as discussed in Section 5.2.
Observation Properties
---
| PPTA-DR2e | PPTA-DR3
Reference epoch (MJD) | 55486 | 55486
Reference terrestrial time standard | BIPM(2018) | BIPM(2020)
Solar system ephemeris | DE436 | DE440
Observing span (yrs) | 22.016 | 14.641
Number of observations | 6396 | 3161
tempo2 Binary model | T2 | T2
Measurements
| PPTA-DR2e | PPTA-DR3
Spin frequency, $f$ (s-1) | $\mathbf{173.6879456649435(2)}$ | 173.6879456649439(4)
Spin frequency derivative, $\dot{f}$ ($10^{-15}\,$s-2) | $\mathbf{-1.7283626(13)}$ |  $-$1.728367(3)
Spin frequency second derivative, $\ddot{f}$ ($10^{-30}\,$s-3) | $\mathbf{-61(15)}$ | –
Right ascension, $\alpha$ (hh:mm:ss, J2000) |  $\mathbf{04}$:$\mathbf{37}$:$\mathbf{15.9284039(2)}$ | 04:37:15.9284042(4)
Declination, $\delta$ (dd:mm:ss, J2000) |  $\mathbf{-47}$:$\mathbf{15}$:$\mathbf{09.303707(2)}$ |  $-$47:15:09.303700(4)
Dispersion measure, DM (pc cm-3) | 2.64539(3)† | 2.64539(3)†
Proper motion in $\alpha$, $\mu_{\alpha}\cos\delta$ (mas yr-1) | $\mathbf{121.4429(5)}$ | 121.4420(6)
Proper motion in $\delta$, $\mu_{\delta}$ (mas yr-1) | $\mathbf{-71.4715(5)}$ |  $-$71.4717(7)
Parallax, $\pi$ (mas) | $\mathbf{6.43(4)}$ | 6.43(5)
Orbital period (days) | $\mathbf{5.74104582(16)}$ | 5.74104635(19)
Epoch of periastron (MJD) | $\mathbf{54530.17223(19)}$ | 54530.1722(2)
Projected semi-major axis (lt-s) | $3.36671463(3)$ | $\mathbf{3.36671466(3)}$
Longitude of periastron (deg) | $\mathbf{1.358(12)}$ | 1.359(13)
Eccentricity ($\times 10^{-5}$) | $1.91819(11)$ | $\mathbf{1.91805(10)}$
Orbital period derivative, $\dot{P_{b}}$ ($\times 10^{-12}$) | $\mathbf{3.7329(15)}$ | 3.733(3)
Periastron advance, $\dot{\omega}$ (deg yr-1) | $\mathbf{0.0134(6)}$ | 0.0155(8)
Companion mass, $M_{c}$ ($M_{\odot}$) | $0.224(4)$ | $\mathbf{0.221(4)}$
Longitude of ascending node, $\Omega$ (deg) | $208.5(9)$ | $\mathbf{208.3(8)}$
Orbital inclination, $i$ (deg) | $137.504(18)$ | $\mathbf{137.506(16)}$
Derivations
| PPTA-DR2e | PPTA-DR3
Pulsar mass (Shapiro delay), $M_{p}$ ($M_{\odot}$) | $1.446(45)$ | $\mathbf{1.418(44)}$
Pulsar distance (Shklovskii effect), $D_{\rm shk}$ (pc) | $\mathbf{156.96(11)}$ | $156.98(15)$
Pulsar distance (parallax), $D_{\pi}$ (pc) | $\mathbf{155.4(9)}$ | $155.4(12)$
Radial velocity, $v_{r}$ (km s-1) |  $\mathbf{-75(18)}$‡ | –
Figure 3: Gaussian probability densities for the pulsar mass, distance, and
orbital inclination measurements from various PPTA data sets, relative to the
broad prior adopted for NICER early analyses and model selection.

###  4.1 Precision measurements of neutron star mass, distance, and orbital
inclination

The timing model parameters of most importance for this work ($M_{c}$, $i$,
and $\dot{P}_{b}$) are measured with high significance, are uncorrelated
(between other timing model parameters and the noise processes), and are
Gaussian distributed. The NICER models use our PPTA-DR3 measurements of
$M_{p}=1.418\pm 0.044$ $M_{\odot}$, $D=156.98\pm 0.15$ pc, and $i=137.506\pm
0.016$ degrees, derived from these key parameters (with a Monte Carlo
simulation), as Gaussian priors.

For the purpose of model selection and early testing in mid-2019, the
inference of NICER observations adopted Gaussian priors (“Broad prior” in
Figure 3) that were centered near the best available measurement at the time
(from PPTA-DR2.5). The standard deviation of the priors were broader for $i$
and $D$ as these terms were thought to be most susceptible to errors in noise
models. These values are listed with measurements from all other data sets in
Table 2 and visualized in Figure 3.

The orbital inclination is assumed to be equal to the inclination of the spin
axis of the neutron star itself. MSPs are formed by mass transfer from a
companion star that works to align the spin angular momentum with the orbital
angular momentum over $\sim 100\,$Myr (Tauris & Savonije, 1999). One NICER
analysis of PSR J0740+6620 assumed a $\pm 5^{\circ}$ uniform prior around the
timing model measurement of $i$ to allow for slight deviations from this
expected alignment (Miller et al., 2021). It was determined, however, that the
data was insensitive to the precise choice of $i$ prior.

##  5 Discussion

Here we discuss the consistency of the parameter measurements from various
data sets, including comparisons with the broader literature, and our
derivation of the radial velocity.

###  5.1 Distance measurements

Knowing MSP distances to high precision is valuable to PTA searches for
continuous gravitational waves. If the distance error is much less than the
wavelength of the gravitational wave, the pulsar term can be predicted,
forming a coherent baseline for the array. We have inferred a high-precision
distance from the Shklovskii effect on the orbital period derivative, $D_{\rm
shk}=156.96\pm 0.11$ pc, after correcting for differential Galactic
acceleration. There is a $1.7\sigma$ discrepancy with the more uncertain
distance inferred from the timing parallax, $D_{\pi}=155.4\pm 0.9$ pc. The
independent astrometric distance from very-long-baseline interferometry (VLBI;
Deller et al., 2008), $D_{\rm VLBI}=156.3\pm 1.3$ pc does not resolve the
tension. A new distance constraint from an analysis of scintillation arcs,
$D_{\rm scint}=158.1\pm 1.5$ pc, also differs from the timing parallax
distance by $\sim 1.5\sigma$ and is more consistent with the $\dot{P}_{b}$
distance (D. Reardon et al. 2024, submitted).

###  5.2 Second spin frequency derivative

Using PPTA-DR2e and the noise model of Goncharov et al. (2021), we measure a
second derivative of the spin frequency, $\ddot{f}=(-61\pm 15)\times
10^{-30}\,$s-3. This parameter could be interpreted as arising from kinematic
effects (see Section 3.4), or from excess low-frequency noise that was not
characterized by the PPTA-DR2 noise models.

For the kinematic interpretation, we assume that the PPTA-DR2 noise model
remains accurate when extrapolated to apply to the extended legacy data
(therefore, $\ddot{f}_{\rm red}=0$). We derive the radial velocity of the
pulsar, $v_{r}=-75\pm 18\,$km s-1 from Equation 4. The uncertainty accounts
for covariance between $\ddot{f}$ and the Gaussian processes in the noise
model, assuming fixed amplitudes and exponents. Although the noise model is
incomplete for PPTA-DR2e (Goncharov et al., 2021), it whitens the residuals at
the longest timescales that are most relevant for obtaining an $\ddot{f}$
measurement. The noise models of Goncharov et al. (2021) were determined
simultaneously with the $\ddot{f}$ in the timing model so that the expected
signature, if present, was not absorbed by the achromatic power-law Gaussian
process.

To demonstrate the alternate $\ddot{f}$ interpretation of excess low-frequency
noise, we instead extrapolate the common-spectrum stochastic process inferred
from the 30 MSPs in PPTA-DR3 (Reardon et al., 2023a). Assuming that this
steeper-spectrum noise process (with $\gamma=-3.87$, $\log_{10}A=-14.5$)
applies to PPTA-DR2e, we measure $\ddot{f}=(-19\pm 45)\times 10^{-30}$ s-3,
and therefore, a kinematic interpretation is not required. The constraint on
radial velocity assuming this alternate noise model is $v_{r}=-24\pm 40\,$km
s-1. Under this interpretation, the measured $\ddot{f}$ from PPTA-DR2e arises
because the Goncharov et al. (2021) noise model underestimates the achromatic
noise PSD at the lowest frequencies.

A recent analysis of scintillation arcs using data from the MeerKAT radio
telescope gives an independent estimate of the pulsar radial velocity by
solving the three-dimensional geometry of the pulsar bow shock,
$v_{r}=-44.5\pm 4.3\,$km s-1 (D. Reardon et al. 2024, submitted), a difference
of $1.6\sigma$ from the value derived from the kinematic $\ddot{f}$
interpretation. If this radial velocity is accurate, then $\ddot{f}=-36.3\pm
3.5\times 10^{-30}\,$s-3 could be included in the timing model of PSR
J0437$-$4715 to improve the accuracy of recovered steep spectrum stochastic
processes, including the putative gravitational-wave background.

While the radial velocity and excess red noise are the most likely
explanations for the measured $\ddot{f}$, it is also possible that some
component is the consequence of reflex motion from a third, widely separated
body in the system. To constrain the mass and orbital separation of a third
body, we use the model presented in Rasio (1994), developed for the triple-
system globular cluster MSP PSR B1620$-$26 (Backer et al., 1993). We chose
this model because it also provides the expected secular perturbations to the
pulsar-white dwarf orbit. In this case, the third body mass $M_{3}$ and
semimajor axis $a_{3}$ can be estimated as (Rasio, 1994)

$\ddot{f}=-\frac{f}{c}\frac{G^{3/2}M_{3}}{a_{3}^{7/2}}\sqrt{M_{p}+M_{c}+M_{3}}\cos(\omega+\phi)\sin
i,$ (5)

where $\phi$ is the position phase of the outer orbit relative to $\omega$ and
the model assumes that the inner and outer orbits are coplanar. In the case
that $M_{3}\ll M_{p}+M_{c}$, we can approximate the relationship between
$M_{3}$ and $a_{3}$ to be

$\displaystyle M_{3}$ $\displaystyle\approx
0.3M_{\earth}\left(\frac{\ddot{f}}{61\times 10^{-30}\,{\rm
s}^{-3}}\right)\left(\frac{f}{173.7\,{\rm Hz}}\right)^{-1}$ (6)
$\displaystyle~{}~{}~{}\times\left(\frac{M_{\rm
tot}}{1.67M_{\odot}}\right)^{-1/2}\left(\frac{a_{3}}{100\,{\rm
AU}}\right)^{7/2},$

where $M_{\rm tot}$ is the total mass of the system and we have assumed
$\cos(\omega+\phi)=-1$. For a circular orbit with radius $100$ AU, the orbital
period would be $\approx 430$ yr.

There is additional observational evidence that an outer companion would need
to be low mass. Optical photometric and spectroscopic analysis of the
companion(s) to PSR J0437$-$4715 is consistent with a single helium white
dwarf star lacking infrared excess that would be indicative of a massive
nondegenerate companion (Durant et al., 2012). However, an $M_{3}\sim 0.03$
$M_{\odot}$ L or T dwarf companion in a $2000$ au, $40000$ yr orbit would be
less luminous in the infrared than the system white dwarf (Charnay et al.,
2018). Such a companion would induce perturbations in $\dot{\omega}$ and
$\dot{x}$ of $\mathcal{O}(10^{-11})$ deg yr-1 and $\mathcal{O}(10^{-20})$
respectively (Rasio, 1994). For any plausible $M_{3}$, these secular
perturbations are orders of magnitude smaller than the current precision of
these measurements.

###  5.3 Consistency between PPTA data sets

We find that most model parameters in PPTA-DR2e and PPTA-DR3 are consistent
within (1$\sigma$) uncertainties, with the exception of spin frequency $f$
(1.3$\sigma$), declination $\delta$ (1.5$\sigma$), proper motion in right
ascension $\mu_{\alpha}$ (1.1$\sigma$), and the highly covariant parameters
$P_{b}$ and $\dot{\omega}$ (2.1$\sigma$). The discrepancy of $P_{b}$ and
$\dot{\omega}$ likely originates from different noise model choices but has no
significant effect on the other parameters of the model. The $\dot{\omega}$
value from PPTA-DR2e is more precise and more consistent with the expectation
from general relativity than PPTA-DR3 (assuming the component masses from
PPTA-DR3 and no contamination from kinematic effects as secular variations are
accounted for with the Kopeikin terms). However, for these analyses, the
$\dot{\omega}$ was not used to constrain the pulsar mass because the
uncertainty is too large.

The Shapiro delay and annual-orbital parallax parameters are more precise with
PPTA-DR3 because of the improved data quality. These effects gives the orbital
inclination and the longitude of ascending node $\Omega$, which have also been
measured using scintillation arcs in the PPTA-DR2 data set (Reardon et al.,
2020). The value of $\Omega$ derived from these scintillation arcs differs by
$2.2\sigma$, which is most likely to be explained by weakness in the
scintillation arc curvature model, which assumes a stable interstellar
scattering screen over many years.

Table 2: Inferred pulsar mass, distance, and orbital inclination for PSR J0437$-$4715, using various PPTA data sets. Data set | Pulsar mass, $M_{p}$ ($M_{\odot}$) | Distance, $D$ (pc) | Orbital inclination, $i$ (deg)
---|---|---|---
PPTA-DR1e† | $1.44\pm 0.07$ | $156.79\pm 0.25$ | $137.56\pm 0.04$
PPTA-DR2 | $1.404\pm 0.050$ | $156.78\pm 0.20$ | $137.490\pm 0.016$
PPTA-DR2e | $1.446\pm 0.045$ | $156.96\pm 0.11$ | $137.504\pm 0.018$
PPTA-DR2.5 | $1.411\pm 0.056$ | $156.69\pm 0.30$ | $137.485\pm 0.038$
PPTA-DR3 | $1.418\pm 0.044$ | $156.98\pm 0.15$ | $137.506\pm 0.016$
Broad NICER prior | $1.41\pm 0.07$ | $157.0\pm 1.0$ | $137.50\pm 0.10$
† from Reardon et al. (2016) |  |  |

##  6 Conclusion

The precise ($\sim 3\%$ uncertainty) pulsar mass for PSR J0437$-$4715 paves
the way for a stringent constraint on the neutron star EoS. We find that the
pulsar has a mass, $M_{p}=1.418\pm 0.044$ $M_{\odot}$, near the center of the
population of neutron star masses (Özel & Freire, 2016). A precise radius
measurement from NICER, enabled by the radio-frequency timing mass, imposes a
tight restriction on the allowed mass-radius relationship predicted by the
EoS, particularly for this central mass range (Choudhury et al., 2024;
Rutherford et al., 2024).

In our attempt to produce the most accurate pulsar mass measurement possible,
we developed the first data sets with a complete noise model for PSR
J0437$-$4715 (PPTA-DR2.5 and PPTA-DR3). We compared the timing model
parameters from multiple generations of PPTA data sets (with different noise
models), and found that the parameters important for characterizing the
physical properties of the system do not significantly vary. This is
encouraging because data sets containing observations prior to 2007 are known
to contain excess unmodeled noise.

Using the longest data set, PPTA-DR2e, we measure a significant second
derivative of the pulsar spin frequency, $\ddot{f}$. We discussed two likely
interpretations of this measurement, in the form of kinematic effects due to a
pulsar radial velocity of $v_{r}=75\pm 18$ km s-1, or excess low-frequency
noise such as from a gravitational-wave background. Additionally, we derived a
constraint on the allowed mass and orbital separation of any third body in the
case that some component of $\ddot{f}$ is due to reflex motion. Future
analyses may consider using a modified Fourier basis (Keith & Niţu, 2023) to
ensure robust radial velocity measurements, particularly if a gravitational-
wave background begins to dominate the achromatic noise in longer data sets.

Future PPTA data sets could also be improved by employing new techniques such
as matrix equation template matching for instrumental calibration, matrix
template matching for ToA measurements (van Straten, 2013; Rogers et al.,
2024), or wide-band timing (Curyło et al., 2023).

## Acknowledgements

Murriyang, the Parkes radio telescope, is part of the Australia Telescope
National Facility (https://ror.org/05qajvd42) which is funded by the
Australian Government for operation as a National Facility managed by CSIRO.
We acknowledge the Wiradjuri people as the Traditional Owners of the
Observatory site. Parts of this research were conducted by the Australian
Research Council Centre of Excellence for Gravitational Wave Discovery
(OzGrav), through project numbers CE170100004 and CE230100016. R.M.S.
acknowledges support through Australian Research Council Future Fellowship
FT190100155. Work at NRL is supported by NASA. L.Z. is supported by ACAMAR
Postdoctoral Fellowship and the National Natural Science Foundation of China
(grant No. 12103069).

## References

  * Agazie et al. (2023a) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJ, 951, L8, doi: 10.3847/2041-8213/acdac6
  * Agazie et al. (2023b) —. 2023b, ApJ, 951, L10, doi: 10.3847/2041-8213/acda88
  * Antoniadis et al. (2023a) Antoniadis, J., Arumugam, P., Arumugam, S., et al. 2023a, A&A, 678, A50, doi: 10.1051/0004-6361/202346844
  * Antoniadis et al. (2023b) —. 2023b, A&A, 678, A49, doi: 10.1051/0004-6361/202346842
  * Arzoumanian et al. (2018) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018, ApJS, 235, 37, doi: 10.3847/1538-4365/aab5b0
  * Backer et al. (1993) Backer, D. C., Foster, R. S., & Sallmen, S. 1993, Nature, 365, 817, doi: 10.1038/365817a0
  * Bassa et al. (2016) Bassa, C. G., Janssen, G. H., Stappers, B. W., et al. 2016, MNRAS, 460, 2207, doi: 10.1093/mnras/stw1134
  * Bell & Bailes (1996) Bell, J. F., & Bailes, M. 1996, ApJ, 456, L33, doi: 10.1086/309862
  * Bell et al. (1993) Bell, J. F., Bailes, M., & Bessell, M. S. 1993, Nature, 364, 603, doi: 10.1038/364603a0
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29, doi: 10.1088/0067-0049/216/2/29
  * Bovy et al. (2012) Bovy, J., Allende Prieto, C., Beers, T. C., et al. 2012, ApJ, 759, 131, doi: 10.1088/0004-637X/759/2/131
  * Britton (2000) Britton, M. C. 2000, ApJ, 532, 1240, doi: 10.1086/308595
  * Charnay et al. (2018) Charnay, B., Bézard, B., Baudino, J. L., et al. 2018, ApJ, 854, 172, doi: 10.3847/1538-4357/aaac7d
  * Choudhury et al. (2024) Choudhury, D., Salmi, T., Vinciguerra, S., et al. 2024, arXiv e-prints, arXiv:2407.06789, doi: 10.48550/arXiv.2407.06789
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, Nature Astronomy, 4, 72, doi: 10.1038/s41550-019-0880-2
  * Curyło et al. (2023) Curyło, M., Pennucci, T. T., Bailes, M., et al. 2023, ApJ, 944, 128, doi: 10.3847/1538-4357/aca535
  * Damour & Deruelle (1986) Damour, T., & Deruelle, N. 1986, Ann. Inst. Henri Poincaré Phys. Théor, 44, 263
  * Damour & Taylor (1992) Damour, T., & Taylor, J. H. 1992, Phys. Rev. D, 45, 1840, doi: 10.1103/PhysRevD.45.1840
  * Deller et al. (2008) Deller, A. T., Verbiest, J. P. W., Tingay, S. J., & Bailes, M. 2008, ApJ, 685, L67, doi: 10.1086/592401
  * Demorest et al. (2010) Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081, doi: 10.1038/nature09466
  * Donner et al. (2020) Donner, J. Y., Verbiest, J. P. W., Tiburzi, C., et al. 2020, A&A, 644, A153, doi: 10.1051/0004-6361/202039517
  * Durant et al. (2012) Durant, M., Kargaltsev, O., Pavlov, G. G., et al. 2012, ApJ, 746, 6, doi: 10.1088/0004-637X/746/1/6
  * Edwards et al. (2006) Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549, doi: 10.1111/j.1365-2966.2006.10870.x
  * Ellis et al. (2019) Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2019, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE, Astrophysics Source Code Library, record ascl:1912.015.  http://ascl.net/1912.015
  * Folkner et al. (2014) Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. 2014, Interplanetary Network Progress Report, 42-196, 1
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, ApJ, 915, L12, doi: 10.3847/2041-8213/ac03b8
  * Gendreau et al. (2012) Gendreau, K. C., Arzoumanian, Z., & Okajima, T. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8443, Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, ed. T. Takahashi, S. S. Murray, & J.-W. A. den Herder, 844313, doi: 10.1117/12.926396
  * Gendreau et al. (2016) Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, T. Takahashi, & M. Bautz, 99051H, doi: 10.1117/12.2231304
  * Goncharov et al. (2021) Goncharov, B., Reardon, D. J., Shannon, R. M., et al. 2021, MNRAS, 502, 478, doi: 10.1093/mnras/staa3411
  * Hobbs et al. (2011) Hobbs, G., Miller, D., Manchester, R. N., et al. 2011, PASA, 28, 202, doi: 10.1071/AS11016
  * Hobbs et al. (2020) Hobbs, G., Manchester, R. N., Dunning, A., et al. 2020, PASA, 37, e012, doi: 10.1017/pasa.2020.2
  * Hobbs et al. (2006) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655, doi: 10.1111/j.1365-2966.2006.10302.x
  * Holmberg & Flynn (2004) Holmberg, J., & Flynn, C. 2004, MNRAS, 352, 440, doi: 10.1111/j.1365-2966.2004.07931.x
  * Jennings et al. (2024) Jennings, R. J., Cordes, J. M., Chatterjee, S., et al. 2024, ApJ, 964, 179, doi: 10.3847/1538-4357/ad2930
  * Johnston et al. (1993) Johnston, S., Lorimer, D. R., Harrison, P. A., et al. 1993, Nature, 361, 613, doi: 10.1038/361613a0
  * Kaplan et al. (2016) Kaplan, D. L., Kupfer, T., Nice, D. J., et al. 2016, ApJ, 826, 86, doi: 10.3847/0004-637X/826/1/86
  * Keith & Niţu (2023) Keith, M. J., & Niţu, I. C. 2023, MNRAS, 523, 4603, doi: 10.1093/mnras/stad1713
  * Keith et al. (2013) Keith, M. J., Coles, W., Shannon, R. M., et al. 2013, MNRAS, 429, 2161, doi: 10.1093/mnras/sts486
  * Kerr et al. (2020) Kerr, M., Reardon, D. J., Hobbs, G., et al. 2020, PASA, 37, e020, doi: 10.1017/pasa.2020.11
  * Kopeikin (1995) Kopeikin, S. M. 1995, ApJ, 439, L5, doi: 10.1086/187731
  * Kopeikin (1996) —. 1996, ApJ, 467, L93, doi: 10.1086/310201
  * Lam et al. (2019) Lam, M. T., McLaughlin, M. A., Arzoumanian, Z., et al. 2019, ApJ, 872, 193, doi: 10.3847/1538-4357/ab01cd
  * Lentati et al. (2014a) Lentati, L., Alexander, P., Hobson, M. P., et al. 2014a, MNRAS, 437, 3004, doi: 10.1093/mnras/stt2122
  * Lentati et al. (2013) —. 2013, Phys. Rev. D, 87, 104021, doi: 10.1103/PhysRevD.87.104021
  * Lentati et al. (2014b) Lentati, L., Hobson, M. P., & Alexander, P. 2014b, MNRAS, 444, 3863, doi: 10.1093/mnras/stu1721
  * Lentati et al. (2016) Lentati, L., Shannon, R. M., Coles, W. A., et al. 2016, MNRAS, 458, 2161, doi: 10.1093/mnras/stw395
  * Liu et al. (2012) Liu, K., Keane, E. F., Lee, K. J., et al. 2012, MNRAS, 420, 361, doi: 10.1111/j.1365-2966.2011.20041.x
  * Liu et al. (2018) Liu, X. J., Bassa, C. G., & Stappers, B. W. 2018, MNRAS, 478, 2359, doi: 10.1093/mnras/sty1202
  * Liu et al. (2019) Liu, X. J., Keith, M. J., Bassa, C. G., & Stappers, B. W. 2019, MNRAS, 488, 2190, doi: 10.1093/mnras/stz1801
  * Manchester et al. (2013) Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, e017, doi: 10.1017/pasa.2012.017
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24, doi: 10.3847/2041-8213/ab50c5
  * Miller et al. (2021) —. 2021, ApJ, 918, L28, doi: 10.3847/2041-8213/ac089b
  * Özel & Freire (2016) Özel, F., & Freire, P. 2016, ARA&A, 54, 401, doi: 10.1146/annurev-astro-081915-023322
  * Özel et al. (2016) Özel, F., Psaltis, D., Güver, T., et al. 2016, ApJ, 820, 28, doi: 10.3847/0004-637X/820/1/28
  * Park et al. (2021) Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. 2021, AJ, 161, 105, doi: 10.3847/1538-3881/abd414
  * Parthasarathy et al. (2021) Parthasarathy, A., Bailes, M., Shannon, R. M., et al. 2021, MNRAS, 502, 407, doi: 10.1093/mnras/stab037
  * Rasio (1994) Rasio, F. A. 1994, ApJ, 427, L107, doi: 10.1086/187376
  * Reardon et al. (2016) Reardon, D. J., Hobbs, G., Coles, W., et al. 2016, MNRAS, 455, 1751, doi: 10.1093/mnras/stv2395
  * Reardon et al. (2020) Reardon, D. J., Coles, W. A., Bailes, M., et al. 2020, ApJ, 904, 104, doi: 10.3847/1538-4357/abbd40
  * Reardon et al. (2021) Reardon, D. J., Shannon, R. M., Cameron, A. D., et al. 2021, MNRAS, 507, 2137, doi: 10.1093/mnras/stab1990
  * Reardon et al. (2023a) Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023a, ApJ, 951, L6, doi: 10.3847/2041-8213/acdd02
  * Reardon et al. (2023b) —. 2023b, ApJ, 951, L7, doi: 10.3847/2041-8213/acdd03
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21, doi: 10.3847/2041-8213/ab481c
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJ, 918, L27, doi: 10.3847/2041-8213/ac0a81
  * Rogers et al. (2024) Rogers, A. F., van Straten, W., Gulyaev, S., et al. 2024, arXiv e-prints, arXiv:2407.20015, doi: 10.48550/arXiv.2407.20015
  * Rutherford et al. (2024) Rutherford, N., Mendes, M., Svensson, I., et al. 2024, arXiv e-prints, arXiv:2407.06790, doi: 10.48550/arXiv.2407.06790
  * Shannon & Cordes (2010) Shannon, R. M., & Cordes, J. M. 2010, ApJ, 725, 1607, doi: 10.1088/0004-637X/725/2/1607
  * Shannon et al. (2014) Shannon, R. M., Osłowski, S., Dai, S., et al. 2014, MNRAS, 443, 1463, doi: 10.1093/mnras/stu1213
  * Shannon et al. (2016) Shannon, R. M., Lentati, L. T., Kerr, M., et al. 2016, ApJ, 828, L1, doi: 10.3847/2041-8205/828/1/L1
  * Shapiro (1964) Shapiro, I. I. 1964, Phys. Rev. Lett., 13, 789, doi: 10.1103/PhysRevLett.13.789
  * Shklovskii (1970) Shklovskii, I. S. 1970, Soviet Ast., 13, 562
  * Tauris & Savonije (1999) Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928, doi: 10.48550/arXiv.astro-ph/9909147
  * Turner et al. (2021) Turner, J. E., McLaughlin, M. A., Cordes, J. M., et al. 2021, ApJ, 917, 10, doi: 10.3847/1538-4357/abfafe
  * van Haasteren & Levin (2013) van Haasteren, R., & Levin, Y. 2013, MNRAS, 428, 1147, doi: 10.1093/mnras/sts097
  * van Straten (2003) van Straten, W. 2003, PhD thesis, Swinburne University of Technology, Australia
  * van Straten (2013) —. 2013, ApJS, 204, 13, doi: 10.1088/0067-0049/204/1/13
  * van Straten et al. (2012) van Straten, W., Demorest, P., & Oslowski, S. 2012, Astronomical Research and Technology, 9, 237, doi: 10.48550/arXiv.1205.6276
  * Verbiest et al. (2008) Verbiest, J. P. W., Bailes, M., van Straten, W., et al. 2008, ApJ, 679, 675, doi: 10.1086/529576
  * Xu et al. (2023) Xu, H., Chen, S., Guo, Y., et al. 2023, Research in Astronomy and Astrophysics, 23, 075024, doi: 10.1088/1674-4527/acdfa5
  * Zic et al. (2023) Zic, A., Reardon, D. J., Kapur, A., et al. 2023, PASA, 40, e049, doi: 10.1017/pasa.2023.36