# Gravitational wave echo of relaxion trapping

Abhishek Banerjee<EMAIL_ADDRESS>Department of Particle
Physics and Astrophysics, Weizmann Institute of Science, Rehovot 7610001,
Israel  Eric Madge<EMAIL_ADDRESS>Department of Particle
Physics and Astrophysics, Weizmann Institute of Science, Rehovot 7610001,
Israel  Gilad Perez<EMAIL_ADDRESS>Department of Particle Physics
and Astrophysics, Weizmann Institute of Science, Rehovot 7610001, Israel
Wolfram Ratzinger<EMAIL_ADDRESS>PRISMA+ Cluster of Excellence and
Mainz Institute for Theoretical Physics, Johannes Gutenberg-Universität Mainz,
55099 Mainz, Germany  Pedro Schwaller<EMAIL_ADDRESS>PRISMA+
Cluster of Excellence and Mainz Institute for Theoretical Physics, Johannes
Gutenberg-Universität Mainz, 55099 Mainz, Germany

###### Abstract

To solve the hierarchy problem, the relaxion must remain trapped in the
correct minimum, even if the electroweak symmetry is restored after reheating.
In this scenario, the relaxion starts rolling again until the backreaction
potential, with its set of local minima, reappears. Depending on the time of
barrier reappearance, Hubble friction alone may be insufficient to retrap the
relaxion in a large portion of the parameter space. Thus, an additional source
of friction is required, which might be provided by coupling to a dark photon.
The dark photon experiences a tachyonic instability as the relaxion rolls,
which slows down the relaxion by backreacting to its motion, and efficiently
creates anisotropies in the dark photon energy-momentum tensor, sourcing
gravitational waves. We calculate the spectrum of the resulting gravitational
wave background from this new mechanism, and evaluate its observability by
current and future experiments. We further investigate the possibility that
the coherently oscillating relaxion constitutes dark matter and present the
corresponding constraints from gravitational waves.

††preprint:

##  I Introduction

Figure 1: Available parameter space for $\xi=100$ (top), and for $\xi=10$
(bottom). The red and orange shaded regions are excluded by the indicated
constraints or combinations thereof. Above the red solid line, the relaxion
decay constant becomes super-Planckian. The grey dashed line encloses the
parameter space in which relaxation can be realized without dark photon
friction, which is discussed in more details in Appendix A. The prospective GW
sensitivity of $\mu$Ares (green) as well as SKA after an observation period of
5 years (turquoise) and 20 years (blue) is indicated by the respective
coloured regions. In the purple coloured region, a subrange of the viable
reappearance temperatures can be excluded using current NANOGrav data from the
11 year data set. The regions bounded by the coloured dotted lines need super-
Planckian decay constants to be probed by the respective experiment. In the
lower panel, the grey shaded region can reproduce our best-fit spectrum (at
$T_{\text{ra}}\sim$20\text{\,}\mathrm{MeV}$$) to the potential stochastic GW
background seen in the recent NANOGrav data. An animated version of the lower
panel of the figure, explicitly showing the dependence on the reappearance
temperature, can be found in the ancillary material of this paper.

The weakness of the gravitational force ensures that gravitational waves (GWs)
can store information about the history of our Universe. On the other hand,
though, it also implies that only dramatic and cosmological events may leave
imprints that can be observed at present or in the near future. GW detectors
are, for instance, sensitive to primordial phase transitions, inflation, and
oscillatory motion of ultralight fields (see, e.g., Refs. [1, 2] for a recent
discussion of potential signals of cosmological origin).

In this work, we introduce a different mechanism that leads to the production
of GWs in the early Universe. In the following, we demonstrate the idea using
the relaxion framework [3]; however, the core scheme may be applicable to
other forms of new physics. At early times, the relaxion is trapped in a local
minimum. However, assuming that the reheating temperature is above the
electroweak (EW) scale,111The requirement of a high reheating temperature,
motivated by a large class of inflation models (see, e.g., Ref. [4] and
references therein), is generic and also needed in most models that explain
the observed baryon abundance (see ,e.g., Ref. [5] and references therein).
we expect the EW symmetry to be restored. Hence, the minimum in which the
relaxion was originally trapped disappears, and the relaxion starts rolling
down its potential. At some temperature below the EW scale, we expect the set
of local minima to reappear. There is a limited region of the model’s
parameter space such that the relaxion would be trapped still within the same
minimum [6, 7]. In this case, even the minimal realization of the relaxation
mechanism could lead to a viable ultra light dark matter (DM) candidate [6],
which may be tested in the future due to the presence of relaxion-Higgs mixing
[8, 9].

Despite the attractiveness of this minimal setup, there is a sizeable region
of parameter space where the post-reheating displacement of the relaxion would
be too large and the relaxion would not be trapped but instead would descend
in an uncontrolled way toward the global minimum. This runaway can be avoided
if an additional source of friction is added to the model. Coupling a scalar
field to a dark $U(1)$ gauge field strength is known to lead to an efficient
energy loss mechanism in the form of tachyonic dark photon field production
[10, 11] (see Refs. [12, 13] for the case of a relaxion or other axionlike
particles (ALPs) in the mass range $\gtrsim\mathcal{O}($\mathrm{MeV}$)$, where
self-friction may be induced, or Ref. [14] for stopping via an instability in
a modified relaxion potential, albeit outside the region of interest of this
work). Dark photon production quickly reaches a quasi-steady-state where the
friction balances the slope of the potential. At each time, the dominantly
produced momentum-mode $k=\xi aH$ is about to exit the tachyonic band, where
$\xi$ is an $\mathcal{O}(10-100)$ parameter, $H$ is the Hubble rate, and $a$
is the scale factor.

The energy density stored in this mode is roughly constant, $\rho_{X}\sim
m_{\phi}^{2}f_{\phi}^{2}$ [6], where $m_{\phi}$ and $f_{\phi}$ are the
relaxion mass and the decay constant, respectively, and is the source of GW
production. The rolling of the relaxion and the dark photon production stop
around the time when the potential barriers reappear, which traps the relaxion
and ends this epoch.222For simplicity, we are agnostic about the details of
the reappearance process of the backreaction potential and assume
instantaneous reappearance. Various considerations regarding barrier
reappearance can also generate gravitational waves, a scenario considered in
Ref. [15].

Given the above mechanism, we can estimate the GW signal as follows. First,
the peak frequency at present time is obtained by redshifting the dominant
$k$-mode at the time of reappearance, cf. Eq. 35,

$f_{\text{peak}}\sim\frac{a_{\text{ra}}}{a_{0}}\xi
H_{\text{ra}}\sim$1\text{\,}\mathrm{\SIUnitSymbolMicro
Hz}$\left(\frac{\xi}{25}\right)\\!\left(\frac{T_{\text{ra}}}{$1\text{\,}\mathrm{GeV}$}\right),$
(1)

where $T_{\text{ra}}$ is the temperature when the potential barriers reappear.
Here and in the following, quantities indexed “0” and “ra” are evaluated today
and at the time when the potential barriers reappear, respectively. Second,
the GW amplitude is roughly given by the square of the energy density stored
in the dark photon just before $T_{\text{ra}}$, cf. Eq. 28,

$\displaystyle\\!\\!\Omega_{\text{GW}}^{\text{peak}}$
$\displaystyle\sim\frac{1}{\rho_{c}^{0}}\,\frac{(c_{\text{eff}}\,\rho_{X}^{\text{ra}}/f_{\text{peak}}^{\text{ra}})^{2}}{M_{\text{Pl}}^{2}}\left(\frac{a_{\text{ra}}}{a_{0}}\right)^{4}$
$\displaystyle\sim${10}^{-10}$\left(\frac{25}{\xi}\right)^{\\!2}\\!\left(\frac{m_{\phi}}{$0.1\text{\,}\mathrm{eV}$}\frac{f_{\phi}}{${10}^{10}\text{\,}\mathrm{GeV}$}\right)^{\\!4}\\!\left(\frac{{$1\text{\,}\mathrm{GeV}$}}{T_{\text{ra}}}\right)^{\\!8}\\!$
(2)
$\displaystyle\sim${10}^{-12}$\left(\frac{25}{\xi}\right)^{\\!2}\\!\left(\frac{$0.1\text{\,}\mathrm{eV}$}{m_{\phi}}\frac{\sin\theta_{\phi
h}}{${10}^{-13}$}\right)^{\\!\\!12}\\!\left(\frac{$\mathrm{TeV}$}{\Lambda}\frac{$\mathrm{GeV}$}{T_{\text{ra}}}\right)^{\\!8}\\!\\!,\\!\\!$

where $\rho_{c}^{0}=3M_{\text{Pl}}^{2}H_{0}^{2}$ is the critical energy
density of the Universe today, $f_{\text{peak}}^{\text{ra}}=\xi H_{\text{ra}}$
is the peak frequency without red-shifting, and $c_{\text{eff}}^{2}\simeq 1/4$
is the efficiency factor for converting dark photon energy into GWs obtained
from the analytic calculation in Appendix C. Note that
$c_{\text{eff}}\sim\mathcal{O}\left(1\right)$ as the relaxion potential energy
at each time is predominantly deposited in a narrow range of exponentially
growing dark photon momentum modes, resulting in large time-dependent
inhomogeneities and an efficient generation of anisotropic stress sourcing
GWs, in a similar manner as in the case of tachyonic preheating [16, 17, 18,
19]. In the last line, we have reexpressed the signal strength in terms of the
relaxion-Higgs mixing $\sin\theta_{h\phi}$ and the cutoff scale $\Lambda$.
This result shows that some of the parameter space of the model may lead to a
visible signal in near-future GW experiments, allowing us to probe parameter
regions that are currently unexplored by other experiments, as discussed
below. In addition, we note that the relation between the physical parameters
of the models and the GW amplitude is given by $\sin\theta_{h\phi}\propto
m_{\phi}\times\left(\Omega_{\text{GW}}^{\text{peak}}\right)^{1/12}\,,$ showing
a rather mild dependence on the actual amplitude.

For convenience, the range of relaxion masses, $m_{\phi}$, and mixing angles
with the Higgs boson, $\sin\theta_{h\phi}$, which can be probed via current or
future GW experiments as well as the corresponding constraints on the
parameter space are summarized in Fig. 1. The green and blue/turquoise
coloured regions can be accessed with $\mu$Ares and the Square-Kilometre Array
(SKA) observatory, respectively, depending on the temperature of barrier
reappearance. In the purple region, the reappearance temperature is restricted
by current data from the North American Nanohertz Observatory for
Gravitational Waves (NANOGrav). In addition to that, in the grey shaded
region, we present the parameter range in which our model can account for the
potential GW signal recently observed in NANOGrav data. The grey dashed line
encloses the region in which the relaxion can be trapped without dark photon
friction. On the other hand, as the figure illustrates, there is a large
fraction of parameter space where an additional source of friction is required
for the viability of the relaxion mechanism and thus motivates us to add a
relaxion-dark photon coupling to the model. A more detailed discussion of the
figure is deferred to Sec. IV.

This paper is organized as follows. We briefly review the relaxion and the
dark photon dynamics in Sec. II. Sections II.1 and II.2 contain a brief
discussion of the interplay of the relaxion-dark photon dynamics and the
possibility of the relaxion being DM, respectively. In Sec. II.3 we review the
constraints on our model and discuss the available parameter space.
Subsequently, the production mechanism of the GW background is studied in Sec.
III. We derive the GW spectrum in Sec. III.1 and briefly describe how the
detectability of the signal is evaluated in Sec. III.2. The results of this
paper are then discussed in Sec. IV. Section V concludes the paper. A brief
discussion of the minimal relaxion scenario (without a dark photon coupling)
is deferred to Appendix A. Further details regarding the calculation of the
dark photon and GW spectra are provided in Appendices B and C, respectively.

##  II Setup

In this work, we consider the relaxion $\phi$ coupled to a dark photon field
$X_{\mu}$,

$-\mathcal{L}\supset
V(H,\phi)+\frac{r_{X}}{4}\,\frac{\phi}{f_{\phi}}\,X_{\mu\nu}\widetilde{X}^{\mu\nu}\,,$
(3)

with the potential of the relaxion field $\phi$ and Higgs doublet $H$ given by

$V(H,\phi)=V_{\text{roll}}(\phi)+\mu_{H}^{2}(\phi)|H|^{2}+\lambda|H|^{4}+V_{\text{br}}(H,\phi)\,,$
(4)

where $\lambda$ is the Higgs quartic coupling and

$\displaystyle V_{\text{roll}}(\phi)$ $\displaystyle=-cg\Lambda^{3}\phi\,,$
(5a) $\displaystyle\mu_{H}^{2}(\phi)$
$\displaystyle=\Lambda^{2}-g\Lambda\phi\,,$ (5b) $\displaystyle
V_{\text{br}}(H,\phi)$
$\displaystyle=-\frac{\Lambda_{\text{br}}^{4}}{v_{H}^{2}}|H|^{2}\cos\frac{\phi}{f_{\phi}}\,.$
(5c)

