# Bridging the micro-Hz gravitational wave gap via Doppler tracking with the
Uranus Orbiter and Probe Mission:
Massive black hole binaries, early universe signals and ultra-light dark
matter

Lorenz Zwick1<EMAIL_ADDRESS>Deniz Soyuer2 Daniel J. D’Orazio1 David
O’Neill1 Andrea Derdzinski3,4 Prasenjit Saha5 Diego Blas6,7 Alexander C.
Jenkins8 Luke Zoltan Kelley9 1Niels Bohr International Academy, The Niels Bohr
Institute, Blegdamsvej 17, DK-2100, Copenhagen, Denmark
2Department of Astrophysics, University of Zurich, Winterthurerstrasse 190,
8057 Zurich, Switzerland
3Department of Life and Physical Sciences, Fisk University, 1000 17th Avenue
N., Nashville, TN 37208, USA
4Department of Physics & Astronomy, Vanderbilt University, 2301 Vanderbilt
Place, Nashville, TN 37235, USA
5Physik-Institut, Universität Zürich, Winterthurerstrasse 190, CH-8057 Zürich,
Switzerland
6Institut de Fisica d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain
7Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluís
Companys 23, 08010 Barcelona, Spain
8Department of Physics and Astronomy, University College London, London WC1E
6BT, UK
9Department of Astronomy, University of California, Berkeley, California
94720, USA

###### Abstract

With the recent announcement by NASA’s Planetary Science and Astrobiology
Decadal Survey 2023-2032, a priority flagship mission to the planet Uranus is
anticipated. Here, we explore the prospects of using the mission’s radio
Doppler tracking equipment to detect gravitational waves (GWs) and other
analogous signals related to dark matter (DM) over the duration of its
interplanetary cruise. By employing a methodology to stack tracking data in
combination with Monte-Carlo Markov-Chain parameter recovery tests, we show
that the mission will be sensitive to GWs over the wide frequency range of
$3\times 10^{-9}$ Hz to $10^{-1}$ Hz, provided that tracking data is taken
consistently over a large fraction of the cruise duration. Thus, the mission
has the potential to fill the gap between pulsar timing and space-based-
interferometry GW observatories. Within this assumption, we forecast the
detection of $\mathcal{\mathcal{O}}(1-100)$ individual massive black hole
binaries using two independent population models. Additionally, we determine
the mission’s sensitivity to both astrophysical and primordial stochastic
gravitational wave backgrounds, as well as its capacity to test, or even
confirm via detection, ultralight DM models. In all these cases, the tracking
of the spacecraft over its interplanetary cruise would enable coverage of
unexplored regions of parameter space, where signals from new phenomena in our
Universe may be lurking.

##  I Introduction

Almost half a century after the launch of the Voyager 2 space probe, the
prospect of a new visit to the ice giant Uranus is finally crystallizing
again: A flagship mission to the planet Uranus (hereafter referred to as
“UOP”, standing for Uranus Orbiter and Probe) has been declared a priority,
according to the announcement by NASA’s Planetary Science and Astrobiology
Decadal Survey 2023-2032.111Survey provided in: Origins, Worlds, and Life: A
Decadal Strategy for Planetary Science and Astrobiology 2023-2032 (2022) The
announcement is timely; numerous publications have underlined the rich
potential of such a mission in terms of planetary science over the past few
years [54, 42, 43, 51, 59, 93, 63, 49, 58]. Yet, any mission to the outer
Solar system must undergo a long cruise in interplanetary space before
reaching its destination. According to a notional, short duration trajectory
of the UOP (see Fig. 1), just the transfer between Jupiter and Uranus is
expected to require around 9 years, during which the planetary science yield
is virtually non-existent.

Several works have highlighted how the UOP’s long cruise may instead provide
unique opportunities for non-Uranian science. Prospects include the possible
measurement of the local dark matter content in the Solar system [89, 109], as
well as the localisation of the hypothetical Planet 9 on the sky [24]. The
opportunity to take such measurements would be afforded by the presence of a
radio link between the spacecraft and the Earth, which would allow one to
closely monitor changes in the spacecraft’s velocity and reconstruct
deviations from the expected mission trajectory. Among all, the most timely
application of interplanetary radio tracking is arguably the detection of
gravitational wave (GW) signals. The idea of detecting GWs in the Doppler
tracking data of planetary missions has a rich history and has been attempted
previously with Pioneer 11 [8], the Galileo–Ulysses–Mars Observer coincidence
experiment [5, 14], and Cassini [30, 13, 9], albeit without any successful
detection candidates due to insufficient sensitivity. As first highlighted in
Ref. [97], the prospective UOP offers an exceptional opportunity compared to
the previous attempts due to a longer spacecraft cruise, a longer detector
arm, and potentially improved radio technology.222Recently, there was another
proposal to use tracking of satellites to Mars to study GWs of frequencies
from $10^{-4}$ Hz to $10^{-1}$ Hz [15]. As compared to our proposal, these
constraints are too weak to be relevant for known sources, while they will be
largely superseded in the future by LISA [29].

In this work, we consider how the collection of tracking data over the entire
interplanetary cruise can increase the sensitivity of the mission to GWs from
many sources, allowing us to probe the yet poorly covered micro-Hz regime.
Additionally, we outline the mission’s capability to greatly expand the range
of ultra-light dark matter (ULDM) models [57, 21] that can be directly
detected, including the possibility of a first detection of DM in the Solar
System. The key for these results is a method developed here to stack
individual short-duration tracking runs, showing how the mission’s GW
frequency sensitivity band and DM mass sensitivity are in principle only
limited by the total observation time and the data cadence, in analogy to the
Pulsar Timing Array analysis for nano-Hz GW detection [hereafter PTA, see,
e.g., 46, 90, 53, 6, 2]. For the UOP, this places the frequency limits to be
between $3\times 10^{-9}$ Hz ($\sim 1/(10$ yr); the interplanetary cruise
time) and $10^{-1}$ Hz ($1/(10$ s); the resolution time), respectively, and
the DM mass limits in the range $10^{-23}$ eV/$c^{2}$ to $10^{-15}$
eV/$c^{2}$.

In summary, we show that, beyond its Uranian science goals, the UOP could
double as a sensitive GW detector in the micro-Hz band, probing the
astrophysics of massive black holes or signals from the early Universe, and as
a dark matter detector for models not accessible today by any other probe.
This should be achievable without any sacrifice to planetary science goals by
tracking the mission consistently over its interplanetary cruise. We also
remark that the mission is expected to fly in concomitance with the Laser
Interferometer Space Antenna (LISA [29]). The opportunities for
complementarity between these two missions were briefly discussed in Ref.
[97]. They deserve a more thorough analysis, which we leave for future work.

This paper is structured as follows: In Section II we introduce the basics of
GW and ULDM detection via Doppler tracking and detail our strategy to track
signals over the entire mission duration. In Section III, we construct GW
sensitivity curves and define three possible mission configurations with
increasing levels of technological development. In Section IV we estimate the
number of detections of individual massive black hole binaries with two
independent population models, as well as the associated stochastic
gravitational wave background. Section V and Section VI include our forecast
constraints on GW signals from the early universe, as well as various DM
models, respectively. Lastly, we discuss the prospects to further improve the
sensitivity of the mission in Section VII and present our concluding remarks
in Section VIII.

Figure 1: One of the many proposed mission plans of a prospective Uranus
mission that involves a Jupiter Gravity Assist and a subsequent cruise to
Uranus. This specific mission timeline corresponds to the shortest proposed
cruise time and is available in the public mission document.

##  II Doppler Tracking of GWs and other signals

###  II.1 Response of Doppler tracking to a passing GW

The passage of a GW between the Earth and the spacecraft modifies the
propagation time of photons, inducing small timing shifts in the tracking
system. These shifts are best described by defining the dimensionless
fractional frequency fluctuation of a two-way Doppler system:

$y_{2}(t)={\Delta\nu}/{\nu_{0}},$ (1)

where $\nu_{0}$ is the tracking system’s carrier frequency, most often in the
Ka-band ($\sim$30 GHz) or the X-band ($\sim$8 GHz), while $\Delta\nu$ is the
measured Doppler shift with respect to the carrier frequency. Following Ref.
[7], we can write down an explicit formula for the frequency fluctuation due
to a passing GW:

$\displaystyle
y_{2}^{\scriptscriptstyle\rm{GW}}(t)=\frac{\mu-1}{2}\bar{\Psi}(t)~{}-~{}$
$\displaystyle\mu\bar{\Psi}\left(t-\frac{\mu+1}{2}T_{2}\right)$
$\displaystyle~{}+~{}$ $\displaystyle\frac{\mu+1}{2}\bar{\Psi}(t-T_{2}),$ (2)

where $T_{2}$ is the two-way light travel time between the Earth and the
spacecraft, $\mu=\hat{\bm{k}}\cdot\hat{\bm{n}}$ is the projection of the unit
wave vector $\hat{\bm{k}}$ of the GW onto the unit vector connecting the Earth
and the spacecraft $\hat{\bm{n}}$, and $\bar{\Psi}$ is the projection of the
GW amplitude onto the Doppler link [105]. The latter is given by
$\bar{\Psi}(t)=(\hat{\bm{n}}\cdot\mathrm{\bf{h}}(t)\cdot\hat{\bm{n}})/(1-\mu^{2})$,
where
$\mathrm{\bf{h}}(t)=\mathrm{{h}}_{+}(t)\,\mathrm{\bf{e}}_{+}+\,\mathrm{{h}}_{\times}(t)\,\mathrm{\bf{e}}_{\times}$
is the GW amplitude (i.e., the strain) and $\mathrm{\bf{e}}_{+,\times}$ are
the usual “plus” and “cross” polarization states of a transverse, traceless
plane GW [82, 22]. While Eq. (2) can be used to describe any incident GW, in
this Section we will focus on coherent, monochromatic GWs with the form:

$\displaystyle\mathrm{h}_{+}(t)=A\frac{1+\cos^{2}\iota}{2}\cos(f_{\scriptscriptstyle\rm
GW}t+\phi_{0}),$ (3a)
$\displaystyle\mathrm{h}_{\times}(t)=A\cos\iota\sin(f_{\scriptscriptstyle\rm
GW}t+\phi_{0})\,,$ (3b)

where $f_{\scriptscriptstyle\rm GW}$ is the frequency of the GW, $A$ its
amplitude, $\iota$ is the binary’s inclination, and $\phi_{0}$ is an initial
phase [see, e.g., 72, 74, 69]. In particular, Eqs. (3) can be used to describe
signals produced by supermassive BH binaries in the early inspiral stages, the
most promising astrophysical source of GWs in the nano-Hz to micro-Hz
frequency band, and will be discussed in detail in Section IV.

###  II.2 Response of Doppler tracking to a SGWB

Following [9], the results from the previous section can be immediately used
to derive a bound on stochastic GW backgrounds (SGWBs). Recall that many
sources of GWs produce a signal that is comprised of several independent waves
that interfere incoherently. These are the kind of signals expected from early
Universe physics [25], and many astrophysical sources (see, e.g. [91] for
examples in the micro-Hz band). SGWBs are modelled by considering the power
spectrum of the sky- and polarization-average of the signal given in Eq. (2):
$S_{y_{2}}(f)$. Following [37], this quantity can be expressed as
$S_{y_{2}}(f)=\bar{R}_{2}(f)S_{\scriptscriptstyle\rm GW}(f)$, where
$S_{\scriptscriptstyle\rm GW}(f)$ is the GW power spectrum, and the transfer
function $\bar{R}_{2}(f)$ can be found in [9]. In the limit of low frequencies
($2\pi fT_{2}\ll 1$), it reads simply $\bar{R}_{2}(f)\approx(8/15)\left(\pi
fT_{2}\right)^{2}$, while at higher frequencies it displays some oscillatory
behaviour. In order to investigate the sensitivity of Doppler tracking to
SGWBs, it is customary to express the constraints in terms of the normalised
energy density of the SGWB for logarithmic frequency unit:

