# Revisiting Stochastic Gravitational-wave Background in the Strong Signal
Case

Zheng-Cheng Liang  MOE Key Laboratory of TianQin Mission, TianQin Research
Center for Gravitational Physics $\&$ School of Physics and Astronomy,
Frontiers Science Center for TianQin, CNSA Research Center for Gravitational
Waves, Sun Yat-sen University (Zhuhai Campus), Zhuhai 519082, China  Zhi-Yuan
Li  MOE Key Laboratory of TianQin Mission, TianQin Research Center for
Gravitational Physics $\&$ School of Physics and Astronomy, Frontiers Science
Center for TianQin, CNSA Research Center for Gravitational Waves, Sun Yat-sen
University (Zhuhai Campus), Zhuhai 519082, China  En-Kun Li  MOE Key
Laboratory of TianQin Mission, TianQin Research Center for Gravitational
Physics $\&$ School of Physics and Astronomy, Frontiers Science Center for
TianQin, CNSA Research Center for Gravitational Waves, Sun Yat-sen University
(Zhuhai Campus), Zhuhai 519082, China  Jian-dong Zhang  MOE Key Laboratory of
TianQin Mission, TianQin Research Center for Gravitational Physics $\&$ School
of Physics and Astronomy, Frontiers Science Center for TianQin, CNSA Research
Center for Gravitational Waves, Sun Yat-sen University (Zhuhai Campus), Zhuhai
519082, China  Yi-Ming Hu<EMAIL_ADDRESS>MOE Key Laboratory of TianQin
Mission, TianQin Research Center for Gravitational Physics $\&$ School of
Physics and Astronomy, Frontiers Science Center for TianQin, CNSA Research
Center for Gravitational Waves, Sun Yat-sen University (Zhuhai Campus), Zhuhai
519082, China

###### Abstract

Weak-signal limit is often used in estimating stochastic gravitational-wave
background (SGWB) intensities. This approximation fails and the signal-to-
noise ratio (SNR) can be much weaker when background signals are loud compared
to the detector noise. In this work, we highlight this limitation for the SGWB
detection using space-borne detector networks. For the TianQin + LISA network,
the SNR estimated under the weak-signal limit might be off by as large as an
order of magnitude. Contour plots of SNR over the parameter spaces are also
presented to indicate regions susceptible to this discrepancy. Our results
suggest that DA and DB type extragalactic double white dwarfs may yield an
SGWB with SNR surpassing 100 after 1 year of operation in the weak-signal-
limit scenario, with a redshift-independent merger rate of about $500\,\,{\rm
Mpc^{-3}\,Myr^{-1}}$. In fact, this value falls significantly below the
necessary threshold. Similar influences arise for first-order phase
transitions, yet pinning down parameter regions remains formidable due to
model uncertainties.

###### pacs:

04.25.dg, 04.40.Nr, 04.70.-s, 04.70.Bw

##  I Introduction

The stochastic  background (SGWB) is a superposition of independent and
unresolvable GW from diverse sources, with the astrophysical component
including binary black hole (BBH), binary neutron star (BNS), double white
dwarf (DWD), etc Christensen (2019); Romano (2019). Alternatively, the
cosmological component traces back to early Universe processes, such as
inflation, first-order phase transitions, and cosmic defects de Araujo  _et
al._ (2000); Martinovic  _et al._ (2021). This diversity enables the SGWB to
cover multiple frequency ranges. A recent significant milestone was made by
collaborations utilizing Pulsar Timing Arrays (PTA), offering compelling
evidence for the existence of SGWB in the nanoHertz (nHz) band Agazie  _et
al._ (2023); Xu  _et al._ (2023); Antoniadis  _et al._ (2023); Reardon  _et
al._ (2023). While ground-based GW detectors have yet to detect the SGWB in
the hundred-hertz frequency band, the upper limit of SGWB intensity in this
range has been estimated based on available observation data Abbott  _et al._
(2021). In the future, with the potential advancements in space-borne GW
detectors, SGWB could be detected within the millihertz (mHz) frequency band
Luo  _et al._ (2016); Amaro-Seoane  _et al._ (2017).

Space-borne GW detectors, like TianQin Luo  _et al._ (2016), Laser
interferometry Space Antenna (LISA) Amaro-Seoane  _et al._ (2017), are
typically composed of three satellites arranged in an equilateral triangle
formation. However, the unequal armlengths resulting from detector motion
present a challenge in canceling laser noise, leading to the requirement of
time delay interferometry (TDI) technique Tinto and Armstrong (1999); Tinto
and Dhurandhar (2021). Among the various TDI combinations, the unequal-arm
Michelson with the noise-orthogonal basis (A,E,T) is commonly employed. The
A/E channels are specifically tailored for detecting GW, providing optimized
detection sensitivity Vallisneri  _et al._ (2008). The T channel exhibits a
substantially lower sensitivity to GW compared to the A/E channels, making it
more suitable for monitoring detector noise rather than efficiently detecting
GW Hogan and Bender (2001).

One often relies on the auto- or cross-correlation of the channel data to
derive the power spectral density (PSD) of SGWB. The null-channel method
utilizes the auto-correlation of the null channel to effectively track and
isolate noises, allowing for the identification of the SGWB component within a
single space-borne detector Tinto  _et al._ (2001); Hogan and Bender (2001);
Adams and Cornish (2014); Boileau  _et al._ (2021a); Muratore  _et al._
(2022); Cheng  _et al._ (2022). While in scenarios where a network of space-
borne detectors is deployed, each responding to a common SGWB but with
uncorrelated noises, distinguishing between the SGWB and detector noise can be
achieved through cross-correlating the channel data from two or more detectors
Hellings and Downs (1983); Michelson (1987); Christensen (1992); Flanagan
(1993); Allen and Romano (1999); Ungarelli and Vecchio (2001); Abbott  _et
al._ (2018); Hu  _et al._ (2024).

Following the null-channel and cross-correlation methods, the expectation of
correlation measurement is primarily determined by the PSD of SGWB and the
overlap reduction function (ORF), which quantifies the reduction in the
correlation to an SGWB resulting from the relative separation and orientation
of the channels. On the other hand, the variance of correlation measurement is
solely influenced by the noise PSD, given the SGWB is significantly lower than
the detector noise (weak-signal limit). Based on the expectation and variance,
the detection signal-to-noise ratio (SNR) can be calculated, taking into
account factors such as correlation time accumulation and frequency band
integration. Building upon the SNR for the weak-signal limit, the power-law
integrated sensitivity (PLIS) and peak-integrated sensitivity (PIS)
sensitivity curve is derived and commonly used to showcase the detection
capability for the SGWB Thrane and Romano (2013); Schmitz (2021).
Additionally, inversely solving the SNR can be employed to constrain the
detectable regions of source parameters, as exemplified by SNR contour plots
Caprini  _et al._ (2016); Kuroyanagi  _et al._ (2018); Caprini  _et al._
(2020).

While the weak-signal limit is applicable in most scenarios, there may exist
SGWB with intensities matching or exceeding the detector noise Amaro-Seoane
_et al._ (2012); Pan and Yang (2020); Huang  _et al._ (2020); Sharma and Harms
(2020); Zhou  _et al._ (2023); Wu and Li (2023); Torres-Orjuela  _et al._
(2023); Mentasti  _et al._ (2023); Song  _et al._ (2024). Therefore, Allen et
al. extended the SNR to handle arbitrarily large SGWB, specifically tailored
for an ORF-invariant channel pair Allen and Romano (1999); Cornish (2002);
Kudoh  _et al._ (2006). Their findings suggested that the existence of SGWBs
can impact the variance of correlation measurements, potentially resulting in
an overestimation of detection capability. This SNR estimator has been applied
in a series of subsequent studies Wang and Han (2021); Brzeminski  _et al._
(2022); Liang  _et al._ (2023a); Hu  _et al._ (2024). However, to the best of
our knowledge, it has not been utilized for constructing PLIS and PIS
sensitivity curves, nor in generating SNR contour plots.

In this paper, we further develop a more robust SNR estimation, adaptable to
multiple channel pairs and time-varying ORF. By incorporating the correction
of the SNR estimation, we reconstruct the PLIS and PIS sensitivity curves.
Focusing on TianQin, the TianQin I+II network, and the TianQin + LISA network,
we identify gaps in sensitivity curves between the scenarios of approximation
and non-approximation. Additionally, to illustrate the gaps in constraining
detectable parameter regions between these two scenarios, we employ the SNR
contour plot.

The structure of the paper is as follows. In Sec. II, we introduce the space
mission discussed in this study. Sec. III and IV are dedicated to SGWB
detection with and without the weak-signal limit. The case study is presented
in Sec. V. Finally, a brief conclusion and discussion is provided in Sec. VI.

##  II Space mission

Several space missions have been proposed for GW detection, such as TianQin
Luo  _et al._ (2016), LISA Amaro-Seoane  _et al._ (2017), among others. This
paper specifically centers on the TianQin and LISA missions.

###  II.1 detector design

As shown in Fig. 1, TianQin is a space-borne GW detector comprising three
identical satellites in Earth’s orbit. These satellites will have a 3.64-day
orbital period and will be positioned at a radius of approximately
$10^{5}\,\,{\rm km}$. When operational, they will form an equilateral triangle
with an arm length ($L_{\rm TQ}$) of about $1.7\times 10^{5}\,\,{\rm km}$.
TianQin will operate on a “three months on + three months off” schedule, with
a proposed second mission, TianQin II, aimed at filling detection gaps Ye  _et
al._ (2019). The orbital planes of TianQin and TianQin II are arranged
perpendicular to each other, creating the TianQin I+II network. The nominal
work scheme leads to no overlap in the operation times of TianQin and TianQin
II; we instead consider a scenario where the observation time can be extended
and correspond to an enhanced operation mode of “four months on + two months
off”.

LISA, on the other hand, is designed to orbit the Sun, trailing Earth by
approximately $20\degree$ and maintaining a fixed $60\degree$ angle with the
ecliptic plane. It will comprise three satellites separated by a distance of
about $2\times 10^{6}\,\,{\rm km}$, forming the armlength ($L_{\rm LISA}$) of
the detector. LISA will share the similar operating period and detection band
as TianQin, which opens up the opportunity to establish the TianQin + LISA
network, allowing for a half-a-year overlap within 1 year.

For convenience, shorthand notations can be used in figures and equations: TQ
for TianQin, TQ II for TianQin II, TT for the TianQin I+II network, and TL for
the TianQin + LISA network.

Figure 1: Schematic diagram of TianQin, TianQin II, and LISA in the ecliptic
plane.

###  II.2 TDI channel

Due to the detector motion, maintaining an equilateral shape for a space-borne
detector is challenging, making it difficult to implement an equal-arm
Michelson to cancel laser noise. In response to this challenge, the TDI
combination has been proposed Tinto and Armstrong (1999). The first-generation
TDI is initially designed for static detectors and includes components such as
the Sagnac, unequal-arm Michelson, Relay, Monitor, and Beacon Tinto and
Dhurandhar (2021). This paper focuses on the unequal-arm Michelson (X,Y,Z),
which can be established using each satellite and its adjacent link.
Furthermore, the corresponding noise-orthogonal unequal-arm Michelson can be
constructed by Vallisneri  _et al._ (2008)

$\displaystyle{\rm A}$ $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2}}({\rm
Z}-{\rm X}),$ $\displaystyle{\rm E}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{6}}({\rm X}-2{\rm Y}+{\rm Z}),$
$\displaystyle{\rm T}$ $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{3}}({\rm
X}+{\rm Y}+{\rm Z}).$ (1)

Following the successful cancel of laser noise, secondary noise, including
optical-path noise and acceleration noise, become more pronounced. Assuming
stationary noise conditions, the PSD for the TDI channel (A, E, T) can be
calculated by

$\displaystyle P_{\rm n_{\rm A/E}}(f)$ $\displaystyle=$
$\displaystyle\frac{2\sin^{2}\big{[}\frac{f}{f_{\ast}}\big{]}}{L^{2}}\bigg{[}\bigg{(}\cos\big{[}\frac{f}{f_{\ast}}\big{]}+2\bigg{)}S_{\rm
p}(f)+2\bigg{(}\cos\big{[}\frac{2f}{f_{\ast}}\big{]}+2\cos\big{[}\frac{f}{f_{\ast}}\big{]}+3\bigg{)}\frac{S_{\rm
a}(f)}{(2\pi f)^{4}}\bigg{]},$ $\displaystyle P_{\rm n_{T}}(f)$
$\displaystyle=$
$\displaystyle\frac{8\sin^{2}\big{[}\frac{f}{f_{\ast}}\big{]}\sin^{2}\big{[}\frac{f}{2f_{\ast}}\big{]}}{L^{2}}\bigg{(}S_{\rm
p}(f)+4\sin^{2}\big{[}\frac{f}{2f_{\ast}}\big{]}\frac{S_{\rm a}(f)}{(2\pi
f)^{4}}\bigg{)},$ (2)

where the characteristic frequency $f_{*}=c/(2\pi L)$, $S_{\rm p}$ and $S_{\rm
a}$ represent the PSDs of optical-path noise and acceleration noise,
respectively. For more comprehensive information on the parameters of TianQin
and LISA, readers are encouraged to refer to Refs. Mei  _et al._ (2020); Babak
_et al._ (2021). Unless otherwise stated, our study primarily focuses on the
case with using the A/E/T channels.

##  III Formalism

###  III.1 Stochastic background

In the transverse-traceless gauge, the metric perturbation of an SGWB can be
expressed as a superposition of plane wave:

$\displaystyle h(t,\vec{x})$ $\displaystyle=$
$\displaystyle\sum_{P=+,\times}\int_{-\infty}^{\infty}{\rm d}f\int_{S^{2}}{\rm
d}\hat{\Omega}_{\hat{k}}\,\widetilde{h}_{P}(f,\hat{k})\textbf{e}^{P}(\hat{k})$
(3) $\displaystyle\times e^{{\rm i}2\pi f[t-\hat{k}\cdot\vec{x}(t)/c]},$

