Post-inflationary axions: a minimal target for axion haloscopes

Marco Gorghettoa and Edward Hardyb

a Department of Particle Physics and Astrophysics, Weizmann Institute of
Science,
Herzl St 234, Rehovot 761001, Israel

b Department of Mathematical Sciences, University of Liverpool,
Liverpool, L69 7ZL, United Kingdom

An axion-like-particle (ALP) in the post-inflationary scenario with domain
wall number $N>1$ can be dark matter if the residual $\mathbb{Z}_{N}$ symmetry
has a small explicit breaking. Although we cannot determine the full dynamics
of the system reliably, we provide evidence that such an ALP can account for
the observed dark matter abundance while having a relatively small decay
constant and consequently a possibly large coupling to photons. In particular,
we determine the number of domain walls per Hubble patch around the time when
they form using numerical simulations and combine this with analytic
expectations about the subsequent dynamics. We show that the strongest
constraint on the decay constant is likely to come from the dark matter ALPs
being produced with large isocurvature fluctuations at small spatial scales.
We also comment on the uncertainties on the dark matter small-scale structure
that might form from these overdensities, in particular pointing out the
importance of quantum pressure in the $N=1$ case.

###### Contents

  1. 1 Introduction
  2. 2 Early universe evolution
  3. 3 Comparison with numerical results
  4. 4 Observational implications
  5. 5 Discussion and future directions
  6. A More details on the pre-inflationary scenario
  7. B The relic abundance from strings
  8. C Simulations

##  1 Introduction

By virtue of being the most robust known solution to the strong CP problem,
the QCD axion is one of the best-motivated extensions to the Standard Model
(SM) [1, 2, 3, 4, 5, 6, 7]. It is also plausible that one or more axion-like-
particles (ALPs), i.e. pseudo-Nambu-Goldstone bosons of spontaneously broken
approximate global U(1) symmetries, exist (we will use the name ‘axion’ to
refer to both ALPs and QCD axions). For instance, current understanding
suggests that axions might be common in many string theory constructions [8,
9, 10, 11, 12]. Additionally, axions are good candidates to be dark matter
(DM) [13, 14, 15]: for masses $\lesssim{\,{\rm eV}}$ and typical interaction
strengths, they are stable on cosmological timescales [16, 17], and a
nonrelativistic relic population is inevitably produced in the early Universe.

Motivated by these appealing features, numerous axion dark matter detection
experiments have been proposed targetting different ranges of axion masses
[18, 19, 20]. Many of these rely on the axion having a coupling to photons

$\mathcal{L}\supset\,\frac{1}{4}g_{a\gamma\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}~{},$
(1)

where $a$ is the axion, $F_{\mu\nu}$ is the electromagnetic field strength,
$\tilde{F}_{\mu\nu}=\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}F^{\rho\sigma}$ is
its dual, and $g_{a\gamma\gamma}$ is the axion-to-photon coupling constant.
For the QCD axion this coupling is

$\displaystyle g_{a\gamma\gamma}$ $\displaystyle=\frac{\alpha_{em}}{2\pi
f_{a}}\left(\frac{E}{N}-1.92\right)$ (2) $\displaystyle=2.0\cdot
10^{-16}{\,{\rm
GeV}}^{-1}\left(\frac{E}{N}-1.92\right)\left(\frac{m_{a}}{10^{-6}{\,{\rm
eV}}}\right)~{},$

where $f_{a}$ is the axion decay constant and $\alpha_{em}$ is the fine-
structure constant, and the second line follows from the fixed relation
between the QCD axion mass and decay constant. $E$ and $N$ are the anomaly
coefficients of the Peccei-Quinn (PQ) symmetry with respect to
electromagnetism and QCD respectively [21, 22, 23]. Meanwhile, the term in eq.
(2) that is independent of $E$ comes from the mixing of the axion with the
neutral pion, ultimately due to the PQ anomaly with respect to QCD and the
fact that quarks are charged. Unless $E$ is very large, the couplings in eq.
(2) are challenging to reach experimentally.111See however e.g. [24, 25, 26]
for ways the QCD axion coupling to photons can be enhanced. This difficulty
motivates understanding the plausible values of the axion-to-photon coupling
in other theories that are more general than the QCD axion but which might not
solve the strong-CP problem.

An ALP’s coupling to photons $g_{a\gamma\gamma}$ does not get a contribution
analogous to the constant part in eq. (2) (because the U(1) symmetry that
gives rise to it is not anomalous under QCD) so might simply vanish.
Nevertheless, it is possible that the new physics associated to the ALP
carries an anomaly with respect to electromagnetism, so $g_{a\gamma\gamma}$
could plausibly be non-zero. Similarly to eq. (2), we parametrise such a
coupling by

$g_{a\gamma\gamma}=\frac{\alpha_{em}}{2\pi f_{a}}\frac{E}{N}~{}.\ $ (3)

In minimal theories $E$ and $N$ are determined by the anomaly coefficients of
the U(1) symmetry with respect to QED and a new confining gauge group, and in
many concrete models $E/N\simeq\mathcal{O}(1)$. For convenience we further
define the UV model-independent ‘coupling’ $\tilde{g}_{a\gamma\gamma}\equiv
g_{a\gamma\gamma}/(E/N)=\alpha_{em}/(2\pi f_{a})$; the relation between this
and the true coupling ${g}_{a\gamma\gamma}$ depends on the underlying theory.
Unlike the QCD axion, an ALP’s decay constant is independent of its mass. As a
result, for fixed values of the mass and $E/N\neq 0$, the axion-to-photon
coupling can be large if $f_{a}$ is small enough. However, if the axion is to
comprise the full dark matter abundance in a standard cosmological history,
$f_{a}$ has a lower bound, ultimately because the energy density in the axion
field available to produce dark matter is proportional to $f_{a}^{2}$.

There are two broad axion cosmological scenarios (see [27, 28] for an
intermediate regime). In the pre-inflationary scenario the global $U(1)$
symmetry that gives rise to the axion is broken during inflation and is never
subsequently restored. In this case, the axion field initially takes a
constant value over the observable Universe, and the axion relic abundance
$\Omega_{a}$ is produced by the misalignment mechanism. If the ALP’s mass is
temperature-independent

$\frac{\Omega_{a}^{\rm mis}}{\Omega_{\rm DM}}\simeq 2.2\cdot
10^{-3}\left(\frac{f_{a}}{10^{12}~{}{\,{\rm
GeV}}}\right)^{2}\left(\frac{m_{a}}{10^{-6}~{}{\,{\rm
eV}}}\right)^{1/2}h(\theta_{0})\,\theta_{0}^{2}~{},$ (4)

where $\Omega_{\rm DM}$ is the observed DM abundance,
$\theta_{0}\in(-\pi,\pi]$ is the initial misalignment angle, and $h$ is a
function that accounts for anharmonicities in the axion potential. For typical
axion potentials $h\left(\theta_{0}\right)\simeq 1$ for $\theta_{0}\lesssim 1$
and increases slowly as $\theta_{0}$ approaches $\pi$. Constraints on the size
of primordial isocurvature perturbations give a lower bound on
$\pi-|\theta_{0}|$ that depends on the Hubble scale during inflation, see
Appendix A. For instance, for a cosine potential $h(\theta_{0})$ increases
only logarithmically with $1/(\pi-|\theta_{0}|)$, and $H_{I}\simeq 1{\,{\rm
eV}}$ requires $\pi-|\theta_{0}|\gtrsim 10^{-18}$ which gives
$h(\pi-10^{-18})\lesssim 200$. From eq. (4), if $\Omega^{\rm
mis}_{a}=\Omega_{\rm DM}$ the axion-to-photon coupling is constrained to be
222An ALP’s misalignment relic abundance, and consequently the maximum
$\tilde{g}_{a\gamma\gamma}$, can be larger if the axion mass is temperature
dependent. This is also possible in theories with multiple axions [29],
different cosmological histories, and different axion potentials [30].

$\tilde{g}_{a\gamma\gamma}\lesssim 5\cdot 10^{-17}{\,{\rm
GeV}}^{-1}\left(\frac{m_{a}}{10^{-6}{\,{\rm
eV}}}\right)^{1/4}h(\theta_{0})^{1/2}\,\theta_{0}~{}.$ (5)

In the post-inflationary scenario the $U(1)$ symmetry is instead unbroken at
the end of inflation or is subsequently restored e.g. by finite temperature
effects. As a result, the field is initially inhomogeneous over our observable
Hubble patch and the axion relic abundance is determined by the decay of a
network of topological defects: global cosmic strings and domain walls. In
this paper, we argue that an ALP in the post-inflationary scenario can be all
the DM with a much smaller $f_{a}$ (i.e., for fixed $E/N$, a larger coupling
to photons) than in the pre-inflationary scenario while being compatible with
observational constraints, the most important of which is expected to come
from bounds on primordial isocurvature perturbations.333A related analysis of
the same theory can be found in [31, 32] focusing on the gravitational waves
from domain walls, see also [33] for a similar study of a different theory
that has an axion. This is possible because, if $N>1$, the string-domain wall
network might be long-lived, which results in the dark matter produced by its
decay being redshifted less. We also analyse the impact of the isocurvature
constraints on possible gravitational wave signals from topological defects
and discuss the properties of, and the uncertainties in, the dark matter
substructure that might form in the early Universe.

We stress that it is difficult to analyse the system’s dynamics from first
principles because they are both non-linear and involve large scale
separations, and as a result the field evolution is presently unknown.
Nevertheless, we make progress by combing results from numerical simulations
with reasonable standard assumptions about the dynamics of the domain wall
network (which we can only very partially check). In particular, in Section 3
we determine number of domain walls per Hubble patch $\mathcal{A}_{\star}$ at
the time when the axion mass becomes cosmologically relevant (when the Hubble
parameter is $H_{\star}$), which is important for the resulting dark matter
abundance, and provide evidence that the axions emitted from the string-domain
wall system are, at most, mildly relativistic. Our results for
$\mathcal{A}_{\star}$, applicable for both ALPs and the QCD axion, suggest
that this quantity increases logarithmically with $f_{a}/H_{\star}$, and mean
that an extrapolation of simulation results is needed. We also discuss the
numerous remaining uncertainties on the evolution of $\mathcal{A}$ after this
time, and on the spatial distribution of the axions produced from the string-
domain wall system.

The structure of the paper is as follows. In Section 2 we introduce the post-
inflationary scenario for ALPs and discuss analytic expectations about the
dynamics and the dark matter relic abundance, and in Section 3 we present
results from numerical simulations of the string-wall system. In Section 4 we
analyse the observational constraints on the theory and discuss its relevance
for upcoming axion dark matter searches, and we also comment on the possible
gravitational waves signals and dark matter substructure. In Section 5 we
describe future directions and extensions.

##  2 Early universe evolution

An axion-like-particle (ALP) is the pseudo-Nambu Goldstone boson, $a$, of an
approximate global U(1) symmetry that is spontaneously broken at a high scale
$v$. Because $a$ is an angular variable, the theory preserves an exact
discrete shift symmetry $a\rightarrow a+2\pi v$. The simplest model in which
this happens consists of a complex scalar $\phi$ with an approximate global
U(1) symmetry and potential

$V_{\phi}(\phi)=\frac{m_{r}^{2}}{2v^{2}}\left(|\phi|^{2}-v^{2}\right)^{2}+\dots~{},$
(6)

that induces a vacuum expectation value $|\langle\phi\rangle|=v$. Expanding
$\phi=\frac{1}{\sqrt{2}}(r+v)e^{ia/v}$, the axion $a$ is a Goldstone mode,
while the radial mode $r$ has (large) mass $m_{r}$, which is often of order
$v$.

We assume that the U(1) symmetry is broken by non-perturbative effects, e.g.
because it is anomalous with respect to a new confining gauge group. Such a
breaking leads to a potential for $a$ that is invariant under $a\rightarrow
a+2\pi f_{a}$. In general $f_{a}$ and $v$ are related by $v=Nf_{a}$, where $N$
has to be an integer otherwise the symmetry $a\rightarrow a+2\pi v$ would be
violated (in KSVZ-like models, $N$ is the number of fermions charged under
U(1)). As a result, the potential induced by a particular non-perturbative
contribution does not necessarily break the U(1) symmetry completely, and can
leave a subgroup $\mathbb{Z}_{N}$ unbroken. For instance, the total axion
potential might be

$V(a)=m_{a}^{2}f_{a}^{2}\left[1-\cos(a/f_{a})\right]+\delta V_{PQ}(a)~{},$ (7)

where the first term comes from the non-perturbative contribution.

In eq. (7), $\delta V_{PQ}(a)$ represents additional sources of U(1) breaking
that violate the discrete symmetry $a\rightarrow a+2\pi f_{a}$ while still
preserving $a\rightarrow a+2\pi v$. Such breaking is expected because it is
believed that there are no exact global symmetries in quantum gravity [34, 35,
36] (implications for the QCD axion have been studied in e.g. [37, 38, 39, 40,
41, 42, 43, 44, 45, 46]).444This happens e.g. via higher dimensional operators
of dimension $n+m$ of the form
$\alpha\phi^{n}(\phi^{\dagger})^{m}/M_{p}^{n+m-4}$, where $M_{p}$ is the
Planck Mass and $n\neq m$. If $\alpha=\mathcal{O}(1)$, $n+m$ needs to be large
in order for the U(1) to be an approximate symmetry. On the other hand, if
$\alpha$ is exponentially suppressed (as expected from some non-perturbative
gravitational processes such as wormholes [34, 35]) $n+m$ can be small. As an
example, although the particular form does not matter, one can consider

$\delta V_{PQ}(a)=\epsilon f_{a}^{4}\cos(a/(Nf_{a})-\theta_{1})~{},$ (8)

with $\theta_{1}=\mathcal{O}(1)$. In the absence of $\delta V_{PQ}$ there
would be $N$ inequivalent degenerate vacua, but the energy densities of these
vacua are shifted by $\simeq\epsilon f_{a}^{4}$ due to by $\delta V_{PQ}$. We
demand $\epsilon<m_{a}^{2}/f_{a}^{2}$ so that the global minimum of the full
potential is a perturbation of one of the $\mathbb{Z}_{N}$ preserving minima
and the axion mass is dominated by the first term in eq. (7) [47]. In the
following we will take $\theta_{1}=0$, but this choice will not have any
effect.

The equations of motion (EoM) from the potentials in eqs. (6) and (7) admit
topologically non-trivial field configurations called strings and domain walls
respectively [48, 49, 50, 51] (see [52] for a review). These form if the field
is inhomogeneous, as is the case in the post-inflationary scenario. We assume
a standard cosmological history, and, as we discuss subsequently, in all
theories of interest the strings and domain walls are destroyed before matter-
radiation equality (MRE). During the preceding era of radiation domination the
Hubble parameter is $H\equiv\dot{R}/R=1/(2t)$, where $R$ is the scale factor
of the Universe and $t$ is the cosmic time. Unless otherwise stated we assume
that the axion potential is independent of the temperature of the Universe.

###  2.1 The scaling regime

After the PQ symmetry spontaneously breaks, $\phi$ has initially random
spatial fluctuations. The axion potential in eq. (7) is cosmologically
irrelevant while $H\gg m_{a}$ so a network of cosmic strings forms. The
strings occur in regions of space in which the axion U(1) is wrapped non-
trivially, and they have (evolving) energy per unit length, i.e. tension,

$\mu=\pi v^{2}\log(m_{r}/H)~{}.$ (9)

The string network is driven to an attractor, scaling, solution on which there
is approximately one Hubble length of string per Hubble volume. To study this
regime, numerical simulations of the nonlinear EoM from the potential in eq.
(6) are required and an extrapolation in $\log(m_{r}/H)$ is needed. In
particular, simulations can only be carried out with a relatively small value
of $\log(m_{r}/H)\lesssim 8\div 9$ so the scale separation between the string
thickness (of order $m_{r}^{-1}$) and the Hubble distance is far from the
values relevant shortly before $H=m_{a}$ [53].

We define the number of strings per Hubble volume
$\xi=\lim_{\mathcal{V}\rightarrow\infty}\ell_{\rm st}t^{2}/\mathcal{V}$ where
$\ell_{\rm st}$ is the length of string in a volume $\mathcal{V}$. Simulations
show that during the scaling regime $\xi$ has a logarithmic increase,
$\xi\simeq 0.24\log(m_{r}/H)$. The energy density in (long) strings is, from
eq. (9), $\rho_{\rm s}=4\xi\mu H^{2}$ and the network releases energy density
into (relativistic) axion waves at a rate $\Gamma_{a}^{\rm
st}\simeq\xi\mu/t^{3}$. Additionally, simulations suggest that at the relevant
large $\log(m_{r}/H)$ most of the axions are emitted with momentum of order
Hubble. The resulting axion number density released until time $t$ is given by
$n_{a}^{\rm st}(t)\simeq\xi\mu H$, which is dominantly produced in the final
few Hubble times prior to $t$, see Appendix B.555See however [54] which
suggests a scale invariant instantaneous emission spectrum, which would reduce
the relic abundance by a factor of $\log_{\star}$ relative to eq. (10).