$\Omega_{\scriptscriptstyle\rm GW}(f)\equiv\frac{8\pi
G}{3H_{0}^{2}}\frac{\mathrm{d}\rho_{\scriptscriptstyle\rm GW}}{\mathrm{d}\log
f}=\frac{8\pi^{2}f^{3}S_{y_{2}}(f)}{3\bar{R}_{2}(f)H_{0}^{2}},$ (4)

where we will consider the Hubble parameter to be $H_{0}\approx 67\,$km/s/Mpc
[4]. In practice, we will impose that $\Omega_{\scriptscriptstyle\rm
GW}(f)<\frac{8\pi^{2}f^{3}S_{y_{2}}(f)}{3\bar{R}_{2}(f)H_{0}^{2}}$, where
$S_{y_{2}}(f)$ represents the power spectrum of the noise of the experiment.

###  II.3 Response of Doppler tracking to ultra-light dark matter

Dark matter (DM) particles of masses $m_{\scriptscriptstyle\rm DM}\lesssim 1$
eV/$c^{2}$ are distributed in the Milky Way at distances smaller than their de
Broglie wavelength (associated with the ‘size’ of the particle). As a result,
DM can be approximated in this range of masses as a classical non-relativistic
massive field $\phi_{\scriptscriptstyle\rm DM}$ [57]. From the virialised
properties of DM in the Milky Way, this field will be distributed in coherent
patches of characteristic length

$l\approx 0.9\times
10^{6}\mathrm{~{}AU}\left(\frac{10^{-3}}{\sigma_{0}}\right)\left(\frac{10^{-20}\mathrm{eV}/c^{2}}{m_{\scriptscriptstyle\rm
DM}}\right),$ (5)

where $\sigma_{0}$ is the DM velocity dispersion, generated by virialization
in the Milky Way. In these regions, the DM field can be approximated as a
monochromatic function oscillating at a frequency

$f_{\scriptscriptstyle\rm DM}\approx 2.4\times
10^{-6}\left(\frac{m_{\scriptscriptstyle\rm
DM}}{10^{-20}\mathrm{eV}/c^{2}}\right)\mathrm{Hz},$ (6)

and coherent over $c^{2}/\sigma_{0}^{2}$ oscillations. As a result, the
gravitational potentials generated by the galactic DM (derived from solving
the Poisson equation, which is quadratic in $\phi_{\scriptscriptstyle\rm DM}$)
inherit the coherent features of the DM field, and oscillate with frequency
$f_{\rm grav}\approx 2f_{\scriptscriptstyle\rm DM}$ in these patches (cf. Eq.
(8) below). When approaching the lowest viable masses, these DM models are
known as fuzzy dark matter [56] or ultra-light dark matter (ULDM) [57].

The Doppler tracking suggested above for GWs can be easily adapted to search
for ULDM in different ways. First, notice that, at leading order, the
fluctuations of the gravitational potential
$g_{ij}\equiv-(1-2\psi)\delta_{ij}$ modify the relative frequency of the radio
signal as [62, 21, 83, 101]

$\Delta\nu/\nu_{0}=\psi(t)-\psi\left(t-T_{2}\right).$ (7)

For ULDM, the fluctuating piece of the gravitational potential $\psi$ behaves
(at leading order in $\sigma_{0}\approx 10^{-3}$) as

$\begin{split}\psi(t)\approx\psi_{g}(t)&\equiv\,6\cdot
10^{-16}\left(\frac{\rho_{\scriptscriptstyle\rm
DM}}{0.3\,\mathrm{GeV}/\mathrm{cm}^{3}}\right)\\\
&\left(\frac{10^{-23}\,\mathrm{eV}/c^{2}}{m_{\scriptscriptstyle\rm
DM}}\right)^{2}\cos(4\pi f_{\scriptscriptstyle\rm
DM}t+2\alpha_{0}),\end{split}$ (8)

where $\rho_{\scriptscriptstyle\rm DM}$ is the dark matter energy density and
$\alpha_{0}$ is a random phase.333Technically, the interference of the
different classical waves that comprise the ULDM implies that the local
density has a degree of randomness over this value, which corresponds to the
average DM energy density [45]. We will assume that the DM energy density in
the Solar System is 1/2 of the averaged one (expected to be of order
$\rho_{\rm DM}\approx 0.3\,\rm{GeV/cm^{3}}$) to provide our estimates. This
immediately allows one to use the bounds on
$y_{2}^{\scriptscriptstyle\mathrm{GW}}(t)$ to constrain the mass of dark
matter candidates (it is enough to use (7) in (1), and connect it to (2)). To
find the precise bound, it is important to realise that from the lack of
polarization and directional dependence, a factor of $\sqrt{15}$ appears when
comparing the bounds on rms of $h(t)$ vs $\psi(t)$ [21].

Second, it is normally expected that DM is not completely ‘dark’, and couples,
though very weakly, to matter or light. As a result, the fluctuations of
$\phi_{\rm DM}$ can affect the tracking data in new ways. For instance, new
forces may affect the trajectory of gravitating probes or the properties of
the tracking beam, see, e.g., [71, 1, 40, 47, 81, 77]. For our current study,
we will focus on two simple, though representative possibilities of _direct
coupling_ of DM to ordinary matter and light. In the first model, dark matter
is coupled to the constituents of the detectors, the spacecraft, and the
tracking beam through a _universal coupling_ , common to all particles of the
standard model of particle physics. This universality is similar to the
coupling of the metric $g_{\mu\nu}$ in general relativity. As a consequence,
in this model the coupling of matter and light to gravity ($g_{\mu\nu}$) and
dark matter ($\phi_{\rm DM}$) happens through the interaction with an
effective metric
$\bar{g}_{\mu\nu}=g_{\mu\nu}(1+2\alpha(\phi_{\scriptscriptstyle\rm DM}))$. In
this situation, matter and light will move in geodesics of $\bar{g}_{\mu\nu}$,
implying that the expression (7) will still be valid, identifying444More
specifically, at leading order in perturbations over the flat background,
$\bar{g}_{ij}=-(1-2\psi-2\alpha(\phi_{\scriptscriptstyle\rm DM}))\delta_{ij}$.
$\psi\equiv\psi_{g}(t)+\alpha(\phi_{\scriptscriptstyle\rm DM})$. The function
$\alpha(\phi)$ is model-dependent. To illustrate the potential of UOP to test
DM models, we will choose $\alpha(\phi_{\scriptscriptstyle\rm
DM})=\phi^{2}_{\scriptscriptstyle\rm DM}/\Lambda^{2}_{2}$, where $\Lambda_{2}$
parameterises the intensity of the coupling. The effect of ULDM is described
in this case through (7), after the identification

$\begin{split}\psi(t)\approx\left(1+40\left(\frac{10^{18}\rm{GeV}}{\Lambda_{2}}\right)^{2}\right)\,\psi_{g}(t).\end{split}$
(9)

To illustrate other possibilities opened from the coupling of the tracking
(light) beam to DM, we will consider a model where DM is an axion-light
particle [28]. In this case, a coupling of the form

$g_{a\gamma\gamma}\phi_{\rm DM}\vec{E}\cdot\vec{B},$ (10)

is expected, where $\vec{E}$ and $\vec{B}$ are the electric and magnetic
fields respectively, and $g_{a\gamma\gamma}$ represents the intensity of the
coupling. Following [17], for the ULDM case one expects frequency oscillations
in the tracking beam of order

$\frac{\Delta\nu}{\nu_{0}}\sim
10^{-16}\left(\frac{g_{a\gamma\gamma}}{10^{-10}\mathrm{GeV}^{-1}}\right)\left(\frac{\mathrm{GHz}}{\nu/2\pi}\right)\sqrt{\frac{\rho_{\scriptscriptstyle\rm{DM}}}{0.3\mathrm{GeV}/\mathrm{cm}^{3}}}$
(11)

at frequency $f_{\scriptscriptstyle\rm DM}$, where we have normalised
$g_{a\gamma\gamma}$ to a value representative of current sensitivity [80, 50].

The previous three effects, summarised by Eqs. (8), (9) and (11) are far from
being comprehensive (for instance, we are not discussing the possibility of
other possibilities for $\alpha(\phi_{\rm DM})$, beyond the quadratic case of
(9), or models where DM is a vector field, as in the case of dark photons).
Still, they represent three families of ULDM models: universal gravitational
coupling; direct coupling to matter; and direct coupling to light, which
illustrate the potential of UOP to explore ULDM models.

###  II.4 Noise in a Doppler tracking system

The response of a Doppler tracking system to GWs and other oscillatory signals
must be compared to its intrinsic noise. Noise sources range from mechanical
(caused, e.g., by flexibility in the antennae), to astrophysical (e.g.,
scintillation due to propagation through the interplanetary plasma). These are
thoroughly compiled and discussed in the extensive review by Ref. [7], and
further discussed in Section VII. In general, the total noise profile can be
described in either the time domain or the frequency domain. In the former,
the size of typical noise fluctuations for every data-point is characterised
by the system’s Allan deviation, $\sigma_{\rm A}$, which is a function of the
chosen averaging time $\tau$. In the latter, noise is best described by its
one-sided power spectrum $S_{n}$, which is a function of frequency $f$. These
two crucial quantities are related by the following equation:

$\displaystyle\sigma_{\rm
A}(\tau)=\sqrt{4\int_{0}^{\infty}S_{n}(f)\frac{\sin(f\tau)^{4}}{(\pi
f\tau)^{2}}\,{\rm d}f}\approx\sqrt{\frac{S_{n}(1/\tau)}{\tau}},$ (12)

where the latter approximation can be verified numerically, and is exact in
the case of white noise [13]. Eq. (12) essentially states that the typical
noise fluctuations given an averaging time $\tau$ are determined by the noise
power spectral density at the reciprocal frequency $f\sim 1/\tau$.

The dependence of the Allan deviation on the averaging time $\tau$ is crucial
when attempting to reconstruct a signal: On the one hand, a longer averaging
time generally leads to reduced noise and improves the quality of each data
point. On the other hand, the signal-to-noise ratio (SNR) also scales as the
square root of the total amount of collected data points $N_{\rm{s}}(\tau)$
which is inherently limited by the minimum sampling rate $1/\tau$ ($N\sim
1/\tau$). These competing effects define a sweet spot for the averaging time,
which compromises between the quantity and quality of tracking data. It can be
found by minimising the function:

$\displaystyle\frac{\sigma_{\rm
A}(\tau)}{\sqrt{N_{\rm{s}}(\tau)}}\sim\sqrt{S_{n}(1/\tau)},$ (13)

where we neglect an arbitrary normalisation given by the total observation
time.

A simple parameterisation of a noise power spectral density is shown in Fig.
2, based on theoretical work by Ref. [13] and extensive analysis of the
Ulysses and Cassini tracking data [14, 30]. The broad features of the power
spectral density can be captured by a simple three power law structure. At
frequencies higher than $f\sim 0.1$ Hz, the noise is dominated by thermal
fluctuations and steeply increases as $\sim f^{2}$. This limit tends to
constrain the minimum averaging time of tracking systems to $\tau\sim 10$ s,
beyond which the Allan deviation rapidly deteriorates. A central region of
optimal sensitivity spans between $\sim 10^{-4}$ Hz and $\sim 10^{-1}$ Hz and
is characterised by a red spectrum $\sim f^{-1/2}$ [30]. The upturn in noise
below $10^{-4}$ Hz is an attempt to model less constrained low-frequency noise
sources, which range from uncertainties in the dynamical modelling of the
spacecraft orbits to the effects of atmospheric turbulence. They are
conservatively estimated to generate a slope of $f^{-2}$, as has been
originally observed in the Ulysses tracking data. However, it is important to
note that these estimates are purely phenomenological. The effects of
achieving a flatter low-frequency slope for the noise, and the resulting
improvements in sensitivity, have been explored in [30]. Ultimately, the
behaviour of the power spectral density below frequencies of $\sim 10^{-4}$ Hz
cannot be studied straightforwardly owing to the rotation of the Earth: The
spacecraft will fall below the horizon of any given tracking station within a
fraction of a day, severing the Doppler link and preventing the continuous
tracking of a coherent signal.