where the wave vector of the gravitational wave is denoted by $\hat{k}$, the
polarization tensor $\textbf{e}^{P}(\hat{k})$ refers to the polarization $P$
of the GW, the speed of light is represented by $c$.

To begin, we constrain our focus on the Gaussian-stationary, unpolarized, and
isotropic SGWB. Therefore, it is reasonable to assume that the SGWB has a
zero-mean Fourier amplitude $\widetilde{h}_{P}(f,\hat{k})$, with the time-
independent variance:

$\langle\widetilde{h}_{P}(f,\hat{k})\widetilde{h}^{*}_{P^{\prime}}(f^{\prime},\hat{k}^{\prime})\rangle=\frac{1}{4}\delta(f-f^{\prime})\delta_{PP^{\prime}}\delta^{2}(\hat{k}-\hat{k}^{\prime})\mathscr{P}_{\rm
h}(|f|,\hat{k}).$ (4)

Here the factor $1/4$ stems from both the definition of one-side PSD and the
average of polarization. The PSD of SGWB incorporates contributions from all
directions in the sky, as reflected in the all-sky integral of the
$\mathscr{P}_{\rm h}$:

$\displaystyle S_{\rm h}(f)$ $\displaystyle=$ $\displaystyle\int_{S^{2}}{\rm
d}\hat{\Omega}_{\hat{k}}\mathscr{P}_{\rm h}(f,\hat{k})$ (5) $\displaystyle=$
$\displaystyle 4\pi\,\mathscr{P}_{\rm h}(f),$

where the second line of this equation is derived under the assumption of
isotropy in the SGWB.

To further characterize the distribution of energy across different
frequencies in the SGWB, a commonly used quantity is the dimensionless energy
spectrum density $\Omega_{\rm gw}$. This quantity serves to capture the ratio
of GW energy density ${\rm d}\rho_{\rm gw}$ within a specific frequency range
[$f$,$f+{\rm d}f$] to the critical energy density $\rho_{\rm gw}$, exhibiting
a direct relationship with $S_{\rm h}$ Allen (1996); Thrane and Romano (2013):

$\Omega_{\rm gw}(f)=\frac{1}{\rho_{\rm c}}\frac{{\rm d}\rho_{\rm gw}}{{\rm
d}(\ln{f})}=\frac{2\pi^{2}}{3H_{0}^{2}}f^{3}S_{\rm h}(f),$ (6)

where the critical energy density is defined as $\rho_{\rm
c}=3H_{0}^{2}c^{2}/(8\pi G)$, with the gravitational constant $G$, and the
Hubble constant $H_{0}$.

###  III.2 Overlap reduction function

The SGWB signal $h_{I}(t)$, observed in detector channel $I$, can be expressed
as the convolution of the metric perturbations $h(t,\vec{x})$ and the detector
response $\mathbb{D}^{ab}(t,\vec{x})$ Romano and Cornish (2017). To simplify
the measurement process without compromising accuracy, one can focus on a
narrow time scale, namely $[t_{0}-T/2,t_{0}+T/2]$, during which the assumption
can be made that the channel response remains unchanged. By employing the
short-term Fourier transform within this time interval, we have:

$\displaystyle h_{I}(t,t_{0})$ $\displaystyle=$
$\displaystyle\mathbb{D}_{I}[t,\vec{x}(t_{0})]*h[t,\vec{x}(t_{0})]$ (7)
$\displaystyle=$ $\displaystyle\sum_{P=+,\times}\int_{-\infty}^{\infty}{\rm
d}f\int_{S^{2}}{\rm
d}\hat{\Omega}_{\hat{k}}F_{I}^{P}(f,\hat{k},t_{0})\widetilde{h}_{P}(f,\hat{k})$
$\displaystyle\times e^{{\rm i}2\pi f[t-\hat{k}\cdot\vec{x}(t_{0})/c]},$

where the position vector of the measurement at time $t$ is denoted by
$\vec{x}$, the channel response is represented as the double contraction of
the channel tensor and the polarization tensor Cornish and Larson (2001).

In terms of the frequency domain SGWB signal

$\displaystyle\widetilde{h}_{I}(f,t_{0})$ $\displaystyle=$
$\displaystyle\sum_{P=+,\times}\int_{S^{2}}\,{\rm{d}}\hat{\Omega}_{\hat{k}}F_{I}^{P}(f,\hat{k},t_{0})\widetilde{h}_{P}(f,\hat{k})$
(8) $\displaystyle\times e^{-{\rm i}2\pi f\hat{k}\cdot\vec{x}(t_{0})/c},$

the PSD of SGWB signals is determined by

$\langle\widetilde{h}_{I}(f,t_{0})\widetilde{h}_{J}^{*}(f^{\prime},t_{0})\rangle=\frac{1}{2}\delta(f-f^{\prime})\Gamma_{IJ}(f,t_{0})S_{\rm
h}(|f|).$ (9)

Here, $I=J$ corresponds to the auto PSD of one channel, while $I\neq J$
indicates the cross PSD of two channels. Combined with Eq. (4), Eq. (5), Eq.
(8), and Eq. (9), the ORF

$\Gamma_{IJ}(f,t_{0})=\frac{1}{4\pi}\int_{S^{2}}{\rm
d}\hat{\Omega}_{\hat{k}}\mathcal{Y}_{IJ}(f,\hat{k},t_{0}),$ (10)

where the antenna pattern involves both the channel response and the
separation vector $\Delta\vec{x}=\vec{x}_{I}-\vec{x}_{J}$ between detectors:

$\displaystyle\mathcal{Y}_{IJ}(f,\hat{k},t_{0})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{P=+,\times}F^{P}_{I}(f,\hat{k},t_{0})F^{P*}_{J}(f,\hat{k},t_{0})$
(11) $\displaystyle\times e^{-{\rm i}2\pi
f\hat{k}\cdot[\vec{x}_{I}(t_{0})-\vec{x}_{J}(t_{0})]/c},$

Note that, for triangular-shaped detectors, it is convenient to designate one
of the vertices as the reference point for the detector position, facilitating
the determination of the separation vector between two detectors Liang  _et
al._ (2023a).

###  III.3 Detection sensitivity

The channel data $s_{I}$ is typically a combination of the SGWB signal $h_{I}$
and the channel noise $n_{I}$. As the SGWB signal and the channel noise are
expected to be uncorrelated, in order to extract the SGWB signal from the
data, a correlation measurement needs to be constructed as follows:

$S_{IJ}=\int_{t_{0}-T/2}^{t_{0}+T/2}{\rm d}t\int_{t_{0}-T/2}^{t_{0}+T/2}{\rm
d}t^{\prime}\left[s_{I}(t)s_{J}(t^{\prime})-\langle
n_{I}(t)n_{J}(t^{\prime})\rangle\right],$ (12)

where the noise PSD is determined by

$\langle
n_{I}(t)n_{J}(t^{\prime})\rangle=\frac{1}{2}\int_{-\infty}^{\infty}{\rm
d}f\,e^{-i2\pi f(t-t^{\prime})}P_{{\rm n}_{IJ}}(|f|).$ (13)

By monitoring the noise using the null-channel method Tinto  _et al._ (2001);
Hogan and Bender (2001), or achieving a zero-value noise PSD through the
cross-correlation method Christensen (1992); Flanagan (1993), it is possible
to identify the SGWB signal. These methods involve different channels from the
same and different detectors. To streamline subsequent analysis, we designate
the detectors as $\mathcal{A}$ and $\mathcal{B}$. Assuming that the employed
channels in each detector share the same auto-correlation ORF
$\mathcal{R}_{\mathcal{A}}$ and noise PSD, the total SNR in the weak-signal
limit Seto (2020); Liang  _et al._ (2022)

$\rho_{\rm tot}=\left[\frac{2\,T_{\rm
tot}}{1+\delta_{\mathcal{A}\mathcal{B}}}\int_{f_{\rm min}}^{f_{\rm max}}{\rm
d}f\,\frac{\bar{\Gamma}^{2}_{\rm tot}(f)\Omega_{\rm
gw}^{2}(f)}{\hat{D}_{\mathcal{A}\mathcal{B}}(f)}\right]^{1/2},$ (14)

where $T_{\rm tot}$ and $[f_{\rm min},f_{\rm max}]$ denote the total
correlation time and the detection band, respectively. The total time-averaged
ORF

$\bar{\Gamma}_{\rm tot}(f)=\sqrt{\frac{1}{T_{\rm tot}}\int_{0}^{T_{\rm
tot}}{\rm d}t_{0}\,\sum_{IJ}|\Gamma_{IJ}(f,t_{0})|^{2}}.$ (15)

The denominator of the integral term

$\hat{D}_{\mathcal{A}\mathcal{B}}(f)=\mathcal{R}_{\mathcal{A}}(f)\mathcal{R}_{\mathcal{B}}(f)\Omega_{{\rm
n}_{\mathcal{A}}}(f)\Omega_{{\rm n}_{\mathcal{B}}}(f),$ (16)

with the noise energy spectrum density

$\Omega_{{\rm
n}_{\mathcal{A}}}(f)=\frac{2\pi^{2}}{3H_{0}^{2}}f^{3}\frac{P_{{\rm
n}_{\mathcal{A}}}(f)}{\mathcal{R}_{\mathcal{A}}(f)}.$ (17)

Figure 2: Normalized ORF for TianQin (red), LISA (orange), the TianQin I+II
network (green), and the TianQin + LISA network (blue).

Figure 2 presents the ORF for TianQin, LISA, the TianQin I+II network, and the
TianQin + LISA network, employing the A/E channels. To account for the
frequency-dependent factor introduced by the TDI, the ORF is normalized by
$C_{\rm TDI}(f)=4\sin\left[2\pi fL/c\right]\sin\left[2\pi
fL^{\prime}/c\right]$, where $L=L^{\prime}$ for a single detector Cornish
(2002); Liang  _et al._ (2023b). It is observed that the total ORF for cross-
correlation between two detectors is much lower than that for the auto-
correlation of each respective detector channel within the mHz frequency band.
This effect becomes notably pronounced for the TianQin + LISA network, where
the large separation results in a considerable decrease in the drop frequency
of ORF. Furthermore, it is worth mentioning that for the A/E channels of the
same detector, the cross-correlation ORF is equal to 0.

The next step involves the analysis of the sensitivity curve. Considering an
SGWB characterized by the energy density spectrum that follows a power-law
distribution:

$\Omega_{\rm gw}(f)=\Omega_{0}(\epsilon)(f/f_{\rm
ref})^{\epsilon}|_{\epsilon=\epsilon_{0}},$ (18)

the SNR threshold $\rho_{\rm thr}$ can be translated into the threshold of
SGWB strength $\Omega_{0}$ based on Eq. (14):

$\Omega_{0}(\epsilon)=\rho_{\rm thr}\left[\frac{2\,T_{\rm
tot}}{1+\delta_{\mathcal{A}\mathcal{B}}}\int_{f_{\rm min}}^{f_{\rm max}}{\rm
d}f\frac{\bar{\Gamma}^{2}_{\rm tot}(f)(f/f_{\rm
ref})^{2\epsilon}}{\hat{D}_{\mathcal{A}\mathcal{B}}(f)}\right]^{-1/2},$ (19)

where the reference frequency $f_{\rm ref}$ is arbitrary. Subsequently, by
considering each frequency with a specific index $\epsilon$, the PLIS curve
can be generated by

$\Omega_{\rm PLIS}(f)={\rm max}_{\epsilon}[\Omega_{0}(\epsilon)(f/f_{\rm
ref})^{\epsilon}].$ (20)

This sensitivity curve, typically plotted on a logarithmic scale,
characterizes the collective envelopes of power-law SGWBs where the SNR aligns
with the preset threshold. Detection expectations for an SGWB hinge on its
power-law energy spectrum density in relation to the PLIS curve. If the
spectrum surpasses the PLIS curve, detection for the SGWB is anticipated;
however, if it falls below the curve, detection feasibility may be limited by
the existing detector configuration.

In addition to following a power-law form, the SGWB can also adhere to other
forms, such as the first-order PT, where the energy spectral density

$\Omega_{\rm gw}(f)=\tilde{\Omega}(\bm{\theta})\mathcal{S}(f,\tilde{f}).$ (21)

Here, $\tilde{\Omega}(\bm{\theta})$ signifies the peak energy spectral density
at the peak frequency $\tilde{f}$, with the spectral function
$\mathcal{S}(f,\tilde{f})$ being contingent upon the cosmological model. To
address this disparity, Schmitz et.al introduced the PIS curve Schmitz (2021):

$\Omega_{\rm PIS}(\tilde{f})=\left[\frac{2\,T_{\rm
tot}}{1+\delta_{\mathcal{A}\mathcal{B}}}\int_{f_{\rm min}}^{f_{\rm max}}{\rm
d}f\frac{\bar{\Gamma}^{2}_{IJ}(f)\mathcal{S}^{2}(f,\tilde{f})}{\hat{D}_{\mathcal{A}\mathcal{B}}(f)}\right]^{-1/2}.$
(22)

Returning to Eq. (14), the SNR can be determined by

$\rho=\frac{\tilde{\Omega}(\bm{\theta})}{\Omega_{\rm PIS}(\tilde{f})}.$ (23)

By selecting a particular set of parameters $\bm{\theta}$, the peak frequency
$\tilde{f}$ and peak energy spectrum density $\tilde{\Omega}(\bm{\theta})$ are
specified. If the peak amplitude exceeds $\rho_{\rm thr}$ times the PIS curve,
the detection of cosmological PT can be confirmed.

##  IV Beyond weak-signal limit