Here, $c$ is an $\mathcal{O}(1)$ number, $g$ is a dimensionless parameter,
$\Lambda$ is the Higgs mass cutoff scale, $\Lambda_{\text{br}}$ is the
backreaction scale, $v_{H}=\langle|H|\rangle=$174\text{\,}\mathrm{GeV}$$ is
the Higgs vacuum expectation value, and $f_{\phi}$ is the decay constant of
the relaxion.

During inflation, the relaxion rolls down the linear slope of its potential
$V_{\text{roll}}$. It thereby scans the Higgs mass parameter $\mu^{2}(\phi)$.
Once $\mu^{2}$ crosses zero, the Higgs acquires a nonvanishing vacuum
expectation value, triggering the breaking of the EW gauge symmetry. The Higgs
then backreacts creating wiggles in the relaxion potential via
$V_{\text{br}}$. Once the Higgs backreaction balances the rolling potential,
the relaxion is trapped in the first minimum it encounters. Choosing
$cg\Lambda^{3}f_{\phi}\sim\Lambda_{\text{br}}^{4}$, we end up with a weak-
scale expectation value for the Higgs boson, solving the hierarchy problem.
The relaxion mass and the relaxion-Higgs mixing angle can then be written as
[7, 9]

$\displaystyle m_{\phi}^{2}\simeq\frac{\Lambda_{\rm
br}^{6}}{f_{\phi}^{2}\Lambda
v_{H}}\,,\quad\sin\theta_{h\phi}\simeq\sqrt{2}\left(\frac{m_{\phi}^{4}f_{\phi}\Lambda^{2}}{v_{H}\,m_{h}^{6}}\right)^{\frac{1}{3}}\\!,$
(6)

in terms of the theory parameters. Here, $m_{h}=$125\text{\,}\mathrm{GeV}$$ is
the Higgs mass.

###  II.1 Relaxion and dark photon evolution

After reheating, the EW symmetry will be restored due to thermal corrections
to the potential, provided that the reheating temperature is above the EW
phase transition temperature. As a consequence, the relaxion will start
rolling again, leading to exponential production of dark photon modes. To see
the interplay, the coupled differential equations describing the evolution of
the spatially homogeneous relaxion and the dark photon modes are given by

$\displaystyle\theta^{\prime\prime}+2aH\theta^{\prime}+\frac{a^{2}}{f_{\phi}^{2}}\frac{\partial
V_{\rm
roll}}{\partial\theta}=-\frac{a^{2}}{f_{\phi}^{2}}\,\frac{r_{X}}{4\,a^{4}}\left<X_{\mu\nu}\tilde{X}^{\mu\nu}\right>,$
(7) $\displaystyle
X^{\prime\prime}_{\lambda}(k,\tau)+(k^{2}-\lambda\,k\,r_{X}\theta^{\prime})X_{\lambda}(k,\tau)=0\,,$
(8)

where $\theta=\phi/f_{\phi}$ and primes denote derivatives with respect to
conformal time $\tau$ with $a\,d\tau=dt$. We have written the dark photon in
terms of its Fourier modes $X_{\lambda}(\mathbf{k},\tau)$ in Coulomb gauge,
$\nabla\cdot\mathbf{X}=0$, as

$\hat{\mathbf{X}}(\mathbf{x},\tau)=\int\\!\frac{d^{3}k}{(2\pi)^{3}}\sum\limits_{\lambda=\pm}X_{\lambda}(k,\tau)\mathbf{\varepsilon}_{\lambda}(\mathbf{k})\hat{a}_{\lambda}(\mathbf{k})+\text{h.c.}$
(9)

Here, $\lambda$ denotes the dark photon helicity, and
$\mathbf{\varepsilon}_{\lambda}$ are the corresponding circular polarization
vectors.

As evident from Eq. 8, dark photon modes with $0<k<\lambda
r_{X}\theta^{\prime}$ are tachyonic and will experience exponential growth
compared to the vacuum fluctuations. The resulting dark photon spectrum then
features anisotropies in its energy-momentum tensor which will act as a source
for GW production, leading to a stochastic GW background [17, 18, 19, 20, 21,
22]. Furthermore, since only modes of the helicity with the same sign as the
relaxion velocity can become tachyonic, the rolling relaxion will produce a
circular polarized dark photon background. In our case, as we choose
$\theta^{\prime}>0$, only the positive-helicity modes are exponentially
produced. For these modes, the solution to the equations of motion is given in
the Wentzel–Kramers–Brillouin (WKB) approximation by

$X_{+}(k,\tau)=\frac{e^{g(k,\tau)}}{\sqrt{2\,\Omega(k,\tau)}}\,,$ (10)

where $\Omega^{2}(k,\tau)=k\,r_{X}|\theta^{\prime}(\tau)|-k^{2}>0$ is the
corresponding tachyonic frequency and
$g(k,\tau)=\int^{\tau}\\!\\!d\tau^{\prime}\,\Omega(k,\tau^{\prime})$. The
approximation holds for $|\Omega^{\prime}/\Omega^{2}|\ll 1$.

At early times, just after reheating, the friction from dark photons can be
neglected. Assuming radiation domination, the relaxion velocity can then be
written as

$\theta^{\prime}(\tau)=\frac{\Lambda_{\rm
br}^{4}}{5f_{\phi}^{2}}\left(\frac{a_{\text{rh}}}{\tau_{\text{rh}}}\right)^{2}\\!\tau^{3}\left[1-\left(\frac{\tau_{\text{rh}}}{\tau}\right)^{5}\right]\,,$
(11)

where we imposed $\theta^{\prime}(\tau_{\text{rh}})=0$. Subsequently, due to
the exponential production, the dark photon friction becomes comparable to the
other terms in the equation of motion, Eq. 7, and asymptotes to the slope of
the rolling potential. The timescale $\tau_{\text{pp}}$ at which particle
production kicks in can be determined by

$\frac{\langle X_{\mu\nu}\tilde{X}^{\mu\nu}\rangle(\tau_{\rm
pp})}{4a^{4}(\tau_{\rm
pp})}\approx\frac{\tilde{k}^{4}e^{2g(\tilde{k},\tilde{\tau})}}{4\pi^{2}a^{4}}\sim\frac{\Lambda_{\text{br}}^{4}}{r_{X}}\,,$
(12)

where $\tilde{k}$ is the mode that dominates the $\langle X\tilde{X}\rangle$
term, given by the saddle point $\left.\partial g(k,\tau)/\partial
k\right|_{\tilde{k}}=0$. After $\tau_{\text{pp}}$, due to the balance between
the potential slope and the backreaction from the dark photon, the relaxion
field velocity becomes proportional to the Hubble rate and evolves as [10, 6]

$\theta^{\prime}(\tau)\approx\frac{\xi}{r_{X}}\,a(\tau)H(\tau)\left(1+\epsilon\log\frac{\tau}{\tau_{\text{pp}}}\right)\approx\frac{\xi}{r_{X}\tau}\,,$
(13)

with a small logarithmic correction ($\epsilon\\!\ll\\!1$). We defined here
the parameter $\xi=\frac{r_{X}|\theta^{\prime}|}{aH}$ at $\tau_{\text{pp}}$.
From Eq. 12 we then obtain $\xi\sim\mathcal{O}(10-100)$ with a mere
logarithmic dependence on the relaxion parameters [10]. The dominating
$k$-mode at each epoch then becomes $\tilde{k}/a\sim
r_{X}\dot{\theta}(t)\sim\xi H$.

Once the Universe has cooled sufficiently, the EW phase transition occurs, and
the wiggles of the backreaction potential reappear. The rolling of the
relaxion between reheating and $T_{\text{ra}}$ leads to a displacement from
the minimum in which it originally settled during inflation by

$\Delta\theta=\int\limits_{\tau_{\text{rh}}}^{\tau_{\text{ra}}}\\!d\tau\,\theta^{\prime}\approx\frac{\xi}{4\,r_{X}}\left[1+\log\frac{H_{\text{pp}}^{2}}{H_{\text{ra}}^{2}}\right].$
(14)

For the relaxion to remain trapped in this minimum, we need to require that
the displacement is less than the distance between the minimum and the next
maximum, $\Delta\theta\lesssim 2\delta$, where
$\delta=\Lambda_{\text{br}}^{2}/(\Lambda v_{H})$ [6, 7]. This sets a lower
bound on the coupling to the dark photons. For smaller couplings, the dark
photon friction is insufficient to prevent the relaxion from rolling into one
of the neighbouring minima. The relaxion then needs to traverse
$\Delta\theta\sim\mathcal{O}(n)$ to end up in the $n$-th minimum, where $n=1$
denotes the minimum in which it stopped during inflation, extending the
parameter space of the theory. However, going beyond the first minimum
requires a careful adjustment of the initial conditions to let the relaxion
stop exactly at the bottom of the $n$-th minimum at reappearance. Otherwise,
the time required for the relaxion to reach the bottom would exceed the age of
the Universe. We thus simply assume $r_{X}=\xi/(2\delta)$ in the following.333
As $\delta=\frac{\Lambda_{\text{br}}^{2}}{\Lambda
v_{H}}<\frac{v_{H}}{\Lambda}$, this implies $r_{X}\gtrsim
10^{2}\left(\frac{\Lambda}{$1\text{\,}\mathrm{TeV}$}\right)\left(\frac{\xi}{$10$}\right)$.
Such large couplings can be obtained in a technically natural way for example
via the clockwork mechanism [23, 24, 25, 26].

###  II.2 Relaxion dark matter

After the reappearance of the Higgs backreaction potential, the displaced
relaxion begins to oscillate around the minimum of its potential, providing a
candidate for ultralight DM as discussed in Ref. [6]. Assuming the maximal
displacement of $\Delta\theta=2\delta$, the relaxion relic abundance is given
by

$\Omega_{\phi}=\frac{2\,m_{\phi}^{2}f_{\phi}^{2}\delta^{2}}{3M_{\text{Pl}}^{2}H_{0}^{2}}\frac{g_{*s}(T_{0})T_{0}^{3}}{g_{*s}(T_{\text{osc}})T_{\text{osc}}^{3}}\,,$
(15)

where $T_{\text{osc}}\sim\min[T_{\text{ra}},\sqrt{m_{\phi}M_{\text{Pl}}}]$ is
the temperature when the relaxion starts to oscillate. Requiring that the
relaxion provides all of DM, i.e., $\Omega_{\phi}h^{2}=0.12$ [27], the
relaxion decay constant can be expressed as

$\displaystyle\\!\\!\\!f_{\phi}\sim$5\text{\times}{10}^{9}\text{\,}\mathrm{GeV}$\left(\frac{\Lambda}{$1\text{\,}\mathrm{TeV}$}\right)^{\\!\frac{2}{5}}\\!\left(\frac{T_{\text{osc}}}{$1\text{\,}\mathrm{GeV}$}\right)^{\\!\frac{9}{10}}\\!\left(\frac{$0.1\text{\,}\mathrm{eV}$}{m_{\phi}}\right)\\!.$
(16)

As we are considering coherently oscillating relaxion DM here, its mass needs
to be less than approximately $10\text{\,}\mathrm{eV}$ in order to be
described by a classical field [28]. For this range of relaxion masses, the
possible decay channels are into two photons and two dark photons,
$\Gamma_{\phi}=\Gamma_{\gamma\gamma}+\Gamma_{XX}$. The decay rate into a dark
photon pair is given by

$\Gamma_{XX}=\frac{r_{X}^{2}}{64\pi}\frac{m_{\phi}^{3}}{f_{\phi}^{2}}\,,$ (17)

while the decay width of $\phi\to\gamma\gamma$ is subdominant compared to that
of $\phi\to XX$ as it is suppressed by the square of the relaxion-Higgs mixing
angle [6], which in turn is bounded from above by $\sin\theta_{h\phi}\lesssim
v_{H}/f_{\phi}$ [9]. The relaxion lifetime hence becomes

$\tau_{\phi}\sim$20\text{\,}\mathrm{Gyr}$\left(\frac{25}{\xi}\right)^{2}\left(\frac{T_{\text{osc}}}{$1\text{\,}\mathrm{GeV}$}\right)^{3}\left(\frac{$0.1\text{\,}\mathrm{eV}$}{m_{\phi}}\right)^{5},$
(18)

where we have chosen $r_{X}=\xi/(2\delta)$. Since the decay of DM into
relativistic particles affects the spectrum of the cosmic microwave background
(CMB) at low-$\ell$ multipoles, the lifetime is constrained as
$\tau_{\phi}>$160\text{\,}\mathrm{Gyr}$$ [29].

As shown in Refs. [30, 31, 32, 22], for $r_{X}(\Delta\theta)_{\rm
sep}\sim\mathcal{O}(10^{2})$, an oscillating ALP may introduce a second phase
of tachyonic dark photon production, which could suppress the DM abundance by
up to two orders of magnitude. This condition is satisfied in parts of the
parameter space where the field displacement is maximal, since there
$r_{X}\sim\mathcal{O}(\xi/\delta)$ and $\xi\sim\mathcal{O}(10-10^{2})$, and
introduces some uncertainty in our estimate of the DM abundance in those
regions.