The scaling regime continues until the time when the Hubble parameter $H$
equals the axion mass, $H=m_{a}$, when the first ($\mathbb{Z}_{N}$-preserving)
contribution to the axion potential in eq. (7) starts to be cosmologically
relevant. We indicate quantities at that moment by $\star$, e.g. $H_{\star}$
is the Hubble parameter at this time, and we define
$\log_{\star}\equiv\log(m_{r}/H_{\star})\simeq\log(f_{a}/m_{a})$, which is of
$\mathcal{O}(100)$ for the $(m_{a},f_{a})$ that we will consider.

The axion waves emitted by strings prior to $H_{\star}$ become nonrelativistic
after $H=H_{\star}$ and contribute to the dark matter relic abundance, see
Section 4.1 of [55] and [56] for more details. Evaluating $n_{a}^{\rm st}$ at
$H=H_{\star}$ and redshifting to today, the resulting dark matter from strings
during the scaling regime, $\Omega_{a}^{\rm st}$, is, for a temperature-
independent $m_{a}$,

$\frac{\Omega_{a}^{\rm st}}{\Omega_{\rm DM}}\simeq
3\left(\frac{\xi_{\star}\log_{\star}}{10^{3}}\right)\left(\frac{Nf_{a}}{10^{12}~{}{\,{\rm
GeV}}}\right)^{2}\left(\frac{m_{a}}{10^{-6}~{}{\,{\rm eV}}}\right)^{1/2}~{}.$
(10)

For the axion masses and decay constants of interest, $10^{-15}~{}{\,{\rm
eV}}\lesssim m_{a}\lesssim{\,{\rm eV}}$ and $10^{8}~{}{\,{\rm GeV}}\lesssim
f_{a}\lesssim 10^{12}~{}{\,{\rm GeV}}$, we have $\log_{\star}\simeq 40\div
80$. Therefore, assuming that the logarithmic growth in $\xi$ is valid also
for $\log(m_{r}/H)\gg 1$, $\xi_{\star}\log_{\star}\simeq 10^{3}$ (if
$\xi_{\star}$ take a different value the abundance can be calculated
accordingly from eq. (10)). Note that for larger values of
$\xi_{\star}\log_{\star}\gtrsim 10^{3}$ or $N\gg 1$, the axion number density
is substantially reduced between the time when the axions are relativistic and
when they become nonrelativistic, because the amplitude of the axion field is
much larger than $f_{a}$ and nonlinear effects from the axion potential become
relevant. As a result, the abundance $\Omega_{a}^{\rm st}$ decreases with
respect to eq. (10), see [55]. For $\xi_{\star}\log_{\star}\lesssim 10^{3}$
the effect is of order 1, comparable to the other uncertainties in eq.
(10).666The non-linear evolution is more important for the QCD axion because
$m_{a}$ has a strong temperature dependence [56].

###  2.2 The expected post-inflationary relic abundance

Once the axion potential is relevant at $H=m_{a}=H_{\star}$, domain walls
form, across which the axion field interpolates between vacua. Each string is
bounded by $N$ domain walls. The domain walls have a tension, i.e. energy per
unit area, given by $\sigma=\beta m_{a}f_{a}^{2}$, where $\beta$ is a
numerical factor ($\beta=8$ for a cosine and is similar for other typical
potentials).

#### $N=1$

If $N=1$ the domain walls connect the same vacuum of the potential in eq. (7),
$a/f_{a}=0$, and the string and domain wall network is unstable and expected
to decay. In the process axion waves are emitted, and subsequently the axion
field oscillates around the global minimum of its potential and behaves as
dark matter. The contribution to the axion abundance produced by the domain
walls and strings at this time is hard to determine reliably. If the network
is destroyed soon after $H=H_{\star}$, it might be relatively unimportant
compared to that from the scaling regime, $\Omega_{a}^{\rm st}$ in eq. (10),
because the energy in strings is enhanced relative to that in domain walls by
the logarithmically boosted string tension (see eq. (14) below). Including
only the contribution from the scaling regime, the value of $f_{a}$ for which
an $N=1$ post-inflationary ALP makes up all the DM is of same order as that in
the pre-inflationary scenario with an initial misalignment angle
$\theta_{0}\sim\pi-10^{-18}$ close to the top of the potential, see eq. (5)
and Figure 4. We do however stress that eq. (10) is a lower bound on the $N=1$
axion relic abundance; the domain wall network might take longer to decay than
expected, in which case it could produce a sizable additional contribution
resulting in smaller required $f_{a}$.

#### $N>1$

In contrast, if $N>1$ the domain walls interpolate between the inequivalent
vacua of the first term in eq. (7). In the absence of $\delta V_{PQ}$ they
would form a topologically stable network after $H=H_{\star}$ [51] which would
come to dominate the energy density in the Universe [57]. However, the energy
bias between the vacua induced by $\delta V_{PQ}$ is believed to lead to the
destruction of the network when $\delta V_{PQ}$ becomes cosmologically
relevant [51, 58, 59, 60]. This is thought to happen when the energy density
in domain walls, $\rho_{\rm w}\simeq\sigma H$, is approximately equal to the
energy density provided by the bias term, $\delta V_{PQ}\simeq\epsilon
f_{a}^{4}$, i.e. the Hubble parameter when the network is destroyed is
$H_{d}\simeq\epsilon f_{a}^{2}/m_{a}$. The energy released by the domain walls
is expected (but not proven) to dominantly go into dark matter axions, with a
small fraction emitted as gravitational waves [31].

If $\epsilon\simeq m_{a}^{2}/f_{a}^{2}$, the energy density difference between
the local minima is of $\mathcal{O}(m_{a}^{2}f_{a}^{2})$ and the
$\mathbb{Z}_{N}$ symmetry is badly broken. As a result, the evolution of the
string-domain wall network, and the axion relic abundance, is likely to be
similar to that in the $N=1$ scenario. However, if $\epsilon\ll
m_{a}^{2}/f_{a}^{2}$ the domain wall network is long-lived and is destroyed
when the Hubble parameter is $H_{d}\ll m_{a}$. Unfortunately, the dynamics in
this case are even harder to study with numerical simulations than the $N=1$
scenario. This is because there is another relevant scale, $m_{a}$, so a
second extrapolation, in $\log(m_{a}/H)$, is needed. Simulations can only
reach values $\log(m_{a}/H)\lesssim 5$, whereas we will see in Section 4.1
that values up to $\log(m_{a}/H_{d})\simeq 60$ (determined by $\delta V_{PQ}$)
are possible. The too-small $\log(m_{r}/H)$ are an issue not only because the
string tension depends on this factor, but also because the string dynamics
are dominated by emission of high energy modes with energy of order $m_{r}$
for $\log(m_{r}/H)\lesssim 9$ [56]. The second extrapolation is perhaps even
more problematic because, as $m_{a}/H$ increases, there is a growing mismatch
between the typical spatial size of the domain walls, which extend over an
Hubble length $\sim H^{-1}$, and the frequency of the axions they are expected
to emit (of order $m_{a}$). It is plausible that this could e.g. decrease the
ability of the domain walls to radiate axions as $m_{a}/H$ increases, which
might even qualitatively change the properties of the network when
$\log(m_{a}/H)\gg 1$.

Nevertheless, we proceed with some standard analytic expectations, which we
will subsequently partially check. It is useful to define the domain-wall area
per Hubble patch
$\mathcal{A}=\lim_{\mathcal{V}\rightarrow\infty}At/\mathcal{V}$, where $A$ is
the total domain wall area in a volume $\mathcal{V}$. $\mathcal{A}$ measures
parametrically the total area of domain walls in one Hubble patch in units of
the Hubble area $H^{-2}$, so the energy density in domain walls is

$\displaystyle\rho_{\rm w}=2\mathcal{A}\sigma H~{}.$ (11)

While $\delta V_{PQ}$ is cosmologically irrelevant, at $H<H_{\star}$ the
domain walls are thought to survive in a domain wall ‘scaling regime’:
$\mathcal{A}$ is expected to be approximately constant in time, and likely not
far from $\mathcal{A}\simeq\mathcal{O}(1)$ [61] (or slightly larger, because
of the logarithmic violation discussed in Section 3). To maintain this, the
network must emit energy density at a rate

$\Gamma_{a}^{\rm w}\simeq\rho_{\rm w}/(2t)~{},$ (12)

i.e. an energy density of order $\rho_{\rm w}(t)$ is released every Hubble
time.777This is because $\Gamma_{a}^{\rm w}(t)=[\dot{\rho}^{\rm free}_{\rm
w}(t_{0})-\dot{\rho}_{\rm w}(t_{0})]_{t_{0}=t}$, where $\rho^{\rm free}_{\rm
w}(t)=\mathcal{A}\sigma/\sqrt{t_{0}t}\propto R(t)^{-1}$ is the energy density
of a system of domain walls that coincides with that in the scaling regime at
$t=t_{0}$, but evolves freely, i.e. the domain walls are fixed in comoving
coordinates.  The momentum of the axions produced is unknown and could be set
by either of the two scales relevant to the domain walls: their size $\simeq
H^{-1}$ or their thickness $\simeq m_{a}^{-1}$. Although the latter is
suggested by results of simulations at small $\log(m_{a}/H_{d})$ and
$\log(m_{r}/H_{d})$, see Section 3, it is not clear whether this remains the
case for physical parameters. Regardless of this uncertainty, because the
emitted energy dominantly goes into at most semi-relativistic modes, the axion
number density produced until time $t$ is approximately
$m_{a}^{-1}\int^{t}_{t_{\star}}dt^{\prime}\Gamma(t^{\prime})(R(t^{\prime})/R(t))^{3}$,
i.e.

$\displaystyle n_{a}^{\rm w}(t)$ $\displaystyle\simeq\int_{\log
t_{\star}}^{\log t}\mathcal{A}\frac{\beta
f_{a}^{2}}{2t^{\prime}}\left(\frac{R(t^{\prime})}{R(t)}\right)^{3}~{}d\log
t^{\prime}$ (13) $\displaystyle=2\mathcal{A}\beta
m_{a}f_{a}^{2}\left[\left(\frac{m_{a}}{H}\right)^{\frac{1}{2}}-1\right]\left(\frac{R_{\star}}{R(t)}\right)^{3}$
(14)

where $R(t)\propto t^{1/2}$ is the scale factor. Because
$R(t^{\prime})^{3}/t^{\prime}\propto t^{\prime 1/2}$, it is clear from eq.
(13) that the number density is mostly produced at the latest times that the
network survives to and is dominated by the axions emitted during the Hubble
time immediately preceding the domain wall network’s destruction. If the
network survives up to $H\ll m_{a}$, the axion number density is
parametrically enhanced, by a factor $(m_{a}/H)^{1/2}$, with respect to that
from misalignment, which is $\mathcal{O}(m_{a}f_{a}^{2})$ at $t=t_{\star}$,
see eq. (14).

At $H=H_{d}$, all the domain walls are destroyed and their energy density
$\rho_{\rm w}=2\mathcal{A}_{d}\sigma H_{d}$ is released, where
$\mathcal{A}_{d}$ denotes the value of $\mathcal{A}$ at this time. The
resulting axion number density is similar that left from the evolution of
domain walls up to $H=H_{d}$, calculated in eq. (14). Therefore, the axion
relic abundance from energy stored in domain walls is approximately

$\displaystyle\frac{\Omega_{a}^{\rm w}}{\Omega_{\rm DM}}$
$\displaystyle\simeq\frac{\beta\mathcal{A}_{d}m_{a}f_{a}^{2}}{T_{\rm
eq}H_{d}^{1/2}M_{\rm Pl}^{3/2}}$ (15) $\displaystyle\simeq
0.1\mathcal{A}_{d}\left(\frac{m_{a}}{H_{d}}\right)^{1/2}\left(\frac{f_{a}}{10^{12}~{}{\,{\rm
GeV}}}\right)^{2}\left(\frac{m_{a}}{10^{-6}~{}{\,{\rm eV}}}\right)^{1/2}~{},$

where we set $\beta=10$ in the second line. For $H_{d}\ll m_{a}$ this
contribution far exceeds that produced by strings before $H=H_{\star}$ given
in eq. (10). Moreover, in Appendix C we show that the abundance produced from
energy in strings after $H=H_{\star}$ is in turn smaller than that released by
strings before $H=H_{\star}$ (fundamentally because $\rho_{s}$ and therefore
$\Gamma_{s}$ decrease faster than $\rho_{\rm w}$ and $\Gamma_{\rm w}$ as $t$
increases). Therefore we expect that eq. (15) approximates the total relic
abundance in this scenario.

If $H_{d}\ll m_{a}$ (i.e. $\epsilon$ is small enough), an axion in the post-
inflationary scenario with $N>1$ can account for the full dark matter
abundance with smaller $f_{a}$ than in the $N=1$ post-inflationary and the
pre-inflationary scenario for any value of $\theta_{0}$ compatible with
isocurvature constraints. The possible coupling of such an axion with
$E/N\simeq 1$ to photons is correspondingly enhanced, in particular

$\tilde{g}_{a\gamma\gamma}\simeq 5\cdot 10^{-16}{\,{\rm
GeV}}^{-1}\mathcal{A}_{d}^{1/2}\left(\frac{m_{a}}{H_{d}}\right)^{1/4}\left(\frac{m_{a}}{10^{-6}~{}{\,{\rm
eV}}}\right)^{1/4}~{}.$ (16)

As mentioned, we expect $H_{d}\sim\epsilon f_{a}^{2}/m_{a}$, but we keep eq.
(16) in terms of $H_{d}$ because we do not have a reliable way of predicting
the relation between $H_{d}$ and $\epsilon$.

###  2.3 The density power spectrum

The number of strings and domain walls per Hubble patch varies by order one
factors between different Hubble patches. Therefore, immediately after the
decay of the string-wall network (when most of the axions are expected to be
produced), the energy density of the axion field is expected to be left with
substantial inhomogeneities over Hubble distances. Being only in the axion
energy density, these are isocurvature fluctuations. Their properties are
conveniently described by the power spectrum $\mathcal{P}_{\delta}(k)$ of the
overdensity field
$\delta(\vec{x})\equiv(\rho_{a}(\vec{x})-\overline{\rho}_{a})/\overline{\rho}_{a}$,
where $\rho_{a}=\frac{1}{2}\dot{a}^{2}+\frac{1}{2}m_{a}^{2}a^{2}$ is the non-
relativistic energy density of the axion when the field is close to its
minimum, $\overline{\rho}_{a}$ is its spatial average, and we defined

$\langle\tilde{X}^{*}(t,\vec{k})\tilde{X}(t,\vec{k}^{\prime})\rangle\equiv(2\pi)^{3}\frac{2\pi^{2}}{k^{3}}\delta^{3}(\vec{k}-\vec{k}^{\prime})\mathcal{P}_{X}(t,k)~{},$
(17)

for a generic field $X$ with Fourier transform $\tilde{X}=\int
d^{3}x\exp(-i\vec{k}\cdot\vec{x})X(\vec{x})$. Unfortunately the power spectrum
of a post-inflationary ALP is not known; below we summarise its possible form
based on analytical arguments and hints from numerical simulations.

Figure 1:  The schematic expected form of the primordial power spectrum
$\mathcal{P}_{\delta}(k)$ of the axion over-density field
$\delta=(\rho_{a}-\overline{\rho}_{a})/\overline{\rho}_{a}$ after the string-
domain wall network decays for $N=1$ and $N>1$, where $k$ is the comoving
momentum. $\mathcal{P}_{\delta}(k)$ is thought to be peaked at a scale set by
the Hubble parameter when the network decays and to have a ‘white noise’
$k^{3}$ tail at smaller momenta. $H_{d}$ and $R_{d}$ are the Hubble parameter
and scale factor when the network is destroyed and $H_{\star}$ and $R_{\star}$
are those when $H=m_{a}$. Free-streaming could suppress the peak of
$\mathcal{P}_{\delta}$ in both scenarios if the dark matter axions are
produced mildly relativistic, and we plot $P_{\delta}$ dashed in this region
to emphasize the especially large uncertainty.