The weak-signal limit is applicable in most scenarios, especially for the
derivation of the sensitivity curves. However, in certain cases, the intensity
of a SGWB could become comparable to or even higher than the detector noise.
For instance, the foreground originating from Galactic double white dwarf
(GDWD) can exceed the noise level of space-borne detector noise Boileau  _et
al._ (2021b), and the foreground from BBHs can dominate over the noise of
next-generation ground-based detectors Regimbau  _et al._ (2017). In such
cases, adopting the SNR beyond the weak-signal limit becomes necessary.

Assuming that all other conditions specified in Eq. (12) hold true, except for
the weak-signal limit. Then for all channel pairs $\\{IJ\\}$ (or equivalently,
$\\{MN\\}$), the expectation and variance of the correlation measurement

$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\langle\sum_{IJ}S_{IJ}\rangle$ (24) $\displaystyle=$
$\displaystyle\sum_{IJ}\int_{0}^{T_{\rm tot}}{\rm
d}t_{0}\int_{-\infty}^{\infty}{\rm d}f\int_{-\infty}^{\infty}{\rm
d}f^{\prime}\,\langle\widetilde{h}_{I}(f,t_{0})\widetilde{h}^{*}_{J}(f^{\prime},t_{0})\rangle\widetilde{Q}_{IJ}(f^{\prime},t_{0})e^{-{\rm
i}2\pi(f-f^{\prime})t_{0}}$ $\displaystyle=$
$\displaystyle\sum_{IJ}\int_{0}^{T_{\rm tot}}{\rm d}t_{0}\int_{0}^{\infty}{\rm
d}f\int_{0}^{\infty}{\rm
d}f^{\prime}\,\delta(f-f^{\prime})\Gamma_{IJ}(f,t_{0})S_{\rm
h}(f)\widetilde{Q}_{IJ}(f^{\prime},t_{0})e^{-{\rm i}2\pi(f-f^{\prime})t_{0}}$
$\displaystyle=$
$\displaystyle\sum_{IJ}\frac{3H_{0}^{2}}{2\pi^{2}}\int_{0}^{T_{\rm tot}}{\rm
d}t_{0}\int_{0}^{\infty}{\rm d}f\,f^{-3}\Gamma_{IJ}(f,t_{0})\Omega_{\rm
gw}(f)\widetilde{Q}_{IJ}(f,t_{0}),$ $\displaystyle\sigma^{2}$ $\displaystyle=$
$\displaystyle\langle\mu^{2}\rangle-\langle\mu\rangle^{2}$ (25)
$\displaystyle=$ $\displaystyle\sum_{IJ,MN}\int_{0}^{T_{\rm tot}}{\rm
d}t_{0}\int_{0}^{T_{\rm tot}}{\rm d}\eta_{0}\int_{-\infty}^{\infty}{\rm
d}f\int_{-\infty}^{\infty}{\rm d}f^{\prime}\int_{-\infty}^{\infty}{\rm
d}\omega\int_{-\infty}^{\infty}{\rm d}\omega^{\prime}$
$\displaystyle\times\bigg{[}\left\langle\big{(}\widetilde{s}_{I}(f,t_{0})\widetilde{s}_{J}^{*}(f^{\prime},t_{0})-\langle\widetilde{n}_{I}(f,t_{0})\widetilde{n}_{J}^{*}(f^{\prime},t_{0})\rangle\big{)}\big{(}\widetilde{s}_{M}(\omega,\eta_{0})\widetilde{s}_{N}^{*}(\omega^{\prime},\eta_{0})-\langle\widetilde{n}_{M}(\omega,\eta_{0})\widetilde{n}_{N}^{*}(\omega^{\prime},\eta_{0})\rangle\big{)}\right\rangle$
$\displaystyle-\langle\widetilde{h}_{I}(f,t_{0})\widetilde{h}_{J}^{*}(f^{\prime},t_{0})\rangle\langle\widetilde{h}_{M}(\omega,\eta_{0})\widetilde{h}_{N}^{*}(\omega^{\prime},\eta_{0})\rangle\bigg{]}\widetilde{Q}_{IJ}(f^{\prime},t_{0})\widetilde{Q}^{*}_{MN}(\omega^{\prime},\eta_{0})e^{-{\rm
i}2\pi(f-f^{\prime})t_{0}}e^{-{\rm i}2\pi(\omega-\omega^{\prime})\eta_{0}}$
$\displaystyle=$ $\displaystyle 2\sum_{IJ,MN}\int_{0}^{T_{\rm tot}}{\rm
d}t_{0}\int_{0}^{\infty}{\rm
d}f\,\bigg{[}\langle\widetilde{s}_{I}(f,t_{0})\widetilde{s}_{M}(f,t_{0})\rangle\langle\widetilde{s}_{J}(f,t_{0})\widetilde{s}_{N}(f,t_{0})\rangle+\langle\widetilde{s}_{I}(f,t_{0})\widetilde{s}_{N}(f,t_{0})\rangle\langle\widetilde{s}_{J}(f,t_{0})\widetilde{s}_{M}(f,t_{0})\rangle\bigg{]}$
$\displaystyle\times\widetilde{Q}_{IJ}(f,t_{0})\widetilde{Q}^{*}_{MN}(f,t_{0})$
$\displaystyle\approx$
$\displaystyle\frac{1+\delta_{\mathcal{A}\mathcal{B}}}{2}\left[\frac{3H_{0}^{2}}{2\pi^{2}}\right]^{2}\sum_{IJ}\int_{0}^{T_{\rm
tot}}{\rm d}t_{0}\int_{0}^{\infty}{\rm
d}f\,f^{-6}D_{\mathcal{A}\mathcal{B}}(f)|\widetilde{Q}_{IJ}(f,t_{0})|^{2},$

where the filter function $\widetilde{Q}$ is arbitrary. Furthermore, for the
last line of Eq. (25) to be valid, it is necessary that the ORF for cross-
correlation is significantly lower than that for auto-correlation, leading to

$D_{\mathcal{A}\mathcal{B}}(f)=\mathcal{R}_{\mathcal{A}}(f)\mathcal{R}_{\mathcal{B}}(f)\left[\Omega_{{\rm
n}_{\mathcal{A}}}(f)+\Omega_{\rm gw}(f)\right]\left[\Omega_{{\rm
n}_{\mathcal{B}}}(f)+\Omega_{\rm gw}(f)\right].$ (26)

Combined with Eq. (24) and Eq. (25), the square of SNR can be obtained by

$\rho^{2}=\frac{\mu^{2}}{\sigma^{2}}=\frac{2}{1+\delta_{\mathcal{A}\mathcal{B}}}\frac{\big{[}\sum_{IJ}\int_{0}^{T_{\rm
tot}}{\rm d}t_{0}\int_{0}^{\infty}{\rm
d}f\,f^{-3}\Gamma_{IJ}(f,t_{0})\Omega_{\rm
gw}(f)\widetilde{Q}_{IJ}(f,t_{0})\big{]}^{2}}{\sum_{IJ}\int_{0}^{T_{\rm
tot}}{\rm d}t_{0}\int_{0}^{\infty}{\rm
d}f\,f^{-6}D_{\mathcal{A}\mathcal{B}}(f)|\widetilde{Q}_{IJ}(f,t_{0})|^{2}}.$
(27)

According to the definition of a positive-definite inner product, which is
applicable to any pair of complex functions $A$ and $B$ Allen and Romano
(1999):

$(A,B):=\sum_{IJ}\int_{0}^{T_{\rm tot}}{\rm d}t_{0}\int_{0}^{\infty}{\rm
d}f\,A_{IJ}(f,t_{0})B_{IJ}^{*}(f,t_{0})D_{\mathcal{A}\mathcal{B}}(f),$ (28)

Eq. (27) transitions to

$\rho_{\rm
opt}^{2}=\frac{2}{1+\delta_{\mathcal{A}\mathcal{B}}}\frac{\big{(}\widetilde{Q}_{IJ}(f,t_{0}),\frac{f^{3}\Gamma^{*}_{IJ}(f,t_{0})\Omega_{\rm
gw}(f)}{D_{IJ}(f)}\big{)}^{2}}{\big{(}\widetilde{Q}_{IJ}(f,t_{0}),\widetilde{Q}_{IJ}(f,t_{0})\big{)}}.$
(29)

Maximizing the SNR yields the solution

$\widetilde{Q}_{IJ}(f,t_{0})=\lambda\frac{f^{3}\Gamma^{*}_{IJ}(f,t_{0})\Omega_{\rm
gw}(f)}{D_{IJ}(f)},$ (30)

with a real constant $\lambda$. Following the preceding derivation, the total
SNR should be modified by

$\rho_{\rm tot}=\sqrt{\frac{2\,T_{\rm
tot}}{1+\delta_{\mathcal{A}\mathcal{B}}}\int_{f_{\rm min}}^{f_{\rm max}}{\rm
d}f\,\frac{\bar{\Gamma}^{2}_{\rm tot}(f)\Omega_{\rm
gw}^{2}(f)}{D_{\mathcal{A}\mathcal{B}}(f)}}.$ (31)

In contrast to Eq. (14), Eq. (31) suggests that the SGWB detection can suffer
negative impacts from the SGWB itself, especially in cases of loud SGWB. To
visualize this impact, plotting sensitivity curves is beneficial.

Although a direct derivation of $\Omega_{0}$ and $\Omega_{\rm PIS}$ directly
from Eq. (18) and Eq. (22) based on Eq. (31) to generate the sensitivity
curves proves challenging, these values can still be determined by inverting
the SNR formula111In this work, we utilize the Newton-Raphson method.. Fig. 3
and Fig. 4 illustrate the PLIS and PIS curves correspondingly, with solid
lines in red, green, and blue denoting TianQin, the TianQin I+II network, and
the TianQin + LISA network 222Note that for the TianQin detector, both auto-
correlations within the same channel and cross-correlations between different
channels are considered. However, in the case of detector networks, only the
cross-correlation between channels from different detectors is taken into
consideration., surpassing the weak-signal limit. Dashed lines including the
weak-signal limit are provided for comparison. Note that the spectral function
for the PIS curves is determined by

$S(f,\tilde{f}_{i})=\frac{3.9(f/\tilde{f}_{i})^{2.9}}{1+2.9(f/\tilde{f}_{i})^{3.8}}.$
(32)

Furthermore, the operation time is set at 1 year, resulting in a correlation
time $T_{\rm tot}$ equal to half a year, four months, and half a year for
TianQin, the TianQin I+II network, and the TianQin + LISA network,
respectively.

While sensitivity curves exhibit minimal gaps with and without the weak-signal
limit when the SNR remains below 10, a notable gap emerges as the SNR reaches
100, especially pronounced for the TianQin + LISA network, with gaps reaching
up to an order of magnitude. This suggests that the energy spectral density
predicted by the weak-signal limit needs to be raised by approximately an
order of magnitude for detectability.

Figure 3: PLIS curve with and without weak-signal limit. Figure 4: PIS curve
with and without weak-signal limit.

##  V Case study

In this section, we will briefly discuss the potential sources that could
generate loud SGWB in the mHz range, including DWDs and first-order PTs.

GDWDs are anticipated to produce foreground comparable to the noise of space-
borne detectors. In terms of our previous work, the PSD of Galactic foreground
can be fit for through the following polynomial Huang  _et al._ (2020); Liang
_et al._ (2022):

$S_{\rm h}(f)=\frac{20}{3}\left[10^{\sum_{i}a_{i}x^{i}}\right]^{2},$ (33)

where $x=\log(f/10^{-3}\,\,{\rm Hz})$, $a_{i}$ represents the polynomial
coefficients that are updated based on different operation times. The specific
values of these coefficients can be found in Table. 1. In terms of Eq. (6),
the energy spectrum density of Galactic foreground can be further derived.

