# Constraints on conformal ultralight dark matter couplings from the European
Pulsar Timing Array

Clemente<EMAIL_ADDRESS>, Adrien<EMAIL_ADDRESS>Enrico Barausse(1,2), Boris Goncharov(3,4,5,6), Diana López Nacir(7,8) , Diego
Blas(9,10), Lijing Shao(11,12), J. Antoniadis(13,14), D. J. Champion(14) , I.
Cognard(15,16), L. Guillemot(15,16), H.Hu(14), M. Keith(17), M. Kramer(14,17),
K. Liu(14,18), D. Perrodin(19), S. A. Sanidas(17), G. Theureau(15,16,20)
0000-0002-0817-2830 0000-0002-4803-2998 0000-0001-6499-6263
0000-0003-3189-5807 0000-0003-4398-1147 0000-0002-1334-8853
0000-0003-4453-3776 0000-0003-1361-7723 0000-0002-1775-9692
0000-0002-9049-8716 0000-0002-3407-8071 0000-0001-5567-5492
0000-0002-4175-2271 0000-0002-2953-7376 0000-0002-1806-2483
0000-0002-5956-5546 0000-0002-3649-276X

###### Abstract

Millisecond pulsars are extremely precise celestial clocks: as they rotate,
the beamed radio waves emitted along the axis of their magnetic field can be
detected with radio telescopes, which allows for tracking subtle changes in
the pulsars’ rotation periods. A possible effect on the period of a pulsar is
given by a potential coupling to dark matter, in cases where it is modeled
with an “ultralight” scalar field. In this paper, we consider a universal
conformal coupling of the dark matter scalar to gravity, which in turn
mediates an effective coupling between pulsars and dark matter. If the dark
matter scalar field is changing in time, as expected in the Milky Way, this
effective coupling produces a periodic modulation of the pulsar rotational
frequency. By studying the time series of observed radio pulses collected by
the European Pulsar Timing Array experiment, we present constraints on the
coupling of dark matter, improving on existing bounds. These bounds can also
be regarded as constraints on the parameters of scalar-tensor theories of the
Fierz-Jordan-Brans-Dicke and Damour-Esposito-Farèse types in the presence of a
(light) mass potential term.

##  1 Introduction

Elucidating the nature of Dark Matter (DM) remains as one of the most pressing
questions of modern physics. The widely accepted Cold Dark Matter (CDM)
paradigm successfully explains numerous aspects of the Universe’s large-scale
structure but encounters difficulties in predicting some observations on
scales smaller than approximately a kiloparsec ($\sim$ kpc). Notably,
observations indicate a flat density profile in the inner regions of galaxies,
contradicting the pure CDM prediction of a steep power-law-like behavior
(cusp-core problem) [1, 2, 3]. Additionally, known challenges arise from the
mismatch between the observed and expected number of satellites of the Milky
Way (MW) (missing satellite problem) [4, 5], and simulations based on
$\Lambda\text{CDM}$ theory suggest that the most massive subhaloes of the MW
would be too dense not to host bright satellites (too-big-to-fail problem)
[6]. Although these issues may be mitigated by considering baryonic feedback
mechanisms, such as supernova feedback [7, 8, 9, 10, 11, 12], another
possibility is to assume that DM is an ultralight scalar field (with mass
$m\sim 10^{-22}~{}\text{eV}$) with negligible self-interactions [13, 14]. In
this scenario, the de Broglie wavelength of the scalar field in galaxies can
reach $\sim\text{kpc}$, suppressing power on smaller scales, while retaining
all the successes of CDM at large scales. Furthermore, the presence of
ultralight scalars is also motivated from a theoretical standpoint by string
theory [15, 16, 17]. In this context, the mass range can be much broader,
which motivates considering very light bosons in a large span of ultra-light
masses (including, but not restricted to, $m\sim 10^{-22}\,$eV) as natural
candidates for theories beyond the Standard Model.

Numerous studies have been conducted to investigate the existence of ultra-
light DM (ULDM). Among them, the integrated Sachs-Wolfe effect on the Cosmic
Microwave Background (CMB) anisotropies excludes masses $m\lesssim
10^{-24}~{}\text{eV}$ [18]. Lyman-$\alpha$ observations strongly suggests a
lower limit of $m\gtrsim 10^{-21}~{}\text{eV}$ for ultra-light candidates
accounting for more than $\sim 30\%$ of the DM [19, 20, 21, 22, 23]. However,
the susceptibility of non-CMB constraints to uncertainties in the modeling of
small-scale structure properties [24, 25] emphasizes the importance of
complementary and independent investigations. In this context, the rotation
curves of well-resolved nearby galaxies also disfavor masses $m\lesssim
10^{-21}\,$eV [26]. In addition, measurements of stellar orbit kinematics in
ultra-faint dwarf (UFD) galaxies may potentially constrain the scalar field
mass to be $m\gtrsim 10^{-19}~{}\text{eV}$, although this remains a topic of
ongoing debate [27, 28]. Dwarf galaxies can also be used to set robust bounds
$m\gtrsim 10^{-22}~{}\text{eV}$ [29, 14]. Given these observations, the
current lore is to consider that ULDM of mass below $m\sim 10^{-22}\,$eV can
not constitute 100% of the dark matter, but masses below these bound are
certainly possible if they constitute a fraction of the dark matter, see Ref.
[30]. This possibility is natural in the case of the axiverse, where several
ULDM candidates coexist at low masses [31] and also for ultra-low mass
particles that may be cosmologically produced to significant values, see e.g.
Ref. [32].

A completely independent method to probe these small masses was suggested in
Ref. [33], where it was pointed out that the oscillating _gravitational
potential_ induced by the presence of ULDM influences the light travel time of
radio signals emitted by pulsars. Pulsar Timing Array (PTA) experiments [34,
35, 36, 37, 38, 39] rely on the exquisite predictability of the millisecond
pulsar (MSP) spin periods behavior. Each time a MSP magnetic field axis points
towards Earth, radio waves are observed by radio telescopes. After measuring
and modelling consecutive pulses in decade-long observational campaigns, PTAs
search for signals of physical effect that affect all of the observed pulsars,
including the ULDM signal. Based on this principle, previous PTA searches have
established 95% upper limits on the local energy density of ULDM, reaching
$\rho\lesssim 0.15$ GeV/cm3 in the mass range $10^{-24.0}~{}\text{eV}\lesssim
m\lesssim 10^{-23.7}~{}\text{eV}$ [40, 41].

The ideas of Ref. [33] rely on the universal gravitational coupling of DM to
ordinary matter. However, ULDM may also be coupled to the Standard Model
fields [42, 43]. Indeed, a natural possibility that respects the weak
equivalence principle is that ULDM may be universally (conformally) coupled to
gravity, or (equivalently, in the Einstein frame) to the Standard Model. This
universal coupling, together with the oscillations of the scalar field in the
MW, would give rise to periodic orbital perturbations in binary systems, which
allows to place constraints on the model [44, 45, 46]. In this context, ULDM
may be regarded as a scalar-tensor theory of the Fierz-Jordan-Brans-Dicke [47,
48, 49, 50] or Damour-Esposito-Farèse type [51, 52], in the presence of a
(light) mass potential term [53]. As a result of the strong gravitational
fields active inside neutron stars, the conformal coupling to gravity gives
rise to an effective (gravity-mediated) interaction between neutron stars (and
thus pulsars) and the scalar ULDM field [54, 55, 56, 57, 51, 58]. This effect,
known as Nordtvedt effect [54], violates the strong equivalence principle, and
it has long been constrained with binary pulsar data [59, 60, 61, 62, 63, 64,
65, 66].

In a recent companion paper [67], following former studies [68, 56, 58], two
of us computed the effect of this effective interaction on the rotational
period of an isolated pulsar. More specifically, Ref. [67] found that the
effective coupling between ULDM and pulsars produces a change in the moment of
inertia (and therefore on the rotational period) of the pulsar. This change is
proportional to the rate of change (in time) of the scalar field, which – as
mentioned above – is expected to oscillate in the MW. Deriving precise
constraints from current observations on the conformal ULDM coupling based on
this effect is the main objective of this work.

This work is structured as follows. In Section 2, we define the Lagrangian of
our theory, and we briefly review some general features of ULDM relevant to
our analysis. Section 3 will be devoted to a detailed description of the
dataset and the procedure used to carry out the analysis. Finally, in Section
4, we will show how a non-minimally coupled ULDM candidate affects the spin
frequency of MSPs. We constrain the coupling strength, resulting in bounds
which are several orders of magnitude better than what obtained from Cassini
tests of General Relativity [69, 44] or from the pulsar in a triple stellar
system [62, 63, 64] in the mass region which PTAs are sensitive to. Our
conclusions are presented in Section 5. The appendices are devoted to
technical details and plots.

##  2 Conformally coupled ultralight dark matter

Light scalar and pseudoscalar fields emerge naturally from string theory and
from theories with pseudo-Goldstone bosons (as the axions introduced to tackle
the strong CP problem) [70, 31, 71]. These fields are, in principle, coupled
to Standard Model fields. Hence it is natural to consider _non-minimally
coupled_ scenarios when exploring their phenomenology. In this work, we
consider a scalar field $\phi$ with mass $m$. The action for this field in the
Einstein frame is given by111Note that our scalar field $\phi$ is not
canonically normalized, but appears in the action multiplied by
$M_{\text{P}}$. This convention makes comparisons with gravitational phenomena
more straightforward.:

$S=M_{\text{P}}^{2}\int\mathrm{d}^{4}x\sqrt{-g}\,\bigg{[}\frac{R}{2}-g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi+m^{2}\phi^{2}\bigg{]}+S_{m}[\psi_{m},\tilde{g}_{\mu\nu}]\;,$
(2.1)

where $M_{\text{P}}$ is the Planck mass and the matter action $S_{m}$ includes
a universal conformal coupling of the scalar field to the matter content
$\psi_{m}$ through the (Jordan) effective metric
$\tilde{g}_{\mu\nu}=\mathcal{A}^{2}(\phi)g_{\mu\nu}$, normalized such that
$\mathcal{A}^{2}(0)=1$. Note that when re-expressed in the Jordan frame, this
direct coupling to matter disappears and is replaced by a conformal coupling
to gravity (i.e., to the Ricci scalar). In other words, this model satisfies
the weak equivalence principle. In particular, the free-fall motion of non-
gravitating objects is universal (i.e., independent of the body composition)
[72]. However, in strongly gravitating objects (such as neutron stars), there
appears an effective (tensor-mediated) coupling between matter and the scalar
field.