###  II.3 Constraints

A successful cosmological relaxation of the Higgs mass requires the inflation
sector to dominate the total energy density,
$3H_{\text{I}}^{2}M_{\text{Pl}}^{2}\gtrsim\Lambda^{4}$, as well as that the
classical motion of the relaxion dominates over quantum fluctuation during
inflation, $(\Delta\phi)_{\rm cl}\gtrsim H_{\rm I}/2\pi$. Here, $H_{\text{I}}$
is the Hubble scale during inflation. Combining these two constraints, we get
an upper bound on the cutoff scale $\Lambda$,

$\displaystyle\Lambda\lesssim\left(\frac{2\pi\sqrt{3}\,M_{\text{Pl}}^{3}\Lambda_{\text{br}}^{4}}{f}\right)^{1/6}\\!\\!.$
(19)

As we are considering a Higgs-dependent backreaction potential, we also
require $\Lambda_{\text{br}}\lesssim v_{H}$ [3, 9].444See Refs. [11, 13] for a
discussion of other stopping mechanisms in which this constraint can be
relaxed. The allowed range of the effective cutoff $\Lambda$ of the theory is

$\displaystyle M_{\text{Pl}}\gtrsim f_{\phi}\gtrsim\Lambda\gtrsim 4\pi
v_{H}\sim$1\text{\,}\mathrm{TeV}$\,.$ (20)

Also, for a generic backreaction potential which does not change the late-time
cosmology, the range of reappearance temperatures is

$\displaystyle v_{H}\gtrsim T_{\text{ra}}\gtrsim T_{\rm
BBN}\sim$10\text{\,}\mathrm{MeV}$\,.$ (21)

For masses below the $\mathrm{eV}$ scale, the relaxion can further mediate
long-range forces. Experiments looking for such interactions (fifth force
experiments, inverse-square-law, and equivalence-principle tests) constrain
the coupling of the relaxion to ordinary matter [33, 34, 35, 36, 37, 38, 39],
which is induced by the relaxion-Higgs mixing angle given in Eq. 6. In a
similar manner, for masses up to the $\mathrm{keV}$ range, the mixing is
constrained from stellar cooling [40, 41, 42, 43], as it induces relaxion-
mediated contributions to the energy loss in stars. Slightly weaker limits on
the mixing angle can furthermore be obtained from bounds on the solar relaxion
flux as constrained by XENON1T and other dark-matter direct detection
experiments [44].

Additional constraints arise when coupling the relaxion to a dark photon field
(see Appendix A for a discussion of the minimal scenario without this
coupling), with the coupling here chosen to saturate the trapping condition,
i.e., $r_{X}=\xi/(2\delta)$. For the dark photon induced friction to trap the
relaxion, reappearance has to occur sufficiently late for the dark photon to
be produced, i.e., $H_{\text{pp}}>H_{\text{ra}}$. This sets a lower bound on
the relaxion mass, $m_{\phi}\gtrsim\sqrt{10}\,\delta H_{\text{ra}}$, for which
the dark photon scenario can be applied.

If this condition is satisfied, throughout its evolution from reheating to
reappearance, the relaxion continuously produces dark photons, depositing
energy density into the latter. At the time of trapping, $t_{\rm ra}$, the
dark photon energy-density can be estimated as

$\rho_{X}(t_{\rm ra})=\frac{1}{a^{4}(t_{\rm ra})}\int^{t_{\rm
ra}}\\!\\!dt^{\prime}a^{4}(t^{\prime})\dot{V}=\frac{1}{2}\,m_{\phi}^{2}f_{\phi}^{2}\,.$
(22)

For the Universe to remain radiation dominated, we require
$\rho_{X}(T_{\text{ra}})$ to be smaller than the radiation energy density at
$T_{\text{ra}}$, i.e. $\rho_{X}(T_{\text{ra}})\lesssim
3M_{\text{Pl}}^{2}H_{\text{ra}}^{2}\sim T_{\text{ra}}^{4}$; otherwise, the
relaxion dominates the total energy density of the Universe. After
$T_{\text{ra}}$, dark photon production stops, and the corresponding energy
density redshifts as that of radiation.

The dark photon is further constrained by its contribution to the energy
density of the Universe at low temperatures, typically parametrized in terms
of the effective number of neutrino species $N_{\text{eff}}$. For temperatures
$T\lesssim m_{e}$, after electrons and positrons have annihilated, the dark
photon energy density (or of any additional relativistic species for that
matter) can be written in terms of its contribution $\Delta N_{\text{eff}}$ to
$N_{\text{eff}}$ as

$\rho_{X}(T)=\frac{7\,\pi^{2}}{120}\Delta
N_{\text{eff}}\left(\frac{4}{11}\right)^{\frac{4}{3}}T^{4},$ (23)

where the $4/11$ factor accounts for the difference between the neutrino and
photon temperature. Using Eq. 22 and the current $95\text{\,}\%$ C.L. limit
$N_{\text{eff}}=3.27\pm 0.29$ from CMB data [27] combined with local
measurements of the Hubble rate [45], as well as
$N_{\text{eff}}^{\text{SM}}=3.046$ in the Standard Model, we obtain a lower
limit on the reappearance temperature as a function of the relaxion mass and
decay constant

$\displaystyle T_{\text{ra}}\gtrsim
2.5\,g_{*s,\text{ra}}^{-1/3}\,\sqrt{m_{\phi}f_{\phi}}\,,$ (24)

where $g_{*s,\text{ra}}$ is the number of entropic degrees of freedom at
reappearance.

If we assume that the relaxion accounts for the full DM abundance, then
plugging Eq. 16 to Eq. 24, we get,

$T_{\text{ra}}\gtrsim$450\text{\,}\mathrm{MeV}$\left(\frac{67}{g_{*s,\text{ra}}}\right)^{\frac{1}{3}}\left(\frac{\Lambda}{$10\text{\,}\mathrm{TeV}$}\right)^{\\!\frac{4}{11}}.$
(25)

Here, we also assume that the relaxion starts to oscillate at $T_{\text{ra}}$,
which is only true for a sufficiently heavy relaxion. As we require
$\Lambda\gtrsim 4\pi v_{H}\sim$1\text{\,}\mathrm{TeV}$$, this sets a lower
bound of $T_{\text{ra}}\gtrsim$240\text{\,}\mathrm{MeV}$$ on the reappearance
temperature for the relaxion DM scenario to be realized.555Note that, for this
value of the reappearance temperature, the relaxion starts to oscillate
directly after reappearance as long as
$m_{\phi}\gtrsim$5\text{\times}{10}^{-10}\text{\,}\mathrm{eV}$$. Upon the same
assumptions, $\Lambda_{\rm br}\lesssim v_{H}$ further sets an upper bound on
the reappearance temperature,

$\displaystyle\\!T_{\text{ra}}\lesssim\min\left[v_{H},\
$80\text{\,}\mathrm{GeV}$\,\left(\frac{${10}^{6}\text{\,}\mathrm{GeV}$}{\Lambda}\right)\left(\frac{96}{g_{\rm*s,ra}}\right)^{\\!\\!\frac{1}{3}}\right].$
(26)

Figure 2: Allowed range of reappearance temperatures $T_{\text{ra}}$ as a
function of the cutoff scale $\Lambda$, while fixing $f_{\phi}$ to the value
reproducing the measured DM abundance, Eq. 16, assuming
$T_{\text{osc}}=T_{\text{ra}}$. The blue shaded region is excluded since
$\Lambda_{\text{br}}>v_{H}$, while in the green shaded regions, the CMB bound
on $\Delta N_{\text{eff}}$, Eq. 25, is violated. In the orange shaded region,
the dark photon energy density further dominates the Universe at reappearance.
The dashed lines are contours of constant
$m_{\phi}f_{\phi}/T_{\text{ra}}^{2}$, which sets the amplitude of the GW
spectrum; cf. Eq. 2.

In Fig. 2, we show the minimal and maximal allowed reappearance temperature
for relaxion DM as a function of the cutoff $\Lambda$ of the theory. Combining
$\Lambda_{\text{br}}<v_{H}$ (blue) and the $N_{\text{eff}}$ constraint
(green), we see that the highest $\Lambda$ for which the relaxion can be
realized as coherently oscillating DM is
$\Lambda\lesssim${10}^{7}\text{\,}\mathrm{GeV}$$, which is in accordance with
the constraints for $T_{\text{ra}}\simeq$6\text{\,}\mathrm{GeV}$$. Because of
the rapid change in the radiative degrees of freedom around the time of the
QCD phase transition, the $\Delta N_{\text{eff}}$ limit on the reappearance
temperature saturates at $T_{\text{ra}}\sim T_{\text{QCD}}$ for
$\Lambda\lesssim$2\text{\,}\mathrm{TeV}$$. We also depict the weaker bound
$\rho_{X}<\rho_{\text{rad}}$ from the total energy density in orange.

Figure 3: Available parameter space (black framed region) for relaxion DM in
the relaxion mass $m_{\phi}$ vs. mixing angle $\sin\theta_{h\phi}$ plane. The
red and orange shaded regions are excluded by the indicated constraints of
combinations thereof. The colored regions inside the viable DM space can be
probed via GWs in $\mu$Ares (green) or SKA (blue/turquoise). The light shading
and solid lines indicate points that can be probed for a subrange of
reappearance temperatures, whereas the darker shaded parts enclosed by dotted
lines are accessible for all valid $T_{\text{ra}}$. An animated version of the
plot scanning the reappearance temperature is enclosed in the ancillary
material of this work.

In Fig. 3, we depict the constraints on the relaxion parameters as a function
of the relaxion mass $m_{\phi}$ and the mixing angle $\sin\theta_{h\phi}$,
where we determined $\Lambda$ and $f_{\phi}$ from the DM abundance using Eqs.
16 and 6. The red and orange shaded regions are excluded by the indicated
constraints, where the shape of the exclusions now partially differs from the
corresponding ones in Fig. 1 as we require here
$\Omega_{\phi}=\Omega_{\text{DM}}$. Scanning over all allowed values of
$T_{\text{ra}}$, the full range of masses and mixing angles for which we can
obtain coherently oscillating relaxion DM is indicated by the black-framed
regions in Fig. 3. We obtain two separated islands of viable parameter space,
one at low masses, $${10}^{-11}\text{\,}\mathrm{eV}$\lesssim
m_{\phi}\lesssim${10}^{-6}\text{\,}\mathrm{eV}$$, with mixing angles around
$${10}^{-23}$\lesssim\sin\theta_{h\phi}\lesssim${10}^{-18}$$, and another
island at high masses, $${10}^{-2}\text{\,}\mathrm{eV}$\lesssim
m_{\phi}\lesssim$1\text{\,}\mathrm{eV}$$, with a narrow range of mixing angles
around $${10}^{-14}$\lesssim\sin\theta_{h\phi}\lesssim${10}^{-11}$$. Note
that, in the high-mass island, the coupling to dark photons is required to
trap the relaxion, whereas in most of the low-mass island, relaxion DM can be
realized without dark photon friction [6]. It shall, moreover, be emphasized
that the low-mass DM regions in the minimal and dark photon scenario are
separated in the reappearance temperature, since relaxion stopping via Hubble
friction requires $m_{\phi}\lesssim\sqrt{8}\,\delta H_{\text{ra}}$, whereas
dark photon production only occurs for $m_{\phi}\gtrsim\sqrt{10}\,\delta
H_{\text{ra}}$. See Appendix A for a detailed discussion on the minimal
scenario.

##  III Gravitational waves

Let us now consider the stochastic GW background generated from anisotropies
in the energy-momentum tensor of the dark photon produced during the
postinflationary evolution of the relaxion. In particular, we focus here on
the GWs sourced during the rolling of the relaxion between reheating and the
EW phase transition. Further gravitational radiation can be generated when the
relaxion oscillates around the minimum of its potential, as explored in Refs.
[20, 21, 22, 46, 47, 48, 49] in the context of general ALP models, or in a
confining phase transition that generates the backreaction potential [15]. In
combination, the various sources may lead to an interesting and rich GW
spectrum with multiple peaks.

The spectrum of a stochastic background of GWs is characterized by its
fractional energy density, i.e. the energy density normalized to the critical
value, $\rho_{c}^{0}=3M_{\text{Pl}}^{2}H_{0}^{2}$, per logarithmic frequency
interval,

$\Omega_{\text{GW}}(f)=\frac{1}{\rho_{c}^{0}}\frac{d\rho_{\text{GW}}}{d\log
f}\,,$ (27)

where the total GW energy density is given by

$\rho_{\text{GW}}=\frac{M_{\text{Pl}}^{2}}{4}\left<\dot{h}_{ij}\dot{h}^{ij}\right>=\int\\!\frac{df}{f}\frac{d\rho_{\text{GW}}}{d\log
f}\,.$ (28)

