Early vs late string networks from a minimal QCD Axion
Marco Gorghettoa,b, Edward Hardyc, Horia Nicolaescud, Alessio Notarie,f,
Michele Redig

aDeutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany
bDepartment of Particle Physics and Astrophysics, Weizmann Institute of
Science,
Herzl St 234, Rehovot 761001, Israel
cRudolf Peierls Centre for Theoretical Physics, University of Oxford, Parks
Road, Oxford OX1 3PU, UK
dDepartment of Mathematical Sciences, University of Liverpool,
Liverpool, L69 7ZL, United Kingdom
eDepartament de Física Fondamental i Institut de Ciéncies del Cosmos,
Universitat de Barcelona,
Martí i Franqués 1, E-08028 Barcelona, Spain
fGalileo Galilei Institute for theoretical physics, Centro Nazionale INFN di
Studi Avanzati
Largo Enrico Fermi 2, I-50125, Firenze, Italy
gINFN Sezione di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, Italy
Department of Physics and Astronomy, University of Florence, Italy

Abstract

> We propose a new regime of minimal QCD axion dark matter that lies between
> the pre- and post-inflationary scenarios, such that the Peccei-Quinn (PQ)
> symmetry is restored only on sufficiently large spatial scales. This leads
> to a novel cosmological evolution, in which strings and domain walls re-
> enter the horizon and annihilate later than in the ordinary post-
> inflationary regime, possibly even after the QCD crossover. Such dynamics
> can occur if the PQ symmetry is restored by inflationary fluctuations, i.e.
> the Hubble parameter during inflation $H_{I}$ is larger than the PQ breaking
> scale $f_{a}$, but it is not thermally restored afterwards. Solving the
> Fokker-Planck equation, we estimate the number of inflationary e-folds
> required for the PQ symmetry to be, on average, restored. Moreover, we show
> that, in the large parts of parameter space where the radial mode is
> displaced from the minimum by de Sitter fluctuations, a string network forms
> due to the radial mode oscillating over the top of its potential after
> inflation. In both cases we identify order one ranges in $H_{I}/f_{a}$ and
> in the quartic coupling $\lambda$ of the PQ potential that lead to the late-
> string dynamics. In this regime the cosmological dark matter abundance can
> be reproduced for axion decay constants as low as the astrophysical
> constraint ${\cal O}(10^{8})$ GeV, corresponding to axion masses up to
> $10^{-2}~{}{\rm eV}$, and with miniclusters with masses as large as
> $\mathcal{O}(10)M_{\odot}$.

###### Contents

  1. 1 Introduction
  2. 2 The QCD axion during, and immediately after, inflation
    1. 2.1 Evolution during inflation
    2. 2.2 Evolution soon after inflation
  3. 3 Formation of strings
    1. 3.1 Inflationary formation
    2. 3.2 Overshoot mechanism
  4. 4 Comparison to simulations
    1. 4.1 Results
  5. 5 Isocurvature perturbations
    1. 5.1 Inflationary Correlations
    2. 5.2 Effects of strings on correlations
    3. 5.3 The effect of overshoot on super-horizon correlations
  6. 6 Landscape of axion scenarios
    1. 6.1 Breaking temperature of PQ symmetry
    2. 6.2 Regimes
    3. 6.3 The late string relic abundance
    4. 6.4 Axion spatial distribution and miniclusters
    5. 6.5 Isocurvature constraints
  7. 7 Discussion
  8. A Approximate correlation functions
    1. A.1 Isocurvature estimates
  9. B More on string production
  10. C Further details of the numerical simulations
    1. C.1 Description of simulations
    2. C.2 The string network
  11. D The effect of overshoot on the density power spectrum

##  1 Introduction

The cosmological dynamics of the QCD axion generally falls into one of two
broad classes depending on when the Peccei-Quinn (PQ) symmetry is last
unbroken [1, 2, 3, 4, 5, 6, 7]. In the first class, which goes under the name
of ‘pre-inflationary’ scenario, the PQ symmetry is broken during the visible
part of inflation and never subsequently restored [8, 9, 10]. Conversely, in
the second class, called the ‘post-inflationary scenario’, the PQ symmetry is
unbroken at some time after inflation and breaks subsequently.

The division into the pre- and post-inflationary regimes is relevant because
these possibilities could be distinguished observationally or experimentally,
and they also have differing requirements on the underlying axion theory. One
of the key reasons for this is that topological defects are present in the
early Universe only in the post-inflationary scenario [11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22]. In this case, when the PQ symmetry breaks, a network
of axion strings forms and quickly approaches an attractor, scaling, solution
with roughly one Hubble length of string per Hubble patch [23, 24, 25, 26, 27,
28]. Subsequently, when the axion mass $m_{a}$ reaches $m_{a}\simeq H$ (with
$H$ the Hubble parameter), domain walls form [29]. If the domain wall number
$N_{\rm DW}$ equals unity, the strings ‘chop up’ and quickly destroy the
resulting network. If $N_{\rm DW}>1$, the domain wall network is stable and an
acceptable cosmological history requires additional ingredients, e.g. an
explicit PQ symmetry breaking energy bias between the $N_{\rm DW}$
vacua.111Cosmic strings are essential for a post-inflationary axion to be
consistent with observations; if there were only domain walls present,
infinitely large domain walls would overclose the Universe both for $N_{\rm
DW}=1$ and also $N_{\rm DW}>1$ even with an energy bias. Due to the axions
produced by defects, the pre- and post-inflationary scenarios lead to almost
complementary mass ranges for which an axion can account for the full observed
dark matter abundance [30, 31], assuming a standard cosmological history.
Additionally, the spatial distribution of dark matter axions on small scales
is different; relatively dense ‘miniclusters’ form in the post-inflationary
scenario but are typically absent in the pre-inflationary case [32, 33, 34,
35] (see however [36]). Meanwhile, pre-inflationary axions are compatible also
with $N_{\rm DW}>1$, but in this case a low scale of inflation $H_{I}\ll{\cal
O}(10^{-5})f_{a}$, where $f_{a}$ is the axion decay constant, is needed to
evade observational constraints on isocurvature.

The originally considered way that the post-inflationary regime can occur is
if the maximal temperature of the Universe $T_{\rm max}$ after inflation [37]
is large enough, $T_{\rm max}\gtrsim f_{a}$. In this case the PQ symmetry is
restored by finite-temperature effects and subsequently breaks through a phase
transition, at some time after the end of inflation (for clarity we call this
the thermal post-inflationary scenario) [38, 39, 40, 41]. In this paper we
instead suppose that $H_{I}\gtrsim f_{a}\gg T_{\rm max}$. As a result, finite
temperature effects do not restore the PQ symmetry after inflation. However,
the PQ symmetry might be effectively restored on average on a given spatial
domain _during_ inflation if inflationary fluctuations are large enough to
efficiently randomise the PQ field in that region [42] (even if in almost all
points the symmetry is locally broken), in which case the effective breaking
of the PQ symmetry happens after the end of inflation when the PQ radial mode
relaxes to its minimum (we refer to this as the stochastic post-inflationary
scenario). Such dynamics have been considered with the radial mode sitting in
the vacuum, in which case only domain walls form, [43] and with the radial
mode fluctuating, which can lead to strings, [44, 45], see also [46, 47]. A
naive expectation is that for a Hubble scale during inflation $H_{I}\gtrsim
2\pi f_{a}$ such restoration is inevitable and the resulting dynamics long
after inflation are indistinguishable from the thermal post-inflationary case.
However, building on the analysis of [44, 45], we will see that the situation
is not this simple.

By studying the evolution of the PQ field during inflation, our first aim is
to determine the values of $H_{I}/f_{a}$ and the quartic $\lambda$ of the PQ
potential for which inflationary fluctuations do actually restore the PQ
symmetry, giving rise to the stochastic post-inflationary scenario. We compare
our analytic predictions to results from numerical simulations of the
evolution of the PQ scalar after inflation, which provide evidence for our
approach. For the same purpose, we also study for the first time
quantitatively (as far as we are aware), a mechanism noted in Ref. [45] by
which oscillations of the radial mode over the top of its potential soon after
the end of inflation can amplify initially small differences between the value
of the PQ scalar at different points in space, leading to a string network in
parts of parameter space that would otherwise not.

Our second aim is to analyse the dynamics of the axion string network at the
boundary of the parameter space where inflationary fluctuations restore the PQ
symmetry, again assuming $T_{\rm max}\ll f_{a}$. In this case, the PQ symmetry
might be restored on average only over large spatial scales, much larger than
the Hubble radius at the end of inflation, while in smaller regions the
symmetry is locally broken. This corresponds to an initially underdense string
network that re-enters the horizon late, possibly even after the QCD
crossover, which we call the late string regime.222This effectively realizes
the scenario of Ref. [48, 49] where the PQ symmetry was assumed to break after
the beginning of visible inflation. String networks with similar features
arising as a result of different dynamics during inflation, often a phase
transition, have been considered in various contexts [50, 51, 52, 53, 54, 55,
56, 57, 58, 59, 60, 48].  Such dynamics occur for order-one ranges of
$H_{I}/f_{a}$ in minimal axion theories with no additional model building, and
affects both the dark matter axion mass and the properties of miniclusters. In
particular, for $N_{\rm DW}=1$ the full observed dark matter abundance can be
obtained for axion decay constants as small as allowed by astrophysical
constraints, $f_{a}\gtrsim 10^{8}$ GeV. We also consider isocurvature
perturbations in such a scenario, which might lead to important bounds
although large uncertainties remain.

The structure of our work is as follows. We start in Section 2 by describing
the dynamics of the PQ field during, and soon after, inflation. In Section 3
we analyse how the inflationary dynamics lead to the formation of a string
network and provide an analytic argument to predict when the network re-enters
the horizon. This analytic analysis is confirmed by numerical simulations in
Section 4. In Section 5 we discuss isocurvature perturbations in this
scenario, in Section 6 we map out the landscape of the minimal axion theory we
focus on and discuss the phenomenology and observational signals. Finally in
Section 7 we conclude and consider directions for future work.

##  2 The QCD axion during, and immediately after, inflation

We consider a complex scalar field $\Phi$ described by the Lagrangian

$\mathcal{L}=|\partial_{\mu}\Phi|^{2}-V\ ,\qquad
V=\lambda\left(|\Phi|^{2}-\frac{f_{a}^{2}}{2}\right)^{2},\qquad\Phi=(\phi_{1}+i\phi_{2})/\sqrt{2}~{}.$
(1)

The axion $a(x)$ is the phase of
$\Phi(x)=\frac{\rho(x)}{\sqrt{2}}e^{ia(x)/f_{a}}$, while the radial mode
$\rho\equiv\sqrt{2}|\Phi|$ has mass $m_{\rho}^{2}=2\lambda f^{2}_{a}$. We
allow the quartic to be somewhat but not too much smaller than
$\mathcal{O}(1)$, e.g. $\lambda\gtrsim 10^{-4}$ (the complications for smaller
$\lambda$ are discussed in Section 3.2). The PQ symmetry breaking scale is set
to be $f_{a}$ as is the case in theories with domain wall number $N_{W}=1$,
which give a minimal viable post-inflationary scenario. The potential $V$ in
eq. (1) arises in the realization of the KSVZ model with one vector-like
fermion [4, 5], however our analysis can be adapted to UV completions with
other potentials and also to $N_{W}>1$.333The axion potential from QCD, to be
included into eq. (1), will be irrelevant during inflation as this will occur
at a high scale, with Hubble parameter of order $f_{a}$, in a regime far from
the stochastic axion scenario [61].

###  2.1 Evolution during inflation

We study the evolution of the PQ field during $N$ e-folds of inflation with,
assumed constant, Hubble parameter $H_{I}$ using the stochastic formalism of
inflation [42] (see also [62]). Throughout this time, scalars with mass
smaller than $\sqrt{2}H_{I}$ experience de Sitter fluctuations similar to
thermal fluctuations with effective temperature $T_{\rm dS}=H_{I}/(2\pi)$ [63,
64]. Once a mode exits the horizon it becomes classical with amplitude
$H_{I}/(2\pi)$, giving a kick to average value of the field. Consequently, the
field value in a coarse-grained region of size $H_{I}^{-3}$ undergoes a random
walk, that is independent in causally disconnected regions, as well as a
motion dictated by the classical potential. The combination of these two
processes can be described through the Fokker-Planck (FP) equation, which
gives the evolution of the probability distribution $P(\Phi)$ of the coarse-
grained field sampled over a region of space that contains many Hubble
volumes.

The general form of the FP equation reads

$\frac{\partial P}{\partial
t}=\frac{H_{I}^{3}}{8\pi^{2}}\sum_{i}\left[\partial_{i}\partial_{i}P+\frac{8\pi^{2}}{3H_{I}^{4}}\partial_{i}(P\partial_{i}V)\right]\,,$
(2)

where $t$ is cosmic time and $\partial_{i}$ are the derivatives with respect
to the fields $\phi_{i}$ that are light during inflation. In general, eq. (2)
admits a time-independent ‘equilibrium’ solution

$P_{\rm eq}(\Phi)\equiv A\,e^{\frac{-8\pi^{2}V(\Phi)}{3H_{I}^{4}}}\,,$ (3)

on which the spread of the field by inflationary fluctuations is balanced by
the classical motion. $A$ is a normalisation constant such that $\int
d\phi_{1}d\phi_{2}P=1$.

For the Lagrangian in eq. (1), the FP equation is two-dimensional with
$\Phi=(\phi_{1}+i\phi_{2})/\sqrt{2}$ and is valid as long as
$m_{\rho}<\sqrt{2}H_{I}$ during inflation (otherwise, the formalism can only
be applied to the axion, in which case string loops can only form by
exponentially suppressed tunnelling processes [65]). On the equilibrium
solution, the complex field is uniformly distributed in the angular direction
and

$P_{\rm eq}(\rho)\propto
e^{-\frac{2\pi^{2}}{3}\alpha^{-4}\left(\rho^{2}/f_{a}^{2}-1\right)^{2}}~{},\qquad\quad\alpha\equiv
H_{I}/(\lambda^{1/4}f_{a})\,.$ (4)

For $\alpha\ll 1$ the radial mode is strongly peaked in the vacuum, i.e. the
expectation value $\langle\rho^{2}\rangle=f_{a}^{2}$. Instead for $\alpha\gg
1$ the radial mode is displaced [47]:

$\langle\rho^{2}\rangle=\sqrt{\frac{3}{2\pi^{3}\lambda}}H_{I}^{2}\,,$ (5)

with the dynamics dominated by the quartic term. In Figure 1 (left) we show
$\langle\rho^{2}\rangle$ as a function of $\alpha$.

The asymptotic solution is approached regardless of the initial conditions.
The timescale, or equivalently the number of e-foldings, on which this happens
is a crucial ingredient for our work. To this aim, following [42, 66], we
define the reduced probability distribution

${\tilde{P}}\equiv e^{\frac{4\pi^{2}V}{3H_{I}^{4}}}P~{},$ (6)

which satisfies

$\frac{\partial\tilde{P}}{\partial
t}=\frac{H_{I}^{3}}{8\pi^{2}}\sum_{i}\left[\partial_{i}^{2}-\left(\left(\partial_{i}v\right)^{2}-\partial^{2}_{i}v\right)\right]\tilde{P}\,,~{}~{}~{}~{}~{}~{}v\equiv\frac{4\pi^{2}}{3H_{I}^{4}}V~{}.$
(7)

It is useful to express this in terms of the radial and angular modes $\rho$
and $\theta\equiv a/f_{a}$ to obtain,

$\frac{\partial\tilde{P}}{\partial
t}=\frac{H_{I}^{3}}{8\pi^{2}}\left[\partial^{2}_{\rho}+\frac{1}{\rho}\partial_{\rho}+\frac{1}{\rho^{2}}\partial_{\theta}^{2}-\left(\left(\partial_{\rho}v\right)^{2}-\frac{1}{\rho}\partial_{\rho}v-{\partial^{2}_{\rho}v}\right)\right]\tilde{P}~{}.$
(8)

We can take advantage of the U(1) symmetry and look for solutions of the form
$\tilde{P}\propto\psi_{m}(t,\rho)e^{im\theta}$, with $m$ integer, given the
periodicity of $\theta$. This leads to