As a first model, we consider the Fierz-Jordan-Brans-Dicke (FJBD) theory [47,
48, 49, 50], in which the conformal coupling is linear:

$\mathcal{A}(\phi)=e^{\alpha\phi}\sim 1+\alpha\,\phi.$ (2.2)

The constant scalar coupling $\alpha$ is constrained by several observations,
notably by the Cassini mission to the level of $\alpha^{2}\lesssim 10^{-5}$
[69] and by the triple system PSR J0337+1715 to the level of
$\alpha^{2}\lesssim 4\times 10^{-6}$ [64]. A simple way to evade this
constraint is to consider masses generating a Yukawa suppression at scales of
order of the typical distances probed by these systems [53]. Since our focus
is on the PTA, and hence to DM masses with Compton wavelength larger than
$\sim 10^{3}$ AU, this constraint applies to the models we explore222The
center of mass of the inner binary in the triple system PSR J0337+1715
completes a rotation around the center of mass of the entire system in about
327 days [64], which corresponds to a distance of $\sim 0.9$ AU..

As a second model, we also consider the Damour-Esposito-Farèse (DEF) gravity
theory [51, 52], where the conformal coupling is quadratic:

$\mathcal{A}(\phi)=e^{\beta\phi^{2}/2}.$ (2.3)

Notice that we have chosen to set the linear coupling of the field in
$\mathcal{A(\phi)}$ to zero, so that FJBD cannot be recovered as a special
case of DEF, rather the two theories are part of a more general theory where
both linear and quadratic couplings are nonzero. The absence of significant
deviations from the General Relativity (GR) predictions in binary pulsar data
requires $\beta\gtrsim-4.3$ (depending on the exact equation of state (EoS)
for the neutron star model) in order to avoid non-perturbative spontaneous
scalarization phenomena [52, 73], while the value $\phi_{0}$ of the scalar
field on cosmological scales is constrained by the Cassini (pulsar in a triple
stellar system) bound to the level $(\beta\phi_{0})^{2}\lesssim 10^{-5}$ and
by the pulsar in a triple stellar system to $(\beta\phi_{0})^{2}\lesssim
4\times 10^{-6}$.

The energy momentum tensor of the scalar field sourcing the metric
$g_{\mu\nu}$ on cosmological scales follows, in the Einstein frame, from
(2.1):

$T_{\mu\nu}=M_{\text{P}}^{2}\left[2\partial_{\mu}\phi\partial_{\nu}\phi-
g_{\mu\nu}\left(\left(\partial\phi\right)^{2}-m^{2}\phi^{2}\right)\right].$
(2.4)

In the Jordan frame $\tilde{T}_{\mu\nu}=\mathcal{A}^{-2}(\phi)T_{\mu\nu}$, it
reduces to

$\tilde{T}_{\mu\nu}=\mathcal{A}^{-2}(\phi)T_{\mu\nu}\simeq\left(1-2\,\alpha\left(\phi\right)\right)T_{\mu\nu},$
(2.5)

where the effective scalar coupling is333In order to avoid confusion, we
stress that $\alpha(\phi)$ is generally different from $\alpha$, but reduces
to it in FJBD theory. Binary pulsars have constrained $|\alpha(\phi)|$ to be
$\lesssim 10^{-2}$ for neutron stars [66].
$\alpha(\phi)=\mathrm{d}\log\mathcal{A}/\mathrm{d}\phi$ and we work in the
limit $\alpha(\phi)\ll 1$. The mass of the scalar field can be as small as
$m\sim 10^{-22}~{}\text{eV}$, and still be a viable very light candidate for
CDM. We will refer to models of masses not far from this limit as ULDM. In
these models, the distance between particles is much smaller than the
corresponding de Broglie wavelength, which implies that they can be treated as
a classical superposition of free waves with dynamical properties generated by
galactic dynamics. For the MW, this superposition has a very small dispersion
velocity ($\sigma_{\phi}\sim 10^{-3}$), and therefore the ULDM field can be
described as a standing wave444Close to the galactic center, the distribution
may condense into a different configuration known as “soliton” or “bose star”
[24]. We will not deal with this situation since the pulsars used by PTA are
far from the galactic center.:

$\displaystyle\phi(\boldsymbol{x},t)=A(\boldsymbol{x})\cos(mt+\theta(\boldsymbol{x})),$
(2.6)

with $A(\boldsymbol{x})$ determined by555Notice that, as compared to field
theory conventions, there are factors of 2 of difference, arising from the
non-canonical normalization of the scalar field in Eq. (2.1).

$\displaystyle\tilde{\rho}_{\phi}\approx\rho_{\phi}=m^{2}M_{\text{P}}^{2}A(\boldsymbol{x})^{2}=\rho\,\hat{\phi}(\boldsymbol{x})^{2},$
(2.7)

where we have used Eq. (2.5) to relate the stress-energy momentum tensor in
the Einstein frame to the one in the Jordan frame, and we have neglected the
terms $\sim(\partial_{i}\phi)^{2}$ which are suppressed by $v_{\phi}^{2}$.
Here, $\rho$ is the average density of the scalar field and
$\hat{\phi}(\boldsymbol{x})$ is a stochastic parameter extracted from the
Rayleigh distribution ($P(\hat{\phi}^{2})=e^{-\hat{\phi}^{2}}$) [74]. This
parametrization takes into account the fact the ULDM configurations are built
by the superposition of several waves of random phases that interfere. For the
average density, we take $\rho=\rho_{\text{DM}}=0.4~{}\text{GeV/cm}^{3}$ as a
benchmark value for the average DM density expected at the Sun position in the
Milky Way [75]. As commented in the introduction, masses $m\lesssim
10^{-22}~{}\text{eV}$ are strongly disfavored to constitute all the dark
matter in our Galaxy. As a result, for these masses one needs to focus on
scenarios where ULDM is a fraction $f_{\text{DM}}\equiv\rho/\rho_{\rm DM}<1$.
The factor $\hat{\phi}$ appears from the interference caused by the wavelike
nature of the scalar field. It reproduces the expected random local value from
the superposition of the waves, and makes the scalar field density
$\tilde{\rho}_{\phi}$ approach $\rho$ when averaging over timescales longer
than the ULDM coherence time

$\tau_{c}\sim\frac{2}{mv^{2}}=2\times
10^{5}~{}\text{yr}\left(\frac{m}{10^{-22}~{}\text{eV}}\right),$ (2.8)

or, equivalently, on a length scale larger than the ULDM coherence length

$l_{c}\sim\frac{1}{mv}\sim v\cdot\tau_{c}\sim
0.4~{}\text{kpc}\left(\frac{10^{-22}~{}\text{eV}}{m}\right).$ (2.9)

We have conveniently normalized the mass to values relevant to PTA
observations.
In the Jordan frame, an oscillating scalar field, such as the one presented in
Eq. (2.6), induces a temporal variation of Newton’s constant. Its measured
value is [51]:

$G=\big{(}8\pi
M_{\text{P}}^{2}\big{)}^{-1}\mathcal{A}^{2}(\phi)\big{(}1+\alpha^{2}(\phi)\big{)},$
(2.10)

where $\alpha(\phi)=\mathcal{A}^{\prime}(\phi)/\mathcal{A}(\phi)$. In order to
use this property of scalar-tensor theories, we follow Ref. [67] and neglect
the scalar field mass on the typical length scale of a pulsar, a valid
approximation given our mass range for ULDM. In turn, a variation of the local
gravitational constant modifies the gravitational mass and the radius of the
neutron star [51]. This dependence is encoded in the sensitivities, which
represent the rate of change of these parameters with respect to changes in
the scalar field [72].

To explore this effect, one can recall that the conservation of angular
momentum $J$ relates the changes in the moment of inertia $I$ of the neutron
star (depending on the local value of $G$) to the observed pulse frequency
$\Omega_{\mathrm{obs}}$ through the relation $J=I\Omega_{\mathrm{obs}}$.
Particular attention has to be paid to the frame used for the definition of
the angular momentum $J$, as the latter is only conserved in the Einstein-
frame (see Ref. [67] for more details). We use the code presented in Ref. [67]
to compute the angular momentum sensitivity, defined by:

$s_{I}=-\frac{1}{2\alpha(\phi)}\frac{\mathrm{d}\ln
I}{\mathrm{d}\phi}\bigg{|}_{N,J}=\frac{1}{2\alpha(\phi)}\frac{\mathrm{d}\ln\Omega_{\mathrm{obs}}}{\mathrm{d}\phi}\bigg{|}_{N,J},$
(2.11)

where the pulsar’s baryon number $N$ and Einstein-frame angular momentum $J$
are kept constant. With this sensitivity at hand, a change in the scalar field
value can be directly related to a change in the frequency of the pulsar, and
consequently to a change in the pulsar’s pulse time of arrival (TOA). We will
use this fact in Section 4 to constrain ULDM models.

##  3 Dataset and methodology

In this work, we analyze the second data release (DR2) [76] of the European
Pulsar Timing Array collaboration (EPTA) [37, 77] . In particular, we use the
EPTA-DR2Full dataset666The dataset can be found at
https://gitlab.in2p3.fr/epta/epta-dr2., consisting of a 24.7 years collection
of TOAs of radio pulses of 25 millisecond pulsars, surveyed with an
approximately biweekly cadence777The cadence is non-uniform. EPTA combines
observations from several telescopes, so sometimes the EPTA-wide cadence can
be much shorter. and observed by five telescopes located in France, Germany,
Italy, the Netherlands, and the United Kingdom. The relation between the time
of emission of a radio pulse and its TOA at the Solar System Barycentre (SSB)
is encoded in a pulsar-specific timing model [78], which takes into account
the pulsar spin down, its position and proper motion, the motion around a
binary companion, etc.. Any deviations between the predicted TOAs and the
actual measurements, referred to as timing residuals, may include
contributions from various sources of noise, including stochastic dispersion
measure fluctuations and irregularities in pulsar rotation [79, 80]. However,
these residuals might also be indicative of processes of astrophysical
significance, c.f. the recent evidence of a stochastic gravitational wave
background (SGWB) in the data of various PTA experiments [81, 82, 83, 84].