Table 1: Coefficients for the polynomial fit for the Galactic foreground. Each successive row corresponds to an increasing operation time $T$. $T$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$ | $a_{6}$
---|---|---|---|---|---|---|---
0.5 yr | -18.7 | -1.23 | -0.801 | 0.832 | -1.96 | 3.09 | -2.38
1 yr | -18.7 | -1.34 | -0.513 | 0.0152 | -1.53 | 4.79 | -5.01
2 yr | -18.7 | -1.39 | -0.610 | 0.577 | 0.00242 | 0.578 | -4.39
4 yr | -18.7 | -1.30 | -0.872 | 0.266 | -5.12 | 15.6 | -15.5
5 yr | -18.7 | -1.32 | -0.322 | -1.68 | -4.49 | 21.6 | -22.6
Table 2: Ingredients contributing to the SGWB spectra arising from first-order PTs. The parameter $\alpha_{\infty}$ determines the threshold for the weakest transition at which a portion of the vacuum energy converts into kinetic energy, driving the expansion of bubbles Ellis  _et al._ (2019). The sound speed is given by $c_{\rm s}=1/\sqrt{3}$, the redshifted, present-day value of the Hubble frequency $h_{*}=1.65\times 10^{-5}\left(\frac{T_{*}}{100{\rm~{}GeV}}\right)\left(\frac{g_{*}}{100{\rm~{}GeV}}\right)^{1/6}$ Hz. $i$ | $p$ | $q$ | $\tilde{\Delta}_{i*}(v_{w})$ | $\frac{\tilde{f}_{i*}}{\beta}$ | $\mathcal{S}_{i}(f,\tilde{f}_{i})$ | $K_{i}(\alpha)$
---|---|---|---|---|---|---
col | 2 | 2 | $\frac{0.44v_{w}^{3}}{1+8.28v_{w}^{3}}$ | $\frac{0.31}{1-0.051v_{w}+0.88v_{w}^{2}}$ | $\left(\frac{f}{\tilde{f}_{i}}\right)^{2.8}\frac{3.8}{1+2.8(f/\tilde{f}_{i})^{3.8}}$ | ${\rm max}\big{[}1-\frac{\alpha_{\infty}}{\alpha},0\big{]}$
sw | 1 | 2 | $0.157v_{w}H_{*}\tau_{\rm sw}$ | $1.16(1-c_{\rm s}/v_{w})$ | $\left(\frac{f}{\tilde{f}_{i}}\right)^{3}\left[\frac{7}{4+3(f/\tilde{f}_{i})^{2}}\right]^{7/2}$ | $\left\\{\begin{gathered}\kappa(\alpha_{N})|_{\alpha_{N}=\alpha},\alpha\leq\alpha_{\infty}\\\ \frac{\alpha_{\infty}}{\alpha}\kappa(\alpha_{N})|_{\alpha_{N}=\alpha_{\infty}},\alpha>\alpha_{\infty}\end{gathered}\right.$
turb | 1 | $\frac{3}{2}$ | $20v_{w}(1-H_{*}\tau_{\rm sw})$ | $1.33(1-c_{\rm s}/v_{w})$ | $\left(\frac{f}{\tilde{f}_{i}}\right)^{3}\frac{[1+(f/\tilde{f}_{i})]^{3/11}}{1+8\pi f/h_{*}}$

Furthermore, it has been highlighted that the SGWB originating from
extragalactic double white dwarfs may exhibit significant intensity levels,
leading to high SNR Liang  _et al._ (2022). Without loss of generality, the
energy spectrum density of SGWB can be described as the cumulative
contribution from each binary system Abbott  _et al._ (2016, 2018):

$\Omega_{\rm gw}(f)=\frac{f}{\rho_{\rm c}}\int\,{\rm
d}\bm{\theta}\int_{0}^{z_{\rm max}}\,{\rm d}z\frac{R_{\rm m}(z)\frac{{\rm
d}E_{\rm gw}}{{\rm d}f_{\rm s}}(f_{\rm s},\bm{\theta})}{(1+z)H(z)},$ (34)

where the population parameter $\bm{\theta}$ primarily encompasses the
component masses of the system. In this paper, we focus on the DA and DB white
dwarfs, for which the mass distribution can be modeled by a Gaussian
distribution with mean $\mu=0.663\,\,M_{\odot}$ and standard deviation
$\sigma=0.177\,\,M_{\odot}$ Rosado (2011). The merger rate is directly linked
to the redshift and can be expressed as follows Abbott  _et al._ (2023):

$R_{\rm m}(z)=R_{0}(1+z)^{\kappa},$ (35)

with the local merger rate $R_{0}$ and index $\kappa$. As suggested in Ref.
Rosado (2011), the $R_{0}$ is assumed to be within the range of 20 to
$500\,\,{\rm Mpc^{-3}\,yr^{-1}}$. The energy spectrum emitted in the source
frame, denoted by ${\rm d}E_{\rm gw}/{\rm d}f_{\rm s}$, is characterized by
the frequency in the source frame $f_{\rm s}=f(1+z)$ Zhu  _et al._ (2013).
Given that the radius of a white dwarf significantly exceeds its Schwarzschild
radius, it is reasonable to assume that below the cut-off frequency where both
stars come into contact, the binary systems consistently remain in the
inspiral phase.

Figure 5: Noise energy spectrum density of TianQin, the TianQin I+II network,
the TianQin + LISA network, along with SGWB energy spectrum density of DWDs
and first-order PTs.

Regarding the first-order PTs, we adopt the model mentioned in our previous
study Liang  _et al._ (2022), where the energy spectrum density is formulated
as

$\Omega_{{\rm
gw}}(f)=\sum_{i}\tilde{\Omega}_{i}(\bm{\theta})\mathcal{S}_{i}(f,\tilde{f}_{i}).$
(36)

Here, the sum accounts for the collision occurring during the PT, as well as
the subsequent acoustic and turbulent phases. The redshifted peak frequency
today for each source is given by Jinno and Takimoto (2017)

$\tilde{f}_{i}=\frac{\tilde{f}_{i*}}{\beta}\left(\frac{\beta}{H_{*}}\right)h_{*},$
(37)

where the peak frequency before redshifting, denoted as $\tilde{f}_{i*}$,
depends on the sound speed $c_{\rm s}$ and the bubble velocity $v_{w}$. The
parameter $\beta/H_{*}$ represents PT duration normalized to the Hubble
parameter, while $h_{*}$ indicates the Hubble frequency. Additional details
are provided in Table. 2.

When GWs are generated with the Hubble rate $H_{*}$, temperature $T_{*}$, and
the total number of relativistic degrees of freedom $g_{*}$, the individual
contributions factoring in the redshift can be succinctly expressed as

$\tilde{\Omega}_{i}(\bm{\theta})=1.67\times
10^{-5}\left(\frac{100}{g_{*}(T_{*})}\right)^{1/3}\left(\frac{H_{*}}{\beta}\right)^{p}\left(K_{i}(\alpha)\frac{\alpha}{1+\alpha}\right)^{q}\tilde{\Delta}_{i*}(v_{w}),$
(38)

where the indices $p,q$, the spectral function
$\mathcal{S}_{i}(f,\tilde{f}_{i})$, the efficiency factor $K_{i}$, and the
peak amplitude $\tilde{\Delta}_{i*}$ are detailed in Table. 2. It is worth
noting that the peak amplitude involves the timescale for acoustic production
Hindmarsh  _et al._ (2017); Ellis  _et al._ (2019); Wang  _et al._ (2020)

$\displaystyle\tau_{\rm sw}$ $\displaystyle=$
$\displaystyle\min\left[{H^{-1}_{*}},R_{*}/U_{\rm f}\right],$ (39)

with the mean bubble separation $R_{*}$ and the average square root of the
fluid velocity $U_{\rm f}$ Espinosa  _et al._ (2010).

Table 3: Detection SNR for Galactic foreground. |  TQ |  TT |  TL
---|---|---|---
With limit |  300 |  120 |  430
Without limit |  170 |  69 |  96

Figure. 5 illustrates the SGWBs from various sources: the purple line
represents the Galactic foreground with a one-year operation time, the gray
region corresponds to the extragalactic background associated with the local
merger rate ranging from 20 to $500\,\,{\rm Mpc^{-3}\,yr^{-1}}$, and the
orange line signifies the cosmological background determined by the parameters
$\\{\alpha,\alpha_{\infty},v_{w},\beta/H_{*},T_{*}\\}=\\{2,1,0.95,50,10^{3}\,{\rm{GeV}}\\}$.
To facilitate comparison between the SGWB intensity and noise levels, the
noise energy spectrum density of TianQin, the TianQin I+II network, and the
TianQin + LISA network are shown in red, green, and blue lines, respectively.
While the extragalactic background can near the detector noise level, both the
Galactic foreground and cosmological background can exceed the noise. Taking
the Galactic foreground as an example, Table. 3 presents the detection SNR
within different detector configurations. It also provides the SNR under weak-
signal conditions for the sake of comparison. In the weak-signal limit, the
TianQin + LISA network achieves the highest SNR for detecting the Galactic
foreground. As moving beyond the weak-signal limit, the SNR drops to
approximately 20% for the TianQin + LISA network, whereas it remains around
60% for TianQin and the TianQin I+II network. Thus, TianQin demonstrates
superior performance in terms of detecting Galactic foreground.

The other two backgrounds exhibit considerable uncertainties in their
parameters compared to the Galactic foreground. We will proceed to utilize SNR
contour plots to determine the detectable parameter regions for them.

In Fig. 6, we present SNR contour plots for the extragalactic background,
showcasing scenarios with and without the weak-signal limit in the left and
right panels for comparison. Not surprisingly, an SNR threshold of 100 leads
to significant variations in the detectable parameter regions between the two
scenarios, especially for the TianQin + LISA network. Failure to achieve an
SNR above 100 within one year of operation would shift the lower limit for
$\kappa$ from around 0 to 2.5, indicating a transition from a scenario with
the weak-signal limit to one without it. This shift suggests that the adverse
effects of loud SGWB self-generation on detection reduce the region of
detectable parameters to a lesser extent than initially anticipated. On the
other hand, assuming a constant merger rate with redshift ($\kappa=0$),
achieving an SNR exceeding 100 after 1 year of operation, under the weak-
signal limit, would estimate the merger rate for EDWDs to be around
$500\,\,{\rm Mpc^{-3}\,yr^{-1}}$. However, this value falls considerably short
of meeting the actual detection requirements.

A similar diminishing effect of parameter constrains can be observed in Fig.
7, where we set the parameters
$\\{\alpha,\alpha_{\infty},v_{w}\\}=\\{2,1,0.95\\}$ for the first-order PTs.
If the SNR surpasses 100 after 1 year of operation, the temperature $T_{*}$ is
expected to be significantly greater than $10\,\,{\rm{GeV}}$, while the PT
duration $\beta/H_{*}$ should be less than $10^{2}$. By leveraging the weak-
signal limit, the detectable parameter region for these two variables can
significantly expand, potentially leading to misestimations. It is important
to highlight that, due to the wide range of variability in fixed parameters,
effectively constraining the parameters of PTs remains a challenging task.

Figure 6: SNR contour plots for the SGWB resulting from EDWD. The left-hand
and right hand panels refer to the scenarios with and without weak-signal
limit, considering a one-year operation time. Figure 7: SNR contour plots for
the SGWB resulting from first-order PTs. The left-hand and right hand panels
refer to the scenarios with and without weak-signal limit, considering a one-
year operation time.

##  VI Conclusions and discussions

In this paper, we have updated the SNR tailored for detecting loud SGWB. By
utilizing the improved SNR estimation, we have plotted PLIS and PIS curves for
TianQin, the TianQin I+II network, and TianQin + LISA network to assess the
detectable SGWB energy spectrum density over an operation time of 1 year. Our
findings indicated that for SNR values below 10, the impact of considering or
neglecting the weak-signal limit on the results is minimal. However, when the
SNR exceeds 100, applying the weak-signal limit can lead to a considerable
overestimation of the detectable energy spectrum density. This overestimation
becomes markedly pronounced for the TianQin + LISA network, where it may be
inflated by approximately an order of magnitude.

Additionally, we have employed SNR contour plots to determine the detectable
parameter regions for loud SGWBs generated by EDWDs and first-order PTs. Our
analysis revealed a significant discrepancy, especially concerning the TianQin
+ LISA network, when comparing results obtained with and without the weak-
signal limit using a selected SNR threshold of 100. Under the weak-signal
limit, the potential for detecting the SGWB from EDWDs of DA and DB types
appears promising within a one-year operation time, maintaining a consistent
merger rate of about $500\,\,{\rm Mpc^{-3}\,Myr^{-1}}$ across different
redshifts. However, this merger rate is notably below the required threshold
as the extragalactic SGWB is strong enough. Similarly, misestimations arise in
constraining the parameters of first-order PTs, where uncertainties in the PT
models only allow for broad constraints on the detectable parameter regions.

We emphasize that our analysis can extend to the statistical inference of SGWB
detection, which encompasses classical (frequentist) inference and Bayesian
inference. The pivotal aspect lies in constructing the likelihood function. By
estimating the expectation and variance in terms of Eq. (24) and Eq. (25), one
can ensure that the likelihood function aligns more closely with the actual
detection process, thereby minimizing bias in parameter estimation. This
approach is not limited to space-borne detectors Adams and Cornish (2014);
Caprini  _et al._ (2019); Pieroni and Barausse (2020); Flauger  _et al._
(2021), but is also relevant for next-generation ground-based detectors
Sachdev  _et al._ (2020); Zhou  _et al._ (2023); Zhong  _et al._ (2023); Pan
and Yang (2023); Bellie  _et al._ (2023); Song  _et al._ (2024). Ongoing work
on statistical inference regarding this matter is underway and will be
presented elsewhere.

###### Acknowledgements.

This work has been supported by the Guangdong Major Project of Basic and
Applied Basic Research (Grant No. 2019B030302001), the National Key Research
and Development Program of China (No. 2020YFC2201400), the Natural Science
Foundation of China (Grants No. 12173104), the Natural Science Foundation of
Guangdong Province of China (Grant No. 2022A1515011862), and the Guangdong
Basic and Applied Basic Research Foundation(Grant No. 2023A15150301). Z.C.L.
is supported by the China Postdoctoral Science Foundation (Grant No.
2023M744094). We also thank Fapeng Huang, Yun Jiang, and Jianwei Mei for
helpful discussions.

## References

  * Christensen (2019) N. Christensen, Rept. Prog. Phys. 82, 016903 (2019), arXiv:1811.08797 [gr-qc] .
  * Romano (2019) J. D. Romano (2019) arXiv:1909.00269 [gr-qc] .
  * de Araujo  _et al._ (2000) J. C. N. de Araujo, O. D. Miranda, and O. D. Aguiar, Phys. Rev. D 61, 124015 (2000), arXiv:astro-ph/0004395 .
  * Martinovic  _et al._ (2021) K. Martinovic, P. M. Meyers, M. Sakellariadou, and N. Christensen, Phys. Rev. D 103, 043023 (2021), arXiv:2011.05697 [gr-qc] .
  * Agazie  _et al._ (2023) G. Agazie _et al._ (NANOGrav), Astrophys. J. Lett. 951, L8 (2023), arXiv:2306.16213 [astro-ph.HE] .
  * Xu  _et al._ (2023) H. Xu _et al._ , Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023) J. Antoniadis _et al._ (EPTA, InPTA:), Astron. Astrophys. 678, A50 (2023), arXiv:2306.16214 [astro-ph.HE] .
  * Reardon  _et al._ (2023) D. J. Reardon _et al._ , Astrophys. J. Lett. 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE] .
  * Abbott  _et al._ (2021) R. Abbott _et al._ (KAGRA, Virgo, LIGO Scientific), Phys. Rev. D 104, 022004 (2021), arXiv:2101.12130 [gr-qc] .
  * Luo  _et al._ (2016) J. Luo _et al._ (TianQin), Class. Quant. Grav. 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Tinto and Armstrong (1999) M. Tinto and J. W. Armstrong, Phys. Rev. D 59, 102003 (1999).
  * Tinto and Dhurandhar (2021) M. Tinto and S. V. Dhurandhar, Living Rev. Rel. 24, 1 (2021).
  * Vallisneri  _et al._ (2008) M. Vallisneri, J. Crowder, and M. Tinto, Class. Quant. Grav. 25, 065005 (2008), arXiv:0710.4369 [gr-qc] .
  * Hogan and Bender (2001) C. J. Hogan and P. L. Bender, Phys. Rev. D 64, 062002 (2001), arXiv:astro-ph/0104266 .
  * Tinto  _et al._ (2001) M. Tinto, J. W. Armstrong, and F. B. Estabrook, Phys. Rev. D 63, 021101 (2001).
  * Adams and Cornish (2014) M. R. Adams and N. J. Cornish, Phys. Rev. D 89, 022001 (2014), arXiv:1307.4116 [gr-qc] .
  * Boileau  _et al._ (2021a) G. Boileau, N. Christensen, R. Meyer, and N. J. Cornish, Phys. Rev. D 103, 103529 (2021a), arXiv:2011.05055 [gr-qc] .
  * Muratore  _et al._ (2022) M. Muratore, D. Vetrugno, S. Vitale, and O. Hartwig, Phys. Rev. D 105, 023009 (2022), arXiv:2108.02738 [gr-qc] .
  * Cheng  _et al._ (2022) J. Cheng, E.-K. Li, Y.-M. Hu, Z.-C. Liang, J.-d. Zhang, and J. Mei, Phys. Rev. D 106, 124027 (2022), arXiv:2208.11615 [gr-qc] .
  * Hellings and Downs (1983) R. Hellings and G. Downs, Astrophys. J. Lett. 265, L39 (1983).
  * Michelson (1987) P. F. Michelson, Monthly Notices of the Royal Astronomical Society 227, 933 (1987).
  * Christensen (1992) N. Christensen, Phys. Rev. D46, 5250 (1992).
  * Flanagan (1993) E. E. Flanagan, Phys. Rev. D48, 2389 (1993), arXiv:astro-ph/9305029 [astro-ph] .
  * Allen and Romano (1999) B. Allen and J. D. Romano, Phys. Rev. D59, 102001 (1999), arXiv:gr-qc/9710117 [gr-qc] .
  * Ungarelli and Vecchio (2001) C. Ungarelli and A. Vecchio, Phys. Rev. D 64, 121501 (2001), arXiv:astro-ph/0106538 .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 120, 091101 (2018), arXiv:1710.05837 [gr-qc] .
  * Hu  _et al._ (2024) Y. Hu, P.-P. Wang, Y.-J. Tan, and C.-G. Shao, Astrophys. J. 961, 116 (2024).
  * Thrane and Romano (2013) E. Thrane and J. D. Romano, Phys. Rev. D 88, 124032 (2013), arXiv:1310.5300 [astro-ph.IM] .
  * Schmitz (2021) K. Schmitz, JHEP 01, 097 (2021), arXiv:2002.04615 [hep-ph] .
  * Caprini  _et al._ (2016) C. Caprini _et al._ , JCAP 04, 001 (2016), arXiv:1512.06239 [astro-ph.CO] .
  * Kuroyanagi  _et al._ (2018) S. Kuroyanagi, T. Chiba, and T. Takahashi, JCAP 11, 038 (2018), arXiv:1807.00786 [astro-ph.CO] .
  * Caprini  _et al._ (2020) C. Caprini _et al._ , JCAP 03, 024 (2020), arXiv:1910.13125 [astro-ph.CO] .
  * Amaro-Seoane  _et al._ (2012) P. Amaro-Seoane _et al._ , Class. Quant. Grav. 29, 124016 (2012), arXiv:1202.0839 [gr-qc] .
  * Pan and Yang (2020) Z. Pan and H. Yang, Class. Quant. Grav. 37, 195020 (2020), arXiv:1910.09637 [astro-ph.CO] .
  * Huang  _et al._ (2020) S.-J. Huang, Y.-M. Hu, V. Korol, P.-C. Li, Z.-C. Liang, Y. Lu, H.-T. Wang, S. Yu, and J. Mei, Phys. Rev. D 102, 063021 (2020), arXiv:2005.07889 [astro-ph.HE] .
  * Sharma and Harms (2020) A. Sharma and J. Harms, Phys. Rev. D 102, 063009 (2020), arXiv:2006.16116 [gr-qc] .
  * Zhou  _et al._ (2023) B. Zhou, L. Reali, E. Berti, M. Çalışkan, C. Creque-Sarbinowski, M. Kamionkowski, and B. S. Sathyaprakash, Phys. Rev. D 108, 064040 (2023), arXiv:2209.01310 [gr-qc] .
  * Wu and Li (2023) J. Wu and J. Li, Phys. Rev. D 108, 124047 (2023), arXiv:2307.05568 [gr-qc] .
  * Torres-Orjuela  _et al._ (2023) A. Torres-Orjuela, S.-J. Huang, Z.-C. Liang, S. Liu, H.-T. Wang, C.-Q. Ye, Y.-M. Hu, and J. Mei, (2023), arXiv:2307.16628 [gr-qc] .
  * Mentasti  _et al._ (2023) G. Mentasti, C. R. Contaldi, and M. Peloso, (2023), arXiv:2312.10792 [gr-qc] .
  * Song  _et al._ (2024) H. Song, D. Liang, Z. Wang, and L. Shao, (2024), arXiv:2401.00984 [gr-qc] .
  * Cornish (2002) N. J. Cornish, Phys. Rev. D 65, 022004 (2002), arXiv:gr-qc/0106058 .
  * Kudoh  _et al._ (2006) H. Kudoh, A. Taruya, T. Hiramatsu, and Y. Himemoto, Phys. Rev. D 73, 064006 (2006), arXiv:gr-qc/0511145 .
  * Wang and Han (2021) G. Wang and W.-B. Han, Phys. Rev. D 104, 104015 (2021), arXiv:2108.11151 [gr-qc] .
  * Brzeminski  _et al._ (2022) D. Brzeminski, A. Hook, and G. Marques-Tavares, JHEP 11, 061 (2022), arXiv:2203.13842 [hep-ph] .
  * Liang  _et al._ (2023a) Z.-C. Liang, Z.-Y. Li, J. Cheng, E.-K. Li, J.-d. Zhang, and Y.-M. Hu, Phys. Rev. D 107, 083033 (2023a), arXiv:2212.02852 [astro-ph.GA] .
  * Ye  _et al._ (2019) B.-B. Ye, X. Zhang, M.-Y. Zhou, Y. Wang, H.-M. Yuan, D. Gu, Y. Ding, J. Zhang, J. Mei, and J. Luo, Int. J. Mod. Phys. D 28, 09 (2019), arXiv:2012.03260 [gr-qc] .
  * Mei  _et al._ (2020) J. Mei _et al._ (TianQin),  (2020), 10.1093/ptep/ptaa114, arXiv:2008.10332 [gr-qc] .
  * Babak  _et al._ (2021) S. Babak, A. Petiteau, and M. Hewitson, (2021), arXiv:2108.01167 [astro-ph.IM] .
  * Allen (1996) B. Allen, in  _Les Houches School of Physics: Astrophysical Sources of Gravitational Radiation_ (1996) pp. 373–417, arXiv:gr-qc/9604033 .
  * Romano and Cornish (2017) J. D. Romano and N. J. Cornish, Living Rev. Rel. 20, 2 (2017), arXiv:1608.06889 [gr-qc] .
  * Cornish and Larson (2001) N. J. Cornish and S. L. Larson, Class. Quant. Grav. 18, 3473 (2001), arXiv:gr-qc/0103075 .
  * Seto (2020) N. Seto, Phys. Rev. D 102, 123547 (2020), arXiv:2010.06877 [gr-qc] .
  * Liang  _et al._ (2022) Z.-C. Liang, Y.-M. Hu, Y. Jiang, J. Cheng, J.-d. Zhang, and J. Mei, Phys. Rev. D 105, 022001 (2022), arXiv:2107.08643 [astro-ph.CO] .
  * Liang  _et al._ (2023b) Z.-C. Liang, Z.-Y. Li, E.-K. Li, J.-d. Zhang, and Y.-M. Hu, (2023b), arXiv:2307.01541 [gr-qc] .
  * Boileau  _et al._ (2021b) G. Boileau, A. Lamberts, N. J. Cornish, and R. Meyer, Mon. Not. Roy. Astron. Soc. 508, 803 (2021b), [Erratum: Mon.Not.Roy.Astron.Soc. 508, 5554–5555 (2021)], arXiv:2105.04283 [gr-qc] .
  * Regimbau  _et al._ (2017) T. Regimbau, M. Evans, N. Christensen, E. Katsavounidis, B. Sathyaprakash, and S. Vitale, Phys. Rev. Lett. 118, 151105 (2017), arXiv:1611.08943 [astro-ph.CO] .
  * Ellis  _et al._ (2019) J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen, JCAP 06, 024 (2019), arXiv:1903.09642 [hep-ph] .
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 131102 (2016), arXiv:1602.03847 [gr-qc] .
  * Rosado (2011) P. A. Rosado, Phys. Rev. D 84, 084004 (2011), arXiv:1106.5795 [gr-qc] .
  * Abbott  _et al._ (2023) R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), Phys. Rev. X 13, 011048 (2023), arXiv:2111.03634 [astro-ph.HE] .
  * Zhu  _et al._ (2013) X.-J. Zhu, E. J. Howell, D. G. Blair, and Z.-H. Zhu, Mon. Not. Roy. Astron. Soc. 431, 882 (2013), arXiv:1209.0595 [gr-qc] .
  * Jinno and Takimoto (2017) R. Jinno and M. Takimoto, Phys. Rev. D 95, 024009 (2017), arXiv:1605.01403 [astro-ph.CO] .
  * Hindmarsh  _et al._ (2017) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. D 96, 103520 (2017), [Erratum: Phys.Rev.D 101, 089902 (2020)], arXiv:1704.05871 [astro-ph.CO] .
  * Wang  _et al._ (2020) X. Wang, F. P. Huang, and X. Zhang, JCAP 05, 045 (2020), arXiv:2003.08892 [hep-ph] .
  * Espinosa  _et al._ (2010) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, JCAP 06, 028 (2010), arXiv:1004.4187 [hep-ph] .
  * Caprini  _et al._ (2019) C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone, and G. Tasinato, JCAP 11, 017 (2019), arXiv:1906.09244 [astro-ph.CO] .
  * Pieroni and Barausse (2020) M. Pieroni and E. Barausse, JCAP 07, 021 (2020), [Erratum: JCAP 09, E01 (2020)], arXiv:2004.01135 [astro-ph.CO] .
  * Flauger  _et al._ (2021) R. Flauger, N. Karnesis, G. Nardini, M. Pieroni, A. Ricciardone, and J. Torrado, JCAP 01, 059 (2021), arXiv:2009.11845 [astro-ph.CO] .
  * Sachdev  _et al._ (2020) S. Sachdev, T. Regimbau, and B. S. Sathyaprakash, Phys. Rev. D 102, 024051 (2020), arXiv:2002.05365 [gr-qc] .
  * Zhong  _et al._ (2023) H. Zhong, R. Ormiston, and V. Mandic, Phys. Rev. D 107, 064048 (2023), [Erratum: Phys.Rev.D 108, 089902 (2023)], arXiv:2209.11877 [gr-qc] .
  * Pan and Yang (2023) Z. Pan and H. Yang, Phys. Rev. D 107, 123036 (2023), arXiv:2301.04529 [gr-qc] .
  * Bellie  _et al._ (2023) D. S. Bellie, S. Banagiri, Z. Doctor, and V. Kalogera, (2023), arXiv:2310.02517 [gr-qc] .

  *[]:
  *[SGWB]: stochastic   () background
  *[GW]: gravitational-wave
  *[BBH]: binary black hole
  *[BNS]: binary neutron star
  *[DWD]: double white dwarf
  *[PTA]: Pulsar Timing Arrays
  *[LISA]: Laser interferometry Space Antenna
  *[TDI]: time delay interferometry
  *[PSD]: power spectral density
  *[ORF]: overlap reduction function
  *[SNR]: signal-to-noise ratio
  *[PLIS]: power-law integrated sensitivity
  *[PIS]: peak-integrated sensitivity
  *[SGWBs]: stochastic   () background
  *[PSDs]: power spectral density
  *[PT]: phase transition
  *[GDWD]: Galactic double white dwarf
  *[BBHs]: binary black hole
  *[DWDs]: double white dwarf
  *[PTs]: phase transition
  *[GDWDs]: Galactic double white dwarf
  *[GWs]: gravitational-wave
  *[EDWDs]: extragalactic double white dwarf
  *[EDWD]: extragalactic double white dwarf