To circumvent this complication we consider a typical tracking run of a single
station to have a maximal duration of approximately $8$ hours, after which the
link may only be reinstated after a clock reset. This amounts to either
switching tracking stations or waiting for the next day.

Figure 2: The blue curve shows the noise power spectral density based on
theoretical work by Ref. [14] and extensive analysis of the Ulysses and
Cassini tracking data [30]. The gray dashed curve is a visualisation of the
trade-off between the total amount of data-points $N_{\rm d}$ and the quality
of each data-point for the purposes of reconstructing the trajectory, where
the absolute values for the latter curve have no significance. The curve is
proportional to $\sqrt{fS_{\rm n}}\times\sqrt{N_{\rm d}}$, where $N_{\rm
d}\sim 1/\tau$ and $\tau$ is the averaging time of the Doppler link. The
lowest point on the gray dashed curve therefore represents the sampling rate
at which optimal SNR is achieved.

For the purposes of this work, we will conservatively assume that the noise of
the Cassini Doppler tracking system is representative of modern interplanetary
spacecraft within the relevant frequency range. We also discuss some more
advanced technological prospects in Section VII.2. With this choice, we are
assuming that noise components at lower frequencies do not bleed into the
sampled range of the power spectral density. Such components take the form of
unmodelled accelerations due to small massive bodies and/or variations in the
solar radiation pressure. Typically, they result in a subdominant contribution
to the Allan deviation at $10\,$s of $\sim 10^{-16}$ [7].

As mentioned above, the noise power spectral density is specified by two
crucial parameters, the optimal sampling rate and its corresponding Allan
deviation. Given the Cassini mission specifications, the two parameters take
the values:

$\displaystyle\tau^{\rm{Cass}}$ $\displaystyle\sim 10\,{\rm
s},~{}~{}~{}~{}\sigma_{\rm A}^{\rm{Cass}}(10\,{\rm s})\sim 10^{-14},$ (14)

which are found by looking for the minimum of Eq. (13). For the purposes of
this work, we fix the optimal sampling rate at 10 s, but allow the Allan
deviation to vary slightly from its Cassini-era value in order to reflect
technological advancements in the last four decades. The effects of these
choices will be further discussed in Section VII.

###  II.5 Strategy for the tracking of low-frequency GWs

Measuring a low-frequency gravitational wave requires a long observation time,
at least of the order $T\sim 1/f_{\rm{GW}}$. Here we devise a strategy to
integrate gravitational wave signals over the total mission duration of
$\sim$10 years. Previous attempts to estimate the sensitivity of Doppler
tracking to such GWs have typically assumed that the noise power spectral
density may be extended to lower frequencies following the phenomenological
slope of $\sim$$f^{-2}$ [13]. Here we take a different approach, which
consists of stacking many independent tracking runs in the time domain with a
maximal duration of 8 hours, taken repeatedly in the entirety of the expected
Uranus mission duration.

The first step for the stacking is to be able to generate noise realisations
for single 8-hour tracking passes as a time series, while preserving the
properties defined by the power spectral density. We construct the
hypothetical noise time series $y_{2}^{\rm{noise}}$ of an entire 10 year
tracking run by performing a discrete inverse Fourier transform of
$\sqrt{fS_{n}}$, where we take the Cassini noise spectral density as shown in
Fig. 2:

$\displaystyle y_{2}^{\rm{noise}}(t)=\frac{\sigma_{\rm A}(10\,{\rm
s})}{\mathcal{N}}\times$
$\displaystyle\sum_{f_{i}\in\left[f_{\rm{min}},f_{\rm{max}}\right]}\sqrt{f_{i}S_{\rm
n}(f_{i})}\left[\sin(f_{i}t+\alpha_{i})+\cos(f_{i}t+\beta_{i})\right].$ (15)

Here $f_{\rm{min}}=(10\,{\rm{yr}})^{-1}$ is the inverse of the cruise time
duration and $f_{\rm{max}}=\tau^{-1}=(10\,{\rm{s}})^{-1}$ is the sampling
rate. The phases $\alpha_{i}$ and $\beta_{i}$ are uniformly sampled from
$[0,2\pi)$, resulting in an approximately Gaussian noise distribution.
Additionally, the normalisation factor $\mathcal{N}$ simply assures that the
standard deviation of the noise takes the correct value $\sigma_{\rm
A}(10\,{\rm s})$, regardless of the specific discretisation.

Figure 3: The nth (with $n=\\{1,..,9\\}$) two-way frequency fluctuation
$y_{2,n}$ due to the total signal (black lines), due to the noise (salmon
lines) and due to the GW signal (dashed lines) for 9 different 8-hour
measurements conducted with a 1-year separation. The 9 equally separated
sections between two-way light travel times of 6000 s and 20000 s, represent
roughly Jupiter’s and Uranus’ semi-major axes, respectively. The horizontal
axis has the units of the resolution time (10s), such that a single 8-hour
observation lasts $8\times 60\times 60{\rm s}/(10{\rm s})=2880$ units. The
observations are then glued together in the time axis, hence the
discontinuities in the waveform of the GW. The change in the shape of the GW
signal is due to the variation of the spacecraft-Earth vector with respect to
the GW source as the spacecraft cruises interplanetary space. The mock
amplitude of the GW ($A=10^{-15}$) is exaggerated to guide the eye (see Eqs.
(3)).

We can simulate tracking data from a virtual mission by averaging evaluations
of Eq. (15) and adding it to the coherent Doppler response to a GW signal over
the time segments during which the Doppler tracking is active. Thus, we model
the mission’s total tracking data as being composed of a sum of the GW
response and the noise realisation, stacked over $N$ consecutive tracking
runs:

$\displaystyle y_{2}^{\rm{tot}}=\bigoplus_{n\leq N}\big{(}y_{2}^{\rm
GW}(t_{n};t)+y_{2}^{\rm noise}(t_{n};t)\big{)},$ (16)

where the tensor sum symbol describes the stacking of several, independent 8 h
passes with arbitrary gaps between them. Here $t_{n}$ corresponds to the
timestamp at which the $n$-th tracking run is initiated, and $0<t<8$ hours
follows a clock that is started anew at the beginning of every tracking run. A
possible realisation of the virtual tracking data, where the gaps have been
removed, is visualised in Fig. 3. We note that from here on, we will assume
that the tracking runs are spaced regularly over the mission duration (see
Section VII). While this is most likely not realistic, investigating the
distribution of tracking runs is only useful once an official mission
trajectory is determined. To summarise, in addition to the Allan deviation,
each realisation of $y_{2}^{\rm{tot}}$ is thus determined by the parameters
characterising both the GW and the spacecraft trajectory. Schematically, they
are given by the following:

$\displaystyle\text{GW }\to$
$\displaystyle\begin{cases}\hat{\bm{k}};&\text{Propagation vector}\\\
f_{\scriptscriptstyle\rm GW};&\text{Frequency}\\\ A;&\text{Amplitude}\\\
\phi_{0};&\text{Phase}\end{cases}$ $\displaystyle\text{Trajectory}\to$
$\displaystyle\begin{cases}\hat{\bm{n}}(t_{n};t);&\text{Link unit vector}\\\
T_{2}(t_{n};t);&\text{Two-way travel time}\end{cases}$

##  III GW sensitivity curves for the prospective Uranus Mission

###  III.1 Detection thresholds and null hypothesis

We shall estimate how well the GW signal can be extracted from the tracking
noise. In a Bayesian statement of the problem, we are looking to evaluate the
posterior probability of a GW with parameters
$\left[\hat{\bm{k}},f_{\scriptscriptstyle\rm GW},A,\phi_{0}\right]$ to be
present in a given realisation of $y_{2}^{\rm{tot}}$. Or in other words, we
wish to evaluate and maximise a likelihood function $\mathcal{L}$ of the form:

$\displaystyle\log\left(\mathcal{L}\right)\propto-\sum\left(y_{2}^{\rm{tot}}-y_{2}^{\rm{trial}}\right)^{2},$
(17)

where $y_{2}^{\rm{trial}}$ is a trial GW response. To this end, we employ a
Markov Chain-Monte Carlo (MCMC) framework, where we use the Monte Carlo
sampling package emcee [44] running 32 parallel walkers in order to recover
the posteriors of the GW parameters
$\left[\hat{\bm{k}},f_{\scriptscriptstyle\rm GW},A,\phi_{0}\right]$. Note
that, since the trajectory parameters
$\left[\hat{\bm{n}}(t_{n};t),T_{2}(t_{n};t)\right]$ only evolve very slowly,
we can essentially treat them as constant over a single 8-hour tracking run.
Therefore, we only have to specify the mission’s trajectory at the timestamps
$t_{n}$:

$\displaystyle\hat{\bm{n}}(t_{n};t)$ $\displaystyle\to\hat{\bm{n}}(t_{n}),$
(18) $\displaystyle T_{2}(t_{n};t)$ $\displaystyle\to T_{2}(t_{n}).$ (19)

This simplification not only drastically expedites the likelihood evaluations,
but also easily allows to efficiently repeat the analysis with different
trajectories. This will become especially convenient once the official mission
trajectory becomes available.

Figure 4: We show the two-dimensional marginalised posterior distributions in
amplitude and frequency for a set of three different injected GW amplitudes,
corresponding to an SNR of 0, 3, and 12 from left to right. Here the SNR is
calculated according to Eq. (24) with respect to the threshold amplitude shown
in Fig. 5. As we can see, increasing SNR produces better parameter estimation
in both frequency and amplitude of the GW, where sources with SNR=3 can
already be constrained. The upper bound on the GW amplitude suffers from a
degeneracy with the source’s position in the sky (not shown) and therefore
stops improving significantly after an SNR$\sim 8$ for this particular mission
trajectory.

Our ultimate goal is to obtain a GW sensitivity curve for the mission which
can be easily compared to the sensitivity curves of other GW observatories. A
possible way to define such a curve is to attempt to maximise the likelihood
of a tracking time series that does not contain any GW signal. The MCMC
walkers will be able to freely sample regions in which the trial GW amplitude
$A^{\rm{trial}}$ is well below the noise. However, they will be strongly
discouraged to sample large values of $A^{\rm{trial}}$, since they are
inconsistent with the absence of a signal. Thus, the posterior probability for
the GW amplitude $A$ will sharply decrease above a certain threshold $A_{s}$:

$\displaystyle\mathcal{P}(A)\sim\begin{cases}{\text{const.}}&{\text{for
}}A<A_{s}\left(\hat{\bm{k}},f_{\scriptscriptstyle\rm GW},\phi_{0}\right),\\\
0&{\text{for }}A>A_{s}\left(\hat{\bm{k}},f_{\scriptscriptstyle\rm
GW},\phi_{0}\right),\end{cases}$ (20)

After marginalising over the initial phase $\phi_{0}$ and wave vector
$\hat{\bm{k}}$, we can specify the threshold to the desired confidence level
and obtain a curve $A_{s}(f_{\scriptscriptstyle\rm{GW}})$, purely as a
function of GW frequency. An example of this procedure is visualised in the
left panel of Fig. 4, where one can clearly see a transition between allowed
and disallowed regions of parameter space in the 2-dimensional posteriors of
the GW frequency and amplitude.

As seen in the bottom panel of Fig. 5, the boundary between the allowed and
disallowed regions is closely approximated by a half-Gaussian. This allows us
to easily find the $1\sigma$ confidence threshold, above which signals are
distinguishable from the null hypothesis. Note that this does not necessarily
imply that the individual signal parameters are well constrained, only that
the overall signal is distinguishable from its absence. This corresponds to
the analogous definition of a sensitivity curve for other GW detectors (see,
e.g., [87]). For the purposes of defining a simple sensitivity curve, we fit
the one sigma boundary with a broken power law, shown in the top panel of Fig.
5.