$\frac{\partial\psi_{m}}{\partial
t}=\frac{H_{I}^{3}}{8\pi^{2}}\left[\partial_{\rho}^{2}+\frac{1}{\rho}\partial_{\rho}-\frac{m^{2}}{\rho^{2}}-V_{\rm
eff}\right]\psi_{m}\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}V_{\rm
eff}\equiv\left(\partial_{\rho}v\right)^{2}-\frac{1}{\rho}\partial_{\rho}v-{\partial^{2}_{\rho}v}~{}.$
(9)

In this form, the right hand-side can be interpreted as the Hamiltonian of a
one-dimensional quantum mechanical particle with potential $V_{\rm eff}$. In
terms of $x\equiv\lambda^{1/4}\rho/H_{I}$ and $\alpha$, the corresponding
operator is

$O_{m}=\frac{\sqrt{\lambda}H_{I}}{8\pi^{2}}\left[\frac{d^{2}}{dx^{2}}+\frac{1}{x}\frac{d}{dx}-\frac{m^{2}}{x^{2}}-\frac{8\pi^{2}}{9}(2x^{2}(\pi^{2}(x^{2}-\alpha^{-2})^{2}-3)+3\alpha^{-2})\right]~{}.$
(10)

$O_{m}$ has always one zero eigenvalue, $\Gamma_{10}$, corresponding to the
equilibrium solution of eq. (7). The other eigenvalues (in units of
$\sqrt{\lambda}H_{I}$) are positive and discrete, and will be denoted by
$\Gamma_{nm}$ with $n>1$ integer. $\Gamma_{nm}$ and the corresponding
eigenfunctions $\psi_{nm}$ can be found numerically.444For example, using
Mathematica’s “NDEigensystem”. The most general solution of the Schröedinger-
like equation (7) is then

$\tilde{P}(t,x,\theta)=\sum_{nm}a_{nm}\psi_{nm}(x)\frac{e^{im\theta}}{\sqrt{2\pi}}e^{-\Gamma_{nm}t}~{},$
(11)

where the coefficients $a_{nm}$ are set by the initial conditions and we take
$\psi_{nm}$ normalized in such a way that
$\int_{0}^{\infty}dx\,x\,\psi_{nm}(x)\,\psi_{n^{\prime}m}(x)=\delta_{nn^{\prime}}$.
Using $\psi_{10}\propto e^{-v}$ (which follows from $O_{m}\tilde{P}_{\rm
eq}=0$ and $\tilde{P}_{\rm eq}=e^{-v}$), the probability distribution is
found,

$P(t,x,\theta)=P_{\rm
eq}(x)+\sum_{(n,m)\neq(1,0)}a_{nm}e^{-v(x)}\psi_{nm}(x)\frac{e^{im\theta}}{\sqrt{2\pi}}e^{-\Gamma_{nm}t}~{}.$
(12)

This shows that $\Gamma_{nm}$ determine the evolution of $P$, which tends to
$P_{\rm eq}$ asymptotically. The late-time behavior is regulated by the
smallest non-zero eigenvalue, corresponding to $(n,m)=(1,1)$. In particular,
$1/\Gamma_{11}$ measures the cosmic time on which equilibrium is approached,
i.e. this happens after a ‘critical’ number of e-folds of expansion

$N_{11}\equiv\frac{H_{I}}{\Gamma_{11}}~{}.$ (13)

Figure 1:  Left: The expectation value of the radial field squared
$\left<\rho^{2}\right>/f_{a}^{2}$ on the equilibrium distribution as a
function of $\alpha=H_{I}/(\lambda^{1/4}f_{a})$. This is at the vacuum for
$\alpha\ll 1$ and displaced by $\alpha^{2}$ for $\alpha\gg 1$. Right: The
combination $\lambda^{1/2}H_{I}/\Gamma_{ij}$ corresponding to the two smallest
non-vanishing eigenvalues of $O_{m}$, $\Gamma_{11}$ and $\Gamma_{20}$, as
function of $\alpha$. The related quantity $N_{11}=H_{I}/\Gamma_{11}$ sets the
number of e-folds on which the probability distribution of the complex scalar
$\Phi$ approaches its equilibrium form during inflation.

Figure 1 (right) shows $\sqrt{\lambda}H_{I}/\Gamma_{nm}$ for the first two
non-vanishing eigenvalues as a function of $\alpha$. The smallest eigenvalue
in units of $H_{I}$ increases substantially as $\alpha$ increases for
$\alpha\lesssim 3$ (with $\lambda$ fixed), corresponding to the fact that
quantum fluctuations become larger so equilibrium is reached faster. In the
limit $\alpha\lesssim 1$,
$\lambda^{1/2}\Gamma_{11}/H_{I}\simeq(H_{I}/(9\lambda^{1/4}f_{a}))^{-2}$ so
$\Gamma_{11}$ is independent of $\lambda$. Meanwhile, for $\alpha\gg 1$ the
lowest eigenvalues are

$\Gamma_{11}=0.108\sqrt{\lambda}H_{I}\,,~{}~{}~{}~{}~{}~{}~{}\Gamma_{20}=0.366\sqrt{\lambda}H_{I}\,,~{}~{}~{}~{}~{}~{}~{}\Gamma_{21}=0.611\sqrt{\lambda}H_{I}\,$
(14)

and in this limit

$N_{11}\simeq\frac{10}{\sqrt{\lambda}}~{},\quad\qquad{\rm for}\ \ \alpha\gg
1~{}.$ (15)

Intuitively, for $\alpha\gg 1$ the variance of $P$ initially grows as
$\sigma^{2}=NH_{I}^{2}/(2\pi)^{2}$ so the equilibrium distribution in eq. (5)
takes $N\propto 1/\sqrt{\lambda}$ e-folds to be reached, independent of
$H_{I}/f_{a}$.555When $V$ is negligible a solution of eq. (2) is a Gaussian
with variance $\sigma^{2}=NH_{I}^{2}/(2\pi)^{2}$ with $N=H_{I}t$. Because
$\Gamma_{11}$ is an increasing function of $\alpha$, eq. (15) gives a lower
bound for the number of e-folds to reach equilibrium, with $\rho$ light.

Importantly for what follows, $N_{11}$ is also roughly the number of e-folds
that it takes for $\Phi$ to be fully randomized; more precisely, on scales
much larger than $d_{11}\equiv e^{N_{11}}/H_{I}$ at the end of inflation,
$a(x)$ acquires random values in the interval $[-\pi,\pi]$, while over smaller
regions it has a nearly uniform value. Equivalently, two points in space that
left each other’s horizons fewer than $N_{11}$ e-folds before the end of
inflation have correlated $\Phi$ and points separated by more than $N_{11}$
e-folds are uncorrelated (we will show this more precisely in Section 5.1).
Consequently, the axion field has non-trivial topology on scales larger than
$d_{11}$. In particular if $N_{11}\gtrsim 50-60$ the field is nearly constant
over our present Hubble patch. Meanwhile, in the opposite regime, at the end
of inflation the field was inhomogeneous over what is now the observable
universe, but it can still be approximately constant over smaller macroscopic
regions. In the first case the PQ symmetry is effectively always broken in our
universe while in the second case we can think that the symmetry is restored
on average at scales larger than $d_{11}$. From eq. (15) the latter scenario
is realized for $\lambda\gtrsim 0.05$, when $\alpha\gg 1$.

Figure 2: The axion field $a(x)/f_{a}\equiv{\rm Arg}(\Phi)$ over a spatial
slice after 6 e-folds of inflation, with $H_{I}/f_{a}=\pi$ and the field
initially in the vacuum $\Phi_{0}=f_{a}/\sqrt{2}$. This is obtained by
neglecting the effect of the PQ potential during inflation, i.e. we take the
real and imaginary parts of $\Phi$ as Gaussian random fields with flat power
spectra (this is a good approximation if $\lambda\ll 1$ but not accurate for
large $\lambda$, in which case $\Phi$ is non-Gaussian and the power spectra
have a non-zero slope, but such an approximation is sufficient for
illustrating the broad picture). Over regions separated by one e-fold of
inflation, the angular field is nearly homogeneous, however over regions
separated by $\simeq 4$ or more e-folds the angular field is random.
Consequently, after inflation we expect the resulting string network to reach
approximately one string per Hubble patch once regions that have size
$e^{4}/H_{I}$ at the end of inflation have re-entered the horizon.

To illustrate the possibility of symmetry restoration only on average at
sufficiently large scales, in Figure 2 we show $a(x)/f_{a}$ across a two-
dimensional spatial slice at the end of $6$ e-folds of inflation. The field is
generated neglecting $V$ (sufficient for our present purposes), so $\phi_{1}$
and $\phi_{2}$ fluctuations are added to an initial homogeneous field value,
arbitrarily set to $\phi_{1}=f_{a}$ and $\phi_{2}=0$, as Gaussian fluctuations
with a flat power spectra of amplitude $H_{I}/(2\pi)$ (i.e. root mean square
fluctuations in the region of $\sqrt{6}H_{I}/(2\pi)$). This reproduces the
expectation from inflation with each point corresponding to a region of size
$H_{I}^{-3}$ at the end of inflation. We use $H_{I}/f_{a}=\pi$, which results
in the radial mode spreading over the top of the potential only after a few
e-folds. Consequently, $a(x)$ is approximately homogeneous over regions
corresponding to $\simeq 3$ e-folds of inflation, however over larger regions
it has substantial inhomogeneities, and the PQ symmetry is effectively
restored. Anticipating our analysis in Section 3, in this example we expect
that the string network that forms after inflation will reach one string per
Hubble patch when the modes that left the horizon $\simeq 4$ e-folds before
the end of inflation have re-entered the horizon.

In the rest of the paper we assume that, over a causal patch $N$ e-folds
before the end of inflation, $\Phi$ has a value $\Phi_{0}$ typical of the
equilibrium distribution, i.e.
$P=\delta(\phi_{1}-\sqrt{2}\Phi_{0})\delta(\phi_{2})$. For $\alpha\ll 1$ this
is reasonable, while for $\alpha\gg 1$ this assumes that equilibrium is
reached before visible inflation. In anticipation of what follows, in Figure 3
we show the probability distribution after $N=25$ inflationary e-folds for
different $\alpha=H_{I}/(\lambda^{1/4}f_{a})$ and $\lambda$, starting from
$P=\delta(\phi_{1}-\sqrt{2}\Phi_{0})\delta(\phi_{2})$, with $\Phi_{0}$ given
by the standard deviation on the equilibrium distribution,
$\Phi_{0}=\sqrt{\left<\rho^{2}\right>/2}$. These results can be thought as the
field distribution over a Hubble patch that re-enters the horizon at around
the QCD crossover (similar plots focusing on $\alpha\simeq 1$ were given in
Ref. [45]). For $\lambda=1$ and $H_{I}=f_{a}$, $\alpha=1$ and $N_{11}\simeq
72$: the radial mode is strongly peaked in the vacuum and the equilibrium
distribution is not reached in 25 e-folds, with the field spreading only part
way around the minimum of its potential; for $\lambda=1$ and $H_{I}=2f_{a}$,
$\alpha=2$ and $N_{11}\simeq 15$ so after 25 e-folds $P$ is part-way to its
equilibrium distribution, on which it is non-zero at $\Phi=0$ with
$\left<\rho^{2}\right>\simeq f_{a}^{2}$ still; and for $H_{I}=5f_{a}$,
$\alpha=5$ and $N_{11}=9$ so the equilibrium distribution is approximately
reached and on this $P$ is unsuppressed at $\Phi=0$ and the expectation value
$\left<\rho^{2}\right>\gtrsim f_{a}^{2}$ is slightly displaced. Meanwhile, for
$\lambda=10^{-3}$, $\alpha\gg 1$ and $N_{11}\simeq 300$ for all $H_{I}$ and
$\Phi\simeq 0$ is not explored in 25 e-folds of inflation in any of these
cases because for larger $H_{I}$ the initial displacement gets larger at the
same rate as the spread of the field: given the initially sharply peaked
distribution, $P$ is approximately Gaussian with variance
$25H_{I}^{2}/(4\pi^{2})\simeq 0.6H_{I}^{2}$, with the radial mode initially
displaced from the origin by $\left<\rho^{2}\right>\simeq 7H_{I}^{2}$.

Figure 3: The probability distribution $P(\Phi)$ after $N=25$ e-folds of
inflation, starting from the initial condition
$P=\delta(\phi_{1}-\sqrt{2}\Phi_{0})\delta(\phi_{2})$ setting
$\Phi_{0}=\sqrt{\left<\rho^{2}\right>/2}$, with $\sqrt{\left<\rho^{2}\right>}$
the expectation value on the equilibrium distribution, shown in Figure 1. We
show this distribution for different Hubble scales during inflation $H_{I}$
and, in the upper and lower panels respectively, $\lambda=1$ and
$\lambda=10^{-3}$. At a fixed $\lambda$, $H_{I}/f_{a}$ controls both the
initial field displacement and the spread of the probability distribution. Red
circles indicates the vacuum $|\Phi|=f_{a}/\sqrt{2}$.

###  2.2 Evolution soon after inflation

We now consider the dynamics of $\Phi$ after the end of inflation when
reheating begins. If $N_{11}\gg 1$, $\Phi$ is roughly constant over several
causal patches, so to a good approximation we initially neglect the spatial
gradients and analyse the evolution patch by patch.666We also neglect any
initial velocity of the field. If $\rho\lesssim f_{a}$ (e.g. for
$H_{I}/f_{a}\lesssim 2$ and $\lambda=1$, see Figure 3), then the evolution is
straightforward: in most patches the field simply rolls along the radial
direction into the nearest minimum of $V$ (except, as we discuss later, for
the possibly few patches starting close to $\rho=0$, that will become the
cores of strings). Conversely if $\alpha=H_{I}/(\lambda^{1/4}f_{a})\gg 1$ the
radial mode takes values $\gg f_{a}$. In this case the field can oscillate
along the radial direction over the top of its potential $V$ many times before
settling down into the minimum, potentially dramatically affecting the field
spatial distribution.

Figure 4:  Solid black line: final value of the radial field, $\rho_{f}$, as a
function of the value at the end of inflation, $\rho_{i}$, assuming a matter
dominated cosmological background during the time when the radial mode relaxes
to the minimum of its potential. Negative values here mean that the field ends
up on the opposite side of the potential (e.g. $\theta=\pi$ when starting at
$\theta=0$). Dashed red line: results obtained assuming a background of
radiation domination.

We focus on a simple scenario where inflation ends with the inflaton
oscillating around the minimum of its potential, which we assume is quadratic
so the Universe is matter dominated at this time with scale factor $R\propto
t^{2/3}$. We also assume that the inflaton decay rate is sufficiently small
that the maximal temperature $T_{\rm max}<f_{a}$; that non-perturbative
processes, e.g. preheating, do not occur; and an instantaneous transition from
inflation to matter domination, so the Hubble parameter is $H_{I}$ at the
start of this era. With these assumptions, for the $H_{I}/f_{a}$, $\lambda$ of
interest, the radial mode settles down to the minimum of its potential while
the Universe is still in matter domination. During these dynamics, in each
patch the equation of motion is approximated by

$\ddot{\rho}+\frac{2}{t}\dot{\rho}+\lambda(\rho^{2}-f_{a}^{2})\rho=0~{}.$ (16)

$\lambda$ can be eliminated by defining $\sqrt{\lambda}t=\tau$, and reheating
starts at $\tau_{0}=2\sqrt{\lambda}/(3H_{I})$. The final value $\rho_{f}$ at
the end of the field’s oscillations is plotted in Figure 4 as a function of
its initial value at the end of inflation $\rho_{i}$, where $\rho_{f}=\pm
f_{a}$ depending on this being on the same/opposite side.777In fact, this
function depends on $H/(\lambda^{1/2}f_{a})$, as can be seen from eq. (16) and
the rescaling of $t$, however the dependence is weak and will not affect our
main results in Figure 6. The result of Figure 4 is shown for
$H_{I}/(\lambda^{1/2}f_{a})=30$. We also plot the function that one would
obtain assuming a background radiation domination, in order to show that the
results depend critically on this assumption. The difference is due to the
larger friction term during matter domination.