Because the EPTA-DR2Full dataset does not yield a strong evidence in favor of
the hallmark Hellings-Downs (HD) [85] inter-pulsar correlation (in contrast to
the 10.3 yr dataset [76, 80, 86, 87]), we only account for possible
contributions from the SGWB via a PTA-wide spatially-uncorrelated but
temporally-correlated noise term, characterized by an amplitude
$A_{\text{GWB}}$ and a spectral index $\gamma_{\text{GWB}}$ (see Table 1 for
more details). Such a signal appears as a precursor to the SGWB [88, 89, 90],
because of the stronger autocorrelation of the Hellings-Downs process.

We utilize Bayesian inference to detect the ULDM-induced deterministic888The
signal induced by ULDM is deterministic, as opposed to the stochastic nature
of e.g. the common red noise process describing the SGWB. signal, while
simultaneously marginalizing over timing model parameters [76] and accounting
for all known sources of noise in the data [80]. Given the model parameters
$\theta$, the likelihood function for the timing residuals, denoted as
$\mathcal{L}(\delta t|\theta)$, is represented as [91, 92, 93, 94, 95]:

$\ln\mathcal{L}(\delta t|\theta)\propto-\frac{1}{2}(\delta
t-\mu)^{\text{T}}C^{-1}(\delta t-\mu).$ (3.1)

In this time-domain Gaussian likelihood, $\delta t$ is a vector with dimension
corresponding to the number of observations. The deterministic ULDM
contribution, which we derive below in Eqs. (4.6) and (4.10), is taken into
account in $\mu$, which includes contributions from both the timing model and
noise processes, as analyzed in Ref. [80]. Temporally-uncorrelated “white”
noise and other sources of uncertainty in TOA measurements are factored in the
diagonal components of the covariance matrix $C$. Off-diagonal elements of the
matrix $C$ could, in principle, contain contributions from temporally
correlated “red” noise; yet they are more commonly incorporated into $\mu$ for
computational efficiency, following the approach described in Refs. [92, 93].
The priors $\pi(\theta)$ of the parameters used for the search are described
in Table 1. Parameter estimations are carried out by evaluating posterior
distributions, denoted as $\mathcal{P}(\theta|\delta
t)\propto\mathcal{L}(\delta t|\theta)\pi(\theta)$, produced with the Parallel-
Tempering-Markov-Chain Monte-Carlo sampler [96] implemented in enterprise [94]
and enterprise_extensions [95], using the PTArcade [97, 98] wrapper and
adapting it to the EPTA dataset.

In this work, we consider the effect of the scalar field on the TOAs, and we
constrain the conformal coupling by looking at its effect on the timing
residuals. Since the duration $T_{\text{obs}}=24.7~{}\text{yr}$ of the EPTA-
DR2Full dataset is much shorter than the coherence time in Eq. (2.8) for the
ULDM considered here, different coherence patches with characteristic
dimension $l_{c}$ will have different $\rho_{\phi}$. However, notice that if
$m$ is sufficiently small so that $l_{c}>R$, where $R$ is the characteristic
radius probed by Galactic rotation curves, we are really observing one single
patch of ULDM. We refer to this case as correlated scenario.

Based on these premises, we thus distinguish three different regimes [40],
according to the interplay between the mass of the ULDM candidate and the
typical interpulsar separation. In the uncorrelated regime, each of the
pulsars timed by the EPTA resides in a different coherent patch. Therefore,
each pulsar is characterized by its own $\hat{\phi}(\boldsymbol{x})$
parameter. As the average interpulsar distance is $\mathcal{O}(\text{kpc})$,
the uncorrelated approximation holds for ULDM masses $m\gtrsim 5\times
10^{-23}~{}\text{eV}$. For masses $m\lesssim 2\times 10^{-24}~{}\text{eV}$
(correlated scenario), the coherence length described by Eq. (2.9) encompasses
the inner Galacto-centric $\sim\,20\,\text{kpc}$, which is the benchmark area
examined by the most accurate measurements of MW rotation curves [75], from
which the local DM abundance is inferred. Therefore, as the kinematic tests of
the DM halo explore the same coherence patch that hosts all the EPTA pulsars,
we can safely absorb the common parameter $\hat{\phi}(\boldsymbol{x})$ into
the measured value of the local ULDM abundance. Equation (2.6) then reads:

$\phi(\boldsymbol{x},t)=\frac{\sqrt{\rho}}{mM_{\text{P}}}\cos(mt+\theta(\boldsymbol{x})),$
(3.2)

with $\rho$ representing the value of the scalar field density $\rho_{\phi}$
in our Galaxy. Finally, for masses $2\times 10^{-24}\text{eV}\lesssim
m\lesssim 5\times 10^{-23}\text{eV}$, one ULDM coherence patch can encompass
all pulsars, but does not reach the typical radius explored by rotation
curves. In this pulsar-correlated scenario, the stochastic parameter
$\hat{\phi}(\boldsymbol{x})$ is common across all the pulsars. However,
estimates of DM density derived from rotation curves average over different
coherence patches. Hence, we maintain $\hat{\phi}(\boldsymbol{x})$ as an
independent parameter and marginalize over it, so that the constraints on
$\rho$ derived from the following analysis can be compared to the density
estimated through rotation curves.

Regardless of the scenario, we always draw one phase parameter
$\theta(\boldsymbol{x})$ per pulsar. This phase encodes the uncertainty on
current pulsar distance measurements [99, 40, 100] (see below Eq. (4.3) in
Sec. 4 for more details).

We focus on the ULDM mass range
$m\in[10^{-24}~{}\text{eV},10^{-21}~{}\text{eV}]$, since this is the interval
where PTA constraints are the most compelling. Notice that the low-mass end
corresponds to a frequency of $f_{\text{low}}\sim 2.4\times
10^{-10}~{}\text{Hz}$, which is far below the inverse of EPTA-DR2Full
observation length $f_{\text{obs}}=1/T_{\text{obs}}\sim 1.3~{}\text{nHz}$. In
this regime, the ULDM-induced signal (see Eq. (4.6) later) can still be
expanded in powers of ($mt$) [42]. The first terms in the expansion are
degenerate with the simultaneous fitting of pulsar spin frequency derivatives
[101, 102, 103]; therefore, the lowest-order term that PTAs are sensitive to
is $(mt)^{3}$. Although this introduces a sensitivity loss - which is also
confirmed e.g. by Fig. 2 \- the ULDM-induced signal amplitude is inversely
proportional to the particle mass (see Eqs. (4.6)- (4.5)). Therefore, we can
still provide significant bounds at $f\lesssim 1/T_{\text{obs}}$ [104]. The
upper end of the ULDM mass range, instead, corresponds to $f_{\text{up}}\sim
2.4\times 10^{-7}~{}\text{Hz}$, which is somewhat lower than the observational
cadence of $f\sim 1/2\,\text{weeks}\sim 8.3\times 10^{-7}~{}\text{Hz}$, after
which the PTA data are not sensitive. Finding no evidence for a signal in the
examined mass range, we present 95% upper limits in the following sections.

##  4 Results

In this section, we apply the theoretical framework laid down in Sec. 2 to
constrain the FJBD and the DEF conformal couplings.

###  4.1 FJBD conformal coupling

Let us focus on FJBD theory. In this case, by inspecting Eq. (2.2) and
recalling the definition of $\alpha(\phi)$, we find $\alpha(\phi)=\alpha$.
Moreover, the numerical analysis carried out in Ref. [67] shows that the
angular momentum sensitivity $s_{I}$ has a very weak dependence on the scalar
coupling $\alpha$ that we neglect in this analysis. Depending on the context,
we conveniently write $s_{I}(\alpha,M)$ as $s_{I}(M)$ or $s_{I}$ to avoid
cluttering the notation.

For small scalar fluctuations such as the ones considered in this article, it
follows from Eq. (2.11) that

$\displaystyle\Omega_{\text{obs}}(t)=\bar{\Omega}\left(1+2\alpha
s_{I}\delta\phi(t)\right)=\bar{\Omega}\left(1+2\alpha
s_{I}\frac{\sqrt{\rho}}{M_{\text{P}}m}\hat{\phi}(\boldsymbol{x})\cos(mt+\theta(\boldsymbol{x}))\right),$
(4.1)

where we used Eq. (2.6) and we denoted the spin frequency of the pulsar in the
absence of the scalar field by $\bar{\Omega}$. Notice that Eq. (4.1) describes
the correct frequency shift only under the assumption that the oscillating
timescale $t_{\text{osc}}\sim 1/m$ is much longer than the timescale on which
a neutron star adjusts its internal structure ($t_{\text{int}}$). This is a
reasonable assumption. For instance, the Vela Pulsar shows a fast core-crust
coupling with a timescale $t_{\text{int}}\sim 10~{}\text{s}$ [105], which
would be larger than $t_{\text{osc}}$ only for ULDM masses $m\gtrsim
10^{-16}~{}\text{eV}$, which are not discussed here (see Sec. 3).

To find the TOA change induced by the scalar field, we write Eq. (2.11) as

$\frac{\delta\Omega_{\text{obs}}}{\bar{\Omega}}=2\alpha
s_{I}\delta\phi=-\frac{\delta P}{\bar{P}},$ (4.2)

where $\bar{P}$ is the pulsar period in the absence of the scalar field. The
total timing residual $\Delta t(t)$ after a time $t$ is the integral of
infinitesimal period variations over time:

$\displaystyle\Delta
t(t)=-\int\mathrm{d}t\frac{\delta\Omega_{\text{obs}}}{\bar{\Omega}}=2\alpha
s_{I}\left.\frac{\sqrt{\rho}}{M_{\text{P}}m^{2}}\hat{\phi}(\boldsymbol{x})\sin(mt+\theta(\boldsymbol{x}))\right|_{t_{\text{start
}}-\frac{d}{c}}^{t_{\text{end }}-\frac{d}{c}},$ (4.3)

where we have highlighted the dependence on the retarded time
$t_{\textit{i}}-d/c$, with $i=\textit{start, end}$ and $d$ referring to the
Earth-pulsar distance. As mentioned before, the pulsar distance can be re-
absorbed in a redefinition of the phase
$\theta(\boldsymbol{x})\rightarrow\theta(\boldsymbol{x})+md/c$. Present
uncertainties on the pulsar distances are $\mathcal{O}(0.1\div
1)~{}\text{kpc}$ [99], implying that this redefinition gives rise to an
effective pulsar-dependent random phase. Therefore, as mentioned in Sec. 3, we
treat $\theta(\boldsymbol{x})$ as a pulsar-specific random parameter, and we
neglect the distance in the retarded time from now on.