Here, $h_{ij}$ denotes the GW metric perturbations, and the dot indicates the
derivative with respect to cosmic time $t$. Switching to conformal time
$\tau$, $dt=a\,d\tau$, where $a$ is the scale factor of the Universe, the
metric reads

$ds^{2}=a(\tau)\left\\{d\tau^{2}-\right[\delta_{ij}+h_{ij}(\mathbf{x},\tau)\,\left]dx^{i}dx^{j}\right\\}\,.$
(29)

During radiation domination, the Einstein equations in the linear regime for
the metric perturbations in momentum space using transverse-traceless gauge
become

$(\partial_{\tau}^{2}+k^{2})\,a(\tau)\,h_{ij}(\mathbf{k},\tau)=\frac{2\,a(\tau)}{M_{\text{Pl}}^{2}}\Pi_{ij}(\mathbf{k},\tau)\,,$
(30)

where $k=|\mathbf{k}|$ is the comoving wave number. The anisotropic stress
tensor $\Pi_{ij}$ relates to the energy-momentum tensor $T_{ij}$, via
$\Pi_{ij}(\mathbf{k},\tau)=\Lambda_{ij}^{ab}(\mathbf{k})T_{ab}(\mathbf{k},\tau)$,
where $\Lambda^{ab}_{ij}=P_{i}^{\,a}P_{j}^{\,b}-\frac{1}{2}P_{ij}P^{ab}$ with
$P_{ij}=\delta_{ij}-k_{i}k_{j}/k^{2}$ being the projector that extracts the
transverse and traceless part [1]. The equations of motion are then solved by
(neglecting the $a^{\prime\prime}$ term which vanishes in a radiation-
dominated universe i.e. for $a\propto\tau$)

$\hat{h}_{ij}(\mathbf{k},\tau)=\frac{2}{M_{\text{Pl}}^{2}}\,\int^{\tau}\\!\\!\\!d\tau^{\prime}\,\frac{a(\tau^{\prime})}{a(\tau)}\,\hat{\Pi}_{ij}(\mathbf{k},\tau)\,\mathcal{G}(k,\tau,\tau^{\prime})\,,$
(31)

where $\mathcal{G}(k,\tau,\tau^{\prime})=\sin[k(\tau-\tau^{\prime})]/k$ is the
causal Green’s function. We denote here the operator form of any quantity $Q$
by $\hat{Q}$.

###  III.1 Gravitational wave production

The GW energy density per logarithmic interval in the comoving momentum $k$ of
a generic stochastic source at conformal time $\tau$ is given by [1]

$\frac{d\rho_{\text{GW}}}{d\log
k}(k,\tau)=\frac{k^{3}}{4\pi^{2}M_{\text{Pl}}^{2}a^{4}(\tau)}\int\limits_{\tau_{i}}^{\tau}\\!\\!d\tau^{\prime}\\!\\!\int\limits_{\tau_{i}}^{\tau}\\!\\!d\tau^{\prime}a(\tau^{\prime})\,a(\tau^{\prime\prime})\,\times\\\
\cos[k(\tau^{\prime}-\tau^{\prime\prime})]\,\Pi^{2}(k,\tau^{\prime},\tau^{\prime\prime})\,,$
(32)

where $\tau_{i}$ is the time at which the GW source starts operating and
$\Pi^{2}(k,\tau^{\prime},\tau^{\prime\prime})$ is the unequal time correlator
of the anisotropic stress, defined as $\langle
0|\hat{\Pi}^{ab}(\mathbf{k},\tau)\hat{\Pi}^{\ast}_{ab}(\mathbf{k}^{\prime},\tau^{\prime})|0\rangle=(2\pi)^{3}\delta(\mathbf{k}-\mathbf{k}^{\prime})\Pi^{2}(k,\tau,\tau^{\prime})$.
In our case, the GWs are generated between reheating and reappearance; hence
$\tau_{i}=\tau_{\text{rh}}$ and $\tau\leq\tau_{\text{ra}}$. As the GWs
produced before the relaxion reaches its terminal velocity will, however, be
subdominant, we can take $\tau_{i}=\tau_{\text{pp}}$, so to first
approximation, the GW signature becomes independent of the temperature to
which the Universe was reheated.

The dark photon anisotropic stress sourcing the GWs can be written in terms of
the dark electric and magnetic fields as

$\hat{\Pi}_{ab}(\mathbf{k},\tau)=-\frac{\Lambda_{ab}^{ij}(\mathbf{k})}{a^{2}(\tau)}\int\\!\frac{d^{3}q}{(2\pi)^{3}}\,\big{[}\hat{E}_{i}(\mathbf{q},\tau)\hat{E}_{j}(\mathbf{k}-\mathbf{q},\tau)\,+\\\
\hat{B}_{i}(\mathbf{q},\tau)\hat{B}_{j}(\mathbf{k}-\mathbf{q},\tau)\big{]}\,.$
(33)

Focusing on the dominant modes which have completed their phase of maximal
tachyonic growth, $q\gtrsim r_{X}|\theta^{\prime}|/2$, we find that GWs with
momentum $k$ are dominantly produced at the time of reappearance. For
frequencies below the peak, both dark photon modes in Eq. 33 have momenta
close to the one that experiences maximal growth at reappearance,
$|\mathbf{q}|\sim|\mathbf{k}-\mathbf{q}|\sim
r_{X}|\theta^{\prime}_{\text{ra}}|/2$, whereas above the peak, one of the
contributing modes must have a larger momentum and therefore have exited the
tachyonic band before reappearance. Details of the calculation are deferred to
Appendix C. The present-day GW spectrum can then be written as

$\Omega_{\text{GW}}(f)=\Omega_{\text{GW}}^{\text{peak}}\,\mathcal{S}_{\xi}\left({f}/{f_{\text{peak}}}\right)$
(34)

with the peak frequency

$f_{\text{peak}}=\frac{k_{\text{ra}}}{2\pi\,a_{0}}=\frac{r_{X}|\theta^{\prime}_{\text{ra}}|}{2\pi\,a_{0}}=\frac{a_{\text{ra}}}{a_{0}}\frac{\xi\,H_{\text{ra}}}{2\pi}\,,$
(35)

where $k_{\text{ra}}$ is the mode that exits the tachyonic band at
reappearance, and the peak amplitude

$\\!\Omega_{\text{GW}}^{\text{peak}}=\frac{1}{\rho_{c}^{0}}\,\frac{a_{\text{ra}}^{4}}{a_{0}^{4}}\,\frac{1}{M_{\text{Pl}}^{2}H_{\text{ra}}^{2}}\,\frac{\Lambda_{\text{br}}^{8}}{4\,r_{X}^{2}}=\frac{1}{\rho_{c}^{0}}\,\frac{a_{\text{ra}}^{4}}{a_{0}^{4}}\,\frac{m_{\phi}^{4}f_{\phi}^{4}}{M_{\text{Pl}}^{2}\xi^{2}H_{\text{ra}}^{2}}\,.\\!$
(36)

We used Eq. 13 here as well as $\xi\simeq 2\delta r_{X}$ and
$\Lambda_{\text{br}}^{4}=m_{\phi}^{2}f_{\phi}^{2}/\delta$ to recast Eqs. 35
and 36 into the estimations in Eqs. 1 and 2. The spectral shape
$\mathcal{S}_{\xi}$ is given by

$\mathcal{S}_{\xi}(x)=\frac{1}{1+\frac{48}{5}(x-1)^{4}+\frac{$19\,965$}{256\,\xi^{2}}\left(x^{-3}+3\,x-4\right)}\,.$
(37)

Note that, similar to the oscillating axion case [20, 21, 46, 22], we obtain a
GW spectrum with an unpolarized low-frequency tail and a chiral spectrum above
the peak.

###  III.2 GW detection

Having predicted the GW spectrum generated by the relaxion dynamics, we can
now evaluate its detectability in current and future experiments. A stochastic
GW background can be detected in a given experiment if its signal-to-noise
ratio (SNR) $\varrho$ exceeds a threshold value $\varrho_{\text{thr}}$. The
SNR is given by [50]

$\varrho^{2}=T_{\text{obs}}\\!\\!\int\limits_{f_{\text{min}}}^{f_{\text{max}}}\\!\\!\\!df\left[\frac{\Omega_{\text{GW}}(f)}{\Omega_{n}(f)}\right]^{2}\,,$
(38)

where $T_{\text{obs}}$ is the period of observation,
$(f_{\text{min}},f_{\text{max}})$ is the frequency range of the detector, and
$\Omega_{n}(f)$ is the detector’s noise spectrum converted to fractional
energy density. For a cross-correlated measurement in a network of detectors,
as, for instance, a pulsar timing array (PTA), the noise spectrum has to be
replaced by the effective noise $\Omega_{\text{eff}}$ of the network (see Ref.
[50] for further details), and the SNR is given by Eq. 38 with an additional
factor of 2.

We present here the projected sensitivities of the planned space-based Laser
Interferometer Space Antenna (LISA) [51, 52], the planned Square-Kilometre
Array (SKA) [53] PTA observatory,666We assume here that the prospective
foreground from supermassive black hole binaries can be subtracted. Further
details on the sensitivities can be found in Ref. [54].  and the proposed
microhertz experiment $\mu$Ares [55]. In addition, we evaluate current limits
from the NANOGrav 11 year dataset [56]. Other GW observatories such as ground-
based interferometers or potential LISA-successor experiments in the decihertz
regime do not cover the frequency range in which a signal is expected in our
scenario.

Recently, NANOGrav has further reported strong evidence for a common-spectrum
stochastic process across the pulsars in their 12.5 year dataset [57], which
might be due to a GW background. However, a GW detection could not be
established yet, due to the lack of conclusive evidence regarding the
interpulsar correlations of this process. Nonetheless, we also present here an
interpretation of this signal as the GW background generated in our model,
fitting our signal to the NANOGrav data based on the procedure outlined in
Ref. [58].

Furthermore, GWs contribute to the total amount of radiation in the Universe
and are therefore subject to constraints on $N_{\text{eff}}$. This sets an
upper bound on the peak amplitude of the spectrum. Using Eq. 23, we obtain

$h^{2}\Omega_{\text{GW}}^{\text{peak}}\leq\frac{a_{\text{rec}}^{4}}{a_{0}^{4}}\frac{7\pi^{2}}{120}\left(\frac{4}{11}\right)^{\frac{4}{3}}\frac{T_{\text{rec}}^{4}}{\rho_{c}^{0}/h^{2}}\frac{\Delta
N_{\text{eff}}}{\mathcal{I}_{\mathcal{S}}(\xi)}\approx${10}^{-6}$,\\!$ (39)

where $\mathcal{I}_{\mathcal{S}}(\xi)=\int\\!d\log
x\,\mathcal{S}_{\xi}(x)\simeq$0.81.5$$ is the integral over the spectral shape
and $T_{\text{rec}}$ is the photon temperature at recombination.

An even stronger bound is obtained from the contribution of the dark photons
to $\Delta N_{\text{eff}}$. Since the GWs are sourced by dark photon
anisotropies, this directly leads to an upper bound on the GW amplitude which
is stronger than the one from the direct contribution of the GWs to $\Delta
N_{\text{eff}}$. Plugging Eq. 24 into Eq. 36, we obtain

$h^{2}\Omega_{\text{GW}}^{\text{peak}}\lesssim${10}^{-8}$\left(\frac{g_{\ast,\text{ra}}}{106.75}\right)^{\\!\frac{1}{3}}\left(\frac{25}{\xi}\right)^{2}\left(\frac{\Delta
N_{\text{eff}}^{X}}{0.52}\right)^{2},$ (40)

where the superscript $X$ indicates that this is the dark photon contribution
to $\Delta N_{\text{eff}}$. A similar constraint on the GW amplitude was
already observed in Ref. [58] in the audible axion scenario. There, it was
noted that the parameter region that provides the best fit to the stochastic
GW background that was potentially observed by the NANOGrav collaboration [57]
is in tension with constraints on $\Delta N_{\text{eff}}$.

##  IV Discussion

Figure 4: GW spectra for the parameter points listed in Table 1. The orange spectrum labeled “$\ast$” corresponds to our best fit to the NANOGrav 12.5 year dataset. The colored regions depict the projected power-law integrated sensitivity curves of LISA (red), $\mu$Ares (green), and SKA (turquoise and blue, assuming 5 and 20 years of observation, respectively) and the current exclusion from the NANOGrav 11 year data set (purple).  |  $m_{\phi}$ (eV) |  $f_{\phi}$ (GeV) |  $T_{\text{ra}}$ (GeV) | $\xi$
---|---|---|---|---
A | $3\text{\times}{10}^{-9}$ | ${10}^{14}$ | $0.02$ | 10
B | $2\text{\times}{10}^{-3}$ | ${10}^{12}$ | $1$ | 25
C1 | $5$ | ${10}^{8}$ | $1$ | 10
C2 | $5\text{\times}{10}^{-3}$ | ${10}^{11}$ | $1$ | 25
C3 | $5\text{\times}{10}^{-6}$ | ${10}^{14}$ | $1$ | 100
D | $0.1$ | ${10}^{12}$ | $30$ | 25
$\ast$ | ${10}^{-9}$ | ${10}^{14}$ | $0.02$ | 10
Table 1: Parameter values for the spectra shown in Fig. 4