The oscillation of $\rho$, which we call _overshooting_ , is important because
it can at least partially randomise $\Phi$ even if this was close to uniform
over $e^{N}$ patches immediately at the end of inflation, as e.g. in the
bottom panel of Figure 3 for $N=25$. To see this, consider the values of the
PQ field at two points in space, $x$ and $0$, such that at the end of
inflation $\arg\left(\Phi(x)\right)\simeq\arg\left(\Phi(0)\right)$ and
$|\Phi(x)|-|\Phi(0)|\ll|\Phi(0)|$. In the absence of overshooting, the value
of $\Phi$ at these points would remain strongly correlated. However, if the
small difference between $|\Phi(x)|$ and $|\Phi(0)|$ is enough that the radial
mode settles down on opposite sides of the potential, i.e. if one of the steps
in Figure 4 is crossed, then $\arg(\Phi(0))$ and $\arg(\Phi(x))$ will differ
by $\pi$ once the radial mode has reached its minimum, effectively amplifying
the initially small relative difference in the field values. We will see that
such overshooting can lead to the formation of strings: this is because of
special initial field values – the boundary points – where the field ends up
exactly on the top of the potential, which will become string cores (this
possibility has previously been noted in Ref. [45]).

To illustrate the effect of overshooting, in Figure 5 we show
$a/f_{a}=\arg(\Phi)$ over a two-dimensional spatial slice for a realisation of
the inflationary initial conditions. As in Figure 2, we generate $\phi_{1}$
and $\phi_{2}$ as Gaussian random fields with flat power spectra (i.e. we
neglect $V$) and fluctuations corresponding to 6 e-folds of inflation. We set
$H_{I}/f_{a}=\pi/2$, and unlike Figure 2, we assume $\Phi_{0}=2.3f_{a}$,
corresponding to $\sqrt{\langle\rho^{2}\rangle}$ in eq. (5) for
$\lambda=10^{-2}$. These values are such that $\Phi$ does not fluctuate over
the top of its potential during inflation. As a result, in the left panel,
which shows $\arg(\Phi)$ immediately at the end of inflation, only small
inhomogeneities are present. However, in the right panel we plot $\arg(\Phi)$
after the evolution to the minimum and overshooting, still assuming
$\lambda=10^{-2}$. Evidently, the angular field acquires large fluctuations on
spatial scales corresponding to $4$ or more e-folds of inflation. These will
be sufficient for an (initially underdense) string network to form even though
over the majority of space $\arg(\Phi)$ is close to either $0$ or $\pi$, with
relatively few points in-between.

Figure 5:  Left: Example of $a/f_{a}\equiv{\rm Arg}(\Phi)$ over a spatial
slice at the end of 6 e-folds of inflation, similar to Figure 2, except for
$H_{I}/f_{a}=\pi/2$ and with the crucial difference that the initial value of
the PQ field is taken to be $\Phi_{0}\simeq 2.3f_{a}$. We generate the
fluctuations of $\phi_{1}$ and $\phi_{2}$ neglecting the effect of the
potential $V$, but note that the $\Phi_{0}$ we choose corresponds to
$\lambda=10^{-2}$ assuming $|\Phi_{0}|$ given by the expectation value in the
equilibrium distribution eq. (5) ($\lambda$ is sufficiently small that the
effect of the potential during inflation is indeed negligible). Only a small
portion of the field has spread over the top of the potential and there are
only relatively small inhomogeneities just after inflation. Right: The same
field configuration as in the left panel, but after accounting for
overshooting, i.e. once the radial mode has settled down to the minimum of its
potential, by solving eq. (16) (neglecting gradients and assuming an
instantaneous transition to matter domination). The radial mode oscillates
over the top of its potential and the initially relatively small
inhomogeneities can change which side of its potential it finishes on. As a
result, the angular field is efficiently randomised. This is efficient enough
that strings form, and subsequently reach ${\cal O}(1)$ string per Hubble
volume when $\simeq 5$ e-folds of inflation re-enter the horizon.

##  3 Formation of strings

As the radial mode settles down to the minimum throughout almost all of space,
cores of strings form from points where the radial mode is at the top of its
potential with non-trivial winding around this point.888We do not consider the
regime in which the radial mode is very heavy but the axion wraps its
fundamental domain many times, in which case the radial mode can be forced to
the top of its potential nucleating strings (as can occur for vector dark
matter produced by inflationary fluctuations) [67]. We define the string
density $\xi$ at a given cosmological time $t$ after the end of inflation as

$\xi=\lim_{{\cal V}\rightarrow\infty}\frac{\ell({\cal V})t^{2}}{{\cal V}}~{},$
(17)

where $\ell({\cal V})$ is the length of strings in a volume ${\cal V}$, such
that $\xi$ measures the string length in Hubble lengths per Hubble volumes. On
the attractor scaling solution $\xi$ is order-one, up to corrections
$\propto\log(m_{\rho}/H)$ [68, 69, 70, 71, 31], while for an underdense
network $\xi\ll 1$. Our key aim is to predict how many e-folds worth of
inflationary fluctuations must re-enter the horizon, which we refer to as $N$
e-folds of inflation re-entering the horizon, for $\xi\simeq 1$ to be reached;
we say that at this time the string network has re-entered the horizon.

We expect that the value of $\xi$, at the time when $N$ e-folds of inflation
have re-entered the horizon, is related to what extent $\Phi$ has spread over
the top of the potential during the last $N$ e-folds of inflation.
Intuitively, if $\Phi$ is completely spread over the top of the potential, as
is the case if the equilibrium distribution is reached, over a given Hubble
patch we expect at least one string. Meanwhile, if the probability
distribution of $\Phi$ is mostly concentrated on one side of the potential, we
expect $\xi\ll 1$, with strings only in a few rare Hubble patches. In the
latter case, $\xi$ will subsequently increase as more e-folds of inflation re-
enter the horizon, until $\xi\simeq 1$ is reached once the critical number of
e-folds have re-entered.

In order to quantify the string density $\xi$ after $N$ e-folds of inflation
have re-entered the horizon, we define the fraction $F$ of the field that has
spread to the other side of the potential as

$F(N)\equiv\frac{1}{2}\left[1-{\rm
Abs}\int_{-\infty}^{\infty}d\phi_{1}\int_{-\infty}^{\infty}d\phi_{2}~{}{\rm
sgn}\left(\phi_{1}\right)~{}B\left(\sqrt{\phi_{1}^{2}+\phi_{2}^{2}}/f_{a}\right)~{}P(\phi_{1},\phi_{2},N)\right]~{}.$
(18)

Here $P$ is the solution the FP equation after $N$ inflationary e-folds
starting from the initial condition
$P(\phi_{1},\phi_{2})=\delta(\phi_{1}-\sqrt{2}|\Phi_{0}|)\delta(\phi_{2})$.999In
a cosmological setting before the final $N$ e-folds of inflation there were
$\simeq 60-N$ visible e-folds of inflation. Consequently, we should account
for a distribution of initial conditions. For most of our subsequent work
however, it suffices to take $|\Phi_{0}|$ to have a single fixed value. In eq.
(18), $B$ accounts for overshoot and is the step-like function plotted in
Figure 4 that takes values $\pm 1$ depending on whether the radial mode
finishes oscillating on the same or the opposite side to it started. $F(N)$
measures the fraction of the field on the less populated side of the
potential. $F=0.5$ corresponds to the field having spread completely over the
potential while $F=0$ corresponds to none of the field having crossed (see
also [72, 73] for related analysis in the context of biased initial
conditions). Note that for $\alpha\ll 1$, the radial mode remains in the
minimum of its potential and $F=0.5$ could be reached by fluctuations only in
the angular direction that would not lead to strings, however this regime is
not relevant to any $H_{I}/f_{a}$, $\lambda$ of phenomenological
interest.101010More precisely, on the equilibrium distribution, the
probability of having points at $\Phi=0$, i.e.
$P(0)\propto\exp(-\frac{2\pi^{2}}{3}\alpha^{-4})$, must be non-negligible, and
provided $\alpha\gtrsim 1.6$ the radial mode does indeed explore the top of
its potential. For smaller values, say $\alpha\lesssim 1$, it might be that
strings do not form even at equilibrium (giving only domain walls as in Ref.
[43]), but this is in the $N_{11}>50$ region of our parameter space anyway,
see eq. (13) and Figure 1. A benefit of using the measure $F$ (as opposed to
other possibilities such as $N_{11}$ defined in eq. (13)) is that $F$ captures
also the overshoot mechanism: $F\simeq 0.5$ can occur either because $P$ has
spread over the top of the potential during the last $N$ e-folds of inflation
(e.g. as in the top-right panel of Figure 3) _or_ because the field has
settled evenly on the two sides of the potential after overshooting despite
being almost homogeneous at the end of inflation (e.g. in the bottom-right
panel), and in both cases we expect a string network to form.

In Section 4 we will present evidence from numerical simulations that $\xi$ at
the time when $N$ e-folds have re-entered the horizon is indeed related to
$F(N)$, for the $N\lesssim 8$ that can be simulated. In particular, we will
see that $\xi\simeq 1$ is obtained for $F\simeq 0.3\div 0.45$. We therefore
postulate that $\xi\simeq 1$ is obtained when $F(N)\simeq 0.4$ even when this
is reached for much larger $N\lesssim 60$, as is the case for a theory for
which a string network re-enters the horizon soon before or even after the QCD
crossover, up to the present day. For a given $H_{I}/f_{a}$, $\lambda$, we
define $N_{s}$ to be the critical number of e-folds for $F$ to reach some
value, which we usually choose to be $0.35$ (analogous results assuming
$F=0.4$ are given in Appendix B to give an indication of the uncertainties).

###  3.1 Inflationary formation

Figure 6:  Left: Estimates of the number of e-folds of inflation $N_{s}$ that
must re-enter the horizon for a network of strings formed from inflationary
fluctuations to reach an order-one number per Hubble patch. We assume $F=0.35$
(see Appendix B for the analogous plot with $F=0.4$). The $N_{s}<25$ region
corresponds to the resulting string network re-entering the horizon before the
QCD crossover. In the dark grey region in the bottom-right corner the radial
mode is heavy, $m_{\rho}>\sqrt{2}H_{I}$, and our discussion does not apply
(the usual isocurvature bounds on $H_{I}/f_{a}$ apply here). Right: Analogous
plot but for string formation via overshoot, again with $F=0.35$. The detailed
shape of the region labelled “$N_{s}>25?$” is particularly sensitive to the
assumed value of $\Phi$ before the visible era of inflation, the assumed
critical $F$ and the cosmological history soon after the end of inflation, so
this is especially uncertain. Note that in reality both production mechanisms
apply, and the split here is to show the parameter space in which each is
dominant. In general, one would compute $F(N)$ defined in eq. (18), which
incorporates both effects, as shown in Figure 8.

We start with ‘inflationary formation’ of strings, occurring when
topologically nontrivial configurations are produced directly by the radial
mode spreading over its potential during inflation, such that strings form
immediately after the end of inflation, in a time of order $1/m_{\rho}$. This
amounts to neglecting overshooting by setting $B=1$ in eq. (18).

To calculate the fraction $F$, we use eq. (18) and solve the FP equation for
$P$. Assuming that the field is initially at $\phi_{1}=\sqrt{2}|\Phi_{0}|$ and
$\phi_{2}=0$, corresponding to $x_{0}=\lambda^{1/4}\sqrt{2}|\Phi_{0}|/H_{I}$
with $\Phi_{0}$ sampled from the equilibrium distribution, we get that the
coefficients for $P$ in eq. (11) are

$a_{nm}=\frac{1}{\sqrt{2\pi}}e^{v(x_{0})}\psi_{nm}(x_{0})~{}.$ (19)

We can now ask how the probability distribution evolves in field space over
time. The fraction $F$ is simply obtained by integrating $P$ over
$\phi_{1}<0$. The $m=0$ wave-functions on the right hand side of eq. (11) do
not contribute to $F$ in eq. (18) given their orthonormality. For $F$
sufficiently close to $1/2$, the lowest eigenvalue ($m=1$) dominates and we
find

$F=1/2+a_{11}b_{11}e^{-N\Gamma_{11}/H_{I}}\,,~{}~{}~{}~{}~{}~{}b_{11}=-\sqrt{\frac{8}{\pi}}\int
dx\,xe^{-v(x)}\psi_{11}(x)\,.$ (20)

This means that the number of inflationary e-folds required for a fraction $F$
to be reached is

$N(F)=\frac{H_{I}}{\Gamma_{11}}\log\left[\frac{2a_{11}b_{11}}{2F-1}\right]\,.$
(21)

Similarly to $N_{11}$, $N(F)$ is proportional to $H_{I}/\Gamma_{11}$ and thus
is controlled by the lowest eigenvalue $\Gamma_{11}$, even though the
coefficient is sensitive to the fraction $F$. Our expectation that $F$
determines the string density and the result $N(F)\simeq N_{11}$ for $F\simeq
0.5$, are consistent with the intuition in Section 2.1 that a string network
reaches $\xi\approx 1$ once, approximately, $N_{11}$ e-folds of inflation have
re-entered the horizon. From eq. (21), $F$ usefully parameterizes the
plausible range of the numerical coefficient, the true values of which can
only be reliably estimated through numerical simulations.111111Alternatively,
we could require that in a certain fraction of space the radial field is
inside a circle of radius $|\rho|<rf_{a}$ around the origin, with $r<1$. Due
to the angular dependence in this case only $m=0$ wave-functions contribute,
so that the second wave-function with $m=0$ would determine $N_{s}$. This
criterion however will depend both on the radius $r$ chosen and the fraction.
Using this measure we define a critical number of e-folds $N_{s}$ for the
network to form, using the threshold $F(N_{s})=0.35$, i.e. assuming that 35%
of points reaches the other side of the potential with respect to $\Phi_{0}$.
In general it differs by $N_{11}$ by a coefficient of about 2.

In Figure 6 (left), we show $H_{I}/f_{a}$ and $\lambda$ for which the network
to re-enters the horizon at the critical number of e-folds $N_{s}$. We show
contours for $N_{s}=25,50$ corresponding to a string network re-entering the
horizon just before the QCD crossover and shortly before the present day
respectively; using $N_{s}=N_{11}=H_{I}/\Gamma_{11}$ gives similar results.
From Figure 1 (left), for $\alpha\lesssim 2$, corresponding to the region
$H_{I}/f_{a}\lesssim 2\lambda^{1/4}$ below the dotted line, the equilibrium
distribution is peaked at the radial mode’s vacuum and in this regime
$\Gamma_{11}/H_{I}$ is approximately independent of $\lambda$. Thus, the
contours are horizontal at $H_{I}/f_{a}\simeq 9/\sqrt{N_{s}}$ and larger
$H_{I}/f_{a}$ simply results in $\Phi$ reaching equilibrium after a smaller
number of e-folds. Conversely, for $\alpha\gg 2$,
$\lambda^{1/2}H_{I}/\Gamma_{11}$ does not change as $H_{I}/f_{a}$ increases at
fixed $\lambda$, because the radial mode is displaced at
$\langle\rho^{2}\rangle\gg f_{a}^{2}$, see Figures 1,3. This leads to vertical
contours at $\lambda\simeq 0.2(N_{s}/25)^{-2}$. As we will see next, in this
region of parameter space the prediction from inflationary fluctuations is
superseded by the effects of overshoot.

###  3.2 Overshoot mechanism

If $P$ does not approach the equilibrium distribution within the visible
number of e-folds of the Universe (e.g. $\lambda\lesssim 0.05$, for $H_{I}\gg
f_{a}$), naively no string network would form.121212See Chapter 14 of Ref.
[74] for a detailed discussion. However, as discussed in Section 2.2, a string
network might still form after inflation by the overshoot mechanism. Yet, it
is not a priori clear if this could lead to the standard scaling solution, or
if it inevitably produces either an underdense network or a collection of
small loops that do not percolate and form a network of long strings. The
results from simulations that we present subsequently show that for $F\simeq
0.35\div 0.45$, a string network with $\xi\simeq 1$ is indeed reached.

In Figure 6 (right) we show contours of $N_{s}$, again defined by
$F(N_{s})=0.35$ as function of $H_{I}/f_{a}$ and $\lambda$ considering only
overshoot. This means that we remove the sign function from the definition of
$F$ in eq. (18), so that it is only the overshoot function that results in
$F\neq 0$. We use the overshoot function in Figure 4 and, as in the case of
inflationary production, we set as initial condition $\Phi_{0}$ sampled from
the equilibrium distribution. In the time-evolution of eq. (18), $P$ is
approximately Gaussian in the parameter space for which overshoot is
important, i.e. the equilibrium distribution is not reached during visible
inflation, with variance $N(H_{I}/2\pi)^{2}$. Additionally, since for the
relevant $H_{I}/f_{a}$ and $\lambda$ the field does not reach the asymptotic
distribution during visible inflation, other initial conditions (i.e.
$\Phi_{0}$) are possible but perhaps less plausible.