If $N=1$, the network decays at $H\simeq H_{\star}$ and therefore the
resulting fluctuations are likely to be on comoving scales approximately set
by $k_{\star}\equiv H_{\star}R_{\star}$. Numerical simulations at small scale
separation, $\log(m_{r}/H_{\star})\simeq 5\div 6$, support this picture, with
$\mathcal{P}_{\delta}(k)$ peaked (and obtaining an $\mathcal{O}(1)$ value) at
comoving momenta $k\simeq\mathcal{O}(10)k_{\star}$ [62, 63, 55]. (This
corresponds to $\delta$ of order one at spatial scales approximately set by
$\lambda_{\star}/\mathcal{O}(10)$, with $\lambda_{\star}=2\pi/k_{\star}$.) At
larger comoving scales, $k\ll k_{\star}$, the perturbations are uncorrelated,
by causality. From eq. (17), this fixes the spectrum to have the ‘white-noise’
form $\mathcal{P}_{\delta}(k)=\left(Ck/k_{\star}\right)^{3}$ for $k\lesssim
k_{\star}/C$, where $C$ is a numerical constant. The axion waves are expected
to be produced marginally relativistic (with momentum
$\mathcal{O}(10)k_{\star}$, i.e. wavelength coinciding with the typical size
of the strings and domain walls – which might also be connected to the density
spectrum being peaked at $\simeq 10k_{\star}$). Therefore they initially free
stream over comoving distances $\simeq(H_{\star}R_{\star})^{-1}$. This might
partly wash out the order-one fluctuations, but should leave the IR ($k^{3}$)
part of the spectrum largely unchanged because the axions will not free-stream
over the corresponding distances [63]. In Figure 1 we show a sketch of the
anticipated power spectrum.

If $N>1$ the spatial distribution of domain walls is expected to lead to order
one fluctuations in the axion energy density immediately after the network’s
destruction. Although not yet established with simulations, it is plausible
that these will be on scales set parametrically by $k_{d}\equiv H_{d}R_{d}$,
where $\mathcal{P}_{\delta}$ will be correspondingly peaked, i.e. on comoving
spatial scales that are larger than $k_{\star}^{-1}$. Similarly to the $N=1$
case, at $k\lesssim k_{d}$ the spectrum acquires a $k^{3}$ white-noise form

$\mathcal{P}_{\delta}(k)\simeq(Ck/k_{d})^{3}\text{ for $k\lesssim
k_{d}/C$}~{},$ (18)

where $C$ is an undetermined constant, possibly not far from order one. As
discussed in the previous Section, the axion momentum distribution is
uncertain in this scenario. If axions are dominantly produced with momentum
$H_{d}$, they are initially highly non-relativistic for $H_{d}\ll m_{a}$ and
their free-streaming will be irrelevant. However, if they are produced with
momentum $m_{a}$, they are initially semi-relativistic and will free-stream
over comoving distances of order $(H_{d}R_{d})^{-1}$.888More precisely, if
produced with momentum $m_{a}$ the axions would free-stream a comoving
distance $\simeq H_{d}^{-1}R_{d}^{-1}(1+\log(t_{\rm eq}/t_{d}))$ where $t_{\rm
eq}$ is the cosmic time at matter-radiation equality [64]. Similarly to the
$N=1$ scenario, fluctuations on comoving scales $k\gtrsim k_{d}$ will be
affected, but the $k^{3}$ part of $\mathcal{P}_{\delta}$ (which is what will
be relevant to the constraints that we discuss shortly) is expected to be left
unchanged. Unfortunately, the reach in $\log(m_{a}/H_{d})$ and
$\log(m_{r}/H_{\star})$ in simulations is not sufficient for us to extrapolate
to obtain the value of $C$ with the physical parameters and the details of the
peak of $\mathcal{P}_{\delta}$ are especially uncertain.

##  3 Comparison with numerical results

###  3.1 The area parameter

The area parameters $\mathcal{A}_{\star}$ and $\mathcal{A}_{d}$ entering in
the abundance of eq. (15) can be extracted from numerical simulations of eqs.
(6) and (7), although these have to be extrapolated in $\log(m_{r}/H)$ and,
for $\mathcal{A}_{d}$, also $\log(m_{a}/H)$. More precisely, we define
$\mathcal{A}$ to be the area of the surfaces where the axion is at the maxima
of the $\mathbb{Z}_{N}$-preserving potential in eq. (7), $a/f_{a}=2\pi n$,
$n=0,\dots,N-1$, which correspond to domain walls when $H\lesssim m_{a}$.
Details of the simulations and additional results may be found in Appendix C.

Figure 2: _Left:_ The time evolution of $\mathcal{A}$, which measures the area
of the surfaces on which $a/f_{a}=2\pi n$, $n=0,\dots,N-1,$ where $N$ is the
domain wall number, during the string scaling regime (i.e. with $m_{a}=0$) in
fat and physical theories. Time is represented by $\log(m_{r}/H)$, and we
factor out $N$ because $\mathcal{A}\propto N$. A clear increase of
$\mathcal{A}$ with $\log(m_{r}/H)$ is evident over a few e-foldings. We also
plot the best fits to the asymptotically-linear function in eq. (19). _Right:_
The evolution of $\mathcal{A}$ with time, parameterised by $\log(m_{a}/H)$,
for different values of the axion mass (i.e. different
$\log(m_{r}/m_{a})\equiv\log(m_{r}/H_{\star})$) and no explicit
$\mathbb{Z}_{N}$ breaking for the fat string system with $N=4$. For comparison
we also show results from simulations with $m_{a}=0$, dashed. As expected,
when $m_{a}\ll H$ the evolution matches that with $m_{a}=0$. At $H\lesssim
m_{a}$ the surfaces measured by $\mathcal{A}$ are domain walls. The increase
in $\mathcal{A}$ while $m_{a}\ll H$ results in a larger domain wall area at
$H\simeq m_{a}$ for larger $\log(m_{r}/H_{\star})$. At $H\lesssim m_{a}$ the
growth in $\mathcal{A}$ ends, although its asymptotic behaviour is challenging
to pin down.

Figure 2 (left) shows the evolution of $\mathcal{A}$ with $m_{a}=0$, i.e.
during the string scaling regime, with time parametrised by $\log(m_{r}/H)$.
In this era $\mathcal{A}$ is proportional to $N$, so we plot $\mathcal{A}/N$.
Results are given for both the so-called ‘fat’ string system (in which
$m_{r}\propto R^{-1}$ in eq. (6)) and the physical system (in which $m_{r}$ is
constant), with initial conditions close to the attractor solution.999The fat
string system displays better numerical properties. For instance, because more
cosmic time is available between two fixed values of $m_{r}/H$, this system
reaches the attractor solution faster. Not surprisingly given the logarithmic
increase of $\xi$ on the attractor, $\mathcal{A}$ grows with $\log(m_{r}/H)$.
For the fat string system $\mathcal{A}$ is approximately proportional to
$\xi$; this can be seen from the plot of $\mathcal{A}(t)/\xi(t)$ in Figure 5
in Appendix C.2. The results for $\mathcal{A}(t)/\xi(t)$ for the physical
system are less clear but are compatible with this ratio approaching a
constant value. This possibility is supported by the fact that
$\mathcal{A}\propto\xi$ is automatic if the average string curvature is
proportional to the Hubble parameter, as is expected by causality during the
scaling regime. On the other hand, even for the fat string system, the
uncertainty on $\mathcal{A}(t)/\xi(t)$ is still large enough that we cannot
exclude a deviation at much later times, $\log(m_{r}/H)\gg 1$.

Given that $\xi(t)$ during the scaling regime is well fitted by a linear
function of $\log\equiv\log(m_{r}/H)$ [56], a natural ansatz for $\mathcal{A}$
is

$\mathcal{A}(t)/N=c_{1}\log+c_{0}+c_{-1}\log^{-1}+c_{-2}\log^{-2}\,,$ (19)

where the last two terms encode possible deviations from linearity that are
irrelevant at late times, but affect the evolution at the small times (i.e.
$\log$s) that are accessible to simulations. This form fits the fat string
data well with $c_{1}=0.35(2)$. It also fits the physical data well, however
the value of $c_{1}$ is more uncertain and we estimate $c_{1}=0.36(10)$, see
Appendix C.2 (a quadratic dependence on $\log$ also gives an acceptable fit to
both $\xi$ and $\mathcal{A}$ for the physical system, see also [56]).
Therefore, at the physically relevant $\log(m_{r}/H_{\star})\simeq 40\div 80$
for an ALP (for $m_{a}$ in the range $10^{-15}~{}{\,{\rm eV}}\div{\,{\rm
eV}}$, c.f. Figure 4), the area parameter is likely to be
$\mathcal{A}_{\star}/N=\mathcal{O}(10\div 100)$. Meanwhile, given the relation
between $m_{a}$ and $f_{a}$, for the QCD axion $\log_{\star}\simeq 65$. As a
result, a linear extrapolation predicts $\mathcal{A}_{\star}\simeq 10$ in this
case (a quadratic extrapolation would instead give $\mathcal{A}_{\star}\simeq
100$). Similarly to $\xi$, these values of $\mathcal{A}_{\star}$ are much
larger than the $\mathcal{O}(1)$ values that occur at small scale separations.

In Figure 2 (right) we show the effect of the $\mathbb{Z}_{N}$-preserving
axion potential $V(a)$ in eq. (7) with $\delta V_{PQ}=0$. In particular, we
plot the evolution of $\mathcal{A}$ for the fat system for different choices
of the axion mass $m_{a}$, which we parametrise by
$\log(m_{r}/m_{a})\equiv\log(m_{r}/H_{\star})$, with time represented by
$\log(m_{a}/H)$. To allow for a larger dynamical range, $m_{a}$ is also taken
to be proportional to $1/R$, so that the ratio $m_{r}/m_{a}$ is constant.
(Although this is not physical, we expect the qualitative features of the
dynamics to be similar, like in the fat string system.) For comparison, we
also show the results obtained without any axion potential, $m_{a}=0$ (dashed
lines, which correspond to the evolution in Figure 2, left). We plot the data
starting from $\log(m_{r}/H)=3$, by which time strings are cleanly-defined
objects. Note that the simulated values of $\log(m_{r}/m_{a})$ are much
smaller than the physically relevant ones. The results indicate that at
$\log(m_{a}/H)\simeq 2$ the logarithmic growth of $\mathcal{A}$ stops. Not
surprisingly, the increase in the area of the surfaces where $a/f_{a}=2\pi n$
with $\log_{\star}$ translates into a larger number of domain walls at this
time (at least for the parameters possible in simulations) and the relation
$\mathcal{A}\propto N$ remains approximately true.101010As shown in Figure 3
of Appendix C.2, the growth of $\xi$ continues unaffected slightly longer.

Figure 2 (right) also shows a hint that $\mathcal{A}$ might decrease at large
$\log(m_{a}/H)$ in the simulations with larger $\log_{\star}$. Although
tentative, such a change is plausible and could indicate the domain walls,
rather than the strings, starting to dominate the dynamics of the system (as
might be expected once $\rho_{\rm w}\gtrsim\rho_{\rm s}$). Unfortunately, we
do not have sufficient numerical reach to determine the asymptotic behaviour,
e.g. it might be that the system is drawn to a new attractor with
$\mathcal{A}$ independent of $\log_{\star}$ but we cannot determine whether
this is the case. We also cannot rule out more dramatic qualitative changes in
the network’s dynamics at $\log_{\star}\gtrsim 5$ and $\log(m_{a}/H)\gg 1$.
Additionally, the non-linear transient of the waves produced during the
scaling regime (expected to happen when $\log_{\star}\gg 1$) could affect the
evolution of the domain walls around the time when $H=H_{\star}$ in a way that
we do not have control of in both the case of an ALP and especially the QCD
axion.

Given these multiple complicating factors, we have little control on the
physically relevant value of $\mathcal{A}_{d}$, and thus leave it unfixed in
what follows. The resulting uncertainty on $f_{a}$ and
$\tilde{g}_{a\gamma\gamma}$ is slightly tempered by the square-root in eq.
(16).

###  3.2 The emission spectrum

As discussed above eq. (13), the momentum of the axions emitted by the string-
domain wall network is important for the axion relic abundance. The momentum
distribution is encoded into the energy spectrum $\partial\rho_{a}/\partial
k$, which can be extracted from numerical simulations as described in Appendix
C. However, while domain walls are present they contaminate the axion field
(and because they have thickness $\simeq m_{a}^{-1}$ they cannot be easily
masked). Instead, it is more convenient to carry out simulations with non-zero
$\mathbb{Z}_{N}$ breaking parameter, and extract the axion energy spectrum
after all the domain walls are destroyed. However, we stress that we can only
determine the evolution and the spectrum at very unphysical values of
$m_{a}/m_{r}$ and $H_{d}/m_{r}$, and the true emission spectrum could differ
from simulation results dramatically.

Figure 3:  _Left:_ The domain wall area parameter $\mathcal{A}$ in fat string
and domain wall theories with $N=4$ and $\log(m_{r}/m_{a})=2.7$ for different
values of the $\mathbb{Z}_{N}$ symmetry breaking parameter
$\tilde{\epsilon}=v\sqrt{\epsilon}/(Nm_{a})$. _Right:_ The axion energy
density spectra $\partial\rho_{a}/\partial\log k$ at the final simulation
time, after the network is destroyed, for the simulations plotted in the left
panel. The $x$-axis is normalised relative to the comoving momentum
$k_{m}\equiv Rm_{a}$ corresponding to the axion mass $m_{a}$. We also indicate
the comoving momentum corresponding to the string core scale $k_{r}/2$. The
coloured ticks indicate $2\pi k_{d}$, where $k_{d}$ is the comoving momentum
that is equal to $H$ at the time when the network is destroyed. For
sufficiently small $\tilde{\epsilon}$ (i.e. a relatively long lived network),
the axions seem to be emitted by domain walls at most mildly relativistic,
with a spectrum peaked at momenta of the order of the axion mass $k\simeq
k_{m}$, but it is challenging to extrapolate this result at the physical
parameter space.

In Figure 3 (left) we plot $\mathcal{A}$ as a function of time in such
simulations for the fat string-domain wall system with different amounts of
$\mathbb{Z}_{N}$ breaking, which we parameterise by
$\tilde{\epsilon}=v\sqrt{\epsilon}/(Nm_{a})$ where $\epsilon$ is defined in
eq. (8). We fix $\log(m_{r}/H_{\star})=2.7$ (i.e. $m_{r}/m_{a}=15$) and $N=4$.
As expected, larger $\tilde{\epsilon}$ results in the network being destroyed
earlier.111111We do not attempt to fit a relation between the
$\tilde{\epsilon}$ and the Hubble parameter when the network is destroyed
$H_{d}$, which would require extremely careful extrapolations in $m_{r}/m_{a}$
and $m_{a}/H_{d}$. In Figure 3 (right) we plot the axion energy spectrum at
the end of the simulations shown in the left panel, by which time all the
strings and domain walls are destroyed, the axion field amplitude is $\ll
f_{a}$ over all of space (apart from a few oscillons [62], which contain a
negligible fraction of the axion energy density), and the spectrum has reached
a fixed shape. Because of the fat strings and domain walls, axions produced
with proper momentum $m_{a}$ or $m_{r}/2$ (the latter corresponding to the
string core scale) always have, and remain at, the same comoving momentum
$k_{m}=m_{a}(t)R(t)$ or $k_{r}/2=m_{r}(t)R(t)/2$ respectively. At the values
of $\log(m_{r}/H)$ that can be simulated, the strings are expected to emit
energy dominantly to axions with UV momentum $m_{r}$ [53, 56] although this
has only been shown with $m_{a}=0$.

In Figure 3 (right), the axion energy spectrum increases as
$\tilde{\epsilon}$, and therefore $H_{d}$, decreases. This indicates that, at
least for the smallest $\tilde{\epsilon}$, a sizable fraction of the energy
density in axion waves is coming from the domain walls (because the relic
abundance from strings is dominantly produced at $H=H_{\star}$ rather than
later). Although $m_{r}/m_{a}=15$ only gives a mild separation between $k_{a}$
and $k_{r}/2$, Figure 3 suggests that with longer lived domain walls an
increased fraction of the axion waves are produced at most borderline rather
that ultra-relativistic, with physical momenta $\simeq m_{a}$ rather than
$m_{r}/2$ (or equivalently, comoving momentum $\simeq k_{m}$ rather than
$k_{r}/2$). Although at these scale separations the emission spectrum from
domain walls is peaked at $k_{m}$ rather than $k_{d}$ (the latter of which is
set by the Hubble parameter when the networks is destroyed),121212We have
checked that results from simulations with different $N$ and $m_{r}/m_{a}$ are
consistent. we cannot exclude that this changes at much larger values of
$\log(m_{r}/m_{a})$ and $\log(m_{a}/H)$, and we leave a full extrapolation for
future work.

##  4 Observational implications

###  4.1 Constraints

There are numerous constraints on $g_{a\gamma\gamma}$ that apply to either any
axion or any axion that is DM. Additionally there are bounds that are specific
to post-inflationary axions.

Given that their energy redshifts differently to dark matter (in particular,
slower), the domain walls must decay before matter-radiation equality (MRE) if
the axion makes up a sizable component of dark matter. In the $N>1$ case, this
requires that the Hubble parameter when the network is destroyed $H_{d}$ is
larger than $H_{\rm eq}$, where ‘eq’ denotes quantities at MRE (similarly, in
$N=1$ theories $H_{\star}>H_{\rm eq}$ is needed). If the axion makes the whole
DM, using eq. (15), this constrains

$f_{a}\gtrsim 10^{7}{\,{\rm
GeV}}\frac{1}{\mathcal{A}_{d}^{1/2}}\left(\frac{10^{-6}~{}{\,{\rm
eV}}}{m_{a}}\right)^{1/2}~{}.$ (20)