In Fig. 4 we show some example GW spectra. The corresponding parameter values
are listed in Table 1. The coupling to dark photons is set to
$r_{X}=\xi/(2\delta)$. In addition, the projected sensitivities of LISA (red),
$\mu$Ares (green), and SKA with an observation period of 5 years (turquoise)
and 20 years (blue), represented by the corresponding power-law integrated
sensitivity curves [50], are depicted as shaded regions. The current
$95\text{\,}\%$ upper limits from the 11 year NANOGrav data set [56] are shown
in purple. The example spectra B – D are detectable with $\mu$Ares, whereas
the benchmark A is accessible via SKA and is excluded by NANOGrav. The
benchmark labeled “$\ast$” shown in orange corresponds to our best fit to the
NANOGrav 12.5 year dataset as discussed further below.

At low frequencies, the GW spectra behave approximately as $f^{3}$, in
accordance with the expectations based on causality arguments using that the
anisotropic stress of a causal source cannot be correlated at scales above the
horizon size at the time of production [59, 1]. At high frequencies, the
spectra fall approximately like $f^{-4}$, allowing a simple distinction from
the much steeper falling GW background generated from oscillating [20, 21, 46,
22] or rotating [60] axionlike fields. It should further be noted that, when
the peak position is fixed, changing $\xi$ barely affects the UV tail, while
the IR tail goes as $\xi^{2}$ (cf. Fig. 9), potentially allowing to one
disentangle the reappearance temperature and $\xi$ in the peak frequency Eq.
35 and thereby facilitating the determination of the relaxion parameters from
a hypothetical observed signal. Larger values of $\xi$ further result in a
flatter peak, although this may be an artifact of our analytic approximation;
cf. Eq. 37.

Figure 5: Values of the peak frequency and amplitude of the GW spectrum which
can be obtained in the relaxion DM scenario. The edges of the polygon
correspond to the minimal and maximal amplitudes which can be obtained for
$\xi=100$ (solid lines) and $\xi=10$ (dashed lines), considering the case when
the relaxion starts to oscillate immediately after barrier reappearance.

The range of peak frequencies and amplitudes that can be attained for
coherently oscillating relaxion DM is displayed in Fig. 5, considering the
case when the relaxion starts oscillating at the time the wiggles reappear.
The polygons are obtained from the allowed range for the cutoff $\Lambda$, for
$T_{\text{ra}}$ between approximately $240\text{\,}\mathrm{MeV}$ and $v_{H}$;
cf. Fig. 2. The solid and dashed lines assume $\xi=100$ and $\xi=10$,
respectively. Peak positions inside the polygons can be realized. The
sensitivities of $\mu$Ares, LISA and PTAs are again indicated as shaded
regions. Note that the shading and solid lines correspond to the sensitivity
for $\xi=100$. For $\xi=10$, the detection reach is degraded to the
correspondingly coloured dashed lines.

For a large part of the peak frequencies and amplitudes that can be realized
with relaxion DM, an observable signal is obtained, although mostly requiring
futuristic space-based interferometers such as $\mu$Ares for observation. For
low values of $\xi\sim 10$, a detection with SKA is possible. The present-day
sensitivity of NANOGrav and expected reach of LISA are, however, insufficient
for a detection. While NANOGrav is able to exclude parts of the parameter
space if the DM assumption is relaxed (cf. benchmark A in Fig. 4), the
sensitivity of LISA will remain insufficient even in this more general case.

We also show the GW sensitivity for relaxion DM as coloured regions in Fig. 3,
using the same colouring scheme as above. The (light) coloured regions
enclosed by the solid green and blue lines here indicate the relaxion masses
and mixing angles for which, at least in a subrange of the reappearance
temperatures in accordance with the constraints, an observable signal in
$\mu$Ares and SKA can be obtained. While $\mu$Ares covers the full low-mass
island as well as the range $m_{\phi}\lesssim$0.1\text{\,}\mathrm{eV}$$ in the
high-mass island, the sensitivity of SKA is limited to slightly lower DM
masses. Note that, as the GW spectrum strongly depends on the reappearance
temperature, a nonobservation in these experiments would in most cases not
rule out the coloured parameter space, as a detection can be evaded by
adjusting the reappearance temperature to shift the signal outside the
experiment’s reach. In the dark-green coloured region bounded by the dotted
lines, however, the generated stochastic GW background is observable in
$\mu$Ares for the full range of allowed reappearance temperatures,
guaranteeing a detectable signal in this region. The temperature dependence of
the relaxion constraints and GW limits can be seen explicitly in the animated
version of the figure that can be found in the ancillary material.

Last but not least, let us now dismiss the assumption that the relaxion
constitutes DM and return to Fig. 1, where we again indicate the parameter
regions in which the GW background can be observed in $\mu$Ares (green),
NANOGrav (purple), or SKA (blue/turquoise) by the respective colouring. The
coloured regions respect the $f_{\phi}<M_{\text{Pl}}$ constraint. If we allow
for super-Planckian decay constant, the regions extend to the dotted lines.

Regarding the GW sensitivity, the reader needs to be aware that the figure
shows the projection of a four-dimensional plot, as $T_{\text{ra}}$ and
$f_{\phi}$ (or $\Lambda$) are not fixed. While red and orange shading marks
the values of $m_{\phi}$ and $\sin\theta_{h\phi}$ for which it is not possible
to evade the respective constraints by adjusting the remaining parameters
(i.e. these coloured regions are definitely excluded), the GW contours (blue,
turquoise, purple and green) correspond to the maximal reach of the respective
experiments. They are obtained by taking the maximal SNR that can be achieved
in each experiment for the given values of $m_{\phi}$ and
$\sin\theta_{h\phi}$. Given the experimental sensitivities we assume here, a
detection via GWs can be evaded in all of the viable parameter space. In
particular, the purple colouring and lines do not indicate that the
corresponding parameter points are excluded by NANOGrav data but that NANOGrav
constrains the reappearance temperature (and the decay constant) in this
region. The same also applies to the projections for SKA and $\mu$Ares.
Furthermore, although the sensitivities overlap in the plot, $\mu$Ares and
PTAs operate in different frequency regimes and are therefore typically
sensitive to different ranges of $T_{\text{ra}}$. Nonetheless, we can conclude
that current and future GW experiments can potentially detect the stochastic
GWs from relaxion-stopping via dark photon emission, and thereby constrain the
parameter space. We again enclose an animated version of the figure in the
ancillary material available online, where the dependence on $T_{\text{ra}}$
(but not on $\Lambda$) is shown.

Figure 6: Best-fit point (cross) to the NANOGrav 12.5 year data as well as
$1\sigma$ (dark) and $2\sigma$ (light) regions, fixing $\xi$ to 10 (blue), 25
(orange), and 100 (green). The dotted lines enclose the region that can be
obtained without being in tension with $\Delta N_{\text{eff}}$ or BBN at the
respective value of $\xi$.

In addition to the current and prospective exclusion range, we also work out
the parameter range in which our model can account for the potential GW signal
recently observed in NANOGrav [57]. Fitting our spectrum to the data using the
same procedure as in Ref. [58], where we keep $\xi$ fixed and only fit the
peak frequency and amplitude, we obtain the best-fit points and the
corresponding $1\sigma$ and $2\sigma$ contours shown in Fig. 6. We further
indicate the minimal peak frequency dictated by the lower bound on the
reappearance temperature, $T_{\text{ra}}\gtrsim T_{\text{BBN}}$, as well as
the maximal peak amplitude consistent with the constraint on $\Delta
N_{\text{eff}}$ by dotted lines. While these bounds exclude an explanation of
the observed stochastic process in terms of our model for large values of
$\xi\sim 100$ (green), we can reach into the $1\sigma$ and $2\sigma$ regions
for intermediate $\xi\sim 25$ (orange), and for $\xi\sim 10$ (blue), we can
account for the NANOGrav best-fit point.

In the $\xi=10$ (lower) panel of Fig. 1, we indicate the values of $m_{\phi}$
and $\sin\theta_{h\phi}$ for which we can attain the best-fit point
($f_{\text{peak}}=$3.3\text{\,}\mathrm{nHz}$$,
$h^{2}\Omega_{\text{GW}}^{\text{peak}}=$1.2\text{\times}{10}^{-9}$$) by the
grey shaded region. Note that, as we fix $\xi$ and $f_{\text{peak}}$, this
also fixes the reappearance temperature to
$T_{\text{ra}}\sim$20\text{\,}\mathrm{MeV}$$, while the peak amplitude then
fixes $f_{\phi}$ as a function of $m_{\phi}$. Hence, in the grey shaded
region, $\Lambda$ can be adjusted within the constraints to obtain the
respective value of the mixing angle. The best-fit spectrum for $\xi=10$ is
also depicted in orange in Fig. 4.

##  V Conclusion

In this work, we have considered the possibility of probing the relaxion,
which was proposed to ameliorate the Higgs hierarchy problem, via
gravitational waves. A coupling to dark photons tames the relaxion excursion
after reheating and thus in turn opens up a large fraction of the parameter
space which was excluded in the minimal scenario without dark photons.
Furthermore, dark photon production after reheating can act as a source for
the generation of a stochastic gravitational wave background. The
gravitational waves are sourced by the anisotropies that are induced by the
tachyonic production of dark photons between the electroweak phase transition
and big bang nucleosynthesis. Hence, instead of the inflationary dynamics
responsible for solving the hierarchy problem, we are probing here the late-
time dynamics of the relaxion.

We have shown that this stochastic gravitational wave background can be probed
by various current (NANOGrav) and future (SKA and $\mu$Ares) gravitational
wave detectors. In addition, we also highlight the parameter range in which
our gravitational wave signal can account for the common-spectrum process
observed in the most recent NANOGrav data [57]. Alongside the existing
theoretical constraints, we have presented the relaxion parameter space which
can be detected or excluded by the gravitational wave observatories in Fig. 1.

We find that the spectral shape of the gravitational wave signal in our model
falls as the fourth power of the frequency above the peak, unlike a steeper
falling gravitational wave signals generated by other axionlike field dynamics
[20, 21, 22, 60], whereas it behaves like $f^{3}$ in the infrared, as expected
based on causality arguments. An observation of the spectrum in the range
around the peak should allow for a determination of the reappearance
temperature, while the amplitude can be used to determine the product of the
relaxion mass and decay constant.

Furthermore, we have shown that the relaxion can constitute dark matter in the
present Universe in the mass range of
$${10}^{-11}\text{\,}\mathrm{eV}$\lesssim
m_{\phi}\lesssim${10}^{-6}\text{\,}\mathrm{eV}$$ and
$${10}^{-2}\text{\,}\mathrm{eV}$\lesssim
m_{\phi}\lesssim$1\text{\,}\mathrm{eV}$$. While this scenario cannot be
constrained with current NANOGrav data, most of the dark matter parameter
space will be accessible by SKA and/or $\mu$Ares in the future. Hence, with
the advent of gravitational wave astronomy, we are now facing promising
prospects for probing the relaxation of the electroweak scale via the
stochastic gravitation wave background generated when stopping the relaxion,
independently of whether the relaxion constitutes dark matter or not.

###### Acknowledgements.

We are grateful to Hyungjin Kim for collaboration at the early stages of this
paper and participating in numerous vital discussions that seeded this work.
We thank Edoardo Vitagliano for pointing out the updated fifth force bounds.
The work of A.B. is supported by the Azrieli Foundation. The work of E.M. is
supported by the Minerva Foundation. The work of G.P. is supported by grants
from BSF-NSF, Friedrich Wilhelm Bessel research award, GIF, ISF, Minerva,
Yeda-Sela-SABRA-WRC, and the Segre Research Award. The work of W.R. and P.S.
is supported by the Deutsche Forschungsgemeinschaft (DFG), Project No.
438947057 and by the Cluster of Excellence “Precision Physics, Fundamental
Interactions, and Structure of Matter” (PRISMA+ EXC 2118/1) funded by the
German Research Foundation (DFG) within the German Excellence Strategy
(Project No. 39083149).

##  Appendix A Minimal relaxion scenario

Here, we discuss the minimal relaxion scenario where a coupling to dark
photons is absent. As discussed in Sec. II, due to the restoration of the EW
symmetry, the relaxion starts rolling again if the reheating temperature is
higher than the EW scale, $v_{H}$. To trap the relaxion in the minimum where
it stopped during inflation, we need to require that it is displaced by less
than $\Delta\theta\lesssim 2\delta$, where $2\delta$ is the separation between
the minimum and the maximum [7]. In the absence of a dark photon coupling, the
Hubble friction is sufficient to trap the relaxion if
$m_{\phi}\lesssim\sqrt{8}\,\delta H_{\text{ra}}$, where $H_{\text{ra}}$ is the
Hubble scale at the time of barrier reappearance [6]. In combination with the
constraints discussed in Eqs. 19, 20 and 21, this trapping condition limits
the mixing angle as,

