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

Clemente Smarra, Adrien Kuntz, Enrico Barausse, Boris Goncharov, Diana López Nacir, Diego Blas, Lijing Shao, J Antoniadis, D.J Champion, Ismaël Cognard, et al.

► To cite this version:

Clemente Smarra, Adrien Kuntz, Enrico Barausse, Boris Goncharov, Diana López Nacir, et al. Constraints on conformal ultralight dark matter couplings from the European Pulsar Timing Array. Physical Review D, 2024, 110 (4), pp.043033. 10.1103/PhysRevD.110.043033. hal-04584540

HAL Id: hal-04584540

https://hal.science/hal-04584540v1

Submitted on 15 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraints on conformal ultralight dark matter couplings from the European Pulsar Timing Array

Clemente Smarra,¹,²,* Adrien Kuntz,¹,²,† Enrico Barausse,¹,² Boris Goncharov

Diana López Nacir,⁷,⁸ Diego Blas,⁹,¹⁰ Lijing Shao,¹¹,¹² J. Antoniadis,¹³,¹⁴

D. J. Champion,¹⁴ I. Cognard,¹⁵,¹⁶ L. Guillemot,¹⁵,¹⁶ H. Hu,¹⁴ M. Keith,¹⁷ M.

Kramer,¹⁴,¹⁷ K. Liu,¹⁴,¹⁸ D. Perrodin,¹⁹ S. A. Sanidas,¹⁷ and G. Theureau,¹⁵,²⁰,²¹

¹SISSA, Via Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste

²IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy

³Gran Sasso Science Institute (GSSI), I-67100 L'Aquila, Italy

⁴INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy

⁵Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Hannover, Germany

⁶Leibniz Universität Hannover, D-30167 Hannover, Germany

⁷Universidad de Buenos Aires, Facultad de Ciencias Exactas y Naturales, Departamento de Física. Buenos Aires, Argentina.

⁸CONICET - Universidad de Buenos Aires, Instituto de Física de Buenos Aires (IFIBA). Buenos Aires, Argentina

⁹Institut de Fisica d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain

¹⁰Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluís Companys 23, 08010 Barcelona, Spain

¹¹Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China

¹²National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China

¹³FORTH Institute of Astrophysics, N. Plastira 100, 70013, Heraklion, Greece

¹⁴Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

¹⁵Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans/CNRS, 45071 Orléans Cedex 02, France

¹⁶Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France

¹⁷Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Manchester M13 9PL

¹⁸Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China

¹⁹Affiliation: INAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, 09047 Selargius (CA), Italy

²⁰Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France

²¹Laboratoire Univers et Théories LUTh, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190 Meudon, France

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.

## I. 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 (~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–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],

<footnote>* csmarra@sissa.it</footnote>

<footnote>† akuntz@sissa.it</footnote>
and simulations based on Λ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–12], another possibility is to assume that DM is an ultralight scalar field (with mass $m \sim 10^{-22}$ eV) with negligible self-interactions [13, 14]. In this scenario, the de Broglie wavelength of the scalar field in galaxies can reach $\sim$ 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–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}$ eV [18]. Lyman-$\alpha$ observations strongly suggests a lower limit of $m \gtrsim 10^{-21}$ eV for ultra-light candidates accounting for more than $\sim$ 30% of the DM [19–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}$ 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}$ eV [14, 29]. 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–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/cm$^3$ in the mass range $10^{-24.0}$ eV $\lesssim m \lesssim 10^{-23.7}$ 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–46]. In this context, ULDM may be regarded as a scalar-tensor theory of the Fierz-Jordan-Brans-Dicke [47–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 [51, 54–58]. This effect, known as Nordtvedt effect [54], violates the strong equivalence principle, and it has long been constrained with binary pulsar data [59–66].

In a recent companion paper [67], following former studies [56, 58, 68], 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 II, we define the Lagrangian of our theory, and we briefly review some general features of ULDM relevant to our analysis. Section III will be devoted to a detailed description of the dataset and the procedure used to carry out the analysis. Finally, in Section IV, 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 [44, 69] or from the pulsar in a triple stellar system [62–64] in the mass region which PTAs are sensitive to. Our conclusions are presented in Section V. The appendices are devoted to technical details and plots.
## II. 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) [31, 70, 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 by <fnref n="1" />:

$$S = M_P^2 \int d^4 x \sqrt{-g} \left[ \frac{R}{2} - g^{\mu\nu} \partial_\mu \phi \partial_\nu \phi + m^2 \phi^2 \right] + S_m[\psi_m, \tilde{g}_{\mu\nu}], \quad (1)$$

where $M_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} = A^2(\phi)g_{\mu\nu}$, normalized such that $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–50], in which the conformal coupling is linear:

$$A(\phi) = e^{\alpha\phi} \sim 1 + \alpha\phi. \quad (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 explore<fnref n="2" />.

As a second model, we also consider the Damour-Esposito-Farëse (DEF) gravity theory [51, 52], where the

conformal coupling is quadratic:

$$A(\phi) = e^{\beta\phi^2/2}. \quad (3)$$

Notice that we have chosen to set the linear coupling of the field in $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 (1):

$$T_{\mu\nu} = M_P^2 \left[ 2\partial_\mu\phi\partial_\nu\phi - g_{\mu\nu} \left( (\partial\phi)^2 - m^2\phi^2 \right) \right]. \quad (4)$$

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

$$\tilde{T}_{\mu\nu} = A^{-2}(\phi)T_{\mu\nu} \approx (1 - 2\alpha(\phi))T_{\mu\nu}, \quad (5)$$

where the effective scalar coupling is<fnref n="3" /> $\alpha(\phi) = d \log A/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}$ 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 wave<fnref n="4" />:

$$\phi(\mathbf{x}, t) = A(\mathbf{x}) \cos(mt + \theta(\mathbf{x})), \quad (6)$$

with $A(\mathbf{x})$ determined by<fnref n="5" />

$$\tilde{\rho}_{\phi} \approx \rho_{\phi} = m^2 M_P^2 A(\mathbf{x})^2 = \rho \hat{\phi}(\mathbf{x})^2, \quad (7)$$

<footnote n="1">¹ Note that our scalar field $\phi$ is not canonically normalized, but appears in the action multiplied by $M_P$. This convention makes comparisons with gravitational phenomena more straightforward.</footnote>

<footnote n="2">² The 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.</footnote>

<footnote n="3">³ In 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].</footnote>

<footnote n="4">⁴ Close 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.</footnote>

<footnote n="5">⁵ Notice 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. (1).</footnote>
where we have used Eq. (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 ~ ($\partial_i \phi)^2$) which are suppressed by $v_\phi^2$. Here, $\rho$ is the average density of the scalar field and $\hat{\phi}(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_{DM} = 0.4$ 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}$ 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_{DM} \equiv \rho/\rho_{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{10^{-22} \text{ eV}}{m} \right), \quad (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). \quad (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. (6), induces a temporal variation of Newton's constant. Its measured value is [51]:

$$ G = (8\pi M_P^2)^{-1} \mathcal{A}^2(\phi)(1 + \alpha^2(\phi)), \quad (10) $$

where $\alpha(\phi) = \mathcal{A}'(\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_{obs}$ through the relation $J = I\Omega_{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(\hat{\phi})} \left. \frac{d \ln I}{d\hat{\phi}} \right|_{N,J} = \frac{1}{2\alpha(\hat{\phi})} \left. \frac{d \ln \Omega_{\text{obs}}}{d\hat{\phi}} \right|_{N,J}, \quad (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 IV to constrain ULDM models.

### III. 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 dataset$^6$, consisting of a 24.7 years collection of TOAs of radio pulses of 25 millisecond pulsars, surveyed with an approximately biweekly cadence$^7$ 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–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_{GWB}$ and a spectral index $\gamma_{GWB}$ (see Table I for more details). Such a signal appears as a precursor to the SGWB [88–90], because of the stronger autocorrelation of the Hellings-Downs process.

We utilize Bayesian inference to detect the ULDM-induced deterministic$^8$ signal, while simultaneously

<footnote>$^6$ The dataset can be found at https://gitlab.in2p3.fr/epta/epta-dr2.</footnote>

<footnote>$^7$ The cadence is non-uniform. EPTA combines observations from several telescopes, so sometimes the EPTA-wide cadence can be much shorter.</footnote>

<footnote>$^8$ The signal induced by ULDM is deterministic, as opposed to the</footnote>
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–95]:

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

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. (19) and (23), 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 I. 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_{\mathrm{obs}} = 24.7$ yr of the EPTA-DR2Full dataset is much shorter than the coherence time in Eq. (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}(\mathbf{x})$ parameter. As the average interpulsar distance is $\mathcal{O}(\mathrm{kpc})$, the *uncorrelated* approximation holds for ULDM masses $m \gtrsim 5 \times 10^{-23}$ eV. For masses $m \lesssim 2 \times 10^{-24}$ eV (*correlated* scenario), the coherence length described by Eq. (9) encompasses the inner Galacto-centric $\sim 20$ 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}(\mathbf{x})$ into the measured value of the local ULDM abundance. Equation (6) then reads:

$$ \hat{\phi}(\mathbf{x}, t) = \frac{\sqrt{\rho}}{m M_{\mathrm{P}}} \cos(mt + \theta(\mathbf{x})), \quad (13) $$

with $\rho$ representing the value of the scalar field density $\rho_\phi$ in our Galaxy. Finally, for masses $2 \times 10^{-24}$ eV $\lesssim m \lesssim 5 \times 10^{-23}$ 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}(\mathbf{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}(\mathbf{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(\mathbf{x})$ per pulsar. This phase encodes the uncertainty on current pulsar distance measurements [40, 99, 100] (see below Eq. (16) in Sec. IV for more details).

We focus on the ULDM mass range $m \in [10^{-24} \mathrm{eV}, 10^{-21} \mathrm{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_{\mathrm{low}} \sim 2.4 \times 10^{-10}$ Hz, which is far below the inverse of EPTA-DR2Full observation length $f_{\mathrm{obs}} = 1/T_{\mathrm{obs}} \sim 1.3$ nHz. In this regime, the ULDM-induced signal (see Eq. (19) 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–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. (19)-(18)). Therefore, we can still provide significant bounds at $f \lesssim 1/T_{\mathrm{obs}}$ [104]. The upper end of the ULDM mass range, instead, corresponds to $f_{\mathrm{up}} \sim 2.4 \times 10^{-7}$ Hz, which is somewhat lower than the observational cadence of $f \sim 1/2$ weeks $\sim 8.3 \times 10^{-7}$ 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.

## IV. RESULTS

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

<footnote>stochastic nature of e.g. the common red noise process describing the SGWB.</footnote>
A. FJBD conformal coupling

Let us focus on FJBD theory. In this case, by inspect-
ing Eq. (2) and recalling the definition of α(φ), we find
α(φ) = α. 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 α
that we neglect in this analysis. Depending on the con-
text, we conveniently write s_I(α, M) as s_I(M) or s_I to
avoid cluttering the notation.

For small scalar fluctuations such as the ones consid-
ered in this article, it follows from Eq. (11) that

$$
\begin{align}
\Omega_{\text{obs}}(t) &= \bar{\Omega} (1 + 2\alpha s_I \delta\phi(t)) \nonumber \\
&= \bar{\Omega} \left( 1 + 2\alpha s_I \frac{\sqrt{\rho}}{M_{\text{P}}m} \hat{\phi}(\mathbf{x}) \cos(mt + \theta(\mathbf{x})) \right), \tag{14}
\end{align}
$$

where we used Eq. (6) and we denoted the spin frequency
of the pulsar *in the absence of* the scalar field by $\bar{\Omega}$.
Notice that Eq. (14) 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$ s [105], which would be larger than $t_{\text{osc}}$ only for ULDM
masses $m \gtrsim 10^{-16}$ eV, which are not discussed here (see
Sec. III).

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

$$
\frac{\delta \Omega_{\text{obs}}}{\Omega} = 2\alpha s_I \delta \phi = - \frac{\delta P}{P}, \quad (15)
$$

where $\tilde{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:

$$
\begin{align}
\Delta t(t) &= - \int dt \frac{\delta \Omega_{\text{obs}}}{\bar{\Omega}} \nonumber \\
&= 2\alpha s_I \frac{\sqrt{\rho}}{M_{\text{P}}m^2} \hat{\phi}(\mathbf{x}) \sin(mt + \theta(\mathbf{x})) \Big|_{t_{\text{start}} - \frac{d}{c}}^{t_{\text{end}} - \frac{d}{c}}, \tag{16}
\end{align}
$$

where we have highlighted the dependence on the re-
tarded time tᵢ - d/c, with i = 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 θ(**x**) → θ(**x**) + md/c. Present uncertainties
on the pulsar distances are O(0.1 ÷ 1) kpc [99], imply-
ing that this redefinition gives rise to an effective pulsar-
dependent random phase. Therefore, as mentioned in
Sec. III, we treat θ(**x**) as a pulsar-specific random pa-
rameter, 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_c}{2m} [\hat{\phi}_E^2 \sin(2mt + \theta_E) - \hat{\phi}_P^2 \sin(2mt + \theta_P)], \quad (17)
$$

with E, P labeling respectively the Earth and the pulsar<fnref n="9" />,
we can define an *effective* amplitude

$$
\Psi = 2\alpha \frac{\sqrt{\rho}}{M_P m}, \quad (18)
$$

such that

$$
\Delta t(t) = \frac{\Psi}{m} s_I \hat{\phi}(\mathbf{x}) \sin(mt + \theta(\mathbf{x})) \Big|_{t_{\text{start}}}^{t_{\text{end}}}. \quad (19)
$$

This form helps to understand what a sensible prior on Ψ
may be. In fact, noticing that Eqs. (17) and (19) have the
same form (differing only for the presence of the Earth
term and some O(1) factors), we can use the same prior
for Ψ and Ψ_c, i.e. Log₁₀ – Uniform [–20, –12]. 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 experi-
  ments as M ± δM, we draw the mass from a normal
distribution centered on M, with uncertainty δM
and truncated for masses below M_min = 1.1 M_☉
and above M_max ~ 2 M_☉. The precise value of
M_max depends on the EoS considered (see Ap-
pendix A for more details);

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

At this point, it is worth obtaining some analytical under-
standing of the results that we expect for α. In particular,
writing the residual induced by an ULDM candidate of
mass m as Δ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), \quad (20)
$$

which, through Eq. (18), yields

$$
\alpha \sim \frac{m M_{\mathrm{P}}}{\sqrt{\rho}} \Psi_{\mathrm{c}}\left(\frac{m}{2}\right) .
\qquad(21)
$$

<footnote n="9">⁹ The stochastic parameter $\hat{\phi}_{\text{P}}$ corresponds to $\hat{\phi}$ in our notation.</footnote>
By substituting the upper limits for $\Psi_c$ found in Fig. 1 of Ref. [40], Eq. (21) 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. (21), assuming for reference $\rho = \rho_{DM} = 0.4$ 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. (17) and Eq. (19) (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_{DM}}$ to take into account the relative energy density of this ULDM candidate $f_{DM}$. This is indeed the quantity constrained by Eq. (18), 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}$ eV and $m \sim 10^{-21.4}$ 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. (1)–(2)) will in general produce both an $\alpha$-dependent residual and a propagation residual, described respectively by Eqs. (16)–(17). 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. (19) is similar to Eq. (12) in Ref. [42], where direct couplings between ULDM and ordinary matter were also studied with PTA data, apart from $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><img src="image_2.png" /><figcaption>FIG. 1: Upper limits on $\log_{10} \alpha$ at 95% credibility compared to the analytical estimate described by Eq. (21), 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_{DM} = 0.4$ GeV/cm$^3$, while the brown dashed lines displays the expected behavior described by Eq. (21). The dotted line shows the degradation of the bounds when choosing $\rho = \min(\bar{\rho}, \rho_{DM})$, optimistically setting $\bar{\rho}$ to the upper bounds presented in Ref. [40]. Smaller $\rho$ yield a stronger degradation of the limits.</figcaption></figure>

## B. DEF theory

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

$$s_I = \frac{1}{2\beta\phi} \left. \frac{d \ln \Omega_{\text{obs}}}{d\phi} \right|_{N,J} \qquad (22)$$

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
<figure><img src="image_3.png" /><figcaption>FIG. 2: Upper limits on $\log_{10}(\alpha\sqrt{f_{DM}})$ 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 I. 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.</figcaption></figure>

then reads:

$$ 
\begin{aligned}
\Delta t(t) &= -\int \frac{\delta\Omega_{\text{obs}}}{\Omega} dt = \\
&= -2\beta s_I(\beta) \int \frac{\rho}{m^2 M_P^2} \hat{\phi}^2(\mathbf{x}) \cos^2(mt + \theta(\mathbf{x})) \\
&= \frac{\Psi}{2m} \beta s_I(\beta) \hat{\phi}^2(\mathbf{x}) \sin(2mt + \theta(\mathbf{x})) \Big|_{t=\text{start}}^{t=\text{end}}, \quad (23)
\end{aligned}
$$

where we used $\Psi = \rho/(m^2 M_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–103]. As $s_I = s_I(\beta)$, Eq. (23) 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_{DM}$ and $\rho = \rho_{DM}$. Moreover, the deterministic residual in Eq. (23) depends on the sign of $\beta$, and not only on its absolute value. While one might naïvely think that a sign flip $\beta \to -\beta$ could be absorbed by a redefinition of the random phase $\theta(\mathbf{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 $|\beta|$ for the correlated, pulsar-correlated and uncorrelated limits and for the selected MPA1 EoS, showing that our analysis implies $|\beta| \lesssim 2.2$ in the range $10^{-24}$ eV $\lesssim m \lesssim 10^{-23.5}$ eV for $\rho = \rho_{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 $|\beta| = 4.3$. As expected from the form of the signal in Eq. (23), 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_{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_{DM}$ for masses $10^{-24}$ eV $\lesssim m \lesssim 10^{-23.7}$ eV and $\rho \gtrsim 0.7 \rho_{DM}$ in the mass range $10^{-23.7}$ eV $\lesssim m \lesssim 10^{-23.4}$ eV. This remark should be taken into account when interpreting the results in Figs. 3-4.
<figure><img src="image_4.png" /><figcaption>FIG. 3: Upper limits on |β| (β < 0) at 95% credibility versus the ULDM mass. The top left panel shows the results for ρ = ρ$_{DM}$ , the top right panel assumes ρ = 0 . 5 ρ$_{DM}$ and the bottom panel displays the results for ρ = 0 . 3 ρ$_{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 I.</figcaption></figure>

V. CONCLUSIONS

Conformally coupled ULDM induces periodic varia-
tions in the gravitational mass and in the radius of pul-
sars [67], with a timescale given by the mass of the par-
ticle. By conservation of angular momentum, this trans-
lates into an oscillating behavior of the pulsar spin fre-
quency, which is accurately measured by the PTA col-
laborations. This effect can be used to set constraints
on the coupling of ULDM to matter, characterized by
the conformal factor A ( φ ) linking the Einstein and Jor-
dan frame metrics. In this work, we have analyzed the
FJBD and the DEF scalar-tensor theories with an ultra-
light scalar mass, under the assumption that the scalar
field constitutes (part of) the DM, thus exploiting differ-
ent functional forms of the conformal factor A ( φ ).

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 [62–64]. Moreover, for masses $m \sim 10^{-23}$ 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
<figure><img src="image_5.png" /></figure>

FIG. 4: Upper limits on $\beta$ ($\beta > 0$) at 95% credibility versus the ULDM mass. The top left panel shows the results for $\rho = \rho_{DM}$, the top right panel assumes $\rho = 0.5 \rho_{DM}$ and the bottom panel displays the results for $\rho = 0.3 \rho_{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_{DM}$, while it is completely unconstraining when $\rho = 0.3 \rho_{DM}$. Priors on the parameters relevant for the search are chosen according to the scheme presented in Table I.

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 (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 Pul-
sar timing", and the EU Horizon 2020 Research and
Innovation Programme under the Marie Skłodowska-
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 Gen-
eralitat de Catalunya. D. Blas acknowledges the sup-
port 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 Na-
tional 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 Agree-
ment number: 101094354) The Nançay radio Observa-
tory is operated by the Paris Observatory, associated
with the French Centre National de la Recherche Sci-
entifique (CNRS), and partially supported by the Re-
gion Centre in France. I. Cognard, L. Guillemot and
G. Theureau acknowledge financial support from “Pro-
gramme 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 ac-
knowledge financial support from Agence Nationale de
la Recherche (ANR-18-CE31-0015), France. Pulsar re-
search at Jodrell Bank Centre for Astrophysics is sup-
ported 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 “GU-
VIRP - Gravity tests in the UltraViolet and InfraRed
with Pulsar timing” and the INAF 2023 Large Grant
“Gravitational Wave Detection using Pulsar Timing Ar-
rays”.
## Appendix A: Parameters of the search

Table I 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* ± *δM*, we sample the mass parameter from a normal prior distribution centered on *M*, with uncertainty *δM* and truncated below *M*<sub>min</sub> = 1.1 *M*<sub>☉</sub> and above *M*<sub>max</sub> = (2.2 *M*<sub>☉</sub>, 2.4 *M*<sub>☉</sub>, 2.2 *M*<sub>☉</sub>) for the (AP4 [110], MPA1 [111], SLy [112]) EoS, respectively. We will label this as TruncNorm(*μ*, *σ*), where *μ* = *M* and *σ* = *δM*. Otherwise, we will just assume a uniform prior between *M*<sub>min</sub> and *M*<sub>max</sub> for the three EoS choices. We also plot the constraints on the FJBD scalar coupling *α* for the different EoS of the pulsar interior presented in Ref. [67].

<table><caption>TABLE I: Parameters employed for the search along with their respective priors. In the correlated limit, the parameter $\hat{\phi}_P^2$ is accounted for by a redefinition of $\Psi$, while in the pulsar-correlated regime $\hat{\phi}_P^2 = \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*<sub>min</sub>, *M*<sub>max</sub>] (see main text for more details). We display the priors on the DEF scalar coupling $\beta$ both for $\beta < 0$ and $\beta > 0$. *U* stands for the uniform distribution and *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.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th><th>Occurrence</th></tr></thead><tbody><tr><td colspan="4">White Noise (σ = E<sub>f</sub><sup>2</sup>σ<sub>TOA</sub><sup>2</sup> + E<sub>q</sub><sup>2</sup>)</td></tr><tr><td>E<sub>f</sub></td><td>EFAC</td><td>U(0, 10)</td><td>1 per pulsar</td></tr><tr><td>E<sub>q</sub></td><td>EQUAD</td><td>Log<sub>10</sub>-U(-10, -5)</td><td>1 per pulsar</td></tr><tr><td colspan="4">Red Noise (RN)</td></tr><tr><td>A<sub>red</sub></td><td>RN power-law amplitude</td><td>Log<sub>10</sub>-U(-20, -11)</td><td>1 per pulsar</td></tr><tr><td>γ<sub>red</sub></td><td>RN power-law spectral index</td><td>U(0, 10)</td><td>1 per pulsar</td></tr><tr><td colspan="4">ULDM (F)</td></tr><tr><td>Ψ</td><td>ULDM signal amplitude</td><td>Log<sub>10</sub>-U(-20, -12)</td><td>1 per PTA</td></tr><tr><td>m [eV]</td><td>ULDM mass</td><td>Log<sub>10</sub>-U(-24, -21)</td><td>1 per PTA</td></tr><tr><td>φ̂<sup>2</sup></td><td>Pulsar factor</td><td>e<sup>-x</sup></td><td>1 per pulsar</td></tr><tr><td>θ</td><td>Pulsar signal phase</td><td>U(0, 2π)</td><td>1 per pulsar</td></tr><tr><td colspan="4">ULDM (D)</td></tr><tr><td>f<sub>DM</sub></td><td>ULDM fraction</td><td>U(0.01, 0.30)</td><td>1 per PTA</td></tr><tr><td>β</td><td>DEF scalar coupling</td><td>U(-4.3, 0)<br>or<br>U(0, 150)</td><td>1 per PTA</td></tr><tr><td>m [eV]</td><td>ULDM mass</td><td>Log<sub>10</sub>-U(-24, -21)</td><td>1 per PTA</td></tr><tr><td>φ̂<sup>2</sup></td><td>Pulsar factor</td><td>e<sup>-x</sup></td><td>1 per pulsar</td></tr><tr><td>θ</td><td>Pulsar signal phase</td><td>U(0, 2π)</td><td>1 per pulsar</td></tr><tr><td colspan="4">Common spatially Uncorrelated Red Noise (CURN)</td></tr><tr><td>A<sub>GWB</sub></td><td>CURN strain amplitude</td><td>Log<sub>10</sub>-U(-18, -11)</td><td>1 per PTA</td></tr><tr><td>γ<sub>GWB</sub></td><td>CURN spectral index</td><td>U(0, 7)</td><td>1 per PTA</td></tr><tr><td colspan="4">Pulsar Masses</td></tr><tr><td>M<sub>J0030</sub> [113]</td><td>PSR J0030+0451 mass</td><td>N(1.44, 0.15)</td><td>1 per PTA</td></tr><tr><td>M<sub>J1012</sub> [114]</td><td>PSR J1012+5307 mass</td><td>N(1.72, 0.16)</td><td>1 per PTA</td></tr><tr><td>M<sub>J1713</sub> [115]</td><td>PSR J1713+0747 mass</td><td>N(1.3, 0.2)</td><td>1 per PTA</td></tr><tr><td>M<sub>J1738</sub> [116]</td><td>PSR J1738+0333 mass</td><td>N(1.47, 0.07)</td><td>1 per PTA</td></tr><tr><td>M<sub>J1909</sub> [117]</td><td>PSR J1909-3744 mass</td><td>N(1.438, 0.024)</td><td>1 per PTA</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>FIG. 5: Upper limits on $log_{10} \alpha$ at 95% credibility versus the ULDM mass, using the prescription for $\rho_{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 I.</figcaption></figure>

[1] Ricardo A. Flores and Joel R. Primack, “Observational and Theoretical Constraints on Singular Dark Matter Halos,” ApJ 427, L1 (1994), arXiv:astro-ph/9402004 [astro-ph].

[2] Ben Moore, “Evidence against dissipation-less dark matter from observations of galaxy haloes,” Nature 370, 629–631 (1994).

[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, A13 (2015).

[4] Anatoly Klypin, Andrey V. Kravtsov, Octavio Valenzuela, and Francisco Prada, “Where are the missing galactic satellites?” The Astrophysical Journal 522, 82–92 (1999).

[5] Ben Moore, Sebastiano Ghigna, Fabio Governato, George Lake, Thomas Quinn, Joachim Stadel, and Paolo Tozzi, “Dark Matter Substructure within Galactic Halos,” ApJ 524, L19–L22 (1999), arXiv:astro-ph/9907411 [astro-ph].

[6] Michael Boylan-Kolchin, James S. Bullock, and Manoj Kaplinghat, “Too big to fail? The puzzling darkness of massive Milky Way subhaloes,” MNRAS 415, L40–L44 (2011), arXiv:1103.0007 [astro-ph.CO].

[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, L72–L78 (1996).

[8] F. Governato, A. Zolotov, A. Pontzen, C. Christensen, S. H. Oh, A. M. Brooks, T. Quinn, S. Shen, and J. Wadsley, “Cuspy no more: how outflows affect the central dark matter and baryon distribution in Λ cold dark matter galaxies,” Monthly Notices of the Royal Astronomical Society 422, 1231–1240 (2012), https://academic.oup.com/mnras/article-pdf/422/2/1231/3466602/mnras0422-1231.pdf.

[9] Alyson M. Brooks, Michael Kuhlen, Adi Zolotov, and Dan Hooper, “A baryonic solution to the missing satellites problem,” The Astrophysical Journal 765, 22 (2013).

[10] T. K. Chan, D. Kereš, J. Oñorbe, P. F. Hopkins, A. L. Muratov, C.-A. Faucher-Giguère, and E. Quataert, “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, 2981–3001 (2015), https://academic.oup.com/mnras/article-pdf/454/3/2981/4038253/stv2165.pdf.

[11] Jose Oñorbe, Michael Boylan-Kolchin, James S. Bullock, Philip F. Hopkins, Dušan Kereš, Claude-André Faucher-Giguère, Eliot Quataert, and Norman Murray, “Forged in fire: cusps, cores and baryons in low-mass dwarf galaxies,” Monthly Notices of the Royal Astronomical Society 454, 2092–2106 (2015), 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, 2573–2590 (2016), https://academic.oup.com/mnras/article-pdf/459/3/2573/8105757/stw713.pdf.

[13] Wayne Hu, Rennan Barkana, and Andrei Gruzinov, “Cold and fuzzy dark matter,” Phys. Rev. Lett. 85, 1158–1161 (2000), arXiv:astro-ph/0003365.

[14] Lam Hui, Jeremiah P. Ostriker, Scott Tremaine, and Edward Witten, “Ultralight scalars as cosmological dark matter,” Physical Review D 95 (2017), 10.1103/physrevd.95.043541.
[15] Michael B. Green, J. H. Schwarz, and Edward Witten, *SUPERSTRING THEORY. VOL. 1: INTRODUCTION*, Cambridge Monographs on Mathematical Physics (1988).

[16] Peter Svrcek and Edward Witten, “Axions in string theory,” Journal of High Energy Physics **2006**, 051 (2006).

[17] Asimina Arvanitaki, Savas Dimopoulos, Sergei Dubovsky, Nemanja Kaloper, and John March-Russell, “String axiverse,” Phys. Rev. D **81**, 123530 (2010).

[18] Renée Hlozek, Daniel Grin, David J. E. Marsh, and Pedro G. Ferreira, “A search for ultralight axions using precision cosmological data,” Phys. Rev. D **91**, 103512 (2015).

[19] Vid Iršič, Matteo Viel, Martin G. Haehnelt, James S. Bolton, and George D. Becker, “First constraints on fuzzy dark matter from lyman-alpha forest data and hydrodynamical simulations,” Phys. Rev. Lett. **119**, 031302 (2017).

[20] Eric Armengaud, Nathalie Palanque-Delabrouille, Christophe Yèche, David J. E. Marsh, and Julien Baur, “Constraining the mass of light bosonic dark matter using SDSS Lyman-alpha forest,” Monthly Notices of the Royal Astronomical Society **471**, 4606–4614 (2017), https://academic.oup.com/mnras/article-pdf/471/4/4606/1963524/stx1870.pdf.

[21] Takeshi Kobayashi, Riccardo Murgia, Andrea De Simone, Vid Iršič, and Matteo Viel, “Lyman-alpha constraints on ultralight scalar dark matter: Implications for the early and late universe,” Physical Review D **96** (2017), 10.1103/phys-rev.d.96.123514.

[22] Matteo Nori, Riccardo Murgia, Vid Iršič, Marco Baldi, and Matteo Viel, “Lyman-alpha forest and non-linear structure characterization in Fuzzy Dark Matter cosmologies,” Monthly Notices of the Royal Astronomical Society **482**, 3227–3243 (2018), https://academic.oup.com/mnras/article-pdf/482/3/3227/26653692/sty2888.pdf.

[23] Keir K. Rogers and Hiranya V. Peiris, “Strong bound on canonical ultralight axion dark matter from the lyman-alpha forest,” Phys. Rev. Lett. **126**, 071302 (2021).

[24] Hsi-Yu Schive, Ming-Hsuan Liao, Tak-Pong Woo, Shing-Kwong Wong, Tzihong Chiueh, Tom Broadhurst, and W-Y. Pauchy Hwang, “Understanding the core-halo relation of quantum wave dark matter from 3d simulations,” Phys. Rev. Lett. **113**, 261302 (2014).

[25] Jiajun Zhang, Hantao Liu, and Ming-Chung Chu, “Cosmological simulation for fuzzy dark matter model,” Frontiers in Astronomy and Space Sciences **5** (2019), 10.3389/fspas.2018.00048.

[26] Nitsan Bar, Diego Blas, Kfir Blum, and Sergey Sibiryakov, “Galactic rotation curves versus ultralight dark matter: Implications of the soliton-host halo relation,” Phys. Rev. D **98**, 083027 (2018), arXiv:1805.00122 [astro-ph.CO].

[27] Kohei Hayashi, Elisa G. M. Ferreira, and Hei Yin Jowett Chan, “Narrowing the mass range of fuzzy dark matter with ultrafaint dwarfs,” The Astrophysical Journal Letters **912**, L3 (2021).

[28] Neal Dalal and Andrey Kravtsov, “Not so fuzzy: excluding fdm with sizes and stellar kinematics of ultra-faint dwarf galaxies,” (2022), arXiv:2203.05750 [astro-ph.CO].

[29] David J. E. Marsh and Jens C. Niemeyer, “Strong Constraints on Fuzzy Dark Matter from Ultrafaint Dwarf Galaxy Eridanus II,” Phys. Rev. Lett. **123**, 051103 (2019), arXiv:1810.08543 [astro-ph.CO].

[30] Elisa G. M. Ferreira, “Ultra-light dark matter,” Astron. Astrophys. Rev. **29**, 7 (2021), arXiv:2005.03254 [astro-ph.CO].

[31] Asimina Arvanitaki, Savas Dimopoulos, Sergei Dubovsky, Nemanja Kaloper, and John March-Russell, “String Axiverse,” Phys. Rev. D **81**, 123530 (2010), arXiv:0905.4720 [hep-th].

[32] Louis Hamaide, Hendrik Müller, and David J. E. Marsh, “Searching for dilaton fields in the Lyman-α forest,” Phys. Rev. D **106**, 123509 (2022), arXiv:2210.03705 [astro-ph.CO].

[33] Andrei Khmelnitsky and Valery Rubakov, “Pulsar timing signal from ultralight scalar dark matter,” Journal of Cosmology and Astroparticle Physics **2014**, 019–019 (2014).

[34] R. S. Foster and D. C. Backer, “Constructing a Pulsar Timing Array,” ApJ **361**, 300 (1990).

[35] R. N. Manchester, G. Hobbs, M. Bailes, W. A. Coles, W. van Straten, M. J. Keith, R. M. Shannon, N. D. R. Bhat, A. Brown, S. G. Burke-Spolaor, D. J. Champion, A. Chaudhary, R. T. Edwards, G. Hampson, A. W. Hotan, A. Jameson, F. A. Jenet, M. J. Kesteven, J. Khoo, J. Kocz, K. Maciesiak, S. Oslowski, V. Ravi, J. R. Reynolds, J. M. Sarkissian, J. P. W. Verbiest, Z. L. Wen, W. E. Wilson, D. Yardley, W. M. Yan, and X. P. You, “The Parkes Pulsar Timing Array Project,” PASA **30**, e017 (2013), arXiv:1210.6130 [astro-ph.IM].

[36] M. A. McLaughlin, “The North American Nanohertz Observatory for Gravitational Waves,” Classical and Quantum Gravity **30**, 224008 (2013), arXiv:1310.0758 [astro-ph.IM].

[37] Michael Kramer and David J. Champion, “The European Pulsar Timing Array and the Large European Array for Pulsars,” Classical and Quantum Gravity **30**, 224009 (2013).

[38] Bhal Chandra Joshi, Prakash Arumugasamy, Manjari Bagchi, Debades Bandyopadhyay, Avishek Basu, Neelam Dhanda Batra, Suryarao Bethapudi, Arpita Choudhary, Kishalay De, L. Dey, A. Gopakumar, Y. Gupta, M. A. Krishnakumar, Yogesh Maan, P. K. Manoharan, Arun Naidu, Rana Nandi, Dhruv Pathak, Mayuresh Surnis, and Abhimanyu Susobhanan, “Precision pulsar timing with the ORT and the GMRT and its applications in pulsar astrophysics,” Journal of Astrophysics and Astronomy **39**, 51 (2018).

[39] M. Bailes, A. Jameson, F. Abbate, E. D. Barr, N. D. R. Bhat, L. Bondonneau, M. Burgay, S. J. Buchner, F. Camilo, D. J. Champion, I. Cognard, P. B. Demorest, P. C. C. Freire, T. Gautam, M. Geyer, J. M. Griessmeier, L. Guillemot, H. Hu, F. Jankowski, S. Johnston, A. Karastergiou, R. Karuppusamy, D. Kaur, M. J. Keith, M. Kramer, J. van Leeuwen, M. E. Lower, Y. Maan, M. A. McLaughlin, B. W. Meyers, S. Oslowski, L. S. Oswald, A. Parthasarathy, T. Pennucci, B. Posselt, A. Possenti, S. M. Ransom, D. J. Reardon, A. Ridolfi, C. T. G. Schollar, M. Serylak, G. Shaifullah, M. Shamohammadi, R. M. Shannon, C. Sobey, X. Song, R. Spiewak, I. H. Stairs, B. W. Stappers, W. van Straten, A. Szary, G. Theureau, V. Venkatraman Krishnan, P. Weltevrede, N. Wex, T. D. Abbott, G. B. Adams, J. P. Burger, R. R. G. Gamatham, M. Gouws, D. M. Horn, B. Hugo, A. F. Joubert, J. R. Manley, K. McAlpine, S. S. Passmoor, A. Peens-Hough, Z. R. Ramudzuli, A. Rust, S. Salie, L. C. Schwardt, R. Siebrits, G. Van Tonder, V. Van Tonder, and M. G. Welz, “The
MeerKAT telescope as a pulsar facility: System verification and early science results from MeerTime," PASA 37, e028 (2020), arXiv:2005.14366 [astro-ph.IM].

[40] EPTA Collaboration, Clemente Smarra, Boris Goncharov, Enrico Barausse, J. Antoniadis, S. Babak, A.-S. Bak Nielsen, C. G. Bassa, A. Berthereau, M. Bonetti, E. Bortolas, P. R. Brook, M. Burgay, R. N. Caballero, A. Chalumeau, D. J. Champion, S. Chanlaridis, S. Chen, I. Cognard, G. Desvignes, M. Falxa, R. D. Ferdman, A. Franchini, J. R. Gair, E. Graikou, J.-M. Grießmeier, L. Guillemot, Y. J. Guo, H. Hu, F. Iraci, D. Izquierdo-Villalba, J. Jang, J. Jawor, G. H. Janssen, A. Jessner, R. Karuppusamy, E. F. Keane, M. J. Keith, M. Kramer, M. A. Krishnakumar, K. Lackeos, K. J. Lee, K. Liu, Y. Liu, A. G. Lyne, J. W. McKee, R. A. Main, M. B. Mickaliger, I. C. Nițu, A. Parthasarathy, B. B. P. Perera, D. Perrodin, A. Petiteau, N. K. Porayko, A. Possenti, H. Quelquejay Leclere, A. Samajdar, S. A. Sanidas, A. Sesana, G. Shaifullah, L. Speri, R. Spiewak, B. W. Stappers, S. C. Susarla, G. Theureau, C. Tiburzi, E. van der Wateren, A. Vecchio, V. Venkatraman Krishnan, J. Wang, L. Wang, and Z. Wu, "Second data release from the european pulsar timing array: Challenging the ultralight dark matter paradigm," Phys. Rev. Lett. 131, 171001 (2023).

[41] Nataliya K. Porayko, Xingjiang Zhu, Yuri Levin, Lam Hui, George Hobbs, Aleksandra Grudskaya, Konstantin Postnov, Matthew Bailes, N.D. Ramesh Bhat, William Coles, Shi Dai, James Dempsey, Michael J. Keith, Matthew Kerr, Michael Kramer, Paul D. Lasky, Richard N. Manchester, Stefan Osłowski, Aditya Parthasarathy, Vikram Ravi, Daniel J. Reardon, Pablo A. Rosado, Christopher J. Russell, Ryan M. Shannon, Renée Spiewak, Willem van Straten, Lawrence Toomey, Jingbo Wang, Linqing Wen, and Xiaopeng You, "Parkes pulsar timing array constraints on ultralight scalar-field dark matter," Physical Review D 98 (2018), 10.1103/physrevd.98.102002.

[42] David E. Kaplan, Andrea Mitridate, and Tanner Trickle, "Constraining fundamental constant variations from ultralight dark matter with pulsar timing arrays," Phys. Rev. D 106, 035032 (2022).

[43] Adeela Afzal, Gabriella Agazie, Akash Anumarlapudi, Anne M. Archibald, Zaven Arzoumanian, Paul T. Baker, Bence Bécsy, Jose Juan Blanco-Pillado, Laura Blecha, Kimberly K. Boddy, Adam Brazier, Paul R. Brook, Sarah Burke-Spolaor, Rand Burnette, Robin Case, Maria Charisi, Shami Chatterjee, Katerina Chatzioioannou, Belinda D. Cheeseboro, Siyuan Chen, Tyler Cohen, James M. Cordes, Neil J. Cornish, Fronefield Crawford, H. Thankful Cromartie, Kathryn Crowter, Curt J. Cutler, Megan E. DeCesar, Dallas DeGan, Paul B. Demorest, Heling Deng, Timothy Dolch, Brendan Drachler, Richard von Eckardstein, Elizabeth C. Ferrara, William Fiore, Emmanuel Fonseca, Gabriel E. Freedman, Nate Garver-Daniels, Peter A. Gentile, Kyle A. Gersbach, Joseph Glaser, Deborah C. Good, Lydia Guertin, Kayhan Gültekin, Jeffrey S. Hazboun, Sophie Hourihane, Kristina Islo, Ross J. Jennings, Aaron D. Johnson, Megan L. Jones, Andrew R. Kaiser, David L. Kaplan, Luke Zoltan Kelley, Matthew Kerr, Joey S. Key, Nima Laal, Michael T. Lam, William G. Lamb, T. Joseph W. Lazio, Vincent S. H. Lee, Natalia Lewandowska, Rafael R. Lino dos Santos, Tyson B. Littenberg, Tingting Liu, Duncan R. Lorimer, Jing Luo, Ryan S. Lynch, Chung-Pei Ma, Dustin R. Madison, Alexander McEwen, James W. McKee, Maura A. McLaughlin, Natasha McMann, Bradley W. Meyers, Patrick M. Meyers, Chiara M. F. Mingarelli, Andrea Mitridate, Jonathan Nay, Priyamvada Natarajan, Cherry Ng, David J. Nice, Stella Koch Ocker, Ken D. Olum, Timothy T. Pennucci, Benetge B. P. Perera, Polina Petrov, Nihan S. Pol, Henri A. Radovan, Scott M. Ransom, Paul S. Ray, Joseph D. Romano, Shashwat C. Sardesai, Ann Schmiedekamp, Carl Schmiedekamp, Kai Schmitz, Tobias Schröder, Levi Schult, Brent J. Shapiro-Albert, Xavier Siemens, Joseph Simon, Magdalena S. Siwek, Ingrid H. Stairs, Daniel R. Stinebring, Kevin Stovall, Peter Stratmann, Jerry P. Sun, Abhimanyu Susobhanan, Joseph K. Swiggum, Jacob Taylor, Stephen R. Taylor, Tanner Trickle, Jacob E. Turner, Caner Unal, Michele Vallisneri, Sonali Verma, Sarah J. Vigeland, Haley M. Wahl, Qiaohong Wang, Caitlin A. Witt, David Wright, Olivia Young, Kathryn M. Zurek, and The NANOGrav Collaboration, "The nanograv 15 yr data set: Search for signals from new physics," The Astrophysical Journal Letters 951, L11 (2023).

[44] Diego Blas, Diana Ló pez Nacir, and Sergey Sibiryakov, "Ultralight dark matter resonates with binary pulsars," Physical Review Letters 118 (2017), 10.1103/physrevlett.118.261102.

[45] Diego Blas, Diana Ló pez Nacir, and Sergey Sibiryakov, "Secular effects of ultralight dark matter on binary pulsars," Physical Review D 101 (2020), 10.1103/physrevd.101.063016.

[46] Pavel Kůs, Diana López Nacir, and Federico R. Urban, "Bayesian sensitivity of binary pulsars to ultra-light dark matter," (2024), arXiv:2402.04099 [astro-ph.HE].

[47] M. Fierz, "On the physical interpretation of P.Jordan's extended theory of gravitation," Helv. Phys. Acta 29, 128-134 (1956).

[48] Pascual Jordan, "The present state of Dirac's cosmological hypothesis," Z. Phys. 157, 112-121 (1959).

[49] C. Brans and R. H. Dicke, "Mach's principle and a relativistic theory of gravitation," Phys. Rev. 124, 925-935 (1961).

[50] R. H. Dicke, "Mach's principle and invariance under transformation of units," Phys. Rev. 125, 2163-2167 (1962).

[51] Thibault Damour and Gilles Esposito-Farese, "Tensor multiscalar theories of gravitation," Class. Quant. Grav. 9, 2093-2176 (1992).

[52] Thibault Damour and Gilles Esposito-Farèse, "Nonperturbative strong-field effects in tensor-scalar theories of gravitation," Phys. Rev. Lett. 70, 2220-2223 (1993).

[53] Justin Alsing, Emanuele Berti, Clifford M. Will, and Helmut Zaglauer, "Gravitational radiation from compact binary systems in the massive Brans-Dicke theory of gravity," Phys. Rev. D 85, 064041 (2012), arXiv:1112.4903 [gr-qc].

[54] Kenneth Nordtvedt, "Equivalence principle for massive bodies. ii. theory," Phys. Rev. 169, 1017-1025 (1968).

[55] D. M. Eardley, "Observable effects of a scalar gravitational field in a binary pulsar," Astrophysical Journal 196, L59-L62 (1975).

[56] Clifford M. Will and Helmut W. Zaglauer, "Gravitational Radiation, Close Binary Systems, and the Brans-Dicke Theory of Gravity," ApJ 346, 366 (1989).
[57] C. M. Will, “Gravitational radiation from binary systems in alternative metric theories of gravity: dipole radiation and the binary pulsar.” ApJ 214, 826–839 (1977).

[58] C. M. Will, *Theory and experiment in gravitational physics* (1993).

[59] Thibault Damour and Joseph H. Taylor, “Strong field tests of relativistic gravity and binary pulsars,” Phys. Rev. D 45, 1840–1868 (1992).

[60] J. M. Weisberg and J. H. Taylor, “Relativistic binary pulsar b1913+16: Thirty years of observations and analysis,” (2004), arXiv:astro-ph/0407149 [astro-ph].

[61] M. Kramer, Ingrid H. Stairs, R.N. Manchester, M.A. McLaughlin, A.G. Lyne, et al., “Tests of general relativity from timing the double pulsar,” Science 314, 97–102 (2006), arXiv:astro-ph/0609417 [astro-ph].

[62] S. M. Ransom, I. H. Stairs, A. M. Archibald, J. W. T. Hessels, D. L. Kaplan, M. H. van Kerkwijk, J. Boyles, A. T. Deller, S. Chatterjee, A. Schechtman-Rook, A. Berndsen, R. S. Lynch, D. R. Lorimer, C. Karako-Argaman, V. M. Kaspi, V. I. Kondratiev, M. A. McLaughlin, J. van Leeuwen, R. Rosen, M. S. E. Roberts, and K. Stovall, “A millisecond pulsar in a stellar triple system,” Nature 505, 520–524 (2014).

[63] Anne M. Archibald, Nina V. Gusinskaia, Jason W. T. Hessels, Adam T. Deller, David L. Kaplan, Duncan R. Lorimer, Ryan S. Lynch, Scott M. Ransom, and Ingrid H. Stairs, “Universality of free fall from the orbital motion of a pulsar in a stellar triple system,” Nature 559, 73–76 (2018).

[64] Voisin, G., Cognard, I., Freire, P. C. C., Wex, N., Guillemot, L., Desvignes, G., Kramer, M., and Theureau, G., “An improved test of the strong equivalence principle with the pulsar in a triple star system,” A& A 638, A24 (2020).

[65] M. Kramer, I. H. Stairs, R. N. Manchester, N. Wex, A. T. Deller, W. A. Coles, M. Ali, M. Burgay, F. Camilo, I. Cognard, T. Damour, G. Desvignes, R. D. Ferdman, P. C. C. Freire, S. Grondin, L. Guillemot, G. B. Hobbs, G. Janssen, R. Karuppusamy, D. R. Lorimer, A. G. Lyne, J. W. McKee, M. McLaughlin, L. E. Münch, B. B. P. Perera, N. Pol, A. Possenti, J. Sarkissian, B. W. Stappers, and G. Theureau, “Strong-field gravity tests with the double pulsar,” Phys. Rev. X 11, 041050 (2021).

[66] Junjie Zhao, Paulo C. C. Freire, Michael Kramer, Lijing Shao, and Norbert Wex, “Closing a spontaneous-scalarization window with binary pulsars,” Class. Quant. Grav. 39, 11LT01 (2022), arXiv:2201.03771 [astro-ph.HE].

[67] Adrien Kuntz and Enrico Barausse, “Angular momentum sensitivities in scalar-tensor theories,” (2024), arXiv:2403.07980 [gr-qc].

[68] Thibault Damour and Gilles Esposito-Farèse, “Tensor-scalar gravity and binary-pulsar experiments,” Phys. Rev. D 54, 1474–1491 (1996).

[69] B. Bertotti, L. Iess, and P. Tortora, “A test of general relativity using radio links with the cassini spacecraft,” Nature 425, 374–376 (2003).

[70] T. Damour and Alexander M. Polyakov, “The String dilaton and a least coupling principle,” Nucl. Phys. B 423, 532–558 (1994), arXiv:hep-th/9401069.

[71] Jihn E. Kim, “Light Pseudoscalars, Particle Physics and Cosmology,” Phys. Rept. 150, 1–177 (1987).

[72] Clifford M. Will, *Theory and Experiment in Gravitational Physics* (Cambridge University Press, 1993).

[73] Lijing Shao, Noah Sennett, Alessandra Buonanno, Michael Kramer, and Norbert Wex, “Constraining nonperturbative strong-field effects in scalar-tensor gravity by combining pulsar timing and laser-interferometer gravitational-wave detectors,” Phys. Rev. X 7, 041025 (2017), arXiv:1704.07561 [gr-qc].

[74] Andrés Castillo, Jorge Martin-Camalich, Jorge Terol-Calvo, Diego Blas, Andrea Caputo, Ricardo Tanausú Génova Santos, Laura Sberna, Michael Peel, and Jose Alberto Rubiño Martín, “Searching for dark-matter waves with PPTA and QUIJOTE pulsar polarimetry,” JCAP 06, 014 (2022), arXiv:2201.03422 [astro-ph.CO].

[75] Fabrizio Nesti and Paolo Salucci, “The Dark Matter halo of the Milky Way, AD 2013,” JCAP 07, 016 (2013), arXiv:1304.5127 [astro-ph.GA].

[76] EPTA Collaboration, Antoniadis, J., Babak, S., Bak Nielsen, A.-S., Bassa, C. G., Berthereau, A., Bonetti, M., Bortolas, E., Brook, P. R., Burgay, M., Caballero, R. N., Chalumeau, A., Champion, D. J., Clianlaridis, S., Chen, S., Cognard, I., Desvignes, G., Falxa, M., Ferdman, R. D., Franchini, A., Gair, J. R., Goncharov, B., Graikou, E., Grießmeier, J.-M., Guillemot, L., Guo, Y. J., Hu, H., Iraci, F., Izquierdo-Villalba, D., Jang, J., Jawor, J., Janssen, G. H., Jessner, A., Karuppusamy, R., Keane, E. F., Keith, M. J., Kramer, M., Krishnakumar, M. A., Lackeos, K., Lee, K. J., Liu, K., Liu, Y., Lyne, A. G., McKee, J. W., Main, R. A., Mickaliger, M. B., Niţu, I. C., Parthasarathy, A., Perera, B. B. P., Perrodin, D., Petiteau, A., Porayko, N. K., Possenti, A., Quelquejay Leclere, H., Samajdar, A., Sanidas, S. A., Sesana, A., Shaifullah, G., Speri, L., Spiewak, R., Stappers, B. W., Susarla, S. C., Theureau, G., Tiburzi, C., van der Wateren, E., Vecchio, A., Venkatraman Krishnan, V., Verbiest, J. P. W., Wang, J., Wang, L., and Wu, Z., “The second data release from the european pulsar timing array - i. the dataset and timing analysis,” A& A 678, A48 (2023).

[77] G. Desvignes, R. N. Caballero, L. Lentati, J. P. W. Verbiest, D. J. Champion, B. W. Stappers, G. H. Janssen, P. Lazarus, S. Oslowski, S. Babak, C. G. Bassa, P. Brem, M. Burgay, I. Cognard, J. R. Gair, E. Graikou, L. Guillemot, J. W. T. Hessels, A. Jessner, C. Jordan, R. Karuppusamy, M. Kramer, A. Lassus, K. Lazaridis, K. J. Lee, K. Liu, A. G. Lyne, J. McKee, C. M. F. Mingarelli, D. Perrodin, A. Petiteau, A. Possenti, M. B. Purver, P. A. Rosado, S. Sanidas, A. Sesana, G. Shaifullah, R. Smits, S. R. Taylor, G. Theureau, C. Tiburzi, R. van Haasteren, and A. Vecchio, “High-precision timing of 42 millisecond pulsars with the European Pulsar Timing Array,” MNRAS 458, 3341–3380 (2016), arXiv:1602.08511 [astro-ph.HE].

[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, 1549–1574 (2006), arXiv:astro-ph/0607664 [astro-ph].

[79] Boris Goncharov, D. J. Reardon, R. M. Shannon, Xing-Jiang Zhu, Eric Thrane, M. Bailes, N. D. R. Bhat, S. Dai, G. Hobbs, M. Kerr, R. N. Manchester, S. Oslowski, A. Parthasarathy, C. J. Russell, R. Spiewak, N. Thyagarajan, and
J. B. Wang, "Identifying and mitigating noise sources in precision pulsar timing data sets," MNRAS 502, 478-493 (2021), arXiv:2010.06109 [astro-ph.HE].

[80] EPTA Collaboration and InPTA Collaboration, Antoniadis, J., Arumugam, P., Arumugam, S., Babak, S., Bagchi, M., Nielsen, A.-S. Bak, Bassa, C. G., Bathula, A., Berthereau, A., Bonetti, M., Bortolas, E., Brook, P. R., Burgay, M., Caballero, R. N., Chalumeau, A., Champion, D. J., Chanlaridis, S., Chen, S., Cognard, I., Dandapat, S., Deb, D., Desai, S., Desvignes, G., Dhanda-Batra, N., Dwivedi, C., Falxa, M., Ferdman, R. D., Franchini, A., Gair, J. R., Goncharov, B., Gopakumar, A., Graikou, E., Grießmeier, J.-M., Guillemot, L., Guo, Y. J., Gupta, Y., Hisano, S., Hu, H., Iraci, F., Izquierdo-Villalba, D., Jang, J., Jawor, J., Janssen, G. H., Jessner, A., Joshi, B. C., Kareem, F., Karuppusamy, R., Keane, E. F., Keith, M. J., Kharbanda, D., Kikunaga, T., Kolhe, N., Kramer, M., Krishnakumar, M. A., Lackeos, K., Lee, K. J., Liu, K., Liu, Y., Lyne, A. G., McKee, J. W., Maan, Y., Main, R. A., Mickaliger, M. B., Nițu, I. C., Nobleson, K., Paladi, A. K., Parthasarathy, A., Perera, B. B. P., Perrodin, D., Petiteau, A., Porayko, N. K., Possenti, A., Prabu, T., Leclere, H. Quelquejay, Rana, P., Samajdar, A., Sanidas, S. A., Sesana, A., Shaifullah, G., Singha, J., Speri, L., Spiewak, R., Srivastava, A., Stappers, B. W., Surnis, M., Susarla, S. C., Susobhanan, A., Takahashi, K., Tarafdar, P., Theureau, G., Tiburzi, C., van der Wateren, E., Vecchio, A., Krishnan, V. Venkatraman, Verbiest, J. P. W., Wang, J., Wang, L., and Wu, Z., "The second data release from the european pulsar timing array - ii. customised pulsar noise models for spatially correlated gravitational waves," A&A 678, A49 (2023).

[81] Gabriella Agazie et al. (NANOGrav), "The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background," Astrophys. J. Lett. 951, L8 (2023), arXiv:2306.16213 [astro-ph.HE].

[82] J. Antoniadis et al. (EPTA, InPTA):, "The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals," Astron. Astrophys. 678, A50 (2023), arXiv:2306.16214 [astro-ph.HE].

[83] Daniel J. Reardon et al., "Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array," Astrophys. J. Lett. 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE].

[84] Heng Xu, Siyuan Chen, Yanjun Guo, Jinchen Jiang, Bojun Wang, Jiangwei Xu, Zihan Xue, R. Nicolas Caballero, Jianping Yuan, Yonghua Xu, Jingbo Wang, Longfei Hao, Jingtao Luo, Kejia Lee, Jinlin Han, Peng Jiang, Zhiqiang Shen, Min Wang, Na Wang, Renxin Xu, Xiangping Wu, Richard Manchester, Lei Qian, Xin Guan, Menglin Huang, Chun Sun, and Yan Zhu, "Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I," Research in Astronomy and Astrophysics 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE].

[85] R. W. Hellings and G. S. Downs, "Upper limits on the isotropic gravitational radiation background from pulsar timing analysis." ApJ 265, L39-L42 (1983).

[86] EPTA Collaboration and InPTA Collaboration, Antoniadis, J., Arumugam, P., Arumugam, S., Babak, S., Bagchi, M., Bak Nielsen, A.-S., Bassa, C. G., Bathula, A., Berthereau, A., Bonetti, M., Bortolas, E., Brook, P. R., Burgay, M., Caballero, R. N., Chalumeau, A., Champion, D. J., Chanlaridis, S., Chen, S., Cognard, I., Dandapat, S., Deb, D., Desai, S., Desvignes, G., Dhanda-Batra, N., Dwivedi, C., Falxa, M., Ferdman, R. D., Franchini, A., Gair, J. R., Goncharov, B., Gopakumar, A., Graikou, E., Grießmeier, J.-M., Guillemot, L., Guo, Y. J., Gupta, Y., Hisano, S., Hu, H., Iraci, F., Izquierdo-Villalba, D., Jang, J., Jawor, J., Janssen, G. H., Jessner, A., Joshi, B. C., Kareem, F., Karuppusamy, R., Keane, E. F., Keith, M. J., Kharbanda, D., Kikunaga, T., Kolhe, N., Kramer, M., Krishnakumar, M. A., Lackeos, K., Lee, K. J., Liu, K., Liu, Y., Lyne, A. G., McKee, J. W., Maan, Y., Main, R. A., Mickaliger, M. B., Nițu, I. C., Nobleson, K., Paladi, A. K., Parthasarathy, A., Perera, B. B. P., Perrodin, D., Petiteau, A., Porayko, N. K., Possenti, A., Prabu, T., Quelquejay Leclere, H., Rana, P., Samajdar, A., Sanidas, S. A., Sesana, A., Shaifullah, G., Singha, J., Speri, L., Spiewak, R., Srivastava, A., Stappers, B. W., Surnis, M., Susarla, S. C., Susobhanan, A., Takahashi, K., Tarafdar, P., Theureau, G., Tiburzi, C., van der Wateren, E., Vecchio, A., Venkatraman Krishnan, V., Verbiest, J. P. W., Wang, J., Wang, L., and Wu, Z., "The second data release from the european pulsar timing array - iii. search for gravitational wave signals," A&A 678, A50 (2023).

[87] J. Antoniadis, P. Arumugam, S. Arumugam, P. Auclair, S. Babak, M. Bagchi, A. S. Bak Nielsen, E. Barausse, C. G. Bassa, A. Bathula, A. Berthereau, M. Bonetti, E. Bortolas, P. R. Brook, M. Burgay, R. N. Caballero, C. Caprini, A. Chalumeau, D. J. Champion, S. Chanlaridis, S. Chen, I. Cognard, M. Crisostomi, S. Dandapat, D. Deb, S. Desai, G. Desvignes, N. Dhanda-Batra, C. Dwivedi, M. Falxa, F. Fastidio, R. D. Ferdman, A. Franchini, J. R. Gair, B. Goncharov, A. Gopakumar, E. Graikou, J. M. Grießmeier, A. Gualandris, L. Guillemot, Y. J. Guo, Y. Gupta, S. Hisano, H. Hu, F. Iraci, D. Izquierdo-Villalba, J. Jang, J. Jawor, G. H. Janssen, A. Jessner, B. C. Joshi, F. Kareem, R. Karuppusamy, E. F. Keane, M. J. Keith, D. Kharbanda, T. Khizriev, T. Kikunaga, N. Kolhe, M. Kramer, M. A. Krishnakumar, K. Lackeos, K. J. Lee, K. Liu, Y. Liu, A. G. Lyne, J. W. McKee, Y. Maan, R. A. Main, M. B. Mickaliger, H. Middleton, A. Neronov, I. C. Nitu, K. Nobleson, A. K. Paladi, A. Parthasarathy, B. B. P. Perera, D. Perrodin, A. Petiteau, N. K. Porayko, A. Possenti, T. Prabu, K. Postnov, H. Quelquejay Leclere, P. Rana, A. Roper Pol, A. Samajdar, S. A. Sanidas, D. Semikoz, A. Sesana, G. Shaifullah, J. Singha, C. Smarra, L. Speri, R. Spiewak, A. Srivastava, B. W. Stappers, D. A. Steer, M. Surnis, S. C. Susarla, A. Susobhanan, K. Takahashi, P. Tarafdar, G. Theureau, C. Tiburzi, R. J. Truant, E. van der Wateren, S. Valtolina, A. Vecchio, V. Venkatraman Krishnan, J. P. W. Verbiest, J. Wang, L. Wang, and Z. Wu, "The second data release from the european pulsar timing array: V. implications for massive black holes, dark matter and the early universe," (2023), arXiv:2306.16227 [astro-ph.CO].

[88] Zaven Arzoumanian, Paul T. Baker, Harsha Blumer, Bence Bécsy, Adam Brazier, Paul R. Brook, Sarah Burke-Spolaor, Shami Chatterjee, Siyuan Chen, James M. Cordes, Neil J. Cornish, Fronefield Crawford, H. Thankful Cromartie, Megan E. Decesar, Paul B. Demorest, Timothy Dolch, Justin A. Ellis, Elizabeth C. Ferrara, William Fiore, Emmanuel Fonseca, Nathan Garver-Daniels, Peter A. Gentile, Deborah C. Good, Jeffrey S. Hazboun, A. Miguel Holgado, Kristina Islo, Ross J. Jennings, Megan L. Jones, Andrew R. Kaiser, David L. Kaplan, Luke Zoltan Kelley, Joey Shapiro Key, Nima Laal, Michael T. Lam, T. Joseph W. Lazio, Duncan R. Lorimer, Jing Luo, Ryan S. Lynch, Dustin R. Madison, Maura A.
McLaughlin, Chiara M. F. Mingarelli, Cherry Ng, David J. Nice, Timothy T. Pennucci, Nihan S. Pol, Scott M. Ransom, Paul S. Ray, Brent J. Shapiro-Albert, Xavier Siemens, Joseph Simon, Renée Spiewak, Ingrid H. Stairs, Daniel R. Stinebring, Kevin Stovall, Jerry P. Sun, Joseph K. Swiggum, Stephen R. Taylor, Jacob E. Turner, Michele Vallisneri, Sarah J. Vigeland, Caitlin A. Witt, and Nanogray Collaboration, "The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background," ApJ 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE].

[89] Boris Goncharov, R. M. Shannon, D. J. Reardon, G. Hobbs, A. Zic, M. Bailes, M. Curylo, S. Dai, M. Kerr, M. E. Lower, R. N. Manchester, R. Mandow, H. Middleton, M. T. Miles, A. Parthasarathy, E. Thrane, N. Thyagarajan, X. Xue, X. J. Zhu, A. D. Cameron, Y. Feng, R. Luo, C. J. Russell, J. Sarkissian, R. Spiewak, S. Wang, J. B. Wang, L. Zhang, and S. Zhang, "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, L19 (2021), arXiv:2107.12112 [astro-ph.HE].

[90] S. Chen, R. N. Caballero, Y. J. Guo, A. Chalumeau, K. Liu, G. Shaifullah, K. J. Lee, S. Babak, G. Desvignes, A. Parthasarathy, H. Hu, E. van der Wateren, J. Antoniadis, A. S. Bak Nielsen, C. G. Bassa, A. Berthereau, M. Burgay, D. J. Champion, I. Cognard, M. Falxa, R. D. Ferdman, P. C. C. Freire, J. R. Gair, E. Graikou, L. Guillemot, J. Jang, G. H. Janssen, R. Karuppusamy, M. J. Keith, M. Kramer, X. J. Liu, A. G. Lyne, R. A. Main, J. W. McKee, M. B. Mickaliger, B. B. P. Perera, D. Perrodin, A. Petiteau, N. K. Porayko, A. Possenti, A. Samajdar, S. A. Sanidas, A. Sesana, L. Speri, B. W. Stappers, G. Theureau, C. Tiburzi, A. Vecchio, J. P. W. Verbiest, J. Wang, L. Wang, and H. Xu, "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, 4970-4993 (2021), arXiv:2110.13184 [astro-ph.HE].

[91] Rutger van Haasteren, Yuri Levin, Patrick McDonald, and Tingting Lu, "On measuring the gravitational-wave background using Pulsar Timing Arrays," Monthly Notices of the Royal Astronomical Society 395, 1005-1014 (2009), arXiv:0809.0791 [astro-ph].

[92] L. Lentati, P. Alexander, M. P. Hobson, F. Feroz, R. van Haasteren, K. J. Lee, and R. M. Shannon, "TEMPONEST: a Bayesian approach to pulsar timing analysis," Monthly Notices of the Royal Astronomical Society 437, 3004-3023 (2014), arXiv:1310.2120 [astro-ph.IM].

[93] S. R. Taylor, L. Lentati, S. Babak, P. Brem, J. R. Gair, A. Sesana, and A. Vecchio, "All correlations must die: Assessing the significance of a stochastic gravitational-wave background in pulsar timing arrays," Phys. Rev. D 95, 042002 (2017), arXiv:1606.09180 [astro-ph.IM].

[94] Justin A. Ellis, Michele Vallisneri, Stephen R. Taylor, and Paul T. Baker, "Enterprise: Enhanced numerical toolbox enabling a robust pulsar inference suite," Zenodo (2020).

[95] Stephen R. Taylor, Paul T. Baker, Jeffrey S. Hazboun, Joseph Simon, and Sarah J. Vigeland, "enterprise_extensions," (2021), v2.3.3.

[96] Justin Ellis and Rutger van Haasteren, "jellis18/ptmcmc sampler: Official release," (2017).

[97] Andrea Mitridate, "Ptarcade," (2023), 10.5281/zenodo.7876430.

[98] Andrea Mitridate, David Wright, Richard von Eckardstein, Tobias Schröder, Jonathan Nay, Ken Olum, Kai Schmitz, and Tanner Trickle, "PTArcade," (2023), arXiv:2306.16377 [hep-ph].

[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, 39 (2012).

[100] Jai chan Hwang, Donghui Jeong, Hyerim Noh, and Clemente Smarra, "Pulsar timing array signature from oscillating metric perturbations due to ultra-light axion," Journal of Cosmology and Astroparticle Physics 2024, 014 (2024).

[101] Jeffrey S. Hazboun, Joseph D. Romano, and Tristan L. Smith, "Realistic sensitivity curves for pulsar timing arrays," Phys. Rev. D 100, 104028 (2019), arXiv:1907.04341 [gr-qc].

[102] R. Blandford, R. Narayan, and R. W. Romani, "Arrival-time analysis for a millisecond pulsar," Journal of Astrophysics and Astronomy 5, 369-388 (1984).

[103] Harikrishnan Ramani, Tanner Trickle, and Kathryn M. Zurek, "Observability of dark matter substructure with pulsar timing correlations," Journal of Cosmology and Astroparticle Physics 2020, 033 (2020).

[104] Caner Unal, Federico R. Urban, and Ely D. Kovetz, "Probing ultralight scalar, vector and tensor dark matter with pulsar timing arrays," (2022), arXiv:2209.02741 [astro-ph.CO].

[105] Mark Abney, Richard I. Epstein, and Angela V. Olinto, "Observational constraints on the internal structure and dynamics of the vela pulsar," The Astrophysical Journal 466, L91 (1996).

[106] Raissa F. P. Mendes, "Possibility of setting a new constraint to scalar-tensor theories," Phys. Rev. D 91, 064024 (2015), arXiv:1412.6789 [gr-qc].

[107] David Anderson and Nicolás Yunes, "Solar System constraints on massless scalar-tensor gravity with positive coupling constant upon cosmological evolution of the scalar field," Phys. Rev. D 96, 064037 (2017), arXiv:1705.06351 [gr-qc].

[108] Gemma Janssen et al., "Gravitational wave astronomy with the SKA," PoS AASKA14, 037 (2015), arXiv:1501.00127 [astro-ph.IM].

[109] Hyungjin Kim, "Gravitational interaction of ultralight dark matter with interferometers," JCAP 12, 018 (2023), arXiv:2306.13348 [hep-ph].

[110] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, "Equation of state of nucleon matter and neutron star structure," Phys. Rev. C 58, 1804-1828 (1998).

[111] H. Müther, M. Prakash, and T.L. Ainsworth, "The nuclear symmetry energy in relativistic brueckner-hartree-fock calculations," Physics Letters B 199, 469-474 (1987).

[112] F. Douchin and P. Haensel, "A unified equation of state of dense matter and neutron star structure," Astronomy & Astrophysics 380, 151-167 (2001).
[113] M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K. Harding, W. C. G. Ho, J. M. Lattimer, R. M. Ludlam, S. Mahmoodifar, S. M. Morsink, P. S. Ray, T. E. Strohmayer, K. S. Wood, T. Enoto, R. Foster, T. Okajima, G. Prigozhin, and Y. Soong, “Psr j0030+0451 mass and radius from nicer data and implications for the properties of neutron star matter,” The Astrophysical Journal Letters **887**, L24 (2019).

[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**, 4031–4042 (2020), 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**, A128 (2011).

[116] J. Antoniadis, M. H. van Kerkwijk, D. Koester, P. C. C. Freire, N. Wex, T. M. Tauris, M. Kramer, and C. G. Bassa, “The relativistic pulsar–white dwarf binary PSR J1738+0333 – I. Mass determination and evolutionary history,” Monthly Notices of the Royal Astronomical Society **423**, 3316–3327 (2012), 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**, L113 (2005).