In analogy with the ULDM search results in Ref. [40], where the timing
residual for the model considered was written as

$\delta
t_{\text{DM}}=\frac{\Psi_{\text{c}}}{2m}[\hat{\phi}^{2}_{\text{E}}\sin{(2mt+\theta_{\text{E}})}-\hat{\phi}^{2}_{\text{P}}\sin{(2mt+\theta_{\text{P}})}],$
(4.4)

with $E,P$ labeling respectively the Earth and the pulsar999The stochastic
parameter $\hat{\phi}_{\text{P}}$ corresponds to $\hat{\phi}$ in our notation.
, we can define an effective amplitude

$\Psi=2\alpha\frac{\sqrt{\rho}}{M_{\text{P}}m},$ (4.5)

such that

$\Delta
t(t)=\left.\frac{\Psi}{m}s_{I}\hat{\phi}(\boldsymbol{x})\sin(mt+\theta(\boldsymbol{x}))\right|_{t_{\text{start
}}}^{t_{\text{end }}}.$ (4.6)

This form helps to understand what a sensible prior on $\Psi$ may be. In fact,
noticing that Eqs. (4.4) and (4.6) have the same form (differing only for the
presence of the Earth term and some $\mathcal{O}(1)$ factors), we can use the
same prior for $\Psi$ and $\Psi_{c}$, i.e.
$\text{Log}_{10}-\text{Uniform}\left[-20,-12\right]$. In other words, the
similarity between the two equations shows that what we are really testing is
whether the PTA data can constrain the presence of a sinusoidal signal.

As previously stated, the sensitivity $s_{I}$ as a function of the pulsar mass
is computed from a fit to the models of Ref. [67]. In particular, we
consistently utilize the pulsar gravitational mass as a parameter of the fit
instead of the inertial one, because the former is the value measured by
experiments. We implement the pulsar masses in the analysis in the following
way:

  * •

if a pulsar mass is determined from other experiments as $M\pm\delta M$, we
draw the mass from a normal distribution centered on $M$, with uncertainty
$\delta M$ and truncated for masses below $M_{\text{min}}=0.2\,M_{\odot}$ and
above $M_{\text{max}}\sim 2\,M_{\odot}$. The precise value of $M_{\text{max}}$
depends on the EoS considered (see Appendix A for more details);

  * •

if we have no determination of the pulsar mass, we draw it from a uniform
distribution (e.g. $M\in\left[0.2\,M_{\odot},2.2M_{\odot}\right]$ for the AP4
model, see Appendix A for more details).

At this point, it is worth obtaining some analytical understanding of the
results that we expect for $\alpha$. In particular, writing the residual
induced by an ULDM candidate of mass $m$ as $\Delta t(m)$, we notice that

$\Delta t(m)\sim\delta
t_{DM}\left(\frac{m}{2}\right)\rightarrow\Psi\sim\Psi_{c}\left(\frac{m}{2}\right),$
(4.7)

which, through Eq. (4.5), yields

$\alpha\sim\frac{mM_{\text{P}}}{\sqrt{\rho}}\Psi_{c}\left(\frac{m}{2}\right).$
(4.8)

By substituting the upper limits for $\Psi_{c}$ found in Fig. 1 of Ref. [40],
Eq. (4.8) gives us an approximate estimate of the upper limits that we expect
to find. Fig 1 shows the comparison between the correlated limit and its
theoretical prediction based on Eq. (4.8), assuming for reference
$\rho=\rho_{\text{DM}}=0.4~{}\text{GeV/cm}^{3}$. As can be seen, our
analytical prediction captures the scaling and the overall shape of the
constraints, but there are deviations caused by the intrinsic difference
between the signals in Eq. (4.4) and Eq. (4.6) (for example, the fact that, in
our scenario, the timing residual depends on the mass of the pulsar through
the sensitivity).

In the following, we present results in terms of $\alpha\sqrt{f_{\text{DM}}}$
to take into account the relative energy density of this ULDM candidate
$f_{\text{DM}}$. This is indeed the quantity constrained by Eq. (4.5), and
allows for rapidly obtaining the relevant bound on $\alpha$ once a value for
the scalar field density $\rho$ is chosen. Fig. 2 displays the upper limits
for the correlated, pulsar correlated and uncorrelated analyses.

We detect the existence of additional signal power above the common red noise
background for masses around $m\sim 10^{-22.7}~{}\text{eV}$ and $m\sim
10^{-21.4}~{}\text{eV}$ across all the three analyses. The first excess is
consistent with what was observed in recent searches [40, 43], while the
second one is thought to be associated with unmodeled perturbations in the
orbital elements of Mercury, whose synodic period matches the detected excess
frequency [41]. While both of them could in principle be interpreted as
evidence of non-minimally coupled ULDM candidates, the fact that they can be
accounted for by different physical models makes us more cautious in drawing
definitive conclusions.

In order our results to be consistent, we also need to ensure that the effect
that we are constraining is not sub-dominant with respect to the TOA induced
by the purely gravitational effect of the ULDM oscillations [33, 40], which we
are neglecting in our analysis. To understand the interplay between our
analysis and the analysis à la Khmelnitsky-Rubakov performed in Ref. [40], it
is sufficient to notice that a non-minimally coupled ULDM candidate of mass
$m$ described by the FJBD action (Eqs. (2.1)–(2.2)) will in general produce
both an $\alpha$-dependent residual and a propagation residual, described
respectively by Eqs. (4.3)–(4.4). While the former has a dependence on the
scalar coupling $\alpha$ and has a characteristic frequency $2\pi f=m$, the
latter only depends on the density of the scalar field and has a
characteristic frequency $2\pi f=2m$ [33, 40]. Therefore, we only need to be
careful selecting a value of $\rho$ which is not already excluded by the
analysis à la Khmelnitsky-Rubakov when extracting a bound on $\alpha$. Fig. 1
shows the effect of this remark on our bounds when we include the constraints
on $\rho$ found in Ref. [40].

Additionally, we notice in passing that the form of the signal in Eq. (4.6) is
similar to Eq. (12) in Ref. [42], where direct couplings between ULDM and
ordinary matter were also studied with PTA data, apart from $\mathcal{O}(1)$
numerical factors. Mapping our signal to Eq. (12) in Ref. [42], we find that
the limits obtained here are in general agreement with their analysis.
Although the two models induce a similar TOA change, it is important to point
out that they are fundamentally different. Indeed, the model presented in Ref.
[42] introduces a direct and particle-dependent coupling of the scalar to
matter, which implies a violation (albeit small) of the weak equivalence
principle. Instead, our model relies on a universal conformal coupling of the
scalar to gravity, and the weak equivalence principle is satisfied (although
the strong equivalence principle is violated).

Figure 1: Upper limits on $\log_{10}\alpha$ at 95% credibility compared to the
analytical estimate described by Eq. (4.8), for the AP4 EoS. The solid line
shows the upper limits on $\log_{10}\alpha$ for the correlated analysis,
assuming that the background DM density is
$\rho_{\text{DM}}=0.4\,\text{GeV/cm}^{3}$, while the brown dashed lines
displays the expected behavior described by Eq. (4.8). The dotted line shows
the degradation of the bounds when choosing
$\rho=\text{min}(\bar{\rho},\rho_{\text{DM}})$, optimistically setting
$\bar{\rho}$ to the upper bounds presented in Ref. [40]. Smaller $\rho$ yield
a stronger degradation of the limits. Figure 2: Upper limits on
$\log_{10}\left(\alpha\sqrt{f_{\text{DM}}}\right)$ at 95% credibility versus
the ULDM mass. We compare results for the correlated, pulsar-correlated and
uncorrelated scenarios in solid, dashed and dotted lines, respectively. The
results are obtained using the AP4 EoS, and the priors on the parameters of
the search are presented in Table 1. Bounds from the Cassini mission and from
the pulsar in a triple stellar system constrain $\alpha^{2}\lesssim 10^{-5}$
and $\alpha^{2}\lesssim 4\times 10^{-6}$, respectively.

###  4.2 DEF theory

To study the DEF theory (quadratic conformal coupling), recall that it is
characterized by $\alpha(\phi)=\beta\phi$, cf. (2.3). Therefore, the angular
momentum sensitivity defined in Eq. (2.11) becomes

$s_{I}=\frac{1}{2\beta\phi}\frac{\mathrm{d}\ln\Omega_{\mathrm{obs}}}{\mathrm{d}\phi}\bigg{|}_{N,J}.$
(4.9)

Here, we cannot neglect the dependence of $s_{I}$ on $\beta$; therefore we
write explicitly $s_{I}=s_{I}(\beta,M)$ or $s_{I}=s_{I}(\beta)$, depending on
the context. The induced timing residual then reads:

$\displaystyle\Delta
t(t)=-\int\frac{\delta\Omega_{\text{obs}}}{\bar{\Omega}}dt=-2\beta
s_{I}(\beta)$
$\displaystyle\int\frac{\rho}{m^{2}M_{\text{P}}^{2}}\hat{\phi}^{2}(\boldsymbol{x})\cos^{2}(mt+\theta(\boldsymbol{x}))$
$\displaystyle=\left.\frac{\Psi}{2m}\beta
s_{I}(\beta)\hat{\phi}^{2}(\boldsymbol{x})\sin(2mt+\theta(\boldsymbol{x}))\right|_{t_{\text{start
}}}^{t_{\text{end }}},$ (4.10)

where we used $\Psi=\rho/(m^{2}M_{\text{P}}^{2})$ and again neglect the
dependence on the retarded time as well as the constant term in the cosine
expansion, as it yields a linear contribution which is absorbed by the pulsar
timing model [101, 102, 103].

As $s_{I}=s_{I}(\beta)$, Eq. (4.10) depends separately on the scalar field
density, parametrized in terms of $\Psi$, and on the DEF scalar coupling
$\beta$. Therefore, unlike the FJBD case, $\rho$ and $\beta$ (or,
equivalently, $\Psi$ and $\beta$) must be two independent parameters of the
search.