$\displaystyle\sin\theta_{h\phi}$
$\displaystyle\gtrsim\frac{m_{\phi}^{3/2}\Lambda_{\rm
min}}{m_{h}^{2}(2\,H_{\text{ra}}^{2})^{1/4}}\,\,\,$ $\displaystyle({\rm
from}\,\,\Lambda\gtrsim\Lambda_{\rm min})\,,$ (41)
$\displaystyle\sin\theta_{h\phi}$
$\displaystyle\lesssim\frac{4\,M_{\text{Pl}}\,m_{\phi}H_{\text{ra}}}{m_{h}^{2}\,v_{H}}\,\,\,$
$\displaystyle({\rm from}\,\,f_{\phi}\lesssim M_{\text{Pl}})\,,$ (42)
$\displaystyle\sin\theta_{h\phi}$
$\displaystyle\lesssim\frac{2^{\frac{5}{4}}\,v_{H}\sqrt{m_{\phi}H_{\text{ra}}}}{m_{h}^{2}}\,\,\,$
$\displaystyle({\rm from}\,\,\Lambda_{\text{br}}\lesssim v_{H})\,,$ (43)

within which the relaxion can be trapped by the Hubble friction. Moreover, the
above equations restrict the relaxion mass for which the minimal scenario can
be realized to
$m_{\phi}\lesssim$2\text{\times}{10}^{-5}\text{\,}\mathrm{eV}$$.

In addition to these constraints, for the relaxion not to overclose the
Universe we need to require $\Omega_{\phi}\lesssim\Omega_{\text{DM}}$. This
further constrains the range of the relaxion mass, and thus the minimal
relaxion scenario can only be realized for

$\displaystyle
m_{\phi}\lesssim$5\text{\times}{10}^{-8}\text{\,}\mathrm{eV}$\,.$ (44)

Figure 7: Available parameter space in the minimal relaxion scenario when the
reheating temperature is higher than the EW scale. The red and orange shaded
regions are excluded by the indicated constraints or combinations thereof. The
blue shaded region is further excluded if the relaxion is not coupled to the
dark photon source term, leading to a field excursion of more than $2\delta$.
The dark photon production cannot occur in the region above the grey dashed
line. Relaxion DM can be realized within the solid and dashed black contours,
respectively, depending on whether a coupling to dark photons is assumed or
not. In the turquoise parameter region, the condition $\Delta\theta\lesssim
2\delta$ further leads to an overproduction of relaxion DM in the minimal
scenario.

Figure 7 shows the available parameter space in the relaxion mass $m_{\phi}$
vs. mixing angle $\sin\theta_{h\phi}$ plane. We assume here the minimal model
without coupling to dark photons and do not require the relaxion to account
for DM. Again, the red and orange shaded regions are excluded by the indicated
constraints. The condition that the relaxion remains trapped in the same
minimum in which it settled before reheating, $\Delta\theta\lesssim 2\delta$,
then additionally excludes the region shaded in blue, with only the white
region remaining as viable parameter space. In the turquoise shaded region,
fulfilling $\Delta\theta\lesssim 2\delta$ without dark photons further leads
to an overclosure of the Universe. The region in which the relaxion can
constitute DM is enclosed by the dashed black contour. If the relaxion is
coupled to a dark photon, we further open up the blue and turquoise shaded
regions of the parameter space. With a dark photon, relaxion DM can be
realized within the solid black contour (same as the black framed regions of
Fig. 3). Note that the dark photon scenario can only be applied in the region
below the grey dashed line. Above the line, the condition
$H_{\text{pp}}>H_{\text{ra}}$ required for dark photon production is not
satisfied (cf. the discussion in Sec. II.3).

##  Appendix B Dark photon spectrum

As discussed in Sec. II.1, since $\theta^{\prime}>0$ in our case, positive-
helicity dark photon modes with $k<k_{\text{ex}}=r_{X}\theta^{\prime}$ are
tachyonic and experience exponential growth, whereas modes outside the
tachyonic band oscillate. In particular, the amplitude of the oscillating
modes does not change in time. At any time, the peak of the dark photon
spectrum is given by the mode that experiences maximal growth with
$k=k_{\text{m}}(\tau)=|r_{X}\theta^{\prime}|/2$. We therefore take the ansatz

$X_{+}(k,\tau)=\begin{cases}\mathcal{A}_{k}\cos\left(k\tau-\xi\right)&\text{for}\
k_{\text{m}}<k<k_{\text{pp}},\\\ 0&\text{otherwise},\end{cases}$ (45)

where we neglect the negative-helicity modes as well as all modes that did not
experience maximal growth yet, since these are exponentially suppressed. We
treat modes with $k_{\text{ex}}<k<k_{\text{m}}$ like the oscillating modes in
order to include the peak of the spectrum in our estimate.

Neglecting the energy produced before particle production, the dark photon
energy density is given by

$\rho_{X}=\int\limits_{k_{\text{m}}}^{k_{\text{pp}}}\\!\frac{dk\,k^{2}}{4\pi^{2}a^{4}}\left(|X^{\prime}|^{2}+k^{2}|X|^{2}\right)=\\!\int\limits_{k_{\text{m}}}^{k_{\text{pp}}}\\!\\!\frac{dk}{k}\,\frac{|\mathcal{A}_{k}|^{2}k^{5}}{4\pi^{2}a^{4}},$
(46)

where $k_{\text{pp}}=k_{\text{ex}}(\tau_{\text{pp}})$ is the mode that exits
the tachyonic band at particle production. On the other hand, we can determine
the energy density spectrum of the dark photons from the relaxion energy
density as

$\displaystyle\rho_{X}$
$\displaystyle=-\\!\int\limits_{\tau_{\text{pp}}}^{\tau}\\!d\eta\,\frac{a^{4}(\eta)}{a^{4}(\tau)}\frac{\partial
V}{\partial\theta}\,\frac{\partial\theta}{\partial\eta}\approx\int\limits_{k_{\text{m}}}^{k_{\text{pp}}}\\!\frac{dk}{k}\,\frac{k_{\text{m}}^{4}}{k^{4}}\frac{\xi\Lambda_{\text{br}}^{4}}{r_{X}}\,,$
(47)

where we assumed that at each time energy is dominantly transferred into the
maximally growing mode, which goes like $k_{\text{m}}\propto 1/\tau$. Hence,
we can rewrite $\mathcal{A}_{k}=\mathcal{A}_{X}k^{-9/2}$ with

$\mathcal{A}_{X}=\frac{\pi}{2}\,\Lambda_{\text{br}}^{2}\,\sqrt{\frac{\xi}{r_{X}}}\,a_{\text{ra}}^{2}k_{\text{ra}}^{2}\,,$
(48)

where $k_{\text{ra}}$ is the mode that exits the tachyonic band at
reappearance, $k_{\text{ra}}=\xi/\tau_{\text{ra}}$. The dark photon spectrum
therefore becomes

$X_{+}(k,\tau)=\begin{cases}\mathcal{A}_{X}k^{-\frac{9}{2}}\cos\left(k\tau-\xi\right)&\text{for}\
k_{\text{m}}<k<k_{\text{pp}},\\\ 0&\text{otherwise}.\end{cases}$ (49)

Figure 8: Simulated (coloured lines) and expected (dashed black line)
amplitude of the dark photon modes at different times. The deep (light)
coloured lines correspond to the positive (negative) helicity. At each time,
the expected peak momentum is indicated by the vertical dashed line in the
corresponding colour, whereas the black-dashed vertical line indicates the
upper bound $k<k_{\text{pp}}$ on the tachyonic dark photon momentum.

To corroborate our estimation, we have simulated the dark photon and relaxion
evolution after reheating, solving the equations of motion 7 and 8
numerically, using $5000$ logarithmically spaced, discretized momenta. In Fig.
8, we present the result of the numeric simulation for the dark photon modes
as a function of the momentum $k$. The simulation assumes
$m_{\phi}=$1\text{\,}\mathrm{meV}$$,
$f_{\phi}=$2.35\text{\times}{10}^{13}\text{\,}\mathrm{GeV}$$,
$\Lambda=$100\text{\,}\mathrm{TeV}$$, and $\xi=77$. The value of $r_{X}$ has
been determined numerically from Eq. 12 [6]. We show the spectrum at
$a/a_{\text{pp}}=4$ (red), $20$ (green), $100$ (orange), and $500$ (blue),
where the full and light coloured lines correspond to the positive and
negative helicity, respectively. As expected, the positive helicity modes
dominate over the negative helicity by far. Furthermore, the amplitude for
positive helicities indeed follows a $k^{-9/2}$ law [cf. Eq. 49, i.e.
$\sqrt{2k}|X_{+}|$ goes as $k^{-4}$] between the momentum $k_{\text{pp}}$ that
exits the tachyonic band at particle production and the peak momentum
$k_{\text{m}}(\tau)$ that experiences the largest growth rate at the
respective time. The peak momenta and $k_{\text{pp}}$ are indicated by the
coloured and black, dashed vertical lines, respectively.

##  Appendix C Calculation of the gravitational wave spectrum

The energy density of a stochastic GW background is given by

$\displaystyle\rho_{\text{GW}}(\mathbf{x},\tau)$
$\displaystyle=\frac{M_{\text{Pl}}^{2}}{4\,a^{2}(\tau)}\left<\hat{h}^{\prime}_{ab}(\mathbf{x},\tau)\hat{h}^{\prime
ab}(\mathbf{x},t)\right>$ (50)
$\displaystyle=\frac{M_{\text{Pl}}^{2}}{8\pi^{2}a^{2}(\tau)}\int\\!dk\,k^{2}\mathcal{P}_{h^{\prime}}(k,\tau)\,,$

where $\mathcal{P}_{h^{\prime}}$ is the power spectral density of
$h^{\prime}_{ab}$ defined by $\langle\hat{h}^{\prime
ab}(\mathbf{k},\tau)\hat{h}^{\prime\ast}_{ab}(\mathbf{k}^{\prime},\tau)\rangle=(2\pi)^{3}\delta(\mathbf{k}-\mathbf{k}^{\prime})\mathcal{P}_{h^{\prime}}(k,\tau)$.
We hence obtain the GW spectrum

$\frac{d\rho_{\text{GW}}}{d\log
k}(k,\tau)=\frac{M_{\text{Pl}}^{2}k^{3}}{8\pi^{2}a^{2}(\tau)}\mathcal{P}_{h^{\prime}}(k,\tau)\,.$
(51)

Using the solution Eq. 31, we obtain

$\mathcal{P}_{h^{\prime}}(k,\tau)=\frac{2}{M_{\text{Pl}}^{4}a^{2}(\tau)}\int\limits_{\tau_{\text{pp}}}^{\tau}\\!\\!d\tau^{\prime}\\!\\!\int\limits_{\tau_{\text{pp}}}^{\tau}\\!\\!d\tau^{\prime\prime}\,\cos[k(\tau^{\prime}-\tau^{\prime\prime})]\,\times\\\
a(\tau^{\prime})\,a(\tau^{\prime\prime})\,\Pi^{2}(k,\tau^{\prime},\tau^{\prime\prime})$
(52)

with the unequal time correlator again defined as $\langle
0|\hat{\Pi}^{ab}(\mathbf{k},\tau)\hat{\Pi}^{\ast}_{ab}(\mathbf{k}^{\prime},\tau^{\prime})|0\rangle=(2\pi)^{3}\delta(\mathbf{k}-\mathbf{k}^{\prime})\Pi^{2}(k,\tau,\tau^{\prime})$,
where we averaged over one period in $\tau$ to get rid of the
$\cos[k(2\tau-\tau^{\prime}-\tau^{\prime\prime})]$ part.

Inserting Eqs. 33 and 10 into Eqs. 52 and 51 and rewriting the dark electric
and magnetic fields in terms of the dark photon modes,

$\displaystyle\hat{\mathbf{E}}(\mathbf{k},\tau)$
$\displaystyle=\sum\limits_{\lambda}X^{\prime}_{\lambda}(k,\tau)\,\mathbf{\varepsilon}_{\lambda}(\mathbf{k})\,\hat{a}_{\lambda}(\mathbf{k})+\text{h.c.},$
(53a) $\displaystyle\hat{\mathbf{B}}(\mathbf{k},\tau)$
$\displaystyle=\sum\limits_{\lambda}\lambda\,k\,X_{\lambda}(k,\tau)\,\mathbf{\varepsilon}_{\lambda}(\mathbf{k})\,\hat{a}_{\lambda}(\mathbf{k})+\text{h.c.},$
(53b)

we obtain the spectrum at the time of reappearance as

$\frac{d\rho_{\text{GW}}}{d\log
k}=\frac{k^{3}}{2\pi^{2}M_{\text{Pl}}^{2}a_{\text{ra}}^{4}}\int\\!\frac{d^{3}q}{(2\pi)^{3}}\,|\Sigma_{ab}(\mathbf{k},\mathbf{q})|^{2}\times\\\
\left(|I_{c}(\mathbf{k},\mathbf{q})|^{2}+|I_{s}(\mathbf{k},\mathbf{q})|^{2}\right)\,,$
(54)