At a fixed $\lambda$, as $H_{I}/f_{a}$ increases the spread of $P$ over $N$
e-folds of inflation increases, and the “steps” of $B$ in Figure 4 do not get
much wider at larger $\Phi_{0}$; thus the angular mode is increasingly
randomised. The contours do however show a complex shape, labelled
‘$N_{s}>25?$’, which is associated to the details of the step function
determining overshoot. Initial field values close to one of the steps can
easily lead to strings because only a small change is needed for the radial
mode to end on the other side, whereas if the initial field value is in the
middle of one of the steps it is relatively harder to form strings. The latter
leads to regions parallel to the line $H/f_{a}\propto\lambda^{1/4}$ (since
$H/\lambda^{1/4}$ controls the equilibrium field displacement used in
$\Phi_{0}$) that correspond to the steps of $B$. The regions end at large
enough $H_{I}/f_{a}$, because at large $H_{I}$ the field’s probability
distribution is always spread enough that it does not fully fit in a single
step of $B$. We stress that the shape of these complex regions depends in
detail on the initial field value; the critical value of $F$, which we do not
know precisely; and the details of the overshoot function assumed, which
depends on the expansion history immediately after inflation (in making this
plot we assume that of Figure 4).131313Additionally, in calculating $F(N)$ to
make this plot, we approximate that the radial mode $\rho$ always has value
$\rho_{0}$ $N$ efolds before the end of inflation. In reality however, over
what the now the observable Universe $\rho$ has some distribution of values
due to the preceding efolds of inflation. If taken into account, this would
smear our some of the sharp features in the Figure. Further details can be
found in Appendix C.1. As a result, these regions of the plot should be taken
as a indication of the values of $H_{I}$, $\lambda$ at which such, complex,
dynamics might occur. The result of overshoot is that string can form with
relatively small $N_{s}$ also in the region $H_{I}/f_{a}\gtrsim
2\lambda^{1/4}$, i.e. above the dotted line, where inflationary production is
inefficient, see Figure 8 (left).

Finally, we note that Ref. [75] finds that if the PQ field reaches values
$|\Phi|\gg 10^{4}f_{a}$ strings can form after inflation as a result of
parametric resonance (see also [76, 77, 78, 79, 80, 81, 82, 83, 84]). Such
field values can occur with $\Phi$ roughly at its equilibrium value eq. (5)
for $\lambda\ll 10^{-2}$ (much smaller than those we consider), or it could
simply be assumed that at the start of the visible era of inflation $\Phi$
took a much larger field value than eq. (5). Similarly to the overshoot
mechanism, parametric resonance leads to strings as a result of the radial
mode oscillating over the top of its potential immediately after inflation.
However, in parametric resonance it is the small fluctuations in the direction
orthogonal to $\Phi$’s large oscillations that are dominantly amplified. Such
dynamics are relevant in complementary parts of parameter space to those we
focus on.

##  4 Comparison to simulations

We now present results from numerical simulations of the string network. We
stress that the purpose of these is to test the general picture described
above and demonstrate that the fraction $F$ defined in eq. (18) sets the
string density of an underdense network, at least approximately. However, as
we discuss shortly we cannot directly study the regime of interest, and we
also modify the physical system in some ways that deviate from the true
dynamics.

In our simulations, the field $\Phi$ is discretized on a lattice with $n^{3}$
grid points. We set the field configuration in the initial conditions to
reproduce that from inflation, by generating a realisation of the Langevin
equation [85]. To do so, we follow a method described in [45], see Appendix C.
At the end of this, each lattice point corresponds to a single region of size
$H_{I}^{-3}$ at the end of inflation. Given that we can only simulate grids
with $n^{3}\sim 2048^{3}$ lattice points, our entire simulation volume
corresponds to $N\sim 8$ e-folds of inflation, much smaller than the regime
$25\lesssim N\lesssim 60$ of our interest.

We subsequently evolve the field forward in time with the discretized
equations of motion of the Lagrangian in eq. (1), using a standard leap-frog
algorithm. We start our simulations at $H=m_{\rho}$ for all the values of
$H_{I}$ and $\lambda$ used in generating the initial conditions. Moreover, we
evolve the mass of the radial mode throughout the simulations as
$m_{\rho}\propto 1/R(t)$ where $R(t)$ is the scale factor, by using a time
dependent $\lambda$. This is known as the “fat-string” trick, which means the
string core size grows in time so that it contains a constant number of
lattice points throughout the simulation (in practice $m_{\rho}\Delta=1$ is
sufficient to resolve the string cores where $\Delta$ is the physical lattice
spacing). Although the equations of motion of fat string system differ from
those of the physical system by order-1 amounts, many features of the two
networks agree. As discussed in Appendix C, these choices maximise the number
of e-folds worth of modes from inflation that we can observe re-entering the
horizon. We evolve the system in radiation domination, $H=1/(2t)$, as is
appropriate to the time when the string network re-enters the horizon in the
late-string regime. Such an expansion history is unlikely to be realistic at
the time when radial mode relaxes to the minimum of its potential after
inflation. However, our aim is to test whether, in the $\xi\ll 1$ (under-
dense) regime, $\xi$ depends on $F$ and this is not expected to be sensitive
to the details of the cosmological expansion (we cannot exclude that the exact
value of $F$ that leads to a given $\xi$ might have some dependence, but this
is beyond the precision we can aim for).

As the system is evolved, a string network automatically forms and interacts
(with strings first being cleanly defined objects at $m_{\rho}(t)/H(t)\simeq
10$). We calculate the length of strings in the simulations, and
correspondingly $\xi$ as defined in eq. (17), using a standard algorithm [86].
The simulations can be run until $HL\simeq 2$, where $L$ is the physical box
size, at which point finite volume systematic errors become significant.

We extract $\xi$ at regular intervals throughout each simulation. For each
such “time-shot” we consider the number of e-folds of inflation $N$ that have
re-entered the horizon at such a time. Then we calculate the corresponding
value of $F$ for each set of initial conditions by numerically solving (with
Mathematica’s “NDSolve”) the FP equation, starting from
$P(\phi_{1},\phi_{2})=\delta(\phi_{1}-\sqrt{2}|\Phi_{0}|)\delta(\phi_{2})$,
for $N$ e-folds (i.e. until $t=N/H_{I}$) and finally inserting the resulting
probability $P$ into eq. (18), including the function $B$ that accounts for
overshoot.141414We could also have extracted $F$ from the initial conditions
of the simulation, which would be equivalent in the limit of large statistics.
As expected, the calculated $F$ increases as each simulation progresses and
more fluctuations re-enter the horizon, and for a given simulation time $F$ is
larger when $H_{I}/f_{a}$ used to generate the initial conditions is larger
(keeping $\lambda$ constant).

###  4.1 Results

Figure 7:  The axion string density $\xi$ in simulations starting from initial
conditions corresponding to inflationary fluctuations with different value of
the Hubble parameter during inflation, $H_{I}$, the quartic coupling,
$\lambda$, and radial mode field value at the start of inflation, $\rho_{0}$.
The results are shown as a function of $F$, defined in eq. (18), which
measures the size of the inflationary fluctuations that have reentered the
horizon at the time when $\xi$ is measured. On the left, data with initial
conditions such that the radial mode stays mostly in the vacuum and strings
form from inflationary fluctuations over the top of the potential. On the
right, initial conditions such that the radial mode does not spread over the
top of the potential during inflation, but is displaced from the minimum and
subsequently overshoots the top of its potential leading to strings forming.
Here $\lambda=0$ is the value used to generate the initial conditions, but the
field was subsequently evolved with non-zero $\lambda$. In both cases, $\xi$
is approximately determined by $F$ regardless of the value of $H_{I}/f_{a}$
and $\lambda$.

We carry out simulations for a range of initial conditions (i.e.
$H_{I}/f_{a}$, $\lambda$, and value of $\Phi_{0}$ at the start of the $\simeq
8$ simulated e-folds of inflation) leading to string networks that are
underdense compared to the attractor, scaling, solution. The main result from
these is presented in Figure 7, where we plot the string density $\xi$ as a
function of the value of $F$ for data sets with differing initial conditions.

The results in the left panel correspond to initial conditions such that $F$
is non-zero mostly because of fluctuations of the radial mode over the top of
its potential during inflation, i.e. there is very little overshooting after
inflation (in particular, $\Phi_{0}=f_{a}/\sqrt{2}$). The key, non-trivial,
point is that data from different sets of initial conditions, which lead to
different relations between $F$ and the simulation time, approximately
overlap. In other words, the values of $\xi$ at $F$ corresponding to the final
time of some of the simulations match those at early simulation times in
others. This suggests that $\xi$ in the underdense regime is indeed well-
predicted by $F$.

Meanwhile, the results in the right panel correspond to initial conditions
with relatively large initial radial mode value,
$\Phi_{0}=\\{3f_{a}/\sqrt{2},5f_{a}/\sqrt{2}\\}$, such that the field is not
spread over the top of the potential during inflation and strings can only
form by the overshoot mechanism. (In such simulations, to a good approximation
we neglect the potential $V$ in generating the initial conditions, i.e. we set
$\lambda=0$, but perform the subsequent evolution after inflation with non-
zero $\lambda$.) These results confirm that the radial mode overshooting
indeed leads to the formation of a string network. Additionally, there is
again a clear trend that larger $F$ leads to more strings although the
relation is somewhat less sharp than in the left panel. We note that the
values of $F$ that correspond to a particular $\xi$ are slightly larger than
in the case of inflationary fluctuations over the potential.

We also note that for all the initial conditions plotted, the growth of $\xi$
as a function of Hubble parameter in simulations is $\xi=C(f_{a}/H)^{n}$ where
the power of $n$ is roughly in the range $1/4$ to $1/3$, universally across
different initial conditions and $C$ varies (see Appendix C and Figure 11 for
details). Such a growth is consistent with an underdense string network re-
entering the horizon (such an increase is distinct from the
$\xi\propto\log(m_{\rho}/H)$ growth on the attractor solution). The value of
the power $n$ is smaller than that for a network of long, non-interacting,
strings which would grow as $\xi\propto 1/R^{2}\propto f_{a}/H$ in radiation
domination (as in our simulations). This is expected, because as well as long
strings, small loops are produced and these annihilate as they enter the
horizon. We have extracted the distribution of loop lengths in our underdense
networks and confirmed that such networks do indeed contain many small loops
compared to the scaling distribution. Long strings are relatively rare, but it
is these that ultimately fix the value of $F$ at which $\xi\simeq 1$ and
scaling is reached.151515Our overshoot simulations with the largest $F$ at the
final simulation time (at $\log(m_{\rho}/H)\simeq 7$) reach $\xi\simeq 0.6$,
which is not quite at the scaling value, $\xi\simeq 0.8$ at this
$\log(m_{\rho}/H)$. However, $\xi$ is still increasing with cosmic time, and
we expect scaling to be reached not much beyond the reach of such simulations.

Of course, there are limitations to our results. For each $\lambda$, there is
only a relatively small range of $H_{I}/f_{a}$ that lead to a network that is
both underdense but dense enough to get reasonable statistics for $\xi$, given
our available simulation volume. Related to this, as mentioned, a key
assumption of our work is that the approximate value of $F$ that leads to
$\xi\simeq 0.5$ after $4\div 7$ e-folds of inflation re-enter the horizon is
similar to when $25-60$ e-folds re-enter. Despite being plausible from our
results, it would be interesting to test this extrapolation in more detail,
although unfortunately exponentially larger lattice sizes would be required
for substantial progress.

##  5 Isocurvature perturbations

The inflationary fluctuations of the complex field $\Phi$ discussed in the
previous sections are uncorrelated with the inflaton fluctuations.
Consequently, if they ultimately lead to different axion relic abundances in
different regions, they result in dark matter isocurvature perturbations. As
we review below, if a network of cosmic strings is formed at early times,
reaching a scaling solution, the fluctuations of $\Phi$ on cosmological scales
might be expected to be erased during the post-inflationary evolution of the
field. According to this lore the standard post-inflationary scenario where
the symmetry is restored during inflation is believed not to be affected by
the corresponding isocurvature constraints. In contrast, in the region of
parameters that leads to the late decay of the string-wall system (i.e. after
the QCD crossover), the initial fluctuations might survive and result in a
contribution to dark matter isocurvature.

We define the axion overdensity field
$\delta_{a}(\vec{x})\equiv(\rho_{a}(\vec{x})-\bar{\rho}_{a})/\bar{\rho}_{a}$
where
$\rho_{a}=\frac{1}{2}\dot{a}^{2}+\frac{1}{2}|\nabla{a}|^{2}+\frac{1}{2}m_{a}^{2}a^{2}$
is the axion dark matter energy density once any strings present have been
destroyed and the axion field is non-relativistic; $\bar{\rho}_{a}$ is its
spatial average. The power spectrum $\Delta_{\delta}^{2}(k)$ of the axion
overdensity field is defined by

$\langle\tilde{\delta}_{\delta}^{*}(\mathbf{k})\tilde{\delta}_{a}(\mathbf{k}^{\prime})\rangle=\frac{2\pi^{2}}{k^{3}}\Delta^{2}_{\delta}(k)\delta^{3}(\mathbf{k}-\mathbf{k}^{\prime})\
.$ (22)

where $\tilde{\delta}_{a}$ is the Fourier transform of $\delta_{a}$. We define
the component of $\Delta_{\delta}^{2}$ that comprises isocurvature
fluctuations to be $\Delta_{\rm iso}^{2}$. Meanwhile, the component of
$\Delta_{\delta}^{2}$ consisting of the almost scale-invariant spectrum of
curvature perturbations that is supposed to be originated by the inflaton
fluctuations is $\Delta_{\mathcal{R}}^{2}(k)=A_{s}(k/k_{\rm CMB})^{n_{s}-1}$
where $A_{s}=2.1\times 10^{-9}$ and $k_{\rm CMB}=0.05~{}{\rm Mpc}^{-1}$.

There are strong observational constraints on isocurvature fluctuations on
cosmological scales. Planck cosmic microwave background (CMB) data [87] bound
the fraction of dark matter isocurvature fluctuations relative to the
curvature perturbations. For a $\Delta_{\rm iso}^{2}(k)$ that is approximately
constant in the interval $k=[0.002,0.1]{\rm Mpc}^{-1}$ one finds $r_{\rm
iso}\equiv\Delta_{\rm iso}^{2}(k_{\rm CMB})/\Delta_{\mathcal{R}}^{2}(k_{\rm
CMB})\lesssim 0.04$.

###  5.1 Inflationary Correlations

We start by analysing the super-horizon correlations in the axion field
immediately after the end of inflation. We first consider the regime
(following Chapter 14 of [74]) $H_{I}/f_{a}\gtrsim\lambda^{1/4}$ and
sufficiently small $\lambda$, such that the radial mode is displaced to field
values $\gg f_{a}$, and that the equilibrium solution is not reached during
visible inflation. As discussed, for such theories the value of the complex
field $\Phi_{0}$ prior to the visible e-folds of inflation cannot be predicted
but we continue to make the reasonable assumption that
$\rho_{0}^{2}\simeq\left<\rho^{2}\right>\sim H^{2}_{I}/\sqrt{\lambda}$, where
$\left<\rho^{2}\right>$ is the expectation value on the equilibrium
distribution, eq. (5). In this case after the visible e-folds of inflation,
$\Phi$ has a nearly constant classical value
$\Phi_{0}=\rho_{0}e^{i\theta_{0}}/\sqrt{2}$ in our Hubble patch (more
generally, this is the case unless $\left|\Phi_{0}\right|^{2}\lesssim
50\left(H_{I}/(2\pi)\right)^{2}\ll\left<\rho^{2}\right>/2$).

One finds that immediately after inflation [47],

$\left<\delta\theta^{2}(\mathbf{x})\right>=\frac{NH_{I}^{2}}{(2\pi)^{2}\langle|\Phi_{0}|^{2}\rangle}~{},$
(23)