This is a crucial observation: in the FJBD case, we constrained a combination
of $\rho$ and $\alpha$, namely $\Psi$, and we rephrased the results into
bounds on $\alpha$ a-posteriori, choosing a reference value of $\rho$ only in
post-processing. Here, instead, we are forced to impose an explicit prior on
$\rho$.
Therefore, we focus on two reference values for the scalar field density;
namely, $\rho=0.5~{}\rho_{\text{DM}}$ and $\rho=\rho_{\text{DM}}$. Moreover,
the deterministic residual in Eq. (4.10) depends on the sign of $\beta$, and
not only on its absolute value. While one might naïvely think that a sign flip
$\beta\rightarrow-\beta$ could be absorbed by a redefinition of the random
phase $\theta(\boldsymbol{x})$, the sensitivity $s_{I}$ does actually depend
on the sign of $\beta$, as shown in Ref. [67] (see e.g. Figs. 5–6 of that
work). Hence, in the following, we present results for both positive and
negative values of $\beta$. As for negative $\beta$, values
$\beta\lesssim-4.3$ would generate non-perturbative strong-field effect
inducing $\mathcal{O}(1)$ variations from GR [52], and are therefore not
considered in the present work (being ruled out by binary pulsars [73]). Fig.
3 displays the upper bounds on $\lvert\beta\rvert$ for the correlated, pulsar-
correlated and uncorrelated limits and for the selected MPA1 EoS, showing that
our analysis implies $\lvert\beta\rvert\lesssim 2.2$ in the range
$10^{-24}~{}\text{eV}\lesssim m\lesssim 10^{-23.5}~{}\text{eV}$ for
$\rho=\rho_{\text{DM}}$. We also plot the constraints obtained for a scalar
field constituting 50% or 30% of DM in Fig. 3. Whenever our analysis is prior
dominated, the upper limits represent the upper end (in absolute value) of our
prior, namely $\lvert\beta\rvert=4.3$. As expected from the form of the signal
in Eq. (4.10), larger values of the scalar field density yield stronger (and
wider) constraints. For positive $\beta$, instead, we focus on $\beta<150$, as
the code provided in Ref. [67] to compute the sensitivities is unstable for
higher values of $\beta$. The results are presented in Fig. 4 for the same
choices of scalar field density. Even in this case, larger values of the
scalar field density translate into more constraining upper bounds, which can
be as low as $\beta\lesssim 20$ for $\rho=\rho_{\text{DM}}$. Again, the upper
bounds are chosen to coincide with the upper end of our prior, namely
$\beta=150$, whenever our analysis is prior dominated.
Finally, let us remind that the analysis in Ref. [40] excludes ultralight
scalar field densities $\rho\gtrsim 0.3~{}\rho_{\text{DM}}$ for masses
$10^{-24}~{}\text{eV}\lesssim m\lesssim 10^{-23.7}~{}\text{eV}$ and
$\rho\gtrsim 0.7~{}\rho_{\text{DM}}$ in the mass range
$10^{-23.7}~{}\text{eV}\lesssim m\lesssim 10^{-23.4}~{}\text{eV}$. This remark
should be taken into account when interpreting the results in Figs. 3-4.

Figure 3: Upper limits on $\lvert\beta\rvert$ ($\beta<0$) at 95% credibility
versus the ULDM mass. The top left panel shows the results for
$\rho=\rho_{\text{DM}}$, the top right panel assumes
$\rho=0.5~{}\rho_{\text{DM}}$ and the bottom panel displays the results for
$\rho=0.3~{}\rho_{\text{DM}}$. We compare results for the correlated, pulsar-
correlated and uncorrelated scenarios in solid, dashed and dotted lines,
respectively. Whenever the bound is prior dominated, the upper limits
represent the upper end (in absolute value) of our prior. Priors on the
parameters relevant for the search are chosen according to the scheme
presented in Table 1.

Figure 4:  Upper limits on $\beta$ ($\beta>0$) at 95% credibility versus the
ULDM mass. The top left panel shows the results for $\rho=\rho_{\text{DM}}$,
the top right panel assumes $\rho=0.5~{}\rho_{\text{DM}}$ and the bottom panel
displays the results for $\rho=0.3~{}\rho_{\text{DM}}$. We compare results for
the correlated, pulsar-correlated and uncorrelated scenarios in solid, dashed
and dotted lines, respectively. Whenever the bound is prior dominated, the
upper limits represent the upper end of our prior. We notice that the pulsar-
correlated case yields valid bounds almost only when the scalar field density
is $\rho=\rho_{\text{DM}}$, while it is completely unconstraining when
$\rho=0.3~{}\rho_{\text{DM}}$. Priors on the parameters relevant for the
search are chosen according to the scheme presented in Table 1.

##  5 Conclusions

Conformally coupled ULDM induces periodic variations in the gravitational mass
and in the radius of pulsars [67], with a timescale given by the mass of the
particle. By conservation of angular momentum, this translates into an
oscillating behavior of the pulsar spin frequency, which is accurately
measured by the PTA collaborations. This effect can be used to set constraints
on the coupling of ULDM to matter, characterized by the conformal factor
$\mathcal{A}(\phi)$ linking the Einstein and Jordan frame metrics. In this
work, we have analyzed the FJBD and the DEF scalar-tensor theories with an
ultralight scalar mass, under the assumption that the scalar field constitutes
(part of) the DM, thus exploiting different functional forms of the conformal
factor $\mathcal{A}(\phi)$.

In the FJBD theory, where $\mathcal{A}(\phi)\sim 1+\alpha\phi$, we find
$\log_{10}\alpha\lesssim-4.5$ across the entire frequency range considered,
vastly overperforming both Cassini bounds [69] and the constraints from the
pulsar in a triple stellar system system [62, 63, 64]. Moreover, for masses
$m\sim 10^{-23}~{}\text{eV}$, our analysis yields the even tighter bound
$\log_{10}\alpha\lesssim-7$. Let us recall, however, that the previously
mentioned bounds have a wider range of applicability than ours, since they
also constrain massless scalar-tensor theories.

In the DEF theory, where $\mathcal{A}\simeq 1+\beta\phi^{2}$, we distinguish
between positive and negative values of $\beta$, which yield a different
expression for the sensitivity $s_{I}$. In the low mass range, we find
$-2\lesssim\beta\lesssim 20$ for a scalar field density
$\rho=\rho_{\text{DM}}$. We also explore how this bound relaxes when the
scalar field density constitutes 50% or 30% of the DM density
$\rho_{\text{DM}}$. Once again, this is competitive with respect to existing
bounds that can be found in the literature [73, 106, 107].

In summary, we have shown that PTA data alone can constrain conformal ULDM
couplings for masses below $\sim 10^{-21}$ eV at levels not yet explored by
other observations. These models include scenarios where the ULDM constitutes
all the dark matter in the Universe, and scenarios where the ULDM is a
fraction of the total dark matter content. All the future improvements of PTA
searches (e.g. those related to SKA [108]) will impact the searches we
performed in this work. Furthermore, it was recently emphasized in [109] that
the effects coming from the interference of the different modes comprising the
ULDM field (recall (2.6) and the velocity dispersion $\sigma_{\phi}\sim
10^{-3}$) will generate a signal at frequencies below $m\sigma_{\phi}$ that
may allow our current analysis to access ULDM models of higher masses.

## Acknowledgments

C. Smarra wishes to thank Cecilia Sgalletta for useful discussions about data
processing. E. Barausse, A. Kuntz and C. Smarra acknowledge support from the
European Union’s H2020 ERC Consolidator Grant “GRavity from Astrophysical to
Microscopic Scales” (Grant No. GRAMS-815673), the PRIN 2022 grant “GUVIRP -
Gravity tests in the UltraViolet and InfraRed with Pulsar timing”, and the EU
Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-
Curie Grant Agreement No. 101007855. The research leading to these results has
received funding from the Spanish Ministry of Science and Innovation
(PID2020-115845GB-I00/AEI/10.13039/501100011033). IFAE is partially funded by
the CERCA program of the Generalitat de Catalunya. D. Blas acknowledges the
support from the Departament de Recerca i Universitats de la Generalitat de
Catalunya al Grup de Recerca i Universitats from Generalitat de Catalunya to
the Grup de Recerca 00649 (Codi: 2021 SGR 00649). D. López Nacir acknowledges
support from University of Buenos Aires and CONICET. LS was supported by the
National SKA Program of China (2020SKA0120300), the Beijing Natural Science
Foundation (1242018), and the Max Planck Partner Group Program funded by the
Max Planck Society. J.Antoniadis acknowledges support from the European
Commission (ARGOS-CDS; Grant Agreement number: 101094354) The Nançay radio
Observatory is operated by the Paris Observatory, associated with the French
Centre National de la Recherche Scientifique (CNRS), and partially supported
by the Region Centre in France. I. Cognard, L. Guillemot and G. Theureau
acknowledge financial support from “Programme National de Cosmologie and
Galaxies” (PNCG), and “Programme National Hautes Energies” (PNHE) funded by
CNRS/INSU-IN2P3-INP, CEA and CNES, France. I. Cognard, L. Guillemot and G.
Theureau acknowledge financial support from Agence Nationale de la Recherche
(ANR-18-CE31-0015), France. Pulsar research at Jodrell Bank Centre for
Astrophysics is supported by an STFC Consolidated Grant (ST/T000414/1;
ST/X001229/1). This work is also supported as part of the “LEGACY” MPG-CAS
collaboration on 1052 low-frequency gravitational wave astronomy. D. Perrodin
acknowledges support from the PRIN 2022 grant “GUVIRP - Gravity tests in the
UltraViolet and InfraRed with Pulsar timing” and the INAF 2023 Large Grant
“Gravitational Wave Detection using Pulsar Timing Arrays”.

##  Appendix A Parameters of the search

Table 1 summarizes the parameters used in the search along with their priors.
We will add a label F if the parameter is used only in the FJBD analysis,
while a label D will signal parameters used only in the DEF analysis. As
stated in the main text, if a pulsar mass is measured from other experiments
to be $M\pm\delta M$, we sample the mass parameter from a normal prior
distribution centered on $M$, with uncertainty $\delta M$ and truncated below
$M_{\text{min}}=0.2\,M_{\odot}$ and above
$M_{\text{max}}=(2.2\,M_{\odot},2.4\,M_{\odot},2.2\,M_{\odot})$ for the (AP4
[110], MPA1 [111], SLy [112]) EoS, respectively. We will label this as
TruncNorm($\mu,\sigma$), where $\mu=M$ and $\sigma=\delta M$. Otherwise, we
will just assume a uniform prior between $M_{\text{min}}$ and $M_{\text{max}}$
for the three EoS choices. We also plot the constraints on the FJBD scalar
coupling $\alpha$ for the different EoS of the pulsar interior presented in
Ref. [67].