where we used
$[\hat{a}_{\lambda}(\mathbf{k}),\hat{a}_{\lambda^{\prime}}^{\dagger}(\mathbf{k}^{\prime})]=(2\pi)^{3}\delta_{\lambda\lambda^{\prime}}\delta(\mathbf{k}-\mathbf{k}^{\prime})$
and the invariance of the integrand upon interchanging
$\mathbf{q}\to\mathbf{k}-\mathbf{q}$. We only considered here the positive
helicity modes and defined
$\Sigma_{ab}(\mathbf{k},\mathbf{q})=\Lambda_{ab}^{ij}(\mathbf{k})\varepsilon_{i}^{+}(\mathbf{q})\varepsilon_{j}^{+}(\mathbf{k}-\mathbf{q})$
as well as

$I_{c/s}(\mathbf{k},\mathbf{q})=-\int\limits_{\tau_{\text{pp}}}^{\tau_{\text{ra}}}\\!\frac{d\tau}{a(\tau)}\,\left\\{\begin{array}[]{l}\cos(k\tau)\\\
\sin(k\tau)\rule{0.0pt}{12.0pt}\end{array}\right\\}\chi(q,l,\tau)$ (55)

with
$\chi(\mathbf{k},\mathbf{q},\tau)=X_{+}^{\prime}(q,\tau)X_{+}^{\prime}(l,\tau)+qX_{+}(q,\tau)\,lX_{+}(l,\tau)$
and $l=|\mathbf{k}-\mathbf{q}|$. Using the dark photon spectrum Eq. 49 we
obtain

$\chi(q,l,\tau)=\frac{\mathcal{A}_{X}^{2}\cos[(q-l)\tau]}{(ql)^{\frac{7}{2}}}$
(56)

for $\min(q,l)>k_{\text{m}}(\tau)$.

Expressing the momenta in terms of $k_{\text{ra}}$, we can rewrite the energy
density as

$\displaystyle\frac{d\rho_{\text{GW}}}{d\log k}$
$\displaystyle=\frac{\Lambda_{\text{br}}^{8}\xi^{2}}{r_{X}^{2}H_{\text{ra}}^{2}M_{\text{Pl}}^{2}}\
x^{3}\int\limits_{\frac{1}{2}}^{\infty}\\!dr\int\limits_{-1}^{1}\\!d\cos\theta\,\frac{|\Sigma_{ab}(x,r,\cos\theta)|^{2}}{512\,r^{5}s^{7}}\left(|\tilde{I}_{c}(x,r,s,\xi)|^{2}+|\tilde{I}_{s}(x,r,s,\xi)|^{2}\right)\
\Theta\left(s-1/2\right)\,,$ (57)

where the remaining integrals only depend on $x=k/k_{\text{ra}}$ and $\xi$. We
defined here $r=q/k_{\text{ra}}$ and $s=l/k_{\text{ra}}$. The polarization
part is given by (see, e.g., Ref. [20] for further details)

$\displaystyle|\Sigma_{ab}|^{2}=\sum\limits_{\lambda=\pm}\left[\frac{1+\lambda\cos\theta}{2}\right]^{2}\left[\frac{1+\lambda\frac{x-r\cos\theta}{s}}{2}\right]^{2}\\!\\!\\!,$
(58)

where $\lambda$ now denotes the GW helicity, and the time integrals evaluate
to

$\displaystyle\tilde{I}_{c/s}\\!=\\!\begin{cases}\\!\operatorname{Ci}\\!\left[(x+r-s)\xi\right]-\operatorname{Ci}\\!\left[\frac{(x+r-s)\xi}{2\,\min(r,s)}\right]\\!+(r\leftrightarrow
s)\\\
\operatorname{Si}\\!\left[(x+r-s)\xi\right]-\,\operatorname{Si}\\!\left[\frac{(x+r-s)\xi}{2\,\min(r,s)}\right]\\!+(r\leftrightarrow
s)\rule{0.0pt}{14.0pt}\end{cases}$ (59)

where $\operatorname{Ci}$ and $\operatorname{Si}$ are the cosine and sine
integral function,

$\displaystyle\operatorname{Ci}(z)=\int\limits_{z}^{\infty}\\!dt\,\frac{\cos
t}{t}\,,\qquad\operatorname{Si}(z)=\int\limits_{0}^{z}\\!dt\,\frac{\sin
t}{t}\,.$ (60)

The GW spectrum can now be obtained by evaluating Eq. 57 numerically. In this
paper, we, however, use an analytic approximation to the spectrum based on the
amplitude at the peak and the asymptotic behaviour in the UV and IR.

Since the GW momentum is given by the sum of two dark photon momenta,
$\mathbf{k}=\mathbf{q}+\mathbf{l}$, and since the time integrals are dominated
by the late-time behaviour, the peak of the GW spectrum will roughly be given
by twice the peak momentum of the dark photon spectrum at reappearance, i.e.
we take $k_{\text{peak}}=k_{\text{ra}}$. As the arguments of the cosine and
sine integrals in Eq. 59 are proportional to
$\xi\sim\mathcal{O}\left(10-100\right)$, we can expand for large $\xi$.
Assuming that the cosine and sine terms remaining in the expansion oscillate
quickly and therefore average to zero, the corresponding amplitude for $x=1$
then evaluates to777The analytic integration actually yields
$\frac{$278\,880\,340$+$44\,998\,983$\log 3}{$1\,254\,113\,280$}$ $\approx
0.262$ instead of $1/4$.

$\displaystyle\frac{d\rho_{\text{GW}}^{\text{peak}}}{d\log
k}=\left[\frac{1}{4\,\xi^{2}}+\mathcal{O}\left(\xi^{-4}\right)\right]\frac{\xi^{2}}{r_{X}^{2}}\frac{\Lambda_{\text{br}}^{8}}{H_{\text{ra}}^{2}M_{\text{Pl}}^{2}}\,.$
(61)

For $k\gg k_{\text{ra}}$, we have
$x+r-s=r(1+\cos\theta)+\mathcal{O}\left(x^{-1}\right)$ and
$x-r+s=2x+\mathcal{O}\left(1\right)$. The $(r\leftrightarrow s)$ part of Eq.
59 is hence suppressed by $1/x$, and only the first two terms contribute.
Again expanding for large $\xi$, we obtain that the spectrum behaves as

$\displaystyle\frac{d\rho_{\text{GW}}^{\text{UV}}}{d\log
k}=\left[\frac{5}{192\,\xi^{2}}+\mathcal{O}\left(\xi^{-4}\right)\right]\frac{\xi^{2}}{r_{X}^{2}}\frac{\Lambda_{\text{br}}^{8}}{H_{\text{ra}}^{2}M_{\text{Pl}}^{2}}\frac{k_{\text{ra}}^{4}}{k^{4}}$
(62)

in the UV. For $k\ll k_{\text{ra}}$, on the other hand, the arguments of the
cosine and sine integrals become proportional to $x\pm r\mp
s=x(1\pm\cos\theta)+\mathcal{O}\left(x^{2}\right)$, so we can expand the
integrals for low arguments, yielding the IR behaviour

$\displaystyle\frac{d\rho_{\text{GW}}^{\text{IR}}}{d\log
k}=\frac{64}{$19\,965$}\frac{\xi^{2}}{r_{X}^{2}}\frac{\Lambda_{\text{br}}^{8}}{H_{\text{ra}}^{2}M_{\text{Pl}}^{2}}\frac{k^{3}}{k_{\text{ra}}^{3}}\,.$
(63)

As can be seen from Eq. 58, at low momentum, $x\ll 1$, where we have $s\sim
r$, both GW helicities contribute equally, whereas at high momentum, where
$x\gg 1$ and hence $s\sim x$, the negative helicity is suppressed by a factor
of $1/k^{2}$ compared to the positive one.

Our approximation to the full spectrum is now obtained by combining these
results inversely. Neglecting GW production after reappearance, the energy
density will subsequently simply redshift as $\rho_{\text{GW}}\sim a^{-4}$, so
the spectrum today as a function of the comoving wave number $k$ becomes

$\displaystyle\Omega_{\text{GW}}(k)=\frac{1}{\rho_{c}^{0}}\left(\frac{a_{\text{ra}}}{a_{0}}\right)^{4}\\!\frac{1}{M_{\text{Pl}}^{2}H_{\text{ra}}^{2}}\frac{\Lambda_{\text{br}}^{8}}{4\,r_{X}^{2}}\left\\{1+\frac{$19\,965$}{256\,\xi^{2}}\left[\left(\frac{k_{\text{ra}}^{3}}{k^{3}}-1\right)+3\left(\frac{k}{k_{\text{ra}}}-1\right)\right]+\frac{48}{5}\left(\frac{k}{k_{\text{ra}}}-1\right)^{4}\right\\}^{-1},$
(64)

where the first term in the curly bracket reproduces our estimate for the peak
amplitude; the last term and the first part of the second term reproduce the
spectrum in the IR and UV, respectively; and the remaining term in the square
bracket renders $k=k_{\text{ra}}$ an extremum of the spectrum. Further noting
that frequency $f$ and comoving wave number $k$ are related via
$kd\tau=\frac{k}{a}dt=2\pi fdt$, we end up with the present-day spectrum in
Eqs. 34, 35, 36 and 37.

Figure 9: Spectral shape of the stochastic GW background for $\xi=10$ (blue),
$30$ (orange), and $100$ (green). The solid lines correspond to the numeric
calculation, whereas the dashed ones are the analytic approximation. The
dotted lines indicate the IR and UV behaviour, and the black cross marks our
estimate of the peak amplitude.

In Fig. 9, we compare the numerical result of the spectral shape based on Eq.
57 (solid lines) to the analytic approximation in Eq. 64 (dashed lines) for
different values of $\xi$. In addition, the UV and IR behaviour Eqs. 62 and 63
as well as the peak amplitude Eq. 61 are indicated by dotted lines and a
cross, respectively. Our approximation agrees well with the numeric result.

Figure 10: Simulation (full colours) and analytic approximation (light
colours) of the GW spectrum for a reappearance temperature of
$T_{\text{ra}}=$150\text{\,}\mathrm{GeV}$$ (blue) and
$750\text{\,}\mathrm{GeV}$ (orange). The vertical dashed lines indicate the
expected peak frequency, whereas the dot-dashed and dotted curves correspond
to the positive and negative helicity contributions to the simulated spectrum.

Finally, Fig. 10 shows the GW spectrum from our numeric simulation (cf. Fig.
8), assuming a reappearance temperature of
$T_{\text{ra}}=$150\text{\,}\mathrm{GeV}$$ (blue) and
$750\text{\,}\mathrm{GeV}$ (orange), corresponding to $a/a_{\text{ra}}=500$
and $a/a_{\text{ra}}=100$, respectively. We see that the GW background is
unpolarized below the peak, whereas above the peak the positive polarization
(dot-dashed) dominates over the negative one (dotted). Our estimates based on
Eq. 64 (light lines) overestimate the simulated spectra by an
$\mathcal{O}\left(1-10\right)$ factor.

