# Gravitational wave detection of DFSZ axion models

Aidi Yang<EMAIL_ADDRESS>Fa Peng Huang<EMAIL_ADDRESS>MOE
Key Laboratory of TianQin Mission, TianQin Research Center for Gravitational
Physics & School of Physics and Astronomy, Frontiers Science Center for
TianQin, Gravitational Wave Research Center of CNSA, Sun Yat-sen University
(Zhuhai Campus), Zhuhai 519082, China

###### Abstract

As a promising dark matter candidate, the axion particle can naturally solve
the strong CP problem of the standard model through the U(1) Peccei-Quinn
symmetry-breaking process. However, detecting this high-energy process poses a
significant challenge for colliders. We precisely calculate the phase
transition gravitational wave signals of this symmetry-breaking process in the
popular DFSZ axion model. By comparing the signals with the expected
sensitivity of Cosmic Explorer and Einstein Telescope, we demonstrate that
Cosmic Explorer could detect this process with a signal-to-noise ratio higher
than 8 for the benchmark points. Furthermore, by performing Fisher analysis,
we find that if signals are observed, the bubble wall velocity will be the
first phase transition parameter to be determined. The future gravitational
wave detectors are expected to provide a new approach to exploring concrete
axion models.

##  I Introduction

The strong CP problem in the standard model and the microscopic nature of dark
matter are long-standing fundamental problems in particle physics and
cosmology. Axion particle Peccei:1977hh ; Peccei:1977ur ; Weinberg:1977ma ;
Wilczek:1977pj ; Kim:1979if ; Shifman:1979if ; Zhitnitsky:1980tq ; Dine:1981rt
; Sikivie:2020zpn  can naturally solve the strong CP problem and also be a
promising dark matter candidate simply through the misalignment mechanism
Preskill:1982cy ; Dine:1982ah ; Sikivie:2006ni . The U(1) Peccei-Quinn (PQ)
symmetry-breaking process is essential to solve the strong CP problem
Peccei:1977hh ; Peccei:1977ur . The scale of the PQ symmetry breaking process
determines the axion mass and its couplings with the standard model particles.
Thus, it is important to explore the PQ symmetry breaking scale and the
symmetry breaking dynamics by experiments. However, the colliders and other
traditional experiments are challenging to detect this high-energy process.
Motivated by the great achievements of gravitational wave (GW) observations,
we consider the possible GW detection of the dynamical process of the U(1) PQ
symmetry breaking since in the early universe the PQ symmetry breaking process
might generate a strong first-order phase transition(SFOPT) Croon:2019iuh ;
Dev:2019njv ; VonHarling:2019rgb ; DelleRose:2019pgi ; Ghoshal:2020vud . In
this case, an observable stochastic gravitational wave background (SGWB) can
be produced. We investigate the possibility of exploring this symmetry-
breaking process through GW experiments. Taking one concrete model of Dine-
Fischler-Srednicki-Zhitnitsky (DFSZ) axion theory as a typical example
Dine:1981rt ; Zhitnitsky:1980tq , we precisely calculate the phase transition
dynamics and its associated phase transition GW signals from bubble collisions
Caprini:2015zlo ; Huber:2008hg ; Kosowsky:1992vn ; Kosowsky:1992rz ;
Lewicki:2020azd , sound waves Hindmarsh:2013xza ; Hindmarsh:2015qta , and
turbulence Caprini:2009yp ; RoperPol:2019wvy . Compared with the expected
sensitivity curves of Cosmic Explorer (CE) LIGOScientific:2016wof ;
Reitze:2019iox ; Hall:2022dik ; Hall:2020dps  and the Einstein Telescope (ET)
Sathyaprakash:2012jk ; Maggiore:2019uih . We demonstrate that the GW signals
from axion models can be detected by CE with a signal-to-noise ratio (SNR)
higher than 8. Furthermore, using the method of Fisher matrix (FM)
Fisher:1922saa , we find that the bubble wall velocity will be the first phase
transition parameter to be fixed if the signals are observed.

For the sake of integrity, we briefly review the axion models in Sec. II and
the detection principles of SGWB in Sec. III. Readers familiar with these two
sections can skip them and proceed directly to Sec. IV. In Sec. IV, we
precisely calculate the phase transition GW signals generated by the PQ
symmetry breaking process in the DFSZ model and discuss their detectability at
future GW detectors, such as CE and ET. Then, we perform Fisher analysis for
the key phase transition parameters at CE in Sec. V. A summary is given in
Sec. VI. More details are shown in the Appendix.

##  II Axion and DFSZ models

###  II.1 Strong CP problem and axion solution

As a $SU(3)_{c}$ non-Abelian gauge theory in 4-dimensional spacetime, the QCD
Lagrangian density contains a CP violation term originating from the instanton
effects

$\mathcal{L}_{\mathrm{QCD}}\supset\bar{\theta}\frac{g_{s}^{2}}{32\pi^{2}}G^{a\mu\nu}\tilde{G}_{a\mu\nu}$
(1)

where $g_{s}$ is the $SU(3)_{c}$ gauge coupling constant, $G^{a\mu\nu}$ is the
gluon field-strength tensor, and $\tilde{G}_{a\mu\nu}$ is the dual of the
gluon field-strength tensor. This CP-violating $\bar{\theta}$ can induce a
neutron electric dipole moment of order Pospelov:1999ha

$d_{n}^{(th)}\simeq 2.4(1.0)\times 10^{-16}\bar{\theta}\,\,\text{e~{}cm}$ (2)

However, current experiments give an upper bound to the electric dipole moment
of the neutron as Baker:2006ts ; Abel:2020pzs

$|d_{n}^{(exp)}|<1.8\times 10^{-26}\,\text{e~{}cm}$ (3)

This requires negligible $\bar{\theta}$, namely,

$\bar{\theta}<10^{-10}$ (4)

Why is the $\bar{\theta}$ parameter in QCD so small? This is the so-called
strong CP problem in the standard model of particle physics.

Axion models provide natural solutions to dynamically explain the strong CP
problem. The strong CP problem can be resolved by introducing a new global
U(1) PQ symmetry. When the U(1) PQ symmetry is spontaneously broken, it gives
rise to a new Nambu-Goldstone boson called axion. The mass and interaction
properties of the axion are closely related to the introduced PQ scalar field
$\sigma$. The interaction between the axion field $a$ and the gluon field
takes the form of $aG\tilde{G}/f_{a}$, where $f_{a}$ is the decay constant of
the axion. The presence of the axion allows for the replacement of the CP-
violating term $\bar{\theta}$ with a dynamic field known as the axion field.
The axion field dynamically adjusts itself to a value that effectively cancels
out the $\bar{\theta}$ by sitting at the minimum of the axion potential.
Within this dynamical mechanism, the strong CP problem is solved by the
emergence of the axion particle.

The original Peccei-Quinn-Weinberg-Wilczek (PQWW) axion model Peccei:1977hh ;
Peccei:1977ur ; Weinberg:1977ma ; Wilczek:1977pj  assumed the PQ symmetry
breaking happens just at the electroweak energy scale, namely,

$f_{a}\simeq 246~{}\rm GeV.$ (5)

It was immediately excluded by experimental data, including beam-dump and LEP
experiments results Asano:1981nh ; Mimasu:2014nea . Later, two types of
invisible axion models are proposed with much higher PQ energy scale to
decrease the axion mass and reduce its couplings to ordinary matter,

$f_{a}\gg 246~{}\rm GeV$ (6)

###  II.2 Invisible axion model

####  II.2.1 Kim-Shifman-Vainshtein-Zakharov model

One type of invisible axion is the Kim-Shifman-Vainshtein-Zakharovb(KSVZ)
model where exotic quark $Q_{E}+Q_{E}^{C}$ with the and a new complex gauge-
singlet scalar $\sigma$ are included Kim:1979if ; Shifman:1979if ; Kim:1986ax
. The relevant Lagrangian density is

$\mathcal{L}_{\mathrm{KSVZ}}\supset\lambda_{K}\sigma
Q_{E}Q_{E}^{C}-\lambda\left(|\sigma|^{2}-\frac{f_{a}^{2}}{2}\right)^{2}$ (7)

All standard model fields are neutral under the $U(1)$ PQ symmetry in its
linear realization. The PQ charge of the scalar field $X_{\sigma}=1$ and the
PQ charge of the exotic quark $X_{Q_{E}}=X_{Q_{E}^{C}}=-1/2$. The Yukawa
coupling constant $\lambda_{K}$ plays a crucial role in this model and can be
treated as a free parameter. This term facilitates the interaction between the
PQ field and the heavy exotic quarks. Like the PQWW model, a global $U(1)$ PQ
symmetry generates a shift in the axion field through a chiral rotation.
However, the potential in this model spontaneously breaks the global $U(1)$ PQ
symmetry. Therefore, the KSVZ model offers a solution to the strong CP
problem.

####  II.2.2 Dine-Fischler-Srednicki-Zhitnitsky model

In this work, we focus on the other type of invisible axion model, namely, the
well-known DFSZ model. All the standard model fields are PQ-charged as the
original PQWW model. This model consists of two Higgs doublets and a complex
scalar field, denoted as $\sigma$. As one component of the complex scalar
field, the DFSZ axion is linked to the standard model through the two Higgs
doublets $H_{u}$ and $H_{d}$. The scalar field transforms according to the
following expression

$\sigma\rightarrow e^{iX_{\sigma}\gamma}\sigma,\quad H_{\mathrm{u}}\rightarrow
e^{iX_{u}\gamma}H_{\mathrm{u}},\quad H_{\mathrm{d}}\rightarrow
e^{iX_{d}\gamma}H_{\mathrm{d}}$ (8)

where

$X_{\mathrm{u}}+X_{\mathrm{d}}=2X_{\sigma}=1$ (9)

The general tree-level potential is Dine:1981rt

$\displaystyle
V\left(\sigma,H_{\mathrm{u}},H_{\mathrm{d}}\right)=\lambda_{\mathrm{u}}\left(\left|H_{\mathrm{u}}\right|^{2}-v_{\mathrm{u}}^{2}\right)^{2}+\lambda_{\mathrm{d}}\left(\left|H_{\mathrm{d}}\right|^{2}-v_{\mathrm{d}}^{2}\right)^{2}$
(10)
$\displaystyle\quad+\lambda\left(|\sigma|^{2}-v_{\sigma}^{2}\right)^{2}+\left(\lambda_{a}\left|H_{\mathrm{u}}\right|^{2}+\lambda_{b}\left|H_{\mathrm{d}}\right|^{2}\right)|\sigma|^{2}$
$\displaystyle\quad+\lambda_{c}\left(H_{\mathrm{u}}^{i}\epsilon_{ij}H_{\mathrm{d}}^{j}\sigma^{n}+\text{
h.c.
}\right)+\lambda_{d}\left|H_{\mathrm{u}}^{i}\epsilon_{ij}H_{\mathrm{d}}^{j}\right|^{2}+\lambda_{e}\left|H_{\mathrm{u}}^{*}H_{\mathrm{d}}\right|^{2}$

where $n=1/X_{\sigma}$ ($X_{\sigma}$ is PQ charge of $\sigma$), $v_{\sigma}$
is the vev of the $\sigma$. Additionally, all couplings in this model are
assumed to be real. The proposed model suggests that $n$ is equal to 2, and
the number parameter of the domain wall is denoted by $N_{DW}$, which is 6. We
set $X_{u}=1$, $X_{d}=0$, and $X_{\sigma}=\frac{1}{2}$. As before, the
potential spontaneously breaks the global U(1) PQ symmetry. Our research
focused on calculating the first DFSZ-type models, specifically the minimal
DFSZ model known as DFSZ-I. In DFSZ-I, only an additional complex scalar field
with $n=2$ is introduced.

###  II.3 axion dark matter and axion detection

Axion particles can be a natural dark matter candidate produced by the
misalignment mechanism. There are various experiments to detect the axion
particles or axion dark matter, such as helioscope searches Sikivie:1983ip ;
Lazarus:1992ry ; Moriyama:1998kd , haloscopes ADMX:2009iij , and the radio
telescope searches Pshirkov:2007st ; Huang:2018lxq ; Hook:2018iia , see recent
review Sikivie:2020zpn  for details. GW might provide new approaches to
explore axion Chiang:2020aui ; Croon:2019iuh ; DelleRose:2019pgi ;
Ahmadvand:2021vxs ; Dev:2019njv ; Ghoshal:2020vud ; VonHarling:2019rgb ;
Yang:2023aak ; Yang:2023vwm ; Xie:2022uvp . To explore the possibility that
the DFSZ axion model can generate an observable SGWB, we continue to review
the detection principles of the SGWB and the sensitivity curves of GW
detectors in the next section. This will lay the foundation for us to
calculate the strength of the phase transition GW signals predicted by the
DFSZ model and evaluate their detectability in Sec.IV.