where $\theta=\rm{Arg}(\Phi)$, $\delta\theta=\theta-\bar{\theta}$,
$\bar{\theta}$ is the spatial average over the entire visible universe and
$N\approx 50$ is the total number of e-folds of visible inflation. If, on
large scales, the correlation in eq. (23) remains unchanged from the end of
inflation to the QCD scale, which neglects the effect of strings and the
relaxation of $\Phi$ to the minimum via overshoot, then eq. (23) results in
isocurvature perturbations in a possible contribution to the axion relic
abundance produced via misalignment. Approximating
$\rho_{a}(\mathbf{x})\propto\theta^{2}(\mathbf{x})$, appropriate for $\theta$
well away from $\pi$, and (for the moment) considering only misalignment
production of axions, for small $\lambda$ we obtain that the isocurvature part
of the power spectrum is approximately flat with an amplitude (at the end of
inflation) [88, 89]

$\Delta_{\rm
iso}^{2}=\frac{H_{I}^{2}}{\pi^{2}\langle|\Phi|^{2}\rangle\theta_{0}^{2}}\simeq\frac{\sqrt{\lambda}}{(\pi\theta_{0})^{2}}~{},$
(24)

where in the last equality we used eq. (5). Comparing with isocurvature
constraints, naively this would exclude such a scenario unless $\lambda$ is
exceedingly small, $\lambda\lesssim 10^{-19}$. However, as we discuss shortly,
both overshoot and the string network that can form as a result of this can
modify this conclusion.

Next we study the case that $H_{I}/f_{a}$ and $\lambda$ are such that at the
end of visible inflation the probability distribution of $\Phi$ is
approximately given by the equilibrium solution, i.e. $N_{11}\lesssim 50$. In
this regime, we can compute correlation functions using the stochastic
formalism of inflation [42]. We are interested in the correlation of the phase
of $\Phi$. Following [42, 66], in the limit $\alpha=H/(\lambda^{1/4}f_{a})\gg
1$ one finds at the end of inflation,161616Correlation functions of arbitrary
functions of the fields can be computed in a similar fashion. The tilt is
determined by the lowest eigenvalue that contributes to the correlation
function and the overlap of wave-functions.

$\langle\cos\left(\theta(\mathbf{x})\right)\cos\left(\theta(0)\right)\rangle\simeq\frac{\kappa_{11}^{2}}{|R_{e}H_{I}\mathbf{x}|^{\frac{2\Gamma_{11}}{H_{I}}}}~{},$
(25)

where $R_{e}$ is the scale factor at the end of inflation and

$\kappa_{11}=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta\cos\theta
e^{i\theta}\int_{0}^{\infty}dx\,x\,\psi_{11}(x)\,\psi_{10}(x)\simeq 0.48~{}.$
(26)

As in the previous case, if we assume that the initial fluctuations in eq.
(25) remain unchanged until late times, they could lead to isocurvature in a
misalignment-like contribution to the axion dark matter abundance. Under this
assumption, we estimate $\Delta_{\rm iso}^{2}$ considering only misalignment
production and approximating that the axion abundance $\rho_{a}\propto
1-\cos\theta$ (which is accurate for $\theta\ll\pi$, has the correct periodic
behaviour, and is easily evaluated). This results in the overdensity power
spectrum

$\Delta_{\rm
iso}^{2}(k)\simeq\kappa^{2}\beta\left(\frac{k}{R_{e}H_{I}}\right)^{\beta}\,,~{}~{}~{}~{}~{}~{}\beta=\frac{2\Gamma_{11}}{H_{I}}~{}\,.$
(27)

In the limit $H_{I}\gg f_{a}$ (with $\lambda\gtrsim 0.01$ so that
$N_{11}\lesssim 50$, c.f. eq. (15)), using eq. (14) the spectrum tilt is
$\beta=0.22\sqrt{\lambda}$. This can be compared with the approximate result
in the mean-field approximation, described in Appendix A, eq. (51), which
numerically gives $\alpha=0.25\sqrt{\lambda}$. The isocurvature power spectrum
in eq. (27) at CMB scales can be written as

$\Delta_{\rm iso}^{2}(k_{\rm CMB})\sim e^{-2N_{\rm vis}/N_{11}}~{},\qquad
N_{\rm vis}\equiv\log(R_{e}H_{I}/k_{\rm CMB})~{},$ (28)

where $N_{\rm vis}$ is the number of e-folds of visible inflation. This
implies that, for $N_{\rm vis}\simeq 50$, under the preceding assumptions the
Planck constraint [87] would be satisfied only for $N_{11}<5$, which in turn
would translate to $\lambda\gtrsim 1$ from eq. (15). However, as we have seen,
for $N_{s}\lesssim 25$ a string network re-enters the horizon and reaches the
scaling regime prior to the QCD crossover, and the effects of these strings
must be taken into account. Conversely, we will argue in Section 6 that in the
regime with strings reentering the horizon after the QCD crossover a
contribution from misalignment as considered in this subsection will be
important.

###  5.2 Effects of strings on correlations

As mentioned, the large-scale correlations in the phase of $\Phi$ induced by
inflation, see eqs. (23) and (25), are not necessarily expected to lead to
significant isocurvature perturbations in the late-time dark matter relic
density. In this subsection we consider how the initial fluctuations might be
suppressed or erased by the effects of a string network (which might re-enter
the horizon before or after the QCD crossover), and subsequently in Section
5.3 we analyse the impact of overshoot.

The first possible effect occurs because a string/domain wall network emits a
significant abundance of dark matter axions. If the network reaches $\xi\simeq
1$ before the QCD crossover, axions are continuously emitted from the network
to sustain the subsequent scaling regime.171717Results from simulations [30,
31] suggest that the abundance from this source is substantially larger than
the would-be abundance from misalignment. Provided the scaling regime is
indeed an attractor that does not depend on the initial conditions, the energy
densities produced by strings (and, once the axion mass is relevant, domain
walls) in different patches should not be correlated, as in the scenario where
the PQ symmetry is restored due to thermal effects. Moreover, in the late
string regime, as we analyse in Section 6, the string/domain wall network that
re-enters the horizon after the QCD crossover gives a substantial contribution
to the axion dark matter abundance; we will assume that this component carries
no isocurvature (although this also deserves a deeper analysis that we leave
for future work). As a result, even if there is a component of the axion
abundance that arises from misalignment and retains long-distance correlations
from inflation, the resulting isocurvature spectrum will be suppressed by
$\Delta_{\rm iso}^{2}\propto(\rho_{\rm mis}/\rho_{\rm string})^{2}$, where
$\rho_{\rm mis}$ and $\rho_{\rm string}$ are the axion dark matter energy
densities from the ‘misalignment-like’ production and the string/domain wall
network respectively.

Secondly, if a network of strings re-enters the horizon before the QCD
crossover, it is not clear whether there is even a sub-dominant contribution
to the dark matter abundance from misalignment that inherits the long-distance
correlations eq. (27). For $\xi\simeq 1$ at the QCD crossover, such
misalignment production cannot be clearly isolated from the string/domain
network, and, relatedly, $\theta$ is not necessarily constant on large scales
during the evolution. The technical reason for this is that in the string
cores $\theta$ is not defined and so $\theta$ is not a conserved quantity on
large scales.181818As example consider a single string that is initially not
in the vacuum. The phase necessarily winds $[-\pi,\pi]$ around the string but
the interval will not be wound uniformly. However, as time progresses,
gradients will smooth out the field configuration driving the system to the
string vacuum thus changing the value of $\theta$ around the string. (Instead,
in absence of strings and with the radial mode in its vacuum everywhere the
zero mode decouples from the equations of motion and $\theta$ is a conserved
quantity on very large scales.) As a result, a string network in scaling could
erase long-distance correlations in $\theta$. However, to our knowledge
whether this indeed happens has not been properly addressed in the literature.
It could, in principle, be clarified by numerical simulations, however this is
technically difficult because only $O(10)$ e-folds of evolution can be
followed in the simulations, and moreover initial conditions that lead to a
string network that reaches scaling quickly have a power spectrum that is
already highly suppressed at large scales. We find that initial conditions
with a flat power spectrum lead to an underdense network with $\xi\ll 1$ (that
will only reach scaling beyond the reach of simulation), and in this case we
see no erasing of large distance correlations. But this could still be
consistent with the picture above, given that most Hubble patches do not yet
contain a string. Given the importance of this question for the axion post-
inflationary scenario, it would be interesting if future, larger, simulations
could provide direct evidence for such effects.

Conversely, in the late string regime, most Hubble patches do not contain a
string when the Universe has temperature $T\simeq\mathrm{GeV}$, at which time
the axion potential becomes relevant. As a result, in this case, there is a
cleanly defined misalignment contribution to the axion relic abundance, which
seems likely to retain the long-distance correlations from inflationary
fluctuations because strings are rare. We do however note that even in this
case, the small string loops that decay prior to the QCD crossover might
affect the long-distance correlations even if the density of long string is
not yet high enough for scaling.

###  5.3 The effect of overshoot on super-horizon correlations

Finally, we discuss the possibility that the relaxation of the radial mode to
the vacuum – via overshooting its potential in a significant fraction of space
– can change the large distance correlations in $\theta$. In particular, for
distances $|\mathbf{x}|$ over which the probability distribution of
$\Phi(\mathbf{x})$ immediately at the end of inflation covers many of the
“steps” in Figure 4, overshooting results in $\cos\theta(\mathbf{x})$ taking
an effectively random value from
$\\{\cos\theta_{0}(\mathbf{x}),-\cos(\theta_{0}(\mathbf{x}))\\}$, where
$\theta_{0}(\mathbf{x})$ is the value before overshooting. Consequently,
$\langle\cos\theta(\mathbf{x})\cos\theta(0)\rangle$ can average towards zero,
even if $\theta_{0}(0)$ and $\theta_{0}(\mathbf{x})$ are highly correlated.

In reality, $\cos(\theta(\mathbf{x}))$ and $\cos(\theta(0))$ are not
completely decorrelated by overshooting, and the extent of the suppression by
overshooting is important when confronting observational constraints. We have
numerically investigated realizations of $\Phi$ with inflationary initial
conditions for $H_{I}$, $\lambda$ such that $\left<\rho^{2}\right>\gg
f_{a}^{2}$. Although we leave a complete study for future work, preliminary
analysis suggests that overshooting does indeed suppress large distance
correlations, with a suppression factor that is consistent with

$\frac{\langle\cos\theta(\mathbf{x})\cos\theta(0)\rangle}{\langle\cos\theta_{0}(\mathbf{x})\cos\theta_{0}(0)\rangle}\sim(1-2F(N))^{2}~{},$
(29)

where $N=\log(|\mathbf{x}|H_{I})$ is the number of e-folds of inflation that
separates $0$ and $\mathbf{x}$ (results from simulations showing a suppression
can be found in Appendix D). Such a dependence is plausible: as expected, the
correlation vanishes if $\Phi$ is evenly distributed over both sides of the
potential after overshooting (i.e. $F\simeq 0.5$), and for a correlation to
remain neither $\Phi(0)$ or $\Phi(\mathbf{x})$ can be “randomised”.

##  6 Landscape of axion scenarios

We now turn to discuss phenomenology. Intuitively, theories such that a string
network re-enters the horizon, i.e. reaches $\xi\simeq 1$, sufficiently early
are indistinguishable from the thermal axion post-inflationary scenario (up to
possible isocurvature, as just discussed). On the other hand, if a string
network re-enters the horizon later, the evolution can be qualitatively
different with important implications for phenomenology.

###  6.1 Breaking temperature of PQ symmetry

In Section 3, we showed how the combination of inflationary fluctuations and
overshoot predicts, for a given $H_{I}$ and $\lambda$, a critical number of
e-folds $N_{s}$ that need to re-enter the horizon for $\xi\simeq 1$ to be
reached. This corresponds to $F(N_{s})$ equal to some critical value, which we
assume to be $\simeq 0.35$.191919In practice, in our plots we approximate
$N_{s}={\rm Min}[N_{s}^{\rm inf.}\,,N_{s}^{\rm o.s.}]$ from the two sources,
which is a good approximation over the majority of parameter space where one
or other formation mechanism dominates. Defining, as in eq. (28), $N_{\rm
vis}$ to be the total number of visible e-folds of inflation, the relevant
modes exited the horizon

$N_{\rm PQ}=N_{\rm vis}-N_{s}~{},$ (30)

from the beginning of visible inflation. $N_{\rm PQ}$ has a simple relation to
the temperature of the Universe when the mode re-enters the horizon. Solving
$H_{0}e^{N_{\rm PQ}}=RH$ in radiation domination one finds,

$T_{\rm PQ}=T_{0}e^{N_{\rm
PQ}+4}\left(\frac{10}{g_{*}(T_{PQ})}\right)^{1/6}~{},$ (31)

where $T_{0}$ is the CMB photon temperature today and $g_{*}(T)$ is the number
of relativistic degrees of freedom. For instance, modes that renter at
$T_{PQ}=1\,,0.15\,,0.001$ GeV exited at $N_{\rm PQ}=25\,,23\,,18$. Eqs. (31)
and (30) allow to relate $N_{s}$ to the temperature when the string network
re-enters the horizon, which is a key quantity for the phenomenology.

It is useful to compare $T_{PQ}$ with the temperature $T_{\star}$ when
$H(T_{\star})=m_{a}(T_{\star})$, around which time the axion potential becomes
cosmologically relevant and the axion starts oscillating. This reads

$T_{\star}=\left[\frac{3\sqrt{10}}{\pi
g_{*}^{1/2}(T_{\star})}\left(\frac{\Lambda_{1}}{\Lambda}\right)^{\frac{\alpha_{m}}{2}}\frac{M_{p}}{f_{a}}\right]^{\frac{2}{\alpha_{m}+4}}\simeq
5.5\,{\rm GeV}\left[\frac{10^{8}\,{\rm
GeV}}{f_{a}}\right]^{\frac{1}{6}}\left[\frac{60}{g_{*}(T_{\star})}\right]^{\frac{1}{12}}\,,$
(32)

where $m_{a}(T)=m_{a}(T=0)(\Lambda_{1}/T)^{\alpha_{m}/2}$, with
$\Lambda_{1}\simeq 150$ MeV, $\Lambda=75$ MeV, $\alpha_{m}\simeq 8$, and
$m_{a}(0)=\Lambda^{2}/f_{a}$.

Note that $N_{\rm vis}$ depends on $H_{I}$ and the, uncertain, details of the
cosmological expansion, e.g. soon after inflation, see eq. (28). For the
scenarios of most interest to us, $H_{I}\sim f_{a}\sim 10^{10}~{}{\rm GeV}$.
We also assume that after inflation the cosmological history is simply
(relatively slow) reheating followed by radiation domination, in which case
$N_{\rm vis}$ depends logarithmically on the reheating temperature $T_{\rm
RH}$, c.f. eq. (28). As a reference value, we take $N_{\rm vis}=50$.

###  6.2 Regimes

Depending on the value of $N_{s}$ different scenarios are realized:

Figure 8:  Estimate of number $N_{s}$ of e-folds of inflation that must re-
enter the horizon for the resulting string network to reach a density of
roughly one string per Hubble patch, i.e. scaling. We account for production
of strings both directly from inflationary fluctuations and from the overshoot
mechanism, with the criteria fraction $F(N_{s})=0.35$ as defined in eq. (18)
(see Appendix B for the analogous plot for $0.4$). For $H/f_{a}\lesssim
2\lambda^{1/4}$ (below the dotted red line) production of strings directly by
inflationary fluctuations dominates, while in the opposite regime the
overshoot mechanism is mostly responsible for the generation of the network.
Above the red line that separates the two regimes the radial mode $\rho$ gets
an expectation value $\gtrsim 2f_{a}$ during inflation and below its
distribution is peaked in the vacuum. In the dark grey region the radial mode
is heavier than Hubble so that the stochastic formalism for $\Phi$ is not
applicable.

####  $N_{s}>N_{\rm vis}$: Broken PQ symmetry

If $N_{s}>N_{\rm vis}$ then, for any practical purpose, the PQ symmetry is
broken throughout the evolution of the universe realizing a pre-inflationary
axion scenario. Because we assume $H_{I}\gtrsim f_{a}$, this regime is
excluded by isocurvature perturbations [90, 91, 92, 93, 94, 88, 89] (unless
any of the following happens: (1) $\lambda\lesssim 10^{-19}$ from eq. (24),
(2) the initial condition for $\Phi$ at the beginning of the visible inflation
is very large, much bigger than the expected value from FP equilibrium
distribution).