Table 1: Parameters employed for the search along with their respective priors. In the correlated limit, the parameter $\hat{\phi}^{2}_{\text{P}}$ is accounted for by a redefinition of $\Psi$, while in the pulsar-correlated regime $\hat{\phi}^{2}_{\text{P}}=\hat{\phi}^{2}$ is a common free parameter. Only the pulsars whose masses have been measured from other experiments are presented, along with the relevant reference. For the other pulsars, we choose uniform priors with support $[M_{\text{min}},M_{\text{max}}]$ (see main text for more details). We display the priors on the DEF scalar coupling $\beta$ both for $\beta<0$ and $\beta>0$. $\mathcal{U}$ stands for the uniform distribution and $\mathcal{N}$ stands for the TruncNorm distribution (see main text). The white noise parameters EFAC (TOA error Excess FACtor) and EQUAD (TOA Error excess in QUADrature) are introduced for every receiver-backend system in every pulsar. Parameter | Description | Prior | Occurrence
---|---|---|---
White Noise
$\left(\sigma=E_{\text{f}}^{2}\sigma^{2}_{\text{TOA}}+E_{\text{q}}^{2}\right)$
$E_{\text{f}}$ | EFAC | $\mathcal{U}(0,10)$ | 1 per pulsar
$E_{\text{q}}$ | EQUAD | Log10-$\mathcal{U}(-10,-5)$ | 1 per pulsar
Red Noise (RN)
$A_{\text{red }}$ | RN power-law amplitude | Log10-$\mathcal{U}(-20,-11)$ | 1 per pulsar
$\gamma_{\text{red }}$ | RN power-law spectral index | $\mathcal{U}(0,10)$ | 1 per pulsar
ULDM (F)
$\Psi$ | ULDM signal amplitude |  $\text{Log}_{10}$-$\mathcal{U}(-20,-12)$ | 1 per PTA
$m~{}[\mathrm{eV}]$ | ULDM mass | Log10-$\mathcal{U}(-24,-21)$ | 1 per PTA
$\hat{\phi}^{2}$ | Pulsar factor | $e^{-x}$ | 1 per pulsar
$\theta$ | Pulsar signal phase | $\mathcal{U}(0,2\pi)$ | 1 per pulsar
ULDM (D)
$f_{\text{DM}}$ | ULDM fraction | $\mathcal{U}(0.01,0.30)$ | 1 per PTA
$\beta$ | DEF scalar coupling |  |  $\mathcal{U}(-4.3,0)$
---
or
$\mathcal{U}(0,150)$
1 per PTA
$m~{}[\mathrm{eV}]$ | ULDM mass | Log10-$\mathcal{U}(-24,-21)$ | 1 per PTA
$\hat{\phi}^{2}$ | Pulsar factor | $e^{-x}$ | 1 per pulsar
$\theta$ | Pulsar signal phase | $\mathcal{U}(0,2\pi)$ | 1 per pulsar
Common spatially Uncorrelated Red Noise (CURN)
$A_{\mathrm{GWB}}$ | CURN strain amplitude | Log10-$\mathcal{U}(-18,-11)$ | 1 per PTA
$\gamma_{\mathrm{GWB}}$ | CURN spectral index | $\mathcal{U}(0,7)$ | 1 per PTA
Pulsar Masses
$\text{M}_{\text{J0030}}$ [113] | PSR J0030+0451 mass | $\mathcal{N}(1.44,0.15)$ | 1 per PTA
$\text{M}_{\text{J1012}}$ [114] | PSR J1012+5307 mass | $\mathcal{N}(1.72,0.16)$ | 1 per PTA
$\text{M}_{\text{J1713}}$ [115] | PSR J1713+0747 mass | $\mathcal{N}(1.3,0.2)$ | 1 per PTA
$\text{M}_{\text{J1738}}$ [116] | PSR J1738+0333 mass | $\mathcal{N}(1.47,0.07)$ | 1 per PTA
$\text{M}_{\text{J1909}}$ [117] | PSR J1909-3744 mass | $\mathcal{N}(1.438,0.024)$ | 1 per PTA

Figure 5: Upper limits on $\log_{10}\alpha$ at 95% credibility versus the ULDM
mass, using the prescription for $\rho_{\text{DM}}$ detailed in the main text.
We compare results for the correlated, pulsar-correlated and uncorrelated
scenarios in solid, dashed and dotted lines, respectively. The results are
obtained using the MPA1 (left panel) and the SLy (right panel) EoS, while the
priors on the parameters of the search are chosen according to Table 1.

