Early vs late string networks from a minimal QCD Axion

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 $O(10^8)$ GeV, corresponding to axion masses up to $10^{-2}~{\rm eV}$, and with miniclusters with masses as large as $O(10)M_\odot$.


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 ≃ H (with H the Hubble parameter), domain walls form [29].If the domain wall number N DW equals unity, the strings 'chop up' and quickly destroy the resulting network.If N 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 DW vacua. 1 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 DW > 1, but in this case a low scale of inflation H I ≪ 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 max after inflation [37] is large enough, T max ≳ 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 ≳ f a ≫ T 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 ≳ 2π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 λ 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 max ≪ 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. 2 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 DW = 1 the full observed dark matter abundance can be obtained for axion decay constants as small as allowed by astrophysical constraints, f a ≳ 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 Φ described by the Lagrangian The axion a(x) is the phase of Φ(x) = ρ(x) √ 2 e ia(x)/fa , while the radial mode ρ ≡ √ 2|Φ| has mass m 2 ρ = 2λf 2 a .We allow the quartic to be somewhat but not too much smaller than O(1), e.g.λ ≳ 10 −4 (the complications for smaller λ 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 postinflationary 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. 3 .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. 4

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 [65,66]).Throughout this time, scalars with mass smaller than √ 2H I experience de Sitter fluctuations similar to thermal fluctuations with effective temperature T dS = H I /(2π) [67,68].Once a mode exits the horizon it becomes classical with amplitude H I /(2π), giving a kick to average value of the field.Consequently, the field value in a coarse-grained region of size H −3 I 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 (Φ) of the coarse-grained field sampled over a region of space that contains many Hubble volumes.
The general form of the FP equation reads where t is cosmic time and ∂ i are the derivatives with respect to the fields ϕ i that are light during inflation.In general, eq. ( 2) admits a time-independent 'equilibrium' solution on which the spread of the field by inflationary fluctuations is balanced by the classical motion.A is a normalisation constant such that dϕ 1 dϕ 2 P = 1.
For the Lagrangian in eq. ( 1), the FP equation is two-dimensional with Φ = (ϕ 1 + iϕ 2 )/ √ 2 and is valid as long as m ρ < √ 2H 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 [69]).On the equilibrium solution, the complex field is uniformly distributed in the angular direction and For α ≪ 1 the radial mode is strongly peaked in the vacuum, i.e. the expectation value ⟨ρ 2 ⟩ = f 2 a .Instead for α ≫ 1 the radial mode is displaced [47]: with the dynamics dominated by the quartic term.In Figure 1 (left) we show ⟨ρ 2 ⟩ as a function of α.
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,70], we define the reduced probability distribution P ≡ e 4π 2 V 3H 4 which satisfies It is useful to express this in terms of the radial and angular modes ρ and θ ≡ a/f a to obtain, We can take advantage of the U(1) symmetry and look for solutions of the form P ∝ ψ m (t, ρ)e imθ , with m integer, given the periodicity of θ.This leads to In this form, the right hand-side can be interpreted as the Hamiltonian of a one-dimensional quantum mechanical particle with potential V eff .In terms of x ≡ λ 1/4 ρ/H I and α, the corresponding operator is O m has always one zero eigenvalue, Γ 10 , corresponding to the equilibrium solution of eq.(7).The other eigenvalues (in units of √ λH I ) are positive and discrete, and will be denoted by Γ nm with n > 1 integer.Γ nm and the corresponding eigenfunctions ψ nm can be found numerically. 5The most general solution of the Schröedinger-like equation ( 7) is then where the coefficients a nm are set by the initial conditions and we take ψ nm normalized in such a way that which follows from O m Peq = 0 and Peq = e −v ), the probability distribution is found, This shows that Γ nm determine the evolution of P , which tends to P eq asymptotically.The latetime behavior is regulated by the smallest non-zero eigenvalue, corresponding to (n, m) = (1, 1).In particular, 1/Γ 11 measures the cosmic time on which equilibrium is approached, i.e. this happens after a 'critical' number of e-folds of expansion and in this limit Intuitively, for α ≫ 1 the variance of P initially grows as σ 2 = N H 2 I /(2π) 2 so the equilibrium distribution in eq. ( 5) takes N ∝ 1/ √ λ e-folds to be reached, independent of H I /f a . 6Because Γ 11 is an increasing function of α, eq. ( 15) gives a lower bound for the number of e-folds to reach equilibrium, with ρ light.
Importantly for what follows, N 11 is also roughly the number of e-folds that it takes for Φ to be fully randomized; more precisely, on scales much larger than d 11 ≡ e N 11 /H I at the end of inflation, a(x) acquires random values in the interval [−π, π], 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 Φ 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 ≳ 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 This is obtained by neglecting the effect of the PQ potential during inflation, i.e. we take the real and imaginary parts of Φ as Gaussian random fields with flat power spectra (this is a good approximation if λ ≪ 1 but not accurate for large λ, in which case Φ 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 ≃ 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.
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 λ ≳ 0.05, when α ≫ 1.
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 ϕ 1 and ϕ 2 fluctuations are added to an initial homogeneous field value, arbitrarily set to ϕ 1 = f a and ϕ 2 = 0, as Gaussian fluctuations with a flat power spectra of amplitude H I /(2π) (i.e.root mean square fluctuations in the region of √ 6H I /(2π)).This reproduces the expectation from inflation with each point corresponding to a region of size H −3 I at the end of inflation.We use H I /f a = π, 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 ≃ 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 ≃ 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, Φ has a value Φ 0 typical of the equilibrium distribution, i.e.P = δ(ϕ 1 − √ 2Φ 0 )δ(ϕ 2 ).For α ≪ 1 this is reasonable, while for α ≫ 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 α = H I /(λ 1/4 f a ) and λ, starting from P = δ(ϕ 1 − √ 2Φ 0 )δ(ϕ 2 ), with Φ 0 given by the standard deviation on the equilibrium distribution, Φ 0 = ⟨ρ 2 ⟩ /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 α ≃ 1 were given in Ref. [45]).For λ = 1 and H I = f a , α = 1 and N 11 ≃ 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 λ = 1 and H I = 2f a , α = 2 and N 11 ≃ 15 so after 25 e-folds P is part-way to its equilibrium distribution, on which it is non-zero at Φ = 0 with ρ 2 ≃ f 2 a still; and for H I = 5f a , α = 5 and N 11 = 9 so the equilibrium distribution is approximately reached and on this P is unsuppressed at Φ = 0 and the expectation value ρ 2 ≳ f 2 a is slightly displaced.Meanwhile, for λ = 10 −3 , α ≫ 1 and N 11 ≃ 300 for all H I and Φ ≃ 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 2 I /(4π 2 ) ≃ 0.6H 2 I , with the radial mode initially displaced from the origin by ρ 2 ≃ 7H 2 I .