##  III  Stochastic gravitational wave background detection

The theoretical framework of the DFSZ axion model in the last section predicts
that it might generate a SFOPT, which in turn can produce an observable SGWB.
Unlike GW signals from individual astrophysical events, SGWB michelson1987 ;
christensen1992 ; Allen:1997ad  arises from the superposition of GWs emitted
by countless sources across the universe, such as binary mergers, supernova
explosions, cosmic strings, and other astrophysical processes. As the SGWB
originates from all directions and consists of GWs with different frequencies,
its signal is typically weak and distributed over a wide frequency range. We
cannot directly measure the properties of individual SGWB sources. To detect
the SGWB, multiple GW detectors are typically employed to observe
simultaneously, and their data is subjected to long-term integration and
cross-correlation analysis. The correlation analysis of the GW detector
network can extract the SGWB signal. In this section, we briefly summarize the
key points for SGWB detection Maggiore:2007ulw .

###  III.1 Signal analysis

In the actual detection process, the output signal of the detector $s_{gw}(t)$
consists of the GW signal $h_{gw}(t)$ and noise $n_{gw}(t)$

$s_{gw}(t)=h_{gw}(t)+n_{gw}(t)$ (11)

Here, it is assumed that the noise $n_{gw}(t)$ is a stationary random process,
and its autocorrelation function depends only on the time difference, not on
the specific time point. The autocorrelation function of the noise $R(\tau)$
is defined as

$R(\tau)=\langle
n_{gw}(t)n_{gw}(t+\tau)\rangle=\int_{-\infty}^{\infty}n_{gw}(t)n_{gw}(t+\tau)dt$
(12)

The noise power spectral density (PSD) function $S_{n}(f)$ describes the
distribution of noise power in the frequency domain. According to the Wiener-
Khinchin theorem, there exists a Fourier transform relationship between the
autocorrelation function and the PSD for a stationary random process. For a
stationary random process,
$\langle\tilde{n}_{gw}^{*}(f)\tilde{n}_{gw}(f^{\prime})\rangle$ is non-zero
only when $f=f^{\prime}$ ($\tilde{n}_{gw}(f)$ is Fourier transform of
$n_{gw}(t)$). Combining these concepts, we obtain the following equation

$\left\langle\tilde{n}_{gw}^{*}(f)\tilde{n}_{gw}\left(f^{\prime}\right)\right\rangle=\delta\left(f-f^{\prime}\right)\frac{1}{2}S_{n}(f)$
(13)

Eq. (13) describes the statistical properties of stationary random noise in
the frequency domain. Since $n_{gw}(t)$ is real,
$\tilde{n}_{gw}(-f)=\tilde{n}_{gw}^{*}(f)$ and $S_{n}(-f)=S_{n}(f)$. Limiting
the time interval to $-T_{t}/2<t<T_{t}/2$, $T_{t}$ is the observation time. We
have

$\delta(f=0)\rightarrow\left[\int_{-T_{t}/2}^{T_{t}/2}dte^{i2\pi
ft}\right]_{f=0}=T_{t}$ (14)

Therefore, it can be derived that $n_{gw}^{2}(t)=\int_{0}^{\infty}dfS_{n}(f)$.
This formula shows how $S_{n}(f)$ is associated with the mean square value of
the noise in the time domain. Other methods are usually used to analyze the
noise signal, such as using the PSD or autocorrelation function of the noise.
The general form of the detector output is given by Eq. (11). To extract the
GW signal from the noise, the cross-correlation between the GW signal and the
detector signal is calculated first

$\frac{1}{T_{t}}\int_{0}^{T_{t}}dts_{gw}(t)h_{gw}(t)=\frac{1}{T_{t}}\int_{0}^{T_{t}}h_{gw}^{2}(t)+\frac{1}{T_{t}}\int_{0}^{T_{t}}dtn_{gw}(t)h_{gw}(t)$
(15)

then

$\frac{1}{T_{t}}\int_{0}^{T_{t}}dth_{gw}^{2}(t)\sim h_{0}^{2}$ (16)

where $h_{0}$ is the characteristic amplitude of $h_{gw}(t)$. Since the noise
$n_{gw}(t)$ and the signal $h_{gw}(t)$ are uncorrelated Maggiore:2007ulw

$\frac{1}{T_{t}}\int_{0}^{T_{t}}dtn_{gw}(t)h_{gw}(t)\sim\left(\frac{T_{0}}{T_{t}}\right)^{1/2}n_{0}h_{0}$
(17)

where $n_{0}$ is the characteristic amplitude of $n_{gw}(t)$. $T_{0}$ is the
oscillation period. Assuming we know the form of $h_{gw}(t)$, then we consider
which filter function can maximize the SNR. This technique is called matched
filtering.

$\hat{s}=\int_{-\infty}^{\infty}dts_{gw}(t)K(t)$ (18)

where $K(t)$ is called the filter function. SNR is defined as $S/N$ with $S$
being the expected value of $\hat{s}$ when the signal is present, and $N$ is
the root mean square value of $\hat{s}$ when the signal is absent. Because
$\langle n_{gw}(t)\rangle=0$, we have

$\displaystyle S$ $\displaystyle=\int_{-\infty}^{\infty}dt\langle
s_{gw}(t)\rangle K(t)$ (19)
$\displaystyle=\int_{-\infty}^{\infty}dth_{gw}(t)K(t)$
$\displaystyle=\int_{-\infty}^{\infty}df\tilde{h}_{gw}(f)\tilde{K}^{*}(f)$

In the absence of a signal, the variance of $N^{2}$ can be obtained by
calculating the convolution of $K(t)$ with $R(\tau)$ which is given by Eq.
(12)

$\displaystyle N^{2}$
$\displaystyle=\int_{-\infty}^{\infty}dtdt^{\prime}K(t)K\left(t^{\prime}\right)\int_{-\infty}^{\infty}df\frac{1}{2}S_{n}(f)e^{2\pi
if\left(t-t^{\prime}\right)}$ (20)
$\displaystyle=\int_{-\infty}^{\infty}df\frac{1}{2}S_{n}(f)|\tilde{K}(f)|^{2}$

Therefore, the SNR for a single detector is

$SNR_{\mathrm{single}}=\frac{\int_{-\infty}^{\infty}df\tilde{h}_{gw}(f)\tilde{K}^{*}(f)}{\left[\int_{-\infty}^{\infty}df\frac{1}{2}S_{n}(f)|\tilde{K}(f)|^{2}\right]^{1/2}}$
(21)

A single detector cannot effectively distinguish between SGWB signals and
detector noise because they both have random properties. However, by using two
detectors and comparing their outputs, the SGWB signal can be extracted
Allen:1997ad . This is achieved by calculating the cross-correlation function
between the outputs of the two detectors. To optimize the SNR, a filter
function $Q(t)$ is introduced Maggiore:2007ulw . The differences in the
locations and orientations of the two detectors affect the correlation of the
GW signals they detect, which can be described by the overlap reduction
function $\Gamma(f)$ Allen:1997ad . To maximize the SNR ratio, the filter
function $\tilde{Q}(f)$ which is the Fourier transform of ${Q}(t)$, needs to
be optimized. The optimized expression for the SNR is

$SNR_{\mathrm{optimized}}=T_{t}^{1/2}\frac{\int_{-\infty}^{\infty}dfS_{h}(f)\Gamma(f)\tilde{Q}(f)}{\left[\int_{-\infty}^{\infty}df|\tilde{Q}(f)|^{2}S_{2n}^{2}(f)\right]^{1/2}}$
(22)

where $S_{h}(f)$ is the PSD of the SGWB, and
$S_{2n}(f)=\left[S_{n,1}(f)S_{n,2}(f)\right]^{1/2}$ is the geometric mean of
the noise spectral densities of the two detectors. The goal is to find the
optimal form of $\tilde{Q}(f)$ that maximizes the SNR expression. This
optimization problem can be solved using the Schwarz inequality, which states
that for any two functions $u(f)$ and $v(f)$

$\left|\int_{-\infty}^{\infty}dfu(f)v^{*}(f)\right|^{2}\leq\int_{-\infty}^{\infty}df|u(f)|^{2}\int_{-\infty}^{\infty}df|v(f)|^{2}$
(23)

with equality holding if and only if $u(f)=\lambda_{s}v(f)$ for some constant
$\lambda_{s}$. We can identify

$u(f)=S_{h}(f)\Gamma(f)\tilde{Q}(f),\quad v(f)=\frac{1}{S_{2n}(f)}$ (24)

Applying the Schwarz inequality, we obtain

$SNR_{\mathrm{optimized}}\leq
T_{t}^{1/2}\left[\int_{-\infty}^{\infty}df\frac{S_{h}^{2}(f)\Gamma^{2}(f)}{S_{2n}^{2}(f)}\right]^{1/2}$
(25)

with equality holding if and only if

$\tilde{Q}(f)=\lambda_{s}\frac{S_{h}(f)\Gamma(f)}{S_{2n}^{2}(f)}$ (26)

This result gives the optimal form of the filter function $\tilde{Q}(f)$ that
maximizes the SNR. To quantify the detectability of the GW signal in a given
model with SFOPT, the SNR needs to be calculated

$SNR_{\mathrm{gw}}=\sqrt{T_{t}\int_{f_{\min}}^{f_{\max}}df\left(\frac{h_{100}^{2}\Omega_{\mathrm{GW}}}{h_{100}^{2}\Omega_{\mathrm{sens}}}\right)^{2}}$
(27)

where $h_{100}=0.72$ eliminates sensitivity to the observed value of the
Hubble constant $H_{0}$

$H_{0}=h_{100}\times 100\mathrm{~{}km}\mathrm{~{}s}^{-1}\mathrm{Mpc}^{-1}$
(28)

$h_{100}^{2}\Omega_{\text{sens }}$ corresponds to the expected sensitivity for
given detectors. To demonstrate the detectability of GW signals, we calculate
their SNR and compare it with the threshold SNR value of the detector. If the
calculated SNR is larger than the threshold, the detector can detect the GW
signal. To better understand the sensitivity of a detector, we need to
consider the sensitivity curve.

###  III.2 Sensitivity curves

The sensitivity curve of a detector characterizes its ability to detect
signals. The three commonly used definitions for sensitivity curves are the
characteristic strain, the amplitude spectral density, and the energy density
Thrane:2013oya .

####  III.2.1  The characteristic strain

The characteristic strain denoted as $h_{c}(f)$ can be derived from the metric
perturbations (GW) $h_{ab}(t,\vec{x})$

$h_{ab}(t,\vec{x})=\sum_{A}\int_{-\infty}^{\infty}df\int_{S^{2}}d\hat{\Omega}h_{A}(f,\hat{\Omega})e^{i2\pi
f(t-\hat{\Omega}\cdot\vec{x})}e_{ab}^{A}(\hat{\Omega})$ (29)

where $h_{A}(f,\hat{k})$ is the Fourier transformation of $h_{ab}(t,\vec{x})$.
$e_{ab}^{A}(\hat{k})$ represents the polarization tensors for GWs, with $A$
taking values of $+$ and $\times$. Additionally, we introduce $\hat{\Omega}$
as a unit vector denoting a specific direction on the two-sphere and
$\hat{k}=2\pi f\hat{\Omega}$. The characteristic strain $h_{c}(f)$ for GWs is
defined as

$\left\langle
h_{A}(f,\hat{k})h_{A^{\prime}}^{*}\left(f^{\prime},\hat{k}^{\prime}\right)\right\rangle=\frac{1}{16\pi}\delta\left(f-f^{\prime}\right)\delta_{AA^{\prime}}\delta^{2}\left(\hat{k},\hat{k}^{\prime}\right)h_{c}(f)^{2}f^{-1}$
(30)

where $\delta_{AA^{\prime}}$ indicates that GWs with different polarization
states are uncorrelated. $\delta^{2}(\hat{k},\hat{k}^{\prime})$ is the delta
function on the sphere which indicates that GWs propagating in different
directions are uncorrelated. The characteristic strain noise amplitude for a
single detector is denoted by $h_{n}(f)$.

####  III.2.2 The strain power spectral density

The strain PSD of a GW, $S_{h}(f)$, is known to have a relationship with
$h_{c}(f)$ Thrane:2013oya

$\sqrt{S_{h}(f)}=h_{c}(f)f^{-1/2}$ (31)

By comparison, an equivalent expression for the strain noise PSD $S_{n}(f)$
can be established

$\sqrt{S_{n}(f)}=h_{\mathrm{n}}(f)f^{-1/2}$ (32)

####  III.2.3 The energy density $\Omega_{\mathrm{gw}}(f)$