## References

  * Caprini and Figueroa [2018] C. Caprini and D. G. Figueroa, Cosmological backgrounds of gravitational waves, Class. Quant. Grav. 35, 163001 (2018), arXiv:1801.04268 [astro-ph.CO] .
  * Caprini  _et al._ [2020] C. Caprini _et al._ , Detecting gravitational waves from cosmological phase transitions with LISA: an update, J. Cosmol. Astropart. Phys. 2020 (03), 024, arXiv:1910.13125 [astro-ph.CO] .
  * Graham  _et al._ [2015] P. W. Graham, D. E. Kaplan, and S. Rajendran, Cosmological Relaxation of the Electroweak Scale, Phys. Rev. Lett. 115, 221801 (2015), arXiv:1504.07551 [hep-ph] .
  * Salopek  _et al._ [1989] D. S. Salopek, J. R. Bond, and J. M. Bardeen, Designing density fluctuation spectra in inflation, Phys. Rev. D 40, 1753 (1989).
  * Bödeker and Buchmüller [2020] D. Bödeker and W. Buchmüller, Baryogenesis from the weak scale to the grand unification scale, (2020), arXiv:2009.07294 [hep-ph] .
  * Banerjee  _et al._ [2019] A. Banerjee, H. Kim, and G. Perez, Coherent relaxion dark matter, Phys. Rev. D 100, 115026 (2019), arXiv:1810.01889 [hep-ph] .
  * Banerjee  _et al._ [2020] A. Banerjee, H. Kim, O. Matsedonskyi, G. Perez, and M. S. Safronova, Probing the relaxed relaxion at the luminosity and precision frontiers, J. High Energy Phys. 2020 (07), 153, arXiv:2004.02899 [hep-ph] .
  * Choi and Im [2016a] K. Choi and S. H. Im, Constraints on relaxion windows, J. High Energy Phys. 2016 (12), 093, arXiv:1610.00680 [hep-ph] .
  * Flacke  _et al._ [2017] T. Flacke, C. Frugiuele, E. Fuchs, R. S. Gupta, and G. Perez, Phenomenology of relaxion-Higgs mixing, J. High Energy Phys. 2017 (06), 050, arXiv:1610.02025 [hep-ph] .
  * Choi  _et al._ [2017] K. Choi, H. Kim, and T. Sekiguchi, Dynamics of the cosmological relaxation after reheating, Phys. Rev. D 95, 075008 (2017), arXiv:1611.08569 [hep-ph] .
  * Hook and Marques-Tavares [2016] A. Hook and G. Marques-Tavares, Relaxation from particle production, J. High Energy Phys. 2016 (12), 101, arXiv:1607.01786 [hep-ph] .
  * Fonseca  _et al._ [2020a] N. Fonseca, E. Morgante, R. Sato, and G. Servant, Axion fragmentation, J. High Energy Phys. 2020 (04), 010, arXiv:1911.08472 [hep-ph] .
  * Fonseca  _et al._ [2020b] N. Fonseca, E. Morgante, R. Sato, and G. Servant, Relaxion fluctuations (self-stopping relaxion) and overview of relaxion stopping mechanisms, J. High Energy Phys. 2020 (05), 080, [Erratum: JHEP 01, 012 (2021)], arXiv:1911.08473 [hep-ph] .
  * Wang [2019] S.-J. Wang, Paper-boat relaxion, Phys. Rev. D 99, 095026 (2019), arXiv:1811.06520 [hep-ph] .
  * Barducci  _et al._ [2021] D. Barducci, E. Bertuzzo, and M. A. Tupia, Gravitational tests of electroweak relaxation, J. High Energy Phys. 2021 (07), 119, arXiv:2011.05795 [astro-ph.CO] .
  * Khlebnikov and Tkachev [1997] S. Y. Khlebnikov and I. I. Tkachev, Relic gravitational waves produced after preheating, Phys. Rev. D 56, 653 (1997), arXiv:hep-ph/9701423 .
  * Easther and Lim [2006] R. Easther and E. A. Lim, Stochastic gravitational wave production after inflation, J. Cosmol. Astropart. Phys. 2006 (04), 010, arXiv:astro-ph/0601617 .
  * Dufaux  _et al._ [2009] J.-F. Dufaux, G. Felder, L. Kofman, and O. Navros, Gravity waves from tachyonic preheating after hybrid inflation, J. Cosmol. Astropart. Phys. 2009 (03), 001, arXiv:0812.2917 [astro-ph] .
  * Adshead  _et al._ [2018] P. Adshead, J. T. Giblin, and Z. J. Weiner, Gravitational waves from gauge preheating, Phys. Rev. D 98, 043525 (2018), arXiv:1805.04550 [astro-ph.CO] .
  * Machado  _et al._ [2019] C. S. Machado, W. Ratzinger, P. Schwaller, and B. A. Stefanek, Audible axions, J. High Energy Phys. 2019 (01), 053, arXiv:1811.01950 [hep-ph] .
  * Machado  _et al._ [2020] C. S. Machado, W. Ratzinger, P. Schwaller, and B. A. Stefanek, Gravitational wave probes of axionlike particles, Phys. Rev. D 102, 075033 (2020), arXiv:1912.01007 [hep-ph] .
  * Ratzinger  _et al._ [2021] W. Ratzinger, P. Schwaller, and B. A. Stefanek, Gravitational waves from an axion-dark photon system: A lattice study, SciPost Phys. 11, 001 (2021), arXiv:2012.11584 [astro-ph.CO] .
  * Choi  _et al._ [2014] K. Choi, H. Kim, and S. Yun, Natural inflation with multiple sub-Planckian axions, Phys. Rev. D 90, 023545 (2014), arXiv:1404.6209 [hep-th] .
  * Choi and Im [2016b] K. Choi and S. H. Im, Realizing the relaxion from multiple axions and its UV completion with high scale supersymmetry, J. High Energy Phys. 2016 (01), 149, arXiv:1511.00132 [hep-ph] .
  * Kaplan and Rattazzi [2016] D. E. Kaplan and R. Rattazzi, Large field excursions and approximate discrete symmetries from a clockwork axion, Phys. Rev. D 93, 085007 (2016), arXiv:1511.01827 [hep-ph] .
  * Giudice and McCullough [2017] G. F. Giudice and M. McCullough, A clockwork theory, J. High Energy Phys. 2017 (02), 036, arXiv:1610.07962 [hep-ph] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Kolb and Turner [1990] E. W. Kolb and M. S. Turner,  _The Early Universe_ , Frontiers in Physics, Vol. 69 (Westview Press, 1990).
  * Audren  _et al._ [2014] B. Audren, J. Lesgourgues, G. Mangano, P. D. Serpico, and T. Tram, Strongest model-independent bound on the lifetime of Dark Matter, J. Cosmol. Astropart. Phys. 2014 (12), 028, arXiv:1407.2418 [astro-ph.CO] .
  * Agrawal  _et al._ [2018] P. Agrawal, G. Marques-Tavares, and W. Xue, Opening up the QCD axion window, J. High Energy Phys. 2018 (03), 049, arXiv:1708.05008 [hep-ph] .
  * Kitajima  _et al._ [2018a] N. Kitajima, T. Sekiguchi, and F. Takahashi, Cosmological abundance of the QCD axion coupled to hidden photons, Phys. Lett. B 781, 684 (2018a), arXiv:1711.06590 [hep-ph] .
  * Kitajima  _et al._ [2021] N. Kitajima, J. Soda, and Y. Urakawa, Nano-Hz Gravitational-Wave Signature from Axion Dark Matter, Phys. Rev. Lett. 126, 121301 (2021), arXiv:2010.10990 [astro-ph.CO] .
  * O’Hare and Vitagliano [2020] C. A. J. O’Hare and E. Vitagliano, Cornering the axion with $CP$-violating interactions, Phys. Rev. D 102, 115026 (2020), arXiv:2010.03889 [hep-ph] .
  * Tan  _et al._ [2020] W.-H. Tan _et al._ , Improvement for Testing the Gravitational Inverse-Square Law at the Submillimeter Range, Phys. Rev. Lett. 124, 051301 (2020).
  * Bergé  _et al._ [2018] J. Bergé, P. Brax, G. Métris, M. Pernot-Borràs, P. Touboul, and J.-P. Uzan, MICROSCOPE Mission: First Constraints on the Violation of the Weak Equivalence Principle by a Light Scalar Dilaton, Phys. Rev. Lett. 120, 141101 (2018), arXiv:1712.00483 [gr-qc] .
  * Chen  _et al._ [2016] Y. J. Chen, W. K. Tham, D. E. Krause, D. López, E. Fischbach, and R. S. Decca, Stronger Limits on Hypothetical Yukawa Interactions in the 30–8000 nm Range, Phys. Rev. Lett. 116, 221102 (2016), arXiv:1410.7267 [hep-ex] .
  * Kapner  _et al._ [2007] D. J. Kapner, T. S. Cook, E. G. Adelberger, J. H. Gundlach, B. R. Heckel, C. D. Hoyle, and H. E. Swanson, Tests of the Gravitational Inverse-Square Law Below the Dark-Energy Length Scale, Phys. Rev. Lett. 98, 021101 (2007), arXiv:hep-ph/0611184 .
  * Schlamminger  _et al._ [2008] S. Schlamminger, K. Y. Choi, T. A. Wagner, J. H. Gundlach, and E. G. Adelberger, Test of the Equivalence Principle Using a Rotating Torsion Balance, Phys. Rev. Lett. 100, 041101 (2008), arXiv:0712.0607 [gr-qc] .
  * Hoskins  _et al._ [1985] J. K. Hoskins, R. D. Newman, R. Spero, and J. Schultz, Experimental tests of the gravitational inverse square law for mass separations from 2 to 105 cm, Phys. Rev. D 32, 3084 (1985).
  * Grifols  _et al._ [1989] J. A. Grifols, E. Masso, and S. Peris, Energy loss From the sun and RED giants: Bounds on short range baryonic and leptonic forces, Mod. Phys. Lett. A 4, 311 (1989).
  * Cadamuro and Redondo [2012] D. Cadamuro and J. Redondo, Cosmological bounds on pseudo Nambu-Goldstone bosons, J. Cosmol. Astropart. Phys. 2012 (02), 032, arXiv:1110.2895 [hep-ph] .
  * Raffelt [2012] G. Raffelt, Limits on a $CP$-violating scalar axion-nucleon interaction, Phys. Rev. D 86, 015001 (2012), arXiv:1205.1776 [hep-ph] .
  * Hardy and Lasenby [2017] E. Hardy and R. Lasenby, Stellar cooling bounds on new light particles: plasma mixing effects, J. High Energy Phys. 2017 (02), 033, arXiv:1611.05852 [hep-ph] .
  * Budnik  _et al._ [2019] R. Budnik, O. Davidi, H. Kim, G. Perez, and N. Priel, Searching for a solar relaxion or scalar particle with XENON1T and LUX, Phys. Rev. D 100, 095021 (2019), arXiv:1909.02568 [hep-ph] .
  * Riess  _et al._ [2018] A. G. Riess _et al._ , New parallaxes of galactic cepheids from spatially scanning the Hubble space telescope: Implications for the Hubble constant, Astrophys. J. 855, 136 (2018), arXiv:1801.01120 [astro-ph.SR] .
  * Salehian  _et al._ [2021] B. Salehian, M. A. Gorji, S. Mukohyama, and H. Firouzjahi, Analytic study of dark photon and gravitational wave production from axion, J. High Energy Phys. 2021 (05), 043, arXiv:2007.08148 [hep-ph] .
  * Soda and Urakawa [2018] J. Soda and Y. Urakawa, Cosmological imprints of string axions in plateau, Eur. Phys. J. C 78, 779 (2018), arXiv:1710.00305 [astro-ph.CO] .
  * Kitajima  _et al._ [2018b] N. Kitajima, J. Soda, and Y. Urakawa, Gravitational wave forest from string axiverse, J. Cosmol. Astropart. Phys. 2018 (10), 008, arXiv:1807.07037 [astro-ph.CO] .
  * Chatrchyan and Jaeckel [2021] A. Chatrchyan and J. Jaeckel, Gravitational waves from the fragmentation of axion-like particle dark matter, J. Cosmol. Astropart. Phys. 2021 (02), 003, arXiv:2004.07844 [hep-ph] .
  * Thrane and Romano [2013] E. Thrane and J. D. Romano, Sensitivity curves for searches for gravitational-wave backgrounds, Phys. Rev. D 88, 124032 (2013), arXiv:1310.5300 [astro-ph.IM] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ (LISA), Laser Interferometer Space Antenna, (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Robson  _et al._ [2019] T. Robson, N. J. Cornish, and C. Liu, The construction and use of LISA sensitivity curves, Class. Quant. Grav. 36, 105011 (2019), arXiv:1803.01944 [astro-ph.HE] .
  * Janssen  _et al._ [2015] G. Janssen _et al._ , Gravitational wave astronomy with the SKA, PoS AASKA14, 037 (2015), arXiv:1501.00127 [astro-ph.IM] .
  * Breitbach  _et al._ [2019] M. Breitbach, J. Kopp, E. Madge, T. Opferkuch, and P. Schwaller, Dark, cold, and noisy: Constraining secluded hidden sectors with gravitational waves, J. Cosmol. Astropart. Phys. 2019 (07), 007, arXiv:1811.11175 [hep-ph] .
  * Sesana  _et al._ [2019] A. Sesana _et al._ , Unveiling the Gravitational Universe at µHz Frequencies, (2019), arXiv:1908.11391 [astro-ph.IM] .
  * Arzoumanian  _et al._ [2018] Z. Arzoumanian _et al._ (NANOGRAV), The NANOGrav 11 year data set: Pulsar-timing constraints on the stochastic gravitational-wave background, Astrophys. J. 859, 47 (2018), arXiv:1801.02617 [astro-ph.HE] .
  * Arzoumanian  _et al._ [2020] Z. Arzoumanian _et al._ (NANOGrav), The NANOGrav 12.5 yr data set: Search for an isotropic stochastic gravitational-wave background, Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Ratzinger and Schwaller [2021] W. Ratzinger and P. Schwaller, Whispers from the dark side: Confronting light new physics with NANOGrav data, SciPost Phys. 10, 047 (2021), arXiv:2009.11875 [astro-ph.CO] .
  * Caprini  _et al._ [2009] C. Caprini, R. Durrer, T. Konstandin, and G. Servant, General properties of the gravitational wave spectrum from phase transitions, Phys. Rev. D 79, 083519 (2009), arXiv:0901.1661 [astro-ph.CO] .
  * Co  _et al._ [2021] R. T. Co, K. Harigaya, and A. Pierce, Gravitational waves and dark photon dark matter from axion rotations, (2021), arXiv:2104.02077 [hep-ph] .