## References

  * [1] R.A. Flores and J.R. Primack, _Observational and Theoretical Constraints on Singular Dark Matter Halos_ , _ApJ_ 427 (1994) L1 [astro-ph/9402004].
  * [2] B. Moore, _Evidence against dissipation-less dark matter from observations of galaxy haloes_ , _Nature_ 370 (1994) 629.
  * [3] Karukes, E. V., Salucci, P. and Gentile, G., _The dark matter distribution in the spiral ngc 3198 out to 0.22 rvir_ , _A &A_ 578 (2015) A13.
  * [4] A. Klypin, A.V. Kravtsov, O. Valenzuela and F. Prada, _Where are the missing galactic satellites?_ , _The Astrophysical Journal_ 522 (1999) 82.
  * [5] B. Moore, S. Ghigna, F. Governato, G. Lake, T. Quinn, J. Stadel et al., _Dark Matter Substructure within Galactic Halos_ , _ApJ_ 524 (1999) L19 [astro-ph/9907411].
  * [6] M. Boylan-Kolchin, J.S. Bullock and M. Kaplinghat, _Too big to fail? The puzzling darkness of massive Milky Way subhaloes_ , _MNRAS_ 415 (2011) L40 [1103.0007].
  * [7] J.F. Navarro, V.R. Eke and C.S. Frenk, _The cores of dwarf galaxy haloes_ , _Monthly Notices of the Royal Astronomical Society_ 283 (1996) L72.
  * [8] F. Governato, A. Zolotov, A. Pontzen, C. Christensen, S.H. Oh, A.M. Brooks et al., _Cuspy no more: how outflows affect the central dark matter and baryon distribution in $\Lambda$ cold dark matter galaxies_, _Monthly Notices of the Royal Astronomical Society_ 422 (2012) 1231 [https://academic.oup.com/mnras/article-pdf/422/2/1231/3466602/mnras0422-1231.pdf].
  * [9] A.M. Brooks, M. Kuhlen, A. Zolotov and D. Hooper, _A baryonic solution to the missing satellites problem_ , _The Astrophysical Journal_ 765 (2013) 22.
  * [10] T.K. Chan, D. Kereš, J. Oñorbe, P.F. Hopkins, A.L. Muratov, C.-A. Faucher-Giguère et al., _The impact of baryonic physics on the structure of dark matter haloes: the view from the FIRE cosmological simulations_ , _Monthly Notices of the Royal Astronomical Society_ 454 (2015) 2981 [https://academic.oup.com/mnras/article-pdf/454/3/2981/4038253/stv2165.pdf].
  * [11] J. Oñorbe, M. Boylan-Kolchin, J.S. Bullock, P.F. Hopkins, D. Kereš, C.-A. Faucher-Giguère et al., _Forged in fire: cusps, cores and baryons in low-mass dwarf galaxies_ , _Monthly Notices of the Royal Astronomical Society_ 454 (2015) 2092 [https://academic.oup.com/mnras/article-pdf/454/2/2092/9503920/stv2072.pdf].
  * [12] J.I. Read, O. Agertz and M.L.M. Collins, _Dark matter cores all the way down_ , _Monthly Notices of the Royal Astronomical Society_ 459 (2016) 2573 [https://academic.oup.com/mnras/article-pdf/459/3/2573/8105757/stw713.pdf].
  * [13] W. Hu, R. Barkana and A. Gruzinov, _Cold and fuzzy dark matter_ , _Phys. Rev. Lett._ 85 (2000) 1158 [astro-ph/0003365].
  * [14] L. Hui, J.P. Ostriker, S. Tremaine and E. Witten, _Ultralight scalars as cosmological dark matter_ , _Physical Review D_ 95 (2017) .
  * [15] M.B. Green, J.H. Schwarz and E. Witten, _SUPERSTRING THEORY. VOL. 1: INTRODUCTION_ , Cambridge Monographs on Mathematical Physics (7, 1988).
  * [16] P. Svrcek and E. Witten, _Axions in string theory_ , _Journal of High Energy Physics_ 2006 (2006) 051.
  * [17] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, _String axiverse_ , _Phys. Rev. D_ 81 (2010) 123530.
  * [18] R. Hlozek, D. Grin, D.J.E. Marsh and P.G. Ferreira, _A search for ultralight axions using precision cosmological data_ , _Phys. Rev. D_ 91 (2015) 103512.
  * [19] V. Iršič, M. Viel, M.G. Haehnelt, J.S. Bolton and G.D. Becker, _First constraints on fuzzy dark matter from lyman-alpha forest data and hydrodynamical simulations_ , _Phys. Rev. Lett._ 119 (2017) 031302.
  * [20] E. Armengaud, N. Palanque-Delabrouille, C. Yèche, D.J.E. Marsh and J. Baur, _Constraining the mass of light bosonic dark matter using SDSS Lyman-alpha forest_ , _Monthly Notices of the Royal Astronomical Society_ 471 (2017) 4606 [https://academic.oup.com/mnras/article-pdf/471/4/4606/19635244/stx1870.pdf].
  * [21] T. Kobayashi, R. Murgia, A.D. Simone, V. Iršič and M. Viel, _Lyman-alpha constraints on ultralight scalar dark matter: Implications for the early and late universe_ , _Physical Review D_ 96 (2017) .
  * [22] M. Nori, R. Murgia, V. Iršič, M. Baldi and M. Viel, _Lyman-alpha forest and non-linear structure characterization in Fuzzy Dark Matter cosmologies_ , _Monthly Notices of the Royal Astronomical Society_ 482 (2018) 3227 [https://academic.oup.com/mnras/article-pdf/482/3/3227/26653692/sty2888.pdf].
  * [23] K.K. Rogers and H.V. Peiris, _Strong bound on canonical ultralight axion dark matter from the lyman-alpha forest_ , _Phys. Rev. Lett._ 126 (2021) 071302.
  * [24] H.-Y. Schive, M.-H. Liao, T.-P. Woo, S.-K. Wong, T. Chiueh, T. Broadhurst et al., _Understanding the core-halo relation of quantum wave dark matter from 3d simulations_ , _Phys. Rev. Lett._ 113 (2014) 261302.
  * [25] J. Zhang, H. Liu and M.-C. Chu, _Cosmological simulation for fuzzy dark matter model_ , _Frontiers in Astronomy and Space Sciences_ 5 (2019) .
  * [26] N. Bar, D. Blas, K. Blum and S. Sibiryakov, _Galactic rotation curves versus ultralight dark matter: Implications of the soliton-host halo relation_ , _Phys. Rev. D_ 98 (2018) 083027 [1805.00122].
  * [27] K. Hayashi, E.G.M. Ferreira and H.Y.J. Chan, _Narrowing the mass range of fuzzy dark matter with ultrafaint dwarfs_ , _The Astrophysical Journal Letters_ 912 (2021) L3.
  * [28] N. Dalal and A. Kravtsov, _Not so fuzzy: excluding fdm with sizes and stellar kinematics of ultra-faint dwarf galaxies_ , 2022.  https://doi.org/10.48550/arXiv.2203.05750.
  * [29] D.J.E. Marsh and J.C. Niemeyer, _Strong Constraints on Fuzzy Dark Matter from Ultrafaint Dwarf Galaxy Eridanus II_ , _Phys. Rev. Lett._ 123 (2019) 051103 [1810.08543].
  * [30] E.G.M. Ferreira, _Ultra-light dark matter_ , _Astron. Astrophys. Rev._ 29 (2021) 7 [2005.03254].
  * [31] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, _String Axiverse_ , _Phys. Rev. D_ 81 (2010) 123530 [0905.4720].
  * [32] L. Hamaide, H. Müller and D.J.E. Marsh, _Searching for dilaton fields in the Lyman- $\alpha$ forest_, _Phys. Rev. D_ 106 (2022) 123509 [2210.03705].
  * [33] A. Khmelnitsky and V. Rubakov, _Pulsar timing signal from ultralight scalar dark matter_ , _Journal of Cosmology and Astroparticle Physics_ 2014 (2014) 019.
  * [34] R.S. Foster and D.C. Backer, _Constructing a Pulsar Timing Array_ , _ApJ_ 361 (1990) 300.
  * [35] R.N. Manchester, G. Hobbs, M. Bailes, W.A. Coles, W. van Straten, M.J. Keith et al., _The Parkes Pulsar Timing Array Project_ , _PASA_ 30 (2013) e017 [1210.6130].
  * [36] M.A. McLaughlin, _The North American Nanohertz Observatory for Gravitational Waves_ , _Classical and Quantum Gravity_ 30 (2013) 224008 [1310.0758].
  * [37] M. Kramer and D.J. Champion, _The European Pulsar Timing Array and the Large European Array for Pulsars_ , _Classical and Quantum Gravity_ 30 (2013) 224009.
  * [38] B.C. Joshi, P. Arumugasamy, M. Bagchi, D. Bandyopadhyay, A. Basu, N. Dhanda Batra et al., _Precision pulsar timing with the ORT and the GMRT and its applications in pulsar astrophysics_ , _Journal of Astrophysics and Astronomy_ 39 (2018) 51.
  * [39] M. Bailes, A. Jameson, F. Abbate, E.D. Barr, N.D.R. Bhat, L. Bondonneau et al., _The MeerKAT telescope as a pulsar facility: System verification and early science results from MeerTime_ , _PASA_ 37 (2020) e028 [2005.14366].
  * [40] EPTA Collaboration, C. Smarra, B. Goncharov, E. Barausse, J. Antoniadis, S. Babak et al., _Second data release from the european pulsar timing array: Challenging the ultralight dark matter paradigm_ , _Phys. Rev. Lett._ 131 (2023) 171001.
  * [41] N.K. Porayko, X. Zhu, Y. Levin, L. Hui, G. Hobbs, A. Grudskaya et al., _Parkes pulsar timing array constraints on ultralight scalar-field dark matter_ , _Physical Review D_ 98 (2018) .
  * [42] D.E. Kaplan, A. Mitridate and T. Trickle, _Constraining fundamental constant variations from ultralight dark matter with pulsar timing arrays_ , _Phys. Rev. D_ 106 (2022) 035032.
  * [43] A. Afzal, G. Agazie, A. Anumarlapudi, A.M. Archibald, Z. Arzoumanian, P.T. Baker et al., _The nanograv 15 yr data set: Search for signals from new physics_ , _The Astrophysical Journal Letters_ 951 (2023) L11.
  * [44] D. Blas, D.L. Nacir and S. Sibiryakov, _Ultralight dark matter resonates with binary pulsars_ , _Physical Review Letters_ 118 (2017) .
  * [45] D. Blas, D.L. Nacir and S. Sibiryakov, _Secular effects of ultralight dark matter on binary pulsars_ , _Physical Review D_ 101 (2020) .
  * [46] P. Kůs, D. López Nacir and F.R. Urban, _Bayesian sensitivity of binary pulsars to ultra-light dark matter_ , 2402.04099.
  * [47] M. Fierz, _On the physical interpretation of P.Jordan’s extended theory of gravitation_ , _Helv. Phys. Acta_ 29 (1956) 128.
  * [48] P. Jordan, _The present state of Dirac’s cosmological hypothesis_ , _Z. Phys._ 157 (1959) 112.
  * [49] C. Brans and R.H. Dicke, _Mach’s principle and a relativistic theory of gravitation_ , _Phys. Rev._ 124 (1961) 925.
  * [50] R.H. Dicke, _Mach’s principle and invariance under transformation of units_ , _Phys. Rev._ 125 (1962) 2163.
  * [51] T. Damour and G. Esposito-Farese, _Tensor multiscalar theories of gravitation_ , _Class. Quant. Grav._ 9 (1992) 2093.
  * [52] T. Damour and G. Esposito-Farèse, _Nonperturbative strong-field effects in tensor-scalar theories of gravitation_ , _Phys. Rev. Lett._ 70 (1993) 2220.
  * [53] J. Alsing, E. Berti, C.M. Will and H. Zaglauer, _Gravitational radiation from compact binary systems in the massive Brans-Dicke theory of gravity_ , _Phys. Rev. D_ 85 (2012) 064041 [1112.4903].
  * [54] K. Nordtvedt, _Equivalence principle for massive bodies. ii. theory_ , _Phys. Rev._ 169 (1968) 1017.
  * [55] D.M. Eardley, _Observable effects of a scalar gravitational field in a binary pulsar_ , _Astrophysical Journal_ 196 (1975) L59.
  * [56] C.M. Will and H.W. Zaglauer, _Gravitational Radiation, Close Binary Systems, and the Brans-Dicke Theory of Gravity_ , _ApJ_ 346 (1989) 366.
  * [57] C.M. Will, _Gravitational radiation from binary systems in alternative metric theories of gravity: dipole radiation and the binary pulsar._ , _ApJ_ 214 (1977) 826.
  * [58] C.M. Will, _Theory and experiment in gravitational physics_ (1993).
  * [59] T. Damour and J.H. Taylor, _Strong field tests of relativistic gravity and binary pulsars_ , _Phys. Rev. D_ 45 (1992) 1840.
  * [60] J.M. Weisberg and J.H. Taylor, _Relativistic binary pulsar b1913+16: Thirty years of observations and analysis_ , 2004.
  * [61] M. Kramer, I.H. Stairs, R. Manchester, M. McLaughlin, A. Lyne et al., _Tests of general relativity from timing the double pulsar_ , _Science_ 314 (2006) 97 [astro-ph/0609417].
  * [62] S.M. Ransom, I.H. Stairs, A.M. Archibald, J.W.T. Hessels, D.L. Kaplan, M.H. van Kerkwijk et al., _A millisecond pulsar in a stellar triple system_ , _Nature_ 505 (2014) 520.
  * [63] A.M. Archibald, N.V. Gusinskaia, J.W.T. Hessels, A.T. Deller, D.L. Kaplan, D.R. Lorimer et al., _Universality of free fall from the orbital motion of a pulsar in a stellar triple system_ , _Nature_ 559 (2018) 73.
  * [64] Voisin, G., Cognard, I., Freire, P. C. C., Wex, N., Guillemot, L., Desvignes, G. et al., _An improved test of the strong equivalence principle with the pulsar in a triple star system_ , _A & A_ 638 (2020) A24.
  * [65] M. Kramer, I.H. Stairs, R.N. Manchester, N. Wex, A.T. Deller, W.A. Coles et al., _Strong-field gravity tests with the double pulsar_ , _Phys. Rev. X_ 11 (2021) 041050.
  * [66] J. Zhao, P.C.C. Freire, M. Kramer, L. Shao and N. Wex, _Closing a spontaneous-scalarization window with binary pulsars_ , _Class. Quant. Grav._ 39 (2022) 11LT01 [2201.03771].
  * [67] A. Kuntz and E. Barausse, _Angular momentum sensitivities in scalar-tensor theories_ , 2024.
  * [68] T. Damour and G. Esposito-Farèse, _Tensor-scalar gravity and binary-pulsar experiments_ , _Phys. Rev. D_ 54 (1996) 1474.
  * [69] B. Bertotti, L. Iess and P. Tortora, _A test of general relativity using radio links with the cassini spacecraft_ , _Nature_ 425 (2003) 374.
  * [70] T. Damour and A.M. Polyakov, _The String dilaton and a least coupling principle_ , _Nucl. Phys. B_ 423 (1994) 532 [hep-th/9401069].
  * [71] J.E. Kim, _Light Pseudoscalars, Particle Physics and Cosmology_ , _Phys. Rept._ 150 (1987) 1.
  * [72] C.M. Will, _Theory and Experiment in Gravitational Physics_ , Cambridge University Press (1993), 10.1017/CBO9780511564246.
  * [73] L. Shao, N. Sennett, A. Buonanno, M. Kramer and N. Wex, _Constraining nonperturbative strong-field effects in scalar-tensor gravity by combining pulsar timing and laser-interferometer gravitational-wave detectors_ , _Phys. Rev. X_ 7 (2017) 041025 [1704.07561].
  * [74] A. Castillo, J. Martin-Camalich, J. Terol-Calvo, D. Blas, A. Caputo, R.T.G. Santos et al., _Searching for dark-matter waves with PPTA and QUIJOTE pulsar polarimetry_ , _JCAP_ 06 (2022) 014 [2201.03422].
  * [75] F. Nesti and P. Salucci, _The Dark Matter halo of the Milky Way, AD 2013_ , _JCAP_ 07 (2013) 016 [1304.5127].
  * [76] EPTA Collaboration, Antoniadis, J., Babak, S., Bak Nielsen, A.-S., Bassa, C. G., Berthereau, A. et al., _The second data release from the european pulsar timing array - i. the dataset and timing analysis_ , _A & A_ 678 (2023) A48.
  * [77] G. Desvignes, R.N. Caballero, L. Lentati, J.P.W. Verbiest, D.J. Champion, B.W. Stappers et al., _High-precision timing of 42 millisecond pulsars with the European Pulsar Timing Array_ , _MNRAS_ 458 (2016) 3341 [1602.08511].
  * [78] R.T. Edwards, G.B. Hobbs and R.N. Manchester, _TEMPO2, a new pulsar timing package - II. The timing model and precision estimates_ , _MNRAS_ 372 (2006) 1549 [astro-ph/0607664].
  * [79] B. Goncharov, D.J. Reardon, R.M. Shannon, X.-J. Zhu, E. Thrane, M. Bailes et al., _Identifying and mitigating noise sources in precision pulsar timing data sets_ , _MNRAS_ 502 (2021) 478 [2010.06109].
  * [80] EPTA Collaboration and InPTA Collaboration, Antoniadis, J., Arumugam, P., Arumugam, S., Babak, S., Bagchi, M. et al., _The second data release from the european pulsar timing array - ii. customised pulsar noise models for spatially correlated gravitational waves_ , _A & A_ 678 (2023) A49.
  * [81] NANOGrav collaboration, _The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background_ , _Astrophys. J. Lett._ 951 (2023) L8 [2306.16213].
  * [82] EPTA, InPTA: collaboration, _The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals_ , _Astron. Astrophys._ 678 (2023) A50 [2306.16214].
  * [83] D.J. Reardon et al., _Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 951 (2023) L6 [2306.16215].
  * [84] H. Xu, S. Chen, Y. Guo, J. Jiang, B. Wang, J. Xu et al., _Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I_ , _Research in Astronomy and Astrophysics_ 23 (2023) 075024 [2306.16216].
  * [85] R.W. Hellings and G.S. Downs, _Upper limits on the isotropic gravitational radiation background from pulsar timing analysis._ , _ApJ_ 265 (1983) L39.
  * [86] EPTA Collaboration and InPTA Collaboration, Antoniadis, J., Arumugam, P., Arumugam, S., Babak, S., Bagchi, M. et al., _The second data release from the european pulsar timing array - iii. search for gravitational wave signals_ , _A & A_ 678 (2023) A50.
  * [87] J. Antoniadis, P. Arumugam, S. Arumugam, P. Auclair, S. Babak, M. Bagchi et al., _The second data release from the european pulsar timing array: V. implications for massive black holes, dark matter and the early universe_ , 2023\.
  * [88] Z. Arzoumanian, P.T. Baker, H. Blumer, B. Bécsy, A. Brazier, P.R. Brook et al., _The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background_ , _ApJ_ 905 (2020) L34 [2009.04496].
  * [89] B. Goncharov, R.M. Shannon, D.J. Reardon, G. Hobbs, A. Zic, M. Bailes et al., _On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _ApJ_ 917 (2021) L19 [2107.12112].
  * [90] S. Chen, R.N. Caballero, Y.J. Guo, A. Chalumeau, K. Liu, G. Shaifullah et al., _Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: inferences in the stochastic gravitational-wave background search_ , _MNRAS_ 508 (2021) 4970 [2110.13184].
  * [91] R. van Haasteren, Y. Levin, P. McDonald and T. Lu, _On measuring the gravitational-wave background using Pulsar Timing Arrays_ , _Monthly Notices of the Royal Astronomical Society_ 395 (2009) 1005 [0809.0791].
  * [92] L. Lentati, P. Alexander, M.P. Hobson, F. Feroz, R. van Haasteren, K.J. Lee et al., _TEMPONEST: a Bayesian approach to pulsar timing analysis_ , _Monthly Notices of the Royal Astronomical Society_ 437 (2014) 3004 [1310.2120].
  * [93] S.R. Taylor, L. Lentati, S. Babak, P. Brem, J.R. Gair, A. Sesana et al., _All correlations must die: Assessing the significance of a stochastic gravitational-wave background in pulsar timing arrays_ , _Phys. Rev. D_ 95 (2017) 042002 [1606.09180].
  * [94] J.A. Ellis, M. Vallisneri, S.R. Taylor and P.T. Baker, “Enterprise: Enhanced numerical toolbox enabling a robust pulsar inference suite.” Zenodo, Sept., 2020\.  10.5281/zenodo.4059815.
  * [95] S.R. Taylor, P.T. Baker, J.S. Hazboun, J. Simon and S.J. Vigeland, _enterprise_extensions_ , 2021.
  * [96] J. Ellis and R. van Haasteren, _jellis18/ptmcmcsampler: Official release_ , Oct., 2017.  10.5281/zenodo.1037579.
  * [97] A. Mitridate, _Ptarcade_ , .
  * [98] A. Mitridate, D. Wright, R. von Eckardstein, T. Schröder, J. Nay, K. Olum et al., _PTArcade_ , 2306.16377.
  * [99] J.P.W. Verbiest, J.M. Weisberg, A.A. Chael, K.J. Lee and D.R. Lorimer, _On pulsar distance measurements and their uncertainties_ , _The Astrophysical Journal_ 755 (2012) 39.
  * [100] J. chan Hwang, D. Jeong, H. Noh and C. Smarra, _Pulsar timing array signature from oscillating metric perturbations due to ultra-light axion_ , _Journal of Cosmology and Astroparticle Physics_ 2024 (2024) 014.
  * [101] J.S. Hazboun, J.D. Romano and T.L. Smith, _Realistic sensitivity curves for pulsar timing arrays_ , _Phys. Rev. D_ 100 (2019) 104028 [1907.04341].
  * [102] R. Blandford, R. Narayan and R.W. Romani, _Arrival-time analysis for a millisecond pulsar._ , _Journal of Astrophysics and Astronomy_ 5 (1984) 369.
  * [103] H. Ramani, T. Trickle and K.M. Zurek, _Observability of dark matter substructure with pulsar timing correlations_ , _Journal of Cosmology and Astroparticle Physics_ 2020 (2020) 033.
  * [104] C. Unal, F.R. Urban and E.D. Kovetz, _Probing ultralight scalar, vector and tensor dark matter with pulsar timing arrays_ , 2022.  https://doi.org/10.48550/arXiv.2209.02741.
  * [105] M. Abney, R.I. Epstein and A.V. Olinto, _Observational constraints on the internal structure and dynamics of the vela pulsar_ , _The Astrophysical Journal_ 466 (1996) L91.
  * [106] R.F.P. Mendes, _Possibility of setting a new constraint to scalar-tensor theories_ , _Phys. Rev. D_ 91 (2015) 064024 [1412.6789].
  * [107] D. Anderson and N. Yunes, _Solar System constraints on massless scalar-tensor gravity with positive coupling constant upon cosmological evolution of the scalar field_ , _Phys. Rev. D_ 96 (2017) 064037 [1705.06351].
  * [108] G. Janssen et al., _Gravitational wave astronomy with the SKA_ , _PoS_ AASKA14 (2015) 037 [1501.00127].
  * [109] H. Kim, _Gravitational interaction of ultralight dark matter with interferometers_ , _JCAP_ 12 (2023) 018 [2306.13348].
  * [110] A. Akmal, V.R. Pandharipande and D.G. Ravenhall, _Equation of state of nucleon matter and neutron star structure_ , _Phys. Rev. C_ 58 (1998) 1804.
  * [111] H. Müther, M. Prakash and T. Ainsworth, _The nuclear symmetry energy in relativistic brueckner-hartree-fock calculations_ , _Physics Letters B_ 199 (1987) 469.
  * [112] F. Douchin and P. Haensel, _A unified equation of state of dense matter and neutron star structure_ , _Astronomy & Astrophysics_ 380 (2001) 151–167.
  * [113] M.C. Miller, F.K. Lamb, A.J. Dittmann, S. Bogdanov, Z. Arzoumanian, K.C. Gendreau et al., _Psr j0030+0451 mass and radius from nicer data and implications for the properties of neutron star matter_ , _The Astrophysical Journal Letters_ 887 (2019) L24.
  * [114] D. Mata Sánchez, A.G. Istrate, M.H. van Kerkwijk, R.P. Breton and D.L. Kaplan, _PSR J1012+5307: a millisecond pulsar with an extremely low-mass white dwarf companion_ , _Monthly Notices of the Royal Astronomical Society_ 494 (2020) 4031 [https://academic.oup.com/mnras/article-pdf/494/3/4031/33149104/staa983.pdf].
  * [115] Chen, W.-C. and Panei, J. A., _The progenitor of binary millisecond radio pulsar psr j1713+0747_ , _A & A_ 527 (2011) A128.
  * [116] J. Antoniadis, M.H. van Kerkwijk, D. Koester, P.C.C. Freire, N. Wex, T.M. Tauris et al., _The relativistic pulsar–white dwarf binary PSR J1738+0333 – I. Mass determination and evolutionary history_ , _Monthly Notices of the Royal Astronomical Society_ 423 (2012) 3316 [https://academic.oup.com/mnras/article-pdf/423/4/3316/4895633/mnras0423-3316.pdf].
  * [117] B.A. Jacoby, A. Hotan, M. Bailes, S. Ord and S.R. Kulkarni, _The mass of a millisecond pulsar_ , _The Astrophysical Journal_ 629 (2005) L113.