Note that if one were to consider only string production directly by
inflationary fluctuations, we would obtain $N_{s}>N_{\rm vis}$ for any
$\lambda\lesssim 0.05$ since the PQ field is approximately constant in our
Hubble patch (at least for our assumed $\Phi_{0}$), see Figure 6 (left) and
[74]. However, we have shown that due to the overshoot mechanism a network may
still form when $\lambda$ is small, Figure 6 (right), so the region of small
$\lambda$ can still be viable, falling into one of the two subsequent classes
rather than the present one.

####  $N_{s}<N_{\rm vis}-25$: Restored PQ symmetry

For such $N_{s}$, at the end of inflation $\Phi$ has fluctuations over our
observable Universe that are large enough to form strings (either directly or
by overshooting). Moreover, the resulting string network re-enters the horizon
before the QCD crossover (because modes that reenter the horizon at around
this time left the horizon $\simeq 25$ e-folds after the start of visible
inflation). This regime leads to predictions that are similar to the axion
post-inflationary scenario.

In this scenario the axion relic abundance receives at least a contribution
from the axions emitted during the string scaling regime. This is uncertain
because the result depends on the energy density spectrum of the axions
emitted by the strings, and simulation results need to be extrapolated over
orders of magnitude. However, the final result likely leads to an over-
abundance of axions with respect to the misalignment, corresponding to
requiring $f_{a}\lesssim 10^{10}\,\mathrm{GeV}$ for the axion dark matter
abundance to not overclose the Universe [30] (there is disagreement in the
upper bound on $f_{a}$ arising from uncertainties in the form of the axion
emission spectrum from strings [95, 96, 97, 98, 99, 100, 31]).

####  $N_{\rm vis}-25<N_{s}<N_{\rm vis}$: Late strings

In the region of parameter space such that $N_{\rm vis}-25<N_{s}<N_{\rm vis}$,
_late strings_ , the PQ symmetry is effectively broken over the observable
Universe, but only on scales that re-enter the horizon after the QCD
crossover. This results in a string network that re-enters the horizon, i.e.
reaches $\xi\simeq 1$, only when the temperature of the Universe $\lesssim{\rm
GeV}$ (more precisely, when $T_{PQ}<T_{\star}$, see eq. (31)). At such times
the axion mass is relevant and the strings are attached to domain walls that
form at $H=m_{a}$, and (assuming $N_{W}=1$) the string-domain wall network is
immediately unstable. Such a regime can lead to very different phenomenology
compared to if scaling is reached while the axion mass is negligible. The
dynamics are effectively those of the meso-inflationary axion scenario where
the PQ is broken after the beginning of visible inflation Ref. [48].

In Figure 8 we show the regions where the different behaviours are expected as
a function $H_{I}/f_{a}$ and $\lambda$, using the criteria $F(N_{s})=0.35$ to
determine $N_{s}$, with the full $F$ in eq. (18) that incorporates both
inflationary production and overshoot (results with $F=0.4$ are shown in
Appendix B; the difference gives a feeling of the uncertainty involved in our
estimates). In green we show the region where the network forms at $T>$ GeV so
that the standard post-inflationary axion phenomenology is expected, where in
particular axions are mostly emitted by the string network. In the blue and
red regions the annihilation of the string/domain wall network takes place at
$T\lesssim 1$ GeV, i.e. after the axion mass is cosmologically relevant. In
these regions the phenomenology is vastly different from the post-inflationary
scenario. In the grey region the symmetry is always broken in the visible
Hubble patch so that this corresponds to a pre-inflationary axion (ruled out
by isocurvature). As in Figure 6, and discussed in detail in Section 3.2, the
shape of the light gray region, ‘$N_{s}>25?$’ is especially uncertain.

###  6.3 The late string relic abundance

In the late string regime, $N_{\rm vis}-25<N_{s}<N_{\rm vis}$, we can cleanly
distinguish two different contributions to the final axion abundance: _(i)_
from the misalignment and _(ii)_ from the annihilation of the string-wall
network that subsequently re-enters the horizon at $T=T_{PQ}$, related to
$N_{s}$ via eq. (31).

The former arises when the axion mass becomes cosmologically relevant at
$H\simeq m_{a}$, which happens before the inhomogeneities generated during
inflation, and the string network, re-enter the horizon. Thus at this time it
is a good approximation to consider the axion field to be homogeneous in every
causal patch, except for the few patches that contain a string (where the
axion instead wraps the fundamental domain $(-\pi,\pi]$). When $H=m_{a}$, the
axion starts oscillating around the minimum of its potential, and, at least in
those patches that do not contain a string and for which the value of the
axion $a/f_{a}\neq\pi$ (which will contain a domain wall), this leads to a
misalignment contribution to the axion dark matter abundance. Over different
Hubble patches, the initial misalignment angle is randomly distributed in the
interval $(-\pi,\pi]$, so the resulting averaged relic abundance is [101]

$\frac{\Omega_{a}^{\rm mis}}{\Omega_{\rm dm}}\simeq
0.4\,\kappa\left[\frac{f_{a}}{10^{12}\,{\rm
GeV}}\right]^{\frac{7}{6}}\left[\frac{60}{g_{*}(T_{\star})}\right]^{\frac{5}{12}}\,,$
(33)

where $\Omega_{\rm dm}$ the observed dark matter abundance and $\kappa$ is a
constant that depends on $\alpha_{m}$, with $\kappa\simeq 5$ for
$\alpha_{m}=8$. This contribution is expected to inherit isocurvature
perturbations from the inflationary power spectrum eq. (27) and, as we discuss
in Section 6.5, its abundance is constrained by CMB observations to be small,
as is the case for $f_{a}\ll 10^{12}~{}\mathrm{GeV}$.

Soon after the time when $H=m_{a}$, the axion relaxes towards the minimum of
its potential in almost all patches, oscillating like matter, except those
that contain a string or domain wall (which span several Hubble patches
connecting the strings). As mentioned, this string/domain-wall network re-
enters the horizon when the Universe’s temperature is $T_{\rm PQ}<T_{\star}$.
For $T_{\rm PQ}\ll T_{\star}$, the energy of the system at these times is
dominated by that in domain walls rather than in strings. The strings and
domain walls are expected to subsequently annihilate away in order-one Hubble
times, producing a contribution to the axion relic abundance.

The dynamics of the string-wall system are highly nonlinear, and a reliable
determination of the relic abundance that it produces would require dedicated
simulations, which are likely to be challenging given the large difference in
scales between the Hubble rate $H_{PQ}$ when strings reenter the horizon and
the other scales, $H_{PQ}\ll m_{a}\ll f_{a}$. Nevertheless, we can make a
reasonable estimate of the axion relic abundance as follows.

  * •

If $T_{\rm PQ}\gtrsim 150$ MeV, the strings and domain walls re-enter the
horizon and decay while the axion mass has not yet saturated to its zero-
temperature value. Assuming that the system instantaneously decays into non-
relativistic axions, the resulting relic abundance is

$\frac{\Omega_{a}^{\rm dw}}{\Omega_{\rm dm}}=\frac{\rho_{w}(T_{\rm
PQ})}{\rho_{r}(T_{\rm eq})}\frac{T_{\rm eq}}{T_{\rm
PQ}}\left[\frac{g_{*}(T_{\rm eq})}{g_{*}(T_{\rm
PQ})}\right]^{\frac{1}{3}}\simeq
1.2\mathcal{A}\left[\frac{f_{a}}{10^{10}\,{\rm
GeV}}\right]\left[\frac{10}{g_{*}(T_{\rm
PQ})}\right]^{\frac{1}{2}}\left[\frac{150\,{\rm MeV}}{T_{\rm
PQ}}\right]^{5}\,,$ (34)

where $\rho_{w}(T_{\rm PQ})\simeq 2\mathcal{A}\sigma(T_{\rm PQ})H_{\rm PQ}$ is
the energy in domain walls at $T=T_{\rm PQ}$ and $\sigma=\beta m_{a}(T_{\rm
PQ})f_{a}^{2}$ (with $\beta\simeq 9$) their tension, $\rho_{r}$ the energy
density in radiation and $T_{\rm eq}$ the temperature at matter-radiation
equality (MRE). $\mathcal{A}$ is an order-one parameter measuring the domain
wall area (in units of $1/H$) at the time when the string-wall system decays.
As a result, for $T_{\rm PQ}\gtrsim 150$ MeV this contribution to the dark
matter abundance is likely smaller than that from the axions that would be
produced by the string scaling regime.

  * •

Conversely, if $T_{\rm PQ}\lesssim 150~{}{\rm MeV}$ the string-wall system
decays when the mass has saturated to its zero temperature value and

$\frac{\Omega_{a}^{\rm dw}}{\Omega_{\rm dm}}\simeq
1.8\mathcal{A}\left[\frac{f_{a}}{10^{8}\,{\rm
GeV}}\right]\left[\frac{10}{g_{*}(T_{\rm
PQ})}\right]^{\frac{1}{2}}\left[\frac{1\,{\rm MeV}}{T_{\rm PQ}}\right]\,.$
(35)

Note that as the network re-enters the horizon at lower temperatures $T_{PQ}$,
the abundance $\Omega_{a}^{\rm dw}$ increases and correspondingly smaller
$f_{a}$ is needed to reproduce the full dark matter abundance. Thus, the axion
can make up the observed abundance for arbitrarily small $f_{a}$, with the
contribution in eq. (35) dominating that from the average-angle misalignment
in eq. (33). Eventually the required $f_{a}$ are in tension with the
astrophysical bounds from star cooling, $f_{a}\gtrsim 10^{8}$ GeV [102, 103,
104, 105, 106, 107, 108, 109], which corresponds to $T_{\rm PQ}\sim 1$ MeV,
around Big Bang Nucleosynthesis.

In Figure 9 we show the different behaviours and relic abundance as a function
of $f_{a}$ and $\lambda$, keeping fixed $H_{I}/f_{a}=2$. This corresponds to a
slice through Figure 8. The dark matter abundance is reproduced on the black
line. In the green region this coincides with the axion post-inflationary
scenario where the abundance, even though currently uncertain, depends only on
$f_{a}$. In the blue region where the annihilating domain walls dominate the
abundance is reproduced for smaller $f_{a}$ up to the astrophysical bound, and
follows eq. (35). In the intermediate regime where the string and domain wall
network re-enters the horizon soon after the QCD crossover, such that the
abundance predicted from domain walls, eq. (34), is smaller than that expected
from the scaling regime, we do not have a complete estimate for the dark
matter abundance, so we leave a gap in the black line.

Figure 9:  Different axion scenarios, and the corresponding dark matter
abundance, as a function of $\lambda$ and $f_{a}$ for $H_{I}=2f_{a}$ (assuming
the string network re-enters the horizon when $F=0.35$, as in Figure 8). In
the green region the string network re-enters the horizon when the temperature
of the Universe $\gtrsim 1$ GeV so that standard post-inflationary
phenomenology is expected. In blue the string-wall network re-enters the
horizon and annihilates at temperatures $\lesssim 1$ GeV, leading to novel
phenomenology. In the grey region the network never forms within our
observable horizon, as in the pre-inflationary scenario. The observed dark
matter abundance is reproduced on the black line, which requires much smaller
$f_{a}$ in the blue region compared to the standard post-inflationary
scenario. The region $f_{a}\lesssim 10^{8}$ GeV is excluded by astrophysical
constraints.

###  6.4 Axion spatial distribution and miniclusters

Similarly to the standard post-inflationary scenario, the collapse of the
string-wall system is expected to leave order-one inhomogeneities in the axion
dark matter spatial distribution at small length scales. The typical comoving
length $\lambda$ of the inhomogeneities is determined by the size of the
domain walls at the time of collapse, approximately fixed by the inverse
Hubble scale, i.e. $\lambda\simeq 1/(R_{PQ}H_{\rm PQ})$ (here $R_{\rm PQ}$ is
the scale factor when $T=T_{PQ}$, and the network re-enters the horizon). This
amounts to a contribution to the power spectrum of the axion overdensity
field, $\Delta_{\delta}^{2}$ defined in eq. (22) (specifically $\Delta_{\rm
iso}^{2}$), that is peaked at the momentum $k/R_{\rm PQ}\simeq H_{\rm PQ}$,
where it acquires an order-one value. For momenta $k/R_{\rm PQ}\ll H_{\rm
PQ}$, this contribution to $\Delta_{\rm iso}^{2}$ is expected to fall $\propto
C[k/(R_{PQ}H_{\rm PQ})]^{3}$ from causality arguments. Here $C$ is an order-
one coefficient that could be in principle extracted from numerical
simulations of the late string-wall decay. Given the large uncertainties on
the string-wall dynamics, the dark matter abundance and the spectrum, we
refrain from giving a more detail description of the axion field, and briefly
mention below some possible observational consequences of this dark matter
substructure, see [110, 49] for related discussions.

If the axion makes up a substantial fraction of the dark matter, the
overdensities collapse at around matter-radiation equality (MRE), much before
canonical structure formation takes place. A spherical overdensity with
magnitude $\delta\gtrsim 1$ undergoes nonlinear collapse when the scale factor
is $R\simeq R_{\rm eq}/\delta$, where $R_{\rm eq}$ is the scale factor at MRE.
This collapse leads to the formation of a compact gravitationally bound
object, known as axion minicluster, with average density $\rho_{b}\simeq
140\delta^{3}(1+\delta)\rho_{c}\simeq\mathcal{O}(100)\rho_{c}$, where
$\rho_{c}$ is the average dark matter density at the time of collapse [32, 33,
34, 35]. Since the collapse happens around matter radiation equality,
$\rho_{b}\simeq\,100\,{\rm eV}^{4}$, about 8 orders of magnitude larger than
average local dark matter density in the neighborhood of the Solar System.

The typical mass of the miniclusters is fixed by the dark matter mass
contained in an order-one fluctuation, and reads

${M_{b}\simeq\frac{4\pi}{3}\left(\frac{2\pi}{xk_{\rm PQ}/R_{\rm
PQ}}\right)^{3}\bar{\rho}(R_{\rm eq})\simeq
24M_{\odot}\left[\frac{10}{g_{*}(T_{\rm
PQ})}\right]^{\frac{1}{2}}\left[\frac{1\,{\rm MeV}}{T_{\rm
PQ}}\right]^{3}x^{-3}}$ (36)

where $x\simeq\mathcal{O}(1)$ parameterizes the typical physical size of the
fluctuations at $R=R_{\rm PQ}$ as $R_{\rm PQ}\lambda=2\pi/(xk_{\rm PQ}/R_{\rm
PQ})$. Owing to the larger size of the fluctuations (related to the larger
Hubble distance at the time of string-wall decay), for $T_{\rm PQ}\ll 1\,{\rm
GeV}$ the miniclusters are significantly more massive than in the usual
$N_{W}=1$ axion post-inflationary scenario, although their density is still of
the order of the Universe’s density at the time of collapse. Finally, the
$k^{3}$ tail leads to the clustering of miniclusters into larger (more massive
and less dense) compact halos at larger scales [111, 112, 113, 114, 115, 116].

###  6.5 Isocurvature constraints

As mentioned, the axion relic abundance from the late decay of domain walls is
expected to have a white noise isocurvature spectrum $\Delta_{\rm
iso}^{2}(k)=C[k/(R_{PQ}H_{\rm PQ})]^{3}$ at momenta $k/R_{PQ}\ll H_{\rm PQ}$,
but are not expected to inherit the correlations from inflationary
fluctuations, eqs. (24) or (27). Such white-noise isocurvature perturbations
are constrained by CMB [87] and Lyman-$\alpha$ observations [117, 118, 119,
120, 121, 122, 123], which bound the fraction of isocurvature fluctuations
$r_{\rm iso}\equiv\Delta_{\rm iso}^{2}(k_{\rm
CMB})/\Delta_{\mathcal{R}}^{2}(k_{\rm CMB})$. In the late string-wall regime
we have

${r^{1/2}_{\rm iso}=\frac{\Omega_{a}^{\rm dw}}{\Omega_{\rm
dm}}\sqrt{\frac{C}{A_{s}}\frac{k_{\rm CMB}^{3}}{H_{\rm PQ}^{3}}}\simeq 4\cdot
10^{-4}~{}C^{\frac{1}{2}}\mathcal{A}\left[\frac{f_{a}}{10^{8}\,{\rm
GeV}}\right]\left[\frac{10}{g_{*}(T_{\rm
PQ})}\right]^{\frac{3}{4}}\left[\frac{1\,{\rm MeV}}{T_{\rm
PQ}}\right]^{\frac{5}{2}}}\,.$ (37)