In cosmology, one usually use the variable $\Omega_{\mathrm{gw}}(f)$ which
signifies the fractional contribution of the GW energy density Allen:1997ad

$\Omega_{\mathrm{gw}}(f)=\frac{1}{\rho_{c}}\frac{d\rho_{\mathrm{gw}}}{d\ln f}$
(33)

where $\rho_{c}$ is the critical energy density. We have the following
conversion formula

$\Omega_{\mathrm{gw}}(f)=\frac{2\pi^{2}}{3H_{0}^{2}}f^{2}h_{c}^{2}(f)=\frac{2\pi^{2}S_{h}(f)}{3H_{0}^{2}}f^{3}$
(34)

The quantity frequently depicted on sensitivity curves for energy density is
the $\Omega_{\mathrm{gw}}(f)h_{100}^{2}$.

##  IV Phase transition gravitational wave of the DFSZ-I axion model and its
SNR at Cosmic Explorer

In the previous section, we reviewed the detection principles of the SGWB and
introduced the sensitivity curves of GW detectors. This lays the foundation
for our discussion on the SFOPT signals generated by the DFSZ-I axion model in
this section. We will calculate in detail the SFOPT signals produced in the
DFSZ-I model and evaluate their SNR and detectability at the future GW
detectors CE and ET. In this section, we use the finite-temperature effective
field theory Coleman:1973jx ; Dolan:1973qd  to investigate the phase
transition dynamics in DFSZ-I. At zero temperature, the tree-level potential
is given by Dine:1981rt

$\displaystyle\begin{split}V_{\text{tree}}&=-\mu_{1}^{2}|H_{u}|^{2}-\mu_{2}^{2}|H_{d}|^{2}+\lambda_{1}|H_{u}|^{4}+\lambda_{2}|H_{d}|^{4}+\lambda_{4}|H_{u}^{\dagger}H_{d}|^{2}\\\
&\quad-\mu_{3}^{2}|\sigma|^{2}+\lambda_{3}|\sigma|^{4}+\lambda_{12}|H_{u}|^{2}|H_{d}|^{2}+\lambda_{13}|\sigma|^{2}|H_{u}|^{2}\\\
&\quad+\lambda_{23}|\sigma|^{2}|H_{d}|^{2}+\left(\lambda_{5}\sigma^{2}\tilde{H}_{u}^{\dagger}H_{d}+\mathrm{h.c.}\right)\;\end{split}$
(35)

where the doublet $\tilde{H}_{u}=i\tau_{2}H_{u}^{*}$ and $\sigma$ is a complex
singlet. At the scale of PQ symmetry breaking, only the $\sigma$ field
acquires a vacuum expectation value (vev). The vev of the $\sigma$ field is
denoted as $\langle\sigma\rangle=\frac{v_{\sigma}}{\sqrt{2}}$. The singlet and
doublets can be expanded as follows

$\sigma=\frac{1}{\sqrt{2}}\left(v_{\sigma}+\sigma^{0}+i\eta_{\sigma}^{0}\right)$
(36) $H_{u}=\left(\begin{array}[]{c}h_{u}^{R}+ih_{u}^{I}\\\
\frac{1}{\sqrt{2}}(h_{u}^{0}+i\eta_{u}^{0})\end{array}\right)$ (37)
$H_{d}=\left(\begin{array}[]{c}\frac{1}{\sqrt{2}}(h_{d}^{0}+i\eta_{d}^{0})\\\
h_{d}^{R}+ih_{d}^{I}\end{array}\right)$ (38)

The free energy density of the scalar field is the effective potential.
Generally, the effective potential is composed of three parts

$V_{\text{eff }}(v_{\sigma},T)\equiv V_{\text{tree
}}(v_{\sigma})+V_{\text{CW}}(v_{\sigma})+V_{\text{T}}(v_{\sigma},T)$ (39)

where $T$ is temperature and $V_{\text{tree }}$ refers to the tree-level
potential. $V_{\text{CW}}(v_{\sigma})$ denotes the one-loop corrections at the
zero temperature. $V_{\text{T}}(v_{\sigma},T)$ represents the thermal
correction including the daisy resummation at finite temperature. The tree-
level potential can be simplified as

$V_{\text{tree
}}(v_{\sigma})=-\frac{1}{2}\mu_{3}^{2}v_{\sigma}^{2}+\frac{1}{4}\lambda_{3}v_{\sigma}^{4}$
(40)

In the $\overline{MS}$ scheme Coleman:1973jx , the zero-temperature Coleman-
Weinberg (CW) potential can be expressed as follows

$V_{\mathrm{CW}}(v_{\sigma})=\frac{1}{64\pi^{2}}\sum_{S}n_{S}m_{S}^{4}(v_{\sigma})\left[\log\frac{m_{S}^{2}(v_{\sigma})}{\mu^{2}}-\frac{3}{2}\right]$
(41)

where $m_{S}(v_{\sigma})$ represents the field-dependent mass and $n_{S}$ is
the degree of freedom ($d.o.f.$). $\mu$ represents the renormalization scale,
and the constant $\frac{3}{2}$ is for scalars. The thermal correction can be
expressed as follows Carrington:1991hz

$V_{\text{T}}(v_{\sigma},T)=\sum_{S}\frac{n_{S}T^{4}}{2\pi^{2}}J_{B}\left[\frac{m_{S}^{2}(v_{\sigma})}{T^{2}}\right]$
(42)

where $J_{B}$ is the thermal function for bosons

$J_{B}=\int_{0}^{\infty}\mathrm{d}xx^{2}\log\left[1-e^{-\sqrt{x^{2}+m_{S}^{2}/T^{2}}}\right]$
(43)

In finite-temperature field theory, when the temperature is high, the infrared
amplitude of the fields becomes very large, leading to the breakdown of
perturbative expansion. This problem is known as infrared divergence. To solve
this problem, we need to perform daisy resummation on the free energy density.
Daisy resummation resolves this issue by resumming the self-energy diagrams of
the fields. There are usually two schemes for daisy resummation: the Arnold-
Espinosa scheme Arnold:1992rz  and the Parwani scheme Parwani:1991gq . In this
work, we choose the Parwani scheme. $\Pi_{S}^{ij}$ is the contribution of
daisy resummation in thermal corrections. It can be expressed as follows
Basler:2018cwe

$\Pi_{S}^{ij}=\frac{T^{2}}{12}\sum_{k=1}^{n_{\text{scalars }}}L^{ijkk}$ (44)

where $n_{\text{scalars }}$ represents the number of scalar fields.
Furthermore, the $L^{ijkk}$ tensor can be defined as follows

$\mathcal{L}_{S}=-L^{i}\Phi_{i}-\frac{1}{2!}L^{ij}\Phi_{i}\Phi_{j}-\frac{1}{3!}L^{ijk}\Phi_{i}\Phi_{j}\Phi_{k}-\frac{1}{4!}L^{ijkl}\Phi_{i}\Phi_{j}\Phi_{k}\Phi_{l}$
(45)

The Lagrangian $\mathcal{L}_{S}$ represents the model containing scalar fields
$\Phi_{i}$, with $i={h_{u}^{R}},{h_{u}^{0}},{\eta_{u}^{0}},{h_{d}^{I}}$,
${h_{d}^{0}},{\eta_{d}^{0}},{\sigma^{0}},{\eta_{\sigma}^{0}},{h_{u}^{I}},{h_{d}^{R}}$.
The term $L^{ijkl}$ denotes the coefficient corresponding to the interaction
term $\Phi_{i}\Phi_{j}\Phi_{k}\Phi_{l}$. The thermal corrections for the
scalar fields are provided below

$\displaystyle\begin{split}&\Pi_{S}^{00}:\frac{20~{}T^{2}}{12}\lambda_{1}+\frac{12~{}T^{2}}{12}\lambda_{12}+\frac{4~{}T^{2}}{12}\lambda_{13}+\frac{4~{}T^{2}}{12}\lambda_{4}\\\
&\Pi_{S}^{11}:\frac{11~{}T^{2}}{12}\lambda_{1}+\frac{6~{}T^{2}}{12}\lambda_{12}+\frac{2~{}T^{2}}{12}\lambda_{13}+\frac{4~{}T^{2}}{12}\lambda_{4}\\\
&\Pi_{S}^{22}:\frac{11~{}T^{2}}{12}\lambda_{1}+\frac{6~{}T^{2}}{12}\lambda_{12}+\frac{2~{}T^{2}}{12}\lambda_{13}+\frac{4~{}T^{2}}{12}\lambda_{4}\\\
&\Pi_{S}^{33}:\frac{20~{}T^{2}}{12}\lambda_{2}+\frac{12~{}T^{2}}{12}\lambda_{12}+\frac{4~{}T^{2}}{12}\lambda_{23}+\frac{4~{}T^{2}}{12}\lambda_{4}\\\
&\Pi_{S}^{44}:\frac{11~{}T^{2}}{12}\lambda_{2}+\frac{6~{}T^{2}}{12}\lambda_{12}+\frac{2~{}T^{2}}{12}\lambda_{23}+\frac{4~{}T^{2}}{12}\lambda_{4}\\\
&\Pi_{S}^{55}:\frac{11~{}T^{2}}{12}\lambda_{2}+\frac{6~{}T^{2}}{12}\lambda_{12}+\frac{2~{}T^{2}}{12}\lambda_{23}+\frac{4~{}T^{2}}{12}\lambda_{4}\\\
&\Pi_{S}^{66}:\frac{3~{}T^{2}}{12}\lambda_{3}+\frac{6~{}T^{2}}{12}\lambda_{23}+\frac{6~{}T^{2}}{12}\lambda_{13}\\\
&\Pi_{S}^{77}:\frac{3~{}T^{2}}{12}\lambda_{3}+\frac{6~{}T^{2}}{12}\lambda_{23}+\frac{6~{}T^{2}}{12}\lambda_{13}\\\
&\Pi_{S}^{88}:\frac{20~{}T^{2}}{12}\lambda_{1}+\frac{12~{}T^{2}}{12}\lambda_{12}+\frac{4~{}T^{2}}{12}\lambda_{13}+\frac{4~{}T^{2}}{12}\lambda_{4}\\\
&\Pi_{S}^{99}:\frac{20~{}T^{2}}{12}\lambda_{2}+\frac{12~{}T^{2}}{12}\lambda_{12}+\frac{4~{}T^{2}}{12}\lambda_{23}+\frac{4~{}T^{2}}{12}\lambda_{4}\\\
\end{split}$ (46)

We calculate the field-dependent mass at finite temperature and subsequently
perform the following substitution in the Parwani method Parwani:1991gq

$m_{S}^{2}\rightarrow\bar{m}_{S}^{2}=m_{S}^{2}+\Pi_{S}^{ij}$ (47)

After obtaining the effective potential, we can use it to calculate some
important phase transition parameters, such as the nucleation temperature, the
characteristic length scale, the time scale of phase transition and phase
transition strength. The characteristic length is the scale of the energy
released during the phase transition, which is regarded as the mean bubble
separation $R$. Typically, we utilize the time scale of phase transition
$\beta$ to estimate $R$, given by

$R=\frac{(8\pi)^{1/3}}{\beta}v_{w}$ (48)

where $\beta$ is an approximation of the inverse time scale that plays a
critical role in determining the shape of bubbles, and $v_{w}$ is the bubble
wall velocity. The bubble wall velocity is a crucial parameter in determining
the spectra of phase transition GW which influences the strength and energy
budget Espinosa:2010hh ; Wang:2023jto ; Wang:2020zlf ; Jiang:2022btc ;
Megevand:2009gh . $\beta$ can be expressed as

$\beta=-\left.\frac{d}{dt}\left(\frac{S_{3}(T)}{T}\right)\right|_{T=T_{n}}=\left.H(T)T\frac{d}{dT}\left(\frac{S_{3}}{T}\right)\right|_{T=T_{n}}$
(49)

And the normalized parameter $\beta_{*}=\frac{\beta}{H(T)}$. $T_{n}$ is the
nucleation temperature. $S_{3}$ represents the bounce action for a three-
dimensional Euclidean space Linde:1980tt ; Linde1983 ,

$S_{3}=4\pi\int_{0}^{\infty}drr^{2}\left[\frac{1}{2}\left(\frac{dv_{\sigma}}{dr}\right)^{2}+V_{\mathrm{eff}}(v_{\sigma},T)\right]$
(50)

The nucleation rate per unit time and volume, represented by $\Gamma$
Linde:1980tt , can be approximated as $\Gamma\approx\Gamma_{0}e^{-S_{3}}$,
where $\Gamma_{0}$ is proportional to $T^{4}$. To calculate the bubble
nucleation rate, we need to obtain the bubble profiles of $v_{\sigma}$ by
solving the bounce equation with the boundary conditions