Evolution soon after inflation
We now consider the dynamics of Φ after the end of inflation when reheating begins.If N 11 ≫ 1, Φ 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. 7If ρ ≲ f a (e.g. for H I /f a ≲ 2 and λ = 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 ρ = 0, that will become the cores of strings).Conversely if α = H I /(λ 1/4 f a ) ≫ 1 the radial mode takes values ≫ 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.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 ∝ t 2/3 .We also assume that the inflaton decay rate is sufficiently small that the maximal temperature T 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 , λ 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 λ can be eliminated by defining √ λt = τ , and reheating starts at τ 0 = 2 √ λ/(3H I ).The motion is one dimensional if we simply allow for negative values of ρ (otherwise θ changes from 0 to π and the field crosses the top of the potential).The final value ρ 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 ρ i . 8We 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 ρ, which we call overshooting, is important because it can at least partially randomise Φ 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 (Φ(x)) ≃ arg (Φ(0)) and |Φ(x)| − |Φ(0)| ≪ |Φ(0)|.In the absence of overshooting, the value of Φ at these points would remain strongly correlated.However, if the small difference between |Φ(x)| and |Φ(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(Φ(0)) and arg(Φ(x)) will differ by π 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(Φ) over a two-dimensional spatial slice for a realisation of the inflationary initial conditions.As in Figure 2, we generate ϕ 1 and ϕ 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 = π/2, and unlike Figure 2, we assume Φ 0 = 2.3f a , corresponding to ⟨ρ 2 ⟩ in eq. ( 5) for λ = 10 −2 .These values are such that Φ does not fluctuate over the top of its potential during inflation.As a result, in the left panel, which shows arg(Φ) immediately at the end of inflation, only small inhomogeneities are present.However, in the right panel we plot arg(Φ) after the evolution to the minimum and overshooting, still assuming λ = 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(Φ) is close to either 0 or π, with relatively few points in-between.

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  2, except for H I /f a = π/2 and with the crucial difference that the initial value of the PQ field is taken to be Φ 0 ≃ 2.3f a .We generate the fluctuations of ϕ 1 and ϕ 2 neglecting the effect of the potential V , but note that the Φ 0 we choose corresponds to λ = 10 −2 assuming |Φ 0 | given by the expectation value in the equilibrium distribution eq.( 5) (λ 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 O(1) string per Hubble volume when ≃ 5 e-folds of inflation re-enter the horizon.
point. 9We define the string density ξ at a given cosmological time t after the end of inflation as where ℓ(V) is the length of strings in a volume V, such that ξ measures the string length in Hubble lengths per Hubble volumes.On the attractor scaling solution ξ is order-one, up to corrections ∝ log(m ρ /H) [31,[72][73][74][75], while for an underdense network ξ ≪ 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 ξ ≃ 1 to be reached; we say that at this time the string network has re-entered the horizon.We expect that the value of ξ, at the time when N e-folds of inflation have re-entered the horizon, is related to what extent Φ has spread over the top of the potential during the last N e-folds of inflation.Intuitively, if Φ 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 Φ is mostly concentrated on one side of the potential, we expect ξ ≪ 1, with strings only in a few rare Hubble patches.In the latter case, ξ will subsequently increase as more e-folds of inflation re-enter the horizon, until ξ ≃ 1 is reached once the critical number of e-folds have re-entered.
In order to quantify the string density ξ 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 Here P is the solution the FP equation after N inflationary e-folds starting from the initial condition . 10 In eq. ( 18), B accounts for overshoot and is the step-like function plotted in Figure 4 that takes values ±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 [76,77] for related analysis in the context of biased initial conditions).Note that for α ≪ 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 , λ of phenomenological interest. 11A 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 ≃ 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 ξ at the time when N efolds have re-entered the horizon is indeed related to F (N ), for the N ≲ 8 that can be simulated.In particular, we will see that ξ ≃ 1 is obtained for 0.3 ≲ F ≲ 0.45.We therefore postulate that ξ ≃ 1 is obtained when F (N ) ≃ 0.4 even when this is reached for much larger N ≲ 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 , λ, 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).

Inflationary formation
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 ρ .This amounts to neglecting overshooting by setting B = 1 in eq. ( 18).
Inflationary production . 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 ρ > √ 2H 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 Φ 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.
To calculate the fraction F , we use eq.( 18) and solve the FP equation for P .Assuming that the field is initially at ϕ 1 = √ 2|Φ 0 | and ϕ 2 = 0, corresponding to x 0 = λ 1/4 √ 2|Φ 0 |/H I with Φ 0 sampled from the equilibrium distribution, we get that the coefficients for P in eq. ( 11) are We can now ask how the probability distribution evolves in field space over time.The fraction F is simply obtained by integrating P over ϕ 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 This means that the number of inflationary e-folds required for a fraction F to be reached is Similarly to N 11 , N (F ) is proportional to H I /Γ 11 and thus is controlled by the lowest eigenvalue Γ 11 , even though the coefficient is sensitive to the fraction F .Our expectation that F determines the string density and the result N (F ) ≃ N 11 for F ≃ 0.5, are consistent with the intuition in Section 2.1 that a string network reaches ξ ≈ 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. 12Using 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 Φ 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 λ 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 /Γ 11 gives similar results.From Figure 1 (left), for α ≲ 2, corresponding to the region H I /f a ≲ 2λ 1/4 below the dotted line, the equilibrium distribution is peaked at the radial mode's vacuum and in this regime Γ 11 /H I is approximately independent of λ.Thus, the contours are horizontal at H I /f a ≃ 9/ √ N s and larger H I /f a simply results in Φ reaching equilibrium after a smaller number of e-folds.Conversely, for α ≫ 2, λ 1/2 H I /Γ 11 does not change as H I /f a increases at fixed λ, because the radial mode is displaced at ⟨ρ 2 ⟩ ≫ f 2 a , see Figures 1,3.This leads to vertical contours at λ ≃ 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.

Overshoot mechanism
If P does not approach the equilibrium distribution within the visible number of e-folds of the Universe (e.g.λ ≲ 0.05, for H I ≫ f a ), naively no string network would form. 13However, 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 0.3 ≲ F ≲ 0.45, a string network with ξ ≃ 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 λ 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 ̸ = 0. We use the overshoot function in Figure 4 and, as in the case of inflationary production, we set as initial condition Φ 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π) 2 .Additionally, since for the relevant H I /f a and λ the field does not reach the asymptotic distribution during visible inflation, other initial conditions (i.e.Φ 0 ) are possible but perhaps less plausible.
At a fixed λ, 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 Φ 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 ∝ λ 1/4 (since H/λ 1/4 controls the equilibrium field displacement used in Φ 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). 14As a result, these regions of the plot should be taken as a indication of the values of H I , λ 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 ≳ 2λ 1/4 , i.e. above the dotted line, where inflationary production is inefficient, see Figure 8 (left).
Finally, we note that Ref. [79] finds that if the PQ field reaches values |Φ| ≫ 10 4 f a strings can form after inflation as a result of parametric resonance (see also [80][81][82][83][84][85][86][87][88]).Such field values can occur with Φ roughly at its equilibrium value eq. ( 5) for λ ≪ 10 −2 (much smaller than those we consider), or it could simply be assumed that at the start of the visible era of inflation Φ 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 Φ's large oscillations that are dominantly amplified.Such dynamics are relevant in complementary parts of parameter space to those we focus on.

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 Φ 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 [89].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 −3 I at the end of inflation.Given that we can only simulate grids with n 3 ∼ 2048 3 lattice points, our entire simulation volume corresponds to N ∼ 8 e-folds of inflation, much smaller than the regime 25 ≲ N ≲ 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 ρ for all the values of H I and λ used in generating the initial conditions.Moreover, we evolve the mass of the radial mode throughout the simulations as m ρ ∝ 1/R(t) where R(t) is the scale factor, by using a time dependent λ.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 ρ ∆ = 1 is sufficient to resolve the string cores where ∆ 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 ξ ≪ 1 (under-dense) regime, ξ 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 ξ 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 ρ (t)/H(t) ≃ 10).We calculate the length of strings in the simulations, and correspondingly ξ as defined in eq. ( 17), using a standard algorithm [90].The simulations can be run until HL ≃ 2, where L is the physical box size, at which point finite volume systematic errors become significant.
We extract ξ 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 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. 15As 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 λ constant).