Lyman-$\alpha$ observations require $r^{1/2}_{\rm iso}\lesssim 4\cdot 10^{-3}$
for a white-noise power spectrum [119, 120], and provide a bound on the latest
decay time in this scenario, which is slightly less stringent than the
combination of those from dark matter overproduction and the astrophysical
bound on $f_{a}$, see [48].

In the late string regime, in addition to the axion relic abundance from
domain walls there is a cleanly separated, sub-dominant, contribution from
misalignment. As discussed in Section 6.3, this is produced around the time of
the QCD crossover and inherits any super-horizon correlations in $\cos\theta$
that remain at this time, and given that the axion is already massive when the
network annihilates it is difficult to imagine that the correlations are
destroyed at late times. Note however that the super-horizon correlations at
the QCD crossover might be suppressed relative to those immediately at the end
of inflation by overshoot, as discussed in Section 5.3. Given that the string
network is expected to re-enter the horizon when $F\simeq 0.3\div 0.45$ is
reached, which must happen long before CMB time, the corresponding $F$ at CMB
scales is likely to be somewhat closer to $0.5$. Therefore, if eq. (29) turns
out to be the true dependence, the resulting suppression can be significant.
Moreover, it is possible that the small loops that are temporarily present
even in an underdense network until they decay might partly erase the super-
horizon correlations prior to the QCD crossover. Considering the most
pessimistic scenario that neither overshoot or the dynamics of loops erases
the initial correlations, the resulting isocurvature power spectrum would be
given by eq. (27), suppressed by the fraction squared of dark matter that
arises from misalignment,

$\Delta_{\rm iso}^{2}(k_{\rm CMB})\sim\frac{\Omega_{\rm mis}^{2}}{\Omega_{\rm
wall}^{2}}e^{-2N_{\rm vis}/N_{11}}~{}.$ (38)

Considering inflationary production, we expect $N_{s}\simeq N_{11}$ up to a
numerical factor that we do not know reliably (e.g. assuming a critical
$F=0.35$, $N_{s}\simeq 1.5N_{11}$, c.f. eq. (21) where
$a_{11}b_{11}\simeq-0.65$). Consequently, in the late string regime with
$25\lesssim N_{s}\lesssim 50$, we expect the suppression from the exponential
factor in eq. (38) to be at least two orders of magnitude. (For overshoot
production $N_{\rm vis}\ll N_{11}$ so the exponential does not lead to a
suppression but, as mentioned in Section 5.3, an extra suppression from the
overshooting dynamics is expected.) Meanwhile, from eq. (33) the abundance
from misalignment is suppressed for small $f_{a}$, with $\Omega_{a}^{\rm
mis}/\Omega_{a}^{\rm dw}\simeq 10^{-5}$ for $f_{a}\simeq
10^{8}~{}\mathrm{GeV}$. Nevertheless, eq. (38) remains a potentially strong
constraint given that $\Delta_{\rm iso}^{2}(k_{\rm CMB})\lesssim 10^{-10}$
from Planck data [87] (so even a fraction $10^{-4}$ of dark matter with with
order one isocurvature perturbations would be excluded). At face value, this
would require $f_{a}$ close to the current astro-physical limit for the late
string regime bound to survive cosmological bounds, although we stress that
our discussion is only tentative and future dedicated study would be useful.

##  7 Discussion

In this work we explored a new regime of the minimal QCD axion that leads to
novel phenomenology, including a modified dark matter abundance, dominantly
produced by domain walls annihilating after the QCD crossover, and larger mass
axion miniclusters than in the standard post-inflationary scenario. Our
discussion applies in particular to the simplest concrete realization of the
post-inflationary axion, the KSVZ model with domain wall number $N_{W}=1$, and
occurs over order-one ranges of $H_{I}/f_{a}$ and $\lambda$.

We focused on the scenario where the PQ symmetry is not restored by thermal
effects, $T_{\rm max}\lesssim f_{a}$. According to the common lore, for
$H_{I}/(2\pi)\gtrsim f_{a}$ the PQ symmetry is restored during inflation while
for $H_{I}/(2\pi)\lesssim f_{a}$ it is broken throughout the cosmological
history of the universe.202020We note that additional couplings between $\Phi$
and the inflaton could change these expectations [124]. Very different
dynamics and phenomenology ought to apply to the two scenarios: in the first
case a string network is expected to form realizing the axion post-
inflationary scenario, while in the latter the axion pre-inflationary scenario
would be realized.

Clearly this picture must be incomplete and at least in the intermediate
region $H_{I}/(2\pi)\simeq f_{a}$ deviations from post or pre-inflationary
axion cosmology should be expected. Moreover, the boundaries of the two
standard regimes are potentially important. As pointed out in Ref. [45], the
dynamics in fact depend on both $H_{I}/f_{a}$ and the quartic coupling
$\lambda$ of the PQ complex field. Roughly the landscape is the following:

  1. 1.

If the quartic is sizable, say $\lambda\gtrsim 0.1$, and $H_{I}\gtrsim 2f_{a}$
inflationary fluctuations randomize the field efficiently so that a string
network forms soon after the end on inflation. Subsequently, the network
rapidly approaches the scaling solution so that the details of the initial
conditions are erased and the usual post-inflationary axion phenomenology is
expected.

  2. 2.

For smaller $\lambda\lesssim 0.1$ the network formed soon after inflation is
underdense (i.e. less than a string per Hubble patch) and will re-enter the
horizon at later and later times during radiation domination until it re-
enters after the axion starts to oscillate, i.e. at $T\lesssim 1$ GeV. Once
this happens the phenomenology becomes very different, effectively realizing
the meso-inflationary scenario of Ref. [48, 49] (but in a minimal theory
without any additional interactions between e.g. $\Phi$ and the inflaton). The
string network cannot reach the scaling regime and it effectively annihilates
as soon as the perturbations re-enter the horizon. In this case, dark matter
axions are dominantly produced by domain walls, allowing the full observed
dark matter abundance to be obtained for smaller $f_{a}$ than in the axion
post-inflationary scenario, down to the astrophysical bound $f_{a}\gtrsim
10^{8}$ GeV (additionally, possible gravitational wave signals and primordial
black holes resulting from similar dynamics have been considered in Refs.
[125, 126]). These features are reminiscent of theories that reach scaling
early, with $N_{W}>1$ and a small energy bias between the vacua, in which case
a domain wall network also survives beyond the QCD crossover [127]. However,
$N_{W}>1$ theories require an energy bias in a narrow range and might not even
be viable without reintroducing a strong CP problem [128].

  3. 3.

For even smaller $\lambda$ the PQ field is roughly homogeneous immediately
after the end of inflation so that no topologically non-trivial configurations
are initially produced. Contrary to different statements in the literature, we
find that a string network can still form during reheating, for some values of
$H_{I}/f_{a}$ and $\lambda$. During inflation indeed the PQ field can be
displaced from the vacuum so that during reheating it can undergo oscillations
(or ‘overshoot’) over the top of its potential. In this case even small
fluctuations allow the field to be randomized, and if this happens efficiently
enough a string network forms, which can re-enter the horizon either at
$T\gtrsim 1$ GeV or at $T\lesssim 1$ GeV (leading to the phenomenology in 1.
and 2., respectively). If a network is not produced instead such theories are
ruled out by isocurvature constraints.

We believe that the qualitative behaviour outlined above is robust. However,
obtaining precise quantitative predictions in terms of $H/f_{a}$ and $\lambda$
is challenging. We have approached this analytically by studying the evolution
of $\Phi$ during and after inflation, which gives a rough estimate of the
temperature when the initially underdense network is expected to reach one
string per Hubble volume, and thus the scaling solution. Numerical simulations
are needed to establish precisely the landscape in the plane
$(H_{I}/f_{a}\,,\lambda)$, but these suffer from the challenges of large scale
separations. In this paper we have also made a first step in this direction,
and the results of our simulations suggest that the string density in an
underdense network is indeed related to the spread of the radial mode over the
top of its potential, determined by the quantity $F$ in eq. (18), and that
overshoot can lead to a network. More work would be useful to better test the
relation between $F$ and the string network density, e.g. starting from a
wider range of initial conditions including initially more underdense
networks. It would also be interesting to study the evolution and eventual
destruction of a network that has not yet re-entered the horizon when the
axion mass becomes cosmologically relevant. From preliminary investigation,
this appears challenging in the simulations with $\simeq 2000^{3}$ lattice
points that we have the computational power for, due to the multiple, widely
separated, scales in the problem ($H_{PQ}\ll m_{a}\ll f_{a}$, where $H_{PQ}$
is the Hubble rate when strings reenter the horizon), but progress might be
possible with very large grids.

From a phenomenological point of view the main novelties arise when the
network annihilates at temperatures $T\lesssim 1~{}\mathrm{GeV}$. As discussed
in Section 5, in this regime misalignment gives a contribution to the axion
relic abundance (sub-leading to that from domain walls discussed above in 2.),
which is likely to induce isocurvature perturbations that are severely
constrained by the CMB. A future precise determination of the resulting
isocurvature power spectrum is therefore crucial to determine the viability of
this regime. This also raises the question of how efficiently super-horizon
correlations (as produced by inflation) are erased in the axion post-
inflationary scenario when scaling is reached early, assuming the PQ symmetry
is not restored by thermal effects. Although it is normally thought that such
a regime is automatically safe from such constraints, as is the case for the
thermal post-inflationary scenario, this relies crucially on the dynamics of
strings in a way that, although expected, is yet to be confirmed in numerical
simulations. Even a tiny contribution to isocurvature could exclude such
models so it is important to determine how efficiently the large-scale
correlations are erased by the string network. We hope to address this
question with numerical simulation in the future.

The viability of this late string regime is of crucial importance because it
opens the possibility of having QCD axion dark matter at low decay constants,
down to those allowed by current astrophysical bounds $f_{a}~{}\simeq 10^{8}$
GeV, even in the most minimal KSVZ setup. This is exciting for direct
detection experiments planning to search for QCD axions with large mass
$m_{a}\lesssim 10^{-2}$ eV, such as IAXO [129, 130] among others [131, 132,
133]. Additionally, for cosmological observations, it opens the possibility of
the simultaneous existence of cold axions that form the bulk of the observed
dark matter and a component of thermal axions (see e.g. [134, 135, 136, 137]
and references therein) detectable as an ‘effective additional neutrino
species’, $\Delta N_{\rm eff}$, via forthcoming CMB experiments such as CMB-S4
[138] and Simons Observatory [139] and in combination with Large Scale
Structure surveys (e.g. Euclid [140]).

A final important caveat is that in the regime $H_{I}\gg f_{a}\lambda^{1/4}$
our results depend on the assumption that, prior to the visible e-folds of
inflation, $\Phi$ acquired a value typical of its equilibrium distribution
with the value of the radial mode of order eq. (5). If the radial mode instead
had a much smaller value at the start of visible inflation, strings would form
more easily than we predict because the radial mode would immediately spread
over the top of its potential. Conversely, if the radial mode had much larger
values at such times, say of order the Planck scale, the dynamics would again
be different to our results. In the absence of a complete theory of inflation
and the cosmology prior to this, the initial value of $\Phi$ will remain an
important theoretical uncertainty in this part of parameter space.

#### Acknowledgements

We wish to thank Andrea Tesi and Giovanni Villadoro for discussions. MR was
supported by MIUR grant PRIN 2017FMJFMW. EH 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. The work of MG is supported
by the Alexander von Humboldt foundation and has been partially funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 491245950.
The work of AN is supported by the grants PID2019-108122GB-C32 from the
Spanish Ministry of Science and Innovation, Unit of Excellence María de Maeztu
2020-2023 of ICCUB (CEX2019-000918-M) and AGAUR2017-SGR-754. AN is grateful to
the Physics Department of the University of Florence for the hospitality
during the course of this work.

##  Appendix A Approximate correlation functions

In this Appendix we present approximate formulas for the power spectrum of the
PQ field at the end of inflation, see also the appendix of [141]. The
derivations below give an alternative intuitive understanding for the
correlation functions.

The power spectrum of a free scalar field $\phi$ with effective mass
$M\lesssim H_{I}$ during inflation reads

$P(k)=\frac{2\pi^{2}}{k^{3}}\Delta^{2}(k)=\frac{H_{I}^{2}}{2k^{3}}\left(\frac{k}{RH_{I}}\right)^{3-2\nu}\,,~{}~{}~{}~{}~{}~{}~{}~{}\nu=\sqrt{\frac{9}{4}-\frac{M^{2}}{H_{I}^{2}}}\approx\frac{3}{2}-\frac{M^{2}}{3H_{I}^{2}}\,.$
(39)

From this we can compute the field’s variance at the end of inflation,

