# High-cadence Timing of Binary Pulsars with CHIME

Chia Min Tan Department of Physics, McGill University, 3600 rue University,
Montréal, QC H3A 2T8, Canada  The Trottier Space Institute at McGill, 3550 rue
University, Montréal, QC H3A 2A7, Canada  International Centre for Radio
Astronomy Research, Curtin University, Bentley, WA 6102, Australia  Emmanuel
Fonseca Department of Physics and Astronomy, West Virginia University,
Morgantown, WV 26506-6315, USA  Center for Gravitational Waves and Cosmology,
Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Kathryn Crowter
Department of Physics and Astronomy, University of British Columbia, 6224
Agricultural Road, Vancouver, BC V6T 1Z1, Canada  Fengqiu Adam Dong Department
of Physics and Astronomy, University of British Columbia, 6224 Agricultural
Road, Vancouver, BC V6T 1Z1, Canada  Victoria M. Kaspi Department of Physics,
McGill University, 3600 rue University, Montréal, QC H3A 2T8, Canada  The
Trottier Space Institute at McGill, 3550 rue University, Montréal, QC H3A 2A7,
Canada  Kiyoshi W. Masui MIT Kavli Institute for Astrophysics and Space
Research, Massachusetts Institute of Technology, 77 Massachusetts Ave,
Cambridge, MA 02139, USA  Department of Physics, Massachusetts Institute of
Technology, 77 Massachusetts Ave, Cambridge, MA 02139, USA  James W. McKee E.
A. Milne Centre for Astrophysics, University of Hull, Cottingham Road,
Kingston-upon-Hull, HU6 7RX, UK  Centre of Excellence for Data Science,
Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham
Road, Kingston-upon-Hull, HU6 7RX, UK  Bradley W. Meyers International Centre
for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia
Scott M. Ransom National Radio Astronomy Observatory, 520 Edgemont Rd,
Charlottesville, VA 22903, USA  Ingrid H. Stairs Department of Physics and
Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver,
BC V6T 1Z1, Canada

(Received XXX; Revised YYY; Accepted ZZZ)

###### Abstract

We performed near-daily observations on the binary pulsars PSR J0218+4232, PSR
J1518+4904 and PSR J2023+2853 with the Canadian Hydrogen Intensity Mapping
Experiment (CHIME). For the first time, we detected the Shapiro time delay in
all three pulsar-binary systems, using only 2–4 years of CHIME/Pulsar timing
data. We measured the pulsar masses to be $1.49^{+0.23}_{-0.20}$ M⊙,
$1.470^{+0.030}_{-0.034}$ M⊙ and $1.50^{+0.49}_{-0.38}$ M⊙ respectively. The
companion mass to PSR J0218+4232 was found to be $0.179^{+0.018}_{-0.016}$ M⊙.
We constrained the mass of the neutron-star companion of PSR J1518+4904 to be
$1.248^{+0.035}_{-0.029}$ M⊙, using the observed apsidal motion as a
constraint on mass estimation. The binary companion to PSR J2023+2853 was
found to have a mass of $0.93^{+0.17}_{-0.14}$ M⊙; in the context of the near-
circular orbit, this mass estimate suggests that the companion to PSR
J2023+2853 is likely a high-mass white dwarf. By comparing the timing model
obtained for PSR J0218+4232 with previous observations, we found a significant
change in the observed orbital period of the system of $\dot{P_{\rm
b}}=0.14(2)\times 10^{-12}$; we determined that this variation arises from
“Shklovskii acceleration” due to relative motion of the binary system, and
used this measurement to estimate a distance of $d=(6.7\pm 1.0)$ kpc to PSR
J0218+4232. This work demonstrates the capability of high-cadence
observations, enabled by the CHIME/Pulsar system, to detect and refine
general-relativistic effects of binary pulsars over short observing
timescales.

instrumentation: interferometers – methods: observational – radio continuum:
general – pulsars: general – techniques: interferometric – telescopes