The value of the $1\sigma$ threshold $A_{s}$ depends on the noise power
spectral density of the tracking runs, the specified mission trajectory, as
well as the total number of tracking runs $N_{\rm t}$. By performing several
numerical MCMC tests, we find, as one would expect, that its overall
normalisation scales linearly with the Allan deviation and scales as the
square root of the total number of samples:

$\displaystyle A_{s}\sim\sigma_{\rm
A}(\tau)\times\frac{\sqrt{\tau}}{\sqrt{N_{\rm t}}}.$ (21)

The effect of the spacecraft’s trajectory is also folded into the specific
realisation of $A_{s}$, as it modulates the GW response function
$y_{2}^{\scriptscriptstyle\rm{GW}}$ over the 10-year observation window. In
particular the relative orientation of the wave vector $\hat{\bm{k}}$ and the
link vector $\hat{\bm{n}}$ can change significantly for elliptical
trajectories. The presence of curvature in the trajectory is crucial, as it
allows one to break the sky localisation degeneracy that would be present for
a single, straight detector arm.

Figure 5: In the top panel we show the two-dimensional marginalised posterior
samples in amplitude and frequency for a run with no injected signal. We
highlight the threshold at which the absence of a signal is excluded at
$1\sigma$ confidence, which is well fit by a broken power law. In the lower
panel, we showcase the sample count in 50 frequency bins, as a function of
amplitude. We normalise the latter with the $1\sigma$ threshold.

###  III.2 Defining the GW sensitivity curve

To reiterate, the recovered $1\sigma$ threshold $A_{s}$, which corresponds to
the definition of a sensitivity curve, results from distributing $N_{\rm t}$
tracking runs over the entire mission duration of $T=10$ years. For a
monochromatic signal, this results in tracking a total of
$n\sim(f_{\scriptscriptstyle\rm{GW}}\times N_{\rm t}\times 8\,\rm hour)$
cycles. This is reminiscent of matched filtering techniques, which enable one
to compare the sensitivity curves of GW observatories to the characteristic
strain of a signal, rather than its instantaneous amplitude [see, e.g., 87].
Thus, we adjust the $1\sigma$ detection threshold by heuristically factoring
out the number of detected cycles:

$\displaystyle
h_{\rm{sens}}(f_{\scriptscriptstyle\rm{GW}})=A_{s}(f_{\scriptscriptstyle\rm{GW}})\times\left(f_{\scriptscriptstyle\rm{GW}}\times
N_{\rm t}\times 8\,{\rm{hr}}\right)^{1/2},$ (22)

where we define $h_{\rm{sens}}(f_{\scriptscriptstyle\rm{GW}})$ the GW
sensitivity curve of the mission, given a trajectory, an Allan deviation and a
set of tracking day timestamps $t_{n}$:

$\displaystyle
h_{\rm{sens}}(f_{\scriptscriptstyle\rm{GW}})=h_{\rm{sens}}(f_{\scriptscriptstyle\rm{GW}};\sigma_{\rm
A}(10\,{\rm{s})},t_{n}).$ (23)

Then, we can also define the SNR of a GW detected with Doppler tracking in
complete analogy with other gravitational wave detectors:

$\displaystyle{\rm{SNR}}=\sqrt{4\times\int\frac{h_{\rm
c}(f^{\prime})^{2}}{h_{\rm{sens}}(f^{\prime})^{2}f^{\prime}}{\rm
d}f^{\prime}},$ (24)

where $h_{\rm c}$ is the GW’s characteristic strain, and the factor 4 comes
from the normalisation of the one-sided-power spectrum. Here, an SNR of $1$
corresponds to a barely detectable GW that can be distinguished from the
absence of signal with $1\sigma$ confidence. As shown in Fig. 4, higher SNRs
result in both better parameter estimation and in reduced false alarm rates,
also in complete analogy with GW detectors. We find that the frequency and
amplitude of sources with an SNR$\sim 3$ can be recovered within $\sim 0.02$
and $\sim 0.5$ dex, respectively. The bounds improve by a factor $\sim 3$ for
SNR $=12$ sources, though a residual degeneracy between the source’s amplitude
and optimal localisation on the sky becomes manifest, inducing a skew in the
recovered posterior distribution.

###  III.3 Baseline, Priority and Optimistic scenarios

The sensitivity curves constructed above span the interesting frequency range
of few $\sim 10^{-9}$ Hz to $\sim 10^{-5}$ Hz with a slope of approximately
$f^{-1/3}$, before a turning point at $f\sim$1/(8 hours). They then rise with
a slope of $\sim f^{1/3}$ until it is truncated by the Nyquist frequency at
$\sim 10^{-1}$ Hz. Crucially, strain sensitivities of the order $\sim
10^{-15}$ in the nano-Hz to milli-Hz band can be achieved without requiring
significant improvements in the Allan deviation, provided that tracking data
is accumulated over a substantial amount of the entire mission duration. This
is in stark contrast with previous analyses, which did not include the
stacking procedure detailed in Section II.5. We note that this amount of
tracking is not typical for deep interplanetary missions. Here, however, we
argue that it should be prioritised in parallel with Allan deviation
improvements in order to achieve the scientific objectives detailed in this
manuscript. Assuming this level of commitment to the tracking, we define three
possible future scenarios different technological improvements are implemented
with respect to Cassini era ranging.

In the baseline scenario, we do not assume any improvement with respect to
Cassini-era tracking technology. We assume the tracking of the spacecraft for
8 hours each day over the mission’s 10-year interplanetary cruise. The Allan
deviation of this scenario is then $\sigma_{\rm A}(10\,\mathrm{s})=10^{-14}$,
at a sampling rate of 10 seconds. The total number of 8-hour tracking runs is
approximately 3650.

In the priority scenario, we assume that optimising currently available
tracking technology becomes one of the explicit targets of the mission before
launch. We assume that the Allan deviation at 10 seconds is reduced to the
optimal Cassini-era value of $3\times 10^{-15}$, which was achieved with a
longer averaging time of $\tau=10^{4}$ s.

In the optimistic scenario, we assume that technological improvements in
tracking technology (see, e.g., Sec. VII) will have acquired an appropriate
readiness level for an interplanetary mission. As a proxy for such
technological improvements, we assume an Allan deviation at 10 seconds of
$10^{-15}$. Additionally, we assume that several tracking stations distributed
across the globe be utilised, assuring 24 hours of tracking per day. We
discuss the prospects to achieve such low levels of noise in Section VII.2.
Here we simply note that by the time the UOP is operational, almost half a
century of technological improvement will have taken place since the launch of
the Cassini mission.

##  IV Detection forecast for Black Hole binaries

Supermassive black hole binaries (SMBHBs) are among the loudest sources in the
low-frequency GW sky. They are expected as a consequence of galaxy mergers,
although how frequently they form and how they evolve in their environments is
still not fully understood [see, e.g., the seminal papers 12, 104, 99, 55,
64]. Recent constraints on a GW background signal by pulsar timing arrays
(PTA) provide evidence of an SMBHB population in the mass range $\sim
10^{8}-10^{10}M_{\odot}$ [2, 6, 85, 108] while the direct detections of
lighter SMBHB mergers below a few $10^{7}M_{\odot}$ is expected in the 2030s
with the advent of space-based interferometers such as LISA [29] and TianQin
[68]. As we will show, the data stacking method in Section III facilitates the
detection of SMBHBs that are evolving beyond the nano-Hz band and into the
milli-Hz regime. As an illustration, we show the hypothetical sensitivity of
the mission to chirping, equal-mass BH binaries as a waterfall plot in the
bottom panel of Fig. 6. The mission’s peak sensitivity lies in the total mass
range of $10^{8}$ to $10^{10}M_{\odot}$. Such heavy, chirping binaries would
in principle be detectable up to redshifts of $z\sim 10$ for an Allan
deviation of $\sigma_{\rm A}=10^{-14}$ and essentially over the whole
cosmological volume for $\sigma_{\rm A}=3\times 10^{-15}$.

Figure 6: In the top panel, we compare the GW sensitivity of the prospective
Uranus mission to the current sensitivity of PTAs [36] and the sensitivity of
the recently adopted LISA mission. The three blue lines correspond to the
scenarios detailed in Section III. In the bottom panel, we show a waterfall
plot of the mission’s hypothetical detection horizon for chirping equal mass
ratio massive black hole binaries, assuming a maximum of 10 years for the
signal duration and the priority scenario, as described in Section III.3. We
compare it to the same plot for PTAs and for LISA (assuming 4 years of maximal
signal duration). Note that monochromatic sources, expected for PTAs and from
our population models in §IV, will have lower SNR for the same mass and
redshift in the bottom panel.

However, a detection forecast requires us to consider the actual population of
sources we expect to fall in this frequency band, and whether the majority
will be detectable as a monochromatic signal rather than a chirp. Here we
adopt two separate approaches to estimate the redshift-dependent merger rate
of such sources.

###  IV.1 Model 1: Millennium simulation

Here we detail a simple prescription that links the massive black hole merger
rate to the much more established halo merger rate, based on the two Millenium
simulations [38]. The latter works provide a convenient fit characterising the
differential halo merger rate:

$\displaystyle\frac{{\rm d}^{2}\Gamma}{{\rm d}\xi{\rm
d}z}=B_{1}\left(\frac{M_{\rm{halo}}}{10^{12}\,{M}_{\odot}}\right)^{b_{1}}\xi^{b_{2}}\exp\left[\left(\frac{\xi}{B_{2}}\right)^{b_{3}}\right](1+z)^{b_{4}},$
(25)

where $\Gamma$ is the total number of mergers that a halo of mass
$M_{\rm{halo}}$ experiences over cosmic time, $\xi\leq 1$ is the halo merger
mass ratio and the best-fit parameters are given by
$(B_{1},B_{2},b_{1},b_{2},b_{3},b_{4})=(0.0104,9.72\times
9.72,0.133,-1.995,0.263,0.0993)$. We link Eq. (25) to the SMBH merger rate
$\dot{N}_{\bullet\bullet}$ by multiplying the halo merger rate with the black
hole mass function:

$\displaystyle\frac{{\rm d}^{3}\dot{N}_{\bullet\bullet}}{{\rm
d}M_{\bullet}{\rm d}\xi{\rm d}z}$ $\displaystyle={P_{\rm
occ}}(M_{\rm{halo}},z)\frac{4\pi cD^{2}_{\rm{com}}(z)}{(1+z)^{3}}$
$\displaystyle\times\frac{{\rm d}n_{\bullet}}{{\rm d}M_{\bullet}}\frac{{\rm
d}^{2}\Gamma}{{\rm d}\xi{\rm d}z}(\xi,z_{\rm{del}}),$ (26)

where $D_{\rm{com}}$ is the comoving distance at redshift $z$ and we must
additionally supply an occupation fraction $P_{\rm occ}$ of black holes in
halos and a delay prescription $z_{\rm{del}}$ between the nominal halo merger
time and the actual black hole merger time. Here we use the SMBH mass function
${{\rm d}n_{\bullet}}/{{\rm d}M_{\bullet}}$ as reported in [92] and adopt a
simple relation between the halo and SMBH mass from [31]:

$\displaystyle M_{\bullet}$
$\displaystyle=\left[\frac{M_{\rm{halo}}(1+z)}{2\times
10^{7}\,M_{\odot}}\right]^{3/2}\,M_{\odot},$ (27)

which is consistent with both simulations and observations of massive galaxies
[70, 10].

We can now establish the number of detectable sources of gravitational
radiation by integrating the differential SMBH merger rate. We define an SNR
threshold $C$ and count all SMBH binaries that exceed it:

$\displaystyle N_{\bullet\bullet}^{\rm{det}}=\int\int\int\int\frac{{\rm
d}^{3}\dot{N}_{\bullet\bullet}}{{\rm d}M_{\bullet}{\rm d}\xi{\rm
d}z}\frac{1}{\dot{f}_{\scriptscriptstyle\rm{GW}}}\Theta\left({\rm{SNR}-C}\right)$
(28) $\displaystyle\times{\rm d}f_{\scriptscriptstyle\rm{GW}}\,{\rm d}\xi{\rm
d}M_{\bullet}{\rm d}z,$

where $\Theta$ is a Heaviside function, and we distribute all binaries in
frequency bins according to their residence times
$\dot{f}_{\scriptscriptstyle\rm{GW}}^{-1}$. Here it is necessary to make two
additional simplifications. Firstly, we model the occupation fraction $P_{\rm
occ}$ as an unspecified constant. For massive black holes $M>10^{7}M_{\odot}$,
simulations [86] and observations of AGN [48] show that the occupation
fraction $P_{\rm occ}$ is a fraction of order unity. Secondly, we also neglect
the effect of time delays between the halo merger and the BH merger, as they
typically only amount to $10^{8}$–$10^{9}\,\rm years$ for the mass range we
are considering [103, 35]. This simplification is justified a posteriori,
given that tracking will only ever be sensitive to binaries at $z\lesssim 1$
for realistic populations,555The seeming discrepancy with the detection
horizon in Fig. 6 derives from the detectability of monochromatic sources vs.
the hypothetical chirping sources drawn for consistent comparison between
bands in Fig. 6. where a delay of 109 years at most corresponds to a small
redshift shift $\Delta z_{\rm{del}}\sim 0.3$.

###  IV.2 Model 2: Holodeck

We present another estimate of the detectable SMBHB population via the
holodeck code [61, and see Ref. [3]]. Importantly, we utilise holodeck to
build population models that generate nano-Hz-frequency GWBs with an amplitude
that is consistent with the GWB for which multiple PTA experiments have
recently found strong evidence [2, 6, 85, 108].

We utilise holodeck’s semi-analytical SMBHB population models. Specifically,
these models combine an observationally constrained galaxy mass function from
Ref. [66], with a galaxy merger rate from Ref. [88] to generate galaxy
mergers. An SMBH-mass versus bulge-mass relation from Ref. [64] is then used
to calculate SMBHB parameters, which are then self-consistently evolved from
galaxy-scales to nano-Hz separations using a phenomenological binary evolution
model developed for holodeck that combines binary decay due to environmental
interactions and GWs [3]. The change in source redshift over the course of
binary evolution is naturally included. Each of these model components depends
on a number of uncertain parameters with measurements and uncertainties
described in the cited studies.

To sample SMBHBs that are consistent with PTA observations, we generate 1080
holodeck populations by sampling over the uncertainties in each model
parameter66629 parameters are marginalised over in this way, sampling from the
PS$\\_$Astro$\\_$AStrong$\\_$AAll-v1.0 parameter space.. We then calculate the
resulting PTA-measured GWB spectrum, fit a power-law spectrum to the lowest
five frequency bins, and select the populations with amplitudes and spectral
indices consistent with recent measurements In particular, we require the GWB
amplitude (characteristic strain at a frequency of $1\,\mathrm{yr}^{-1}$) that
is $A_{\mathrm{yr}}\in[10^{-15},10^{-14}]$, and a power-law index (in power-
spectral-density of timing residuals) that is $\gamma\in[-5.0,-2.5]$ [see,
e.g., 100]. This procedure is repeated for each UOP sensitivity scenario and
threshold SNR, drawing independent populations for each combination. Out of
the 1080 populations, $130\pm 10$ pass the above selection cuts. Those
populations are then integrated over with a cut in detection SNR similarly to
Eq. (28) to determine detectable SMBHB systems.

###  IV.3 Detection forecast

The results of our detection forecast are summarised in Table 1 and Fig. 7.
Regarding the former, we compute the number of expected detections of
individual binary SMBHs at three separate SNR thresholds (1, 3, 8), for the
three mission scenarios of baseline, priority, and optimistic (§III.3). For
the holodeck results we quote the median, $25\%$, and $75\%$ quantiles derived
from the $\sim 130$ GWB-selected populations described above. For the
Millennium results, we quote expectation values in terms of the occupation
fraction of BHs in halos. The expected number of detections is broadly
consistent for the two adopted population models; the Millennium expectation
values fall within the inner $50\%$ range of holodeck values for
$P_{\rm{occ}}\gtrsim 0.4$. Even in the baseline scenario, which does not
assume any technological improvement with respect to Cassini-era technology,
we would typically expect the detection of $\mathcal{O}(1)$ binaries at ${\rm
SNR}=3$. A mild improvement in the Allan deviation, e.g., to $3\times
10^{-15}$, boosts the number of detections to $\mathcal{O}(10)$. Overall, the
results shown in Table 1 suggest that a mild improvement in Allan deviation
over the baseline mission would guarantee the detection of several well-
constrained monochromatic signals from SMBH binaries, in a mass range that is
yet to be explored. Indeed, as seen in Fig. 7, the majority of detected
sources would lie around a mass of $10^{9}M_{\odot}$, being detectable out to
redshifts of $z\lesssim 1$. The majority of these are detected as almost
monochromatic sources at frequencies of tens to hundreds of nHz, i.e., above
the PTA band (see bottom panel of Fig. 7). A more sensitive mission can pick
up more monochromatic signals at lower frequencies, where a larger number of
binaries reside.

The main difference between the Millennium and the holodeck population models
lies in the redshift distribution of sources, as can be seen in the middle
panel of Fig. 7. As seen in Fig. 8, we also note that in the priority
scenario, the UOP will be sensitive to the stochastic background of SMBH
binaries over the entire frequency range of $10^{-8}$ Hz to $10^{-6}$,
provided that the latter follows the scaling of $f^{-2/3}$ that is reported
by, e.g., Ref, [2].

As a final comment on the topic of SMBHBs, we briefly mention that tracking
the UOP would provide stringent bounds on the existence of an intermediate-
mass companion to the main SMBH of the Milky Way, i.e., SgrA*. Current
constraints arise from a range of techniques, see, e.g., [106, 78]. Quite
relevant for this work are the models where the GWs emitted fall in the micro-
Hz band. In [78], it was shown that several of these models generate GWs with
strains at the Solar system of up to $10^{-12}$. From Fig. 6, it can be easily
extrapolated that all hidden companions that would produce strains of order
$\sim 10^{-15}$ over the frequency range of $10^{-8}$ Hz to $10^{-5}$ Hz would
be detectable.

Combining all of these prospects, we have demonstrated that Doppler tracking
with the UOP is uniquely suited to detect GWs in the frequency regime between
PTAs and space-based interferometry. Sources in this frequency range comprise
heavy SMBH binaries of several $10^{9}M_{\odot}$ approaching coalescence,
inspiraling pre-LISA systems and intermediate-mass black holes in the Milky
Way. Such detections would unveil a more complete picture of SMBH formation,
growth, and interaction, complementing both existing and planned GW detectors
in the low-frequency regime. In the event of a concurrent detection with PTAs
or LISA, a UOP detection may also aid in the difficult task of localising
SMBHB sources [e.g., 97].

Figure 7: The differential distribution of detections (i.e., the probability
density function) shown with the dashed curves and the cumulative distribution
function with the solid curves, for both the Millennium and holodeck (here
chosen to be the median percentile realisation) population models, normalised
with their respective maximum values. The top panel shows the distribution in
mass decades $dN/d\log_{10}(M)$, while the middle panel shows the distribution
in redshift $dN/dz$. The curves are displayed for the optimistic scenario with
an SNR threshold of 3, but apply with little variation for all scenarios and
SNR thresholds. The lowest panel shows the distribution of sources in
$dN/d\log_{10}(f_{\rm GW})$, showing how a more sensitive mission is capable
of detecting more monochromatic signals at lower frequencies, shifting the
PDFs towards the left. The highest SNR sources are therefore always detected
at $\sim 100$ nHz, above the PTA sensitivity band. Note that all the displayed
curves are normalised.

Scenario 1: Baseline Scenario 2: Priority Scenario 3: Optimistic Pop. SNR Nr.
Det Millennium 1 12$\times P_{\rm occ}$ Millennium 3 2$\times P_{\rm occ}$
Millennium 8 0 holodeck 1 $6.8^{+20.0}_{-2.0}$ holodeck 3 $1.0^{+1.8}_{-0.3}$
holodeck 8 $0.15^{+0.28}_{-0.06}$ Pop. SNR Nr. Det Millennium 1 71$\times
P_{\rm occ}$ Millennium 3 12$\times P_{\rm occ}$ Millennium 8 2$\times P_{\rm
occ}$ holodeck 1 $58.2^{+193.1}_{-21.4}$ holodeck 3 $6.5^{+14.8}_{-1.7}$
holodeck 8 $1.2^{+3.1}_{-0.4}$ Pop. SNR Nr. Det Millennium 1 1049$\times
P_{\rm occ}$ Millennium 3 197$\times P_{\rm occ}$ Millennium 8 42$\times
P_{\rm occ}$ holodeck 1 $1375^{+2582}_{-584}$ holodeck 3 $190^{+497}_{-68}$
holodeck 8 $30.1^{+67.2}_{-10.6}$ .

Table 1: Three tables summarising the expected number of detections of SMBHBs
for the UOP scenarios detailed in Section III.3, over the entire mission
duration. The detection estimates are performed with the two population models
explained in Section IV and reported for three different SNR thresholds (see
Fig. 4 for a visualisation of the typical parameter posteriors given these
SNRs). The Millennium results scale with the occupation fraction of SMBHs in
massive galaxies, $P_{\rm occ}$. The holodeck numbers are quoted as median
$\pm$ $25\%$ and $75\%$ quantile values.

##  V Forecast constraints on early Universe signals

As well as filling the gap in our sensitivity to SMBHBs, our proposed GW
searches are also sensitive to potential cosmological signals from the early
Universe. There is a huge range of mechanisms for generating GWs at early
times [25], each of which, if detected, would give valuable new insights into
fundamental physics. In the context of Doppler tracking with the UOP however,
the most interesting prospects are for signals that are peaked in a relatively
small range of frequencies, which risk going undetected by other GW
experiments if there is not sufficient coverage of the micro-Hz frequency
band.

The quintessential example of such a peaked cosmological GW signal is that
generated by a first-order phase transition (FOPT) in the early Universe.
These transitions, in which a fundamental quantum field escapes from a
metastable state by nucleating ‘bubbles’ of a new phase, are a generic
prediction of many extensions to the Standard Model of particle physics [25,
26, 27], notably beyond the reach of even the most powerful particle
accelerators on Earth. Collisions between nucleated bubbles, and the
subsequent acoustic and turbulent motion of the thermal plasma, can produce
strong SGWB signals that peak at a frequency set by the energy scale of the
transition, with higher frequencies corresponding to higher energies and thus
earlier cosmic epochs [60]. This peak frequency is commonly approximated as

$f_{*}\approx
19\,\mu\mathrm{Hz}\times\frac{k_{\mathrm{B}}T_{*}}{100\,\mathrm{GeV}}\frac{\beta/H_{*}}{v_{w}/c}\left(\frac{g_{*}}{106.75}\right)^{1/6},$
(29)

where $T_{*}$ is the temperature of the thermal plasma at the time when the
GWs are generated, $\beta$ is the duration of the transition, measured in
units of the Hubble rate $H_{*}$ at that epoch, $v_{w}$ is the terminal
expansion velocity of the bubble walls, and $g_{*}$ is the number of
relativistic degrees of freedom in the plasma, normalised here to the Standard
Model value. Focusing for simplicity on the contribution from sound waves in
the plasma [52], which is dominant in many scenarios, the SGWB spectrum can be
modelled as a broken power law,

$\Omega_{\mathrm{GW}}(f)\approx\Omega_{\mathrm{GW}}(f_{*})\times(f/f*)^{3}\left[\frac{7}{4+3(f/f_{*})^{2}}\right]^{7/2},$
(30)