Observational limits on the fraction of primordial isocurvature density
perturbations lead to less certain, but probably stronger, constraints on both
$N=1$ [65, 66, 67, 55] and $N>1$ post-inflationary theories. These require any
primordial isocurvature fluctuations (such as the axion’s ones), with power
spectrum $\mathcal{P}_{\rm iso}(k)$, to be suppressed on observable scales
with respect to the initially small adiabatic fluctuations, described by the
almost scale-invariant power spectrum of curvature perturbations,
$\mathcal{P}_{\mathcal{R}}(k)=A_{s}(k/k_{p})^{n_{s}-1}$, where
$k_{p}/R_{0}=0.05\,{\rm Mpc}^{-1}$ is a pivot scale and $A_{s}=2.2\times
10^{-9}$, with $n_{s}\simeq 1$. In particular, at masses $m\gtrsim
10^{-22}\,{\,{\rm eV}}$, the strongest such bound comes from Lyman-$\alpha$
observations (looking at momenta $k/R_{0}=k_{L\alpha}/R_{0}\equiv 10\,{\rm
Mpc}^{-1}$) and limits the fraction of isocurvature perturbations as $f_{\rm
iso}^{2}(k)\equiv\mathcal{P}_{\rm iso}(k)/\mathcal{P}_{\mathcal{R}}(k)\lesssim
0.004$ at $k=k_{p}$ assuming $\mathcal{P}_{\rm iso}(k)\propto k^{3}$ for
$k\in(k_{p},k_{L\alpha})$ [66].

To have a chance of being compatible with isocurvature constraints, the axion
modes on this scale must be in the $k^{3}$ part of the power spectrum, so the
resulting constraint should only depend on the numerical coefficient $C$ and
$m_{a}$. For $N=1$, using eq. (18), $f^{2}_{\rm iso}=(\Omega_{a}/\Omega_{\rm
DM})^{2}(Ck_{p}/k_{\star})^{3}/A_{s}$ and therefore this isocurvature bound
corresponds to

$m_{a}\gtrsim 2\cdot 10^{-17}C^{2}\left(\frac{\Omega_{a}}{\Omega_{\rm
DM}}\right)^{4/3}\,{\,{\rm eV}}~{}.$ (21)

where, based on numerical simulations at small scale separation, a value of
order $C=\mathcal{O}(10^{-1})$ is plausible [55].

Also for $N>1$ the constrained scales must lie in the $k^{3}$ part of
$\mathcal{P}_{\delta}$. If the network decays at $H_{d}\ll m_{a}$ and turns
out to indeed have a peak at $k\simeq k_{d}$, the resulting bound is stronger
than in the $N=1$ case because the order one perturbations are on much larger
spatial scales $k_{d}\ll k_{\star}$, and $\mathcal{P}_{\delta}$ at observable
scales is less suppressed by the $k^{3}$ form. Under this reasonable
assumption, $H_{d}$ plays exactly the role that $m_{a}$ does in the $N=1$
case, so the constraint in [66] can be simply translated to

$H_{d}\gtrsim 2\cdot 10^{-17}C^{2}\left(\frac{\Omega_{a}}{\Omega_{\rm
DM}}\right)^{4/3}\,{\,{\rm eV}}~{}.$ (22)

Combining eq. (22) with the result for the relic abundance eq. (15) gives the
bound from isocurvature perturbations

$f_{a}\gtrsim 6\times 10^{9}~{}{\,{\rm
GeV}}\frac{1}{\mathcal{A}_{d}^{1/2}}\left(\frac{10^{-6}~{}{\,{\rm
eV}}}{m_{a}}\right)^{1/2}C^{1/2}\left(\frac{\Omega_{a}}{\Omega_{\rm
DM}}\right)^{5/6}~{},$ (23)

and a corresponding upper bound on $\tilde{g}_{a\gamma\gamma}$. As discussed,
we cannot reliably extract the value of $C$ in the $N>1$ scenario from
numerical simulations. However, based on results for the $N=1$ string-wall
network, a reasonable expectation is $C=\mathcal{O}(10^{-1}\div 1)$, but not
too much smaller [55]. Note that the bound on $f_{a}$ (and therefore
$\tilde{g}_{a\gamma\gamma}$) is weaker for larger $\mathcal{A}_{d}$ because
this provides the correct dark matter abundance for larger $H_{d}$.

If the axions are produced with initial momentum of order $m_{a}$, there is a
constraint from the axions behaving as warm dark matter because the free-
streaming of the dark matter particles would suppress the formation of small-
scale structures [68, 69, 32]. In theories of thermally produced warm dark
matter, the observational bound on the dark matter mass is $m_{\rm WDM}\gtrsim
10{\,{\rm keV}}$, i.e. the dark matter must be at least borderline non-
relativistic when the Universe has temperature $T=10{\,{\rm keV}}$ [70, 71,
72]. The free streaming length, and therefore the resulting cut-off in the
linear matter power spectrum, only depends on the dark matter equation of
state [73]. Therefore, a limit on $N>1$ theories can be obtained by demanding
that, if the axions are produced semi-relativistic, $H_{d}\gtrsim
H_{T=10{\,{\rm keV}}}$ where $H_{T=10{\,{\rm keV}}}$ is the value of the
Hubble parameter when $T=10{\,{\rm keV}}$. This leads to the constraint on an
axion that makes up all the dark matter of 131313There is uncertainty on this
bound both from the unknown details of the initial power spectrum, e.g. how
relativistic the axions are at $H_{d}$ and also from the fact we did not do a
full analysis of the effects of warm, low mass, axions and instead
reinterpreted existing bounds. With a reliable prediction of the energy
spectrum of the axions produced by the domain walls, one could obtain an
accurate constraint by following the approach of [73].

$f_{a}\gtrsim 10^{9}~{}{\,{\rm
GeV}}\frac{1}{\mathcal{A}_{d}^{1/2}}\left(\frac{10^{-6}~{}{\,{\rm
eV}}}{m_{a}}\right)^{1/2}~{},$ (24)

and a corresponding one on $\tilde{g}_{a\gamma\gamma}$. The resulting bound is
approximately an order of magnitude weaker than the isocurvature one with
$C\simeq\mathcal{O}(0.1)$. If $\Omega_{a}<\Omega_{\rm DM}$ the warm dark
matter bound is less certain. Fitting the numerical results of [74] suggests
that, approximately, the minimum temperature at which the axions must be non-
relativistic (which coincides with the destruction temperature $T_{d}$ with
our present assumptions) is $T_{d}\propto\left(\Omega_{a}/\Omega_{\rm
DM}\right)^{1.5}$ in which case the allowed $f_{a}$ in eq. (24) would scale as
$\propto\left(\Omega_{a}/\Omega_{\rm DM}\right)^{1.25}$. (Note that, as shown
in Ref. [75], the combination of the isocurvature bound and the constraint on
warm dark matter leads to the limit $m_{a}\gtrsim 10^{-18}~{}{\,{\rm eV}}$ in
post-inflationary theories regardless of the scale $k$ at which
$\mathcal{P}_{\delta}$ is of order one.)

We note that the explicit $\mathbb{Z}_{N}$ symmetry breaking parameter
$\epsilon$ enters in all of the constraints only via $H_{d}$, so the
uncertainty on $H_{d}$ as a function of the $\epsilon$ does not propagate to
eqs. (20), (23), and (24).

Figure 4:  _Blue colours:_ The axion-to-photon coupling $g_{a\gamma\gamma}$ as
a function of the axion mass $m_{a}$ in different minimal scenarios, requiring
that an ALP makes up the full DM relic abundance. We assume that this coupling
is given in terms of $f_{a}$ as $g_{a\gamma\gamma}\equiv\alpha_{em}/(2\pi
f_{a})$, i.e. $E/N=1$ in eq. (3); in a theory with $E/N\neq 1$ the
corresponding $g_{a\gamma\gamma}$ will be shifted. Dashed lines represent the
pre-inflationary scenario for different misalignment angles, $\theta_{0}$.
$\theta_{\rm max}$ is the largest such angle compatible with isocurvature
constraints if the Hubble scale during inflation $H_{I}=1{\,{\rm eV}}$
(smaller $\theta_{0}$ are required if $H_{I}$ is larger). The blue shaded
region represents the post-inflationary scenario with $N>1$ as the parameter
$\epsilon$, which controls the $\mathbb{Z}_{N}$-breaking potential, varies.
The area labelled ‘excluded by isocurvature’ is likely to be in tension with
isocurvature constraints, although these depend on the axion’s poorly known
density power spectrum after the network decays; the results are plotted
assuming the coefficient of the white-noise tail in eq. (18) to be $C=0.1$.
The region ‘$H_{d}<H_{\rm eq}$’ is excluded because the domain walls would be
destroyed after matter-radiation equality. Both these constraints depend on
the domain wall area parameter at the time of the network’s decay
$\mathcal{A}_{d}$, which is uncertain (we plot the results for
$\mathcal{A}_{d}=1$ and $10$, dashed and solid). The lower boundary of the
$N>1$ region very approximately coincides with the $N=1$ scenario. _Red/orange
colours:_ Existing experimental and observational bounds and the reach of some
proposed future haloscopes (orange) [76].

###  4.2 Relevance for experimental searches

The ALP-to-photon couplings in the pre- and post-inflationary scenarios are
shown in Figure 4 imposing that $\Omega_{a}=\Omega_{\rm DM}$. To enable
comparison with existing and projected constraints we assume $E/N=1$, so
$g_{a\gamma\gamma}=\tilde{g}_{a\gamma\gamma}\equiv\alpha_{em}/(2\pi
f_{a})$.141414The current constraints on $g_{a\gamma\gamma}$ in Figure 4 come
from various sources [77, 78, 79, 80, 81, 82, 83, 84, 85] (we do not attempt
to be exhaustive and simply present the most important ones in each mass
range). We also plot the bounds from black hole super-radiance, which do not
require an axion-photon coupling [86, 87] and instead are relevant only for
$f_{a}$ not too small, because a too large axion self-interaction hinders the
super-radiance process. The axion-photon coupling for the QCD axion, given in
eq. (2), is also shown; we set the upper and lower edges of the band by
$E/N=44/3$ and $2$ respectively, but we stress that couplings outside the
range plotted are possible (the band is grayed for $m_{a}\gtrsim
10^{-2}~{}{\,{\rm eV}}$ because astrophysical observations rule out many such
QCD axion theories [88, 89, 90, 91, 92]).

In Figure 4, dashed lines show the values of the initial misalignment angle
$\theta_{0}$ for which an ALP in the pre-inflationary scenario makes up the
full relic density. As can be seen, such an axion with $E/N\simeq 1$ can be
all DM and compatible with existing constraints for almost any
$m_{a}\in(10^{-22}~{}{\,{\rm eV}},10^{2}~{}{\,{\rm eV}})$ [93, 94, 95, 96,
82]. Meanwhile, an ALP in the post-inflationary $N>1$ scenario can make up the
full dark matter abundance over the dark blue area as the value of the
explicit symmetry breaking parameter $\epsilon$, which determines $H_{d}$,
varies. Points with smaller $f_{a}$ (i.e. larger $\tilde{g}_{a\gamma\gamma}$)
correspond to smaller $\epsilon$, i.e. smaller $H_{d}$. The upper edge of this
region is set by the lowest $H_{d}$ allowed by isocurvature constraints, which
depend on the uncertain parameters $\mathcal{A}_{d}$ and $C$ (we also plot the
weaker bound from requiring the domain wall network is destroyed before MRE).
The lower edge of the $N>1$ post-inflationary region corresponds to
$\epsilon\simeq m_{a}^{2}/f_{a}^{2}$ so $H_{d}\simeq m_{a}$ and the domain
walls are destroyed soon after forming. In this case the DM abundance, and
therefore $\tilde{g}_{a\gamma\gamma}$, is of the same order as that in the
$N=1$ post-inflationary scenario, which we calculate considering only the
abundance produced by strings eq. (10) (if domain walls give a large
contribution to the relic abundance the prediction would be at larger
$g_{a\gamma\gamma}$).

###  4.3 Possible implications for gravitational wave searches

In a post-inflationary $N>1$ axion theory the evolution of the domain wall
network between $H_{\star}$ and $H_{d}$ produces a stochastic GW background
[97, 61, 31, 32]. Due to the uncertainties in the dynamics of the network, the
rate of production of these GWs is unknown. However, standard assumptions
making use of the quadrupole approximation predict that the domain walls emit
energy density in GWs at a rate $\Gamma_{\rm GW}$ roughly given by (see e.g.
[98, 61, 31])

$\displaystyle\Gamma_{\rm GW}/H$
$\displaystyle\simeq\mathcal{A}^{2}G\sigma^{2}$ (25)
$\displaystyle\simeq\mathcal{A}^{2}Gm_{a}^{2}f_{a}^{4}~{}.$

If eq. (25) is accurate, the energy lost into GWs is much smaller than the
total energy the network emits (in eq. (12)) over all of the allowed parameter
space in Figure 4, and their backreaction is therefore not expected to affect
the dynamics of the domain walls.

Because the domain walls have typical curvature of order Hubble, the
instantaneous emission spectrum of the GWs is expected to be peaked at the
Hubble parameter at the time of production. This is indeed the case in
numerical simulations at small scale separations [61]. If true with the
physical parameters, the GW spectrum today $\Omega_{\rm
GW}(f)\equiv\frac{1}{\rho_{\rm c}}\partial\rho_{\rm GW}/\partial\log f$, where
$f$ is the present-day GW frequency, will be peaked at the frequency
corresponding to the Hubble parameter at the time the network is destroyed
$f^{\rm peak}\simeq 10^{-16}~{}{\rm Hz}\left(H_{d}/H_{\rm eq}\right)^{1/2}$,
i.e.

$f^{\rm peak}\simeq 7\cdot 10^{-7}~{}{\rm
Hz}~{}\mathcal{A}_{d}\left(\frac{\Omega_{a}}{\Omega_{\rm
DM}}\right)^{-1}\left(\frac{f_{a}}{10^{12}~{}{\,{\rm
GeV}}}\right)^{2}\left(\frac{m_{a}}{10^{-6}~{}{\,{\rm eV}}}\right)~{}.$ (26)

Assuming that eq. (25) reproduces the actual emission, the amplitude of the GW
spectrum at its peak is

$\displaystyle\Omega_{\rm GW}|_{\rm peak}$ $\displaystyle\simeq
10^{-31}\mathcal{A}_{d}^{2}\left(\frac{f_{a}}{10^{12}~{}{\,{\rm
GeV}}}\right)^{4}\left(\frac{m_{a}}{H_{d}}\right)^{2}$ (27)
$\displaystyle\simeq 5\cdot
10^{-28}\frac{1}{\mathcal{A}_{d}^{2}}\left(\frac{\Omega_{a}}{\Omega_{\rm
DM}}\right)^{4}\left(\frac{10^{12}~{}{\,{\rm
GeV}}}{f_{a}}\right)^{4}\left(\frac{10^{-6}~{}{\,{\rm
eV}}}{m_{a}}\right)^{2}.$

Eq. (26) and the second line of eq. (27) follow from the expression for the
axion relic abundance in eq. (15) and we allow $\Omega_{a}\neq\Omega_{\rm DM}$
for generality. As is evident from eq. (27), the strongest GW signals occur
for the smallest allowed $f_{a}$. This is because, for fixed $m_{a}$ and
$\Omega_{a}$, smaller $f_{a}$ corresponds to smaller $H_{d}$, so the GWs are
redshifted less.

The possibility of observing these GWs is in tension with the observational
bound on isocurvature perturbations described in Section 4.1: combing eq. (27)
with eq. (23), the isocurvature bound implies

$\Omega_{\rm GW}|_{\rm peak}\lesssim 4\cdot
10^{-19}\frac{1}{C^{2}}\left(\frac{\Omega_{a}}{\Omega_{\rm
DM}}\right)^{2/3}~{},$ (28)