1920

# Dark photon limits from magnetic fields and astrophysical plasmas

Giacomo Marocco
Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Parks
Road,
Oxford OX1 3PU, United Kingdom

###### Abstract

Dark photons may kinetically mix with our photon and modify Maxwell’s
equations. We examine their impact on the Lorentz force acting on the plasma
in astrophysical magnetospheres. Solving the relevant magnetohydrodynamic
equations, we show that any dark photon effect is negligible. The existing
bound on dark photons with a mass below $10^{-14}$ eV, coming from Pioneer 10
observations of Jupiter’s magnetosphere, relies on such an effect and is thus
unreliable. To rectify this, we use recent Juno analyses to isolate any dark
photon-induced modification of the Jovian magnetic field and hence place a
numerically similar limit.

##  I Introduction

There may exist a light gauge boson called the dark photon (DP) $X_{\mu}$ that
kinetically mixes with our photon $A_{\mu}$ Holdom (1986). The Standard Moded
field content is uncharged under the corresponding gauge group, making such a
boson difficult to detect as it mediates no long-range force, even if it is
very light.

The defining Lagrangian of the DP may be written as

$\begin{split}\mathcal{L}&\supset-\frac{1}{4}\left(F_{\mu\nu}F^{\mu\nu}+X_{\mu\nu}X^{\mu\nu}\right)\\\
&+\frac{1}{2}m^{2}X^{\mu}X_{\mu}+J^{\mu}_{\mathrm{EM}}A_{\mu}+\frac{\sin\alpha}{2}F^{\mu\nu}X_{\mu\nu},\end{split}$
(1)