$\langle\phi^{2}(\mathbf{x})\rangle=\frac{H_{I}^{2}}{(2\pi)^{2}}\int_{H_{I}R_{e}e^{-N_{I}}}^{H_{I}R_{e}}\frac{dk}{k}\left(\frac{k}{R_{e}H_{I}}\right)^{\frac{2M^{2}}{3H^{2}}}\approx\left\\{\begin{array}[]{lcl}\displaystyle\frac{3H_{I}^{4}}{8\pi^{2}M^{2}}&&N_{I}>\frac{H_{I}^{2}}{M^{2}}\,,\\\
\displaystyle
N_{I}\frac{H_{I}^{2}}{(2\pi)^{2}}&&N_{I}<\frac{H_{I}^{2}}{M^{2}}\\\
\end{array}\right.\,,$ (40)

where $R_{e}$ is the scale factor at the end of inflation and $N_{I}$ is the
total number of e-folds of inflation. Eq. (40) shows that for small $M\ll
H_{I}/\sqrt{N_{I}}$ it is a good approximation to neglect the potential and
the variance undergoes a random walk, while for $M\simeq H_{I}/\sqrt{N_{I}}$
the variance reaches a constant value within the observable number $N_{I}$ of
e-folds. The fact that $\langle\phi^{2}\rangle$ becomes constant is due to the
balance between the stochastic inflationary fluctuations and the classical
force relaxing the field to the minimum. The convergence of the integral in
the IR is due to the suppression of long modes, with momentum smaller than
$H_{I}R_{e}e^{-N_{I}}$, that thus cease to contribute to
$\langle\phi^{2}\rangle$.

Let us first consider the quartic potential $V=\lambda\phi^{4}/4$. The
effective mass is $M^{2}=3\lambda\langle\phi^{2}\rangle$ and, from eq. (40),
in the stationary regime we thus have

$\langle\phi^{2}\rangle\approx\frac{H_{I}^{2}}{\sqrt{8\pi^{2}\lambda}}\,.$
(41)

We can extend eq. (40) to several scalar fields, obtaining that the
expectation value of $\phi_{i}\phi_{j}$ is

$\langle\phi_{i}\phi_{j}\rangle=(M^{2})^{-1}_{ij}\frac{3H_{I}^{4}}{8\pi^{2}}\,.$
(42)

For the PQ field $V=\lambda/4(\phi_{1}^{2}+\phi_{2}^{2}-f_{a}^{2})^{2}$ so
that,

$M^{2}=\lambda\begin{pmatrix}3\langle\phi_{1}^{2}\rangle+\langle\phi_{2}^{2}\rangle&2\langle\phi_{1}\phi_{2}\rangle-
f_{a}^{2}\\\
2\langle\phi_{1}\phi_{2}\rangle&\langle\phi_{1}^{2}\rangle+3\langle\phi_{2}^{2}\rangle-
f_{a}^{2}\\\ \end{pmatrix}\,.$ (43)

In the limit $H_{I}\gg\lambda^{1/4}f_{a}$, from eq. (42)

$\lambda(3\langle\phi_{1}^{2}\rangle+\langle\phi_{2}^{2}\rangle)\langle\phi_{1}^{2}\rangle=\frac{3H_{I}^{4}}{8\pi^{2}}\,,~{}~{}~{}~{}~{}\lambda(3\langle\phi_{2}^{2}\rangle+\langle\phi_{1}^{2}\rangle)\langle\phi_{2}^{2}\rangle=\frac{3H_{I}^{4}}{8\pi^{2}}\,,~{}~{}~{}~{}~{}\langle\phi_{1}\phi_{2}\rangle=0\,,$
(44)

with solution

$\langle\phi_{1}^{2}\rangle=\langle\phi_{2}^{2}\rangle=\frac{\sqrt{3}H_{I}^{2}}{4\pi\sqrt{2\lambda}}$
(45)

By plugging eq. (45) in eq. (43), and these in eq. (39), we can also compute
the power spectra,

$\Delta^{2}_{\phi_{1}}(k)=\Delta^{2}_{\phi_{2}}(k)\approx\frac{H_{I}^{2}}{(2\pi)^{2}}\left(\frac{k}{R_{e}H_{I}}\right)^{\gamma}\,,~{}~{}~{}~{}~{}\gamma=\sqrt{\frac{2\lambda}{3\pi^{2}}}\,.$
(46)

The results above are equivalent to using the Hartree-Fock approximation, see
[42]. The time scale, or the number of e-folds, for the fields to reach the
equilibrium distribution can be estimated by studying at the evolution of the
expectation value $\langle\phi_{i}\phi_{j}\rangle$ during inflation. This is
determined by the equation,

$\frac{\partial}{\partial
t}\langle\phi_{i}\phi_{j}\rangle=\frac{H_{I}^{3}}{4\pi^{2}}\delta_{ij}-\frac{2M^{2}_{ik}}{3H_{I}}\langle\phi_{k}\phi_{j}\rangle\,,$
(47)

where the first term corresponds to the stochastic noise due to inflationary
fluctuations while the second is the classical drift. Eq. (47) is an
approximation to the FP equation that assumes that the fields are gaussianly
distributed, see [42] for the derivation. The time independent solution is
given by eq. (45). Indeed, taking $\langle\phi_{i}\phi_{j}\rangle=0$ and
$\langle\phi_{1}^{2}\rangle=\langle\phi_{2}^{2}\rangle$, eq. (47) simplifies
to,

$\frac{\partial}{\partial
t}\langle\phi_{1}^{2}\rangle=\frac{H_{I}^{3}}{4\pi^{2}}-\lambda\frac{2\langle\phi_{1}^{2}\rangle(4\langle\phi_{1}^{2}\rangle-
f_{a}^{2})}{3H_{I}}$ (48)

and at equilibrium one recovers eq. (45). The stationary solution is reached
equating the random walk $\frac{H_{I}^{2}}{4\pi^{2}}N$ to the equilibrium
value $\langle\phi_{1}^{2}\rangle$ in eq. (45), which gives a number of
e-folds $N\sim\sqrt{\frac{3\pi^{2}}{2\lambda}}$.

###  A.1 Isocurvature estimates

We can also compute correlation functions relevant for isocurvature, as
discussed in Section 5. Assuming that the system has reached the asymptotic
distribution we have,

$\langle\cos\theta(x)\cos\theta(0)\rangle=\langle\frac{\phi_{1}}{\sqrt{\phi_{1}^{2}+\phi_{2}^{2}}}(x)\frac{\phi_{1}}{\sqrt{\phi_{1}^{2}+\phi_{2}^{2}}}(0)\rangle\approx\frac{1}{\langle\phi_{1}^{2}+\phi_{2}^{2}\rangle}\langle\phi_{1}(x)\phi_{1}(0)\rangle\,,$
(49)

In Fourier space,

$\Delta^{2}_{\cos\theta}(k)=\frac{k^{3}}{2\pi^{2}}\int
d^{3}xe^{ikx}\langle\cos\theta(x)\cos\theta(0)\rangle=\sqrt{\frac{\lambda}{6\pi^{2}}}\left(\frac{k}{R_{e}H_{I}}\right)^{\sqrt{\frac{2\lambda}{3\pi^{2}}}}$
(50)

where we used the expression for $\langle\phi^{2}_{1,2}\rangle$ in equilibrium
in eq. (45) and for the power spectrum for
$\langle\phi_{1}(x)\phi_{1}(0)\rangle$. The suppression at the CMB scales is
thus,

$\Delta^{2}_{\cos\theta}(k_{\rm
CMB})\approx\sqrt{\frac{\lambda}{6\pi^{2}}}e^{-2N_{\rm
vis}\sqrt{\frac{\lambda}{6\pi^{2}}}}$ (51)

This has the same form of eq. (27) with
$N_{11}=\Gamma_{11}/H_{I}=\sqrt{6\pi^{2}/\lambda}\approx 8/\sqrt{\lambda}$.

##  Appendix B More on string production

In this Appendix we present a few additional results on our predictions of the
number of e-folds of inflation that have to re-enter the horizon for a string
network with $\xi\simeq 1$ to be reached (at which point the string scaling
regime starts if $T_{PQ}\gtrsim 1$ GeV, or the string-wall network annihilates
if $T_{PQ}\lesssim 1$ GeV). In particular, in Figure 10 we show the same
results as in Figure 6 in the main text, but with critical fraction $F=0.4$
rather than $0.35$, to give an indication of the uncertainty on our results
(as shown in Figure 7, our results from numerical simulations are not precise
enough to predict the critical $F$ to reach $\xi=1$ to better than the range
$0.3\div 0.45$; the critical value of $F$ might also vary depending on whether
the string network is dominantly produced by inflationary fluctuations or the
overshoot mechanism). In the left panel we show $N_{s}$ considering only
inflationary production of strings. Given that in this case $N_{s}$ has only a
logarithmic dependence on $F$, c.f. eq. (21), the shift in the contours of
$N_{s}=25,~{}50$ are relatively mild, with changes of roughly less than order-
one to the values of $\lambda$ and $H_{I}/f_{a}$ that, respectively, set the
approximately vertical and horizontal contour regions. Meanwhile, in the right
panel we show contours of $N_{s}$ considering only overshoot. The changes
relative to Figure 6 are slightly larger in this case, with the structures on
the left of the plot (corresponding to the radial mode starting in the middle
of one of the ‘overshoot steps’) getting somewhat larger. However, we stress
that the exact details of the parameter space are especially uncertain in this
regime given the unknown value of $\Phi_{0}$, so the additional uncertainty
from $F$ is not a major issue.

Figure 10: Number of e-folds of inflation $N_{s}$ that must re-enter the
horizon for string network close to scaling to be reached, analogous to Figure
6 in the main text, except requiring a critical fraction $F=0.4$ rather than
$0.35$. On the left, considering string production only from inflationary
fluctuations and on the right only via the overshoot mechanism.

##  Appendix C Further details of the numerical simulations

###  C.1 Description of simulations

We generate initial conditions for our simulations (made of $n^{3}$ lattice
points) by the following procedure, equivalent to evolving the Langevin
equation that the field solves during inflation, [45]:

  * •

We start with a homogeneous field $\Phi=\Phi_{0}$ on all lattice points. At
this point, the entire simulation volume represents a size $H_{I}^{-3}$ region
$N=\log_{2}(n)$ e-folds before the end of inflation (with the value of the
field coarse grained over $H_{I}^{-1}$).

  * •

To reproduce one doubling time of inflation (during which the initial
$H_{I}^{-3}$ volume grows to $8H_{I}^{-3}$) we split the homogeneous region
into 8 patches. Inflationary fluctuations are accounted for by giving
$\phi_{1}$ and $\phi_{2}$ (the canonically normalised real and imaginary parts
of $\Phi$) in each of the new regions a ‘kick’ drawn from a Gaussian
distribution with variance $\log(2)H_{I}^{2}/(4\pi^{2})$ and mean $0$.
Finally, we evolve $\Phi$ in each of the regions according to the classical
equations of motion (for a homogeneous field) for a time $\log(2)H_{I}^{-1}$.

  * •

This procedure is repeated, with each homogeneous patch successively split
into $8$, kicked and evolved until each grid point corresponds to an
individual region of size $H_{I}^{3}$.

This generates an approximation to a realisation of the field configuration
produced by inflation. As well as giving a distribution with a probability
distribution of $\Phi$ that approximates the solution to the FP equation, such
an approach captures higher order correlation functions, e.g. non-
Gaussianities, that are present when the potential is non-negligible.

Subsequently the system is evolved forward in time from $t=H_{I}^{-1}$ using
the classical equations of motion from the Lagrangian in eq. (1) discretised
with a leapfrog algorithm that is accurate to second order in the timestep. We
parameterize the cosmic time in simulations by $\log(m_{\rho}/H)$. As
mentioned in the main text, the dynamics of the string network are only
captured without significant systematics provided there is at least $\sim 1$
lattice point per string core (which is of size $m_{\rho}^{-1}$) and
$HL\gtrsim 2$ (where $L$ is the physical size of the box). With the Lagrangian
of eq. (1), the string core size decreases in comoving coordinates $\propto
R^{-1}$ (where $R$ is the scale factor). Consequently, to maximise the range
in $m_{\rho}/H$ of such simulations, the string cores must initially contain
$\sqrt{n}$ lattice points. In other words, $m_{\rho}\Delta=1/\sqrt{n}$
initially, where $n$ the number of lattice points. This is unfortunate for our
purposes, because half of the e-folds worth of inflationary fluctuations are
already inside the string cores at the start of the simulation, and we could
only study $\frac{1}{2}\log n$ e-folds of inflation re-entering the horizon.
Instead, we simulate the so-called fat string system, i.e. with $\lambda$
changing with time in such a way that $m_{\rho}\propto R^{-1}$, and we set the
lattice spacing $\Delta$ as $m_{\rho}\Delta=1$; we start the simulations at
$H=m_{\rho}$. In this way we can study the effects of almost the full
generated $\log n$ e-folds worth of inflationary fluctuations reentering the
horizon (although, as mentioned in the main text, the equations of motion
differ from the physical system).

In practice, for initial conditions that give small $\xi$, statistical
uncertainties on $\xi$ become large prior to the maximum time allowed by the
systematic uncertainties discussed above (even with $\simeq 10$ repeated
runs). Our data sets include sufficiently many simulations that the relative
statistical uncertainty on $\xi$ is within $10\%$ (this is calculated via the
standard error on the mean over repeated runs). We also only present results
from timeshots $\log(m_{\rho}/H)\gtrsim 3$; prior to this the strings have a
thickness comparable to the horizon size and are not cleanly defined objects
(in the case of data sets with a large radial mode expectation value during
inflation, we additionally only consider timeshots once the radial mode has
settled down to the minimum of its potential in most of the simulation
volume).

To produce the $x$ axis of Figure 7, we calculate the value of $F$, defined by
eq. (18), as a function of the simulation time. As described in the main text,
the relation between simulation time and $F$ depends on the values of
$H_{I}/f_{a}$, $\lambda$ (used to generate the initial conditions), as well as
the initial radial mode value $\rho_{0}$. One additional complication arises
in the case of large initial field displacement (Figure 7 right). In this
regime, in fact, the value of the radial mode $\rho$ $N$ doubling times before
the end of inflation is not unique in what corresponds to the comoving
simulation volume. Instead, it has a distribution that occurs because the
radial mode has evolved from the delta-function initial condition (peaked at
$\rho_{0}$) to a non-trivial distribution during the preceding $11-N$ doubling
times of inflation ($11$ is the total number of doubling times in the
simulation, and appears given that our simulations with $2048^{3}$ lattice
points lead to a total of $\log_{2}(2048)=11$ doublings). Thus, to compare to
$\xi$ evaluated over the full simulation box (when $N$ doubling times of
inflation have re-entered the horizon), we need to calculate $F$ taking into
account this distribution. Although the results obtained for $F$ as a function
of simulation time are not numerically very different from those obtained
assuming a fixed value the radial mode $N$ doubling times before the end of
inflation, this procedure results in $F$ being an increasing function of
simulation time, as is expected physically. To produce the plot in Figure 6 we
neglected this issue, and instead assumed that the radial mode value at $N$
e-folds before the end of inflation is unique.

In more detail, we calculate an average $F(N)$:

$\bar{F}(N)=\int_{0}^{\infty}d\rho_{i}~{}P\left(\rho_{i}\right)F\left(N\right)|_{\rho=\rho_{i}}~{},$
(52)

where $P\left(\rho\right)=\int_{0}^{2\pi}d\theta~{}\rho
P(\phi_{1},\phi_{2})|_{\phi_{1}^{2}+\phi_{2}^{2}=\rho^{2}}$ is the probability
distribution of $\rho$ evaluated $N$ e-folds before the end of inflation (due
to the preceding e-folds of inflation, in our case corresponding to $11-N$
doublings) and $F\left(N\right)|_{\rho=\rho_{i}}$ is the function defined by
eq. (18) setting $\rho=\rho_{i}$ $N$ e-folds before the end of inflation. For
simulations with $\lambda=0$ in the initial conditions, $P$ is a Gaussian and
we evaluate eq. (52) directly. These initial conditions lead to a field
displaced by $\rho_{0}\gg f_{a}$ after inflation and should produce strings
via overshoot, thus the issue discussed above is important to account for. For
$\lambda\neq 0$ (and the values of $H_{I}/f_{a}$ used) instead the radial mode
is not significantly displaced and we simply approximate that the probability
distribution of $\rho$ is strongly peaked in the vacuum, which amounts to not
carrying out the average in eq. (52), and is a good approximation for all
$(H_{I}/f_{a},\lambda)$ plotted (additionally, it is numerically impractical
to evaluate the required integrals for each timeshot for each set of initial
conditions; we have however checked for individual data points that this
approximation does not significantly affect the inferred $\bar{F}$).

###  C.2 The string network

Figure 11:  The same simulation results as in Figure 7, but with the string
density $\xi$ plotted as a function of simulation time, parameterized by
$\log(m_{\rho}/H)$. In the left panel, results obtained from initial
conditions such that strings form as a result of fluctuations directly over
the top of the potential. Right panel, results with initial conditions such
that strings form as a result of the radial mode overshooting the potential.

To clarify the properties of underdense string networks, we plot in Figure 11
the string density $\xi$ as a function of simulation time (parameterised by
$\log(m_{\rho}/H)$) rather than $F$. At a given simulation time the string
density varies by a factor of $\sim 2$ among initial conditions with different
$H_{I}/f_{a}$ and $\lambda$, and the rate of increase in $\xi$ also varies. It
is therefore non-trivial that that the values of $\xi$ between different
initial conditions coincide once $\xi$ is expressed as a function of $F$, as
they do in Figure 7. It would however be interesting if simulations carried
out on large grids could allow for underdense networks that have a greater
range of $\xi$ at a given $\log(m_{\rho}/H)$ to test how well our proposed
dependence on $F$ holds. As mentioned in the main text, the increase in $\xi$
with $\log(m_{\rho}/H)$ is exponential rather than linear. This is
characteristic of the network being underdense compared to the critical $\xi$
such that string annihilation balances the rate at which strings reenter the
horizon, as opposed to the logarithmic growth of $\xi$ on the attractor as a
result of this critical point slowly moving with $\log(m_{\rho}/H)$.

##  Appendix D The effect of overshoot on the density power spectrum

In this Appendix, we present a preliminary investigation of the effects of the
relaxation of the PQ field to the minimum via overshoot on correlation
functions at super-horizon scales, as discussed in Section 5.3. First, in
Figure 12 we plot results for the string density $\xi$ and the power spectrum
of $1-\cos(\theta)$, $\Delta^{2}_{1-\cos\theta}$, at different cosmological
times (parameterized by $\log(m_{\rho}/H)$) for numerical simulations carried
out similarly to described in Appendix C (results are shown averaged over 5
simulation runs with different random realizations of the initial conditions).
As initial conditions, we take $\phi_{1,2}$ to have a flat power spectrum with
total variance $4f_{a}^{2}$ centered around $\Phi=f_{a}/\sqrt{2}$, with the
fluctuations cutoff at comoving momentum $k$ such that $k/(R_{i}f_{a})\simeq
0.4$ where $R_{i}$ is the scale factor at the initial simulation time (when
$H=f_{a}$). Such initial conditions lead to an initial under-dense string
network with $\xi\ll 1$. We see that at the initial simulation time
$\Delta^{2}_{1-\cos\theta}$ is almost flat (this is because a flat power
spectrum for $\phi_{1},\phi_{2}$ corresponds to a slighly tilted one for
$1-\cos\theta$). However, a suppression a large scales quickly develops,
already present by $\log(m_{\rho}/H)\simeq 2$, which we interpret as being due
to the radial mode overshooting its potential (the string density at such
times is still small, with most Hubble patches not containing a string). Note
that the simulation is evolved in a background of radiation domination, so for
the initial conditions we choose a sizable fraction of the radial mode