††software:  DSPSR (van Straten & Bailes, 2011), PSRDADA
(http://psrdada.sourceforge.net), PSRCHIVE (Hotan et al., 2004; van Straten et
al., 2012), PRESTO (https://github.com/scottransom/presto), TEMPO
(https://tempo.sourceforge.net), TEMPO2 (Hobbs et al., 2006)

##  1 Introduction

Timing studies of radio pulsars typically depend on the overall time span of
observations to measure parameters that quantify apparent or intrinsic spin
variations. Until recently, the cadence of pulsar-timing observations was
limited to $\sim$monthly observing rates due to competitive telescope
resources. The recent rise of high-cadence pulsar observing programs, such as
the UTMOST (Jankowski et al., 2019; Lower et al., 2020) and CHIME/Pulsar
(CHIME/Pulsar Collaboration et al., 2021) experiments, are producing uniquely
dense data sets for a majority of the known pulsar population.

The CHIME/Pulsar system observes up to 10 sources simultaneously at any given
time, allowing for near-daily observations of a selected set of pulsars. This
circumstance has resulted in several studies that leverage high-cadence
observations to make unique measurements on a growing number of sources. For
example, high-cadence observations of the binary PSR J2108+4516 reveal
extreme, intra-day variations in the electromagnetic properties of the
broadband pulsar signal due to a turbulent wind and/or circumstellar disk of
its high-mass companion star (Andersen et al., 2023). Moreover, high-cadence
CHIME/Pulsar data recently contributed to the refined timing analysis of the
binary PSR J0740+6620 for improved mass constraints (Fonseca et al., 2021),
which directly led to the estimation of its radius using X-ray data obtained
by the NICER telescope (Riley et al., 2021; Miller et al., 2021). High-cadence
CHIME/Pulsar observations are supporting the analysis of pulsars discovered by
the PALFA (Parent et al., 2022) and GBNCC (Stovall et al., 2014) surveys, as
well as contributing to the data set developed by NANOGrav for improved
detection of the nHz-frequency background of gravitational radiation (e.g.,
Lam, 2018). In all of these aforementioned works, the high-cadence nature of
CHIME/Pulsar data affords new views into timing-based phenomena over observing
timescales that cannot be achieved with other Northern-hemisphere instruments.

In this study, we perform high-precision timing analyses of three binary
pulsars and explore the impact of high-cadence observations obtained with the
CHIME/Pulsar backend. Two of these pulsars – PSRs J0218+4232 and J1518+4904 –
have long been known and well-studied using archival data that span decades,
and thus allow for independent checks of model accuracy. The remaining source,
PSR J2023+2853, is a binary pulsar recently discovered in a survey of the
Galactic plane using the Five-hundred Meter Aperture Synthesis Telescope
(FAST; Han et al., 2021); to our knowledge, our study below presents the first
coherent timing solution of PSR J2023+2853 and demonstrates unambiguous period
variations due to orbital motion. In Section 2, we outline the logistics of
our observing program and processing of data. In Section 3, we outline the
procedures used to construct timing models for all three binary pulsars. In
Section 4, we present the results obtained from high-precision timing
analyses.

##  2 Observations & Reduction

All data presented in this work were acquired with the CHIME telescope, a
static interferometer that digitizes the sky as it drifts within its field of
view (CHIME Collaboration et al., 2022). The CHIME telescope digitizes raw
voltages from 1,024 dual-polarization antennae, sensitive to the 400–800 MHz
range, that are mounted on four half-cylinder reflectors that span a $80\times
100\ m^{2}$ area; these digitized samples are then beamformed within the “FX”
correlator to yield up to 10 baseband streams. Beamformed baseband data are
phased to user-specified celestial positions in order to track the known radio
sources.

These baseband streams are transmitted to the CHIME/Pulsar backend, a ten-
node, real-time computing cluster that enables CHIME to make measurements of
radio pulsars and other radio-transient phenomena under user-specified
acquisition modes (CHIME/Pulsar Collaboration et al., 2021). Each CHIME/Pulsar
node receives one baseband stream, complex-sampled at a rate of 2.56 $\mu$s
and for 1,024 frequency channels that span the 400–800 MHz range.

Timing data for our study were obtained using the CHIME/Pulsar backend in its
“fold” mode. For all fold-mode observations, input voltages are coherently
dedispersed (Hankins & Rickett, 1975), and then folded using an existing
timing model to average individual pulses into 10-second “subintegrations”.
Coherent dedispersion is achieved by de-convolving the ISM impulse response
function in the frequency domain from the input baseband prior to detection
and folding of the radio signal, which removes intra-channel smearing due to
pulse dispersion (e.g., Hankins & Rickett, 1975). We used the DSPSR software
suite (van Straten & Bailes, 2011) for real-time coherent dedispersion and
folding with a graphics processing unit mounted onto each CHIME/Pulsar node.

We initially took 10 observations of PSR J2023+2853 in filterbank mode with
the CHIME/Pulsar backend, as no timing model that incorporates the binary
orbital parameter existed at the time we began conducting our measurements. We
folded the filterbank data to the spin period and dispersion measure presented
in Han et al. (2021) and found that the apparent spin period of the pulsar
changes between observations, indicating a potential binary companion. We
modelled the variation in the apparent spin period assuming that they arose
from Doppler variations due to unmodeled orbital motion, and obtained an
initial binary solution that describes a near-circular orbit with binary
period of 0.7 days and projected semi-major axis of 4 lt-s. We then continued
to observe PSR J2023+2853 in “fold” mode using a timing model that includes
these orbital parameters.

The resultant fold-mode data yield time-averaged profiles for all four
components of the Stokes polarization vector measured over 10-s
subintegrations, a user-defined number of profile bins ($n_{\rm bin}$), and
1,024 frequency channels. We set the profile resolution to be $n_{\rm
bin}=512$ for PSR J0218+4232, $n_{\rm bin}=1,024$ for PSR J1518+4904 and
$n_{\rm bin}=256$ for PSR J2023+2853. In the case of PSR J0218+4232, its 2.3
ms spin period and the fixed 2.56-$\mu$s baseband sample rate jointly lead an
upper limit in allowed profile resolution, such that the largest binary value
that can be used to evaluate Stokes profiles for fold-mode observations is
512. The profile resolution of PSR J2023+2853 is set at $n_{\rm bin}=256$,
which is the default number set for any CHIME/Pulsar observation. By the time
we decided to include PSR J2023+2853 in this work, there were approximately 1
year of observations. Hence we decided to not change $n_{\rm bin}$ to keep the
dataset consistent.

###  2.1 Offline Excision and Reduction

During post-processing, we used the PSRCHIVE package (van Straten et al.,
2012) and related utilities for offline cleaning and downsampling of raw
timing data. We removed frequency channels that contain persistent radio-
frequency interference (RFI) using the fold mode data cleaning tool clfd
(Morello et al., 2019), combined with a mask of frequency channels that
produce unusable data for the particular day. These masks arise due to
variable unavailability of a small fraction of FX-correlator nodes that each
produce channelized timeseries for four frequency channels.

We further reduced the RFI-excised timing data by downsampling all spectra
from 1024 frequency sub-bands to a smaller number, depending on the S/N of the
pulsar, and fully integrating over time. We downsampled the number of
frequency sub-bands to $n_{\rm sub}=32$ for PSR J0218+4232, $n_{\rm sub}=8$
for PSR J1518+4904, and $n_{\rm sub}=1$ for PSR J2023+2853. The frequency-
resolved spectra of PSR J0218+4232 and PSR J1518+4904 provide opportunities
for the modeling of dispersion variations, which we enact and discuss below.

###  2.2 Generation and Pruning of TOAs

We used the common technique of cross-correlating Stokes-$I$ profiles with a
noiseless “template” profile in the Fourier domain for determining pulse times
of arrivals (TOAs) with high precision (Taylor, 1992). Standard templates were
generated by aligning and averaging at least 100 days of data using an
adequate a timing solution to form a single, representative profile for each
pulsar considered in this work. These profiles were then fully averaged over
all subintegrations and frequency channels, and de-noised using a wavelet-
smoothing algorithm available in PSRCHIVE, psrsmooth to produce the standard
templates.

As a final step in data preparation, we cleaned all TOAs data sets using two
common methods. First, we applied a signal-to-noise (S/N) threshold to each
TOA, such that a TOA derived from a Stokes-I profile with S/N $<$ 4 was
discarded from the timing analysis presented below. Any remaining TOAs that
significantly deviated from timing-model predictions – due to RFI and/or
instrumental corruption whose statistics bypassed the S/N thresholds and RFI-
cleaining algorithms – were manually removed. These TOA-cleaning steps removed
6.8%, 3.4%, and 10.7% of TOAs from the data sets for PSR J0218+4232, PSR
J1518+4904, and PSR J2023+2853, respectively.

##  3 Timing Methods

An updated timing model based on CHIME/Pulsar data was produced for each
pulsar through the weighted least-squared fitting of the TOAs produced (See
Section 2.2) using the TEMPO111https://tempo.sourceforge.net software package.
For all timing models, we explored the significance of various astrophysical
quantities: astrometric positions and their proper motion; spin frequency and
its first derivative; the Keplerian binary parameters and their first-order
variations; and the Shapiro delay parameters. The timing models for PSR
J1518+4904 and PSR J0218+4232 include a temporal fit of the variation in the
electromagnetic dispersion measure (DM) across the timespan. This method is
not used for PSR J2023+2853 as we only obtained frequency-averaged TOAs that
do not allow for DM modelling over time, primarily due to the faintness of the
detection on individual observation epochs.

We also estimated ad hoc adjustment factors for our raw TOA uncertainties that
produced reduced goodness-of-fit ($\chi^{2}$) statistics of $\sim$ 1\. These
factors, commonly referred to as “EFAC” in pulsar-timing literature,
characterize the amount of increase in uncertainties needed for TOA residuals
to exhibit Gaussian noise largely due to imperfect estimation of the profile
template shape (e.g., Agazie et al., 2023). We found that EFAC values for our
data sets varied from 1.2–1.25, indicating that the TOA-estimation algorithm
described in Section 2 largely determined the arrival times and their
uncertainties in a robust manner. We did not further explore other separable
sources of TOA noise (e.g., “ECORR”) for our sources as such computational
efforts are complicated by the large volume of the CHIME/Pulsar data sets. As
shown in Fonseca et al. (2021), the use of “wideband” TOA algorithms (e.g.,
Pennucci et al., 2014) can reduce the TOA volume by over an order of magnitude
and will therefore be crucial for an efficient analysis of noise properties,
though we reserve such an effort for future work.

###  3.1 Models for Orbital Motion

We used the ELL1 model (Lange et al., 2001) to describe the near-circular
orbital motion of PSR J0218+4232 and PSR J2023+2853. For PSR J1518+4904, we
used the DD and DDGR models (Damour & Deruelle, 1985, 1986) to parameterize
periodic and secular orbital variations due to strong-field gravitation in
terms of the mass of the companion star ($m_{\rm c}$) and the total system
mass ($m_{\rm tot}=m_{\rm p}+m_{\rm c}$, where $m_{\rm p}$ is the mass of the
pulsar).

Both the DDGR and ELL1 models characterize the Keplerian orbits in terms of
the orbital period ($P_{\rm b}$) and radial projection of the semi-major axis
along the line of sight ($x=a\sin i/c$, where $i$ is the orbital inclination
angle), though differ in their treatment of the orbital shape; the DDGR model
uses orbital eccentricity ($e$), the argument of periastron ($\omega$), and
epoch of periastron passage ($T_{0}$) to describe ellipticity, while the ELL1
model instead uses the “Laplace-Lagrange” parameters ($\epsilon_{1}$,
$\epsilon_{2}$) and the epoch of passage through the ascending-node longitude
($T_{\rm asc}$) to describe O($e$) deviations from circular orbits. The ELL1
eccentricity parameters are related in the following manner:

$\displaystyle e$ $\displaystyle=\sqrt{\epsilon_{1}^{2}+\epsilon_{2}^{2}},$
(1) $\displaystyle\omega$
$\displaystyle=\tan^{-1}\bigg{(}\frac{\epsilon_{1}}{\epsilon_{2}}\bigg{)},$
(2) $\displaystyle T_{0}$ $\displaystyle=T_{\rm
asc}+\bigg{(}\frac{\omega}{2\pi}\bigg{)}P_{\rm b}.$ (3)

In practice, the ELL1 model is most applicable for binary pulsars with timing
precision and orbital parameters that mutually satisfy the criterion
$\sigma_{\rm TOA,rms}>xe^{2}$, where $xe^{2}$ is the amplitude of the
O($e^{2}$) correction to the ELL1 model.

Some binary pulsars eventually exhibit deviations from Keplerian motion due to
various particularly extreme properties. Given the short timescale of our
observations, the most likely source of “post-Keplerian” (PK) variations is
strong-field gravitation, which gives rise to: apsidal motion
($\dot{\omega}$); orbital decay ($\dot{P}_{\rm b}$); time dilation and
gravitational redshift ($\gamma$); and the “range” ($r$) and “shape” ($s$) of
the Shapiro delay. According to general relativity, these PK effects are
functions of the mass and geometric information of the binary system (e.g.,
Damour & Taylor, 1992):

$\displaystyle\dot{\omega}={}$ $\displaystyle 3\bigg{(}\frac{P_{\rm
b}}{2\pi}\bigg{)}^{-5/3}\frac{(T_{\odot}m_{\rm tot})^{2/3}}{1-e^{2}},$ (4)
$\displaystyle\dot{P}_{\rm b}={}$
$\displaystyle-\frac{192\pi}{5}\bigg{(}\frac{P_{\rm
b}}{2\pi}\bigg{)}^{-5/3}\frac{T_{\odot}^{5/3}m_{\rm p}m_{\rm c}}{m_{\rm
tot}^{1/3}(1-e^{2})^{7/2}}$
$\displaystyle\times\bigg{(}1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\bigg{)},$
(5) $\displaystyle\gamma={}$ $\displaystyle e\bigg{(}\frac{P_{\rm
b}}{2\pi}\bigg{)}^{1/3}T_{\odot}^{2/3}\frac{m_{\rm c}(m_{\rm tot}+m_{\rm
c})}{m_{\rm tot}^{4/3}},$ (6) $\displaystyle r={}$ $\displaystyle
T_{\odot}m_{\rm c},$ (7) $\displaystyle s\equiv{}$ $\displaystyle\sin
i=x\bigg{(}\frac{P_{\rm b}}{2\pi}\bigg{)}^{-2/3}\frac{m_{\rm
tot}^{2/3}}{T_{\odot}^{1/3}m_{\rm c}},$ (8)

where $T_{\odot}=GM_{\odot}/c^{3}=4.925490947$ $\mu$s converts the masses to
possess units of solar mass. The DD model directly measures the left-hand side
of Equations 4–8 in a theory-independent manner, while the DDGR model assumes
Equations 4–8 to directly measure $m_{\rm c}$ and $m_{\rm tot}$.

###  3.2 Methods for Modeling ISM Variations

We modeled the DM variation of each pulsar using the DMX routine of TEMPO,
which estimates a “local” value of DM within contiguous, non-overlapping time
bins that span each data set. The measured DMs over different epochs are then
subtracted from the residuals to model the other parameters. For PSR
J0218+4232, we set the width of each DMX time bin to be no larger than 5 days;
for PSR J1518+4904, we set the width of DMX time bin to be no larger than 7
days; and for PSR J2023+2853, we set the DMX width to be no larger than 10
days. These widths were selected based on the variations in the observed
trends, but were otherwise chosen arbitrarily.

We did not explicitly model or estimate the presence of pulse scatter-
broadening in any of the three pulsars subject to our study. This choice was
made for two reasons: the use of TEMPO and channelized TOAs does not allow for
simultaneous estimation of DM and scattering-timescale parameters; and there
is a lack of a scatter-broadening signature in their CHIME/Pulsar data sets
that is distinct from intrinsic profile evolution across the CHIME band. We
nonetheless implicitly modeled the composite effects of scattering-broadening
and intrinsic profile evolution by fitting for “frequency-dependent” (FD)
variations in our channelized TOAs. These variations were modeled using the FD
time delay $\Delta t_{\rm FD}=\sum_{i}c_{i}\ln^{i}\nu$, where $c_{i}$ is the
FD coefficient and $\nu$ is the electromagnetic frequency (NANOGrav
Collaboration et al., 2015).

When simultaneously estimated with all other timing parameters with TEMPO, we
found that at least one FD coefficient was sufficient for PSR J0218+4232 and
the other two pulsars yielded insignificant FD parameters. Only two FD
coefficients appear to be statistically significant for PSR J0218+4232, but
the use of two coefficients leads to clear degeneracy with the DMX parameters.
This degeneracy likely arises due to several eras of intra-day variations in
DM not being robustly modeling using our 5-day DMX bin widths. We therefore
used only one FD coefficient in our final timing model for PSR J0218+4232, and
manually estimated an EFAC using this “DMX+FD” timing model in order to
correct our TOA uncertainties. The resultant EFAC value estimated from this
model likely includes systematic biases from DM mis-estimation and any
measurable variations in scatter-broadening, both of which cannot be modeled
by the $c_{i}$ parameters as they are presumed to be constant across each data
set. However, the EFAC for PSR J0218+4232 is only $\sim$ 25% larger than
unity, indicating that the ISM variations are largely captured by the DMX+FD
timing model.

DMX measurements are useful for evaluating ISM fluctuations across a wide
range of length scales, and are often required for any analysis of sub-$\mu$s
timing effects like the Shapiro delay and nHz-frequency gravitational
radiation (Jones et al., 2017; Valentina Sosa Fiscella et al., 2023). The DM
variations observed in PSR J0218+4232 are particularly well-resolved as shown
in Figure 1, such that a lack of DMX evaluation would contribute substantial
“red” noise to our timing measurements and dominate the best-fit
uncertainties. We reserve an analysis of the DMX timeseries for PSR J0218+4232
and other sources for a forthcoming census of low-frequency DM measurements
with CHIME/Pulsar (McKee et al., in preparation).

###  3.3 Goodness-of-Fit Grids of Shapiro Delay Parameters

With the timing model in hand, we used the Bayesian $\chi^{2}$-gridding method
developed by Splaver et al. (2002) to better determine the statistics of the
Shapiro delay parameters for each pulsar. In this method, we iterated over the
Shapiro delay parameters in a $100\times 100$ grid of ($m_{\rm c}$, $\cos i$)
values. At each ($m_{\rm c}$, $\cos i$) grid point, the timing model of each
pulsar is refitted while fixing the values of the Shapiro delay parameters. A
2-D probability density function (PDF), $p(m_{\rm c},\cos i|{\rm data})$, is
then generated by comparing the reduced $\chi^{2}$ values obtained at
different Shapiro delay parameter value pairs.

The two-dimensional PDFs obtained for each pulsar were then used to constrain
the Shapiro-delay parameters and the mass of the pulsar. For the estimates of
$m_{\rm p}$, we used the Keplerian mass function ($m_{\rm f}$), where

$\displaystyle m_{\rm f}$ $\displaystyle=\frac{{(}m_{\rm c}\sin
i{)}^{3}}{m_{\rm tot}^{2}}{=\frac{4\pi^{2}}{T_{\odot}}\frac{x^{3}}{P_{\rm
b}^{2}}},$ (9)

in order to enact the transformation in PDF variables:

$\displaystyle p(m_{\rm c}|{\rm data})$
$\displaystyle=\int_{0}^{\infty}p(m_{\rm c},\cos i|{\rm data})d(\cos i),$ (10)
$\displaystyle p(m_{\rm p}|{\rm data})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}$ $\displaystyle p(m_{\rm
c},\cos i|{\rm data})\delta(\Delta m_{\rm p})d(\cos i)d(m_{\rm c}),$ (11)
$\displaystyle p(\cos i|{\rm data})$ $\displaystyle=\int_{0}^{\infty}p(m_{\rm
c},\cos i|{\rm data})d(m_{\rm c}),$ (12)

where $\delta(x)$ is the Dirac delta function, and $\Delta m_{\rm p}=m_{\rm
p}-(\sqrt{(m_{\rm c}\sin i)^{3}/m_{\rm f}}-m_{\rm c})$. We computed the
cumulative distribution functions (CDFs) for each of the three above
parameters to derive the median, 68.3% and 95.4% credible-interval values for
robust estimates of the relativistic parameters and their statistical
uncertainties.

###  3.4 Constrained Shapiro-delay Grid for PSR J1518+4904

PSR J1518+4904 exhibits apsidal motion that has long been known to be
consistent with the predictions of general relativity (Nice et al., 1996;
Janssen et al., 2008), as given in Equation 4. Since $\dot{\omega}$ is a
function of $m_{\rm tot}$, then its statistical significance can be used as a
weight on grid-based estimates of ($m_{\rm c}$, $\cos i$) and further delimit
regions of preferred values, as has been done for other high-eccentricity
systems. We therefore performed an additional $\chi^{2}$-grid calculation for
CHIME/Pulsar timing data of PSR J1518+4904where, at each grid point, we held
values of $m_{\rm c}$, $m_{\rm tot}$, and $\dot{\omega}$ fixed, using Equation
9 to compute and fix $\sin i$.

We chose to perform grid calculations for PSR J1518+4904 in the ($m_{\rm c}$,
$m_{\rm tot}$) phase space due to the high statistical significance of
$\dot{\omega}$. After extensive analysis, we found that a grid over the
($m_{\rm c}$, $m_{\rm tot}$) phase space produces better constraints on the
probability densities compared with that from the traditional ($m_{\rm c}$,
$\cos i$) space. These calculations ultimately produced a map of probability
density $p(m_{\rm c},m_{\rm tot}|{\rm data})$. Using this map, we derived
estimates of the Shapiro-delay parameters by marginalizing over the
appropriate coordinate axes:

$\displaystyle p(m_{\rm c}|{\rm data})$
$\displaystyle=\int_{0}^{\infty}p(m_{\rm c},m_{\rm tot}|{\rm
data})\delta(\Delta\dot{\omega})d(m_{\rm tot}),$ (13) $\displaystyle p(m_{\rm
p}|{\rm data})$ $\displaystyle=\int_{0}^{\infty}\int_{0}^{\infty}p(m_{\rm
c},m_{\rm tot}|{\rm data})$ $\displaystyle\phantom{0}{}\times\delta(\Delta
m_{\rm p})d(m_{\rm tot})d(m_{\rm c}),$ (14) $\displaystyle p(\cos i|{\rm
data})$ $\displaystyle=\int_{0}^{\infty}\int_{0}^{\infty}p(m_{\rm c},m_{\rm
tot}|{\rm data})$ $\displaystyle\phantom{0}{}\times\delta(\Delta\cos
i)d(m_{\rm tot})d(m_{\rm c}),$ (15)

where $\Delta\cos i=\cos i-\sqrt{1-(m_{\rm f}m_{\rm tot}^{2}/m_{\rm
c}^{3})^{2/3}}$.

Figure 1: Frequency-resolved and frequency-averaged-timing residuals (top and
middle panels, respectively), and dispersion measure variation (bottom panels)
of PSR J0218+4232 and PSR J1518+4904 over the observing span, as measured
using the DMX method. While there is evidence for an annual variation in the
DM of PSR J1518+4904, the ecliptic latitude of PSR J1518+4904,
$\beta=63^{\circ}$, is too large for electron-density variations from the Sun
to appear prominently in our DM timeseries. Figure 2: Frequency-resolved and-
frequency averaged timing residuals (top and middle panels, respectively), and
dispersion measure variation (bottom panel) of PSR J2023+2853 over the
observing span, as measured using the DMX method.

Figure 3: The probability density maps of Shapiro-delay parameters for PSRs
J0218$+$4232, J1518$+$4904, and J2023$+$2853\. The red contours in each map
represent regions of probability that contain 68.3%, 95.4%, and 99.3%
credibility. In the outer panels, red-solid lines denote median values while
the red-dashed lines represent 68.3% credible intervals. For PSR J1518+4904,
the black curves in each map represent the range of allowed masses inferred by
the observed $\dot{\omega}$ when assuming that general relativity describes
the apsidal motion; the precision in our measurement of $\dot{\omega}$ is so
high that the two curves appear as a single line in each map. No credible
intervals are presented in the outer panels for PSR J1518+4904 as the majority
of probability density is not fully encapsulated in the map bounds. Figure 4:
The constrained probability density of the neutron-star masses in the PSR
J1518+4904 system, derived from posterior PDFs of the Shapiro-delay parameters
computed while using the observed apsidal motion as a constraint on the
likelihood function. The red solid line is the median value while the red-
dashed lines represent 68.3% credible intervals of the masses. Table 1: Best-
fit parameters and derived quantities for binary MSPs.  Global Parameters
---
Pulsar name | PSR J0218+4232 | PSR J1518+4904 | PSR J2023+2853
Reference epoch (MJD) | 58917 | 59134.411 | 59674
Observing timespan (MJD) | 58652-59930 | 58552-59929 | 59376-59930
Solar System Planetary Ephemeris | DE440 | DE440 | DE440
Binary model | ELL1 | DDGR | ELL1
Clock Standard | TT(BIPM2021) | TT(BIPM2021) | TT(BIPM2021)
Number of Sub-bands | 32 | 8 | 8
Width of DMX time bin (days) | 5 | 7 | 10
TOA Uncertainty Adjustment Factor (EFAC) | 1.2 | 1.25 | 1.25
Timing Solution & Best-fit Metrics
Right ascension (J2000), $\alpha$ (h:m:s) | 02:18:06.362483(15) | 15:18:16.797784(12) | 20:23:21:06320(7)
Declination (J2000), $\delta$ (d:m:s) | 42:32:17.3417(3) | 49:04:34.08325(12) | 28:53:41.4423(7)
Proper motion in Right ascension, $\mu_{\alpha}\cos\delta$ ($\rm mas\,yr^{-1}$) | 5.39(11) | $-0.72(12)$ | $-2.6(19)$
Proper motion in Declination, $\mu_{\delta}$ ($\rm mas\,yr^{-1}$) | 1.4(2) | $-8.57(11)$ | $-6(4)$
Parallax, $\varpi$ (mas) | — | 4(2) | —
Pulse frequency, $\nu$ (s-1) | 430.461056366647(6) | 24.4289797496997(2) | 88.269785050593(8)
First pulse frequency derivative, $\dot{\nu}$ ($10^{-15}$ s-2) | $-14.34202(16)$ | $-0.016197(10)$ | $-0.212(3)$
Dispersion measure, DM ($\mathrm{pc\,cm^{-3}}$) | 61.233749 | 11.611711 | 22.75
Orbital period, $P_{\rm b}$ (days) | 2.02884608487(5) | 8.63400496116(15) | 0.71823040745(4)
Projected semi-major axis, $x$ (lt-s) | 1.98443196(5) | 20.03942440(11) | 4.0022194(4)
Epoch of periastron passage, $T_{0}$ (MJD) | — | 59125.99829199(5) | —
Orbital eccentricity, $e$ | — | 0.249484383(9) | —
Longitude of periastron, $\omega$ (deg) | — | 342.745426(2) | —
Epoch of ascending node passage, $T_{ASC}$ (MJD) | 58915.445027883(14) | — | 59390.14830429(2)
First Laplace-Lagrange parameter, $\epsilon_{1}$ (10-5) | 0.39(2) | — | 1.11(2)
Second Laplace-Lagrange parameter, $\epsilon_{2}$ (10-5) | 0.50(2) | — | 0.73(2)
First order frequency dependent variation in channelised TOAs | $-1.54(15)\times 10^{-5}$ | — | —
Number of TOAs | 29889 | 7412 | 2623
$\chi^{2}$ | 30192.25 | 7458.33 | 2690.70
Degrees of freedom | 29675 | 7241 | 2560
Reduced $\chi^{2}$ | 1.02 | 1.03 | 1.05
RMS timing residual ($\mu s$) | 6.458 | 4.906 | 12.395
Gridded Estimates of the Masses and Orbital Inclination
Mass of the companion, $m_{\rm c}$ (M⊙) | $0.179^{+0.018}_{-0.016}$ | $1.248^{+0.035}_{-0.029}$ | $0.93^{+0.17}_{-0.14}$
Total mass of the system, $m_{\rm tot}$ (M⊙) | — | $2.7186(7)$ | —
Cosine of the inclination angle, $cos\ i$ | $0.086^{+0.020}_{-0.018}$ | — | $0.123^{+0.036}_{-0.032}$
Derived Quantities
Mass function, $m_{\rm f}$ (M⊙) | 0.0020384176(2) | 0.115908122(2) | 0.13343133(4)
Mass of the pulsar, $m_{\rm p}$ (M⊙) | $1.49^{+0.23}_{-0.20}$ | $1.470^{+0.030}_{-0.034}$ | $1.50^{+0.49}_{-0.38}$
Inclination angle of the binary system, $i$, (degrees) | $85.1^{+1.0}_{-1.2}$ | $49.6^{+1.6}_{-1.8}$ | $83.0^{+1.8}_{-2.0}$
Rate of advance of periastron, $\dot{\omega}$ (deg yr-1) | — | 0.011373(2) | —
Time-dilation and gravitational-redshift parameter, $\gamma$ | — | $0.00463^{+0.00017}_{-0.00014}$ | —
Relativistic orbital decay, $\dot{P_{\rm b}}$ ($10^{-12}$) | — | $-0.001176(5)$ | —

##  4 Results and Discussion

We obtained best-fit timing models for all three binary pulsars using TEMPO
and the methods described in Section 3. The best-fit timing residuals and
available DMX data for all binary pulsars are presented in Figure 1. A summary
of each timing model, derived parameters, and best-fit statistics is presented
in Table 1, while fixing the Shapiro delay values to the median values
obtained from the $\chi^{2}$ grids described in Section 3.3.

###  4.1 New Shapiro-delay Measurements

For the first time, we detected the Shapiro time delay in all three binary
pulsar systems. The near-daily observing cadence has clearly led to desirable
coverage of the orbital phase over a 2–3-yr timescale; the lack of detection
in prior work on PSRs J0218$+$4232 (Navarro et al., 1995; Desvignes et al.,
2016; Perera et al., 2019) and J1518$+$4904 (Janssen et al., 2008) was likely
impacted by low-density orbital phase coverage despite their several-decade
timespan.

We consider the Shapiro delay in PSR J1518+4904 to be detected since total
degeneracy of the ($m_{\rm c}$, $\cos i$) parameters – regions where
probability density is non-zero as $\cos i$ tends to unity – is statistically
disfavored within the range of plausible companion-mass values. However, the
Shapiro delay in the PSR J1518+4904 system is nonetheless weakly detected as
the marginalized constraints on $m_{\rm c}$ and $\cos i$ produce a weak
constraint on $m_{\rm p}$. As a standalone effect, the current Shapiro delay
in PSR J1518+4904 therefore offers no meaningful constraining power to the
mass and geometry of the binary system, other than to restrict likely values
of $i$. However, the use of $\dot{\omega}$ as a statistical weight in the
$\chi^{2}$-grid calculation leads to substantial improvements in the
constraints of physical parameters. We found that $m_{\rm
p}=1.470^{+0.030}_{-0.034}$ M⊙ and $m_{\rm c}=1.248^{+0.035}_{-0.029}$ M⊙
(68.3% credibility), under the assumption that general relativity describes
the observed apsidal motion.

Our estimates of $m_{\rm p}$ for PSRs J0218$+$4232 and J2023$+$2853 – $m_{\rm
p}=1.48^{+0.24}_{-0.20}\textrm{ M}_{\odot}$ and $m_{\rm
p}=1.42^{+0.40}_{-0.32}\textrm{ M}_{\odot}$, respectively – are broadly
consistent with values observed for the recycled-pulsar population, though
additional observations will further strengthen the constraints for
delineating whether these pulsars lie on the high end of the neutron-star mass
spectrum. The estimates of the neutron-star masses found in the PSR J1518+4904
system are consistent with masses found in other DNS systems (e.g., Schwab et
al., 2010).

###  4.2 Discussion on individual sources

####  4.2.1 PSR J0218+4232

The properties of the PSR J0218+4232 system, first discovered by Navarro et
al. (1995), were most recently modelled by Desvignes et al. (2016) and later
Perera et al. (2019) as part of the International Pulsar Timing Array (IPTA)
project. No apparent detection of the Shapiro delay was made in either of
their extended data sets. By contrast, we were able to successfully detect the
Shapiro-delay signature in the PSR J0218+4232 system due to the near-daily
cadence of our CHIME/Pulsar data set. The measurement of $\cos i$ is
consistent with constraints placed on orbital geometry and spin orientation of
PSR J0218+4232 derived from low-frequency polarimetry (Stairs et al., 1999),
and suggests that the misalignment between the axes of orbital angular
momentum and pulsar spin are nearly orthogonal.

The statistical significance of the Shapiro delay remains high even when
incorporating hundreds of degrees of freedom from the DMX model into our
timing analysis. However, when included as a free parameter, the best-fit
measured proper motion in declination was found to be inconsistent with the
measurements by both Desvignes et al. (2016) and Perera et al. (2019), as well
as a Very Long Baseline Interferometry (VLBI) measurement by (Du et al.,
2014), both in magnitude and sign. Moreover, we found that fits of the timing
signature due to parallax yielded negative values with apparent statistical
significance. We re-modeled the properties of PSR J0218+4232 after fixing the
parallax to zero and proper motion to the model-independent values determined
by Du et al. (2014) of $\mu_{\alpha}\cos\delta=5.35$ mas yr-1 and
$\mu_{\delta}=-3.74$ mas yr-1, and found that the best-fit DM variations of
the pulsar changed significantly. We suspect that these discrepancies arise
due to several reasons: the sensitivity of our data set on PSR J0218+4232 to
intra-day DM variations; our observations being scheduled such that the pulsar
is observed on the same local sidereal time every day; and that the binary
orbit of PSR J0218+4232 is slightly more than exactly 2 days, resulting in the
variation of the orbital Doppler shift being absorbed into the DM-variation
modelling. This apparent DM variation is found to be degenerate with the
astrometric parameters and manifests itself as an offset to the measured
proper motion in declination. These issues do not arise with the other pulsars
analyzed in our study.

We nonetheless found that the difference in the proper motion used to model
the binary system resulted in an insignificant difference in the measured
Shapiro Delay parameters. We found $m_{\textrm{c}}$ $=0.182^{+0.019}_{-0.017}$
M⊙ and $\cos\ i=0.0911^{+0.020}_{-0.019}$ after fixing the proper motion of
the system to be the values obtained by Du et al. (2014). These values remain
consistent at the 68.3% credibility level with values derived from a timing
model that sets the proper motion of the system as free parameters.

The CHIME/Pulsar timing model provides a high-cadence snapshot of the PSR
J0218+4232 orbit that may have evolved over the decades since its discovery.
We therefore compared our values of the orbital elements with those most
recently published for the second data release of IPTA (Perera et al., 2019);
the IPTA data set for PSR J0218+4232 spans $\sim$17.5 years and consists of
TOAs collected with a variety of European observatories. As the timing model
from Perera et al. (2019) does not include any Shapiro-delay modelling, we
first downloaded the publicly available data from their work and refit the
release model for PSR J0218+4232 after incorporating the Shapiro-delay
parameters as degrees of freedom. This refitting was done to ensure robust
comparisons of models, and to avoid potential biasing of orbital parameters
due to the presence of an unmodeled Shapiro delay (Freire & Wex, 2010).

Using this new IPTA model, we found that $x$ has remained the same between the
two orbital models to within $1\sigma$, where $\sigma$ is the larger of the
two statistical uncertainties in the fit parameter. We inferred an approximate
rate of change in $x$ to be $\dot{x}\approx\Delta x/\Delta T_{\rm
asc}=-5(4)\times 10^{-15}$ s/s. This limit on $\dot{x}$ is consistent with the
maximum value of $\dot{x}$ due to the changing of orientation due to proper
motion: given our estimate of $i$ and that $(\dot{x})_{\rm max}=x\mu|\cot{i}|$
(Nice et al., 2001), we found that $|\dot{x}|_{\rm max}\sim 10^{-16}$ s/s.

By contrast, we found that there is a significant change in $P_{\rm b}$
between our CHIME/Pulsar timing model and the Shapiro-delay incorporated model
derived from the IPTA data set from Perera et al. (2019). We computed the
change in $P_{\rm b}$ between the two epochs to be $\dot{P_{\rm
b}}\approx\Delta P_{\rm b}/\Delta T_{\rm asc}=0.14(2)\times 10^{-12}$, where
the uncertainty is dominated by those estimated for $P_{\rm b}$.

The fact that $\dot{P_{\rm b}}$ is positive suggests this variation is not
dominated by decay due to the emission of gravitational waves from the system.
Instead, this rate of change most likely arises from variations in the Doppler
shift due the relative motion between PSR J0218+4232 and Earth. This motion
results in contributions of $\dot{P}_{\rm b}$ to the total, observed value,
such that ($\dot{P}_{\rm b})_{\rm obs}$ = $(\dot{P}_{\rm b})_{\rm GR}$ \+
$(\dot{P}_{\rm b})_{\rm DR}$ \+ $(\dot{P}_{\rm b})_{z}$ \+ $(\dot{P}_{\rm
b})_{\mu}$, where $(\dot{P}_{\rm b})_{\rm GR}$ is general-relativistic orbital
decay rate presented in Equation 3.1 and the additional terms are defined as
follows (Shklovskii, 1970; Kuijken & Gilmore, 1989; Damour & Taylor, 1992;
Nice et al., 1996):

$\displaystyle(\dot{P}_{\rm b})_{\rm DR}$ $\displaystyle=-P_{\rm b}\cos
b\bigg{(}\frac{\Theta_{0}^{2}}{cR_{0}}\bigg{)}\bigg{(}\cos
l+\frac{\kappa^{2}}{\sin^{2}l+\kappa^{2}}\bigg{)},$ (16)
$\displaystyle(\dot{P}_{\rm b})_{\rm z}$ $\displaystyle=-1.08\times
10^{-19}\frac{P_{\rm b}}{c}$
$\displaystyle\times\bigg{(}\frac{1.25z}{[z^{2}+0.0324]^{1/2}}+0.58z\bigg{)}\sin
b,$ (17) $\displaystyle(\dot{P}_{\rm b})_{\mu}$
$\displaystyle=\frac{\mu^{2}d}{c}P_{\rm b}.$ (18)

The various terms in Equations 16–18 are defined as follows: $b$ and $l$ are
the Galactic latitude and longitude, respectively; $\Theta_{0}$ and $R_{0}$
are the Galactrocentric circular-speed and distance parameters for the Solar-
System barycenter, respectively; $\kappa=(d/R_{0})\cos b-\cos l$, where $d$ is
the distance to the pulsar; $z=d\sin b$ is the projected vertical distance of
the pulsar-binary system from the Galactic plane; and
$\mu=\sqrt{\mu_{\alpha}^{2}+\mu_{\delta}^{2}}$ is the magnitude of proper
motion. We assumed $\Theta_{0}$ = 236.9(4.2) km s-1 and $R_{0}$ = 8.178(26)
kpc as recently determined by Gravity Collaboration et al. (2019).

By modelling $(\dot{P}_{\rm b})_{\rm obs}$ in terms of Equations 3.1 and
16–18, we found that large values of $d$ are needed to produce our derived
value. We performed a Monte Carlo analysis of $d$ by randomly sampling the
{$m_{\rm p}$, $m_{\rm c}$, $\Theta_{0}$, $R_{0}$, $\mu$, $(\dot{P}_{\rm
b})_{\rm obs}$} parameters, based on their best available uncertainties that
we assumed to be Gaussian222While our Shapiro-delay analysis shows that the
posterior PDFs for $m_{\rm p}$ and $m_{\rm c}$ are non-Gaussian, the
$(\dot{P}_{\rm b})_{\rm GR}$ term is subdominant to all other sources of
variation in $P_{\rm b}$ that therefore has no bearing on our Monte-Carlo
estimate of $d$., and computing $d$ using Equations 3.1 and 16–18 for each
combination of parameter values. Our estimate of $(\dot{P}_{\rm b})_{\rm obs}$
from this Monte Carlo analysis corresponds to $d=(6.7\pm 1.0)$ kpc, where the
statistics reflect the mean and standard deviation of the resultant
distribution in $d$. This distance is statistically consistent with the VLBI
distance determined by Du et al. (2014), $d=6.3$ kpc, where the combined
effects give a predicted estimate of $\dot{P_{\rm b}}\approx 0.12\times
10^{-12}$ for this large distance. Conversely, a distance to the pulsar of
$d=3.15$ kpc, as suggested by Verbiest & Lorimer (2014), would yield an
expected $\dot{P_{\rm b}}\approx 0.06\times 10^{-12}$. This latter value of
$d$ is inconsistent with our estimate of $(\dot{P}_{\rm b})_{\rm obs}$ at the
$4\sigma$ level. Otherwise, the spectroscopic analysis of the white dwarf
companion (Bassa et al., 2003) places the binary pulsar system at a distance
of $d=4$ kpc for a white dwarf model that is consistent with the measured
mass. This gives an expected $\dot{P_{\rm b}}\approx 0.08\times 10^{-12}$,
which differs by $3\sigma$ from our measured $\dot{P_{\rm b}}$ value.

While the difference between our measured $\dot{P_{\rm b}}$ value and expected
$\dot{P_{\rm b}}$ from several distance estimates is tantalising, the specific
form of Equation 4.2.1 assumes a model of vertical acceleration due to the
surface-mass distribution of the Galactic disk estimated by Kuijken & Gilmore
(1989). We separately utilized an alternative model of the Galactic surface-
mass density developed by Holmberg & Flynn (2000, 2004), which predicts that

$\displaystyle(\dot{P}_{\rm b})_{\rm z}$ $\displaystyle=-\bigg{[}2.27|z|_{\rm
kpc}+3.68(1-\exp[-4.61|z|_{\rm kpc}])\bigg{]}$
$\displaystyle\phantom{0}{}\times 10^{-11}\bigg{(}\frac{P_{\rm b}|\sin
b|}{c}\bigg{)}$ (19)

and has been used in recent pulsar-timing analyses (e.g., Lazaridis et al.,
2009; Zhu et al., 2015). Both Equations 4.2.1 and 4.2.1 predict that
$(\dot{P}_{b})_{z}\approx-10^{-14}$ s/s for a distance range $3.15<d<6.3$ kpc,
which are subdominant to the “Shklovskii acceleration” defined in Equation 18.
The preference for a large distance to PSR J0218+4232 therefore remains,
regardless of the choice in a model of $(\dot{P}_{b})_{z}$. However, these
models likely lose accuracy for large $z$; further development of viable
Galactic-acceleration models (e.g., Chakrabarti et al., 2021) will therefore
improve the robustness of our timing-based constraint on $d$ for PSR
J0218+4232.

####  4.2.2 PSR J1518+4904

The most recent analysis of the PSR J1518+4904 binary system was conducted by
Janssen et al. (2008), using a data set that spanned over 12 years. The
properties obtained by Janssen et al. (2008) are broadly similar to the model
we produced in our analysis of CHIME/Pulsar timing data. However, they did not
meaningfully detect the Shapiro delay from their observations. Janssen et al.
(2008) instead derived an upper limit on the inclination of the orbit at
$i<47\deg$, based on their non-detection of the Shapiro delay. Using this
constraint, they estimated the masses of the pulsar and companion to be
$m_{\textrm{p}}$ $<1.17$ M⊙ and $m_{\textrm{c}}$ $>1.55$ M⊙.

In our work, we were able to significantly detect the signature of Shapiro
delay from PSR J1518+4904, with the grid-based inclination of the orbit
estimated to be $i=49.6^{+1.6}_{-1.8}$ deg. The derived inclination angle is
marginally higher than the upper limit placed by Janssen et al. (2008). The
low inclination angle suggests that the Shapiro-delay signature from the PSR
J1518+4904 system is weak, and is only detectable due to a combination of the
significant apsidal motion and high-cadence sampling of our data over much
shorter time span of less than 4 years.

The measured masses of the PSR J1518+4904 system also show that the pulsar is
the more massive object with $m_{\textrm{p}}=1.470^{+0.030}_{-0.034}$ M⊙,
whilst the companion neutron star has a smaller mass of
$m_{\textrm{p}}=1.248^{+0.035}_{-0.029}$ M⊙. The rotational properties of the
pulsar suggest that it is a partially recycled pulsar and would be the first
born neutron star in the system. The first born neutron star being the heavier
object in a double neutron star system is consistent with most others that are
observed (See Tauris et al., 2017; Özel & Freire, 2016, and references
therein.) Tauris et al. (2017) studied the formation of such systems and
suggest that the mass discrepancy between the neutron star is unlikely due to
accretion of mass by the first born neutron star, but rather due to the
progenitor of the second neutron star being stripped of its mass more
significantly due to the first neutron star prior to the second supernova.

We compared our best-fit Keplerian elements with those published by Janssen et
al. (2008), in order to detect or constrain variations between the two data
sets. No significant changes in $x$ or $P_{\rm b}$ beyond 2-$\sigma$ were
found after accounting for clock-correction differences between the two timing
solutions. The limits we derived from this analysis – $|\dot{x}|<3\times
10^{-14}$ and $|\dot{P}_{\rm b}|<6\times 10^{-12}$ – are consistent with those
determined from the longer data set studied by Janssen et al. (2008).

The accuracy and precision in our measurement of $\dot{\omega}$ are nearly
identical to the values published by Janssen et al. (2008). It is worth noting
that this result arises from our analysis despite the CHIME/Pulsar data set
being a factor of $\sim$3.5 shorter in timespan than that analyzed by Janssen
et al. (2008). This circumstance demonstrates that the CHIME/Pulsar instrument
can produce high-accuracy snapshots of pulsar orbits on much faster timescales
than possible elsewhere, owing to the nearly-daily cadence of our
observations, and thus allow for the quickened resolution of orbital
variations over time. As originally noted by Janssen et al. (2008), the
statistical uncertainty in the observed apsidal motion is comparable in
magnitude to variations that arise evolving orientation of the binary system
due to proper motion (e.g., Kopeikin, 1996). We therefore consider our
measured $\dot{\omega}$ to arise purely from general-relativistic orbital
motion. A combination of the CHIME/Pulsar data set with that published by
Janssen et al. (2008), which will lead to substantially improved measurement
of $\dot{\omega}$ and likely other relativistic variations, will be the
subject of future work.

####  4.2.3 PSR J2023+2853

PSR J2023+2853 is a pulsar recently discovered by the FAST Galactic Plane
Pulsar Snapshot survey (Han et al., 2021). Han et al. (2021) noted that PSR
J2023+2853 is detected with significant acceleration, suggesting that the
pulsar has a binary companion. Our independent observations revealed that the
pulsar is in a 0.7-day long orbit, with a detectable Shapiro-delay signature
that revealed a companion with m${}_{\textrm{c}}$ $=0.90^{+0.14}_{-0.12}$ M⊙.
Considering the low eccentricity of the orbit and the mass of the companion,
it is likely to be a high mass carbon-oxygen white dwarf.

The relatively long rotation period of PSR J2023+2853 ($P>10$ ms) compared to
typical millisecond pulsars, combined with the relatively heavy mass of the
companion ($m_{\rm c}>0.5$ M⊙), suggests that it is part of the class of
binary system known as the intermediate-mass binary pulsars (Camilo et al.,
2001, IMBPs). IMBPs are thought to be formed from accretion of a main sequence
companion of 4-8 M⊙, in which the pulsar itself is partially recycled. PSR
J2023+2853 has also a larger eccentricity compared to typical pulsar binary
systems of similar orbital period, as well as a low scale height of
$0.13<z<0.17$ kpc, for an estimated distance to the pulsar of 1.6-2.0 kpc as
predicted by the Galactic electron density models (Cordes & Lazio, 2002; Yao
et al., 2017). These properties are consistent with the system being an IMBP.

An optical source with $r=23.55,g-r=1.43$, SDSS J202320.88+285345.3, is found
in the Sloan Digital Sky Survey (SDSS, York et al., 2000; Abdurro’uf et al.,
2022), around 5” from the timing position of PSR J2023+2853. However, the
photometric properties of the source, specifically the SDSS $u-g,g-r$ and
$r-i$ colors, are inconsistent with it being a white dwarf, using the white
dwarf cooling model described by Bergeron et al. (1995); Holberg & Bergeron
(2006); Kowalski & Saumon (2006); Tremblay et al. (2011); Bergeron et al.
(2011); Blouin et al. (2018); Bédard et al.
(2020)333https://www.astro.umontreal.ca/~bergeron/CoolingModels/, even
considering the Galactic extinction along the line of sight, which gives
$r=24.17,g-r=1.23$ (Schlegel et al., 1998; Schlafly & Finkbeiner, 2011). This
suggests that the optical source is not associated with PSR J2023+2853. No
sources are detected with 5” of PSR J2023+2853 in the Two Micron All Sky
Survey (Skrutskie et al., 2006, 2MASS) and the Gaia mission (Gaia
Collaboration et al., 2016, 2022) The non detection of the WD companion is
expected, as it should have a r-magnitude of 25 in SDSS, smaller than the
detection threshold of the surveys above, based on the r-magnitude of the
detection made on a similar WD in another pulsar-binary system, PSR
J1658$+$3630 (Tan et al., 2020).

##  5 Conclusions

By leveraging the near-daily observations conducted by CHIME/Pulsar, we were
able to, for the first time, model the Shapiro time delay of the binary
pulsars PSRs J0218+4232, J1518+4904 and J2023+2853. In the case of J1518+4904,
we assumed that the observed aspidal motion was described by general
relativity so that its significance could be used to constrain the component
masses of the system. These measurements allowed us to constrain to the masses
of the pulsars in the binary systems to be $1.48^{+0.24}_{-0.20}$ M⊙,
$1.470^{+0.030}_{-0.034}$ M⊙ and $1.50^{+0.49}_{-0.38}$ M⊙ respectively. The
measured pulsar masses of all three system are well within the known mass
distribution of neutron stars (Özel & Freire, 2016).

We were also able to measure the mass of the companion neutron star to PSR
J1518+4904 to be $m_{\rm c}=1.248^{+0.035}_{-0.029}$ M⊙, and for the first
time, model the binary properties of PSR J2023+2853, revealing a relatively
high mass companion of $m_{\rm c}=0.93^{+0.17}_{-0.14}$ M⊙. The measured mass,
together with a near-circular orbit, suggest that the companion is likely a
high-mass carbon-oxygen white dwarf, and that the system is an IMBP. The low
mass of the companion to PSR J0218+4232, a helium white dwarf of $m_{\rm
c}=0.179^{+0.019}_{-0.016}\textrm{ M}_{\odot}$, is consistent with predictions
from models of MSP formation through long-term mass transfer within low-mass
X-ray binary systems (e.g., Tauris & Savonije, 1999).

We also obtained an estimate in the change in orbital period of PSR J0218+4232
between the epoch of the observations from Perera et al. (2019) and our
observations. We found that the primary contribution of the orbital period
change is due to variations in the Doppler shift from Shklovskii acceleration,
and that the change corresponds to a distance to the pulsar of $d$ = ($6.7\pm
1.0$) kpc, consistent with the measurement of $d$ = 6.3 kpc from the VLBI
observation conducted by Du et al. (2014). No significant changes were
observed in the CHIME/Pulsar estimates of the orbital elements of PSR
J1518+4904.

Our success in obtaining measurements of the masses of these pulsars through
Shapiro-delay modelling, using the near-daily observations from CHIME/Pulsar,
has prompted us to observe other binary pulsar systems on a near-daily basis.
These high-cadence observations will help constrain the masses of the objects
in these systems in order to provide an even larger sample of pulsar masses
and relativistic measurments. We also continue to observe the three pulsar
binary systems in this study, in order to further constrain their properties.

Future work on these sources will also investigate the DM variations and noise
properties of the high-cadence CHIME/Pulsar data sets, which could slightly
modify the precision of our mass and geometric measurements. Moreover,
subsequent studies will compute and analyze “wideband” TOAs in order to
maximize arrival-time precision and reduce TOA data volumes for efficient
modeling of timing phenonema.

We acknowledge that CHIME is located on the traditional, ancestral, and
unceded territory of the Syilx/Okanagan people. We thank Marten van Kerkwijk
for valuable discussion. We are grateful to the staff of the Dominion Radio
Astrophysical Observatory, which is operated by the National Research Council
of Canada. CHIME is funded by a grant from the Canada Foundation for
Innovation (CFI) 2012 Leading Edge Fund (Project 31170) and by contributions
from the provinces of British Columbia, Québec and Ontario. The CHIME/FRB
Project, which enabled development in common with the CHIME/Pulsar instrument,
is funded by a grant from the CFI 2015 Innovation Fund (Project 33213) and by
contributions from the provinces of British Columbia and Québec, and by the
Dunlap Institute for Astronomy and Astrophysics at the University of Toronto.
Additional support was provided by the Canadian Institute for Advanced
Research (CIFAR), McGill University and the McGill Space Institute thanks to
the Trottier Family Foundation, and the University of British Columbia. The
CHIME/Pulsar instrument hardware was funded by NSERC RTI-1 grant EQPEQ
458893-2014. This research was enabled in part by support provided by the BC
Digital Research Infrastructure Group and the Digital Research Alliance of
Canada (alliance-can.ca).

## References

  * Abdurro’uf et al. (2022) Abdurro’uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35, doi: 10.3847/1538-4365/ac4414
  * Agazie et al. (2023) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023, ApJ, 951, L10, doi: 10.3847/2041-8213/acda88
  * Andersen et al. (2023) Andersen, B. C., Fonseca, E., McKee, J. W., et al. 2023, ApJ, 943, 57, doi: 10.3847/1538-4357/aca485
  * Bassa et al. (2003) Bassa, C. G., van Kerkwijk, M. H., & Kulkarni, S. R. 2003, A&A, 403, 1067, doi: 10.1051/0004-6361:20030384
  * Bédard et al. (2020) Bédard, A., Bergeron, P., Brassard, P., & Fontaine, G. 2020, ApJ, 901, 93, doi: 10.3847/1538-4357/abafbe
  * Bergeron et al. (1995) Bergeron, P., Wesemael, F., & Beauchamp, A. 1995, Publications of the Astronomical Society of the Pacific, 107, 1047, doi: 10.1086/133661
  * Bergeron et al. (2011) Bergeron, P., Wesemael, F., Dufour, P., et al. 2011, ApJ, 737, 28, doi: 10.1088/0004-637X/737/1/28
  * Blouin et al. (2018) Blouin, S., Dufour, P., & Allard, N. F. 2018, ApJ, 863, 184, doi: 10.3847/1538-4357/aad4a9
  * Camilo et al. (2001) Camilo, F., Lyne, A. G., Manchester, R. N., et al. 2001, ApJ, 548, L187, doi: 10.1086/319120
  * Chakrabarti et al. (2021) Chakrabarti, S., Chang, P., Lam, M. T., Vigeland, S. J., & Quillen, A. C. 2021, ApJ, 907, L26, doi: 10.3847/2041-8213/abd635
  * CHIME Collaboration et al. (2022) CHIME Collaboration, Amiri, M., Bandura, K., et al. 2022, ApJS, 261, 29, doi: 10.3847/1538-4365/ac6fd9
  * CHIME/Pulsar Collaboration et al. (2021) CHIME/Pulsar Collaboration, Amiri, M., Bandura, K. M., et al. 2021, ApJS, 255, 5, doi: 10.3847/1538-4365/abfdcb
  * Cordes & Lazio (2002) Cordes, J. M., & Lazio, T. J. W. 2002, arXiv e-prints, astro, doi: 10.48550/arXiv.astro-ph/0207156
  * Damour & Deruelle (1985) Damour, T., & Deruelle, N. 1985, Ann. Inst. Henri Poincaré Phys. Théor, 43, 107
  * Damour & Deruelle (1986) —. 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
  * Desvignes et al. (2016) Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341, doi: 10.1093/mnras/stw483
  * Du et al. (2014) Du, Y., Yang, J., Campbell, R. M., et al. 2014, ApJ, 782, L38, doi: 10.1088/2041-8205/782/2/L38
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, ApJ, 915, L12, doi: 10.3847/2041-8213/ac03b8
  * Freire & Wex (2010) Freire, P. C. C., & Wex, N. 2010, MNRAS, 409, 199, doi: 10.1111/j.1365-2966.2010.17319.x
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211, doi: 10.48550/arXiv.2208.00211
  * Gravity Collaboration et al. (2019) Gravity Collaboration, Abuter, R., Amorim, A., et al. 2019, A&A”, 625, L10
  * Han et al. (2021) Han, J. L., Wang, C., Wang, P. F., et al. 2021, Research in Astronomy and Astrophysics, 21, 107, doi: 10.1088/1674-4527/21/5/107
  * Hankins & Rickett (1975) Hankins, T. H., & Rickett, B. J. 1975, Methods in Computational Physics, 14, 55, doi: 10.1016/B978-0-12-460814-6.50007-3
  * 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
  * Holberg & Bergeron (2006) Holberg, J. B., & Bergeron, P. 2006, AJ, 132, 1221, doi: 10.1086/505938
  * Holmberg & Flynn (2000) Holmberg, J., & Flynn, C. 2000, MNRAS, 313, 209, doi: 10.1046/j.1365-8711.2000.02905.x
  * Holmberg & Flynn (2004) —. 2004, MNRAS, 352, 440, doi: 10.1111/j.1365-2966.2004.07931.x
  * Hotan et al. (2004) Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302, doi: 10.1071/AS04022
  * Jankowski et al. (2019) Jankowski, F., Bailes, M., van Straten, W., et al. 2019, MNRAS, 484, 3691, doi: 10.1093/mnras/sty3390
  * Janssen et al. (2008) Janssen, G. H., Stappers, B. W., Kramer, M., et al. 2008, A&A, 490, 753, doi: 10.1051/0004-6361:200810076
  * Jones et al. (2017) Jones, M. L., McLaughlin, M. A., Lam, M. T., et al. 2017, ApJ, 841, 125, doi: 10.3847/1538-4357/aa73df
  * Kopeikin (1996) Kopeikin, S. M. 1996, ApJ, 467, L93, doi: 10.1086/310201
  * Kowalski & Saumon (2006) Kowalski, P. M., & Saumon, D. 2006, ApJ, 651, L137, doi: 10.1086/509723
  * Kuijken & Gilmore (1989) Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 605, doi: 10.1093/mnras/239.2.605
  * Lam (2018) Lam, M. T. 2018, ApJ, 868, 33, doi: 10.3847/1538-4357/aae533
  * Lange et al. (2001) Lange, C., Camilo, F., Wex, N., et al. 2001, MNRAS, 326, 274, doi: 10.1046/j.1365-8711.2001.04606.x
  * Lazaridis et al. (2009) Lazaridis, K., Wex, N., Jessner, A., et al. 2009, MNRAS, 400, 805, doi: 10.1111/j.1365-2966.2009.15481.x
  * Lower et al. (2020) Lower, M. E., Bailes, M., Shannon, R. M., et al. 2020, MNRAS, 494, 228, doi: 10.1093/mnras/staa615
  * McKee et al. (in preparation) McKee et al., J. W. in preparation
  * Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJ, 918, L28, doi: 10.3847/2041-8213/ac089b
  * Morello et al. (2019) Morello, V., Barr, E. D., Cooper, S., et al. 2019, MNRAS, 483, 3673, doi: 10.1093/mnras/sty3328
  * NANOGrav Collaboration et al. (2015) NANOGrav Collaboration, Arzoumanian, Z., Brazier, A., et al. 2015, ApJ, 813, 65, doi: 10.1088/0004-637X/813/1/65
  * Navarro et al. (1995) Navarro, J., de Bruyn, A. G., Frail, D. A., Kulkarni, S. R., & Lyne, A. G. 1995, ApJ, 455, L55, doi: 10.1086/309816
  * Nice et al. (1996) Nice, D. J., Sayer, R. W., & Taylor, J. H. 1996, ApJ, 466, L87, doi: 10.1086/310178
  * Nice et al. (2001) Nice, D. J., Splaver, E. M., & Stairs, I. H. 2001, ApJ, 549, 516, doi: 10.1086/319079
  * Özel & Freire (2016) Özel, F., & Freire, P. 2016, ARA&A, 54, 401, doi: 10.1146/annurev-astro-081915-023322
  * Parent et al. (2022) Parent, E., Sewalls, H., Freire, P. C. C., et al. 2022, ApJ, 924, 135, doi: 10.3847/1538-4357/ac375d
  * Pennucci et al. (2014) Pennucci, T. T., Demorest, P. B., & Ransom, S. M. 2014, ApJ, 790, 93, doi: 10.1088/0004-637X/790/2/93
  * Perera et al. (2019) Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019, MNRAS, 490, 4666, doi: 10.1093/mnras/stz2857
  * 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
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Schwab et al. (2010) Schwab, J., Podsiadlowski, P., & Rappaport, S. 2010, ApJ, 719, 722, doi: 10.1088/0004-637X/719/1/722
  * Shklovskii (1970) Shklovskii, I. S. 1970, Soviet Ast., 13, 562
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Splaver et al. (2002) Splaver, E. M., Nice, D. J., Arzoumanian, Z., et al. 2002, ApJ, 581, 509, doi: 10.1086/344202
  * Stairs et al. (1999) Stairs, I. H., Thorsett, S. E., & Camilo, F. 1999, ApJS, 123, 627, doi: 10.1086/313245
  * Stovall et al. (2014) Stovall, K., Lynch, R. S., Ransom, S. M., et al. 2014, ApJ, 791, 67, doi: 10.1088/0004-637X/791/1/67
  * Tan et al. (2020) Tan, C. M., Bassa, C. G., Cooper, S., et al. 2020, MNRAS, 492, 5878, doi: 10.1093/mnras/staa113
  * Tauris & Savonije (1999) Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928, doi: 10.48550/arXiv.astro-ph/9909147
  * Tauris et al. (2017) Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170, doi: 10.3847/1538-4357/aa7e89
  * Taylor (1992) Taylor, J. H. 1992, Philosophical Transactions of the Royal Society of London Series A, 341, 117, doi: 10.1098/rsta.1992.0088
  * Tremblay et al. (2011) Tremblay, P. E., Bergeron, P., & Gianninas, A. 2011, ApJ, 730, 128, doi: 10.1088/0004-637X/730/2/128
  * Valentina Sosa Fiscella et al. (2023) Valentina Sosa Fiscella, S., Lam, M. T., Arzoumanian, Z., et al. 2023, arXiv e-prints, arXiv:2307.13248, doi: 10.48550/arXiv.2307.13248
  * van Straten & Bailes (2011) van Straten, W., & Bailes, M. 2011, PASA, 28, 1, doi: 10.1071/AS10021
  * van Straten et al. (2012) van Straten, W., Demorest, P., & Oslowski, S. 2012, Astronomical Research and Technology, 9, 237
  * Verbiest & Lorimer (2014) Verbiest, J. P. W., & Lorimer, D. R. 2014, MNRAS, 444, 1859, doi: 10.1093/mnras/stu1560
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29, doi: 10.3847/1538-4357/835/1/29
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513
  * Zhu et al. (2015) Zhu, W. W., Stairs, I. H., Demorest, P. B., et al. 2015, ApJ, 809, 41, doi: 10.1088/0004-637X/809/1/41


# CP tests for high-$p_{T}$ multi-leptons

Shaouly Bar-Shalom  Physics Department, Technion–Institute of Technology,
Haifa, Israel

###### Abstract

We propose a generic and model-independent search strategy for probing tree-
level CP-violation in inclusive multi-lepton signals and we extend the
standard expression for tree-level CP-violation in scattering processes at the
LHC to include the cases where the initial state in not self-conjugate. We
then use TeV-scale 4-fermion operators of the form $tu\ell\ell$ and
$tc\ell\ell$ with complex Wilson coefficients as an illustrative example and
show that it may generate ${\cal O}(10\%)$ CP asymmetries that should be
accessible at the LHC with an integrated luminosity of ${\cal O}(1000)$ fb-1.

##  1 Introduction

CP-violation (CPV) may be the key to a deeper understanding of particle
physics and the evolution of the universe. It has far-reaching implications
for cosmology since it is needed to explain the observed baryon asymmetry of
the universe, while CPV within the SM is insufficient for that. Furthermore,
CP is not a symmetry of nature, so that on general grounds, one expects any
generic new physics to entail Beyond the SM (BSM) CP-odd phase(s), e.g., SUSY
models, Mutli-Higgs models, Leptoquarks, Vector-like fermions etc …

In this talk, which is based on [1], we re-examine the formulation of tree-
level CP-violating effects in scattering processes at the LHC and we introduce
a new term to the ”master” CPV expression, which takes into account ”fake” CP-
violating effects that arise when the initial state in not self-conjugate.
Based on that, we motivate an inclusive search for CPV in tri-lepton and four-
leptons events at the LHC:

$\displaystyle pp\to\ell^{\prime-}\ell^{+}\ell^{-}+X_{3}\leavevmode\nobreak\
,$ $\displaystyle
pp\to\ell^{\prime+}\ell^{-}\ell^{+}+\bar{X}_{3}\leavevmode\nobreak\ ,$
$\displaystyle
pp\to\ell^{\prime+}\ell^{\prime-}\ell^{+}\ell^{-}+X_{4}\leavevmode\nobreak\ ,$
(1)

where $\ell,\ell^{\prime}=e,\mu,\tau$ (preferably $\ell\neq\ell^{\prime}$, see
below) and $X_{3}$, $\bar{X}_{3}$ and $X_{4}$ contain in general jets and
missing energy. These include the $e^{\pm}\mu^{+}\mu^{-}$ and
$\mu^{\pm}e^{+}e^{-}$ final states for $\ell,\ell^{\prime}=e,\mu$ and
similarly for the pairs $\ell,\ell^{\prime}=e,\tau$ and
$\ell,\ell^{\prime}=\mu,\tau$ and also the 3-flavor $e\mu\tau$ final state.

Indeed, the multi-leptons signals of Eq. 1 are excellent test ground for new
physics (NP), being the ”end-signals” of many interesting processes such as
$pp\to t\bar{t}V,\leavevmode\nobreak\ t\bar{t}H,tV,tt\bar{t}\bar{t},VV,VVV...$
($V=Z,W$). In particular, they are rich and clean signals in the hadronic
environment of the LHC and are sensitive to many types of underlying NP, e.g.,
lepton-flavor violation, lepton universality violation, lepton-number
violation (same-sign leptons), CP violation etc … Moreover, High-$p_{T}$ (TeV
energies) charged-leptons are still relatively unexplored and are relatively
easily identifiable objects. All that and more were explored recently in
detail [2, 3, 4, 5, 6].

##  2 CPV in multi-leptons events

We construct below generic tests of CP in multi-lepton processes, focusing on
the tri-lepton events of Eq. 1, but our derivation below applies also to the
4-leptons signals of Eq. 1.

We would like to emphasize two important points regarding CPV test in multi-
leptons events at the LHC:

  * •

Sizeable, say ${\cal O}(1\%)$ manifestation of CPV in these multi-leptons
events will be an unambiguous indication of NP, since the SM background to CP
in these processes is expected to be negligible (un-observably small); the CPV
effects from the SM’s CP-odd phase of the CKM matrix can arise only from EW
processes at a higher loop order.

  * •

CP studies at the LHC are more complicated, since the initial state is not
self-conjugate. We thus developed a rigorous formula for testing CPV in an
LHC-like environment, which includes a modification to the CP “master”
formula.

###  2.1 Constructing observables for tree-level CPV

Consider the underlying hard process (parton level) for tri-lepton production
at the LHC (and its CC channel):

$ab\to\ell^{\prime-}\ell^{+}\ell^{-}\leavevmode\nobreak\ \leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\
\bar{a}\bar{b}\to\ell^{\prime-}\ell^{+}\ell^{-}$ (2)

For CPV, at least two interfering amplitudes are required ($M_{1}$ and
$M_{2}$):

$\displaystyle{\cal{M}}_{ab\to\ell^{\prime-}\ell^{+}\ell^{-}}$
$\displaystyle=$ $\displaystyle
M_{1}e^{i\left(\phi_{1}+\delta_{1}\right)}+M_{2}e^{i\left(\phi_{2}+\delta_{2}\right)}\leavevmode\nobreak\
,$ (3)

so that for the charged-conjugate channel we have

$\displaystyle\bar{\cal{M}}_{\bar{a}\bar{b}\to\ell^{\prime+}\ell^{-}\ell^{+}}$
$\displaystyle=$ $\displaystyle
M_{1}e^{i\left(-\phi_{1}+\delta_{1}\right)}+M_{2}e^{i\left(-\phi_{2}+\delta_{2}\right)}\leavevmode\nobreak\
,$ (4)

where $\phi_{1,2}$ and $\delta_{1,2}$ are CP-odd and CP-even phases,
respectively. We will use below the relative phases
$\Delta\phi=\phi_{1}-\phi_{2}$, $\Delta\delta=\delta_{1}-\delta_{2}$.

It is useful to classify CP observables according to their transformation
properties under ”naive time reversal” ($T_{N}$), where $T_{N}$: time $\to$
-time.[7] This is presented in Table 1. In particular, a $T_{N}$-odd CP-odd
observable requires only a CP-odd phase, which could arise already at the
tree-level, while a $T_{N}$-even CP-odd observable needs also a non-vanishing
CP-even phase from final state interactions (FSI) which is typically of a
higher order.

Table 1: Classification of CP-odd observables under naive time reversal $T_{N}$, for the above tri-lepton processes, where $\Delta\delta\equiv\delta_{1}-\delta_{2}$, $\Delta\phi\equiv\phi_{1}-\phi_{2}$ and $\delta_{i},\phi_{i}$ are defined within the amplitudes ${\cal{M}}_{ab\to\ell^{\prime-}\ell^{+}\ell^{-}}$ and $\bar{\cal{M}}_{\bar{a}\bar{b}\to\ell^{\prime+}\ell^{-}\ell^{+}}$ in Eqs. 3 and 4.  Type |  $T_{N}$-odd (CP-odd) |  $T_{N}$-odd (CP-even)
---|---|---
CP-asymmetry | $A_{CP}\propto\cos\Delta\delta\sin\Delta\phi$ | $A_{CP}\propto\sin\Delta\delta\sin\Delta\phi$
Required phases | only CP-odd | Both CP-odd & CP-even
Sensitivity | tree-level CPV | CP-even phase from FSI
|  | (higher order effect)
Expected size | ${\cal O}(10\%)$ | ${\cal O}(0.1\%)$

Therefore, in order to probe CPV in the tri-lepton channels of Eq. 2, we have
constructed in [1] CP-odd observables based on the following $T_{N}$-odd
triple products (TP): 111Useful TP correlations for CP studies in scattering
processes at the LHC, which involve leptons with jets momenta, e.g., in
$t\bar{t}$ and vector-bosons production, have been also suggested in [8, 9,
10, 11, 12, 13].

$\displaystyle{\cal O}_{\tt
CP}=\vec{p}_{\ell^{\prime-}}\cdot\left(\vec{p}_{\ell^{+}}\times\vec{p}_{\ell^{-}}\right)\leavevmode\nobreak\
\leavevmode\nobreak\ ,\leavevmode\nobreak\ \leavevmode\nobreak\
\overline{{\cal O}_{\tt
CP}}=\vec{p}_{\ell^{\prime+}}\cdot\left(\vec{p}_{\ell^{-}}\times\vec{p}_{\ell^{+}}\right)$
(5)

which are odd under $P$ and under $T_{N}$. Under $C$ and $CP$ they transform
as: $C({\cal O}_{\tt CP})=+\overline{{\cal O}_{\tt CP}}$, $C(\overline{{\cal
O}_{\tt CP}})=+{\cal O}_{\tt CP}$, $CP({\cal O}_{\tt CP})=-\overline{{\cal
O}_{\tt CP}}$ and $CP(\overline{{\cal O}_{\tt CP}})=-{\cal O}_{\tt CP}$.

Using the ${\cal O}_{\tt CP}$’s in Eq. 5, we define the following $T_{N}$-odd
(and also $P$-violating) asymmetries:

$\displaystyle A_{T}$ $\displaystyle\equiv$ $\displaystyle\frac{N\left({\cal
O}_{\tt CP}>0\right)-N\left({\cal O}_{\tt CP}<0\right)}{N\left({\cal O}_{\tt
CP}>0\right)+N\left({\cal O}_{\tt CP}<0\right)}\leavevmode\nobreak\ ,$ (6)
$\displaystyle\bar{A}_{T}$ $\displaystyle\equiv$
$\displaystyle\frac{N\left(-\overline{{\cal O}_{\tt
CP}}>0\right)-N\left(-\overline{{\cal O}_{\tt
CP}}<0\right)}{N\left(-\overline{{\cal O}_{\tt
CP}}>0\right)+N\left(-\overline{{\cal O}_{\tt
CP}}<0\right)}\leavevmode\nobreak\ ,$ (7)

where $N\left({\cal O}_{\tt CP}>0\right)$ is the number of events for which
${\rm sign}({\cal O}_{\tt CP})>0$ is measured etc.

The asymmetries $A_{T}$ and $\bar{A}_{T}$ are sensitive to the CP-odd phase
BUT are not proper CP-asymmetries, since they are not eigenstate of CP, i.e.,
$CP(A_{T})=\bar{A}_{T}$. In particular, we find that
$A_{T}\propto\sin(\Delta\delta+\Delta\phi)$ and
$\bar{A}_{T}\propto\sin(\Delta\delta-\Delta\phi)$, so that in general
$A_{T}\neq 0$ and/or $\bar{A}_{T}\neq 0$ could be generated without CPV, i.e.,
if $\Delta\phi=0$ & $\Delta\delta\neq 0$. Therefore, in order to isolate the
pure CPV effects, we need to combine the information from both $A_{T}$ and
$\bar{A}_{T}$ as follows:

$\displaystyle A_{CP}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(A_{T}-\bar{A}_{T}\right)\leavevmode\nobreak\ .$
(8)

###  2.2 The multi-leptons cross-section & CP-asymmetries

The differential (hard) cross-sections for the tri-lepton production processes
under investigation, $ab\to\ell^{\prime-}\ell^{+}\ell^{-}+X$ and , can be
schematically written as:

$\displaystyle
d\hat{\sigma}(\ell^{\prime+}\ell^{-}\ell^{+})=W+U\cdot\cos(\Delta\delta+\Delta\phi)+V\cdot{\cal
O}_{\tt CP}\cdot\sin(\Delta\delta+\Delta\phi)\leavevmode\nobreak\ ,$ (9)

and similarly for the CC channel
$\bar{a}\bar{b}\to\ell^{\prime+}\ell^{-}\ell^{+}+\bar{X}$ with
$\Delta\phi\to-\Delta\phi$ and ${\cal O}_{\tt CP}\to\overline{{\cal O}_{\tt
CP}}$ in Eq. 9. Note that $W\propto|M_{1}|^{2},|M_{2}|^{2}$, $U\propto{\rm
Re}\left(M_{1}M_{2}^{*}\right)$ and the 3rd term in Eq. 9 (which is
$\propto{\cal O}_{\tt CP}$) arises from ${\rm Im}\left(M_{1}M_{2}^{*}\right)$
and it is where the tree-level CPV resides.

Using Eq. 9, we obtain for the $T_{N}$-odd quantities $A_{T}$ and
$\bar{A}_{T}$ in Eqs. 6 and 7:

$\displaystyle A_{T}={\cal
I}_{ab}\sin(\Delta\delta+\Delta\phi)\leavevmode\nobreak\ \leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \bar{A}_{T}={\cal
I}_{\bar{a}\bar{b}}\sin(\Delta\delta-\Delta\phi)$ (10)

with

$\displaystyle{\cal I}_{ab}$ $\displaystyle\propto$
$\displaystyle\frac{\int_{R}d\Phi\cdot f_{a}f_{b}\cdot V\cdot{\tt sign}({\cal
O}_{\tt CP})}{\int_{R}d\Phi\cdot
f_{a}f_{b}\cdot\left(W+U\cdot\cos(\Delta\delta+\Delta\phi)\right)}\leavevmode\nobreak\
,$ (11)

where $d\Phi$ is the phase-space volume element, $R$ is the phase-space region
of integration and $f_{a},f_{b}$ are the PDF’s of the incoming particles
$a,b$; similarly, for the CC channel, ${\cal I}_{\bar{a}\bar{b}}$ is obtained
by replacing $f_{a}f_{b}\to f_{\bar{a}}f_{\bar{b}}$ and ${\cal O}_{\tt
CP}\to\overline{\cal O}_{\tt CP}$.

Thus, the CP-asymmetry $A_{CP}$ of Eq. 8 is given by:

$\displaystyle A_{CP}=\frac{{\cal I}_{ab}+{\cal
I}_{\bar{a}\bar{b}}}{2}\cos\Delta\delta\sin\Delta\phi+\frac{{\cal
I}_{ab}-{\cal
I}_{\bar{a}\bar{b}}}{2}\sin\Delta\delta\cos\Delta\phi\leavevmode\nobreak\ ,$
(12)

where the 1st term is the conventional $T_{N}$-odd CP-odd term which is
sensitive to CPV at the tree-level (i.e., it is $\propto\cos\Delta\delta$ and
therefore does not require a CP-even phase from FSI), while the (new) 2nd term
represent a ”fake” CP effect (i.e., it does not require a CP-odd phase), which
arises when the initial state is not self-conjugate.

To summarize this section, we have constructed three $T_{N}$-odd asymmetries:
$A_{T}$, $\bar{A}_{T}$ and $A_{CP}$. We find that $A_{CP}$, which is supposed
to be a genuine CP-asymmetry by construction, may, also be “contaminated” by
CP-even phases when the initial state is not CP-symmetric, as can be the case
at the LHC. However, Eq. 12 implies that at the tree-level where
$\Delta\delta=0$ (i.e., in the absence of FSI), all three asymmetries are
$\propto$ CP-odd phase and are thus good measures of CPV, regardless of the
initial state properties.

##  3 The new physics framework

We use an effective field theory (EFT) approach to describe the underlying NP
responsible for CPV in tri-lepton production at the LHC and demonstrate our
strategy for the CPV search in multi-lepton events using the following two
scalar and tensor 4-Fermi effective operators:222Note that the typical bounds
on the natural scale of the operators in Eqs. 13 are $\Lambda\lower
2.15277pt\hbox{$\>\buildrel>\over{\sim}\>$}{\cal O}(1)$ TeV, see [3].

$\displaystyle{\cal
O}_{S}=\left(\bar{\ell}_{R}\ell_{R}\right)\left(\bar{t}_{R}u_{R}\right)\leavevmode\nobreak\
\leavevmode\nobreak\ ,\leavevmode\nobreak\ \leavevmode\nobreak\ {\cal
O}_{T}=\left(\bar{\ell}_{R}\sigma_{\mu\nu}\ell_{R}\right)\left(\bar{t}_{R}\sigma_{\mu\nu}u_{R}\right)$
(13)

and similarly for the $tc\ell\ell$ 4-Fermi terms, where

$\displaystyle{\cal
L}_{tu\ell\ell}^{dim6}=\frac{1}{\Lambda^{2}}\left(f_{S}{\cal O}_{S}+f_{T}{\cal
O}_{T}\right)\leavevmode\nobreak\ ,$ (14)

and ${\cal L}={\cal L}_{SM}+{\cal L}_{tu\ell\ell}^{dim6}$ (we assume that the
two operators have a common underlying scale).

These 4-Fermi interactions can be generated by tree-level exchanges of heavy
scalars and tensors in the underlying heavy theory. Interesting examples are
the scalar leptoquarks $S_{1}$ and $R_{2}$, which transform as $(3,1,-1/3)$
and $(3,2,7/6)$, respectively, under the
$\textrm{SU(3)}\times\textrm{SU}(2)\times\textrm{U}(1)$ SM gauge group. In
particular, tree-level exchanges of $S_{1}$ and $R_{2}$ among the lepton-quark
pairs induce the dimension six scalar and tensor 4-Fermi operators in Eqs. 13,
with a specific proportion of the corresponding Wilson coefficients, obeying
the relation:

$\displaystyle|f_{T}|=\frac{1}{4}|f_{S}|\leavevmode\nobreak\ ,$ (15)

which we have used as our benchmark relation in the numerical study of the CPV
effect below.

The 4-Fermi $tu\ell\ell$ and $tc\ell\ell$ terms with complex Wilson
coefficients can drive tree-level CPV in tri-lepton events at the LHC via the
single top-quark production channels (see Fig. 1):

$\displaystyle ug,gg\to\ell^{+}\ell^{-}t,\leavevmode\nobreak\
\ell^{+}\ell^{-}t+j$ (16)

and the charged conjugate channel (i.e., with $\bar{t}$ in the final state),
where $j={\rm light\leavevmode\nobreak\ jet}$. This channel, followed by the
top decay $t\to Wb\to\ell^{\prime}\nu_{\ell^{\prime}}b$, has interesting
implications also for generic BSM searches of new heavy states around the TeV
scale, which can generate top-leptons 4-Fermi contact terms.[3]

Figure 1: Representative Feynman diagrams for the lowest order single top-
quark + di-lepton production channels $pp\to t\ell^{+}\ell^{-}$ and $pp\to
t\ell^{+}\ell^{-}+j$ ($j$ is a light jet), via the $tu\ell^{+}\ell^{-}$
4-Fermi interaction (marked by a heavy dot).

##  4 Collider signals of CPV in tri-lepton events

As discussed above, we consider CPV in the single-top production channel
depicted in Fig. 1, followed by the top decay to a charged lepton and focusing
henceforward on the $e\mu\mu$ final state:

$\displaystyle pp\to t\mu^{+}\mu^{-}+X\to
e^{+}\mu^{+}\mu^{-}+X\leavevmode\nobreak\ ,$ $\displaystyle
pp\to\bar{t}\mu^{+}\mu^{-}+\bar{X}\to
e^{-}\mu^{+}\mu^{-}+\bar{X}\leavevmode\nobreak\ ,$ (17)

where the NP and CPV phases arise from the dim.6 $tu\mu\mu$ and/or $tc\mu\mu$
4-Fermi interactions.

The dominant SM background to the tri-lepton production channel is from $WZ$
production, $pp\to WZ$, followed by the $W$ and $Z$ decays to charged leptons.
The SM contributions to $pp\to e^{+}\mu^{+}\mu^{-}+X$ from other channels such
as $pp\to t\bar{t},t\bar{t}V,tVV,V+jets$ ($V=Z,W$) are much smaller, see [1].
In particular, there is no interference between the $tu\mu\mu$ NP diagrams
(see Fig. 1) and the SM contributions to the tri-leptons signal, so that the
cross-section is of the form:

$\sigma_{e\mu\mu}(m_{\mu\mu}^{min})=\sigma^{SM}_{e\mu\mu}(m_{\mu\mu}^{min})+\frac{f^{2}}{\Lambda^{4}}\sigma^{NP}_{e\mu\mu}(m_{\mu\mu}^{min})\leavevmode\nobreak\
,$ (18)

where we have defined the $m_{\mu\mu}^{min}$-dependent cumulative cross-
section, selecting events with $m_{\mu\mu}>m_{\mu\mu}^{min}$:

$\displaystyle\sigma_{e\mu\mu}(m_{\mu\mu}^{min})\equiv\sigma_{e\mu\mu}(m_{\mu\mu}>m_{\mu\mu}^{min})=\int_{m_{\mu\mu}\geq
m_{\mu\mu}^{min}}dm_{\mu\mu}\frac{d\sigma_{e\mu\mu}}{dm_{\mu\mu}}\leavevmode\nobreak\
,$ (19)

In particular, as shown below, $m_{\mu\mu}$, which is the invariant mass of
the two muons in the final state, is a useful discriminating parameter.333In
the more general case of $pp\to\ell^{\prime\pm}\ell^{+}\ell^{-}+X$, we have
$m_{\mu\mu}\to m_{\ell\ell}$ and $m_{\ell\ell}$ would be the invariant mass of
the ”none-top” opposite sign same-flavor (OSSF) di-leptons from the underlying
hard process, i.e., of the di-leptons produced from the $tt\ell\ell$ vertex
and not from the top-quark decays.

The CP-violation in our case is, therefore, a pure NP effect, since it arises
from the imaginary part of the interference between the two dim.6 scalar and
the tensor operators, if at least one of the corresponding Wilson coefficients
is complex. In particular, the numerator of the CP-asymmetry $A_{CP}$ (and of
the $T_{N}$-odd asymmetries $A_{T}$ and $\bar{A}_{T}$) is proportional to the
CP-violating part of the cross-section for $pp\to t\mu^{+}\mu^{-}\to
e^{+}\mu^{+}\mu^{-}+X$, which is:

$\displaystyle
d\hat{\sigma}(CPV)\propto\epsilon\left(p_{u_{i}},p_{e^{+}},p_{\mu^{+}},p_{\mu^{-}}\right)\cdot{\rm
Im}\left(f_{S}f_{T}^{\star}\right)\leavevmode\nobreak\ ,$ (20)

and similarly for the CC channel $pp\to\bar{t}\mu^{+}\mu^{-}\to
e^{-}\mu^{-}\mu^{+}+\bar{X}$ by replacing
$\epsilon\left(p_{u_{i}},p_{e^{+}},p_{\mu^{+}},p_{\mu^{-}}\right)$ with
$\epsilon\left(p_{\bar{u}_{i}},p_{e^{-}},p_{\mu^{-}},p_{\mu^{+}}\right)$,
where
$\epsilon\left(p_{1},p_{1},p_{3},p_{4}\right)=\epsilon_{\alpha\beta\gamma\delta}p_{1}^{\alpha}p_{2}^{\beta}p_{3}^{\gamma}p_{4}^{\delta}$
and $\epsilon_{\alpha\beta\gamma\delta}$ is the Levi-Chivita tensor. Guided by
the relation between the scalar and tensor couplings in Eq. 15, which is
motivated by the matching of the EFT framework to the leptoquarks exchanges
scenarios, we set below $|f_{S}|=1$, $|f_{T}|=0.25$ with a maximal CP-odd
phase for the $tu\mu\mu$ and $tc\mu\mu$ operators, so that the CPV coupling in
Eq. 20 is set to:

${\rm Im}\left(f_{S}\cdot f_{T}^{\star}\right)=0.25\leavevmode\nobreak\ .$
(21)

##  5 Results

All three asymmetries $A_{CP}$, $A_{T}$ and $\bar{A}_{T}$ are sensitive to the
di-muons invariant mass, since the SM (which is dominant at low $m_{\mu\mu}$)
contributes to the denominators while the CPV NP (which is dominant for high
$m_{\mu\mu}$) contributes to the numerators. We therefore exploit this
$m_{\mu\mu}$ dependence to gain sensitivity to CPV.

Our main results are summarized in Table 2 and Fig. 2: in Fig. 2 we show the
dependence of $A_{CP}$ on $m^{min}_{\ell\ell}$ and in Table 2 we give the
resulting CP-violating and $T_{N}$-odd asymmetries for
$m^{min}_{\ell\ell}=400$ GeV. Results are shown for both the $ug$-fusion and
$cg$-fusion production channels, for $\Lambda=1,2$ TeV and for ${\rm
Im}\left(f_{S}f_{T}^{\star}\right)=0.25$. The SM background from $pp\to
ZW^{\pm}+X$ is included.

Table 2: The expected $T_{N}$-odd and CP asymmetries in tri-lepton events, $pp\to\ell^{\prime\pm}\ell^{+}\ell^{-}+X$, see text.  |  $ug$-fusion: $\Lambda=1(2)$ TeV |  $cg$-fusion: $\Lambda=1(2)$ TeV
---|---|---
$A_{CP}$ | 11.1(7.9)% | 3.9(0.7)%
$A_{T}$ | 16.4(13.5)% | 3.1(0.5)%
$\bar{A}_{T}$ | -5.8(-2.3)% | -4.7(-1.0)%
Figure 2: The expected CP-asymmetry $A_{CP}$, as a function of the invariant
mass cut on the same-flavor di-leptons, see text. The error bars represent the
expected statistical uncertainty with an integrated luminosity of $1000(3000)$
fb-1 for the ug-fusion(cg-fusion) case.

##  6 Summary

Multi-leptons signals provide an excellent and rich testing ground for NP
associated with flavor physics, lepton flavor universality, CP-violation etc…

In this work we have constructed useful CP-asymmetries for measuring CPV
effects in multi-lepton events, introducing a new modification to the classic
formula for tree-level CPV in scattering and decay processes, which takes into
account the effect of an asymmetric initial state and which is particularly
useful for CP studies at the LHC.

Our asymmetries have several new and unique features, that make them
particularly useful for searching for CPV at high-energy colliders:

  * •

They use only multi-lepton final states as probes, which makes them
experimentally highly distinctive.

  * •

They are based on simple kinematic observables that only require the
reconstruction of the relatively easily-identifiable charged-lepton momenta.

  * •

They can be generated by tree-level CPV underlying physics, making them very
sensitive to NP.

  * •

They are generic, since they can probe a wide range of underlying NP.

We find that ${\cal O}(10\%)$ CP-asymmetries may be generated with new CPV
TeV-scale NP, whereas the expected SM background for CPV in multi-lepton
events is at the sub-% level. Finally, with the NP considered, we expect
${\cal O}(10000)$ high-$p_{T}$ tri-lepton events with an integrated luminosity
of ${\cal O}(1000)$ fb-1, if TeV-scale NP is generating the $tu\ell\ell$
4-Fermi interactions. A more detailed study relevant to this talk is given in
[1].

## Acknowledgments

I wish to thank the organizers of the Moriond 2023 conference for their
invitation to present this talk and for the welcoming atmosphere in this
meeting.

## References

## References

  * [1] Y. Afik, S. Bar-Shalom, K. Pal, A. Soni, J. Wudka, arxiv: 2212.09433 (hep-ph), 2022.
  * [2] Y. Afik, S. Bar-Shalom, K. Pal, A. Soni, J. Wudka, Nucl. Phys. B 980, 115849 (2022).
  * [3] Y. Afik, S. Bar-Shalom, A. Soni, J. Wudka, Phys. Rev. D 103, 075031 (2021).
  * [4] Y. Afik, S. Bar-Shalom, J. Cohen, A. Soni, J. Wudka, Phys. Lett. B 811, 135908 (2020).
  * [5] Y. Afik, S. Bar-Shalom, J. Cohen, Y. Rozen, Phys. Lett. B 807, 135541 (2020).
  * [6] Y. Afik, J. Cohen, E. Gozani, E. Kajomovitz, Y. Rozen, JHEP 08, 056 (2018).
  * [7] D. Atwood, S. Bar-Shalom, G. Eilam, A. Soni, Phys. Rept. 347, 1 (2001).
  * [8] O. Antipin and G. Valencia, Phys. Rev. D 79, 013013 (2009).
  * [9] S.K. Gupta, A.S. Mete and G. Valencia, Phys. Rev. D 80, 034013 (2009).
  * [10] T. Han, Y. Li, Phys. Lett. B 683, 278 (2010).
  * [11] S.K. Gupta and G. Valencia, Phys. Rev. D 81, 034013 (2010).
  * [12] A. Hayreter and G. Valencia, Phys. Rev. D 93, 014020 (2016).
  * [13] A. Tiwari and S.K. Gupta, Adv. High Energy Phys. 2021, 6676930 (2021).
  * [14] W. Dekens, J. de Vries, M. Jung and K.K. Vos, JHEP 01, 069 (2019).
  * [15] I. Bigaran and R.R. Volkas, Phys. Rev. D 105, 015002 (2022).