where the peak intensity is

$\displaystyle\begin{split}\Omega_{\mathrm{GW}}(f_{*})\approx 5.7\times
10^{-6}&\times\frac{v_{w}/c}{\beta/H_{*}}\left(\frac{\kappa\alpha}{1+\alpha}\right)^{2}\left(\frac{g_{*}}{106.75}\right)^{-1/3}\\\
&\times\left[1-(1+2\tau_{\mathrm{sw}}H_{*})^{-1/2}\right],\end{split}$ (31)

with $\kappa$ an efficiency parameter and $\tau_{\mathrm{sw}}$ the lifetime of
the sound-wave source, both of which can be modelled as functions of the
parameters $\alpha$, $\beta$, and $v_{w}$ [26, 27].

Following [18, 19], we carry out a scan over the 4D parameter space
$(T_{*},\alpha,\beta,v_{w})$ to assess the sensitivity of our proposed
searches to FOPT signals. As shown in Fig. 9, there is a broad swathe of
parameter space for these models that is accessible with Doppler tracking
searches in the micro-Hz band that would go undetected by pulsar timing arrays
in the nano-Hz band and LISA in the milli-Hz band. These constraints are
complementary to those that have been forecasted for ‘binary resonance’
searches for micro-Hz GWs using Lunar laser ranging data [18, 19], allowing us
to potentially detect or rule out a much larger family of FOPT models.

Figure 8: Forecast sensitivity to the SGWB spectrum $\Omega_{\mathrm{GW}}(f)$
for our baseline, priority, and optimistic scenarios (navy solid, dashed, and
dotted curves), along with existing sensitivities of the NANOGrav PTA (green)
[3] and Doppler tracking with Cassini (blue), as well as forecast
sensitivities of LISA (brown) [87] and binary resonance searches with lunar
laser ranging data (LR, light blue). The black horizontal line ($N_{\rm eff}$
constraint) shows the indirect constraint on the _integrated_ SGWB intensity
from measurements of the Hubble rate during radiation domination. The pink
shaded region indicates the astrophysical SGWB from unresolved SMBHBs (with
amplitude inferred from NANOGrav and fixed $\propto f^{2/3}$ spectral tilt, as
appropriate in the absence of environmental or relativistic effects). The
yellow shaded region indicates a potential cosmological signal from a first-
order phase transition in the early Universe.  Figure 9: Forecast constraints
on cosmological SGWB signals from early Universe phase transitions. These are
shown using 2D projections of the 4D space of parameters describing the
transition: the temperature of the Universe at the time of the transition,
$T_{*}$, the fractional energy density released by the transition, $\alpha$,
the rate of the transition, $\beta$, and the terminal velocity of the bubble
walls, $v_{w}$.

##  VI Forecast on ULDM constraints.

The expected constraints on $h_{c}$ from the top panel of Fig. 6 already give
an idea of the reach of UOP to detect ULDM (compare Eqs. (1) and (7)). This
suggests that for masses $m_{\scriptscriptstyle\rm DM}\gtrsim 10^{-22}$
eV/$c^{2}$, UOP may represent the most sensitive probe for large regions in
parameter space of dark matter models, even in the baseline scenario. Notably,
this mass range is particularly relevant for ULDM phenomenology [57, 41]. To
illustrate this reach, let us consider the equivalent $h_{c}$ that would be
produced by the three families of models we introduced in Section II.3, and
that could be summarised by Eqs. (8), (9) and (11). Recall that $f_{\rm
grav}\approx 2f_{\scriptscriptstyle\rm DM}$, where $f_{\scriptscriptstyle\rm
DM}$ can be found in Eq. (6) for the purely gravitational coupling and the
quadratic direct coupling $\alpha(\phi)=\phi^{2}/\Lambda^{2}$. For the linear
coupling, such as the one in Eq. (10), $f=f_{\rm DM}$.

The purely gravitational case represented by Eq. (8) (which is independent of
other couplings) generates the black line in Fig. 10 as a function of the ULDM
mass (horizontal axes). We notice that the priority UOP scenario may reach a
sensitivity that would directly detect dark matter in the Solar System for
masses below $\sim 10^{-23}$ eV. While these low masses are already in tension
with other astrophysical observations (see, e.g., [11, 41]), this is a very
interesting result since, as compared to other bounds, it would be based on
_direct_ sensitivity to the ULDM field.

Regarding direct coupling, from Eq. (9) and the bounds on the gravitational
effect from Fig. 10, it is easy to understand the constraining power on
$\Lambda_{2}$ for different masses. In the upper panel of Fig. 10 we show the
equivalent $h_{c}$ predicted for two values of $\Lambda_{2}$: one
corresponding to the typical minimum value UOP could constrain (around
$\Lambda_{2}=10^{-18}\,\rm GeV$) and another one that improves on the best
bounds from Cassini ($\Lambda_{2}=10^{-1}6\,\rm GeV$). We see clearly that in
both cases, even the baseline UOP will produce the most stringent results for
ULDM in a large range of masses. To make this more explicit, in Fig. 10, lower
panel, we show the constraints that can be achieved on $\Lambda_{2}$ as a
function of the DM mass for different configurations and compare them to
Cassini and PTA [21] data (see also [20] and the recent [65, 95] for related
bounds). It is clear that the UOP will enormously extend the searches for
these models, with the potential to generate the first direct detection of
dark matter.

Finally, we also show in Fig. 10 the equivalent $h_{c}$ from the effect from
Eq. (11) representing axion-like particles. Recall that this effect also
arises from an oscillating part of the ULDM with frequency $f_{\rm DM}$. We
have used a value of the coupling $g_{a\gamma\gamma}$ characteristic of
current bounds at the relevant masses, see [80]. From this plot, we conclude
that the priority or optimistic cases of the UOP mission would also have a
unique legacy in the search for axion-like particles.

Figure 10: Upper plot: Expected sensitivity of the UOP for ULDM. The effect of
ULDM is represented as a corresponding $h_{c}$ that can be directly compared
with Fig. 6. We plot the gravitational case (Eq. (8)), two cases where ULDM
couples to matter (Eq. (9)), and a case where ULDM couples to light (Eqs.
(11)). Lower plot: Expected sensitivity of FUM for ULDM for the model with
direct coupling of matter to $\alpha(\phi)=\phi^{2}/\Lambda_{2}^{2}$. We
compare with bounds from Cassini [9] and the EPTA [94] (see also [83]).

##  VII Beyond the Baseline: Reducing the Allan deviation

###  VII.1 Noise sources

Predicting realistic improvements in the Allan deviation requires considering
recent advancements in a wide range of technologies. Here we base some
qualitative considerations on the thorough discussion in Ref. [7]. Cassini-era
noise was dominated by 3 major components; the antenna mechanical noise, the
plasma interference noise, and the tropospheric scintillation noise. The first
is simply a mechanical limitation, which can be addressed by implementing a
complementary smaller and stiffer antenna in combination with the main dish,
as well as via three-point antenna calibration, which carries the potential to
reduce $\sigma_{\rm A}$ of Ka-band tracking missions by one order of
magnitude.777Private conversation with Sami W. Asmar, NASA JPL.

The noise of plasma scintillation, the variation of the refractive index of
the interplanetary medium, is localised around $3\times 10^{-3}$ Hz and then
steadily drops down for lower frequencies [73]. It is astrophysical in origin,
and therefore cannot be easily addressed with technological improvements. One
possibility is to prioritise tracking measurements at optimal Sun-Earth-
spacecraft configuration, as well as by upgrading the Doppler link to higher
frequencies, as discussed in Section VII.2.

Tropospheric noise is thought to dominate the noise power spectral density
frequency below $10^{-4}$ Hz. Since the Cassini era, corrective measures based
on water-vapor-radiometers have been shown to be able to reduce the Allan
deviation up to a factor of 10, down to $1.5\times 10^{-15}-3\times 10^{-15}$
[7]. Furthermore, drastic improvements in tropospheric noise would be
achievable by tracking the spacecraft with high-altitude facilities, either
ground-based or with balloons. Finally, an ulterior opportunity to reduce
tropospheric noise would be using multiple measurement points via a radio
telescope array [107].

Glitches and unmodelled accelerations, e.g., from solar winds or imprecisions
in the planetary ephemeris, may also leave residual imprints on the Doppler
time series. However, these outliers will not influence the coherent stacking
of the GW signal in the experiment. A more quantitative treatment of their
influence can be addressed by a more detailed simulation of the mission’s
trajectory, once the satellite specifications are known. Connected to this,
one must also consider the noise from asteroids in the Solar system. A simple
estimate using Brownian-motion argument [67] suggests a noise level in the
solar system barycenter of order $10^{-15}$ in the 10 nHz range due to
asteroids smaller than 80 km. Furthermore, the results in [39] indicate a
lower level of $10^{-18}$ for the potential generated by asteroids on a
satellite at $\sim$AU from the Sun (with typical frequencies around and below
$\mu$Hz), which would be irrelevant to the scientific outcome we outlined.
However, a more detailed analysis confirming this expectation and considering
a realistic mission trajectory and possible mitigation measures using updated
asteroid ephemerides (see [102] for the current strategy followed by NANOGrav)
would be opportune.

###  VII.2 Breakthrough: Optical links?

A true breakthrough in improving the Allan deviation to $\lesssim 10^{-15}$
could be achieved by upgrading the Doppler link to optical frequencies. Here
we briefly review theoretical limits and recent experiments. There are three
distinct techniques in optical ranging.

  * •

The passive-reflector method has reflecting cubes on the target, which reflect
light back to the source. This has been very successful for lunar ranging [75]
but (as we will see below) is not feasible at interplanetary distances.

  * •

The synchronous transponder method has the target receive a signal and then
send another signal back with a fixed delay. This has been recently
demonstrated with the Hayabusa2 mission to $\sim 15$ lunar distances [79],
with ground observations even during the day.

  * •

An asynchronous transponder ranging both sides sends signals at pre-designated
intervals, and the ranging is computed indirectly [32]. The Deep Space Optical
Communications (DSOC) system on the Psyche mission [84, 76] is an example.

We can get an idea of what is and is not feasible using order-of-magnitude
arguments. Similar arguments appear in [33]. Consider a signal at wavelength
$\lambda$ from a transmitting telescope of area $A_{T}$ to a receiving
telescope of area $A_{R}$ at distance $R$. A diffraction-limited beam will
have a solid angle of $\sim\lambda^{2}/A_{T}$. Of the photons transmitted from
$A_{T}$, the fraction received at $A_{R}$ will be $\sim A_{R}/R^{2}$ divided
by the solid angle. Thus a fraction

$\frac{A_{T}A_{R}}{\lambda^{2}R^{2}}$ (32)

will be received. For a transmitting power $P$, the rate of photons received
is then

$\frac{P}{hc\lambda}\frac{A_{T}A_{R}}{R^{2}}.$ (33)

The maximum data rate for a given power and area is thus
$\propto\lambda^{-1}$.

For some orders of magnitude, let us put
$\lambda=$1\text{\times}{10}^{-6}\text{\,}\mathrm{m}$$,
$A_{T}=A_{R}=$1\text{\,}\mathrm{m}^{2}$$ and
$R=$3\text{\times}{10}^{9}\text{\,}\mathrm{m}$$. This makes the photon
fraction $\sim 10^{-13}$ and the photon rate
$\sim$5\text{\times}{10}^{4}\text{\,}\mathrm{W}^{-1}$$.

With a passive reflector the photon fraction (32) gets squared through the
round trip. At lunar distances
($R\simeq$3\text{\times}{10}^{9}\text{\,}\mathrm{m}$$) this is tolerable, but
at interplanetary distances, it would become hopeless. Hence an active system
is needed.

Ref. [98] discusses optical telecom prospects for ice-giant missions briefly,
noting the positives and negatives, which also follow from the above
arguments.

  * $+$