$\frac{d^{2}v_{\sigma}}{dr^{2}}+\frac{2}{r}\frac{dv_{\sigma}}{dr}=\frac{\partial
V_{\mathrm{eff}}}{\partial
v_{\sigma}},\quad\left.\frac{dv_{\sigma}}{dr}\right|_{r=0}=0,\quad
v_{\sigma}\left|{}_{r=\infty}=v_{\sigma\text{false }}\right.$ (51)

where $v_{\sigma\text{false }}$ is the scalar field value $v_{\sigma}$ in the
false vacuum. The overshooting/undershooting approach Coleman:1977py ;
Callan:1977pt  is commonly employed to solve the equation for the single-field
bounce. $T_{n}$ is the point at which the nucleation rate first surpasses the
Hubble rate. More precisely, it is defined as follows

$\frac{\Gamma\left(T_{n}\right)}{H\left(T_{n}\right)^{4}}=1$ (52)

The Hubble rate $H(T)$ during the radiation-dominated era is

$H(T)=\left[\left(\frac{\pi^{2}}{90}\right)g_{*}\right]^{1/2}\frac{T^{2}}{M_{p}}$
(53)

where $M_{p}$ is reduced Planck mass and $g_{*}$ represents the effective
number of relativistic $d.o.f.$. Within the bag model framework, the phase
transition strength parameter $\alpha$ is typically defined as Wang:2020jrd

$\alpha=\left.\frac{\Delta V_{\mathrm{eff}}-T\frac{\partial\Delta
V_{\mathrm{eff}}}{\partial T}}{\rho_{R}}\right|_{T=T_{n}}$ (54)

To derive $\alpha$, we can utilize the expressions for the effective potential
$V_{\mathrm{eff}}$ and its derivative concerning $T_{n}$, respectively. The
effective energy density $\rho_{R}$ can be expressed as
$\rho_{R}=\frac{\pi^{2}g_{*}T_{n}^{4}}{30}$. $\Delta V_{\mathrm{eff}}$ refers
to the difference in effective potential between the symmetric phase and the
broken phase at $T_{n}$.

The GW spectrum, which is generated by three primary processes: bubble
collisions, sound waves in the plasma, and turbulence, can be accurately
calculated by using these key phase transition parameters derived from the
effective potential.

  * •

Bubbles collisions
According to the envelope approximation, the shape of the uncollided bubble
walls is the key factor determining the generation of GWs. Based on the
numerical simulation, we obtain Caprini:2015zlo ; Huber:2008hg ;
Kosowsky:1992vn ; Kosowsky:1992rz ; Lewicki:2020azd