Results
We carry out simulations for a range of initial conditions (i.e.H I /f a , λ, and value of Φ 0 at the start of the ≃ 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 ξ 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, Φ 0 = f a / √ 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 ξ at F corresponding to the final time of some of the simulations match those at early simulation times in others.This suggests that ξ 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, Φ 0 = {3f a / √ 2, 5f a / √ 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 λ = 0, but perform the subsequent evolution after inflation with non-zero λ.)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  18), which measures the size of the inflationary fluctuations that have reentered the horizon at the time when ξ 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 λ = 0 is the value used to generate the initial conditions, but the field was subsequently evolved with non-zero λ.In both cases, ξ is approximately determined by F regardless of the value of H I /f a and λ.
the left panel.We note that the values of F that correspond to a particular ξ 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 ξ as a function of Hubble parameter in simulations is ξ = 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 ξ ∝ log(m ρ /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 ξ ∝ 1/R 2 ∝ 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 ξ ≃ 1 and scaling is reached. 16f course, there are limitations to our results.For each λ, 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 ξ, 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 ξ ≃ 0.5 after ≃ 4 − 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.

Isocurvature perturbations
The inflationary fluctuations of the complex field Φ 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 Φ 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 δ a (⃗ x) ≡ (ρ a (⃗ x) − ρa )/ρ a where ρ a = 1 2 ȧ2 + 1 2 |∇a| 2 + 1 2 m 2 a a 2 is the axion dark matter energy density once any strings present have been destroyed and the axion field is non-relativistic; ρa is its spatial average.The power spectrum ∆ 2 δ (k) of the axion overdensity field is defined by where δa is the Fourier transform of δ a .We define the component of ∆ 2 δ that comprises isocurvature fluctuations to be ∆ 2 iso .Meanwhile, the component of ∆ 2 δ consisting of the almost scale-invariant spectrum of curvature perturbations that is supposed to be originated by the inflaton fluctuations is There are strong observational constraints on isocurvature fluctuations on cosmological scales.Planck cosmic microwave background (CMB) data [91] bound the fraction of dark matter isocurvature fluctuations relative to the curvature perturbations.For a ∆ 2 iso (k) that is approximately constant in the interval k = [0.002,0.1]Mpc −1 one finds r iso ≡ ∆ 2 iso (k CMB )/∆ 2 R (k CMB ) ≲ 0.04.

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 [78]) H I /f a ≳ λ 1/4 and sufficiently small λ, such that the radial mode is displaced to field values ≫ f a , and that the equilibrium solution is not reached during visible inflation.As discussed, for such theories the value of the complex field Φ 0 prior to the visible e-folds of inflation cannot be predicted but we continue to make the reasonable assumption that ρ 2 0 ≃ ρ 2 ∼ H 2 I / √ λ, where ρ 2 is the expectation value on the equilibrium distribution, eq. ( 5).In this case after the visible e-folds of inflation, Φ has a nearly constant classical value Φ 0 = ρ 0 e iθ 0 / √ 2 in our Hubble patch (more generally, this is the case unless One finds that immediately after inflation [47], where θ = Arg(Φ), δθ = θ − θ, θ is the spatial average over the entire visible universe and N ≈ 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 Φ 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 ρ a (x) ∝ θ 2 (x), appropriate for θ well away from π, and (for the moment) considering only misalignment production of axions, for small λ we obtain that the isocurvature part of the power spectrum is approximately flat with an amplitude (at the end of inflation) [92,93] ∆ 2 iso = where in the last equality we used eq.( 5).Comparing with isocurvature constraints, naively this would exclude such a scenario unless λ is exceedingly small, λ ≲ 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 λ are such that at the end of visible inflation the probability distribution of Φ is approximately given by the equilibrium solution, i.e.N 11 ≲ 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 Φ.Following [42,70], in the limit α = H/(λ 1/4 f a ) ≫ 1 one finds at the end of inflation,17 where R e is the scale factor at the end of inflation and 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 ∆ 2 iso considering only misalignment production and approximating that the axion abundance ρ a ∝ 1 − cos θ (which is accurate for θ ≪ π, has the correct periodic behaviour, and is easily evaluated).This results in the overdensity power spectrum In the limit H I ≫ f a (with λ ≳ 0.01 so that N 11 ≲ 50, c.f. eq. ( 15)), using eq.( 14) the spectrum tilt is β = 0.22 √ λ.This can be compared with the approximate result in the mean-field approximation, described in Appendix A, eq. ( 52), which numerically gives α = 0.25 √ λ.The isocurvature power spectrum in eq. ( 27) at CMB scales can be written as where N vis is the number of e-folds of visible inflation.This implies that, for N vis ≃ 50, under the preceding assumptions the Planck constraint [91] would be satisfied only for N 11 < 5, which in turn would translate to λ ≳ 1 from eq. ( 15).However, as we have seen, for N s ≲ 25 a string network reenters 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.

Effects of strings on correlations
As mentioned, the large-scale correlations in the phase of Φ 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 ξ ≃ 1 before the QCD crossover, axions are continuously emitted from the network to sustain the subsequent scaling regime. 18Provided 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 ∆ 2 iso ∝ (ρ mis /ρ string ) 2 , where ρ mis and ρ 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 ξ ≃ 1 at the QCD crossover, such misalignment production cannot be clearly isolated from the string/domain network, and, relatedly, θ is not necessarily constant on large scales during the evolution.The technical reason for this is that in the string cores θ is not defined and so θ is not a conserved quantity on large scales. 19(Instead, in absence of strings and with the radial mode in its vacuum everywhere the zero mode decouples from the equations of motion and θ is a conserved quantity on very large scales.)As a result, a string network in scaling could erase long-distance correlations in θ.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 ξ ≪ 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 ≃ 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.

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 θ.In particular, for distances |x| over which the probability distribution of Φ(x) immediately at the end of inflation covers many of the "steps" in Figure 4, overshooting results in cos θ(x) taking an effectively random value from {cos θ 0 (x), − cos(θ 0 (x))}, where θ 0 (x) is the value before overshooting.Consequently, ⟨cos θ(x) cos θ(0)⟩ can average towards zero, even if θ 0 (0) and θ 0 (x) are highly correlated.
In reality, cos(θ(x)) and cos(θ(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 Φ with inflationary initial conditions for H I , λ such that ρ 2 ≫ f 2 a .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 where N = log(|x|H I ) is the number of e-folds of inflation that separates 0 and x (results from simulations showing a suppression can be found in Appendix D).Such a dependence is plausible: as expected, the correlation vanishes if Φ is evenly distributed over both sides of the potential after overshooting (i.e.F ≃ 0.5), and for a correlation to remain neither Φ(0) or Φ(x) can be "randomised".

Landscape of axion scenarios
We now turn to discuss phenomenology.Intuitively, theories such that a string network re-enters the horizon, i.e. reaches ξ ≃ 1, sufficiently early are indistinguishable from the thermal axion postinflationary 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.

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 λ, a critical number of e-folds N s that need to re-enter the horizon for ξ ≃ 1 to be reached.This corresponds to F (N s ) equal to some critical value, which we assume to be ≃ 0.35. 20 Defining, as in eq. ( 28), N vis to be the total number of visible e-folds of inflation, the relevant modes exited the horizon 20 In practice, in our plots we approximate Ns = Min[N inf.s , N o.s.

s
] from the two sources, which is a good approximation over the majority of parameter space where one or other formation mechanism dominates.
from the beginning of visible inflation.N PQ has a simple relation to the temperature of the Universe when the mode re-enters the horizon.Solving H 0 e N PQ = RH in radiation domination one finds, 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 P Q = 1 , 0.15 , 0.001 GeV exited at N 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 P Q with the temperature T ⋆ when H(T ⋆ ) = m a (T ⋆ ), around which time the axion potential becomes cosmologically relevant and the axion starts oscillating.This reads where m a (T ) = m a (T = 0)(Λ 1 /T ) αm/2 , with Λ 1 ≃ 150 MeV, Λ = 75 MeV, α m ≃ 8, and m a (0) = Λ 2 /f a .Note that N 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 ∼ f a ∼ 10 10 GeV.We also assume that after inflation the cosmological history is simply (relatively slow) reheating followed by radiation domination, in which case N vis depends logarithmically on the reheating temperature T RH , c.f. eq. ( 28).As a reference value, we take N vis = 50.

Regimes
Depending on the value of N s different scenarios are realized: If N s > N 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 ≳ f a , this regime is excluded by isocurvature perturbations [92][93][94][95][96][97][98] (unless any of the following happens: (1) λ ≲ 10 −19 from eq. ( 24), (2) the initial condition for Φ 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 vis for any λ ≲ 0.05 since the PQ field is approximately constant in our Hubble patch (at least for our assumed Φ 0 ), see Figure 6 (left) and [78].However, we have shown that due to the overshoot mechanism a network may still form when λ is small, Figure 6 (right), so the region of small λ can still be viable, falling into one of the two subsequent classes rather than the present one.
N s < N vis − 25: Restored PQ symmetry For such N s , at the end of inflation Φ 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 ≃ 25 e-folds after the start of visible inflation).This regime leads to predictions that are similar to the axion post-inflationary scenario.
. 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 ≲ 2λ 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 ρ gets an expectation value ≳ 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 Φ is not applicable.
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 ≲ 10 10 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 [31,[99][100][101][102][103][104]).The gravitational waves (GWs) emitted by strings are likely to be unobservable in this scenario for f a ≲ 10 10 GeV [30].
In the region of parameter space such that N vis − 25 < N s < N 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 ξ ≃ 1, only when the temperature of the Universe ≲ GeV (more precisely, when T P Q < T ⋆ , 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 λ, 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 dark white regions (labelled by '25 < N s < 50' and 'N s > 25?' respectively) the annihilation of the string/domain wall network takes place at T ≲ 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.

The late string relic abundance
In the late string regime, N vis − 25 < N s < N 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 P Q , related to N s via eq.( 31).
The former arises when the axion mass becomes cosmologically relevant at H ≃ 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 (−π, π]).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 ̸ = π (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 (−π, π], so the resulting averaged relic abundance is [ where Ω dm the observed dark matter abundance and κ is a constant that depends on α m , with κ ≃ 5 for α 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 ≪ 10 12 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 PQ < T ⋆ .For T PQ ≪ T ⋆ , 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 P Q when strings reenter the horizon and the other scales, H P Q ≪ m a ≪ f a .Nevertheless, we can make a reasonable estimate of the axion relic abundance as follows.
• If T PQ ≳ 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 where ρ w (T PQ ) ≃ 2Aσ(T PQ )H PQ is the energy in domain walls at T = T PQ and σ = βm a (T PQ )f 2 a (with β ≃ 9) their tension, ρ r the energy density in radiation and T eq the temperature at matterradiation equality (MRE).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 PQ ≳ 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 PQ ≲ 150 MeV the string-wall system decays when the mass has saturated to its zero temperature value and This estimate assumes that the axions emitted are non-relativistic otherwise the contribution is more suppressed.Note that as the network re-enters the horizon at lower temperatures T P Q , the abundance Ω dw a 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 ≳ 10 8 GeV [106][107][108][109][110][111][112][113], which corresponds to T PQ ∼ 1 MeV, around Big Bang Nucleosynthesis.
In Figure 9 we show the different behaviours and relic abundance as a function of f a and λ, 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.
The dynamics of the late string-wall decay could also source stochastic GWs.Despite the uncertain dynamics, we can estimate the corresponding spectrum with the standard quadrupole approximation, similarly to theories with N W > 1, see e.g.[114].At T = T P Q the walls decay, emitting GW energy density at momentum k ≃ H P Q at a rate dρ gw /dt ≃ Gρ 2 w /H P Q for approximately one Hubble time.This leads to a peak amplitude today . Different axion scenarios, and the corresponding dark matter abundance, as a function of λ 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 ≳ 1 GeV so that standard post-inflationary phenomenology is expected.In blue the string-wall network re-enters the horizon and annihilates at temperatures ≲ 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 ≲ 10 8 GeV is excluded by astrophysical constraints.
where we used eq.( 35) and assumed T P Q ≲ 150 MeV.The corresponding frequency today is f peak ≃ 2 • 10 −11 Hz (T PQ /1 MeV).From eq. ( 36), the GW spectrum is the largest for the smallest T PQ , and is a few orders of magnitude smaller the sensitivity of near-future detectors (Ω gw ≃ 10 −15 ) at the relevant frequency range (f ≳ 10 −9 Hz).We stress however that eq. ( 36) is only an approximate parametric estimate, and the true GW emission could be substantially different.

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 λ 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. λ ≃ 1/(R P Q H PQ ) (here R PQ is the scale factor when T = T P Q , and the network re-enters the horizon).This amounts to a contribution to the power spectrum of the axion overdensity field, ∆ 2 δ defined in eq. ( 22) (specifically ∆ 2 iso ), that is peaked at the momentum k/R PQ ≃ H PQ , where it acquires an order-one value.For momenta k/R PQ ≪ H PQ , this contribution to ∆ 2 iso is expected to fall ∝ C[k/(R P Q H 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 [49,115] 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 δ ≳ 1 undergoes nonlinear collapse when the scale factor is R ≃ R eq /δ, where R 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 ρ b ≃ 140δ 3 (1 + δ)ρ c ≃ O(100)ρ c , where ρ c is the average dark matter density at the time of collapse [32][33][34][35].Since the collapse happens around matter radiation equality, ρ b ≃ 100 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 where x ≃ O(1) parameterizes the typical physical size of the fluctuations at R = R PQ as R PQ λ = 2π/(xk PQ /R PQ ).Owing to the larger size of the fluctuations (related to the larger Hubble distance at the time of string-wall decay), for T PQ ≪ 1 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 [116][117][118][119][120][121].

Isocurvature constraints
As mentioned, the axion relic abundance from the late decay of domain walls is expected to have a white noise isocurvature spectrum ∆ 2 iso (k are not expected to inherit the correlations from inflationary fluctuations, eqs.(24) or (27).Such white-noise isocurvature perturbations are constrained by CMB [91] and Lyman-α observations [122][123][124][125][126][127][128], which bound the fraction of isocurvature fluctuations r iso ≡ ∆ 2 iso (k CMB )/∆ 2 R (k CMB ).In the late string-wall regime we have Lyman-α observations require r 1/2 iso ≲ 4 • 10 −3 for a white-noise power spectrum [124,125], 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 θ 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 superhorizon 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 0.3 ≲ F ≲ 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, Considering inflationary production, we expect N s ≃ N 11 up to a numerical factor that we do not know reliably (e.g.assuming a critical F = 0.35, N s ≃ 1.5N 11 , c.f. eq. ( 21) where a 11 b 11 ≃ −0.65).
Consequently, in the late string regime with 25 ≲ N s ≲ 50, we expect the suppression from the exponential factor in eq. ( 39) to be at least two orders of magnitude.(For overshoot production N vis ≪ 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 Ω mis a /Ω dw a ≃ 10 −5 for f a ≃ 10 8 GeV.Nevertheless, eq. ( 39) remains a potentially strong constraint given that ∆ 2 iso (k CMB ) ≲ 10 −10 from Planck data [91] (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.

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 λ.
We focused on the scenario where the PQ symmetry is not restored by thermal effects, T max ≲ f a .According to the common lore, for H I /(2π) ≳ f a the PQ symmetry is restored during inflation while for H I /(2π) ≲ f a it is broken throughout the cosmological history of the universe. 21Very 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 preinflationary scenario would be realized.
Clearly this picture must be incomplete and at least in the intermediate region H I /(2π) ≃ 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 λ of the PQ complex field.Roughly the landscape is the following: 1.If the quartic is sizable, say λ ≳ 0.1, and H I ≳ 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. For smaller λ ≲ 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 ≲ 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.Φ 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 ≳ 10 8 GeV (additionally, possible gravitational wave signals and primordial black holes resulting from similar dynamics have been considered in Refs.[130,131]).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 [132].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 [133].
3. For even smaller λ 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 λ.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 ≳ 1 GeV or at T ≲ 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 λ is challenging.We have approached this analytically by studying the evolution of Φ 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 , λ), 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 ≃ 2000 3 lattice points that we have the computational power for, due to the multiple, widely separated, scales in the problem (H P Q ≪ m a ≪ f a , where H P Q 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 ≲ 1 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 postinflationary 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 ≃ 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 ≲ 10 −2 eV, such as IAXO [134,135] among others [136][137][138].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.[139][140][141][142][143] and references therein) detectable as an 'effective additional neutrino species', ∆N eff , via forthcoming CMB experiments such as CMB-S4 [144] and Simons Observatory [145] and in combination with Large Scale Structure surveys (e.g.Euclid [146]).
A final important caveat is that in the regime H I ≫ f a λ 1/4 our results depend on the assumption that, prior to the visible e-folds of inflation, Φ 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 Φ will remain an important theoretical uncertainty in this part of parameter space.
The power spectrum of a free scalar field ϕ with effective mass M ≲ H I during inflation reads From this we can compute the field's variance at the end of inflation, 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. (41) shows that for small M ≪ H I / √ N I it is a good approximation to neglect the potential and the variance undergoes a random walk, while for M ≃ H I / √ N I the variance reaches a constant value within the observable number N I of e-folds.The fact that ⟨ϕ 2 ⟩ 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 ⟨ϕ 2 ⟩.
Let us first consider the quartic potential V = λϕ 4 /4.The effective mass is M 2 = 3λ⟨ϕ 2 ⟩ and, from eq. ( 41), in the stationary regime we thus have We can extend eq. ( 41) to several scalar fields, obtaining that the expectation value of ϕ i ϕ j is For the PQ field V = λ/4(ϕ 2 1 + ϕ 2 2 − f 2 a ) 2 so that, In the limit H I ≫ λ 1/4 f a , from eq. ( 43) with solution By plugging eq. ( 46) in eq. ( 44), and these in eq. ( 40), we can also compute the power spectra, 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 ⟨ϕ i ϕ j ⟩ during inflation.This is determined by the equation, where the first term corresponds to the stochastic noise due to inflationary fluctuations while the second is the classical drift.Eq. ( 48) 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. ( 46).Indeed, taking ⟨ϕ i ϕ j ⟩ = 0 and ⟨ϕ 2 1 ⟩ = ⟨ϕ 2 2 ⟩, eq. ( 48) simplifies to, and at equilibrium one recovers eq. ( 46).The stationary solution is reached equating the random walk N to the equilibrium value ⟨ϕ 2 1 ⟩ in eq. ( 46), which gives a number of e-folds N ∼ 3π 2 2λ .

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, In Fourier space, where we used the expression for ⟨ϕ 2  1,2 ⟩ in equilibrium in eq. ( 46) and for the power spectrum for ⟨ϕ 1 (x)ϕ 1 (0)⟩.The suppression at the CMB scales is thus, This has the same form of eq. ( 27) with

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 ξ ≃ 1 to be reached (at which point the string scaling regime starts if T P Q ≳ 1 GeV, or the string-wall network annihilates if T P Q ≲ 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 Inflationary production . 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.reach ξ = 1 to better than the range 0.3 ≲ F ≲ 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 λ 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 Φ 0 , so the additional uncertainty from F is not a major issue.

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 Φ = Φ 0 on all lattice points.At this point, the entire simulation volume represents a size H −3 I region N = log 2 (n) e-folds before the end of inflation (with the value of the field coarse grained over H −1 I ).• To reproduce one doubling time of inflation (during which the initial H −3 I volume grows to 8H −3 I ) we split the homogeneous region into 8 patches.Inflationary fluctuations are accounted for by giving ϕ 1 and ϕ 2 (the canonically normalised real and imaginary parts of Φ) in each of the new regions a 'kick' drawn from a Gaussian distribution with variance log(2)H 2 I /(4π 2 ) and mean 0.
Finally, we evolve Φ in each of the regions according to the classical equations of motion (for a homogeneous field) for a time log(2)H −1 I .• 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 3 I .
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 Φ 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 −1 I 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 ρ /H).As mentioned in the main text, the dynamics of the string network are only captured without significant systematics provided there is at least ∼ 1 lattice point per string core (which is of size m −1 ρ ) and HL ≳ 2 (where L is the physical size of the box).With the Lagrangian of eq. ( 1), the string core size decreases in comoving coordinates ∝ R −1 (where R is the scale factor).Consequently, to maximise the range in m ρ /H of such simulations, the string cores must initially contain √ n lattice points.In other words, m ρ ∆ = 1/ √ 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 1  2 log n e-folds of inflation re-entering the horizon.Instead, we simulate the so-called fat string system, i.e. with λ changing with time in such a way that m ρ ∝ R −1 , and we set the lattice spacing ∆ as m ρ ∆ = 1; we start the simulations at H = m ρ .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 ξ, statistical uncertainties on ξ become large prior to the maximum time allowed by the systematic uncertainties discussed above (even with ≃ 10 repeated runs).Our data sets include sufficiently many simulations that the relative statistical uncertainty on ξ 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 ρ /H) ≳ 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 , λ (used to generate the initial conditions), as well as the initial radial mode value ρ 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 ρ 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 ρ 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 ξ 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  7, but with the string density ξ plotted as a function of simulation time, parameterized by log(m ρ /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.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 ): where P (ρ) = 2π 0 dθ ρP (ϕ 1 , ϕ 2 )| ϕ 2 1 +ϕ 2 2 =ρ 2 is the probability distribution of ρ 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 (N ) | ρ=ρ i is the function defined by eq. ( 18) setting ρ = ρ i N e-folds before the end of inflation.For simulations with λ = 0 in the initial conditions, P is a Gaussian and we evaluate eq. ( 53) directly.These initial conditions lead to a field displaced by ρ 0 ≫ f a after inflation and should produce strings via overshoot, thus the issue discussed above is important to account for.For λ ̸ = 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 ρ is strongly peaked in the vacuum, which amounts to not carrying out the average in eq. ( 53), and is a good approximation for all (H I /f a , λ) 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 F ).

C.2 The string network
To clarify the properties of underdense string networks, we plot in Figure 11 the string density ξ as a function of simulation time (parameterised by log(m ρ /H)) rather than F .At a given simulation time the string density varies by a factor of ∼ 2 among initial conditions with different H I /f a and λ, and the rate of increase in ξ also varies.It is therefore non-trivial that that the values of ξ between different initial conditions coincide once ξ 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 ξ at a given log(m ρ /H) to test how well our proposed dependence on F holds.As mentioned in the main text, the increase in ξ with log(m ρ /H) is exponential rather than linear.This is characteristic of the network being underdense compared to the critical ξ such that string annihilation balances the rate at which strings reenter the horizon, as opposed to the logarithmic growth of ξ on the attractor as a result of this critical point slowly moving with log(m ρ /H).

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 ξ and the power spectrum of 1−cos(θ), ∆ 2  1−cos θ , at different cosmological times (parameterized by log(m ρ /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 ϕ 1,2 to have a flat power spectrum with total variance 4f 2 a centered around Φ = f a / √ 2, with the fluctuations cutoff at comoving momentum k such that k/(R i f a ) ≃ 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 ξ ≪ 1.We see that at the initial simulation time ∆ 2  1−cos θ is almost flat (this is because a flat power spectrum for ϕ 1 , ϕ 2 corresponds to a slighly tilted one for 1 − cos θ).However, a suppression a large scales quickly develops, already present by log(m ρ /H) ≃ 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 overshoots, c.f. Figure 4.The suppression seems to tend to the k 3 form expected if the axion field was completely randomised; a more detailed investigation would be useful.As an additional test, we have evolved the same initial conditions without the spatial gradients in the Hamiltonian.In such a system we again a similar suppression to ∆ 2 1−cos θ on large scales, providing further evidence that this is due to overshoot.
As a further test, we generate a similar realisation of inflationary initial conditions with a flat initial power spectrum with total variance of (4f a ) 2 centered around Φ 0 ≃ 4f a .Rather than evolving as in our field theory simulations, in this case we simply map the radial mode to its vacuum neglecting spatial gradients (and keeping the angular variable constant), i.e. using the step-like function plotted in Figure 4.The form of ∆ 2 1−cos θ before and after this mapping of the radial mode is plotted in Figure 13, assuming both the step-like function appropriate to matter domination and for radiation domination.A suppression of the super-horizon correlations by overshooting is evident.In the matter dominated case there is evidence that the the destruction of large-distance correlations is not complete, i.e. ∆ 2  1−cos θ tends to a non-zero constant as k → 0, corresponding to the fact that the right hand side of eq. ( 29) is not zero if F ̸ = 0.5 (a less clear hint of a similar remaining correlation is also visible in the radiation dominated case).The numerical value that ∆ 2  1−cos θ tends to for k → 0 is roughly compatible with eq. ( 29), given that F ≃ 0.25, 0.49 for matter domination and radiation domination respectively.1−cos θ (right) in a numerical simulation of the string network, starting from initial conditions that lead to an initially under-dense network.We interpret the suppression of ∆ 2 1−cos θ at large spatial scales, which occurs at early times in the simulation, as being due to overshoot randomising the phase of the complex scalar Φ. Results are shown immediately after inflation (blue).Additionally, in orange we show results after mapping the radial mode to its vacuum using the overshoot function plotted in Figure 4 assuming matter domination (left) and radiation domination (right).Note that these are not results from a full numerical evolution of the theories equations of motion, and instead account only for overshoot and neglect spatial gradients.The suppression of ∆ 2  1−cos θ at large scales is therefore due to the effects of overshoot.Numerically, in the case of matter domination the split of the radial mode on the two sides of its potential after overshooting is roughly 25% − 75%, and in this case a non-zero constant contribution to ∆ 2  1−cos θ survives to large scales.In the case of radiation domination, overshoot is more efficient at randomising Φ: there is a 49% − 51% split of the radial mode between the two sides of the potential and ∆ 2  1−cos θ is close to the k 3 white noise prediction for a completely randomized field.

Figure 1 .
Figure1.Left: The expectation value of the radial field squared ρ 2 /f 2 a on the equilibrium distribution as a function of α = H I /(λ 1/4 f a ).This is at the vacuum for α ≪ 1 and displaced by α 2 for α ≫ 1. Right: The combination λ 1/2 H I /Γ ij corresponding to the two smallest non-vanishing eigenvalues of O m , Γ 11 and Γ 20 , as function of α.The related quantity N 11 = H I /Γ 11 sets the number of e-folds on which the probability distribution of the complex scalar Φ approaches its equilibrium form during inflation.

Figure 2 .
Figure2.The axion field a(x)/f a ≡ Arg(Φ) over a spatial slice after 6 e-folds of inflation, with H I /f a = π and the field initially in the vacuum Φ 0 = f a / √ 2. This is obtained by neglecting the effect of the PQ potential during inflation, i.e. we take the real and imaginary parts of Φ as Gaussian random fields with flat power spectra (this is a good approximation if λ ≪ 1 but not accurate for large λ, in which case Φ 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 ≃ 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.

Figure 3 .
Figure 3.The probability distribution P (Φ) after N = 25 e-folds of inflation, starting from the initial condition P = δ(ϕ 1 − √ 2Φ 0 )δ(ϕ 2 ) setting Φ 0 = ⟨ρ 2 ⟩ /2, with ⟨ρ 2 ⟩ 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, λ = 1 and λ = 10 −3 .At a fixed λ, H I /f a controls both the initial field displacement and the spread of the probability distribution.Red circles indicates the vacuum |Φ| = f a / √ 2.

Figure 4 .
Figure 4. Solid black line: final value of the radial field, ρ f , as a function of the value at the end of inflation, ρ 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.θ = π when starting at θ = 0).Dashed red line: results obtained assuming a background of radiation domination.

Figure 5 .
Figure5.Left: Example of a/f a ≡ Arg(Φ) over a spatial slice at the end of 6 e-folds of inflation, similar to Figure2, except for H I /f a = π/2 and with the crucial difference that the initial value of the PQ field is taken to be Φ 0 ≃ 2.3f a .We generate the fluctuations of ϕ 1 and ϕ 2 neglecting the effect of the potential V , but note that the Φ 0 we choose corresponds to λ = 10 −2 assuming |Φ 0 | given by the expectation value in the equilibrium distribution eq.(5) (λ 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 O(1) string per Hubble volume when ≃ 5 e-folds of inflation re-enter the horizon.

Figure 7 .
Figure 7.The axion string density ξ in simulations starting from initial conditions corresponding to inflationary fluctuations with different value of the Hubble parameter during inflation, H I , the quartic coupling, λ, and radial mode field value at the start of inflation, ρ 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 ξ 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 λ = 0 is the value used to generate the initial conditions, but the field was subsequently evolved with non-zero λ.In both cases, ξ is approximately determined by F regardless of the value of H I /f a and λ.

(Figure 11 .
Figure 11.The same simulation results as in Figure7, but with the string density ξ plotted as a function of simulation time, parameterized by log(m ρ /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.

Figure 12 .
Figure 12.The evolution of the string density ξ (left) and the power spectrum ∆2  1−cos θ (right) in a numerical simulation of the string network, starting from initial conditions that lead to an initially under-dense network.We interpret the suppression of ∆ 2 1−cos θ at large spatial scales, which occurs at early times in the simulation, as being due to overshoot randomising the phase of the complex scalar Φ.

Figure 13 .
Figure13.The power spectrum of 1 − cos(θ) for a realization of inflationary initial conditions consisting of ϕ 1,2 having a flat power spectrum of (Gaussian) fluctuations of total variance of (4f a ) 2 (equivalent to a Hubble parameter during inflation of H I = 4πf a ) centred around Φ 0 ≃ 4f a , as a function of comoving momentum k/R.Results are shown immediately after inflation (blue).Additionally, in orange we show results after mapping the radial mode to its vacuum using the overshoot function plotted in Figure4assuming matter domination (left) and radiation domination (right).Note that these are not results from a full numerical evolution of the theories equations of motion, and instead account only for overshoot and neglect spatial gradients.The suppression of ∆2  1−cos θ at large scales is therefore due to the effects of overshoot.Numerically, in the case of matter domination the split of the radial mode on the two sides of its potential after overshooting is roughly 25% − 75%, and in this case a non-zero constant contribution to ∆2  1−cos θ survives to large scales.In the case of radiation domination, overshoot is more efficient at randomising Φ: there is a 49% − 51% split of the radial mode between the two sides of the potential and ∆2  1−cos θ is close to the k 3 white noise prediction for a completely randomized field.