where $X_{\mu\nu}=\partial_{\mu}X_{\nu}-\partial_{\nu}X_{\mu}$, $\sin\alpha$
parametrises the kinetic mixing, and $m$ is related to the DP mass. We remain
agnostic as to the necessary UV completion of such a model, limiting ourselves
to low energy scales where this approach is valid. To get to a basis with
diagonal kinetic terms, we let $A_{\mu}\rightarrow\frac{A_{\mu}}{\cos\alpha}$
and $X_{\mu}\rightarrow X_{\mu}+\sin\alpha A_{\mu}$, so that the Lagrangian
reads

$\begin{split}\mathcal{L}&\supset-\frac{1}{4}\left(F_{\mu\nu}F^{\mu\nu}+X_{\mu\nu}X^{\mu\nu}\right)\\\
&+\tilde{e}J^{\mu}_{\mathrm{EM}}A_{\mu}+\frac{m^{2}}{2}\left(X^{\mu}X_{\mu}+2\chi
X^{\mu}A_{\mu}+\chi^{2}A^{\mu}A_{\mu}\right),\end{split}$ (2)

where $\chi\equiv\tan\alpha$ and $\tilde{e}\equiv e(1+\chi^{2})^{1/2}$. The
Lagrangian still possesses the $U(1)$ gauge invariance of electromagnetism,
despite appearances, in contrast to the Proca Lagrangian. The mass matrix is
singular, containing a zero eigenvalue as dictated by gauge invariance; the
other eigenvalue, the mass of the DP, is $m_{X}^{2}=(1+\chi^{2})m^{2}$. The
model parameter space is then simply spanned by $\\{m_{X},\epsilon\\}$.

There exist many constraints on the DP parameter space, see, e.g., Jaeckel and
Ringwald (2010); Fabbrichesi  _et al._ (2020) for a review of these. For
ultralight DPs (below around ${10}^{-10}\text{\,}\mathrm{e}\mathrm{V}$), there
are very strong constraints if the DP constitutes dark matter. If this is not
the case, the bounds are weaker. A constraint comes from a Cavendish-type
experiment Williams  _et al._ (1971) exploring deviations from Coulomb’s law
Bartlett and Lögl (1988). The lack of spectral distortions measured by COBE
Fixsen  _et al._ (1996) in the CMB places cosmological constraints Caputo  _et
al._ (2020); García  _et al._ (2020). Observations of the magnetic fields
around the Earth Goldhaber and Nieto (1971) also place limits Bartlett and
Lögl (1988); Fischbach  _et al._ (1994); Kloor  _et al._ (1994). Finally,
Kloor et al. Kloor  _et al._ (1994) claim that measurements of Jupiter’s
magnetosphere, as reported in Smith  _et al._ (1974a), may bound dark photons
at still lower masses.

We critically examine Kloor et al.’s recasting of the photon mass bound coming
from Jupiter’s magnetic field in terms of the DP model. We highlight the
distinction between bounds that come from the observation of the forces on the
plasma in a magnetosphere, and those that come from the the measured magnitude
of the magnetic field ascribed to a dipole source. The former, in the
particular case of Jupiter, has been used to place a robust bound on the mass
of the photon Davis  _et al._ (1975), but does not constrain DPs, as we will
show. Furthermore, the latter is not enough to place constraints without
including the external (non-intrinsic) contributions to the Jovian magnetic
field, which formerly were not taken into account.

The secure bound on the photon mass is derived from the non-observation of the
effect of a massive photon on the plasma in the Jovian magnetosphere. As any
photon mass increases, the size of the Lorentz force exerted by the magnetic
field on a plasma increases, and at some point is in conflict with the
observed flow. We derive the DP-altered magnetohydrodynamic (MHD) equations
and show that no such effect occurs and hence no recasting of the photon mass
bound is possible. We thus turn to measurements from the Juno mission Bolton
_et al._ (2017) to obtain the size of any anomalous magnetic field that cannot
be ascribed to the Jovian magnetosphere. This allows a reliable bound to be
placed on very low mass DPs.

Incidentaly, we demonstrate that the most stringent bounds on the photon mass
Ryutov (2007) from observations of the approximate equilibrium of the solar
wind around Pluto Ness and Burlaga (2001); Burlaga  _et al._ (2003) do not
constrain DPs either.

The layout of the paper is as follows. In section II, we present the modified
form of Maxwell’s equations, as well as the MHD equations in the quasistatic
limit. In section III, we present the bounds coming from careful observation
of the Earth’s magnetic field. In section IV, we review the knowledge of
Jupiter’s magnetosphere that allows us to obtain an upper limit on the size of
external magnetic fields around Jupiter, and thus to set a DP limit. In
section V, we solve the relevant MHD equations to show that the forces on
plasmas in the magnetosphere are largely unaffected by any DP, and hence their
behaviour cannot be used to constrain DPs. We show in section VI that the same
conclusions are reached regarding the observation of currents in the far
heliosphere around Pluto. Finally, we conlude in section VII. Appendix A
contains a derivation of the magnetic field created by a point-like MDM when
Maxwell’s equations are modified by the inclusion of a DP.

##  II DP-modifed Maxwell and MHD equations

In this section, we derive the equations governing the magnetic fields and
currents within magnetospheres. In order to make contact with the common form
of Maxwell’s equations, we will work in non-relativistic notation. We write
the regular electric and magnetic fields as $\mathbf{E}$ and $\mathbf{B}$,
respectively, which may be written in terms of the scalar potential $\phi$ and
vector potential $\mathbf{A}$ as

$\displaystyle\mathbf{E}$ $\displaystyle=-\nabla\phi-\partial_{t}\mathbf{A},$
(3a) $\displaystyle\mathbf{B}$ $\displaystyle=\nabla\times\mathbf{A}.$ (3b)

Similarly, we write the dark electric and magnetic fields as
$\boldsymbol{\eta}$ and $\boldsymbol{\beta}$, respectively. These are given in
terms of their potentials as

$\displaystyle\boldsymbol{\eta}$
$\displaystyle=-\nabla\psi-\partial_{t}\mathbf{X},$ (4a)
$\displaystyle\boldsymbol{\beta}$ $\displaystyle=\nabla\times\mathbf{X}.$ (4b)

From the Lagrangian (2), we may derive the dark photon-modified Maxwell’s
equations

$\displaystyle\nabla\cdot\mathbf{E}+\chi^{2}m^{2}\phi$
$\displaystyle=\tilde{e}\rho+\chi m^{2}\psi,$ (5a)
$\displaystyle\nabla\times\mathbf{B}+\chi^{2}m^{2}\mathbf{A}$
$\displaystyle=\tilde{e}\mathbf{J}+\partial_{t}\mathbf{E}+\chi
m^{2}\mathbf{X},$ (5b) $\displaystyle\nabla\times\mathbf{E}$
$\displaystyle=-\partial_{t}\mathbf{B},$ (5c)
$\displaystyle\nabla\cdot\mathbf{B}$ $\displaystyle=0.$ (5d)

Here, $\rho$ is the charge density and $\mathbf{J}$ is the charge current,
both in units of $\tilde{e}$. Note that the homogeneous equations remain
unchanged compared to the case with no dark photon. The dark photon also
satisfies its set of Dark Maxwell’s equations

$\displaystyle\nabla\cdot\boldsymbol{\eta}+m^{2}\psi$ $\displaystyle=\chi
m^{2}\phi$ (6a) $\displaystyle\nabla\times\boldsymbol{\beta}+m^{2}\mathbf{X}$
$\displaystyle=\partial_{t}\boldsymbol{\eta}+\chi m^{2}\mathbf{A}$ (6b)
$\displaystyle\nabla\times\boldsymbol{\eta}$
$\displaystyle=-\partial_{t}\boldsymbol{\beta}$ (6c)
$\displaystyle\nabla\cdot\beta$ $\displaystyle=0.$ (6d)

Additionally, the Lorentz force law is only modified by a renormalisation of
the electric charge Goldhaber and Nieto (1971), so we may write

$\mathbf{F}=\tilde{e}(\mathbf{E}+\mathbf{v}\times\mathbf{B}),$ (7)

where $\mathbf{v}$ is the velocity of a charged particle. Since the Lorentz
force law (7) is only modified by its overall scale, we may carry over all of
the equations of kinetic plasma theory with the substitution
$e\rightarrow\tilde{e}$.

To obtain a closed set of equations, we must further specify the equations of
motion for the plasma. By the observation at the end of the previous
paragraph, we get Ohm’s law:

$\tilde{e}\mathbf{J}=\sigma(\mathbf{E}+\mathbf{v}\times\mathbf{B}),$ (8)

where $\sigma$ is the plasma conductivity. For plasmas in the solar system, we
may take the large $\sigma$ limit; we further take the quasistatic limit to
find the set of MHD equations that we will use:

$\displaystyle(-\nabla^{2}+\chi^{2}m^{2})\mathbf{B}$
$\displaystyle=\tilde{e}\nabla\times\mathbf{J}+\chi m^{2}\boldsymbol{\beta},$
(9a) $\displaystyle\nabla\times\mathbf{E}$
$\displaystyle=-\partial_{t}\mathbf{B},$ (9b)
$\displaystyle\mathbf{E}+\mathbf{v}\times\mathbf{B}$ $\displaystyle=0,$ (9c)
$\displaystyle(-\nabla^{2}+m^{2})\boldsymbol{\beta}$ $\displaystyle=\chi
m^{2}\mathbf{B},$ (9d) $\displaystyle\nabla\times\boldsymbol{\eta}$
$\displaystyle=-\partial_{t}\boldsymbol{\beta},$ (9e)
$\displaystyle\partial_{t}\rho+\nabla\cdot(\rho\mathbf{v})$ $\displaystyle=0,$
(9f) $\displaystyle\rho\frac{d\mathbf{v}}{t}$ $\displaystyle=-\nabla
P+\tilde{e}\mathbf{J}\times\mathbf{B}.$ (9g)

Here $P$ is the fluid pressure, $\mathbf{v}$ the fluid velocity, $\mathbf{J}$
is still the current density in units of $\tilde{e}$, but we have switched
notation and $\rho$ is now the fluid density as the charge density does not
appear any longer.

##  III Earth’s magnetic field

The sources of magnetic fields around the Earth are well enough known that
bounds have been set on the photon’s mass from measurements of the magnetic
field caused by the magnetic dipole moment of the Earth Goldhaber and Nieto
(1968); Fischbach  _et al._ (1994). These bounds are based on the technique of
Schrödinger Schrödinger (1943). As the process for bounding dark photons will
be similar, let us briefly recapitulate the logic of the massive photon case.

If the photon has a mass $\mu$, the magnetic field caused by a dipole
$\mathbf{D}=D\hat{z}$ is

$\mathbf{B}(\mathbf{x})=\frac{De^{-\mu r}}{4\pi r^{3}}\Big{[}\left(1+\mu
r+\frac{1}{3}(\mu
r)^{2}\right)(3\hat{z}\cdot\hat{\mathbf{x}}\hat{\mathbf{x}}-\hat{z})-\frac{2}{3}(\mu
r)^{2}\hat{z}\Big{]}.$ (10)

At a constant radius, the terms in the first bracket of this equation simply
effect a rescaling of the magnetic field compared to the massless case. The
last term is a new effect though, and from the point of view of massless
magnetostatics it appears like an external field antiparallel to the dipole
direction. The ratio $S_{\mu}$ of the external field caused by the mass of the
photon to the equatorial field strength is

$S_{\mu}=\frac{2}{3}\frac{(\mu R)^{2}}{1+\mu R+(\mu R)^{2}/3},$ (11)

where the measurement is taken at a radius $R$. There are a number of true
external fields, which for the Earth are known. Subtracting these off leaves a
bound on $S_{\mu}$ which translates into a bound on the photon mass
$\mu\lesssim$8\text{\times}{10}^{-16}\text{\,}\mathrm{e}\mathrm{V}$$ Fischbach
_et al._ (1994). It is worth emphasising that the true external fields must be
known to obtain a bound from this method. This is true of the Earth; it was
not true of Jupiter at the time of early analyses that have been used in the
literature, as stated in Davis  _et al._ (1975).

The case of a dark photon follows in the same way, as analysed in Bartlett and
Lögl (1988); Fischbach  _et al._ (1994) (see Appendix A for an alternate
derivation of the field to all orders in $\chi$). The resulting ratio $S$ of
the new external field to the field strength at the equator is

$S=\frac{2}{3}\frac{(\chi
m_{X}R)^{2}e^{-m_{X}R}(1+\chi^{2})}{1+\chi^{2}e^{-m_{X}R}\left(1+m_{X}R+\frac{1}{3}(m_{X}R)^{2}\right)}.$
(12)

The above ratio is bounded by $S\leq$3.9\text{\times}{10}^{-4}\text{\,}$$ at a
distance
$R=$6.38\text{\times}{10}^{6}\text{\,}\mathrm{m}$=$3.24\text{\times}{10}^{13}\text{\,}\mathrm{e}\mathrm{V}^{-1}$$
Fischbach  _et al._ (1994). The bounds that result are an order of magnitude
more constraining than any other terrestrial experiment, and extend lower in
mass than even those from COBE.

The limit that we reviewed here relies on the knowledge of the magnetic fields
that arise from sources external to the Earth’s magnetic dipole. Without this,
it is impossible to ascribe any anomalous magnetic field to the presence of a
dark photon, or equivalently a photon mass, rather than an unaccounted
Standard Model process.

Figure 1: Constraints on the dark photon parameter space. The existing
constraints shown are: Earth, from the lack of anomalous external magnetic
fields Fischbach  _et al._ (1994); Cavendish, coming from the accuracy of
Coulomb’s lawBartlett and Lögl (1988); COBE, from a lack of spectral
distortions in the CMB Caputo  _et al._ (2020); García  _et al._ (2020). The
constraints we place from the size of any unaccounted external magnetic fields
on Jupiter Connerney  _et al._ (2020) are labelled as Jupiter (this work).