independent of $m_{a}$ and $\mathcal{A}_{d}$. For
$C\gtrsim\mathcal{O}(10^{-2})$ we have $\Omega_{\rm GW}|_{\rm peak}\lesssim
10^{-15}$, and if the axion is only a small fraction of DM, this limit is not
relaxed.151515Meanwhile, the warm dark matter bound eq. (24), which applies if
the axions are produced mildly relativistic, leads to $\Omega_{\rm GW}|_{\rm
peak}\lesssim 5\cdot 10^{-16}\mathcal{A}_{d}^{-1}\left(\Omega_{a}/\Omega_{\rm
DM}\right)^{-1}$. This is weaker than eq. (28) for $C\gtrsim 10^{-2}$ (see
also [32] for a related analysis). Near-future detectors are likely to only
reach sensitivities $\Omega_{\rm GW}\simeq 10^{-15}$ for any frequency
predicted by eq. (26).161616Pulsar timing arrays will, at best, reach
$\Omega_{\rm GW}\simeq 10^{-15}$ for frequencies around $\simeq 10^{-8}~{}{\rm
Hz}$, and space based interferometers could have similar reach at higher
frequencies. Because GWs redshift like radiation, the energy density they
contain can be constrained by CMB observations [99], but current bounds are in
the region of $\Omega_{\rm GW}\simeq 10^{-6}$ [100] and could be improved by
only a few orders of magnitude (e.g. by the future EUCLID satellite). Finally,
at frequencies $\lesssim 10^{-14}~{}{\rm Hz}$ measurements of the CMB
polarisation data are relevant (see e.g. [101]). Current limits are
$\Omega_{\rm GW}\simeq 10^{-12}$ for GWs of frequency $f_{\rm eq}\simeq
10^{-16}~{}\rm Hz$ [102, 103] (corresponding to $H_{d}=H_{\rm eq}$) and at
higher frequencies these get weaker $\propto(f_{\rm eq}/f)^{2}$ (planned
experiments e.g. LiteBIRD [104], could have a factor $\sim 100$ better
sensitivity). Therefore, unless the unknown constant $C$ (which parametrises
the white-noise tail of the axion density power spectrum, see eq. (18)) is
much smaller than the naive expectation, the allowed GW signals from ALP
domain walls are barely detectable regardless of $m_{a}$ and $f_{a}$. We do
however stress that eq. (25) is just a rough parametric expectation, and an
analysis of the GW emission from domain walls in numerical simulations taking
into account the several required extrapolations would certainly be
worthwhile, and left for future work. We also note that it is possible that
the network stops emitting energy to axions efficiently enough to maintain an
approximately constant $\mathcal{A}$ at scale separations beyond the reach of
simulations, completely altering the network’s dynamics and the resulting GW
prediction.

GWs are also produced by the axion strings prior to $H_{\star}$ [55] and the
resulting stochastic background is detectable, compatibly with constraints
from isocurvature perturbations, dark radiation, and axion dark matter
overproduction, if $f_{a}\gtrsim 10^{14}~{}{\,{\rm GeV}}$, and
$10^{-20}~{}{\,{\rm eV}}<m_{a}<10^{-18}~{}{\,{\rm eV}}$. For such high $f_{a}$
however typical $g_{a\gamma\gamma}$ are suppressed and the relic axions cannot
be directly detected (except for specific models e.g. based on [105, 106],
where $g_{a\gamma\gamma}\ggg 1/f_{a}$).

###  4.4 Comments on ALP dark matter substructure

As discussed in Section 2.3, the decay of the string-domain wall network leads
to order one fluctuations in the axion dark matter overdensity field at small
spatial scales, which might lead to interesting dark matter substructure [107,
108, 109, 110, 63], see [111, 112, 113, 114, 115, 116, 117] for analysis in
the case of the QCD axion. Although their properties are extremely uncertain,
we now briefly comment on the evolution of an ALP’s initial fluctuations
assuming that there are order-one overdensities on scales $H_{d}^{-1}$. As we
will see, quantum pressure and axion self-interactions can play an important,
previously neglected, role.

After the string-wall decays, the axion follows the nonrelativistic
Schrödinger–Poisson EoM, which can be rewritten as the continuity and Euler
equation of the density and velocity fields, $\rho$ and $v$. At the small
spatial scales we are interested in,
$\delta\equiv(\rho_{a}-\overline{\rho}_{a})/\overline{\rho}_{a}=\mathcal{O}(1)$.
Nevertheless, we can obtain a qualitative understanding of the evolution of a
generic mode (with comoving momentum $k$) by considering the linearised form
of the EoM

$\ddot{\delta}+2H\dot{\delta}-\left[4\pi
G\rho_{a}-\frac{k^{2}}{R^{2}}\left(\frac{k^{2}}{4m^{2}R^{2}}+\frac{g\rho_{a}}{m_{a}^{2}}\right)\right]\delta=0\,,$
(29)

where $g\equiv\lambda/m_{a}^{2}\simeq-1/(8f_{a}^{2})$ parametrises the axion
quartic coupling $V\supset\lambda a^{4}/3$, which we assume is attractive as
is the case for a cosine potential (higher order couplings are present but
negligible). The square bracket in eq. (29) is the Laplacian of the total
potential $\Phi+\Phi_{Q}+\Phi_{\rm self}$ acting on the field, where
$\Phi_{Q}\equiv-\left(\nabla^{2}\sqrt{\rho_{a}}\right)/(2m_{a}^{2}R^{2}\sqrt{\rho_{a}})$
and $\Phi_{\rm self}\equiv(g\rho_{a})/m_{a}^{2}$ are the ‘quantum’ pressure
and self-interaction potentials respectively, while $\Phi$ is the
gravitational one, satisfying $\nabla^{2}\Phi=4\pi GR^{2}\rho_{a}$. Both
$\Phi$ and $\Phi_{\rm self}$ are attractive and tend to increase $\delta$
(given that $g<0$), while quantum pressure is repulsive and makes $\delta$
oscillate, preventing perturbations from collapsing. Note that eq. (29)
reduces to the usual equation describing cold dark matter perturbations (the
Meszaros equation) in the limits $m\rightarrow\infty$ and $g\rightarrow 0$. We
refer to [118, 119, 120, 121, 122] for the derivation of eq. (29).

The EoM (29) takes a simplified form when expressed in terms of $y\equiv
R/R_{\rm eq}$,

$\delta^{\prime\prime}+\frac{3y+2}{2y(y+1)}\delta^{\prime}-\frac{3}{2y(y+1)}\left[1-\left(\frac{k}{k_{J}^{\rm
eq}}\right)^{4}\frac{1}{y}+\left(\frac{k}{k_{\rm self}^{\rm
eq}}\right)^{2}\frac{1}{y^{2}}\right]\delta=0\,,$ (30)

where $\delta^{\prime}\equiv d\delta/dy$ and $k_{J}$ and $k_{\rm self}$ are
two comoving momentum scales associated with quantum pressure (‘quantum’ Jeans
scale [123]) and self-interactions, defined by

$k_{J}\equiv R(16\pi G\rho_{a}m_{a}^{2})^{1/4}~{},\qquad\qquad k_{\rm
self}\equiv Rm_{a}(4\pi G/|g|)^{1/2}~{},$ (31)

respectively, and the superscript ‘eq’ refers to quantities calculated at
$R=R_{\rm eq}$. For very IR modes ($k/k_{J}^{\rm eq}\ll 1$ and $k/k_{\rm
self}^{\rm eq}\ll 1$, which includes the normal adiabatic modes at
cosmological scales) quantum pressure and self-interactions are negligible,
and eq. (30) shows that such perturbations increase linearly around MRE as
$\delta\propto 1+(3/2)R/R_{\rm eq}$ while order-one fluctuations collapse into
bound objects at $R=R_{\rm eq}/\delta$. However, for large enough $k$ quantum
pressure and self-interactions affect the evolution by preventing the modes
from collapsing, or enhancing their growth.

Interestingly, at MRE the quantum Jeans scale $k_{J}^{\rm eq}$ coincides with
$k_{\star}$, modulo an order one factor:

$\frac{k_{\star}}{k_{J}^{\rm eq}}\simeq 0.53~{},$ (32)

similarly to vector bosons produced by inflationary fluctuations [124]. This
relation holds regardless of the value of $m_{a}$ and even Newton’s constant
(but it is modified if the axion mass is temperature dependent).171717In eq.
(32) we assume that the number of relativistic degrees of freedom $g_{\star}$
at $H_{\star}$ is given by the standard model high temperature value, which is
appropriate for $m_{a}\gtrsim 10^{-5}{\,{\rm eV}}$. However, the numerical
factor in eq. (32) only depends on $g_{\star}^{1/12}$, so this has a
negligible impact on our conclusions (e.g. of order 100 beyond Standard Model
degrees of freedom around the TeV scale would not affect our conclusions).

#### $\boldsymbol{N=1}$

In theories with $N=1$ the density power spectrum $\mathcal{P}_{\delta}$ is
expected to be peaked at $k_{\rm peak}=xk_{\star}$ with
$\mathcal{P}_{\delta}(k_{\rm peak})\sim\mathcal{O}(1)$, where $x$ is a
numerical factor that is highly uncertain (simulations at
$\log(m_{r}/H_{\star})\simeq 5$ suggest $x\sim\mathcal{O}(10)$ [55] although
this could be modified prior to MRE by free-streaming or, as we discuss
shortly, self-interactions). Therefore, given eq. (32), the collapse of order-
one fluctuations will be affected by quantum pressure unless $x\ll 1$, which
would be surprising. In particular if $x\gtrsim 2$ quantum pressure will
prevent fluctuations collapsing at MRE, and compact objects are not expected
to form from these fluctuations.

Meanwhile, for all the values of $m_{a}$ and $f_{a}$ of interest
$k_{\star}/k_{\rm self}^{\rm eq}\ll 1$, so at MRE the self-interactions do not
affect the modes at $k\simeq k_{\star}$ or smaller. However, self-interactions
could be important deep in radiation domination,181818We thank Asimina
Arvanitaki for pointing this out. when the corresponding terms in eqs. (29)
and eq. (30) can be larger than the quantum pressure term because $\Phi_{g}$
is proportional to the density. By solving eq. (30) starting at
$R_{0}=R_{\star}$ with $\delta^{\prime}(R_{0})=0$, one can see that indeed
perturbative modes at momenta $k\gtrsim 10^{-2}k_{\star}$ grow substantially
during radiation domination due to the self-interactions. This could
significantly modify the initial spectrum of density fluctuations before the
overdensities collapse. Given the large uncertainties on the initial spectrum,
we leave a full study for future work.

#### $\boldsymbol{N>1}$

Conversely, for $N>1$ the string network decays at a time when
$H_{d}<H_{\star}$ and the order one fluctuations are expected to lie on larger
comoving scales, $xk_{d}$. The numerical value of $x$, although unknown, might
be $\mathcal{O}(10)$, similarly to the $N=1$ scenario. At MRE, the ratio
between the size of these fluctuations and the quantum Jeans scale reads

$\displaystyle\frac{k_{\rm peak}}{k_{J}^{\rm eq}}$ $\displaystyle\simeq
x\left(\frac{H_{d}}{m_{a}}\right)^{1/2}~{}.$ (33)

Thus, for $H_{d}\ll m_{a}$ (i.e. well within the blue region in the parameter
space in Figure 4), quantum pressure is not expected to be relevant in the
dynamics around MRE. Additionally, by solving eq. (30) starting at
$R_{0}=R_{d}$, self-interactions are also negligible throughout for modes
$k\simeq k_{d}$ or smaller. Thus, the order-one perturbations are expected to
collapse into compact objects (miniclusters) with density
$\sim\overline{\rho}(R_{\rm eq})$, i.e. set by the average dark matter density
when they collapse. The mass of the miniclusters will be similar to the mass
contained within an initially order one fluctuation

$\displaystyle M_{d}$
$\displaystyle\simeq\frac{4\pi}{3}\left(\frac{2\pi}{k_{d}}\right)^{3}R_{\rm
eq}^{3}\overline{\rho}(R_{\rm eq})x^{-3}$ (34) $\displaystyle\simeq 3.5\cdot
10^{-10}M_{\odot}A_{d}^{-3}\left(\frac{10^{12}~{}{\,{\rm
GeV}}}{f_{a}}\right)^{6}\left(\frac{10^{-6}~{}{\,{\rm
eV}}}{m_{a}}\right)^{3}x^{-3}~{}.$

If $x\sim\mathcal{O}(10)$, $M_{d}$ varies from roughly $10^{-22}M_{\odot}$ to
$10^{2}M_{\odot}$ over the allowed parameter space in Figure 4. Larger mass
substructure will form at larger scales from the $k^{3}$ part of the initial
density power spectrum, as these modes grow and collapse during matter
domination, however the resulting clumps have lower density and are more
likely to be subsequently destroyed by e.g. tidal disruption [124].

##  5 Discussion and future directions

We have shown that it is plausible that a post-inflationary axion-like-
particle with $N>1$ can comprise the full dark matter abundance with values of
$g_{a\gamma\gamma}$ that are fairly easily discoverable by several proposed
haloscopes and are also relevant for searches using astrophysical
observations, see Figure 4.191919Among others, a post-inflationary $N>1$ ALP
could be discovered by the experiments TOORAD [125], BREAD [126], BRASS [127],
ORGAN [128], MADMAX [129], ALPHA [130], the next stages of ADMX, KLASH [131],
DM-Radio [132], WISPLC [133], the broadband mode of ABRACADABRA [134], and
searches with polaritons [135] (projections are not plotted for all of these).
Meanwhile, parts of the parameter space are constrained by observations by the
telescopes VIMOS [85] and MUSE [84] and extra-galactic background light (EBL),
and the parameter space is close to bounds from measurements of the
polarisation of photons emitted from magnetic white dwarf stars [80] (‘MWD’).
Even larger $\tilde{g}_{a\gamma\gamma}$ are possible if the axion makes up
only a fraction of the dark matter, but the discoverability in this case
depends on the details of a particular experiment. Meanwhile, the DM
substructure that is likely to form in $N>1$ theories with $H_{d}\ll m_{a}$
could provide complementary observational prospects. For instance, the
miniclusters could be detected through gravitational interactions (e.g. by
searches with pulsar timing arrays [136]) or via an ALP-photons coupling (e.g.
by radio signals from collisions between miniclusters and neutron stars [137,
138, 139, 140]).202020We note that the same theory can also lead to
interesting signals if the string network survives until the time when the
cosmic microwave background is formed and the axion makes up only a small
fraction of DM [141].

Although large parts of our analysis have rested on analytic expectations, we
have made some progress in studying the string-domain wall system with
numerical simulations. In particular, we have provided evidence that
$\mathcal{A}_{\star}$ increases proportionally to $\log(m_{r}/H_{\star})$, see
Figure 2. As a result, when extrapolated to the physical point
($\log(m_{r}/H_{\star})=\mathcal{O}(100)$), $\mathcal{A}_{\star}$ is likely to
take a value that is at least an order of magnitude larger than occurs in
simulations. Such results are also directly applicable to the QCD axion, for
which $\log(m_{r}/H_{\star})\simeq 65$ and $\mathcal{A}_{\star}/N\simeq 10$,
although the non-linear evolution of the axion waves at around $H=H_{\star}$
in this case could affect the domain wall dynamics.

Of course many uncertainties remain. The evolution of the string-domain wall
network and the area parameter at $H\ll H_{\star}$ in $N>1$ theories is
uncertain, and as a result the DM abundance is not known. Similarly, the
density power spectrum $\mathcal{P}_{\delta}$ of the relic axions is not
reliably known. Although not a problem in our present work, the relation
between the PQ symmetry breaking Lagrangian parameter $\epsilon$ and $H_{d}$
is also unknown apart from an analytic guess. This is particularly important
for the QCD axion in the postinflationary scenario, because it determines
whether the PQ breaking potential must be fine tuned, see e.g. [51, 97, 142,
60, 143].

The theory we focused on is quite minimal from a high energy perspective. All
that is required is an axion potential with $N>1$ degenerate minima, a small
breaking of the remaining discrete $\mathbb{Z}_{N}$, and the PQ symmetry to be
unbroken at the end of, or after, inflation. The first of these is typical for
many sources of axion mass, e.g. if it comes from a strongly coupled hidden
sector [144]. The second is not only easily accommodated but is thought to be
inevitable [34, 35, 36]. Finally, for the $f_{a}$ we have considered, the PQ
symmetry can easily be unbroken at the end of or after inflation: either due
to quantum fluctuations during inflation, finite temperature effects after
reheating, or a coupling between the inflaton and the radial mode [145, 146,
147, 148, 55, 149]. Additionally, it has recently been shown that, over most
of the post-inflationary ALP parameter space in Figure 4, an ALP with a
coupling to photons of the form eq. (3) with $E/N\simeq\mathcal{O}(1)$ is only
possible if the Standard Model gauge group does not arise from a grand unified
theory with a QCD axion unless the Lagrangian is fine-tuned [150].
Consequently, as well as being relatively easily discoverable, a post-
inflationary ALP could give remarkable information not only about the
cosmological history of the Universe (because the PQ symmetry was unbroken
after inflation) but also the high energy completion of the Standard Model.

## Acknowledgements

We thank Asimina Arvanitaki for useful discussions and Giovanni Villadoro for
discussions and collaboration on related work. EH thanks Jed Thompson, David
Cyncynates, and Olivier Simon for useful discussions, and acknowledges the UK
Science and Technology Facilities Council for support through the Quantum
Sensors for the Hidden Sector collaboration under the grant ST/T006145/1 and
UK Research and Innovation Future Leader Fellowship MR/V024566/1. We
acknowledge use of the University of Liverpool Barkla HPC cluster and of the
CINECA Marconi Skylake partition.

##  Appendix A More details on the pre-inflationary scenario

Here we summarise some standard results about ALP dark matter in the pre-
inflationary scenario, focusing on the regime where the ALP is initially close
to the top of a temperature independent cosine potential, i.e.
$\pi-|\theta_{0}|\ll 1$ (a similar analysis can be found in e.g. [118]). Such
an initial condition does not correspond to a fine-tuning of a Lagrangian
parameter. Rather, it is a ‘cosmological tuning’ in that our observable
Universe corresponds to a rare part of the pre-inflationary, presumably
random, axion field distribution.