The data rates can be orders of magnitude higher. Ref. [96] demonstrated
uplink and downlink rates of 19 and 38 Mbps respectively from lunar distances.
The Psyche mission to the asteroid belt was designed for rates of 250 Mbps
from deep space [23], and a demonstration of the data rate (showing a cat
pursuing a laser pointer) was widely reported in the media and social media.

  * $-$

The narrowness of the optical beam requires accurate pointing, including the
effect of light travel time.

  * $-$

To be weather-proof, the system requires redundant ground stations (or balloon
or orbiting stations).

Laboratory experiments of asynchronous ranging [16] give precisions of better
than 1 mm, which amounts to $3\times 10^{-16}$ of the Uranus distance.
Reference [34] discuss the error budget for interplanetary ranging, and
conclude that there is potential for (sub-)mm range accuracy. That is,
accuracies even better than the optimistic scenario considered earlier appear
feasible with optical ranging.

##  VIII Conclusion: Mission Outlook

We have demonstrated that the prospective flagship mission to Uranus (UOP) has
rich scientific opportunities that go beyond its nominal planetary science
goals. Doppler tracking data may be accumulated over the entire interplanetary
cruise phase, and subsequently stacked according to the methodology showcased
in Section II.5. This results in the capacity to detect imprints from
deterministic GWs from individual SMBH binaries (Section IV), astrophysical as
well as early universe stochastic GW backgrounds (Sections IV and V), and
finally gravitational and non-gravitational couplings between the spacecraft
and dark matter (Section VI). Such varied and timely prospects would greatly
expand the scientific yield of the mission, broadening its relevancy to
communities beyond the planetary sciences. So the question beckons: what is
required to actually achieve what we have shown to be possible in this paper?
We note that our analysis of the mission’s capacity to distinguish various
signals should be repeated with more sophisticated noise models, including,
e.g., non-Gaussianities and the time dependence of different noise
contributions along the cruise phase. A more thorough treatment of the
mission’s trajectory is also crucial, in particular considering the influence
of unmodelled accelerations due to uncertainties in the masses of planets,
asteroids, and solar radiation pressure. Additionally, a more sophisticated
data analysis pipeline should take into account the possibility of multiple
signals to be present at the same time, including the possibility of
stochastic backgrounds overshadowing individual signals. Nevertheless, taking
these limitations at face value, we have demonstrated how two crucial factors
determine the magnitude of the achievable scientific goals. They are:

  * •

The total amount of tracking data collected during the cruise phase, which
could range from a few data points to constant tracking at the optimal allowed
cadence of 0.1 Hz for the entire duration of the cruise phase.

  * •

The technological feasibility of reducing the Allan deviation with respect to
its baseline value of $\sim 10^{-14}$ at $10\,$s, as observed for the Cassini
mission data more than three decades ago, in 1992.

We maintain that these two points should become part of the mission
requirements, in light of the numerous scientific opportunities demonstrated
in this work. With them, the UOP mission has the potential for groundbreaking
discoveries about black holes living in galactic centers, the study of
particle physics beyond the standard model and the reach of current particle
accelerators, and the nature of dark matter.

###### Acknowledgements.

L.Z., D.S, D.J.D, D.O., A.D. and P.S. acknowledge the support of the
International Space Science Institute (ISSI) in Bern, Switzerland. L.Z. and
D.O. acknowledge support from ERC Starting Grant No. 121817–BlackHoleMergs.
D.J.D. and D.O. received funding from the Danish Independent Research Fund
through Sapere Aude Starting Grant No. 121587. A.D. acknowledges support from
NSF grant AST-2319441. We thank Mark Hofstadter, Joseph Lazio and Daniele
Durante for useful discussions. D.S. thanks Timothée Schaeffer for continuing
moral support and inspiration. L.Z. also acknowledges G.M.’s curiosity. The
research leading to these results has received funding from the Spanish
Ministry of Science and Innovation
(PID2020-115845GB-I00/AEI/10.13039/501100011033). IFAE is partially funded by
the CERCA program of the Generalitat de Catalunya. D.B. acknowledges the
support from the Departament de Recerca i Universitats de la Generalitat de
Catalunya al Grup de Recerca i Universitats from Generalitat de Catalunya to
the Grup de Recerca 00649 (Codi: 2021 SGR 00649). A.C.J. is supported by the
Science and Technology Facilities Council (STFC) through the UKRI Quantum
Technologies for Fundamental Physics Programme [grant number ST/T005904/1].
This work was partly enabled by the UCL Cosmoparticle Initiative.