##  IV Jupiter’s magnetic field

A number of missions have measured the magnetic field of Jupiter. Pioneer 10
took data in 1973 Smith  _et al._ (1974b), and a year later the Pioneer 11
magnetometer Acuna and Ness (1975) obtained the most accurate measurement of
the intrinsic, “internal” Jovian dipole Acuna and Ness (1976) until the Juno
mission Connerney  _et al._ (2018). The value obtained by the latter for the
size of the dipole moment is 4.170$\,R_{J}\,$G, where
$R_{J}=$71\,492\text{\,}\mathrm{k}\mathrm{m}$$ is the radius of Jupiter. The
currents in the Jovian magnetosphere cause “external” magnetic fields of order
$100\text{\,}\mathrm{n}\mathrm{T}$, and so at distances of order 10 $R_{J}$
the internal and external magnetic fields become comparable. Hence, in order
to see any deviation from the predicted form of a MDM-sourced magnetic field
at such radii, detailed knowledge of the external fields is required.

Voyager 1 measured not only the magnetic field around Jupiter Ness  _et al._
(1979a), but critically also the currents and particles in its magnetosphere
Bridge  _et al._ (1979a); Krimigis  _et al._ (1979a); Warwick  _et al._
(1979). Voyager 2 obtained similar results for the field strength Ness  _et
al._ (1979b) and the magnetosphere’s contents Bridge  _et al._ (1979b);
Krimigis  _et al._ (1979b); Pearce  _et al._ (1979). By fitting an empirical
model Connerney  _et al._ (1981) to the magnetic field observed by Juno, the
Jovian magnetosphere was shown to contain a large current sheet near the
equatorial plane extending out from around 8 $R_{J}$ Connerney  _et al._
(2020). Using only Pioneer 10 data taken at radii of around 10 $R_{J}$, as
done in Kloor  _et al._ (1994), it is impossible to set a reliable bound, as
one must account for the external distortion of the magnetic field produced by
the internal Jovian dynamo.

Data from the Juno magnetometer Connerney  _et al._ (2017) was analysed by
Connerney et al. Connerney  _et al._ (2020) to find the impact of the Jovian
magnetodisc on the magnetic field. The resulting upper bound on the size of
any unaccounted for external field, projected along the direction of the
dipole and at a distance
$R=5R_{J}=$1.8\text{\times}{10}^{15}\text{\,}\mathrm{e}\mathrm{V}^{-1}$$, is
$10\text{\,}\mathrm{n}\mathrm{T}$. In the absence of error bars in Connerney
et al.’s analysis, we choose to double the size of this field to obtain a more
conservative bound. The size of the intrinsic dipole at this radius is also
needed, for which we use the analysis of Connerney  _et al._ (2018) to obtain
a value of $3\text{\times}{10}^{3}\text{\,}\mathrm{n}\mathrm{T}$. We thus
bound the ratio (12) as $S\lesssim$6\text{\times}{10}^{-3}\text{\,}$$, as
compared to ${10}^{-1}\text{\,}\mathrm{u}$sed in Kloor  _et al._ (1994). The
resulting limits are shown in Fig. 1.

One may worry that an empirical fit of the external magnetic field is not
sufficient to extract a bound. To understand why, let us review the how this
analysis work. Such a fit proceeds in two stages. First, the coefficients of
the multipole moments of the internal field are measured by fitting to
magnetic field close to Jupiter’s surface, where external fields are highly
subdominant 111This analysis in fact shows that Jupiter is not well described
as a dipole source, and the effect of higher order moments is significant,
which also must be subtracted off.. Then, using this internal model, the field
at larger distances is fit to a particular magnetodisc model. In doing this
fit, the effect of a DP may erroneously be included into the magnetodisc
parameters if there is a degeneracy. If the time-variation of the magnetic
field at fixed locations were available (as is the case on Earth, for
example), then one could hope to isolate the time-independent effect of a DP
from the external sources. However, since this is not the case, the external
magnetic field ought to be determined from direct measurements of the plasma
environment surrounding Jupiter. Models of the magnetosphere’s charged
particle distribution derived from the above fit in fact show consistency with
the measured data by Voyager 1 & 2 to within a factor of two Divine and
Garrett (1983), lending support to the conservative bound we have placed. A
Juno era consistency check would further support such a conclusion.

##  V MHD around Jupiter

The magnetic fields around Jupiter were analysed to obtain a bound on the
photon mass $\mu\leq$6\text{\times}{10}^{-16}\text{\,}\mathrm{e}\mathrm{V}$$
Davis  _et al._ (1975), even without detailed knowledge of the surrounding
currents. The procedure used to get this bound is different to that of the
previous section, as emphasised in Ryutov (1997). Let us briefly review how
this was done.

For a massive photon, Ampère’s law is modified so that a dipolar magnetic
field rescales the azimuthal current by a factor $1+(\mu r)^{2}$ compared to
the case of a massless photon Ryutov (2007). Hence the force
$\mathbf{F}=\mathbf{J}\times\mathbf{B}$ exerted on the plasma also increases
by the same factor. Thus if the photon mass is too large, there can be no
other balancing forces to achieve the observed equilibrium, and so such a
scenario may be ruled out. This modified relation between currents and the
magnetic field is the source of the above limit.

How does a dark photon modify the currents in Jupiter’s magnetic field? We
follow Ryutov Ryutov (2007) and use a purely kinematic model for the standard
electromagnetic interactions. We consider a situation in which the magnetic
field on some surface of radius $r_{0}$ is known to be a dipole, so that we
may write

$\mathbf{B}=B_{0}\cos\theta\hat{r}.$ (13)

on this surface. Then, as we move away, the vanishing divergence of
$\mathbf{B}$ implies

$\mathbf{B}=B_{0}r_{0}^{2}\frac{\cos\theta}{\mathbf{r}^{2}}\hat{r}.$ (14)

We now would like to solve equation (9d) to find the dark magnetic field in
the presence of such a source. The particular solution to this equation may be
written in Cartesian co-ordinates as

$\boldsymbol{\beta}(\mathbf{x})=\frac{\chi m^{2}}{4\pi}\int
d^{3}\mathbf{\tilde{x}}\frac{e^{-m|\mathbf{\tilde{x}}|}}{|\mathbf{\tilde{x}}|}\mathbf{B}(\mathbf{x}+\mathbf{\tilde{x}}).$
(15)

We will solve the integral perturbatively in two different limits. As in the
case of a massive photon, we ignore the homogeneous solution as unlike in
regular magnetostatics there is no spatially homogeneous solution, and we are
just left with unphysical exponential solutions.

###  V.1 Large mass limit

We first proceed in the limit $|\mathbf{x}|m\gg 1$. The strategy is to now
expand $\mathbf{B}(\mathbf{x}+\mathbf{\tilde{x}})$ around the point
$\mathbf{x}$ and then evaluate the integral (15). The exponential damping
factor will tame the divergent pieces in the Taylor series of $\mathbf{B}$ as
long as we are in the specified limit. The result is

$\boldsymbol{\beta}(\mathbf{x})=\chi
B_{0}r_{0}^{2}\left(\boldsymbol{\beta}(\mathbf{x})^{(0)}+\frac{\boldsymbol{\beta}(\mathbf{x})^{(2)}}{(mr)^{2}}+\frac{\boldsymbol{\beta}(\mathbf{x})^{(4)}}{(mr)^{4}}+\ldots\right),$
(16)

where

$\begin{split}\boldsymbol{\beta}(\mathbf{x})^{(0)}&=\frac{\cos\theta}{r^{2}}\hat{r},\\\
\boldsymbol{\beta}(\mathbf{x})^{(2)}&=-2\left(\frac{\cos\theta}{r^{2}}\hat{r}+\frac{\sin\theta}{r^{2}}\hat{\theta}\right),\\\
\boldsymbol{\beta}(\mathbf{x})^{(4)}&=-8\left(\frac{\cos\theta}{r^{2}}\hat{r}+2\frac{\sin\theta}{r^{2}}\hat{\theta}\right).\end{split}$
(17)

In order to find how this dark magnetic field affects the charge current, we
plug the result (16) into (9a). Note that the zeroth order piece in (16)
cancels what looks like a photon mass term in (9a). This means that the only
effect of a dark magnetic field is higher order in $mr$ than one would expect
from a naive comparison between the Proca Lagrangian and the Lagrangian we
consider. This may be seen explicitly in the solution for $\mathbf{J}$ which
is

$\mathbf{J}=B_{0}r_{0}^{2}\frac{\sin\theta}{r^{2}}\left(1+\chi^{2}+\frac{4\chi^{2}}{(mr)^{2}}+\ldots\right).$
(18)

If the change in the current induced by the dark magnetic field were to have a
noticeable effect, we must have at least $\chi>1$, since this result only
holds if $mr\gg 1$. However, $\chi>1$ is outside of the realm of applicability
of the tree-level model of dark photon electromagnetism that we consider, and
so no consistent constraints can be obtained.

###  V.2 Small mass limit

Let us now work in the other limit where the characteristic scale of
observation $R$ is much less than the Compton wavelength of the dark photon
$mR\ll 1$. The solution to (9d) may then be approximated as

$\boldsymbol{\beta}\simeq\chi(mR)^{2}\mathbf{B},$ (19)

and so to leading order in $mR$ (9a) is just

$(-\nabla^{2}+\chi^{2}m^{2})\mathbf{B}=\tilde{e}\nabla\times\mathbf{J}+\ldots,$
(20)

which is of the same form as that arising from the Proca Lagrangian with a
photon mass $\mu=\chi m$. We may thus recycle the known results to show that
the force exerted on the plasma by the current is

$\tilde{e}\mathbf{J}\times\mathbf{B}=\frac{B_{0}^{2}r_{0}^{2}\sin\theta\cos\theta}{4\pi
r^{5}}\left(1+(\chi mR)^{2}+\ldots\right).$ (21)

This would be incompatible with observations if the extra term in brackets
were somewhat larger than one; in the limit $mR\ll 1$, this requires a $\chi$
larger than one, which is again outside of the scope of classical physics. No
bound can be placed on dark photons from current forces around Jupiter in this
limit of the dark photon mass.

We have now shown that the observation of a steady state in the neighbourhood
of Jupiter places no consistent constraint on the dark photon coupling, in
contrast to the case of a massive photon.

##  VI MHD in the far heliosphere

The strongest reliable constraint on the photon mass
$\mu\lesssim$10\text{\times}{10}^{-18}\text{\,}\mathrm{e}\mathrm{V}$$ comes
from consideration of the solar wind in the far heliosphere Ryutov (2007). Can
one derive a constraint on dark photons in a similar manner?

The starting point is the model of solar winds first put forward by Parker
Parker (1958). The components of the magnetic field when we no longer neglect
the azimuthal component of the magnetic field that comes from the Sun’s
rotation, as is necessary at large distances $r\gg r_{0}$, are given by

$\begin{split}B_{r}&=B_{0}r_{0}^{2}\frac{\cos\theta}{r^{2}},\\\
B_{\theta}&=0,\\\
B_{\phi}&=B_{0}r_{0}^{2}\frac{\omega}{v_{m}}(r-r_{0})\frac{\sin\theta}{r^{2}},\end{split}$
(22)

where $\omega$ is the angular velocity of the sun and $v_{m}$ is the outward
velocity of the solar wind. In the far heliosphere, the magnetic field is
almost purely azimuthal. We now follow the same procedure as section V to see
what effect there is on the current.

###  VI.1 Large mass limit

Carrying out the same expansion as in section V.1, we write

$\boldsymbol{\beta}(\mathbf{x})=\boldsymbol{\beta}(\mathbf{x})^{(0)}+\frac{\boldsymbol{\beta}(\mathbf{x})^{(2)}}{(mr)^{2}}+\ldots,$
(23)

and solve (9d) perturbatively to find

$\begin{split}\boldsymbol{\beta}(\mathbf{x})^{(0)}&=\chi\mathbf{B},\\\
\boldsymbol{\beta}(\mathbf{x})^{(2)}&=-2\chi\mathbf{B}.\end{split}$ (24)

Thus (9a) gives

$\tilde{e}\nabla\times\mathbf{J}=\frac{2\mathbf{B}}{r^{2}}(1+\chi^{2}+\ldots),$
(25)

which implies that $\mathbf{J}$ has $r$\- and $\theta$-components of rough
size

$J\sim\frac{B_{\phi}}{r}(1+\chi^{2})^{1/2},$ (26)

with an associated force causing equatorial compression and radial
deceleration of magnitude

$F\sim\frac{B_{\phi}^{2}}{r}(1+\chi^{2})^{1/2}.$ (27)

If this force were bigger than the force $\rho v^{2}/r$, then the velocity
would not be radial and constant, as is observed. We thus need $F<\sigma\rho
v^{2}/r$ where $\sigma$ is some fudge factor to account for our crudeness;
this is the same limit used by Ryutov (2007) to derive a limit on the photon
mass. Translating this into dark photon parameters, we find the constraint

$1+\chi^{2}\lesssim 10^{3/2};$ (28)

the bound is trivially satisfied for any perturbative coupling.

###  VI.2 Small mass limit

As in section V.2, to leading order in the small mass limit, we again have

$(-\nabla^{2}+\chi^{2}m^{2})\mathbf{B}=\tilde{e}\nabla\times\mathbf{J}+\ldots.$
(29)

This is solved by

$\tilde{e}\mathbf{J}=2\frac{B_{\phi}}{r}\left(1+(\chi mr)^{2}\right)\hat{r},$
(30)

and again we see the dark photon correction is always negligible.

##  VII Conclusions