For $\pi-|\theta_{0}|\ll 1$, numerical results for the function $h$ that
enters the axion relic abundance in eq. (4) can be approximated by (to an
accuracy of better than $10\%$ for $\pi-|\theta_{0}|\lesssim 10^{-3}$)

$h(\theta_{0})\simeq\log^{\,\gamma}\left(\frac{1}{\pi-|\theta_{0}|}\right)~{},\qquad\gamma\simeq
1.4\,,$ (35)

where $\gamma$ is fit numerically. Therefore, in this limit, to obtain full DM
abundance requires

$\pi-|\theta_{0}|\simeq\exp\left(-15\left(\frac{10^{12}~{}{\,{\rm
GeV}}}{f_{a}}\right)^{2\gamma}\left(\frac{10^{-6}~{}{\,{\rm
eV}}}{m_{a}}\right)^{\gamma/2}\right)~{}.$ (36)

Therefore for a fixed $m_{a}$ the axion decay constant decreases only
logarithmically with $\pi-|\theta_{0}|$ as

$f_{a}\propto\frac{1}{\log(1/(\pi-|\theta_{0}|))^{\gamma/2}}~{}.$ (37)

The smallest allowed $\pi-|\theta_{0}|$ in eq. (37) is determined by
isocurvature bounds. Any scalar field that has mass smaller than $H_{I}$
during inflation obtains fluctuations during inflation that are of the order
of $H_{I}/2\pi$. For the axion these are isocurvature and translate to
fluctuations $\delta\theta_{0}$ in the misalignment angle of

$\delta\theta_{0}\simeq\frac{H_{I}}{2\pi f_{a}}~{}.$ (38)

This induces isocurvature perturbations in the axion dark matter energy
density $\rho_{a}$ at all scales of

$\displaystyle\frac{\delta\rho_{a}}{\rho_{a}}$
$\displaystyle\simeq\frac{1}{\rho_{a}}\frac{\partial\rho_{a}}{\partial\theta_{0}}\delta\theta_{0}$
(39) $\displaystyle\simeq\frac{H_{I}}{2\pi
f_{a}}\frac{1}{\rho_{a}}\frac{\partial\rho_{a}}{\partial\theta_{0}}~{},$

and an isocurvature contribution to $\mathcal{P}_{\delta}$ (defined in eq.
(17)) that is independent of $k$. Observational constraints require
$\delta\rho_{a}/\rho_{a}<\delta_{\rm max}\simeq 10^{-5.5}$ [102], so, in the
large misalignment limit,

$\frac{\gamma}{(\pi-|\theta_{0}|)\log(1/(\pi-|\theta_{0}|))}\frac{H_{I}}{2\pi
f_{a}}\lesssim\delta_{\rm max}~{},$ (40)

which, in combination with eq. (36), leads to a lower bound on $f_{a}$ for a
fixed $m_{a}$ if $\Omega_{a}=\Omega_{\rm DM}$. For $H_{I}\simeq{\,{\rm eV}}$
values $\pi-|\theta_{0}|\gtrsim 10^{-18}$ are allowed. In principle somewhat
smaller values of $H_{I}$ are compatible with cosmological observations
provided reheating is fast enough that the Universe is radiation dominated
before big bang nucleosynthesis. However, if $H_{I}\lesssim m_{a}$ our
assumption that the axion starts oscillating at $H\simeq m_{a}$ during
radiation domination fails and the expression for the relic abundance will be
modified in a way that depends on the particular cosmological history.

For completeness, we note that in the opposite regime $\theta_{0}\lesssim 1$,
$\rho_{a}\propto m_{a}^{2}f_{a}^{2}\theta_{0}^{2}$, so if the axion makes up
all the dark matter the isocurvature constraint is $f_{a}\theta_{0}\gtrsim
H_{I}/(2\pi\delta_{\rm max})$. Thus, using eq. (4) with $h(\theta_{0})=1$, the
isocurvature bound is $m_{a}\lesssim 2.5\cdot 10^{-2}~{}{\,{\rm
eV}}\left({\,{\rm eV}}/H_{I}\right)^{1/4}$.

##  Appendix B The relic abundance from strings

Prior to the time when $H=H_{\star}$ the string network should emit most of
the axions at momentum $x_{0}H$, with $x_{0}$ a numerical factor that is
expected to be $\simeq 10$ based on simulation results. Consequently, the
axion number density at a time $t>t_{\star}$ as a result of emission prior to
$H=H_{\star}$ is

$\displaystyle n_{a}^{\rm st}(t)$ $\displaystyle\simeq\int_{\log t_{0}}^{\log
t_{\star}}\frac{2\xi}{x_{0}}\log\left(\frac{m_{r}}{H(t^{\prime})}\right)~{}\frac{\pi
f_{a}^{2}}{t^{\prime}}\left(\frac{R(t^{\prime})}{R(t)}\right)^{3}~{}d\log
t^{\prime}$ (41) $\displaystyle\simeq\frac{8\xi_{\star}}{x_{0}}\pi
m_{a}f_{a}^{2}\log_{\star}\left(\frac{R_{\star}}{R(t)}\right)^{3}~{},$

where we used $t_{0}\ll t_{\star}$ in the second line. Eq. (41) leads to an
axion relic abundance given by eq. (10). On the other hand, between the times
when $H=H_{\star}$ and $H=H_{d}$ the majority of the energy is expected to be
emitted into axions that are non-relativistic or at most mildly relativistic.
As a result, at a time $t_{\star}<t<t_{d}\equiv 1/(2H_{d})$ the axion number
density from emission from strings after $H=H_{\star}$ is

$\displaystyle n_{a}^{\rm st}(t)$ $\displaystyle\simeq\int_{\log
t_{\star}}^{\log
t}2\xi\log\left(\frac{m_{r}}{H(t^{\prime})}\right)~{}\frac{\pi
f_{a}^{2}}{m_{a}t^{\prime
2}}\left(\frac{R(t^{\prime})}{R(t)}\right)^{3}~{}d\log t^{\prime}$ (42)
$\displaystyle\simeq 16\xi_{\star}\pi
m_{a}f_{a}^{2}\log_{\star}\left(1-\left(\frac{H(t)}{m_{a}}\right)^{1/2}\right)\left(\frac{R_{\star}}{R(t)}\right)^{3}~{}.$

Because $t_{\star}=1/(2m_{a})$, from eqs. (41) and (42) the majority of the
axions are produced soon before and soon after the time when $H=H_{\star}$.
Comparing with eq. (13), for $H_{d}\ll m_{a}$, $n_{a}^{\rm st}$ is much
smaller than the number density emitted by domain walls, $n_{a}^{\rm w}$,
during the same period. In particular, this is true once $H_{d}/m_{a}$ is
small enough to overcome the logarithmic enhancement in eqs. (41) and (42).

##  Appendix C Simulations

###  C.1 Details of simulations

We solve the equations of motion of a complex scalar field $\phi$ with
Lagrangian

$\mathcal{L}=|\partial_{\mu}\phi|^{2}-V(\phi)~{},$ (43)

in a Friedmann-Robertson-Walker background with scale factor $R(t)\propto
t^{1/2}$, where $t$ is the cosmic time and the Hubble parameter $H=1/(2t)$.
The potential $V$ is of the general form of eq. (6), and we choose the
particular example

$V\left(\phi\right)=\frac{m_{r}^{2}}{2v^{2}}\left(|\phi|^{2}-\frac{v^{2}}{2}\right)^{2}+\frac{m_{a}^{2}v^{2}}{N^{2}}\left(1-\frac{\sqrt{2}|\phi|}{v}\cos\left(N\frac{a}{v}\right)\right)+{\tilde{\epsilon}}^{2}m_{a}^{2}v^{2}\left(1-\frac{\sqrt{2}|\phi|}{v}\cos\left(\frac{a}{v}\right)\right)~{},$
(44)

where $\phi=\frac{v+r}{\sqrt{2}}e^{ia/v}$. We follow [61] in multiplying the
cosine in eq. (44) by factor of $\sqrt{2}|\phi|/v$ to avoid numerical
instabilities (this has no effect away from the string cores). The magnitude
of the final, $\mathbb{Z}_{N}$ breaking, term is parametrised by
$\tilde{\epsilon}$, related to $\epsilon$ defined in the main text by
$\tilde{\epsilon}=v\sqrt{\epsilon}/(Nm_{a})$.212121The choice that the minima
of the $\mathbb{Z}_{N}$ preserving and breaking contributions to the potential
coincide does not affect the dynamics. The EoM depend on $m_{a}$ only through
the ratio $m_{r}/m_{a}$, so we report different values of the axion mass
always in terms of this. We solve the EoM on a discrete lattice of constant
comoving size using a standard finite difference algorithm, typically on grids
containing $\sim 3000^{3}$ points.

The maximum ratios $m_{r}/H$ and $m_{a}/H$ that can be reached in simulations
are limited by the combination of three requirements. First, the physical
lattice spacing $\Delta$ must be sufficiently small that the string and domain
wall cores are resolved, i.e. $\Delta m_{r}\lesssim 1$ and $\Delta
m_{a}\lesssim 1$ at all times. We pick $\Delta m_{r}=1$, which we have checked
is sufficient to not introduce the systematic uncertainties [53]. Second,
while the network is present there must be at least a few Hubble patches
within the simulation, i.e. $HL\gtrsim 1$ where $L$ is the physical box size.
We run simulations to $HL=2$ (in Appendix C.3 we show this is safe from finite
volume systematics). As a result,
$\log(m_{r}/H_{\star})\lesssim\log(m_{r}/H)\lesssim 7$ for the grid sizes that
we use. Finally, $m_{a}/m_{r}$ must be sufficiently smaller than one so that
the potential in eq. (44) leads to domain walls across which the radial mode
stays close to its vacuum expectation value [151] , rather than interpolating
across the top of its potential (such configurations are not the physical
domain walls). In practice we find $m_{r}/m_{a}\geq 5$ is sufficient, which in
turn fixes $\log(m_{a}/H)\lesssim 5$.

We focus on the ‘fat string and domain wall’ system in which $m_{r}(t)\propto
m_{a}(t)\propto t^{-1/2}$, and only consider the evolution of the physical
theory (with $m_{r}$ constant) in simulations of the string scaling regime
(i.e. prior to the time when $H=H_{\star}$). In the fat string system, the
number of lattice points inside string cores and domain walls remains constant
in time. As mention in the main text, the fat string and domain wall system
reaches the attractor faster. Moreover, modes with physical momentum
$k/R=m_{r}$ or $k/R=m_{a}$ correspond to the same comoving momentum regardless
of when they are produced, so energy emitted with momentum close to the string
core scale ($m_{r}$) does not contaminate the spectrum at smaller momenta,
during the scaling regime and the domain wall evolution.

Notice that the fat domain walls mean that the energy density emitted by the
domain wall network per Hubble time evolves as $\Gamma/H\simeq\sigma H\propto
t^{-3/2}$ rather than the physical $\propto t^{-1}$, see eq. (12). However,
because the energy density in axion waves redshifts $\propto R^{-4}$, the
relic abundance is still dominantly produced at late times. Additionally, the
energy density in domain walls $\simeq Hm_{a}f_{a}^{2}$ still dominates that
in strings $\simeq\pi f_{a}^{2}\log(m_{r}/H)H^{2}$ when $H\ll m_{a}$.

We identify strings and domain walls using the algorithms in [151] and [152,
61] respectively and calculate the scaling parameters $\xi$ and $\mathcal{A}$
from their definitions given in Section 2. As in the main text, when plotting
results, it is convenient to measure the time in simulations in terms of
$\log(m_{r}/H)$ or $\log(m_{a}/H)$. Similarly to [53, 56]
$\partial\rho_{a}/\partial k$ can be written in terms of the Fourier
transforms $\tilde{a}({\bf k})$ and $\tilde{\dot{a}}({\bf k})$ of the axion
field $a({\bf x})$ and its time derivative $\dot{a}({\bf x})$ respectively:

$\frac{\partial\rho_{a}}{\partial{k/R}}=\frac{k^{2}/R^{2}}{(2\pi)^{3}\mathcal{V}}\int
d\Omega\frac{1}{2}\left(|\tilde{\dot{a}}({\bf k})|^{2}+\left(\frac{{\bf
k}^{2}}{R^{2}}+m_{a}^{2}\right)|\dot{a}({\bf k})|^{2}\right)~{},$ (45)

where $k$ is the magnitude of the axion comoving momentum ${\bf k}$,
$\mathcal{V}$ is the simulation volume, and $\Omega$ is the solid angle.

###  C.2 The scaling parameters before network destruction

In this Appendix we show simulation results for the string length $\xi$ and
the domain wall area parameter $\mathcal{A}$ during the string scaling regime
and after $H_{\star}$, as discussed in Section 2.2.

Figure 5: Left: The string length scaling parameters $\xi$ as a function of
time represented by $\log(m_{a}/H)$ for different values of the axion mass,
parametrized by $\log(m_{r}/H_{\star})=\log(m_{r}/m_{a})$. We set the number
of domain wall to $N=4$ and no explicit $\mathbb{Z}_{N}$ symmetry breaking and
study the fat string system. Statistical uncertainties from multiple
simulation runs are shaded. Dashed lines show the evolution of $\xi$ with
$m_{a}=0$. While $\log(m_{a}/H)\lesssim 3.5$, $\xi$ grows as in the $m_{a}=0$
case unaffected by the axion mass and subsequently $\xi$ deviates, affected by
the axion potential and the (stable) network of domain walls that this leads
to. Right: The ratio $\mathcal{A}/\xi$ for the fat and physical string systems
with $m_{a}=0$ and $N=4$ and, for the fat system, also with non-zero $m_{a}$,
again parameterised by $\log(m_{r}/H_{\star})=\log(m_{r}/m_{a})$, and with no
explicit $\mathbb{Z}_{N}$ symmetry breaking.

In Figure 5 (left) we plot the evolution of $\xi$ in the fat system for
different values of the axion mass, normalised as $m_{a}/m_{r}$ (i.e.
different $\log(m_{r}/H_{\star})$) with $\tilde{\epsilon}=0$ (meaning no
$\mathbb{Z}_{N}$ breaking) and $N=4$. This is analogous to the plot of
$\mathcal{A}$ in Figure 2 (left). For comparison, we also show the evolution
in the absence of the axion mass, corresponding to the string scaling regime.
For sufficiently small $\log(m_{a}/H)$, $\xi$ has the well-known growth
proportional to $\log(m_{r}/H)$ [151, 153, 53, 154]. $\xi$ deviates from the
$m_{a}=0$ evolution at $\log(m_{r}/H_{\star})\gtrsim 3$; indicating that the
domain walls that bound the strings at such times are affecting the dynamics.
We do not attempt to extrapolate $\xi$’s evolution to large $\log(m_{a}/H)$
and $\log(m_{a}/H_{\star})$. In Figure 5 (right) we plot the ratio
$\mathcal{A}/\xi$ as a function of $\log(m_{r}/H)$ for the fat and physical
string systems with $m_{a}=0$ (black lines) and also with $m_{a}\neq 0$
(coloured lines) for the fat strings, all with $N=4$. For fat strings this
ratio seems approximately constant while $m_{a}\ll H$. Meanwhile, for the
physical system with $m_{a}=0$ the evolution is less clear but is consistent
with a slow approach to a constant ratio. The ratio in simulations with
$m_{a}\neq 0$ deviates from the $m_{a}=0$ evolution.

Finally, we give more details on the fit of $\mathcal{A}$ prior to the time
when $H=H_{\star}$, discussed in the main text. For the fat string system, the
form in eq. (19) fits the data well, with $c_{1}$ consistent regardless of the
$\log(m_{r}/H)$ in the range $3\div 5.5$ at which the fit is started. A
quadratic ansatz of the kind
$\mathcal{A}(t)=c_{2}\log^{2}+c_{1}\log+c_{0}+c_{-1}\log^{-1}+c_{-2}\log^{-2}$
results in $c_{2}=-0.002\pm.02$ compatible with zero but with a large
uncertainty. The linear form in eq. (19) also fits the physical data well, but
the best fit $c_{1}$ varies depending on the $\log(m_{r}/H)$ from which the
fit is started. Considering the change in $c_{1}$ for initial $\log(m_{r}/H)$
between $3$ and $5$, we estimate $c_{1}=0.36(10)$. We also note that in the
physical case a quadratic fit of the previous form gives $c_{2}\simeq
0.15(10)$, slightly incompatible with zero.

###  C.3 Systematic uncertainties

Figure 6: Comparison between the domain wall area parameter $\mathcal{A}$ and
string length parameter $\xi$ as a function of time from simulations in
different box sizes, but which are otherwise identical (with $m_{r}/m_{a}=5$,
$\tilde{\epsilon}=0$ and $N=4$). This different simulations reach $HL=2$ at
different times (indicated on the lower axis). The agreement between each
simulation and the largest simulation (green), which contains many Hubble
patches at those points, indicates that $\mathcal{A}$ and $\xi$ are unaffected
by finite volume systematics until at least $HL=2$.