$h_{100}^{2}\Omega_{\mathrm{co}}(f)\simeq 1.67\times
10^{-5}\left(\frac{H_{*}R}{(8\pi)^{1/3}}\right)^{2}\left(\frac{\kappa_{\sigma}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{1/3}\frac{0.11v_{w}}{0.42+v_{w}^{2}}\frac{3.8\left(f/f_{\mathrm{co}}\right)^{2.8}}{1+2.8\left(f/f_{\mathrm{co}}\right)^{3.8}}$
(55)

where $\kappa_{\phi}$ represents the fraction of vacuum energy converted into
the gradient energy of the scalar field, and $H_{*}$ is the Hubble parameter
at $T_{n}$. The peak frequency of the GW spectrum generated by bubble
collisions is given by

$f_{\mathrm{co}}\simeq 1.65\times
10^{-5}\mathrm{~{}Hz}\frac{(8\pi)^{1/3}}{H_{*}R}\left(\frac{0.62v_{w}}{1.8-0.1v_{w}+v_{w}^{2}}\right)\left(\frac{T_{n}}{100~{}\mathrm{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}$
(56)

  * •

Sound wave
The expansion of bubbles leads to the formation of sound waves in the
surrounding fluid. Based on the numerical simulation, we obtain
Hindmarsh:2013xza ; Hindmarsh:2015qta

$\displaystyle h_{100}^{2}\Omega_{sw}(f)\simeq$ $\displaystyle 2.65\times
10^{-6}{\beta}^{-1}\left(\frac{\kappa_{sw}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{1/3}v_{w}\left(f/f_{sw}\right)^{3}\left(\frac{7}{4+3\left(f/f_{sw}\right)^{2}}\right)^{7/2}$
(57)

where $\kappa_{sw}$ represents the fraction of vacuum energy that is converted
into the kinetic energy of the fluid’s bulk motion and the peak frequency of
the GW spectrum generated by sound wave is given by

$f_{\mathrm{sw}}\simeq 2.6\times
10^{-5}\mathrm{~{}Hz}\frac{1}{H_{*}R}\left(\frac{T_{n}}{100~{}\mathrm{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}$
(58)

  * •

Turbulence
As the fluid propagates, the complex physical process of turbulence can be
triggered. Based on the numerical simulation, we obtain Caprini:2009yp ;
RoperPol:2019wvy

$h_{100}^{2}\Omega_{\mathrm{turb}}(f)\simeq 1.14\times
10^{-4}H_{*}R\left(\frac{\kappa_{\mathrm{turb}}\alpha}{1+\alpha}\right)^{3/2}\left(\frac{100}{g_{*}}\right)^{1/3}\frac{\left(f/f_{\text{turb
}}\right)^{3}}{\left(1+f/f_{\text{turb }}\right)^{11/3}\left(1+8\pi
f/H_{*}\right)}$ (59)

where $\kappa_{\mathrm{turb}}$ represents how effectively vacuum energy is
converted into turbulent flow

$\kappa_{\text{turb }}=\tilde{\epsilon}\kappa_{sw}$ (60)

We chose $\tilde{\epsilon}=0.1$ and the peak frequency of the GW spectrum
generated by turbulence is given by

$f_{\text{turb }}\simeq 7.91\times
10^{-5}\mathrm{~{}Hz}\frac{1}{H_{*}R}\left(\frac{T_{*}}{100~{}\mathrm{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}$
(61)

Table 1: Four sets benchmark parameters of DFSZ-I model  |  $\mu_{3}$ [GeV] | $\lambda_{3}$ | $\lambda_{13}$ | $\lambda_{23}$ |  $\mu$ [GeV] |  $T_{n}$ [GeV] |  $v_{\sigma}(T_{n})$ [GeV] | $\alpha$ | $\beta^{*}$
---|---|---|---|---|---|---|---|---|---
$BP_{1}$ | $7.13\times 10^{5}$ | $1.19\times 10^{-4}$ | 1.91 | 0.70 | $1.12\times 10^{6}$ | $1.34\times 10^{7}$ | $3.31\times 10^{7}$ | 0.87 | 330
$BP_{2}$ | $6.13\times 10^{6}$ | $1.12\times 10^{-3}$ | 1.79 | 0.59 | $2.69\times 10^{6}$ | $1.02\times 10^{7}$ | $2.55\times 10^{7}$ | 0.81 | 586
$BP_{3}$ | $7.13\times 10^{5}$ | $1.66\times 10^{-3}$ | 0.91 | 0.70 | $3.15\times 10^{6}$ | $9.76\times 10^{6}$ | $4.06\times 10^{7}$ | 1.63 | 579
$BP_{4}$ | $6.13\times 10^{5}$ | $1.20\times 10^{-4}$ | 1.90 | 0.60 | $1.35\times 10^{6}$ | $1.36\times 10^{7}$ | $3.50\times 10^{7}$ | 0.67 | 936

To find the parameter space that can generate a strong phase transition, we
performed numerical calculations on the phase transition dynamics by using the
CosmoTransitions package Wainwright:2011kj . We aim for the phase transition
to occur at a high energy scale to generate strong GW signals. So during the
calculations, we imposed the following optimization conditions on the
parameters: the critical temperature $T_{c}$ should be as high as possible and
at least greater than the electroweak scale $100~{}\mathrm{GeV}$; $v_{\sigma}$
should significantly exceed $246~{}\mathrm{GeV}$ since the vev is related to
the energy scale of the phase transition; and the ratio of the scalar field
mass to the nucleation temperature should be much smaller than 1 to ensure the
applicability of perturbation expansion in Eq. (43). By systematically
scanning the parameter space and obtaining results that maximize these
conditions, we obtained several parameter combinations that have a high-energy
scale. We present four benchmark parameter sets ($BP_{1}$, $BP_{2}$, $BP_{3}$,
and $BP_{4}$) for the DFSZ-I model in Tab. 1. These sets are selected based on
two main criteria: (1) consistency with the current experimental constraints
on the axion-photon coupling, $g_{a\gamma\gamma}$, and (2) the ability to
generate a SFOPT with high-energy scale. These parameter sets will be used to
calculate the GW spectra. Then be compared with the sensitivity curves of the
CE and ET detectors to determine whether these detectors can observe the GW
signals predicted by this model. During the calculations, we find some
parameters have a relatively smaller impact on the resulting phase transition
GW signals. For simplicity, these model parameters $\lambda_{1}$,
$\lambda_{4}$, and $\lambda_{12}$ are fixed at 0.0001, $\lambda_{2}$ is set to
0.1, and $\mu_{1}$ and $\mu_{2}$ are fixed at $6\times 10^{5}$ GeV. The
remaining model parameters, including $\mu_{3}$, $\lambda_{3}$, $\mu$,
$\lambda_{13}$, and $\lambda_{23}$, are varied to generate different phase
transition parameters.

We obtained phase transitions with energy scales up to
$f_{PQ}=v_{\sigma}(T_{n})=40~{}\mathrm{PeV}$, where $f_{PQ}$ is the PQ
symmetry-breaking scale. Since the axion mass is inversely proportional to
$f_{PQ}$ Gorghetto:2018ocs , the high-energy scale result indicates that it
may be possible for invisible axions to exist in the DFSZ-I model.

Table 2: Axion mass, $g_{a\gamma\gamma}$ for different $f_{PQ}$ values. |  $f_{PQ}$ [GeV] |  $m_{a}$ [eV] |  $g_{a\gamma\gamma}$ [GeV-1]
---|---|---|---
$BP_{1}$ | $3.31\times 10^{7}$ | $0.17$ | $5.17\times 10^{-10}$
$BP_{2}$ | $2.55\times 10^{7}$ | $0.22$ | $6.71\times 10^{-10}$
$BP_{3}$ | $4.06\times 10^{7}$ | $0.14$ | $4.22\times 10^{-10}$
$BP_{4}$ | $3.50\times 10^{7}$ | $0.16$ | $4.89\times 10^{-10}$
Figure 1: The $m_{a}$ and $g_{a\gamma\gamma}$ obtained from all four benchmark
points in the ALP-photon coupling exclusion diagram from OHare2020 .

$m_{a}=5.691(51)\left(\frac{10^{9}~{}\mathrm{GeV}}{f_{PQ}}\right)~{}\mathrm{meV}$
(62)

The numbers in parentheses represent the uncertainty in the last two digits of
the main value. Given the $m_{a}$, we can calculate the axion-photon coupling
constant $g_{a\gamma\gamma}$. The calculation as follows
GrillidiCortona:2015jxo

$g_{a\gamma\gamma}=\frac{\alpha_{c}}{2\pi
f_{PQ}}\left(\frac{E}{N}-1.92(4)\right)=\left(0.203(3)\frac{E}{N}-0.39(1)\right)\frac{m_{a}}{\mathrm{GeV}^{2}}$
(63)

where $\alpha_{c}$ is the fine structure constant, $E/N$ is the ratio of the
electromagnetic anomaly coefficient to the colour anomaly coefficient. For the
DFSZ-I model, $E/N$ is $8/3$ Zhitnitsky:1980tq ; Dine:1981rt .

Figure 1 shows the ALP-photon coupling exclusion plot. The coloured regions
represent the excluded parameter space from various ALP-photon coupling
experiments, such as CAST, ADMX, etc. The yellow band denotes the
theoretically predicted range for axions. The four elliptical points
$BP_{1-4}$ in the plot represent four benchmark parameter sets in Tab.1, which
are shown in red, blue, green, and black, respectively. The corresponding
$m_{a}$ and $g_{a\gamma\gamma}$ values for these points are listed in Table 2.
The $m_{a}$ for $BP_{1}$, $BP_{2}$, $BP_{3}$, and $BP_{4}$ is
$0.17~{}\mathrm{eV}$, $0.22~{}\mathrm{eV}$, $0.14~{}\mathrm{eV}$, and
$0.16~{}\mathrm{eV}$, respectively. All four points fall within the expected
yellow region showing that these points satisfy the constraints from existing
experiments. This result identifies parameter sets that are consistent with
experimental constraints and the ability to achieve high-energy scale SFOPT.

Figure 2: GW spectra for DFSZ-I model. The light-green and orange curves
depict the sensitivity curves with SNR=1 for ET and CE, respectively. The red,
blue, green, and black curves represent the phase transition GW signals for
$BP_{1}$, $BP_{2}$, $BP_{3}$, and $BP_{4}$ with $v_{w}=0.9$. The relevant
parameters are summarized in Tab. 1.

Figure 2 shows the GW spectrum for four benchmark parameter sets. The light-
green and orange curves represent the sensitivity curves for ET and CE with
SNR=1 (more details of the sensitivity curves for detectors will be introduced
in the next section), respectively, over a running time of $1\times 10^{8}$ s.
The red, blue, green, and black curves represent the GW signals for $BP_{1}$,
$BP_{2}$, $BP_{3}$, and $BP_{4}$ with the $v_{w}=0.9$. The SNR values using
the definition of Eq. (27) for each benchmark point are

$SNR:\quad BP_{1}(93.26)\quad BP_{2}(14.34)\quad BP_{3}(29.12)\quad
BP_{4}(0.77)$ (64)

Table 3 shows that the SNR of $BP_{4}$ is only 0.77, which is significantly
below the threshold value SNR=8 Hall:2020dps  for the CE detector, making its
signal difficult to capture by the CE detector. $BP_{1}$ has an SNR of 93.26,
$BP_{2}$ has an SNR of 14.34 and $BP_{3}$ has an SNR of 29.12. The SNR values
of these three curves exceed the threshold, indicating that the CE detector
can detect them. The detectability of $BP_{1}$, $BP_{2}$, and $BP_{3}$
indicates that the DFSZ-I model can generate observable axion phase transition
GW signals.

Figure 3: GW spectra for DFSZ-I model with $v_{w}=0.92$ and $v_{w}=0.8$. The
light-green and orange curves depict the sensitivity curves with SNR=1 for ET
and CE. $BP_{1}$ (red solid line) and $BP_{1}$ (red dash-dotted line)
correspond to the expected GW signals of $v_{w}=0.92$ and $v_{w}=0.8$,
respectively. $BP_{3}$ (green solid line) and $BP_{3}$ (green dash-dotted
line) correspond to the expected GW signals of $v_{w}=0.92$ and $v_{w}=0.8$,
respectively. The relevant parameters are summarized in Tab. 1.

To investigate the impact of the unknown parameter $v_{w}$ on the GW spectrum,
we consider two representative values which may produce stronger GW signals
during the phase transition process: $v_{w}=0.92$ and $v_{w}=0.8$. Fig. 3
shows the GW spectrum for two benchmark parameter sets $BP_{1}$ and $BP_{3}$
that have higher SNR with different $v_{w}$. $BP_{1}$ (red solid line) and
$BP_{1}$ (red dash-dotted line) correspond to the expected GW signals of
$v_{w}=0.92$ and $v_{w}=0.8$, respectively. $BP_{3}$ (green solid line) and
$BP_{3}$ (green dash-dotted line) correspond to the expected GW signals of
$v_{w}=0.92$ and $v_{w}=0.8$, respectively. We find that there are differences
in the frequency and strength of the GW spectra under different $v_{w}$ values
in Fig. 3. In order to assess the constraining power of CE on these phase
transition parameters more precisely, we need to quantify the uncertainties of
each parameter. FM analysis is a powerful statistical tool which is widely
used for estimating the precision of parameter estimates. In the next section,
we will use FM analysis to calculate the relative uncertainties of the phase
transition parameters in the DFSZ-I model and evaluate the constraining power
of CE on these parameters.

Table 3: SNR values for two benchmark parameter sets at three different bubble wall velocities $v_{w}$. | $v_{w}=0.92$ | $v_{w}=0.9$ | $v_{w}=0.8$
---|---|---|---
$BP_{1}$ | 102.26 | 93.26 | 56.83
$BP_{3}$ | 31.98 | 29.12 | 17.60

##  V Fisher matrix analysis of the phase transition parameter at Cosmic
Explorer

FM is a powerful statistical tool which is widely used for estimating the
precision of parameter estimates in various fields, including GW astronomy
Fisher:1922saa . For a probability distribution $f(x;\Theta),$ of an
observable random variable $x$ that depends on a set of unknown parameters
$\Theta=(\Theta_{1},\dots,\Theta_{k})$, the FM is defined as the expected
value of the second-order partial derivatives of the log-likelihood function
with respect to the parameters. One of the critical properties of the FM is
that the inverse of the FM corresponds to the covariance matrix of the
observable parameters as the probability distribution approximates a Gaussian
distribution Hamimeche:2008ai . The diagonal elements of the inverse FM are
the variances (squared uncertainties) of the parameter estimates. Therefore,
smaller values of these diagonal elements indicate more precise parameter
estimates with lower uncertainties. We will use FM analysis to calculate the
relative uncertainties of phase transition parameters, which enables us to
assess the constraining power of detectors such as the CE and ET on these
parameters. To apply FM analysis, we first need to construct a likelihood
model for CE.

###  V.1 Likelihood model

In our study, the observational data is the GW signal detected by CE, and the
theoretical model parameters are the dynamical parameters of the phase
transition generated by the DFSZ-I axion model. In the following discussion,
we will introduce in detail how to construct the likelihood model using CE, in
preparation for the subsequent FM analysis.

The GW strain in the frequency domain, $\tilde{h}_{gw}\left(f_{n}\right)$, is
given by the discrete Fourier transform of the time-domain strain, $h_{gw}(t)$
Gowling:2021gcy

$\tilde{h}_{gw}\left(f_{n}\right)=\frac{1}{\sqrt{N_{t}}}\sum_{m=0}^{N_{t}-1}h_{gw}(t_{m})\exp\left(-2\pi
if_{n}t_{m}\right)$ (65)

where $N_{t}=T_{tot}/T_{s}$ is the total number of samples, $T_{tot}=1\times
10^{8}\mathrm{~{}s}$ is the total observation time for the CE mission, and
$T_{s}$ is the data sampling period. The time $t_{m}$ is given by
$t_{m}=mT_{s}$, where $m$ is an integer ranging from $-N_{t}/2$ to $N_{t}/2$.
The frequency $f_{n}$ is defined as

$f_{n}=n/T_{\text{tot }}$ (66)

where $n$ is an integer ranging from $-N_{t}/2$ to $N_{t}/2$. The GW energy
density spectrum, $\Omega_{\mathrm{gw}}\left(f_{n}\right)$, is related to the
strain in the frequency domain by

$\Omega_{\mathrm{gw}}\left(f_{n}\right)=\left(\frac{4\pi^{2}}{3H_{0}^{2}}\right)f_{n}^{2}\left|\tilde{h}_{gw}\left(f_{n}\right)\right|^{2}$
(67)

The PSD of the CE detector can be expressed as $E_{n}$ and the variances of
the Fourier amplitudes are denoted as $A_{n}$. To optimize computational
efficiency, we employ a simplified approach for data grouping known as
frequency binning. This involves dividing the frequency range into $N_{b}$
bins, where the spacing between consecutive bins follows a logarithmic scale.
Consequently, each bin contains a specific number of frequencies, denoted as
$n_{b}$

$n_{b}=\left[\left(f_{b}-f_{b-1}\right)T_{\text{tot }}\right],0\leq b\leq
N_{b}$ (68)

where $f_{b}$ and $f_{b-1}$ are the upper and lower frequency boundaries of
the $b$-th bin. The relationship between $A_{n}$ and $E_{n}$ can be expressed
as

$\bar{E}_{n}=\frac{A_{b}}{n_{b}}\sum_{n\in I_{b}}\frac{E_{n}}{A_{n}}$ (69)

with

$\frac{1}{A_{b}}=\frac{1}{n_{b}}\sum_{n\in{I_{b}}}\frac{1}{A_{n}}$ (70)

$A_{b}$ is the mean value of the distribution. Let $I_{b}$ denote a group of
integers numbers. Consider a set of independent normal random variables
$A_{1},\dots,A_{n}$, the variable

$\sum_{n\in I_{b}}\frac{E_{n}}{A_{n}}=E_{1}/A_{1}+...+E_{n}/A_{n}$ (71)

approximates a Gamma distribution with shape parameter $2n_{b}$ when $n_{b}$
is sufficiently large, according to the Central Limit Theorem. This
distribution arises from scaling the Chi-square distribution

$p\left(\bar{E}_{b}\mid
A_{b}\right)=\prod_{b=1}^{N_{b}}\frac{1}{\left(n_{b}-1\right)!}\frac{n_{b}}{A_{b}}\left(n_{b}\frac{\bar{E}_{b}}{A_{b}}\right)^{n_{b}-1}\exp\left(-n_{b}\frac{\bar{E}_{b}}{A_{b}}\right)$
(72)

Under the Gaussian approximation, we have

$p\left(\bar{E}_{b}\mid
A_{b}\right)=\prod_{b=1}^{N_{\mathrm{b}}}\left(\frac{n_{b}}{2\pi
A_{b}^{2}}\right)^{\frac{1}{2}}\exp\left(-\frac{1}{2}\frac{n_{b}\left(\bar{E}_{b}-A_{b}\right)^{2}}{A_{b}^{2}}\right)$
(73)

From Eq. (91) in the Appendix, the natural logarithm of the likelihood
function with the Gaussian approximation is

$l_{G}=\ln(p)=-\frac{1}{2}\sum_{b=1}^{N_{\mathrm{b}}}\frac{n_{b}\left(\bar{E}_{b}-A_{b}\right)^{2}}{A_{b}^{2}}-\sum_{b=1}^{N_{\mathrm{b}}}\ln
2A_{b}+\text{ const. }$ (74)

The FM can be obtained from Eq. (97) is

$\displaystyle F_{ij}^{G}=$ $\displaystyle=E\left[\frac{\partial
l_{G}}{\partial\Theta_{i}}\frac{\partial l_{G}}{\partial\Theta_{j}}\right]$
(75) $\displaystyle=E\left[\frac{\partial l_{G}}{\partial A_{b}}\frac{\partial
l_{G}}{\partial A_{b}}\frac{\partial A_{b}}{\partial\Theta_{i}}\frac{\partial
A_{b}}{\partial\Theta_{j}}\right]$
$\displaystyle=E\left[\frac{\partial^{2}}{\partial^{2}A_{b}}\left(-\frac{1}{2}\sum_{b=1}^{N_{b}}\frac{n_{b}\left(\bar{E}_{b}-A_{b}\right)}{A_{b}^{2}}-\sum_{n=1}^{N_{b}}\ln
2A_{b}\right)\right]\frac{\partial A_{b}}{\partial\Theta_{i}}\frac{\partial
A_{b}}{\partial\Theta_{j}}$
$\displaystyle=\sum_{b=1}^{N_{\mathrm{b}}}\frac{n_{b}}{A_{b}^{2}}\frac{\partial
A_{b}}{\partial\Theta_{i}}\frac{\partial A_{b}}{\partial\Theta_{j}}$
$\displaystyle=T_{\text{obs }}\sum_{b=1}^{N_{\mathrm{b}}}\frac{\Delta
f_{b}}{\Omega_{\mathrm{t}}^{2}}\frac{\partial\Omega_{\mathrm{t}}}{\partial\Theta_{i}}\frac{\partial\Omega_{\mathrm{t}}}{\partial\Theta_{j}}$

To establish the theoretical model for the data, we will employ the total
power spectrum $\Omega_{\mathrm{t}}$ (defined in Eq. (84)) instead of the
spectral densities $A_{b}$. The covariance matrix is the inverse of the FM,
and the square roots of the diagonal entries correspond to the relative
uncertainties of the parameters.

###  V.2  Cosmic Explorer noise

The CE and the ET are third-generation ground-based GW detectors, anticipated
to commence operations post-2030 LIGOScientific:2016wof ; Reitze:2019iox ;
Hall:2022dik ; Hall:2020dps ; Sathyaprakash:2012jk ; Maggiore:2019uih . In
comparison to second-generation detectors, their sensitivity is expected to
increase by an order of magnitude, enabling the detection of frequencies as
low as several hertz. This sensitivity enhancement will allow for the
observation of $\mathcal{O}(10^{5}\text{-}10^{6})$ compact binary coalescence
events (CBCs) within a year-long observation period, not limited to binary
black hole (BBH) and binary neutron star (BNS) mergers. As shown in Fig. 2, CE
can detect GW signals from SFOPT in axion models. Therefore, we choose CE as
an example to perform FM analysis. To better understand the noise model of CE,
we need to carefully analyze various noise sources and their impact on the
detector’s sensitivity.

####  V.2.1 shot noise

The first noise source we consider is the shot noise of the laser, stemming
from the quantum nature of photons. In GW detectors, shot noise limits the
sensitivity at high frequencies. One method to mitigate the impact of shot
noise on experimental results is increasing the number of photons received by
the detector, achievable by amplifying the power of the beam, $P_{\text{bs}}$.
The initial beam power, $P_{0}$, refers to the power of the original beam from
the laser source incident on the interferometer. A beam splitter is an optical
component that divides an incoming beam into two parts, and $C$ represents the
multiplier by which the beam power is increased before reaching the beam
splitter. Therefore, $P_{\text{bs}}\equiv CP_{0}=1650~{}\mathrm{kw}$ denotes
the power measured at the beam splitter after power amplification.

The shot noise is formulated as Maggiore:2007ulw

$\left.S_{n}^{1/2}(f)\right|_{\text{shot
}}=\frac{1}{8\mathcal{F}L}\left(\frac{4\pi\hbar\lambda_{\mathrm{L}}c}{\eta
P_{\mathrm{bs}}}\right)^{1/2}\sqrt{1+\left(f/f_{p}\right)^{2}}$ (76)

where $\mathcal{F}$ denotes the Finesse of the Fabry-Perot cavities, a measure
of the optical quality and reflectivity of the mirrors within the Fabry-Perot
interferometer, set at $\mathcal{F}=200$. $L$ is arm lengths of
$40~{}\mathrm{km}$ and $\eta=0.96$ is the efficiency of the photodetector. The
pole frequency, $f_{p}=500~{}\mathrm{Hz}$, is a crucial parameter describing
the frequency response of a GW detector, especially those utilizing Fabry-
Perot cavities. The symbol $\lambda_{\mathrm{L}}$ denotes the laser wavelength
used in the interferometer, which is approximately $1550~{}\mathrm{nm}$.

####  V.2.2 radiation pressure

In GW detection, enhancing the experimental signal through increased laser
beam power is a common practice Caves:1980pha . However, while this approach
improves detector sensitivity, it also introduces additional noise in the form
of radiation pressure. This pressure is the movement of mirrors as the beam
traverses its path and reflects, exerting a force on the mirrors in the
system. At high frequencies, shot noise typically dominates over radiation
pressure noise, while at low frequencies, radiation pressure noise becomes
more significant. The strain sensitivity due to radiation pressure typically
takes the following form Maggiore:2007ulw

$\left.S_{n}^{1/2}(f)\right|_{\mathrm{rad}}=\frac{16\sqrt{2}\mathcal{F}}{ML(2\pi
f)^{2}}\sqrt{\frac{\hbar}{2\pi}\frac{P_{\mathrm{bs}}}{\lambda_{L}c}}\frac{1}{\sqrt{1+\left(f/f_{p}\right)^{2}}}$
(77)

where $M=320~{}\text{kg}$ represents the mass of the mirror.

Combining the shot noise and radiation pressure, one obtains the optical
readout noise with its PSD given by

$\left.S_{n}(f)\right|_{\text{opt }}=\left.S_{n}(f)\right|_{\text{shot
}}+\left.S_{n}(f)\right|_{\text{rad }}$ (78)

So that the energy density is

$\Omega_{\mathrm{opt}}=\left(\frac{4\pi^{2}}{3H_{0}^{2}}\right)f^{3}S_{n}(f)|_{\text{opt
}}$ (79)

####  V.2.3 Seismic noise

Seismic noise arises from the continuous movement of the ground, with
amplitudes around several micrometres. The strain sensitivity due to seismic
noise typically can be represented as Maggiore:2007ulw

$\frac{x(f)}{L}\simeq\frac{A}{L}\left(\frac{1\mathrm{~{}Hz}}{f}\right)^{\nu}\mathrm{~{}m}\cdot\mathrm{Hz}^{-1/2}$
(80)

Above $1\mathrm{~{}Hz}$, the exponent $\nu\simeq 2$, with the magnitude of $A$
around $10^{-7}$. Dividing $x(f)$ by the length $L$, we obtain the noise
strain sensitivity. The energy density is

$\Omega_{\mathrm{se}}=\left(\frac{4\pi^{2}}{3H_{0}^{2}}\right)f^{3}\left(\frac{x(f)}{L}\right)^{2}$
(81)

Seismic noise predominantly governs the lower frequency domain, while at
higher frequencies, shot noise and radiation pressure noise become
increasingly significant. The total power spectrum energy density of the CE
noise model is

$\Omega_{\mathrm{n}}=\Omega_{\mathrm{opt}}+\Omega_{\mathrm{se}}$ (82)

####  V.2.4 Compact binary coalescence

The GW power spectrum of CBCs foreground, such as BBH and BNS
LIGOScientific:2017zlf ; Rosado:2011kv ; Sachdev:2020bkk  is often defined as

$\Omega_{\mathrm{cbc}}(f)=\Omega_{\mathrm{ref}}\left(\frac{f}{f_{\mathrm{ref}}}\right)^{2/3}$
(83)

where $\Omega_{\mathrm{ref}}$ represents the energy spectral density at a
reference frequency $f_{\mathrm{ref}}$. Due to the uncertain spatial
distribution of BBH and BNS sources, $\Omega_{\mathrm{ref}}$ is uncertain,
typically ranging from $10^{-7}$ to $10^{-10}$ Sachdev:2020bkk . For this
analysis, we assume $\Omega_{\mathrm{ref}}=10^{-10}$. For BBH mergers, which
typically occur within a frequency range from a few hertz to several hundred
hertz, we set a reasonable reference frequency
$f_{\mathrm{ref}}=200~{}\mathrm{Hz}$. This frequency is indicative of the
coalescence phase during the BBH merger process, where the GW signal is at its
strongest. For BNS mergers, which happen at slightly higher frequencies
compared to BBH mergers, a higher reference frequency of $1000~{}\mathrm{Hz}$
is selected.

The total energy density is

$\Omega_{\mathrm{t}}=\Omega_{\mathrm{n}}+\Omega_{\mathrm{cbc}}+\Omega_{\mathrm{gw}}$
(84)

###  V.3 The relative uncertainties

After obtaining Eq. (84) for the CE detector, we can substitute it into Eq.
(75) to calculate the relative uncertainty of specific GW model parameters. We
choose the sound wave mechanism, which dominates the GW spectrums. The
considered parameters include the phase transition strength, Hubble-scaled
mean bubble separation $R_{H}=H_{*}R$, and nucleation temperature. These
parameters can be encapsulated in a vector $\vec{\Theta}$, defined as
$\vec{\Theta}=\left(v_{\mathrm{w}},\alpha,R_{H},\left(T_{\mathrm{n}}/~{}\mathrm{GeV}\right)\right)$.
To calculate the relative uncertainties of the GW spectrum generated by the
sound wave mechanism, we rewrite Eq. (57) in the following form to include the
$R_{H}$

$h_{100}^{2}\Omega_{sw}(f)\cong 1.64\times
10^{-6}\frac{R_{H}^{2}}{\bar{U}_{f}}\left(\frac{\kappa_{v}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{1/3}\left(\frac{f}{f_{sw}}\right)^{3}\left(\frac{7}{4+3\left(f/f_{sw}\right)^{2}}\right)^{7/2}$
(85)

$f_{sw}\cong 2.6\times
10^{-5}\mathrm{~{}Hz}\frac{1}{R_{H}}\left(\frac{T_{n}}{100~{}\mathrm{Gev}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}$
(86)

$\bar{U}_{f}^{2}\approx\frac{3}{4}\frac{\kappa_{v}\alpha}{1+\alpha}$ (87)

Figure 4: Relative uncertainty with phase transition dynamics parameters wall
speed $v_{w}$, transition strength $\alpha$, Hubble-scaled mean bubble spacing
$R_{H}$, and nucleation temperature $T_{n}$.

We have selected a set of phase transition parameters around $BP_{1}$ from
Tab.1 that can be detected by CE. We set $T_{n}$ to be $1.34\times
10^{7}~{}\mathrm{GeV}$ and $R_{H}=0.1$. The parameter space is sampled by
varying $v_{w}$ from 0.92 to 0.94 and $\alpha$ from 0.58 to 0.65.

Figure 4 contains four subplots, each showing the relative uncertainties of
$v_{w}$, $\alpha$, $R_{H}$, and $T_{n}$ within the given ranges of $v_{w}$ and
$\alpha$: $\Delta v_{w}/v_{w}$, $\Delta\alpha/\alpha$, $\Delta R_{H}/R_{H}$,
and $\Delta T_{n}/T_{n}$. A colour bar on the right side of the figure
indicates the magnitude of the relative uncertainties. From dark blue to light
green, the relative uncertainties gradually increase. It can be seen from the
figure that the relative uncertainty of $v_{w}$ is the smallest, as the colour
in the $v_{w}$ subplot is the most blue-shifted. These results indicate that
detectors like CE will be most sensitive to and place the strongest
constraints on measurements of the $v_{w}$. The results obtained through the
FM analysis are consistent with the findings presented in Fig. 3.

##  VI Discussions and Conclusion

We have investigated GW detection of the DFSZ axion model at CE. Taking the
popular DFSZ-I axion model as a typical example, we precisely calculate the
phase transition dynamics of the U(1) PQ symmetry-breaking process and its
associated phase transition GW. After scanning different model parameter
combinations in the DFSZ-I model, we find that the U(1) PQ symmetry-breaking
can be the cosmological SFOPT at high energy scales up to $40~{}\mathrm{PeV}$
that current experiments are unable to reach. We calculate the axion mass and
axion-photon coupling constant predicted by the DFSZ model and compare them
with current experimental constraints to find the benchmark points. We find
that the parameters can not only generate a SFOPT with a high-energy scale but
also are consistent with the current experimental constraints on the axion-
photon coupling. Then, we quantified the detectability of the DFSZ-I model at
the CE and ET by precisely calculating the SNR. For $v_{w}=0.9$, three signals
achieve $SNR>8$ at CE, indicating that they can be observed by CE.

Considering the uncertainty of GW spectra from the bubble wall velocity, we
provide GW spectra under different bubble wall velocities as examples. The
results show that the GW spectrum is sensitive to the value of the bubble wall
velocity. Based on this, to further quantify the uncertainties of phase
transition parameters at future GW detectors, we perform FM analysis on the
key phase transition parameters, including $v_{w}$, $\alpha$, $R_{H}$, and
$T_{n}$. The results show that the relative uncertainty of $v_{w}$ is smaller
than the relative uncertainties of other parameters. Thus, if the predicted
phase transition GW spectra are observed at detectors like CE, bubble wall
velocity $v_{w}$ will be the first parameter to be measured or most strongly
constrained.

Looking ahead, the construction and operation of next-generation GW detectors
such as ET and CE are expected to improve the ability to detect phase
transition GW. This study demonstrates that CE detectors have the potential to
detect GWs from the U(1) PQ symmetry-breaking process, providing a novel
approach to probe the existence of axion. More comprehensive studies in
various axion models are left in our future studies.

###### Acknowledgements.

We thank Seokhoon Yun for the helpful discussions on axion constraints and
Siyu Jiang for making useful suggestions on the manuscript. This work is
supported by the National Natural Science Foundation of China (NNSFC) under
Grant No.12205387. This work is partly supported in part by Guangdong Major
Project of Basic and Applied Basic Research (Grant No. 2019B030302001).

##  Appendix A The field-dependent mass

To determine the field-dependent mass of each physical particle, we need to
diagonalize the mass-squared matrix. The elements of the mass-squared matrix
$M$ for the bosons can be expressed as follows

$\displaystyle M_{00}$
$\displaystyle=v_{\sigma}^{2}\lambda_{13}-2\mu_{1}^{2},$ $\displaystyle
M_{11}$ $\displaystyle=\frac{v_{\sigma}^{2}\lambda_{13}}{2}-\mu_{1}^{2},$
$\displaystyle M_{22}$
$\displaystyle=\frac{v_{\sigma}^{2}\lambda_{13}}{2}-\mu_{1}^{2}$
$\displaystyle M_{33}$
$\displaystyle=v_{\sigma}^{2}\lambda_{23}-2\mu_{2}^{2},$ $\displaystyle
M_{44}$ $\displaystyle=\frac{v_{\sigma}^{2}\lambda_{23}}{2}-\mu_{2}^{2},$
$\displaystyle M_{55}$
$\displaystyle=\frac{v_{\sigma}^{2}\lambda_{23}}{2}-\mu_{2}^{2}$
$\displaystyle M_{66}$ $\displaystyle=3v_{\sigma}^{2}\lambda_{3}-\mu_{3}^{2},$
$\displaystyle M_{77}$ $\displaystyle=v_{\sigma}^{2}\lambda_{3}-\mu_{3}^{2},$
$\displaystyle M_{88}$ $\displaystyle=v_{\sigma}^{2}\lambda_{13}-2\mu_{1}^{2}$
$\displaystyle M_{99}$
$\displaystyle=v_{\sigma}^{2}\lambda_{23}-2\mu_{2}^{2},$ $\displaystyle
M_{03}$ $\displaystyle=-v_{\sigma}^{2}\lambda_{5},$ $\displaystyle M_{30}$
$\displaystyle=-v_{\sigma}^{2}\lambda_{5}$ $\displaystyle M_{14}$
$\displaystyle=\frac{v_{\sigma}^{2}\lambda_{5}}{2},$ $\displaystyle M_{41}$
$\displaystyle=\frac{v_{\sigma}^{2}\lambda_{5}}{2},$ $\displaystyle M_{25}$
$\displaystyle=-\frac{v_{\sigma}^{2}\lambda_{5}}{2}$ $\displaystyle M_{52}$
$\displaystyle=-\frac{v_{\sigma}^{2}\lambda_{5}}{2},$ $\displaystyle M_{89}$
$\displaystyle=\frac{v_{\sigma}^{2}\lambda_{5}}{2},$ $\displaystyle M_{98}$
$\displaystyle=\frac{v_{\sigma}^{2}\lambda_{5}}{2}$

Let

$G_{44}=M_{44}+\Pi_{S}^{44},\quad G_{55}=M_{55}+\Pi_{S}^{55},\quad
G_{66}=M_{66}+\Pi_{S}^{66}$ (88) $G_{77}=M_{77}+\Pi_{S}^{77},\quad
G_{88}=M_{88}+\Pi_{S}^{88},\quad G_{99}=M_{99}+\Pi_{S}^{99}$ (89)

The mass square matrix $M$ with the daisy corrections is

$M=\left(\begin{array}[]{cccccccccc}M_{00}+\Pi_{S}^{00}&0&0&M_{03}&0&0&0&0&0&0\\\
0&M_{11}+\Pi_{S}^{11}&0&0&M_{14}&0&0&0&0&0\\\
0&0&M_{22}+\Pi_{S}^{22}&0&0&M_{25}&0&0&0&0\\\
M_{30}&0&0&M_{33}+\Pi_{S}^{33}&0&0&0&0&0&0\\\ 0&M_{41}&0&0&G_{44}&0&0&0&0&0\\\
0&0&M_{52}&0&0&G_{55}&0&0&0&0\\\ 0&0&0&0&0&0&G_{66}&0&0&0\\\
0&0&0&0&0&0&0&G_{77}&0&0\\\ 0&0&0&0&0&0&0&0&G_{88}&M_{89}\\\
0&0&0&0&0&0&0&0&M_{98}&G_{99}\end{array}\right)$ (90)

##  Appendix B Efficiency coefficients

TABLE 4 shows the relationship between efficiency coefficients for different
velocities.

Table 4: Efficiency coefficients for different velocities. velocities | efficiency coefficients
---|---
$v_{w}\ll c_{s}$ | $\kappa_{A}\simeq v_{w}^{6/5}\frac{6.9\alpha}{1.36-0.037\sqrt{\alpha}+\alpha}$
$v_{w}=c_{s}$ | $\kappa_{B}\simeq\frac{\alpha^{2/5}}{0.017+\left(0.997+\alpha\right)^{2/5}}$
$v_{w}=c_{J}$ |  $\kappa_{C}\simeq\frac{\sqrt{\alpha}}{0.135+\sqrt{0.98+\alpha}}\quad$ and $\quad c_{J}=\frac{\sqrt{\frac{2}{3}\alpha+\alpha^{2}}+\sqrt{1/3}}{1+\alpha}$
$v_{w}\rightarrow 1$ | $\kappa_{D}\simeq\frac{\alpha}{0.73+0.083\sqrt{\alpha}+\alpha}$
$v_{w}\lesssim c_{s}$ | $\kappa\left(v_{w}\lesssim c_{s}\right)\simeq\frac{c_{s}^{11/5}\kappa_{A}\kappa_{B}}{\left(c_{s}^{11/5}-v_{w}^{11/5}\right)\kappa_{B}+v_{w}c_{s}^{6/5}\kappa_{A}}$
$v_{w}\gtrsim c_{J}$ | $\kappa\left(v_{w}\gtrsim c_{J}\right)\simeq\frac{\left(c_{J}-1\right)^{3}c_{J}^{5/2}v_{w}^{-5/2}\kappa_{C}\kappa_{D}}{\left[\left(c_{J}-1\right)^{3}-\left(v_{w}-1\right)^{3}\right]c_{J}^{5/2}\kappa_{C}+\left(v_{w}-1\right)^{3}\kappa_{D}}$
$v_{w}=c_{s}$ | $\delta\kappa\simeq-0.9\log\frac{\sqrt{\alpha}}{1+\sqrt{\alpha}}$
|
$\kappa\left(c_{s}<v_{w}<c_{J}\right)\simeq\kappa_{B}+\left(v_{w}-c_{s}\right)\delta\kappa+\frac{\left(v_{w}-c_{s}\right)^{3}}{\left(c_{J}-c_{s}\right)^{3}}\left[\kappa_{C}-\kappa_{B}-\left(c_{J}-c_{s}\right)\delta\kappa\right]$

##  Appendix C Fisher information

Once the probability density distribution function $f(\Theta;x)$ is known, the
logarithm of the likelihood function can be expressed as

$l(\Theta;x)=\ln[f(\Theta;x)]$ (91)

The score function is defined as the partial derivative of the natural
logarithm of the likelihood function

$\frac{\partial}{\partial\Theta}l(\Theta;x)$ (92)

Fisher information (FI) represents the second cross-moment of the score vector

$F(\Theta)=\mathrm{E}\left[\left(\frac{\partial}{\partial\Theta}l(\Theta;x)\right)\left(\frac{\partial}{\partial\Theta}l(\Theta;x)\right)^{T}\right]$
(93)

Given mild regularity conditions, the score function’s expected value equals
zero

$\mathrm{E}\left[\frac{\partial}{\partial\Theta}l(\Theta;x)\right]=0$ (94)

As a consequence, the FI is defined as the variance of the score

$\displaystyle F(\Theta)$
$\displaystyle=\mathrm{E}\left[\left(\frac{\partial}{\partial\Theta}l(\Theta;x)\right)\left(\frac{\partial}{\partial\Theta}l(\Theta;x)\right)^{T}\right]$
(95)
$\displaystyle=\mathrm{E}\left[\left\\{\frac{\partial}{\partial\Theta}l(\Theta;x)-\mathrm{E}\left[\frac{\partial}{\partial\Theta}l(\Theta;x)\right]\right\\}\left\\{\frac{\partial}{\partial\Theta}l(\Theta;x)-\mathrm{E}\left[\frac{\partial}{\partial\Theta}l(\Theta;x)\right]\right\\}\right]$
$\displaystyle=\operatorname{Var}\left[\frac{\partial}{\partial\Theta}l(\Theta;x)\right]$

It is assumed that the observations in the sample are assumed to be
independently and identically distributed. The maximum likelihood estimator,
represented by $\hat{\Theta}$, is approximately equal to the inverse of the
FI.

$\operatorname{Var}[\hat{\Theta}]\approx F(\Theta)^{-1}$ (96)

In matrix form, $\Theta$ is a vector with N parameters. The FI is a matrix
$F_{ij}$

$F_{ij}=E\left[\frac{\partial l(\Theta;x)}{\partial\Theta_{i}}\frac{\partial
l(\Theta;x)}{\partial\Theta_{j}}\right]$ (97)

and the inverse of FM is the covariance matrix $C_{ij}$

$C_{ij}=F_{ij}^{-1}$ (98)

##  Appendix D Notation and corresponding dimension of signal analysis.

TABLE 5 shows the notation of signal analysis and their dimension in natural
units.

Table 5: Notation of signal analysis variable | definition | Dimension[E]
---|---|---
$h_{ab}(t,\vec{x})$ | metric perturbation | 0
$h_{A}(f,\hat{k})$ | Fourier coefficients of metric perturbation | -1
$S_{h}(f)$ | strain PSD of a gravitational-wave background | -1
$\Omega_{\mathrm{gw}}(f)$ | fractional energy density spectrum of a GW background | 0
$h_{c}(f)$ | characteristic strain for GWs | 0
$h_{gw}(t)$ | detector response to GWs | 0
$n_{gw}(t)$ | detector response to noise | 0
$s_{gw}(t)$ | detector response | 0
$\tilde{h}_{gw}\left(f_{n}\right)$ | discrete Fourier transform of $h_{gw}(t)$ | 0
$\tilde{h}_{gw}(f)$ | Fourier transform of $h_{gw}(t)$ | -1
$\tilde{n}_{gw}(f)$ | Fourier transform of $n_{gw}(t)$ | -1
$S_{n}(f)$ | strain noise PSD for a single detector | -1
$h_{n}(f)$ | characteristic strain noise amplitude for a single detector | 0

## References

  * (1) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440-1443 (1977) doi:10.1103/PhysRevLett.38.1440
  * (2) R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791-1797 (1977) doi:10.1103/PhysRevD.16.1791
  * (3) S. Weinberg, Phys. Rev. Lett. 40, 223-226 (1978) doi:10.1103/PhysRevLett.40.223
  * (4) F. Wilczek, Phys. Rev. Lett. 40, 279-282 (1978) doi:10.1103/PhysRevLett.40.279
  * (5) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979) doi:10.1103/PhysRevLett.43.103
  * (6) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 166, 493-506 (1980) doi:10.1016/0550-3213(80)90209-6
  * (7) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980)
  * (8) M. Dine, W. Fischler and M. Srednicki, Phys. Lett. B 104, 199-202 (1981) doi:10.1016/0370-2693(81)90590-6
  * (9) P. Sikivie, Rev. Mod. Phys. 93, no.1, 015004 (2021) doi:10.1103/RevModPhys.93.015004 [arXiv:2003.02206 [hep-ph]].
  * (10) J. Preskill, M. B. Wise and F. Wilczek, Phys. Lett. B 120, 127-132 (1983) doi:10.1016/0370-2693(83)90637-8
  * (11) M. Dine and W. Fischler, Phys. Lett. B 120, 137-141 (1983) doi:10.1016/0370-2693(83)90639-1
  * (12) P. Sikivie, Lect. Notes Phys. 741, 19-50 (2008) doi:10.1007/978-3-540-73518-2_2 [arXiv:astro-ph/0610440 [astro-ph]].
  * (13) D. Croon, R. Houtz and V. Sanz, JHEP 07, 146 (2019) doi:10.1007/JHEP07(2019)146 [arXiv:1904.10967 [hep-ph]].
  * (14) L. Delle Rose, G. Panico, M. Redi and A. Tesi, JHEP 04, 025 (2020) doi:10.1007/JHEP04(2020)025 [arXiv:1912.06139 [hep-ph]].
  * (15) P. S. B. Dev, F. Ferrer, Y. Zhang and Y. Zhang, JCAP 11, 006 (2019) doi:10.1088/1475-7516/2019/11/006 [arXiv:1905.00891 [hep-ph]].
  * (16) A. Ghoshal and A. Salvio, JHEP 12, 049 (2020) doi:10.1007/JHEP12(2020)049 [arXiv:2007.00005 [hep-ph]].
  * (17) B. Von Harling, A. Pomarol, O. Pujolàs and F. Rompineve, JHEP 04, 195 (2020) doi:10.1007/JHEP04(2020)195 [arXiv:1912.07587 [hep-ph]].
  * (18) C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau, P. Schwaller and G. Servant, et al. JCAP 04, 001 (2016) doi:10.1088/1475-7516/2016/04/001 [arXiv:1512.06239 [astro-ph.CO]].
  * (19) S. J. Huber and T. Konstandin, JCAP 09, 022 (2008) doi:10.1088/1475-7516/2008/09/022 [arXiv:0806.1828 [hep-ph]].
  * (20) A. Kosowsky and M. S. Turner, Phys. Rev. D 47, 4372-4391 (1993) doi:10.1103/PhysRevD.47.4372 [arXiv:astro-ph/9211004 [astro-ph]].
  * (21) A. Kosowsky, M. S. Turner and R. Watkins, Phys. Rev. Lett. 69, 2026-2029 (1992) doi:10.1103/PhysRevLett.69.2026
  * (22) M. Lewicki and V. Vaskonen, Eur. Phys. J. C 81, no.5, 437 (2021) [erratum: Eur. Phys. J. C 81, no.12, 1077 (2021)] doi:10.1140/epjc/s10052-021-09232-3 [arXiv:2012.07826 [astro-ph.CO]].
  * (23) M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, Phys. Rev. D 92, no.12, 123009 (2015) doi:10.1103/PhysRevD.92.123009 [arXiv:1504.03291 [astro-ph.CO]].
  * (24) M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, Phys. Rev. Lett. 112, 041301 (2014) doi:10.1103/PhysRevLett.112.041301 [arXiv:1304.2433 [hep-ph]].
  * (25) C. Caprini, R. Durrer and G. Servant, JCAP 12, 024 (2009) doi:10.1088/1475-7516/2009/12/024 [arXiv:0909.0622 [astro-ph.CO]].
  * (26) A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili and A. Kosowsky, Phys. Rev. D 102, no.8, 083512 (2020) doi:10.1103/PhysRevD.102.083512 [arXiv:1903.08585 [astro-ph.CO]].
  * (27) B. P. Abbott et al. [LIGO Scientific], Class. Quant. Grav. 34, no.4, 044001 (2017) doi:10.1088/1361-6382/aa51f4 [arXiv:1607.08697 [astro-ph.IM]].
  * (28) D. Reitze, R. X. Adhikari, S. Ballmer, B. Barish, L. Barsotti, G. Billingsley, D. A. Brown, Y. Chen, D. Coyne and R. Eisenstein, et al. Bull. Am. Astron. Soc. 51, no.7, 035 (2019) [arXiv:1907.04833 [astro-ph.IM]].
  * (29) E. D. Hall, Galaxies 10, no.4, 90 (2022) doi:10.3390/galaxies10040090
  * (30) E. D. Hall, K. Kuns, J. R. Smith, Y. Bai, C. Wipf, S. Biscans, R. X. Adhikari, K. Arai, S. Ballmer and L. Barsotti, et al. Phys. Rev. D 103, no.12, 122004 (2021) doi:10.1103/PhysRevD.103.122004 [arXiv:2012.03608 [gr-qc]].
  * (31) B. Sathyaprakash, M. Abernathy, F. Acernese, P. Ajith, B. Allen, P. Amaro-Seoane, N. Andersson, S. Aoudia, K. Arun and P. Astone, et al. Class. Quant. Grav. 29, 124013 (2012) [erratum: Class. Quant. Grav. 30, 079501 (2013)] doi:10.1088/0264-9381/29/12/124013 [arXiv:1206.0331 [gr-qc]].
  * (32) M. Maggiore, C. Van Den Broeck, N. Bartolo, E. Belgacem, D. Bertacca, M. A. Bizouard, M. Branchesi, S. Clesse, S. Foffa and J. García-Bellido, et al. JCAP 03, 050 (2020) doi:10.1088/1475-7516/2020/03/050 [arXiv:1912.02622 [astro-ph.CO]].
  * (33) R. A. Fisher, Phil. Trans. Roy. Soc. Lond. A 222, 309-368 (1922) doi:10.1098/rsta.1922.0009
  * (34) M. Pospelov and A. Ritz, Phys. Rev. Lett. 83, 2526-2529 (1999) doi:10.1103/PhysRevLett.83.2526 [arXiv:hep-ph/9904483 [hep-ph]].
  * (35) C. A. Baker, D. D. Doyle, P. Geltenbort, K. Green, M. G. D. van der Grinten, P. G. Harris, P. Iaydjiev, S. N. Ivanov, D. J. R. May and J. M. Pendlebury, et al. Phys. Rev. Lett. 97, 131801 (2006) doi:10.1103/PhysRevLett.97.131801 [arXiv:hep-ex/0602020 [hep-ex]].
  * (36) C. Abel, S. Afach, N. J. Ayres, C. A. Baker, G. Ban, G. Bison, K. Bodek, V. Bondar, M. Burghoff and E. Chanel, et al. Phys. Rev. Lett. 124, no.8, 081803 (2020) doi:10.1103/PhysRevLett.124.081803 [arXiv:2001.11966 [hep-ex]].
  * (37) Y. Asano, E. Kikutani, S. Kurokawa, T. Miyachi, M. Miyajima, Y. Nagashima, T. Shinkawa, S. Sugimoto and Y. Yoshimura, Phys. Lett. B 107, 159 (1981) doi:10.1016/0370-2693(81)91172-2
  * (38) K. Mimasu and V. Sanz, JHEP 06, 173 (2015) doi:10.1007/JHEP06(2015)173 [arXiv:1409.4792 [hep-ph]].
  * (39) J. E. Kim, Phys. Rept. 150, 1-177 (1987) doi:10.1016/0370-1573(87)90017-2
  * (40) P. Sikivie, Phys. Rev. Lett. 51, 1415-1417 (1983) [erratum: Phys. Rev. Lett. 52, 695 (1984)] doi:10.1103/PhysRevLett.51.1415
  * (41) D. M. Lazarus, G. C. Smith, R. Cameron, A. C. Melissinos, G. Ruoso, Y. K. Semertzidis and F. A. Nezrick, Phys. Rev. Lett. 69, 2333-2336 (1992) doi:10.1103/PhysRevLett.69.2333
  * (42) S. Moriyama, M. Minowa, T. Namba, Y. Inoue, Y. Takasu and A. Yamamoto, Phys. Lett. B 434, 147 (1998) doi:10.1016/S0370-2693(98)00766-7 [arXiv:hep-ex/9805026 [hep-ex]].
  * (43) S. J. Asztalos et al. [ADMX], Phys. Rev. Lett. 104, 041301 (2010) doi:10.1103/PhysRevLett.104.041301 [arXiv:0910.5914 [astro-ph.CO]].
  * (44) M. S. Pshirkov and S. B. Popov, J. Exp. Theor. Phys. 108, 384-388 (2009) doi:10.1134/S1063776109030030 [arXiv:0711.1264 [astro-ph]].
  * (45) F. P. Huang, K. Kadota, T. Sekiguchi and H. Tashiro, Phys. Rev. D 97, no.12, 123001 (2018) doi:10.1103/PhysRevD.97.123001 [arXiv:1803.08230 [hep-ph]].
  * (46) A. Hook, Y. Kahn, B. R. Safdi and Z. Sun, Phys. Rev. Lett. 121, no.24, 241102 (2018) doi:10.1103/PhysRevLett.121.241102 [arXiv:1804.03145 [hep-ph]].
  * (47) C. W. Chiang and B. Q. Lu, JCAP 05, 049 (2021) doi:10.1088/1475-7516/2021/05/049 [arXiv:2012.14071 [hep-ph]].
  * (48) M. Ahmadvand, JHEP 10, 109 (2021) doi:10.1007/JHEP10(2021)109 [arXiv:2108.00958 [hep-ph]].
  * (49) J. Yang, N. Xie and F. P. Huang, [arXiv:2306.17113 [hep-ph]].
  * (50) J. Yang and F. P. Huang, Phys. Rev. D 108, no.10, 103002 (2023) doi:10.1103/PhysRevD.108.103002 [arXiv:2306.12375 [hep-ph]].
  * (51) N. Xie and F. P. Huang, Sci. China Phys. Mech. Astron. 67, no.1, 210411 (2024) doi:10.1007/s11433-023-2172-7 [arXiv:2207.11145 [hep-ph]].
  * (52) P. F. Michelson, “On detecting stochastic background gravitational radiation with terrestrial detectors,” Monthly Notices of the Royal Astronomical Society, vol. 227, no. 4, pp. 933–941, 1987.
  * (53) N. Christensen, “Measuring the stochastic gravitational-radiation background with laser-interferometric antennas,” Physical Review D, vol. 46, no. 12, p. 5250, 1992.
  * (54) B. Allen and J. D. Romano, Phys. Rev. D 59, 102001 (1999) doi:10.1103/PhysRevD.59.102001 [arXiv:gr-qc/9710117 [gr-qc]].
  * (55) M. Maggiore, Oxford University Press, 2007, ISBN 978-0-19-171766-6, 978-0-19-852074-0 doi:10.1093/acprof:oso/9780198570745.001.0001
  * (56) E. Thrane and J. D. Romano, Phys. Rev. D 88, no.12, 124032 (2013) doi:10.1103/PhysRevD.88.124032 [arXiv:1310.5300 [astro-ph.IM]].
  * (57) S. R. Coleman and E. J. Weinberg, Phys. Rev. D 7, 1888-1910 (1973) doi:10.1103/PhysRevD.7.1888
  * (58) L. Dolan and R. Jackiw, Phys. Rev. D 9, 3320-3341 (1974) doi:10.1103/PhysRevD.9.3320
  * (59) M. E. Carrington, Phys. Rev. D 45, 2933-2944 (1992) doi:10.1103/PhysRevD.45.2933
  * (60) P. B. Arnold and O. Espinosa, Phys. Rev. D 47, 3546 (1993) [erratum: Phys. Rev. D 50, 6662 (1994)] doi:10.1103/PhysRevD.47.3546 [arXiv:hep-ph/9212235 [hep-ph]].
  * (61) R. R. Parwani, Phys. Rev. D 45, 4695 (1992) [erratum: Phys. Rev. D 48, 5965 (1993)] doi:10.1103/PhysRevD.45.4695 [arXiv:hep-ph/9204216 [hep-ph]].
  * (62) P. Basler and M. Mühlleitner, Comput. Phys. Commun. 237, 62-85 (2019) doi:10.1016/j.cpc.2018.11.006 [arXiv:1803.02846 [hep-ph]].
  * (63) J. R. Espinosa, T. Konstandin, J. M. No and G. Servant, JCAP 06, 028 (2010) doi:10.1088/1475-7516/2010/06/028 [arXiv:1004.4187 [hep-ph]].
  * (64) X. Wang, C. Tian and F. P. Huang, JCAP 07, 006 (2023) doi:10.1088/1475-7516/2023/07/006 [arXiv:2301.12328 [hep-ph]].
  * (65) X. Wang, F. P. Huang and X. Zhang, [arXiv:2011.12903 [hep-ph]].
  * (66) S. Jiang, F. P. Huang and X. Wang, Phys. Rev. D 107, no.9, 095005 (2023) doi:10.1103/PhysRevD.107.095005 [arXiv:2211.13142 [hep-ph]].
  * (67) A. Megevand and A. D. Sanchez, Nucl. Phys. B 825, 151-176 (2010) doi:10.1016/j.nuclphysb.2009.09.019 [arXiv:0908.3663 [hep-ph]].
  * (68) A. D. Linde, Phys. Lett. B 100, 37-40 (1981) doi:10.1016/0370-2693(81)90281-1
  * (69) A. D. Linde, “Decay of the false vacuum at finite temperature,” Nuclear Physics B, vol. 216, no. 2, pp. 421–445, 1983.
  * (70) S. R. Coleman, Phys. Rev. D 15, 2929-2936 (1977) [erratum: Phys. Rev. D 16, 1248 (1977)] doi:10.1103/PhysRevD.16.1248
  * (71) C. G. Callan, Jr. and S. R. Coleman, Phys. Rev. D 16, 1762-1768 (1977) doi:10.1103/PhysRevD.16.1762
  * (72) X. Wang, F. P. Huang and X. Zhang, JCAP 05, 045 (2020) doi:10.1088/1475-7516/2020/05/045 [arXiv:2003.08892 [hep-ph]].
  * (73) C. L. Wainwright, Comput. Phys. Commun. 183, 2006-2013 (2012) doi:10.1016/j.cpc.2012.04.004 [arXiv:1109.4189 [hep-ph]].
  * (74) M. Gorghetto and G. Villadoro, JHEP 03, 033 (2019) doi:10.1007/JHEP03(2019)033 [arXiv:1812.01008 [hep-ph]].
  * (75) G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, JHEP 01, 034 (2016) doi:10.1007/JHEP01(2016)034 [arXiv:1511.02867 [hep-ph]].
  * (76) C. O’Hare, cajohare/AxionLimits: AxionLimits, https://cajohare.github.io/AxionLimits/, [Online; accessed 16-April-2023].
  * (77) S. Hamimeche and A. Lewis, Phys. Rev. D 77, 103013 (2008) doi:10.1103/PhysRevD.77.103013 [arXiv:0801.0554 [astro-ph]].
  * (78) C. Gowling and M. Hindmarsh, JCAP 10, 039 (2021) doi:10.1088/1475-7516/2021/10/039 [arXiv:2106.05984 [astro-ph.CO]].
  * (79) C. Caves, Phys. Rev. Lett. 45, no.2, 75-79 (1980) doi:10.1103/PhysRevLett.45.75
  * (80) B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 120, no.9, 091101 (2018) doi:10.1103/PhysRevLett.120.091101 [arXiv:1710.05837 [gr-qc]].
  * (81) P. A. Rosado, Phys. Rev. D 84, 084004 (2011) doi:10.1103/PhysRevD.84.084004 [arXiv:1106.5795 [gr-qc]].
  * (82) S. Sachdev, T. Regimbau and B. S. Sathyaprakash, Phys. Rev. D 102, no.2, 024051 (2020) doi:10.1103/PhysRevD.102.024051 [arXiv:2002.05365 [gr-qc]].