The magnetic field produced by a MDM receives noticeable corrections when the
photon kinematically mixes with a massive DP. Using this fact and knowledge
about the true external fields, we placed constraints on dark photons whose
Compton wavelength is of the order of the size of the Jupiter. These bounds
extend further down in DP mass than any other (apart from those that make
assumptions about the nature of dark matter). Specifically, we have
$\chi\gtrsim 10^{-1}$ for
$m_{X}\sim${10}^{-15}\text{\,}\mathrm{e}\mathrm{V}$$. This made use of a fit
of the magnetic field to an empirical model of the Jovian magnetodisc. This
bound could be made more compelling by checking for good consistency between
this model and the direct measurements of the surrounding plasma with
forthcoming Juno data. A full comparison of the functional form of the
magnetic field around Jupiter with that predicted by Maxwell theory – taking
into account the effect of the surrounding plasma – would likely give stronger
bounds, but is beyond the scope of this work.

We have shown that bounds on the photon mass coming from the observation of
current flow patterns in astrophysical magnetospheres do not constrain DPs, in
contrast to massive photons. For large DP masses, where one might hope to have
some effect, the leading order terms cancel once one consistently solves the
relevant MHD equations. DPs with small masses fare no better, causing no
significant alteration of the force felt by the plasma.

## Acknowledgements

I thank John Wheater for insightful conversations and suggestions, and Subir
Sarkar and John Wheater for comments on a draft. I am grateful to Cressida for
supporting me during this work.

##  Appendix A The magnetic field from a magnetic dipole moment

We would like to calculate the magnetic field arising from a point magnetic
dipole moment, in analogy with (10). The relevant equations we need to solve
are

$\displaystyle(-\nabla^{2}+\chi^{2}m^{2})\mathbf{A}$
$\displaystyle=\tilde{e}\mathbf{J}+\chi m^{2}\mathbf{X},$ (31a)
$\displaystyle(-\nabla^{2}+m^{2})\mathbf{X}$ $\displaystyle=\chi
m^{2}\mathbf{A},$ (31b)

where we work here in Coulomb gauge $\nabla\cdot\mathbf{A}=0$ and note that in
the quasistatic limit this also implies $\nabla\cdot\mathbf{X}=0$. Upon
Fourier transforming (31), these differential equations become algebraic and
we may solve for the vector potential to find

$\mathbf{A}(x)=\tilde{e}\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}\frac{\mathbf{\tilde{J}}(\mathbf{k})e^{i\mathbf{k}\cdot\mathbf{x}}}{k^{2}+\chi^{2}m^{2}\frac{k^{2}}{k^{2}+m^{2}}},$
(32)

where $\tilde{\mathbf{J}}(\mathbf{k})$ is the Fourier transform of
$\tilde{e}\mathbf{J}(\mathbf{x})$. Approximating the magnetic dipole as coming
form a point-like source at the origin, we have
$\mathbf{J}(\mathbf{x})=-\mathbf{D}\times\nabla\delta(\mathbf{x})$, where
$\mathbf{D}=D\hat{z}$ is the magnetic dipole moment. We then have
$\tilde{\mathbf{J}}=-i\mathbf{D}\times\mathbf{k}$.

We now need to evaluate integrals of the form

$\mathbf{I}=\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}\frac{\mathbf{k}e^{i\mathbf{k}\cdot\mathbf{x}}}{k^{2}+\chi^{2}m^{2}\frac{k^{2}}{k^{2}+m^{2}}},$
(33)

which is most readily done in spherical co-ordinates. If we rotate our axes so
that the polar axis of the momentum integral aligns with $\hat{\mathbf{x}}$,
we see that the $x-$ and $y$-components have vanishing azimuthal angle $\phi$
integrals. The polar integral needed is found to be

$\displaystyle\int_{0}^{\pi}d\theta\,e^{ikr\cos\theta}\cos\theta\sin\theta$
$\displaystyle=2i\frac{-kr\cos(kr)+\sin(kr)}{(kr)^{2}}.$ (34a)

After integrating $\phi$ from 0 to $2\pi$, we may extend the remaining even
integral over the whole real axis and evaluate it by the Residue Theorem as

$\frac{i}{r^{2}}\int^{\infty}_{-\infty}\frac{dk}{(2\pi)^{2}}\,k\frac{-kr\cos(kr)+\sin(kr)}{k^{2}+\chi^{2}m^{2}\frac{k^{2}}{k^{2}+m^{2}}}=\frac{i}{4\pi
r^{2}}\left(\frac{1+\chi^{2}e^{-mr\sqrt{1+\chi^{2}}}}{1+\chi^{2}}+\frac{mr\chi^{2}e^{-mr\sqrt{1+\chi^{2}}}}{\sqrt{1+\chi^{2}}}\right).$
(35)

If we now rotate back to the original co-ordinate system, we find that
$\mathbf{I}$ points in the same direction as $\mathbf{x}$. Recalling (32),
which was $\mathbf{A}=-i\mathbf{D}\times\mathbf{I}$, we find

$\mathbf{A}(\mathbf{x})=\frac{D}{4\pi
r^{2}}\left(\frac{1+\chi^{2}e^{-mr\sqrt{1+\chi^{2}}}}{1+\chi^{2}}+\frac{mr\chi^{2}e^{-mr\sqrt{1+\chi^{2}}}}{\sqrt{1+\chi^{2}}}\right)\hat{\phi},$
(36)

for $\mathbf{D}$ in the $z$-direction. Taking the curl of this equation gives
us the magnetic field sourced by a magnetic dipole moment when a dark photon
coupling is present:

$\begin{split}\mathbf{B}(\mathbf{x})=&\frac{D}{4\pi r^{3}(1+\chi^{2})^{2}}\\\
\cdot\Big{\\{}&\Big{[}1+\chi^{2}e^{-m_{X}r}\left(1+m_{X}r+\frac{1}{3}(m_{X}r)^{2}\right)\Big{]}\Big{(}3\mathbf{\hat{x}}\cdot\hat{z}\mathbf{\hat{x}}-\hat{z}\Big{)}\\\
&-\frac{2}{3}(\chi
m_{X}r)^{2}e^{-m_{X}r}(1+\chi^{2})\hat{z}\Big{\\}}.\end{split}$ (37)

## References

  * Holdom (1986) B. Holdom, Physics Letters B 166, 196 (1986).
  * Jaeckel and Ringwald (2010) J. Jaeckel and A. Ringwald, Annual Review of Nuclear and Particle Science 60, 405 (2010), arXiv:1002.0329 .
  * Fabbrichesi  _et al._ (2020) M. Fabbrichesi, E. Gabrielli, and G. Lanfranchi,  (2020), 10.1007/978-3-030-62519-1, arXiv:2005.01515 [hep-ph] .
  * Williams  _et al._ (1971) E. R. Williams, J. E. Faller, and H. A. Hill, Phys. Rev. Lett. 26, 721 (1971).
  * Bartlett and Lögl (1988) D. F. Bartlett and S. Lögl, Physical Review Letters 61, 2285 (1988).
  * Fixsen  _et al._ (1996) D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, The Astrophysical Journal 473, 576–587 (1996).
  * Caputo  _et al._ (2020) A. Caputo, H. Liu, S. Mishra-Sharma, and J. T. Ruderman, Physical Review Letters 125 (2020), 10.1103/physrevlett.125.221303.
  * García  _et al._ (2020) A. A. García, K. Bondarenko, S. Ploeckinger, J. Pradler, and A. Sokolenko, Journal of Cosmology and Astroparticle Physics 2020, 011–011 (2020).
  * Goldhaber and Nieto (1971) A. S. Goldhaber and M. M. Nieto, Reviews of Modern Physics 43, 277 (1971).
  * Fischbach  _et al._ (1994) E. Fischbach, H. Kloor, R. A. Langel, A. T. Y. Liu, and M. Peredo, Phys. Rev. Lett. 73, 514 (1994).
  * Kloor  _et al._ (1994) H. Kloor, E. Fischbach, C. Talmadge, and G. L. Greene, Physical Review D 49, 2098 (1994).
  * Smith  _et al._ (1974a) E. J. Smith, L. Davis, D. E. Jones, D. S. Colburn, P. J. Coleman, P. Dyal, and C. P. Sonett, Science 183, 305 (1974a).
  * Davis  _et al._ (1975) L. Davis, A. S. Goldhaber, and M. M. Nieto, Physical Review Letters 35, 1402 (1975).
  * Bolton  _et al._ (2017) S. J. Bolton, J. Lunine, D. Stevenson, J. E. P. Connerney, S. Levin, T. C. Owen, F. Bagenal, D. Gautier, A. P. Ingersoll, G. S. Orton, T. Guillot, W. Hubbard, J. Bloxham, A. Coradini, S. K. Stephens, P. Mokashi, R. Thorne, and R. Thorpe, Space Science Reviews 213, 5 (2017).
  * Ryutov (2007) D. D. Ryutov, Plasma Physics and Controlled Fusion 49, B429 (2007).
  * Ness and Burlaga (2001) N. Ness and L. Burlaga, Journal of Geophysical Research 106, 15803 (2001).
  * Burlaga  _et al._ (2003) L. F. Burlaga, N. F. Ness, Y. Wang, and N. R. Sheeley, AIP Conference Proceedings 679, 39 (2003).
  * Goldhaber and Nieto (1968) A. S. Goldhaber and M. M. Nieto, Phys. Rev. Lett. 21, 567 (1968).
  * Schrödinger (1943) E. Schrödinger, Proceedings of the Royal Irish Academy. Section A: Mathematical and Physical Sciences 49, 135 (1943).
  * Connerney  _et al._ (2020) J. E. P. Connerney, S. Timmins, M. Herceg, and J. L. Joergensen, Journal of Geophysical Research: Space Physics 125, e2020JA028138 (2020).
  * Smith  _et al._ (1974b) E. Smith, L. Davis, D. E. Jones, P. Coleman, D. Colburn, P. Dyal, C. Sonett, and A. M. Frandsen, Journal of Geophysical Research 79, 3501 (1974b).
  * Acuna and Ness (1975) M. A. Acuna and N. F. Ness, Space Science Instrumentation 1, 177 (1975).
  * Acuna and Ness (1976) M. Acuna and N. Ness, Journal of Geophysical Research 81, 2917 (1976).
  * Connerney  _et al._ (2018) J. E. P. Connerney, S. Kotsiaros, R. J. Oliversen, J. R. Espley, J. L. Joergensen, P. S. Joergensen, J. M. G. Merayo, M. Herceg, J. Bloxham, K. M. Moore, S. J. Bolton, and S. M. Levin, Geophysical Research Letters 45, 2590 (2018).
  * Ness  _et al._ (1979a) N. F. Ness, M. H. Acuna, R. P. Lepping, L. F. Burlaga, K. W. Behannon, and F. M. Neubauer, Science 204, 982 (1979a).
  * Bridge  _et al._ (1979a) H. S. Bridge, J. W. Belcher, A. J. Lazarus, J. D. Sullivan, R. L. Mcnutt, F. Bagenal, J. D. Scudder, E. C. Sittler, G. L. Siscoe, V. M. Vasyliunas, C. K. Goertz, and C. M. Yeates, Science 204, 987 (1979a).
  * Krimigis  _et al._ (1979a) S. M. Krimigis, T. P. Armstrong, W. I. Axford, C. O. Bostrom, C. Y. Fan, G. Gloeckler, L. J. Lanzerotti, E. P. Keath, R. D. Zwickl, J. F. Carbary, and D. C. Hamilton, Science 204, 998 (1979a).
  * Warwick  _et al._ (1979) J. W. Warwick, J. B. Pearce, A. C. Riddle, J. K. Alexander, M. D. Desch, M. L. Kaiser, J. R. Thieman, T. D. Carr, S. Gulkis, A. Boischot, C. C. Harvey, and B. M. Pedersen, Science 204, 995 (1979).
  * Ness  _et al._ (1979b) N. F. Ness, M. H. Acuna, R. P. Lepping, L. F. Burlaga, K. W. Behannon, and F. M. Neubauer, Science 206, 966 (1979b).
  * Bridge  _et al._ (1979b) H. S. Bridge, J. W. Belcher, A. J. Lazarus, J. D. Sullivan, F. Bagenal, R. L. Mcnutt, K. W. Ogilvie, J. D. Scudder, E. C. Sittler, V. M. Vasyliunas, and C. K. Goertz, Science 206, 972 (1979b).
  * Krimigis  _et al._ (1979b) S. M. Krimigis, T. P. Armstrong, W. I. Axford, C. O. Bostrom, C. Y. Fan, G. Gloeckler, L. J. Lanzerotti, E. P. Keath, R. D. Zwickl, J. F. Carbary, and D. C. Hamilton, Science 206, 977 (1979b).
  * Pearce  _et al._ (1979) J. B. Pearce, A. C. Riddle, J. W. Warwick, J. K. Alexander, M. D. Desch, M. L. Kaiser, J. R. Thieman, T. D. Carr, S. Gulkis, A. Boischot, Y. Leblanc, B. M. Pedersen, and D. H. Staelin, Science 206, 991 (1979).
  * Connerney  _et al._ (1981) J. E. P. Connerney, M. H. Acuña, and N. F. Ness, Journal of Geophysical Research: Space Physics 86, 8370 (1981).
  * Connerney  _et al._ (2017) J. E. P. Connerney, M. Benn, J. B. Bjarno, T. Denver, J. Espley, J. L. Jorgensen, P. S. Jorgensen, P. Lawton, A. Malinnikova, J. M. Merayo, S. Murphy, J. Odom, R. Oliversen, R. Schnurr, D. Sheppard, and E. J. Smith, Space Science Reviews 213, 39 (2017).
  * Divine and Garrett (1983) N. Divine and H. B. Garrett, Journal of Geophysical Research: Space Physics 88, 6889 (1983).
  * Ryutov (1997) D. D. Ryutov, Plasma Physics and Controlled Fusion 39, A73 (1997).
  * Parker (1958) E. N. Parker, The Astrophysical Journal 128, 664 (1958).