Systematic uncertainties can arise in simulations in several ways: the finite
lattice spacing, the finite volume, and the finite timestep. We have checked
that our choices of simulation parameters are sufficient for all the
observables that we extract from simulations to be unaffected. Since the
analysis is very similar to that in [53], in this Appendix we simply
demonstrate that our choice to end simulations when $HL=2$ is sufficient for
finite volume simulations to not affect $\mathcal{A}$ (this is particularly
important because the hint of a decrease in $\mathcal{A}$ in Figure 2 right
happens at the final times). To do so, in Figure 6 we plot $\mathcal{A}$ and
$\xi$ as a function of time for identical theories, but in boxes of different
sizes, so that $HL=2$ happens at different $\log(m_{r}/H)$. After averaging
over multiple runs, we see that there is no deviation in the mean
$\mathcal{A}$ and $\xi$ in the smaller simulations until well after they have
$HL=2$. We have checked this remains the case regardless of $m_{a}/m_{r}$, $N$
and $\tilde{\epsilon}$.

## References

  * [1] R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [2] S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978) 223–226.
  * [3] F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * [4] J. E. Kim, “Weak Interaction Singlet and Strong CP Invariance,” Phys. Rev. Lett. 43 (1979) 103.
  * [5] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, “Can Confinement Ensure Natural CP Invariance of Strong Interactions?,” Nucl. Phys. B 166 (1980) 493–506.
  * [6] A. R. Zhitnitsky, “On Possible Suppression of the Axion Hadron Interactions. (In Russian),” Sov. J. Nucl. Phys. 31 (1980) 260.
  * [7] M. Dine, W. Fischler, and M. Srednicki, “A Simple Solution to the Strong CP Problem with a Harmless Axion,” Phys. Lett. B 104 (1981) 199–202.
  * [8] P. Svrcek and E. Witten, “Axions In String Theory,” JHEP 06 (2006) 051, arXiv:hep-th/0605206.
  * [9] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, “String Axiverse,” Phys. Rev. D 81 (2010) 123530, arXiv:0905.4720 [hep-th].
  * [10] M. Cicoli, M. Goodsell, and A. Ringwald, “The type IIB string axiverse and its low-energy phenomenology,” JHEP 10 (2012) 146, arXiv:1206.0819 [hep-th].
  * [11] B. S. Acharya, K. Bobkov, and P. Kumar, “An M Theory Solution to the Strong CP Problem and Constraints on the Axiverse,” JHEP 11 (2010) 105, arXiv:1004.5138 [hep-th].
  * [12] M. Demirtas, N. Gendler, C. Long, L. McAllister, and J. Moritz, “PQ Axiverse,” arXiv:2112.04503 [hep-th].
  * [13] J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B 120 (1983) 127–132.
  * [14] L. F. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B 120 (1983) 133–136.
  * [15] M. Dine and W. Fischler, “The Not So Harmless Axion,” Phys. Lett. B 120 (1983) 137–141.
  * [16] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, “WISPy Cold Dark Matter,” JCAP 06 (2012) 013, arXiv:1201.5902 [hep-ph].
  * [17] G. Alonso-Álvarez, R. S. Gupta, J. Jaeckel, and M. Spannowsky, “On the Wondrous Stability of ALP Dark Matter,” JCAP 03 (2020) 052, arXiv:1911.07885 [hep-ph].
  * [18] P. Sikivie, “Experimental Tests of the Invisible Axion,” Phys. Rev. Lett. 51 (1983) 1415–1417. [Erratum: Phys.Rev.Lett. 52, 695 (1984)].
  * [19] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, “Experimental Searches for the Axion and Axion-Like Particles,” Ann. Rev. Nucl. Part. Sci. 65 (2015) 485–514, arXiv:1602.00039 [hep-ex].
  * [20] I. G. Irastorza and J. Redondo, “New experimental approaches in the search for axion-like particles,” Prog. Part. Nucl. Phys. 102 (2018) 89–159, arXiv:1801.08127 [hep-ph].
  * [21] H. Georgi, D. B. Kaplan, and L. Randall, “Manifesting the Invisible Axion at Low-energies,” Phys. Lett. B 169 (1986) 73–78.
  * [22] G. Grilli di Cortona, E. Hardy, J. Pardo Vega, and G. Villadoro, “The QCD axion, precisely,” JHEP 01 (2016) 034, arXiv:1511.02867 [hep-ph].
  * [23] M. Gorghetto and G. Villadoro, “Topological Susceptibility and QCD Axion Mass: QED and NNLO corrections,” JHEP 03 (2019) 033, arXiv:1812.01008 [hep-ph].
  * [24] M. Farina, D. Pappadopulo, F. Rompineve, and A. Tesi, “The photo-philic QCD axion,” JHEP 01 (2017) 095, arXiv:1611.09855 [hep-ph].
  * [25] L. Darmé, L. Di Luzio, M. Giannotti, and E. Nardi, “Selective enhancement of the QCD axion couplings,” Phys. Rev. D 103 no. 1, (2021) 015034, arXiv:2010.15846 [hep-ph].
  * [26] A. V. Sokolov and A. Ringwald, “Photophilic hadronic axion from heavy magnetic monopoles,” JHEP 06 (2021) 123, arXiv:2104.02574 [hep-ph].
  * [27] M. Redi and A. Tesi, “The meso-inflationary QCD axion,” arXiv:2211.06421 [hep-ph].
  * [28] K. Harigaya and L.-T. Wang, “More axions from diluted domain walls,” arXiv:2211.08289 [hep-ph].
  * [29] D. Cyncynates, T. Giurgica-Tiron, O. Simon, and J. O. Thompson, “Resonant nonlinear pairs in the axiverse and their late-time direct and astrophysical signatures,” Phys. Rev. D 105 no. 5, (2022) 055005, arXiv:2109.09755 [hep-ph].
  * [30] G. Alonso-Álvarez and J. Jaeckel, “Exploring axionlike particles beyond the canonical setup,” Phys. Rev. D 98 no. 2, (2018) 023539, arXiv:1712.07500 [hep-ph].
  * [31] G. B. Gelmini, A. Simpson, and E. Vitagliano, “Gravitational waves from axionlike particle cosmic string-wall networks,” Phys. Rev. D 104 no. 6, (2021) 061301, arXiv:2103.07625 [hep-ph].
  * [32] G. B. Gelmini, A. Simpson, and E. Vitagliano, “Catastrogenesis: DM, GWs, and PBHs from ALP string-wall networks,” arXiv:2207.07126 [hep-ph].
  * [33] R. Zambujal Ferreira, A. Notari, O. Pujolàs, and F. Rompineve, “High Quality QCD Axion at Gravitational Wave Observatories,” arXiv:2107.07542 [hep-ph].
  * [34] L. F. Abbott and M. B. Wise, “Wormholes and Global Symmetries,” Nucl. Phys. B 325 (1989) 687–704.
  * [35] T. Banks and N. Seiberg, “Symmetries and Strings in Field Theory and Gravity,” Phys. Rev. D 83 (2011) 084019, arXiv:1011.5120 [hep-th].
  * [36] D. Harlow and H. Ooguri, “Symmetries in quantum field theory and quantum gravity,” Commun. Math. Phys. 383 no. 3, (2021) 1669–1804, arXiv:1810.05338 [hep-th].
  * [37] M. Kamionkowski and J. March-Russell, “Planck scale physics and the Peccei-Quinn mechanism,” Phys. Lett. B 282 (1992) 137–141, arXiv:hep-th/9202003.
  * [38] S. M. Barr and D. Seckel, “Planck scale corrections to axion models,” Phys. Rev. D 46 (1992) 539–549.
  * [39] S. Ghigna, M. Lusignoli, and M. Roncadelli, “Instability of the invisible axion,” Phys. Lett. B 283 (1992) 278–281.
  * [40] B. Rai and G. Senjanovic, “Gravity and domain wall problem,” Phys. Rev. D 49 (1994) 2729–2733, arXiv:hep-ph/9301240.
  * [41] R. Holman, S. D. H. Hsu, T. W. Kephart, E. W. Kolb, R. Watkins, and L. M. Widrow, “Solutions to the strong CP problem in a world with gravity,” Phys. Lett. B 282 (1992) 132–136, arXiv:hep-ph/9203206.
  * [42] M. Dine, “Problems of naturalness: Some lessons from string theory,” in Conference on Topics in Quantum Gravity.  7, 1992.  arXiv:hep-th/9207045.
  * [43] B. A. Dobrescu, “The Strong CP problem versus Planck scale physics,” Phys. Rev. D 55 (1997) 5826–5833, arXiv:hep-ph/9609221.
  * [44] P. Cox, T. Gherghetta, and M. D. Nguyen, “A Holographic Perspective on the Axion Quality Problem,” JHEP 01 (2020) 188, arXiv:1911.09385 [hep-ph].
  * [45] S. Fichet and P. Saraswat, “Approximate Symmetries and Gravity,” JHEP 01 (2020) 088, arXiv:1909.02002 [hep-th].
  * [46] W. Yin, “Scale and quality of Peccei-Quinn symmetry and weak gravity conjectures,” JHEP 10 (2020) 032, arXiv:2007.13320 [hep-ph].
  * [47] A. Banerjee, J. Eby, and G. Perez, “From axion quality and naturalness problems to a high-quality $\mathbb{Z}_{4\mathcal{N}}$ QCD relaxion,” arXiv:2210.05690 [hep-ph].
  * [48] T. W. B. Kibble, “Topology of Cosmic Domains and Strings,” J. Phys. A 9 (1976) 1387–1398.
  * [49] A. Vilenkin, “Cosmic Strings,” Phys. Rev. D 24 (1981) 2082–2089.
  * [50] A. Vilenkin and A. E. Everett, “Cosmic Strings and Domain Walls in Models with Goldstone and PseudoGoldstone Bosons,” Phys. Rev. Lett. 48 (1982) 1867–1870.
  * [51] P. Sikivie, “Of Axions, Domain Walls and the Early Universe,” Phys. Rev. Lett. 48 (1982) 1156–1159.
  * [52] P. Sikivie, “Axion Cosmology,” Lect. Notes Phys. 741 (2008) 19–50, arXiv:astro-ph/0610440.
  * [53] M. Gorghetto, E. Hardy, and G. Villadoro, “Axions from Strings: the Attractive Solution,” JHEP 07 (2018) 151, arXiv:1806.04677 [hep-ph].
  * [54] M. Buschmann, J. W. Foster, A. Hook, A. Peterson, D. E. Willcox, W. Zhang, and B. R. Safdi, “Dark matter from axion strings with adaptive mesh refinement,” Nature Commun. 13 no. 1, (2022) 1049, arXiv:2108.05368 [hep-ph].
  * [55] M. Gorghetto, E. Hardy, and H. Nicolaescu, “Observing invisible axions with gravitational waves,” JCAP 06 (2021) 034, arXiv:2101.11007 [hep-ph].
  * [56] M. Gorghetto, E. Hardy, and G. Villadoro, “More axions from strings,” SciPost Phys. 10 no. 2, (2021) 050, arXiv:2007.04990 [hep-ph].
  * [57] Y. B. Zeldovich, I. Y. Kobzarev, and L. B. Okun, “Cosmological Consequences of the Spontaneous Breakdown of Discrete Symmetry,” Zh. Eksp. Teor. Fiz. 67 (1974) 3–11.
  * [58] G. B. Gelmini, M. Gleiser, and E. W. Kolb, “Cosmology of Biased Discrete Symmetry Breaking,” Phys. Rev. D 39 (1989) 1558.
  * [59] S. E. Larsson, S. Sarkar, and P. L. White, “Evading the cosmological domain wall problem,” Phys. Rev. D 55 (1997) 5129–5135, arXiv:hep-ph/9608319.
  * [60] A. Ringwald and K. Saikawa, “Axion dark matter in the post-inflationary Peccei-Quinn symmetry breaking scenario,” Phys. Rev. D 93 no. 8, (2016) 085031, arXiv:1512.06436 [hep-ph]. [Addendum: Phys.Rev.D 94, 049908 (2016)].
  * [61] T. Hiramatsu, M. Kawasaki, K. Saikawa, and T. Sekiguchi, “Axion cosmology with long-lived domain walls,” JCAP 01 (2013) 001, arXiv:1207.3166 [hep-ph].
  * [62] A. Vaquero, J. Redondo, and J. Stadler, “Early seeds of axion miniclusters,” JCAP 04 (2019) 012, arXiv:1809.09241 [astro-ph.CO].
  * [63] C. A. J. O’Hare, G. Pierobon, J. Redondo, and Y. Y. Y. Wong, “Simulations of axionlike particles in the postinflationary scenario,” Phys. Rev. D 105 no. 5, (2022) 055025, arXiv:2112.05117 [hep-ph].
  * [64] E. W. Kolb and M. S. Turner, The Early Universe, vol. 69.  1990\.
  * [65] M. Feix, J. Frank, A. Pargner, R. Reischke, B. M. Schäfer, and T. Schwetz, “Isocurvature bounds on axion-like particle dark matter in the post-inflationary scenario,” JCAP 05 (2019) 021, arXiv:1903.06194 [astro-ph.CO].
  * [66] V. Iršič, H. Xiao, and M. McQuinn, “Early structure formation constraints on the ultralight axion in the postinflation scenario,” Phys. Rev. D 101 no. 12, (2020) 123518, arXiv:1911.11150 [astro-ph.CO].
  * [67] M. Feix, S. Hagstotz, A. Pargner, R. Reischke, B. M. Schäfer, and T. Schwetz, “Post-inflationary axion isocurvature perturbations facing CMB and large-scale structure,” JCAP 11 (2020) 046, arXiv:2004.02926 [astro-ph.CO].
  * [68] P. Bode, J. P. Ostriker, and N. Turok, “Halo formation in warm dark matter models,” Astrophys. J. 556 (2001) 93–107, arXiv:astro-ph/0010389.
  * [69] M. Viel, G. D. Becker, J. S. Bolton, and M. G. Haehnelt, “Warm dark matter as a solution to the small scale crisis: New constraints from high redshift Lyman-$\alpha$ forest data,” Phys. Rev. D 88 (2013) 043502, arXiv:1306.2314 [astro-ph.CO].
  * [70] V. Iršič et al., “New Constraints on the free-streaming of warm dark matter from intermediate and small scale Lyman-$\alpha$ forest data,” Phys. Rev. D 96 no. 2, (2017) 023522, arXiv:1702.01764 [astro-ph.CO].
  * [71] S. Baumholzer, V. Brdar, and E. Morgante, “Structure Formation Limits on Axion-Like Dark Matter,” JCAP 05 (2021) 004, arXiv:2012.09181 [hep-ph].
  * [72] A. Dekker, S. Ando, C. A. Correa, and K. C. Y. Ng, “Warm Dark Matter Constraints Using Milky-Way Satellite Observations and Subhalo Evolution Modeling,” arXiv:2111.13137 [astro-ph.CO].
  * [73] G. Ballesteros, M. A. G. Garcia, and M. Pierre, “How warm are non-thermal relics? Lyman-$\alpha$ bounds on out-of-equilibrium dark matter,” JCAP 03 (2021) 101, arXiv:2011.13458 [hep-ph].
  * [74] R. Diamanti, S. Ando, S. Gariazzo, O. Mena, and C. Weniger, “Cold dark matter plus not-so-clumpy dark relics,” JCAP 06 (2017) 008, arXiv:1701.03128 [astro-ph.CO].
  * [75] M. A. Amin and M. Mirbabayi, “A lower bound on dark matter mass,” arXiv:2211.09775 [hep-ph].
  * [76] C. O’Hare, “cajohare/axionlimits: Axionlimits.” https://cajohare.github.io/AxionLimits/, July, 2020.
  * [77] J. S. Reynés, J. H. Matthews, C. S. Reynolds, H. R. Russell, R. N. Smith, and M. C. D. Marsh, “New constraints on light axion-like particles using Chandra transmission grating spectroscopy of the powerful cluster-hosted quasar H1821+643,” Mon. Not. Roy. Astron. Soc. 510 no. 1, (2021) 1264–1277, arXiv:2109.03261 [astro-ph.HE].
  * [78] C. Dessert, J. W. Foster, and B. R. Safdi, “X-ray Searches for Axions from Super Star Clusters,” Phys. Rev. Lett. 125 no. 26, (2020) 261102, arXiv:2008.03305 [hep-ph].
  * [79] M. Meyer and T. Petrushevska, “Search for Axionlike-Particle-Induced Prompt $\gamma$-Ray Emission from Extragalactic Core-Collapse Supernovae with the $Fermi$ Large Area Telescope,” Phys. Rev. Lett. 124 no. 23, (2020) 231101, arXiv:2006.06722 [astro-ph.HE]. [Erratum: Phys.Rev.Lett. 125, 119901 (2020)].
  * [80] C. Dessert, D. Dunsky, and B. R. Safdi, “Upper limit on the axion-photon coupling from magnetic white dwarf polarization,” arXiv:2203.04319 [hep-ph].
  * [81] CAST Collaboration, V. Anastassopoulos et al., “New CAST Limit on the Axion-Photon Interaction,” Nature Phys. 13 (2017) 584–590, arXiv:1705.02290 [hep-ex].
  * [82] D. Cadamuro and J. Redondo, “Cosmological bounds on pseudo Nambu-Goldstone bosons,” JCAP 02 (2012) 032, arXiv:1110.2895 [hep-ph].
  * [83] D. Wadekar and Z. Wang, “Strong constraints on decay and annihilation of dark matter from heating of gas-rich dwarf galaxies,” arXiv:2111.08025 [hep-ph].
  * [84] M. Regis, M. Taoso, D. Vaz, J. Brinchmann, S. L. Zoutendijk, N. F. Bouché, and M. Steinmetz, “Searching for light in the darkness: Bounds on ALP dark matter with the optical MUSE-faint survey,” Phys. Lett. B 814 (2021) 136075, arXiv:2009.01310 [astro-ph.CO].
  * [85] D. Grin, G. Covone, J.-P. Kneib, M. Kamionkowski, A. Blain, and E. Jullo, “A Telescope Search for Decaying Relic Axions,” Phys. Rev. D 75 (2007) 105018, arXiv:astro-ph/0611502.
  * [86] A. Arvanitaki, M. Baryakhtar, and X. Huang, “Discovering the QCD Axion with Black Holes and Gravitational Waves,” Phys. Rev. D 91 no. 8, (2015) 084011, arXiv:1411.2263 [hep-ph].
  * [87] M. Baryakhtar, M. Galanis, R. Lasenby, and O. Simon, “Black hole superradiance of self-interacting scalar fields,” Phys. Rev. D 103 no. 9, (2021) 095019, arXiv:2011.11646 [hep-ph].
  * [88] G. G. Raffelt, “Astrophysical axion bounds,” Lect. Notes Phys. 741 (2008) 51–71, arXiv:hep-ph/0611350.
  * [89] J. H. Chang, R. Essig, and S. D. McDermott, “Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle,” JHEP 09 (2018) 051, arXiv:1803.00993 [hep-ph].
  * [90] P. Carenza, T. Fischer, M. Giannotti, G. Guo, G. Martínez-Pinedo, and A. Mirizzi, “Improved axion emissivity from a supernova via nucleon-nucleon bremsstrahlung,” JCAP 10 no. 10, (2019) 016, arXiv:1906.11844 [hep-ph]. [Erratum: JCAP 05, E01 (2020)].
  * [91] N. Viaux, M. Catelan, P. B. Stetson, G. Raffelt, J. Redondo, A. A. R. Valcarce, and A. Weiss, “Neutrino and axion bounds from the globular cluster M5 (NGC 5904),” Phys. Rev. Lett. 111 (2013) 231301, arXiv:1311.1669 [astro-ph.SR].
  * [92] M. M. Miller Bertolami, B. E. Melendez, L. G. Althaus, and J. Isern, “Revisiting the axion bounds from the Galactic white dwarf luminosity function,” JCAP 10 (2014) 069, arXiv:1406.7712 [hep-ph].
  * [93] V. Iršič, M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker, “First constraints on fuzzy dark matter from Lyman-$\alpha$ forest data and hydrodynamical simulations,” Phys. Rev. Lett. 119 no. 3, (2017) 031302, arXiv:1703.04683 [astro-ph.CO].
  * [94] V. Lora, J. Magana, A. Bernal, F. J. Sanchez-Salcedo, and E. K. Grebel, “On the mass of ultra-light bosonic dark matter from galactic dynamics,” JCAP 02 (2012) 011, arXiv:1110.2684 [astro-ph.GA].
  * [95] A. X. González-Morales, D. J. E. Marsh, J. Peñarrubia, and L. A. Ureña López, “Unbiased constraints on ultralight axion mass from dwarf spheroidal galaxies,” Mon. Not. Roy. Astron. Soc. 472 no. 2, (2017) 1346–1360, arXiv:1609.05856 [astro-ph.CO].
  * [96] D. J. E. Marsh and J. C. Niemeyer, “Strong Constraints on Fuzzy Dark Matter from Ultrafaint Dwarf Galaxy Eridanus II,” Phys. Rev. Lett. 123 no. 5, (2019) 051103, arXiv:1810.08543 [astro-ph.CO].
  * [97] S. Chang, C. Hagmann, and P. Sikivie, “Studies of the motion and decay of axion walls bounded by strings,” Phys. Rev. D 59 (1999) 023505, arXiv:hep-ph/9807374.
  * [98] M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments.  Oxford Master Series in Physics. Oxford University Press, 2007.
  * [99] T. L. Smith, E. Pierpaoli, and M. Kamionkowski, “A new cosmic microwave background constraint to primordial gravitational waves,” Phys. Rev. Lett. 97 (2006) 021301, arXiv:astro-ph/0603144.
  * [100] L. Pagano, L. Salvati, and A. Melchiorri, “New constraints on primordial gravitational waves from Planck 2015,” Phys. Lett. B 760 (2016) 823–825, arXiv:1508.02393 [astro-ph.CO].
  * [101] M. Kamionkowski and A. Kosowsky, “The Cosmic microwave background and particle physics,” Ann. Rev. Nucl. Part. Sci. 49 (1999) 77–123, arXiv:astro-ph/9904108.
  * [102] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [103] BICEP2, Keck Array Collaboration, P. A. R. Ade et al., “BICEP2 / Keck Array x: Constraints on Primordial Gravitational Waves using Planck, WMAP, and New BICEP2/Keck Observations through the 2015 Season,” Phys. Rev. Lett. 121 (2018) 221301, arXiv:1810.05216 [astro-ph.CO].
  * [104] T. Matsumura et al., “Mission design of LiteBIRD,” J. Low Temp. Phys. 176 (2014) 733, arXiv:1311.2847 [astro-ph.IM].
  * [105] K. Choi and S. H. Im, “Realizing the relaxion from multiple axions and its UV completion with high scale supersymmetry,” JHEP 01 (2016) 149, arXiv:1511.00132 [hep-ph].
  * [106] D. E. Kaplan and R. Rattazzi, “Large field excursions and approximate discrete symmetries from a clockwork axion,” Phys. Rev. D 93 no. 8, (2016) 085007, arXiv:1511.01827 [hep-ph].
  * [107] C. J. Hogan and M. J. Rees, “AXION MINICLUSTERS,” Phys. Lett. B 205 (1988) 228–230.
  * [108] E. W. Kolb and I. I. Tkachev, “Axion miniclusters and Bose stars,” Phys. Rev. Lett. 71 (1993) 3051–3054, arXiv:hep-ph/9303313.
  * [109] K. M. Zurek, C. J. Hogan, and T. R. Quinn, “Astrophysical Effects of Scalar Dark Matter Miniclusters,” Phys. Rev. D 75 (2007) 043511, arXiv:astro-ph/0607341.
  * [110] E. Hardy, “Miniclusters in the Axiverse,” JHEP 02 (2017) 046, arXiv:1609.00208 [hep-ph].
  * [111] M. Fairbairn, D. J. E. Marsh, J. Quevillon, and S. Rozier, “Structure formation and microlensing with axion miniclusters,” Phys. Rev. D 97 no. 8, (2018) 083502, arXiv:1707.03310 [astro-ph.CO].
  * [112] L. Visinelli and J. Redondo, “Axion Miniclusters in Modified Cosmological Histories,” Phys. Rev. D 101 no. 2, (2020) 023008, arXiv:1808.01879 [astro-ph.CO].
  * [113] B. Eggemeier, J. Redondo, K. Dolag, J. C. Niemeyer, and A. Vaquero, “First Simulations of Axion Minicluster Halos,” Phys. Rev. Lett. 125 no. 4, (2020) 041301, arXiv:1911.09417 [astro-ph.CO].
  * [114] D. Ellis, D. J. E. Marsh, and C. Behrens, “Axion Miniclusters Made Easy,” Phys. Rev. D 103 no. 8, (2021) 083525, arXiv:2006.08637 [astro-ph.CO].
  * [115] D. Ellis, D. J. E. Marsh, B. Eggemeier, J. Niemeyer, J. Redondo, and K. Dolag, “Structure of Axion Miniclusters,” arXiv:2204.13187 [hep-ph].
  * [116] V. Dandoy, T. Schwetz, and E. Todarello, “A self-consistent wave description of axion miniclusters and their survival in the galaxy,” arXiv:2206.04619 [astro-ph.CO].
  * [117] X. Shen, H. Xiao, P. F. Hopkins, and K. M. Zurek, “Disruption of Dark Matter Minihaloes in the Milky Way environment: Implications for Axion Miniclusters and Early Matter Domination,” arXiv:2207.11276 [astro-ph.GA].
  * [118] A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, “Large-misalignment mechanism for the formation of compact axion structures: Signatures from the QCD axion to fuzzy dark matter,” Phys. Rev. D 101 no. 8, (2020) 083014, arXiv:1909.11665 [astro-ph.CO].
  * [119] P.-H. Chavanis, “Collapse of a self-gravitating Bose-Einstein condensate with attractive self-interaction,” Phys. Rev. D 94 no. 8, (2016) 083007, arXiv:1604.05904 [astro-ph.CO].
  * [120] P.-H. Chavanis, “Growth of perturbations in an expanding universe with Bose-Einstein condensate dark matter,” Astron. Astrophys. 537 (2012) A127, arXiv:1103.2698 [astro-ph.CO].
  * [121] J.-c. Hwang and H. Noh, “Axion as a Cold Dark Matter candidate,” Phys. Lett. B 680 (2009) 1–3, arXiv:0902.4738 [astro-ph.CO].
  * [122] A. Suárez and P.-H. Chavanis, “Hydrodynamic representation of the Klein-Gordon-Einstein equations in the weak field limit: General formalism and perturbations analysis,” Phys. Rev. D 92 (2015) 023510, arXiv:1503.07437 [gr-qc].
  * [123] L. Hui, “Wave Dark Matter,” Ann. Rev. Astron. Astrophys. 59 (2021) 247–289, arXiv:2101.11735 [astro-ph.CO].
  * [124] M. Gorghetto, E. Hardy, J. March-Russell, N. Song, and S. M. West, “Dark Photon Stars: Formation and Role as Dark Matter Substructure,” arXiv:2203.10100 [hep-ph].
  * [125] D. J. E. Marsh, K.-C. Fong, E. W. Lentz, L. Smejkal, and M. N. Ali, “Proposal to Detect Dark Matter using Axionic Topological Antiferromagnets,” Phys. Rev. Lett. 123 no. 12, (2019) 121601, arXiv:1807.08810 [hep-ph].
  * [126] BREAD Collaboration, J. Liu et al., “Broadband solenoidal haloscope for terahertz axion detection,” arXiv:2111.12103 [physics.ins-det].
  * [127] “Brass website.” http://wwwiexp.desy.de/groups/astroparticle/brass/brassweb.htm.
  * [128] B. T. McAllister, G. Flower, E. N. Ivanov, M. Goryachev, J. Bourhill, and M. E. Tobar, “The ORGAN Experiment: An axion haloscope above 15 GHz,” Phys. Dark Univ. 18 (2017) 67–72, arXiv:1706.00209 [physics.ins-det].
  * [129] MADMAX Working Group Collaboration, A. Caldwell, G. Dvali, B. Majorovits, A. Millar, G. Raffelt, J. Redondo, O. Reimann, F. Simon, and F. Steffen, “Dielectric Haloscopes: A New Way to Detect Axion Dark Matter,” Phys. Rev. Lett. 118 no. 9, (2017) 091801, arXiv:1611.05865 [physics.ins-det].
  * [130] M. Lawson, A. J. Millar, M. Pancaldi, E. Vitagliano, and F. Wilczek, “Tunable axion plasma haloscopes,” Phys. Rev. Lett. 123 no. 14, (2019) 141802, arXiv:1904.11872 [hep-ph].
  * [131] D. Alesini, D. Babusci, D. Di Gioacchino, C. Gatti, G. Lamanna, and C. Ligi, “The KLASH Proposal,” arXiv:1707.06010 [physics.ins-det].
  * [132] M. Silva-Feaver et al., “Design Overview of DM Radio Pathfinder Experiment,” IEEE Trans. Appl. Supercond. 27 no. 4, (2017) 1400204, arXiv:1610.09344 [astro-ph.IM].
  * [133] Z. Zhang, O. Ghosh, and D. Horns, “WISPLC: Search for Dark Matter with LC Circuit,” arXiv:2111.04541 [hep-ex].
  * [134] J. L. Ouellet et al., “First Results from ABRACADABRA-10 cm: A Search for Sub-$\mu$eV Axion Dark Matter,” Phys. Rev. Lett. 122 no. 12, (2019) 121802, arXiv:1810.12257 [hep-ex].
  * [135] A. Mitridate, T. Trickle, Z. Zhang, and K. M. Zurek, “Detectability of Axion Dark Matter with Phonon Polaritons and Magnons,” Phys. Rev. D 102 no. 9, (2020) 095005, arXiv:2005.10256 [hep-ph].
  * [136] H. Ramani, T. Trickle, and K. M. Zurek, “Observability of Dark Matter Substructure with Pulsar Timing Correlations,” JCAP 12 (2020) 033, arXiv:2005.03030 [astro-ph.CO].
  * [137] I. I. Tkachev, “Fast Radio Bursts and Axion Miniclusters,” JETP Lett. 101 no. 1, (2015) 1–6, arXiv:1411.3900 [astro-ph.HE].
  * [138] M. S. Pshirkov, “May axion clusters be sources of fast radio bursts?,” Int. J. Mod. Phys. D 26 no. 07, (2017) 1750068, arXiv:1609.09658 [astro-ph.HE].
  * [139] J. H. Buckley, P. S. B. Dev, F. Ferrer, and F. P. Huang, “Fast radio bursts from axion stars moving through pulsar magnetospheres,” Phys. Rev. D 103 no. 4, (2021) 043015, arXiv:2004.06486 [astro-ph.HE].
  * [140] T. D. P. Edwards, B. J. Kavanagh, L. Visinelli, and C. Weniger, “Transient Radio Signatures from Neutron Star Encounters with QCD Axion Miniclusters,” Phys. Rev. Lett. 127 no. 13, (2021) 131103, arXiv:2011.05378 [hep-ph].
  * [141] P. Agrawal, A. Hook, J. Huang, and G. Marques-Tavares, “Axion string signatures: a cosmological plasma collider,” JHEP 01 (2022) 103, arXiv:2010.15848 [hep-ph].
  * [142] M. Kawasaki, K. Saikawa, and T. Sekiguchi, “Axion dark matter from topological defects,” Phys. Rev. D 91 no. 6, (2015) 065014, arXiv:1412.0789 [hep-ph].
  * [143] K. A. Beyer and S. Sarkar, “Ruling out light axions: the writing is on the wall,” arXiv:2211.14635 [hep-ph].
  * [144] H. Georgi, J. E. Kim, and H.-P. Nilles, “Hidden sector gaugino condensation and the model independent axion,” Phys. Lett. B 437 (1998) 325–330, arXiv:hep-ph/9805510.
  * [145] L. A. Kofman and A. D. Linde, “Generation of Density Perturbations in the Inflationary Cosmology,” Nucl. Phys. B 282 (1987) 555.
  * [146] E. T. Vishniac, K. A. Olive, and D. Seckel, “Cosmic Strings and Inflation,” Nucl. Phys. B 289 (1987) 717–734.
  * [147] J. Yokoyama, “INFLATION CAN SAVE COSMIC STRINGS,” Phys. Rev. Lett. 63 (1989) 712.
  * [148] H. M. Hodges and J. R. Primack, “Strings, texture and inflation,” Phys. Rev. D 43 (1991) 3155–3172.
  * [149] Y. Bao, J. Fan, and L. Li, “Opening up window of post-inflationary QCD axion,” arXiv:2209.09908 [hep-ph].
  * [150] P. Agrawal, M. Nee, and M. Reig, “Axion Couplings in Grand Unified Theories,” arXiv:2206.07053 [hep-ph].
  * [151] L. Fleury and G. D. Moore, “Axion dark matter: strings and their cores,” JCAP 01 (2016) 004, arXiv:1509.00026 [hep-ph].
  * [152] W. H. Press, B. S. Ryden, and D. N. Spergel, “Dynamical Evolution of Domain Walls in an Expanding Universe,” apj 347 (Dec., 1989) .
  * [153] V. B. . Klaer and G. D. Moore, “The dark-matter axion mass,” JCAP 11 (2017) 049, arXiv:1708.07521 [hep-ph].
  * [154] V. B. Klaer and G. D. Moore, “Global cosmic string networks as a function of tension,” JCAP 06 (2020) 021, arXiv:1912.08058 [hep-ph].