## References

  * Abbott et al. [2022] Abbott, R., et al. 2022, Phys. Rev. D, 105, 063030, doi: 10.1103/PhysRevD.105.063030
  * Agazie et al. [2023a] Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJ, 951, L8, doi: 10.3847/2041-8213/acdac6
  * Agazie et al. [2023b] —. 2023b, arXiv e-prints, arXiv:2306.16220, doi: 10.48550/arXiv.2306.16220
  * Aghanim et al. [2020] Aghanim, N., et al. 2020, Astron. Astrophys., 641, A6, doi: 10.1051/0004-6361/201833910
  * Anderson et al. [1992] Anderson, J. D., Armstrong, J. W., Campbell, J. K., et al. 1992, Space Sci. Rev., 60, 591, doi: 10.1007/BF00216869
  * Antoniadis et al. [2023] Antoniadis, J., Arumugam, P., Arumugam, S., et al. 2023, arXiv e-prints, arXiv:2306.16214, doi: 10.48550/arXiv.2306.16214
  * Armstrong [2006] Armstrong, J. W. 2006, Living Reviews in Relativity, 9, 1, doi: 10.12942/lrr-2006-1
  * Armstrong et al. [1987] Armstrong, J. W., Estabrook, F. B., & Wahlquist, H. D. 1987, Astrophys. J. , 318, 536, doi: 10.1086/165390
  * Armstrong et al. [2003] Armstrong, J. W., Iess, L., Tortora, P., & Bertotti, B. 2003, Astrophys. J. , 599, 806, doi: 10.1086/379505
  * Bansal et al. [2023] Bansal, A., Ichiki, K., Tashiro, H., & Matsuoka, Y. 2023, MNRAS, 523, 3840, doi: 10.1093/mnras/stad1608
  * Bar et al. [2018] Bar, N., Blas, D., Blum, K., & Sibiryakov, S. 2018, Phys. Rev. D, 98, 083027, doi: 10.1103/PhysRevD.98.083027
  * Begelman et al. [1980] Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature (London), 287, 307, doi: 10.1038/287307a0
  * Bertotti et al. [1999] Bertotti, B., Vecchio, A., & Iess, L. 1999, Phys. Rev. D, 59, 082001, doi: 10.1103/PhysRevD.59.082001
  * Bertotti et al. [1992] Bertotti, B., Ambrosini, R., Asmar, S. W., et al. 1992, A&AS, 92, 431
  * Bi et al. [2024] Bi, X., Guo, Z., Zou, X., et al. 2024, arXiv e-prints, arXiv:2402.06096, doi: 10.48550/arXiv.2402.06096
  * Birnbaum et al. [2010] Birnbaum, K. M., Chen, Y., & Hemmati, H. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7587, Free-Space Laser Communication Technologies XXII, ed. H. Hemmati, 75870A, doi: 10.1117/12.843315
  * Blas et al. [2020a] Blas, D., Caputo, A., Ivanov, M. M., & Sberna, L. 2020a, Phys. Dark Univ., 27, 100428, doi: 10.1016/j.dark.2019.100428
  * Blas & Jenkins [2022a] Blas, D., & Jenkins, A. C. 2022a, Phys. Rev. D, 105, 064021, doi: 10.1103/PhysRevD.105.064021
  * Blas & Jenkins [2022b] —. 2022b, Phys. Rev. Lett., 128, 101103, doi: 10.1103/PhysRevLett.128.101103
  * Blas et al. [2020b] Blas, D., López Nacir, D., & Sibiryakov, S. 2020b, Phys. Rev. D, 101, 063016, doi: 10.1103/PhysRevD.101.063016
  * Blas et al. [2017] Blas, D., Nacir, D. L., & Sibiryakov, S. 2017, Phys. Rev. Lett. , 118, 261102, doi: 10.1103/PhysRevLett.118.261102
  * Bondi et al. [1959] Bondi, H., Pirani, F. A. E., & Robinson, I. 1959, Proceedings of the Royal Society of London Series A, 251, 519, doi: 10.1098/rspa.1959.0124
  * Buck et al. [2019] Buck, B. R., Allen, G. D., Duerr, E. K., et al. 2019, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10978, Advanced Photon Counting Techniques XIII, ed. M. A. Itzler, J. C. Bienfang, & K. A. McIntosh, 1097809, doi: 10.1117/12.2520617
  * Bucko et al. [2023] Bucko, J., Soyuer, D., & Zwick, L. 2023, MNRAS, 524, L32, doi: 10.1093/mnrasl/slad067
  * Caprini & Figueroa [2018] Caprini, C., & Figueroa, D. G. 2018, Classical and Quantum Gravity, 35, 163001, doi: 10.1088/1361-6382/aac608
  * Caprini et al. [2016] Caprini, C., et al. 2016, JCAP, 04, 001, doi: 10.1088/1475-7516/2016/04/001
  * Caprini et al. [2020] —. 2020, JCAP, 03, 024, doi: 10.1088/1475-7516/2020/03/024
  * Chadha-Day et al. [2022] Chadha-Day, F., Ellis, J., & Marsh, D. J. E. 2022, Sci. Adv., 8, abj3618, doi: 10.1126/sciadv.abj3618
  * Colpi et al. [2024] Colpi, M., Danzmann, K., Hewitson, M., et al. 2024, arXiv e-prints, arXiv:2402.07571, doi: 10.48550/arXiv.2402.07571
  * Comoretto et al. [1992] Comoretto, G., Bertotti, B., Iess, L., & Ambrosini, R. 1992, Nuovo Cimento C Geophysics Space Physics C, 15, 1193, doi: 10.1007/BF02506712
  * Croton [2009] Croton, D. J. 2009, MNRAS, 394, 1109, doi: 10.1111/j.1365-2966.2009.14429.x
  * Degnan [2002] Degnan, J. J. 2002, Journal of Geodynamics, 34, 551, doi: 10.1016/S0264-3707(02)00044-3
  * Degnan [2023] —. 2023, Photonics for Solar Energy Systems IX, 10, 98, doi: 10.3390/photonics10020098
  * Dirkx et al. [2019] Dirkx, D., Prochazka, I., Bauer, S., et al. 2019, Journal of Geodesy, 93, 2405, doi: 10.1007/s00190-018-1171-x
  * Dosopoulou & Antonini [2017] Dosopoulou, F., & Antonini, F. 2017, Astrophys. J. , 840, 31, doi: 10.3847/1538-4357/aa6b58
  * EPTA Collaboration et al. [2023] EPTA Collaboration, InPTA Collaboration, Antoniadis, J., et al. 2023, A&A, 678, A50, doi: 10.1051/0004-6361/202346844
  * Estabrook & Wahlquist [1975] Estabrook, F. B., & Wahlquist, H. D. 1975, General Relativity and Gravitation, 6, 439, doi: 10.1007/BF00762449
  * Fakhouri et al. [2010] Fakhouri, O., Ma, C.-P., & Boylan-Kolchin, M. 2010, MNRAS, 406, 2267, doi: 10.1111/j.1365-2966.2010.16859.x
  * Fedderke et al. [2021] Fedderke, M. A., Graham, P. W., & Rajendran, S. 2021, Phys. Rev. D, 103, 103017, doi: 10.1103/PhysRevD.103.103017
  * Fedderke & Mathur [2023] Fedderke, M. A., & Mathur, A. 2023, Phys. Rev. D, 107, 043004, doi: 10.1103/PhysRevD.107.043004
  * Ferreira [2021] Ferreira, E. G. M. 2021, Astron. Astrophys. Rev., 29, 7, doi: 10.1007/s00159-021-00135-6
  * Fletcher et al. [2020a] Fletcher, L. N., Simon, A. A., Hofstadter, M. D., et al. 2020a, arXiv e-prints, arXiv:2008.12125.  https://arxiv.org/abs/2008.12125
  * Fletcher et al. [2020b] Fletcher, L. N., Helled, R., Roussos, E., et al. 2020b, Planet. Space Sci., 191, 105030, doi: 10.1016/j.pss.2020.105030
  * Foreman-Mackey et al. [2013] Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Foster et al. [2018] Foster, J. W., Rodd, N. L., & Safdi, B. R. 2018, Phys. Rev. D, 97, 123006, doi: 10.1103/PhysRevD.97.123006
  * Foster & Backer [1990] Foster, R. S., & Backer, D. C. 1990, Astrophys. J. , 361, 300, doi: 10.1086/169195
  * Fukusumi et al. [2023] Fukusumi, K., Morisaki, S., & Suyama, T. 2023, Phys. Rev. D, 108, 095054, doi: 10.1103/PhysRevD.108.095054
  * Georgakakis et al. [2019] Georgakakis, A., Comparat, J., Merloni, A., et al. 2019, MNRAS, 487, 275, doi: 10.1093/mnras/sty3454
  * Girija [2023] Girija, A. P. 2023, Acta Astronautica, 202, 104, doi: 10.1016/j.actaastro.2022.10.005
  * Graham et al. [2015] Graham, P. W., Irastorza, I. G., Lamoreaux, S. K., Lindner, A., & van Bibber, K. A. 2015, Ann. Rev. Nucl. Part. Sci., 65, 485, doi: 10.1146/annurev-nucl-102014-022120
  * Helled & Fortney [2020] Helled, R., & Fortney, J. J. 2020, arXiv e-prints, arXiv:2007.10783.  https://arxiv.org/abs/2007.10783
  * Hindmarsh et al. [2014] Hindmarsh, M., Huber, S. J., Rummukainen, K., & Weir, D. J. 2014, Phys. Rev. Lett., 112, 041301, doi: 10.1103/PhysRevLett.112.041301
  * Hobbs et al. [2010] Hobbs, G., Archibald, A., Arzoumanian, Z., et al. 2010, Classical and Quantum Gravity, 27, 084013, doi: 10.1088/0264-9381/27/8/084013
  * Hofstadter et al. [2019] Hofstadter, M., Simon, A., Atreya, S., et al. 2019, Planet. Space Sci., 177, 104680, doi: 10.1016/j.pss.2019.06.004
  * Hopkins et al. [2005] Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2005, Astrophys. J. , 630, 705, doi: 10.1086/432438
  * Hu et al. [2000] Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett. , 85, 1158, doi: 10.1103/PhysRevLett.85.1158
  * Hui et al. [2017] Hui, L., Ostriker, J. P., Tremaine, S., & Witten, E. 2017, Phys. Rev. D, 95, 043541, doi: 10.1103/PhysRevD.95.043541
  * Iorio et al. [2023] Iorio, L., Girija, A. P., & Durante, D. 2023, MNRAS, 523, 3595, doi: 10.1093/mnras/stad1446
  * Jarmak et al. [2020] Jarmak, S., Leonard, E., Akins, A., et al. 2020, Acta Astronautica, 170, 6, doi: 10.1016/j.actaastro.2020.01.030
  * Kamionkowski et al. [1994] Kamionkowski, M., Kosowsky, A., & Turner, M. S. 1994, Phys. Rev. D, 49, 2837, doi: 10.1103/PhysRevD.49.2837
  * Kelley et. al [In Prep.] Kelley et. al. In Prep.
  * Khmelnitsky & Rubakov [2014] Khmelnitsky, A., & Rubakov, V. 2014, JCAP, 02, 019, doi: 10.1088/1475-7516/2014/02/019
  * Kollmann et al. [2020] Kollmann, P., Cohen, I., Allen, R. C., et al. 2020, Space Sci. Rev., 216, 78, doi: 10.1007/s11214-020-00696-5
  * Kormendy & Ho [2013] Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Kůs et al. [2024] Kůs, P., López Nacir, D., & Urban, F. R. 2024.  https://arxiv.org/abs/2402.04099
  * Leja et al. [2020] Leja, J., Speagle, J. S., Johnson, B. D., et al. 2020, Astrophys. J. , 893, 111, doi: 10.3847/1538-4357/ab7e27
  * Loeb [2024] Loeb, A. 2024.  https://arxiv.org/abs/2405.05410
  * Luo et al. [2016] Luo, J., Chen, L.-S., Duan, H.-Z., et al. 2016, Classical and Quantum Gravity, 33, 035010, doi: 10.1088/0264-9381/33/3/035010
  * Maggiore [2018] Maggiore, M. 2018, Gravitational Waves: Volume 2: Astrophysics and Cosmology, doi: 10.1093/oso/9780198570899.001.0001
  * Marasco et al. [2021] Marasco, A., Cresci, G., Posti, L., et al. 2021, MNRAS, 507, 4274, doi: 10.1093/mnras/stab2317
  * Miller et al. [2022] Miller, A. L., Badaracco, F., & Palomba, C. 2022, Phys. Rev. D, 105, 103035, doi: 10.1103/PhysRevD.105.103035
  * Misner et al. [1973] Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation
  * Molera Calvés et al. [2014] Molera Calvés, G., Pogrebenko, S. V., Cimò, G., et al. 2014, A&A, 564, A4, doi: 10.1051/0004-6361/201322925
  * Moore et al. [2015] Moore, C. J., Cole, R. H., & Berry, C. P. L. 2015, Classical and Quantum Gravity, 32, 015014, doi: 10.1088/0264-9381/32/1/015014
  * Müller et al. [2019] Müller, J., Murphy, T. W., Schreiber, U., et al. 2019, Journal of Geodesy, 93, 2195, doi: 10.1007/s00190-019-01296-0
  * Mumolo et al. [2021] Mumolo, J. M., Pestana, N., Clark, H., et al. 2021, Evaluation of asynchronous geiger-mode avalanche photodiode arrays for deep-space optical communications, DRAFT VERSION, Pasadena, CA: Jet Propulsion Laboratory, National Aeronautics and Space Administration, 2021, doi: 2014/54760
  * Nagano et al. [2021] Nagano, K., Nakatsuka, H., Morisaki, S., et al. 2021, Phys. Rev. D, 104, 062008, doi: 10.1103/PhysRevD.104.062008
  * Naoz et al. [2020] Naoz, S., Will, C. M., Ramirez-Ruiz, E., et al. 2020, Astrophys. J. Lett., 888, L8, doi: 10.3847/2041-8213/ab5e3b
  * Noda et al. [2023] Noda, H., Senshu, H., Otsubo, T., et al. 2023, Advances in Space Research, 71, 4196, doi: 10.1016/j.asr.2022.12.057
  * O’Hare [2020] O’Hare, C. 2020, cajohare/AxionLimits: AxionLimits, v1.0, https://cajohare.github.io/AxionLimits/, Zenodo, doi: 10.5281/zenodo.3932430
  * Oshima et al. [2023] Oshima, Y., Fujimoto, H., Kume, J., et al. 2023, Phys. Rev. D, 108, 072005, doi: 10.1103/PhysRevD.108.072005
  * Pirani [1956] Pirani, F. A. E. 1956, Acta Physica Polonica, 15, 389
  * Porayko et al. [2018] Porayko, N. K., et al. 2018, Phys. Rev. D, 98, 102002, doi: 10.1103/PhysRevD.98.102002
  * Pretorius et al. [2021] Pretorius, H., Aull, B., Clark, H., et al. 2021, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11721, Advanced Photon Counting Techniques XV, ed. M. A. Itzler, J. C. Bienfang, & K. A. McIntosh, 117210S, doi: 10.1117/12.2588469
  * Reardon et al. [2023] Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6, doi: 10.3847/2041-8213/acdd02
  * Ricarte & Natarajan [2018] Ricarte, A., & Natarajan, P. 2018, MNRAS, 474, 1995, doi: 10.1093/mnras/stx2851
  * Robson et al. [2019] Robson, T., Cornish, N. J., & Liu, C. 2019, Classical and Quantum Gravity, 36, 105011, doi: 10.1088/1361-6382/ab1101
  * Rodriguez-Gomez et al. [2015] Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, MNRAS, 449, 49, doi: 10.1093/mnras/stv264
  * Sereno & Jetzer [2006] Sereno, M., & Jetzer, P. 2006, MNRAS, 371, 626, doi: 10.1111/j.1365-2966.2006.10670.x
  * Sesana et al. [2008] Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192, doi: 10.1111/j.1365-2966.2008.13682.x
  * Sesana et al. [2021] Sesana, A., et al. 2021, Exper. Astron., 51, 1333, doi: 10.1007/s10686-021-09709-9
  * Shankar et al. [2013] Shankar, F., Weinberg, D. H., & Miralda-Escudé, J. 2013, MNRAS, 428, 421, doi: 10.1093/mnras/sts026
  * Simon et al. [2020] Simon, A. A., Fletcher, L. N., Arridge, C., et al. 2020, Space Sci. Rev., 216, 17, doi: 10.1007/s11214-020-0639-1
  * Smarra et al. [2023] Smarra, C., et al. 2023, Phys. Rev. Lett., 131, 171001, doi: 10.1103/PhysRevLett.131.171001
  * Smarra et al. [2024] —. 2024.  https://arxiv.org/abs/2405.01633
  * Sodnik et al. [2014] Sodnik, Z., Smit, H., Sans, M., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8971, Free-Space Laser Communication and Atmospheric Propagation XXVI, ed. H. Hemmati & D. M. Boroson, 89710W, doi: 10.1117/12.2045510
  * Soyuer et al. [2021] Soyuer, D., Zwick, L., D’Orazio, D. J., & Saha, P. 2021, MNRAS, 503, L73, doi: 10.1093/mnrasl/slab025
  * Spilker [2020] Spilker, T. R. 2020, Philosophical Transactions of the Royal Society of London Series A, 378, 20190488, doi: 10.1098/rsta.2019.0488
  * Springel et al. [2005] Springel, V., Di Matteo, T., & Hernquist, L. 2005, ApJ, 620, L79, doi: 10.1086/428772
  * The International Pulsar Timing Array Collaboration et al. [2023] The International Pulsar Timing Array Collaboration, Agazie, G., Antoniadis, J., et al. 2023, arXiv e-prints, arXiv:2309.00693, doi: 10.48550/arXiv.2309.00693
  * Unal et al. [2022] Unal, C., Urban, F. R., & Kovetz, E. D. 2022.  https://arxiv.org/abs/2209.02741
  * Vallisneri et al. [2020] Vallisneri, M., et al. 2020, doi: 10.3847/1538-4357/ab7b67
  * Vasiliev et al. [2015] Vasiliev, E., Antonini, F., & Merritt, D. 2015, Astrophys. J. , 810, 49, doi: 10.1088/0004-637X/810/1/49
  * Volonteri et al. [2003] Volonteri, M., Haardt, F., & Madau, P. 2003, Astrophys. J. , 582, 559, doi: 10.1086/344675
  * Wahlquist [1987] Wahlquist, H. 1987, General Relativity and Gravitation, 19, 1101, doi: 10.1007/BF00759146
  * Will et al. [2023] Will, C. M., Naoz, S., Hees, A., et al. 2023, Astrophys. J., 959, 58, doi: 10.3847/1538-4357/ad09b3
  * Williams et al. [1998] Williams, S., Bock, Y., & Fang, P. 1998, J. Geophys. Res., 103, 27,051, doi: 10.1029/98JB02794
  * Xu et al. [2023] Xu, H., Chen, S., Guo, Y., et al. 2023, Research in Astronomy and Astrophysics, 23, 075024, doi: 10.1088/1674-4527/acdfa5
  * Zwick et al. [2022] Zwick, L., Soyuer, D., & Bucko, J. 2022, A&A, 664, A188, doi: 10.1051/0004-6361/202243741