More Axions from Strings

We study the contribution to the QCD axion dark matter abundance that is produced by string defects during the so-called scaling regime. Clear evidence of scaling violations is found, the most conservative extrapolation of which strongly suggests a large number of axions from strings. In this regime, nonlinearities at around the QCD scale are shown to play an important role in determining the final abundance. The overall result is a lower bound on the QCD axion mass in the post-inflationary scenario that is substantially stronger than the naive one from misalignment.

scaling violations are indeed manifest. 1 During the scaling regime axions are radiated from the strings, and if the properties of the network throughout this time are understood with sufficient accuracy the axions produced during this phase could also be reconstructed reliably. We should note that a huge extrapolation is still required to connect the ratio of scales that can be computed directly (slightly more than three orders of magnitude) to the physical ratio previously mentioned. However, the presence of the attractor, the fact that the scaling violations are only logarithmic and, as we will show, the fact that the final abundance is mostly determined by the qualitative features of the network, will allow us to perform such an extrapolation with some confidence.
The main inputs required for this programme are the total energy radiated from strings into axions during the scaling regime and the shape of the instantaneous axion spectrum emitted. Using energy conservation and the presence of the scaling law, the first quantity can be linked to one of the main parameters of the scaling solution: the average number of strings per Hubble patch ξ, which is, as we will discuss, a slowly varying function of time. Meanwhile, the spectrum is contained between an infrared (IR) cutoff set by the Hubble scale and an ultraviolet (UV) one set by the string thickness. The absence of any other scales in the problem suggests that, between these two cutoffs, the spectrum should be described by a single power law. The associated spectral index q determines whether the spectrum is IR or UV dominated, i.e. whether the energy of the radiation is distributed over a large number of soft axions (for q > 1) or a small number of hard ones (for q < 1).
Although the spectrum is mostly UV dominated in the range of parameters that can be reached by present simulations [7], we find clear evidence of a non-trivial running of the spectral index, which is more compatible with an IR dominated spectrum once extrapolated to the physical parameters.
These results imply that by the time the axion mass turns on the amplitude of background axion radiation produced by strings at previous times is large. In fact the occupation number of axions emitted by strings would be so large that nonlinear effects of the axion potential cannot be neglected, even considering the axion radiation in isolation, without topological defects.
We study the effects of these nonlinear dynamics in some detail. Their main consequence is a partial reduction of the number density of axions from strings, which however continues to dominate over the naive estimates based on the misalignment mechanism alone, or equivalently, over the results obtained by simulations of the full network of strings and domain walls carried out at the (currently available) unphysical values of the string thickness.
The article is structured as follows. We present our discussion of the most important points of the analysis and the key results in the main text, in Sections 2, 3 and 4. Meanwhile we give all the details of the various analyses, further studies, spin-off results, checks and generalization of the formulas, and interpretations in the Appendices. In particular, in Section 2 we present the results of simulations of the scaling dynamics and the axions produced by strings. In Section 3 we provide both analytical and numerical analysis of the effects of nonlinearities on the axion abundance from strings. In Section 4 we discuss the physical implications of the results and the assumptions and uncertainties behind them. In Appendix A we give details about the numerical simulations. In Appendix B we provide additional analysis of the properties of the string network during the scaling regime, including studies of string velocities, the decoupling of the heavy modes, the axion and radial mode spectra, as well as the systematics. In Appendix C we discuss how logarithmic effects are also visible in the dynamics of single loops in isolation. In Appendix D we identify when and how the scaling regime ends as the axion potential turns on. In Appendix E we give more details and results of both the analytical and numerical analysis of the nonlinear regime during the QCD crossover. In Appendix F we study the effects of the presence of topological defects during the QCD crossover on the evolution of the axion radiation produced during the scaling regime. Finally, in Appendix G we comment on the compatibility of our results with the existing literature.

Axions from Strings: The Scaling Regime
When the PQ symmetry is broken a network of axion strings forms [14][15][16] and this rapidly approaches an attractor solution [17][18][19] during the subsequent evolution of the Universe (extensive evidence for this was given in ref. [7]). The attractor is independent of the network's initial properties, allowing predictions to be made that are independent of the details of the PQ breaking phase transition and of the very early history of the Universe (i.e. at times much earlier than that of the QCD crossover).
The dynamics of the string network is highly nonlinear, and while models have been proposed to describe the main features of the attractor [9][10][11][12][13] they typically rely on a series of (unproven) assumptions. Instead we study the properties of the string network using numerical simulations. In these we integrate the classical equation of motion of the complex scalar field φ that gives rise to the axion numerically, assuming a radiation dominated Universe. 2 For simplicity we choose the Lagrangian which leads to spontaneous PQ symmetry breaking at the scale v. The axion field a(x) is related to the phase of the complex scalar field as φ(x) = v+r(x) √ 2 e ia(x)/v , while the radial mode r(x) is a heavier field of mass m r associated to the restoration of the PQ phase.
The scale v can be trivially reabsorbed in a rescaling of φ, while the scale m r provides the normalization of the physical space and time scales over which the dynamics unfolds. While the clock of the UV physics associated to the radial mode ticks with intervals set by 1/m r , the more phenomenologically relevant clock associated to the IR axion physics ticks at a much slower pace set by the scale 1/H, which keeps slowing down as the Universe expands. For this reason it is more useful to study the dynamics in terms of Hubble e-foldings log(m r /H) = log(t/t 0 ) ("log" for short), where t is the Friedmann-Robertson-Walker time (with metric ds 2 = dt 2 − R 2 (t)dx 2 , R(t) ∝ t 1/2 ) and t 0 is the reference time at m r = H. With appropriate random initial conditions, strings automatically form in simulations and their dynamics are fully captured. Regions of space containing string cores are identified from the variation of the axion field around small loops of lattice points. Further details on our implementation and algorithms are given in Appendix A.
Limits on computational power allow us to evolve grid sizes of at most 4500 3 lattice points. Meanwhile, the lattice spacing ∆ should be such that m r ∆ 1 and the physical length of the box L such that LH 1 to avoid introducing significant systematic uncertainties. As a result, simulations can only access relative small values of log(m r /H) ≈ log(L/∆) 8. In contrast, the vast majority of the axions present when the axion mass becomes cosmologically relevant come from the emission of the string network in the prior few Hubble times. This happens shortly before the time of the QCD crossover when log(m r /H) 60 ÷ 70. Therefore properties of the string network measured in simulations must be extrapolated if reliable, physically relevant, predictions are to be obtained.
A simulation trick, often used in the literature, is to make m r vary with time as m r (t) = m r (t 0 ) t 0 /t (the so-called "fat-string" trick). In this way the string thickness 1/m r (t) stays constant in comoving coordinates. The maximum log(m r (t)/H) = 1/2 log(t/t 0 ) that can be simulated remains unchanged, but this is reached over a much longer physical time, which allows far better convergence to the attractor regime. Although the simulations performed with this trick might lead to different quantitative answers, it is expected (and so far confirmed) that the qualitative behavior is the same. We performed most simulations with both m r constant ("physical") and with the "fat" trick. While we only use the data from the physical simulations to extract the relevant parameters, the results obtained with the fat trick make some features of the attractor solution more manifest and our interpretation of the string dynamics more robust.

String Density
The energy density stored in the string network can be written as where ξ, the number of strings per Hubble patch, counts the total length of the strings inside a Hubble volume in units of Hubble length, namely ξ ≡ lim L→∞ (L) t 2 /L 3 , while µ eff represents the effective tension of the strings, i.e. their energy per unit length. At late times, the latter is approximately equal to the tension of a long straight string in one Hubble patch µ = πv 2 log(m r /H), where v is the PQ breaking scale, which we take equal to the QCD axion decay constant f a from now on (we will discuss how to adapt our results to the more general case v = N f a in Section 3.3).
Such an approximation captures ρ s 's leading dependence on H and the UV parameters of the theory (f a and m r ). Corrections from the boost factors and the curvature of the strings are discussed in Appendix B.3.
The dynamics of strings are well known to be logarithmically sensitive to the evolving scale ratio m r /H. As mentioned above, the string tension is itself a linear function of this logarithm, and con- sequently the effective coupling of large wavelength axions with long strings scales as 1/ log(m r /H) (see e.g. ref. [20]). It is therefore not surprising that the dynamics of the string network, and in particular the parameters of the attractor, might depend non trivially on log(m r /H). This is indeed the case for the parameter ξ, which was observed to "run" in ref. [7] (see also refs. [21][22][23][24][25][26] for further supporting evidence), increasing logarithmically with time.
The growth of ξ is manifest in Fig. 1, which shows ξ as a function of log(m r /H). Each color refers to a set of simulations with different initial string density (initially overdense simulations show first a drop and then a universal increase). The error bars refer to the statistical errors. 3 Simulations ending before log = 7 are data taken in ref. [7] with grids up to 1250 3 , and the remainder are new data collected with bigger grids, up to 4500 3 . When we analyze other properties of the scaling solution we choose the initial conditions that reach the attractor behavior the earliest, indicated with black data points in Fig. 1. 4 Because of the manifest logarithmic increase, the value of ξ at late times could be much larger than that measured directly in simulations. In ref. [7] it was shown that the data is compatible with a linear logarithmic growth. Here we extend that analysis including all the data sets with different initial conditions and with bigger grids, in total comprising about 1000 simulations of which 100 are with grids larger than 4000 3 . We test the linear logarithmic increase with the 3 These take into account both the total number of simulations and the number of independent Hubble patches in each simulation. For this reason the error bars increase toward the end of simulations where fewer Hubble patches are available. 4 These are roughly those with the least overdense initial conditions. following fit ansatz (see Appendix B.1 for more details): where the coefficients c −1,−2 are taken with different values for each data set to account for differing initial conditions, while the coefficients c 1,0 , which survive in the large log limit, are taken universal across all data sets. As explained in [7] the string network starts showing scaling behaviors after log = 4 (when strings can begin efficiently emitting axions with sub-horizon wavelengths), which we choose as our starting point for the fit. 5 The result of the fit is represented by the colored curves in Fig. 1. The ansatz in eq. (3) reproduces all the data for a variety of initial conditions very well over almost 4 e-foldings in time.
The O(1/ log) corrections are relevant only at the smallest values of the logs in the fit, while they become almost irrelevant by the end of the simulations.
The fit value of the slope c 1 = 0.24(2) is definitely nonzero, confirming a non-vanishing universal increase. A straight extrapolation to log = 70 would give ξ = 15 (2). The current precision however does not allow us to exclude an even steeper growth. In fact, a fit with an extra quadratic term (i.e. c 2 log 2 ) gives analogously good results with a positive quadratic coefficient c 2 , which would lead to even bigger values of ξ at large logs. Simulations with the fat trick, which had more time to converge to the attractor, show an even more manifest linear log growth (see Appendix B.1). In particular the data set with initial conditions that reached the attractor the earliest in Fig. 6 leaves very little room for any nonlinear function to be a good fit. This suggests that ξ has a linear behavior in both the physical and fat systems, as opposed to a steeper growth.
Because of the decoupling of the axion field at large values of the log, continued growth of ξ beyond the reach of simulations would be compatible with the expectation that the global string network tends to approximate the Nambu-Goto string one (and the local string one) in the limit log → ∞. Indeed, old Nambu-Goto simulations gave values of ξ NG between 10 and 20 [18,19,27], while more recent local string ones [28] give ξ loc = 15(1). This is a hint that ξ for the global string network will not saturate at least prior to log = 60 ÷ 70 (extrapolating the linear growth).
An enhanced value of ξ was also observed in refs. [22,26] where a large value of the effective string tension was achieved by means of a clever modification of the physics at the string core m r .
However, we should point out that the asymptotic evolution of the string network parameter ξ for axion strings has not yet been fully established. It is still unknown whether the decoupling of the axion from the string dynamics really completes within a finite range of logs or keeps going with an infinite running. As we will see further below, the axion spectrum extracted from field theoretic simulations still shows nontrivial changes in the dynamics that could qualitatively affect the asymptotic behavior of the network. On the other hand, Nambu-Goto simulations could also miss the asymptotic behavior of the network, as they lack the back-reaction of the bulk fields and Kalb-Ramond effective descriptions might not capture the physics of string reconnections and backreaction of UV modes properly. In fact even for local string networks, which are expected to already be in the Nambu-Goto limit, a nontrivial logarithmic evolution of ξ might be present [29].
To summarize, while we cannot exclude the possibility that the observed growth of ξ saturates at larger values of the log, no indication of this is observed in the simulated range (it is particularly clear that the data for the fat system is incompatible with any reasonable function that plateaus soon after log = 8), which suggests that such a saturation could potentially happen only at much later times, if at all. 6 Instead, all approaches seem to agree on a growth of ξ to the range O(10) for log ∼ O(100), which is probably the most plausible and safe extrapolation. For our purposes we will assume the nominal value from our fit ξ = 15 for log = 60 ÷ 70, taking into account that this estimate might receive O(1) corrections.
Another quantity that characterizes the string network is the distribution of string velocities. We study this property in Appendix B.2 where we show that, in agreement with other studies [21,23,30], the strings are mildly relativistic with an average boost factor γ ∼ 1.3 ÷ 1.4. While this value appears to be approximately constant over the simulation time, the distribution of velocities shows a nontrivial evolution, with a subleading portion of the string network reaching increasingly higher boost factors as the log increases. This property is also compatible with the interpretation that the system is evolving towards the Nambu-Goto string network behavior, for which the formation of kinks and cusps explores arbitrary high boosts, and loops oscillate many times instead of shrinking and disappearing after one oscillation (more details are given in Appendices B.2 and C). As a consequence of the increasing Lorentz contraction from higher boosts, finite lattice spacing effects become more severe at larger values of the log. Such effects can be seen in a variety of observables (in particular in those that are more UV sensitive, see Appendix B for examples), and decrease the potential dynamical gain from simulations with bigger grids.

Axion Spectrum
In an expanding universe, eq. (2) and the conservation of energy imply that the string network continuously releases energy at a rate Γ ξµ eff /t 3 8πH 3 f 2 a ξ log(m r /H) (see e.g. [7] for more details). As shown in ref. [7], although most of this energy is emitted into axions, in simulations a non-negligible portion goes into radial modes (between 10% and 20%). Thanks to our new data with larger final logs, and by analyzing the radial excitation spectrum, we find that a significant part of the energy in radial modes is actually produced at the time the network enters the scaling regime and the subsequent emission of radial modes becomes less and less important (see the discussion in the Appendix B.4). This is compatible with the expectation that UV modes decouple from the network evolution at large values of the log (see Appendix B.3). 7 We will therefore assume that at late times the emission of radial modes is negligible and all the energy is released into axions.
The total energy density in axion radiation at late times is therefore ρ a 4/3µ eff ξH 2 log, where  Figure 2: The spectral index of the instantaneous axion emission q as a function of time (represented by log(m r /H)) for the physical (left) and fat string systems (right). The darker/lighter shaded regions correspond to the results of linear/quadratic fits to the data of the simulations (in black). The clear increase in q implies that the axion spectrum is turning IR dominated (i.e. q > 1), a regime that it will reach long before the physically relevant log(m r /H) = 60 ÷ 70.
the last log factor arises from the convolution of the emission rate over time. 8 As explained at length in ref. [7], the contribution of such radiation to the final axion abundance strongly depends on how the energy is distributed over axions of different momenta. A particularly useful quantity is the normalized instantaneous spectrum F (k/H) = ∂ log(Γ)/∂(k/H), which tracks the momentum distribution of axions produced at each moment in time by the string network. As mentioned in the Introduction, F is expected to be approximately a single power law F ∼ 1/k q between the IR scale set by Hubble and the UV one set by the string core. Depending on whether the spectral index q is greater or smaller than unity, most of the axion energy density emitted is thus contained either in a large number of soft axions or in a smaller number of hard ones, with obvious implications for the resulting number density. For example, if F is single power law 1/k q with compact support k ∈ [x 0 H, m r /2], the axion number density turns out to be n a = 8µ eff ξHν(q)/x 0 where the function ν(q) rapidly interpolates between 1 − 1/q for q > 1 and (H/m r ) 1−q for q < 1. It is therefore clear that the spectral index q plays a crucial role in the determination of the axion abundance produced by the string network.
We extract q from simulations using both the physical theory and the fat string trick, with the latter having a cleaner final spectrum with less residual dependence on the initial conditions. We fit q in the range 30 < k/H < m r /4 over which it indeed shows a constant power law behavior. The fitting interval has been chosen somewhat smaller than that over which the network emits axions in order to further reduce possible systematics from finite volume and grid size effects. In Appendix B.4 we show that our results remain consistent as this range is changed, we discuss more properties of the spectra and give details of the simulations used.
The value of q as a function of log(m r /H) is shown in Fig. 2. The data points represent the average of q over many simulations and the error bars measure the associated statistical errors. 9 Although the spectral index is less than unity over the whole simulated range, a nontrivial growth is evident, corresponding to a spectrum that is becoming more IR dominated. The behavior is fit well by a linear function (i.e. q(log) = q 0 + q 1 log) in both the fat and the physical systems (the dark shaded region in Fig. 2). Fits with an extra quadratic term (+q 2 log 2 ) give compatible results (the lighter shaded region in Fig. 2), although with larger uncertainties. This implies that the linear logarithmic growth will continue for, at the very least, a few more e-foldings.
Hence the data in Fig. 2 strongly suggests that the spectrum becomes IR dominated (q > 1) within one or two e-foldings beyond the simulation reach. 10 Note however that the data shown in Fig. 2 represent averages over many simulations: while at early times (log 6) all the simulations that comprise our data sets have q < 1, at late times (log 7.5) a portion already shows an IR dominated instantaneous spectrum with q > 1. This strengthens our confidence that the spectrum indeed turns IR dominated at slightly larger values of log. Further suggestive evidence can be found in Figs. 14 and 15 of Appendix B.4.1, in which the shape of the instantaneous spectrum F at different times is plotted.
This nontrivial log dependence of the emitted axion spectrum correlates with all the other evidence of evolution of the attractor's parameters, in particular with the reduction of UV mode emission. The most conservative extrapolation of the data in Fig. 2 is to values of q larger than unity at late times. Fortunately, as we will explain in the next Section, as long as q > 1 the final axion abundance only has a very weak dependence on its precise value. For this reason we will not attempt to perform a real extrapolation of q from the data in Fig. 2, but we will just assume that at log > 60 its value is definitely larger than unity (say, q > 2).
To summarize, we performed dedicated high-statistics large-grid simulations of the axion string network, providing strong evidence for nontrivial evolution of the network's scaling parameters towards the expected behavior of Nambu-Goto-like strings. In particular, both the string density and the axion spectrum vary in a way that, once extrapolated to the physical parameter region relevant for QCD axions, can make the relic axion component produced by strings orders of magnitude larger than the naive one inferred directly from simulations.
The possibility that topological defects, in particular strings, might provide the dominant contribution of relic axions (much larger than the naive misalignment one) was already argued long ago [32][33][34][35], by assuming that at late time the axion string network's dynamics was well approximated by the Nambu-Goto one, and in particular q > 1. Our results in this Section represent the first clear evidence from full field theory simulations in support of this picture and provide a more detailed characterization of how this limit is approached.

From Strings to Freedom
The scaling regime discussed in the previous Section ends at temperatures of order the QCD scale, when the axion potential becomes relevant and the PQ symmetry is explicitly broken. At this time each string develops N domain walls (where N = v/f a is the QCD anomaly coefficient). For N = 1 (we will discuss the case N > 1 in Section 3.3) there are no conserved quantum numbers left, and the network of strings and walls subsequently decays into axions.
As mentioned in the Introduction, a huge hierarchy of scales forbids a direct numerical study of the system at these times. Given the observed evolution of the properties of the string network (which dramatically changes the dynamics at large scale separations already during the scaling regime) we cannot trust results for the string/wall system dynamics from simulations that are carried out so far away from the physical point. Instead, we focus solely on the contribution of axions produced before the axion potential becomes relevant (i.e. on axions emitted while the system was still in the scaling regime), which requires far fewer theoretical assumptions and extrapolations. To do so, we will study the nonlinear evolution of these axions through the QCD transition in isolation, ignoring the presence of strings and walls and the additional axions they decay into. This allows us to perform direct numerical simulations without the need for any extrapolations. The price to pay is that we miss the component of axions that is produced from the decay of strings and domain walls, which will presumably contribute further to the abundance. In this way we obtain only a lower bound on the final abundance. One may worry that the strings and walls, and the axions produced from them afterwards, could interfere with the evolution of the preexisting axions that we are trying to reconstruct. However, barring an unlikely highly-efficient absorption of background axions by topological defects, their presence is not expected to alter our lower bound considerably, and at worst might weaken it by an order one factor (which, in any case, is not more than other sources of uncertainties that we will discuss at the end of the Section and in Section 4). This fact is further supported by a study in Appendix F.2 where we performed dedicated simulations to analyze the evolution of the axion radiation (as predicted by the scaling regime at log ∼ 60 ÷ 70) when strings and domain walls are included.
Away from topological defects the Hamiltonian density describing the propagation of the axion field is where, as suggested by the dilute instanton gas approximation [36] and supported by recent lattice simulations [37][38][39][40][41] (see also ref. [42] for a recent review), we assume that the axion potential at early times is described by a single cosine potential and the axion mass has a power dependence on the temperature m a ∝ T −α/2 ∝ t α/4 , with α 8 the preferred value. 11 Naively one might think that the axions produced by strings propagate freely like radiation until their momenta (which is typically of order a few H) become of the same order as the axion mass m a , after which they would start propagating as nonrelativistic matter. Throughout this whole process the comoving number density would be conserved. This is true if the axions remain weakly coupled for the whole time. Indeed the axion couplings are suppressed by either H/f a or m a /f a , most of the axions have small momenta of order H ≪ f a and the effective coupling to strings is also suppressed by 1/ log 1. However, as we will see below, the large quantity of axion radiation produced during the scaling regime implies that the average value of the field a 2 1/2 f a , and nonlinear effects have an important effect on the axion number density. This whole process can be studied directly through numerical simulations of the axion field alone. The initial conditions are taken from the axion spectrum emitted by strings during the scaling regime extrapolated to the time t = t , which we define as the moment when H(t ) = m a (t ), since we find that the axion spectrum is still unaffected by the potential at this point (see Appendix D). More axions will be emitted afterwards, however, since their spectrum is unknown, we conservatively do not include them in the initial conditions, and therefore not in our lower bound. Before presenting the results of our simulations we first describe what our expectations are for the effects of nonlinearities. In particular, we derive an analytic formula for the final axion abundance that agrees surprisingly well with the numerical results and correctly reproduces the dependence on the relevant parameters.

Analytic Description
As mentioned in Section 2.2, the energy density of the axions produced by the string network up until t = t is ρ a (t ) ≈ ξ µ eff, H 2 log (from now on the pedex " " on a quantity indicates that it is computed at t = t , log ≡ log(m r /H )), where the last log factor arises from the convolution of axion energy densities emitted over the course of the scaling regime. Using the results of Section 2 on the evolution of the network (in particular the fact that q > 1 long before t * ) the overall energy density ρ a, is distributed with a scale invariant spectrum (up to logarithmic corrections), i.e. ∂ρ a /∂k ∝ 1/k, between the IR cutoff at k ∼ k IR = x 0 H (with x 0 = O(10)) and the redshifted UV scale at k ∼ √ H m r . We refer to Appendix E for the derivation of this result, and to eq. (23) for the explicit form of ∂ρ a /∂k.
The evolution of high frequency modes with k k IR is dominated by the gradient term even long after t = t . Therefore, the nonlinearities arising from the axion potential are negligible for the entire evolution of these modes. As a result, we have to focus only on the IR part of the spectrum, the contribution of which to the energy density is ρ IR ≈ 8ξ µ eff, H 2 ∼ 8πξ log H 2 f 2 a (more precisely we define ρ IR as the integral of the axion spectrum over momenta k < c m m a , with c m = O(1) coefficient, since for higher modes the potential term is subleading). Given the extrapolated values of ξ and log from Section 2, at t the IR axion energy density ρ IR ∼ O(10 4 )H 2 f 2 a is much larger than the contribution from the axion potential This means that at t = t most of the energy density is still contained in the gradient part of the Hamiltonian ( 1 2ȧ 2 + 1 2 (∇a) 2 ). Several implications follow from this fact. First, since the gradient term dominates the Hamiltonian evolution of the field, even the modes with k < m a , which in the linear regime would behave nonrelativistically, will not feel the presence of the potential term and so continue evolving as a free relativistic field after t , until ρ V becomes comparable to ρ IR . Moreover, since the typical gradient of the field is set by H , in order for the gradient term of the Hamiltonian density to account for ρ IR the amplitude of the IR modes needs to be much larger than f a , i.e. a 2 /f 2 a ∼ O(ξ log ). 12 This means that at large ξ log the axion field is mostly a superposition of waves, with wavelengths of order Hubble, that wind and unwind the fundamental axion domain (−πf a , πf a ) several times in a topologically trivial way. Points in space with a/f a ∼ π mod 2π correspond to the core of domain walls with the topology of a sphere. For ξ log 1 there will be multiple domain walls nested inside each others, with a deformed onion-like structure. The presence of these domain walls however does not play any role as long as ρ V < ρ IR since the field continues to evolve freely.
During this period the field keeps redshifting relativistically, the amplitude of the field decreases, ρ IR = ρ IR, (t /t) 2 and the comoving number density of axions remains constant. Meanwhile, as the temperature continues to drop, approaching the QCD transition, the axion mass and ρ V increase rapidly. Eventually, at t = t defined as the time when ρ IR (t ) = c V ρ V (t ) (with c V an order one constant), the presence of the axion potential becomes important and the dynamics turn completely nonlinear. This corresponds to the moment when the domain walls start to be resolved, i.e. when the thickness of each domain wall (∼ 1/m a (t )) has shrunk below the average distance between two walls (∼ k −1 f a / a 2 1/2 ∼ f a /ρ 1/2 IR (t )). Soon after, domain walls, being topologically trivial, annihilate into axions. Except for few loci where oscillons can potentially form (which anyway can only take away a negligible portion of the total energy density, as shown in the next Section), the axion field amplitude continues to drop, rapidly falling below f a . Nonlinearities fade away and conservation of the comoving axion number density is restored.
We can thus assume that during the nonlinear transient (at t ∼ t ) the axion energy density ρ IR ∼ ρ V ∼ m 2 a (t )f 2 a is promptly converted into nonrelativistic axions. The corresponding number density is n str a (t ) = c n ρ IR (t )/m a (t ) c n m a (t )f 2 a , where c n is an order one coefficient taking into account transient effects, extra contributions from higher modes, etc. The value of m a (t ) can be extracted from the definition of t above and for α 1 (i.e. neglecting redshifting effects of ρ IR with respect to the much faster axion mass growth) it is parametrically given by m a (t ) (ξ log ) 1/2 H . We therefore expect that, up to order one factors, the axion number density after the nonlinear regime is n str a (t ) (ξ log ) 1/2 H f 2 a , i.e. it is enhanced by a factor (ξ log ) 1/2 with respect to the misalignment contribution. Note that the enhancement, while substantial, is parametrically smaller than the naive one obtained by assuming that the axion field remains linear throughout the QCD transition, which would be ξ log .
The main effects of the nonlinearities can be simply summarized as follows: the large energy density stored in the axion gradient term delays the moment when the axion mass and potential become relevant. In the meantime the axion mass is growing fast, so that, by the time the potential becomes relevant and the axions nonrelativistic, more energy is required to produce each axion and the comoving number density is suppressed.
The estimate above can be improved by keeping all the order one factors, taking into account the actual shape of the spectrum and the effects of redshifting from t to t . The full computation 12 See eq. (25) in Appendix E for an explicit derivation based on the spectrum. is discussed in Appendix E.1 and gives where n mis,θ 0 =1 a (t ) is the axion number density from misalignment with θ 0 = 1 redshifted to t , the prefactor A ξ log ,x 0 is a function of all the parameters (including the order one coefficients c m , c V , c n ) but with only a mild logarithmic dependence on ξ log , and x 0 (the full form is given in eq. (36)). The dependence on q is further suppressed by 1/ log , as shown in Appendix E.1.
The result in eq. (5) assumes that ξ log 1 and involve several approximations, parametrized by some unknown order one coefficients. These crudely describe the number of IR modes involved in the nonlinear dynamics (c m ), the relative importance of the potential versus the gradient energy in IR modes when nonlinearities become relevant (c V ) and the conversion factor of energy density into number density during the true nonlinear transient (c n ). While these numbers can only be fixed through numerical simulations, the full dependence on ξ log as well as the subleading ones on α and x 0 are genuine prediction of our analysis. As we will see next, they are nicely reproduced by the numerical simulations.

Comparison with Simulations
The dynamics discussed above can be checked by numerically integrating the axion equations of motion from the Hamiltonian density in eq. (4). We start the simulations at t = t with initial conditions set by the axion field radiation that would be produced during the scaling regime for different values of ξ log and x 0 . The form of the spectrum is characterized by the position of the IR cutoff (x 0 ) and the spectral index of the instantaneous spectrum (q), while the overall size is controlled by the parameter ξ log . We carry out simulations with different values of α, which fixes the temperature dependence of the axion mass. More details are given in Appendix E.2.
For sufficiently large ξ log the numerical simulations show that the system indeed continues to evolve as in the absence of a potential after t = t , redshifting as radiation and with a conserved comoving number density. More details and plots are given in Appendix E.3. The larger ξ log is, the longer the period of relativistic redshift lasts. This regime ends, as expected, with a nonlinear transient, after which the mean square field amplitude rapidly drops below f a (see Fig. 23).
At this point the field settles down around the minimum of its potential at a = 0, oscillating with an amplitude that is much smaller than f a almost everywhere. Consequently, the system becomes linear again except in a few localized regions of size m −1 a (t) where the field continues oscillating with an amplitude of the order πf a . These objects, remnants of the large initial field amplitude (with a > f a at t = t ), are known as oscillons or axitons [44,45]. Oscillons are heavy and slowly decay radiating their energy density into axion waves with momentum of order m a . Their lifetime is long enough that they persist until the end of our simulations. However, only a very small portion of the energy density remains trapped in oscillons, so that their presence is irrelevant for the computation of the final axion abundance. More details about the oscillons can be found in Appendix E.3. Everywhere else the axion field is in the linear regime by the end of the simulations. We can therefore calculate the total axion spectrum ∂ρ a /∂k and number density n str a = dk(∂ρ a /∂k)/ω k (ω k = k 2 + m 2 a ). We do so screening away the regions occupied by oscillons, and we use the difference with the unscreened results to estimate the uncertainty introduced by the presence of these objects. As anticipated the difference is small, which confirms that only a negligible portion of the energy density is trapped in oscillons. Moreover, as expected, after the screening the conservation of the comoving number density further improves. Additional discussion and plots are given in Appendix E.2. Thanks to the rapid growth of the axion mass, the nonlinear regime is reached not long after t and the system soon becomes linear again, after a short transient, as the field relaxes below f a . For this reason, in the range of ξ log and f a under consideration (ξ log 10 5 , f a 10 11 GeV), the system reenters the linear regime (and our simulations end) at temperatures that have dropped by, at most, a factor of four from that at t . This is still above the QCD transition, in a regime where the axion potential used in eq. (4) should hold.
The agreement between the numerical simulations and our analytic description is not only qualitative but also quantitative. We compare the two using the ratio Q = n str a /n mis,θ 0 =1 a of eq. (5), between the number density of axions from strings and the reference one from misalignment (with initial misalignment angle θ 0 = 1). Since both n str a and n mis,θ 0 =1 a are conserved per comoving volume at late times, Q asymptotes to a constant value. The results for Q for α = 8 and x 0 = 5, 10, 30 are plotted in Fig. 3 as a function of ξ log , and the comparisons for the other values of α are reported in Appendix E.2. The three parameters of the analytical formula c m , c V , c n have been fixed with a global fit of Q including all the simulations with different values of α, x 0 and ξ log . The agreement between the theoretical estimate and the simulation data is remarkable given that: We also indicate the comoving momenta corresponding to the physical momenta that are equal to the axion mass at t = t (i.e. k com = m a (t )(H /H ) 1/2 , filled dots), which is parametrically the time when the nonlinear transient occurs.
1) all data is fit with just three universal parameters which indeed turn out to be of order one, and 2) the dependence on ξ log , which is a prediction (not the result of a fit), agrees very well over multiple orders of magnitude. The only slight deviation is at low values of ξ log where the approximations used in the analytic formula are not in fact valid. More details about the dependence on the input spectrum, the values of the fitted input parameters and the dependence on the other parameters can be found in Appendix E.2. Here we simply note that, as anticipated, the dependence on the spectral index q of the spectrum from the scaling regime is negligible as long as q is away from unity. Although the numerical simulations are capable of covering the parameter space relevant for the QCD axion (discussed in Section 4), our analytic formula, in addition to providing a better understanding of the physics behind the nonlinear effects, would allow us to interpolate and extrapolate the simulation results to other values of the parameters if needed.
We will discuss the phenomenological implications of our results in Section 4; first we analyze the effects of nonlinearities on the shape of the final axion spectrum in more detail. As shown in Fig. 25 in Appendix E.3, the spectrum continues to redshift almost unaltered after t until it reaches the nonlinear regime at around t = t . At this point the energy contained in modes k m a (t ), is converted into massive nonrelativistic axions. For this to happen axions with k m a (t ) need to combine with each other to generate on-shell axions with mass m a (t ), and the comoving number density of this component cannot be conserved. In other words, nonlinearities remove the IR part of the spectrum via 3-to-1, 5-to-3, etc. processes. The smaller k-modes are those with the larger occupation number and they therefore suffer stronger nonlinear effects. The resulting spectrum after the end of the nonlinear transient will therefore be peaked at physical momenta that were of order m a (t ) at t = t , which is significantly higher than the would-be peak at x 0 H (at t = t ) had the nonlinearity been absent. In particular, the value at the peak grows with ξ log . This is shown in Fig. 4 where we plot the spectrum as a function of the comoving momentum k com ≡ k(H /H) 1/2 for the three values of ξ log = 10 2 , 10 3 , 10 4 , at the initial time t = t and at the final simulation time.
The deformation of the spectrum above could have important implications for the properties of the small scale structures produced by the axion inhomogeneities known as mini-clusters [46] (see also [24,25,47,48] for recent studies). Figure 4 also shows the role of oscillons. These only affect the spectrum at momenta of order m a , indicated by empty dots (see Appendix E.3 for details). Since the largest contribution to the number density comes from the peak of the spectrum, once m a (t) is sufficiently above this (as is the case at late enough times) the screening of oscillons does not significantly affect the measured axion number density. This matches the results for the number density evaluated directly, described above.
We finish this Section by briefly discussing the possible effects of the presence of strings and domain walls during the QCD transition, which have been omitted so far. We first note that at t = t the energy density in the string network is comparable to that we considered from the IR part of the axion radiation (ρ IR ), and it is mostly localized along the strings themselves, so the dynamics of the field away from the strings should be largely unaffected by their presence. After t domain walls start to form but their energy density is bounded by the axion potential, and becomes relevant only much later, when the axion field has relaxed to values a ∼ πf a everywhere. Hence away from strings we do not expect the dynamics of the axion field to be significantly different from those we computed, at least until t ∼ t . At this point the nonlinear transient starts. The difference with respect to our simplified case is that, as well as our topologically trivial domain walls, extra walls surrounded by strings are also present. If the extra string-wall network decays during the transient, then as we saw before the total energy density (that in strings walls and radiation) is expected to convert into axions with a conserved comoving number density of order (ρ tot (t )/m a (t )). If for some reason 13 the string-wall network were to survive for longer, away from them the field would still evolve as calculated above. Therefore, we would expect that the results given above should represent, up to O(1) factors, a lower bound on the axion abundance regardless.
It would be quite surprising if the extra string-wall system were able to wipe away the bulk axions with a high enough efficiency to suppress their big contribution to the final abundance significantly. To further exclude this possibility we performed dedicated simulations where, as well as the axion radiation predicted by the scaling regime, we also included the strings (and the domain walls that form from them) during the mass turn on. In these simulations, the background axion radiation is as it would be with the physical parameters (i.e. with the spectrum and energy density expected at log ∼ 60 ÷ 70). However, the string-domain wall network is evolved with the currently allowed log(m r /H) 7, so ξ log for the string system is much smaller than the physically relevant value. As expected, the presence and decay of strings and domain walls does not significantly alter the evolution of the preexisting background radiation, and thus does not decrease the final abundance. Since in such simulations ξ log for the string network is small, and the emission from the decay of strings is UV dominated, the inclusion of strings also does not noticeably increase the final abundance. We refer to Appendix F for more details and the explicit results of these simulations.
From this study we learn two important lessons. Calculations of the axion abundance from brute force simulations of the whole evolution of the string-domain wall system can easily miss the dominant source of axion emission, underestimating the final relic abundance by more than one order of magnitude. Moreover, the explicit inclusion of strings in the late evolution of the field does not play a role unless their contribution starts becoming comparable to that from radiation during the scaling regime, at which point a tuned cancellation among the two sources would be surprising.

The case N > 1
We now discuss the generalization of our results to the case N = v/f a > 1. First notice that in the equations of motion from the Lagrangian in eq. (1) the scale v can be removed by rescaling the complex scalar field φ. This means that the string dynamics during the scaling regime do not depend on v. The way v enters observables is just fixed by dimensional analysis, and in particular all energy densities, number densities and the string tension are proportional to v 2 . Therefore, the axion spectrum produced during the scaling regime in the general case can be recovered by simply multiplying the results of Section 2 by N 2 .
On the other hand, the axion potential produced by QCD in eq. (4) involves the scale f a . In all our computations in Section 3, the scale v only enters through the axion spectrum via the scaling solution used as an input, where it appears in combination with ξ log . All the results in Section 3 can therefore be generalized by simply substituting ξ log with N 2 ξ log (e.g. in eq. (5) and in Figs. 3 and 4).
The effect of v > f a is therefore to increase the energy density of the axions produced by strings (as a result of the enhanced string tension), increasing the field amplitude and therefore the effects of nonlinearities. Roughly, the final number density of axions will be enhanced by an O(N ) factor, and the peak in the final spectrum will be UV shifted by a similar amount.

Results and Phenomenological Implications
We can now extract some phenomenological implications from the results of the previous Sections. In particular, given the extrapolated values of the axion spectrum from the string scaling regime, Fig. 3 and Section 3.3 provide a lower bound on the axion number density (in terms of the easily computed misalignment result). This can be translated into corresponding bounds on the axion mass and its decay constant requiring that such an abundance does not exceed the current observed dark matter value. As reference values we choose ξ = 15, x 0 = 10, q > 2, α = 8, log = 64, which for N = 1 (as in the minimal KSVZ model [49,50]) imply m a 0.5 meV , f a 10 10 GeV , while for N = 6 (as occurs in e.g. the DSFZ model [51,52]) they imply m a 3.5 meV , f a 2 · 10 9 GeV .
For comparison, the naive axion number density from misalignment in the post-inflationary scenario is obtained by averaging the misalignment relic abundance with a flat θ 0 distribution in the interval [0, 2π]. This gives n mis,avg a 5n mis,θ 0 =1 a , which corresponds to m a 0.028 meV and f a 2 · 10 11 GeV, more than an order of magnitude weaker than our bound.
We do not think that it would be fair to associate an error to the figures in eqs. (6) and (7): shifts of O(1) could be expected, but we would be surprised if these bounds relax by significantly more than a factor of two. To provide a better feeling for the main sources of uncertainty, and the choices of parameters used, we will now go through all the assumptions underlying the numbers above: ξ : We fixed the value of ξ = 15 from the best fit of the scaling solution described in Section 2.
As discussed at length this number assumes that the linear-log behavior observed in simulations extends beyond the simulation range by another order of magnitude. While such an assumption can be questioned, other independent studies support a similar enhanced value. These include refs. [22,26], which partially reproduce the possible effects of an enhanced string tension and find ξ 5; Nambu-Goto simulations, which seem to prefer values between 10 and 20 [18,19,27]; and recent local string simulations, which give ξ = 15(1) [28]. 14 Since the final abundance approximately scales as ξ 1/2 , a factor of 2 shift in ξ would affect the final bound by less than 40%, within our target precision. Substantially larger deviations from our central value seem unlikely. q: The main assumption behind the result above is associated to the spectral index q being larger than unity. Although present simulations cannot provide a proof, our analysis in Section 2 shows that q > 1 is by far the most conservative extrapolation of the results from simulations. This extrapolation is also supported by theoretical arguments about the expectation that the string network approaches the Nambu-Goto dynamics at large values of log . With q > 1 the instantaneous axion spectrum emitted by strings is IR dominated. The corresponding integrated spectrum (which determines the final abundance) is therefore fixed and only very weakly dependent on the actual value of q. As a result, the actual extrapolated value of q does not lead to large uncertainties in our estimate (see Appendix E.3).
log : We set the reference value of log = 64, corresponding to fixing f a 10 10 GeV, m r f a and the value of α = 8 (discussed below). Much smaller values for log are in principle allowed (m r keV from astrophysical and fifth force experimental constraints) although these are much less plausible given that the smallness of m r would not come for free.
x 0 : We set the position of the IR cutoff of the spectrum x 0 = 10 from the results of the simulations at log 8 (see Fig. 14). Within the available range of simulations x 0 is consistent with being constant, although we cannot exclude a slow evolution which would change its value at large log . One might indeed expect that x 0 could increase with ξ 1/2 , as the average interseparation between strings is reduced. 15 The study of the evolution of x 0 from simulations is more challenging than that of q since it is much more sensitive to finite volume effects and it requires a better understanding and modelling of the shape of the IR peak. Fortunately, as discussed in the previous Section and explicitly shown in Appendix E.1, the final abundance is only logarithmically sensitive to x 0 , so even a substantial change at large values of log has a limited impact on the final result. This can also be seen in Fig. 3.
α: We set the index that controls the temperature dependence of the axion potential to the value α = 8, which corresponds to the prediction of the dilute instanton gas approximation at weak coupling. Although this computation is probably out of its regime of validity at the temperatures we are interested in, the same value of α seems to be supported by the most recent lattice QCD simulations. Waiting for an independent check we adopt this value as the reference one and provide the results for generic α in Appendix E. Similarly the results in refs. [37,41] suggest that a simple cosine is a good approximation to the axion potential for the temperatures relevant to the nonlinear regime. 16 -Extra strings-domain walls contribution: The last source of systematic error comes from neglecting the extra contributions from strings and domain walls present after t . As discussed at length in the previous Section, we expect these to add further to the axion abundance, hence our lower bound. If the extra contribution is subdominant then our bounds would turn into central values for the abundance. If the extra contribution dominates, they will become strict inequalities. We cannot exclude a partial destructive interference between these extra contributions and the axion spectrum produced at earlier times. However it would be highly unlikely that it could weaken the bounds in eqs. (6) and (7) beyond an O(1) factor, i.e. by more than the size of the other uncertainties.
When combined with astrophysical constraints, the bounds in eqs. (6) and (7) restrict the allowed parameter space for the QCD axion in the post-inflationary scenario quite substantially. In particular, they motivate efforts to further explore a region of parameter space that could in principle be probed by astrophysics, as well as axion dark matter [53][54][55][56] and non dark matter [57][58][59] experiments. In Fig. 5 we show our bound for the QCD axion mass in the post-inflationary scenario, together with constraints on the axion-photon coupling g aγγ 17 from currently running experiments and the parameter space that could be probed by proposed experiments. 18 Interestingly, the allowed window is almost complementary to that of the pre-inflationary scenario. The upper bound on the mass in this case is f a 3.7 · 10 9 GeV and comes from requiring 15 We thank Javier Redondo for a discussion on this point. 16 We also note that the uncertainty introduced by the number of degrees of freedom in thermal equilibrium at Figure 5: The axion parameter space in terms of its mass and coupling to photons g aγγ . Solid green lines indicate the allowed parameter space in the post-inflationary scenario for the minimal KSVZ model and the DFSZ model, given our constraints from dark matter overproduction (eqs. (6) and (7)), while the dashed green lines indicate our estimate of the uncertainty on these results. The green shaded band indicates the post-inflationary parameter space allowed by the bound for more general axion models. Vertical green lines show our lower bounds on the axion mass at m a = 0.5 meV and m a = 3.5 meV for N = 1 and N = 6 respectively. We also indicate, in red, the allowed axion masses in the pre-inflationary scenario (the corresponding g aγγ lie in the partially transparent grey band), the upper limit of which m a 1.5 · 10 −3 eV is set by isocurvature fluctuations. Existing experimental bounds and observational constraints (solid lines) on g aγγ as a function of the axion mass and the projected sensitivity of proposed experiments (dotted) are also shown. The limit on DFSZ models from white dwarfs and red giants ("WD/RG") is indicated for tan β = 1 [74] , while the supernova-1987a limit on such models ("SN1987a") spans the blurred region as tan β varies (the corresponding constraint on KSVZ models is m a < 15 meV) [69]. The post-inflationary region that we identify could also be probed by future experiments sensitive to the axion's couplings to matter [53,58]. In combination with the bound from supernova, the region of viable QCD axion masses in the post-inflationary scenario is restricted to m a ≈ 0.5 ÷ 20 meV.
that the Hubble parameter during inflation is small enough to avoid observational constraints on isocurvature from Planck [82], but at the same time above m a so as not to deplete the misalignment abundance during inflation. In fact, in the overlapping region both the pre-inflationary and postinflationary scenarios predict nontrivial small scale structures from the axion self-interactions [83], although the details are expected to differ as a consequence of the different origins of field inhomogeneities.

A The String Network on the Lattice
In this Appendix we summarize the methodology behind our numerical simulations of the scaling regime. In these we evolve the equations of motion of the Lagrangian in eq. (1), where∇ is the gradient with respect to the comoving coordinates, on a discrete lattice with a finite time-step 19 (we fix v = f a as in Section 2). We assume radiation domination in a spatially flat Friedmann-Robertson-Walker background, so the scale factor grows as R(t)/R(t 0 ) = (t/t 0 ) 1/2 , t = t and t = t , and by the changes between these times, is relatively small, certainty within our target precision. 17 Defined as L ⊃ − 1 4 g aγγ aFF in the low energy theory, where F is the electromagnetic field strength. 18 The existing experimental and observational bounds shown are from ADMX [60,61], earlier cavity experiments "UF/RBF" [62,63], HAYSTACK [64], CAST [65], observations of horizontal branch starts "HB" [66], supernova 1987a "SN1987a" [67][68][69] (see however [70]) and red giants and white dwarf stars "RG/WD" [71][72][73][74]. The constraints on DSFZ axions from supernova 1987a, red giants and white dwarfs are model dependent via the mixing angle β. For those from white dwarfs and red giants we plot the limit from [74] for tan β = 1. The limit from supernova 1987a is m a < 0.02 eV for tan β = 1 and this barely weakens for smaller tan β but it strengthens slightly (up to the edge of the blurred region) for large tan β [69]. The proposed experiments shown are ABRACADABRA [75], superconducting radio frequency cavities "SRF" [76] (see also [77]), CULTASK [78], MADMAX: [79,80], tunable plasma haloscopes "Plasma" [55], TOORAD [54], phase measurements in cavities "phase" [59], absorption by gapped polaritons "polaritons" [56] and IAXO [57]. 19 Many of the details of our implementation follow those described in Appendix A of [7]. For example, it is most convenient to work in terms of the rescaled field ψ = R(t)φ/f a , so that the Hubble term in eq. (8) is canceled. Rather than reviewing all such technicalities, here we focus on the key features and the differences in our present work.
We carry out simulations of the physical string system, for which m r is constant, and also the so-called fat string system in which m r (t) = m r (t 0 ) (t 0 /t) 1/2 . The core-size of strings is characterized by the length scale associated to the region where |φ| < f a / √ 2 and this is set by m −1 r . Consequently, for physical strings the number of lattice points per string core decreases through a simulation, while for the fat string system it remains constant.
As discussed in Section 2, for a given grid size the maximum log(m r /H) that a simulation can reach is limited by the simultaneous requirements that systematic errors from the finite lattice resolution and from the finite box size do not become too large. The former constrains the maximum value of m r ∆, while the latter imposes a lower bound on HL, where L is the physical box length, defined in Section 2. The corresponding maximum log(m r /H) is the same for simulations of fat and physical strings, however the fat string system evolves for a longer cosmic time before this is reached. The maximum gridsize is limited by the available computational resources, and we carry out simulations with up to 4500 3 lattice points. 20 The relatively large values of log(m r /H) accessible with such grids are vital in identifying the evolution of q described in the main text.
The numerical values of m r ∆ and HL that can be used without introducing significant errors must be determined by direct testing in simulations. For log 6.5, m r ∆ = 1 is sufficient for most observables of interest [7], however the bigger values of log in our present simulations necessitates that these are re-analyzed. We carry out a study of the finite lattice spacing effects from m r ∆ in Appendix B, where we show that some observables are indeed increasingly sensitive as log increases. To maximise the accessible value of log, we ran simulations until HL = 1.5. For this value ξ still coincides with the infinite volume limit (see Appendix B.1). The IR part of the spectrum starts being slightly distorted, but not in the range of momenta used for the extraction of q, which still coincides with the infinite volume limit (see Appendix B.4.2). With these choices of m r ∆ and HL, our simulations reach log 8.

A.1 Selecting the Initial Conditions
For simulations to show the properties and log dependence of the attractor solution as clearly and accurately as possible, the initial conditions need to be fixed as close as possible to the scaling solution. If this is not done the network will go through a transient period as it approaches the attractor, decreasing the range of log over which its properties can be reliably studied. Indeed, the dynamics during the transient will differ from those in the scaling regime, and ξ and q might not show their true asymptotic evolution.
One requirement to be on scaling is connected to the initial density of strings. Simulations with too small a density will fail to reproduce the right properties associated to string interactions responsible for maintaining the attractor regime. Meanwhile, too large densities will lead to an enhanced string interaction rate and an overproduction of radiation with respect to the scaling 20 To achieve this we use MPI parallelization across multiple (up to 48) cluster nodes. regime. A possible criterion to identify the optimal initial conditions is to choose those with the highest density of strings that do not show a clear initial drop of ξ before the observed universal asymptotic growth.
Another source of systematic noise is associated to initial excitations of the strings core. For example, such excitations will be triggered if the initial configuration contains strings with coresizes that are significantly different to those on the attractor, which are parametrically set by m −1 r . As the network evolves the strings cores relax to the properties they have on the attractor regime emitting UV radiation that pollutes the axion spectrum (mostly around the frequency m r /2, as a result of the parametric resonance with the radial modes -see Appendix B.4). Although at late times (log 60 ÷ 70) such radiation is completely negligible because of the huge redshift, the effect can be sizable in the limited extent of simulations.
The initial conditions are more important when studying the physical string system than the fat string one for two reasons. First, thanks to the longer cosmological time range, the fat string system reaches the attractor in a fraction of the total time (i.e. at smaller values of the log) even with untuned initial conditions, and the radiation left over from early times is diluted fairly efficiently by redshifting. In contrast, for physical strings the transient can easily last for the entire span of the simulation. Second, in the fat string system m r /2 corresponds to a fixed comoving momentum. Therefore, the radial and axion modes emitted by the string cores only affect the UV part of the spectrum, outside the region of interest. Meanwhile, for physical strings these modes are redshifted towards smaller frequencies, polluting part of the spectrum used to extract the instantaneous emission F by creating large oscillations (we will see this in more detail in Appendix B.4).
We use initial conditions that contain a fixed (adjustable) density of strings. For the fat string system, these are obtained by evolving eq. (8) starting from a random field configuration until the required total string length inside the box is reached. The field at this time is then used as the initial condition for the main simulation. The strings produced by such a procedure do not generally have the correct core size, since the Hubble parameter at the end of the initial simulation does not match that at the start of the main simulation. However, as mentioned above, the subsequent readjustment of the string cores only affect the UV part of the spectrum and has no consequences for the study of the attractor properties.
For simulations of the physical string system we modify this procedure slightly to overcome the issue with the spectrum discussed above. We generate initial conditions for these by evolving eq. (8) with R ∝ t and m r ∝ R −1 , until the desired total string length is reached. This choice has the advantage that both the Hubble parameter H =Ṙ/R and the string core size m −1 r are constant in comoving coordinates, i.e. H −1 /R and m −1 r /R do not change. We chose m −1 r /R equal to m −1 r /R(t 0 ), i.e. the core-size at the beginning of the actual simulation. Consequently, the strings in the initial conditions have the right core-size, no matter what time the preliminary evolution ends. Moreover, the simulations to generate the initial conditions can run for an arbitrarily long time, since the comoving Hubble radius and the comoving core-size do not change. For small H −1 /R this evolution corresponds to a system with large Hubble friction, which acts as a relaxation period smoothing out fluctuations of the initially random field and diluting preexisting radiation. Meanwhile strings form and their core-sizes relax. We choose H −1 /R = H −1 (t 0 )/R(t 0 ).
With this value the total string length in the box decreases fairly slowly, so the relaxation period lasts a significant amount of time.
For the fat string system we start the main simulations at log(m r (t 0 )/H(t 0 )) = 0. For the physical system we found that cleaner initial conditions were obtained by choosing log(m r /H(t 0 )) = 2. When studying the evolution of ξ we varied the initial string density. Meanwhile, when analyzing the spectrum and energies we fixed the initial string density close to the scaling solution. With our method of generating initial conditions, this happens for ξ(t 0 ) = 10 −2 and ξ(t 0 ) = 0.2 for fat and physical strings respectively.

A.2 String Length and Boost Factors
To identify strings and calculate their length, we adopt the algorithm proposed in Appendix A.2 of [21]. This involves counting the plaquettes that are pierced by a string, and converting the result to a length using a statistical correction factor. In doing so it is assumed that the strings are equally distributed in all directions. 21 We calculate the boost factor γ in two ways, the first as in Appendix A.2 of ref. [21] and the second as in ref. [30]. Briefly, the first method estimates the string velocity from the relativistic contraction of the string core, extracted from the derivative of the field on the gridpoints near the center of the string. The second instead measures the speed at which the points x such that φ(t, x) = 0 change in time. Both methods give the (local) γ-factor at each gridpoint where a string is identified. The frequency distribution function dξγ dγ of γ-factors throughout the string network, defined in Appendix B.2, can be calculated easily. The average γ-factor of the network is defined as the mean over all the gridpoints where a string is identified. We checked that both methods give approximately consistent results, however the method of ref. [30] leads to superluminal γ-factors on some grid points, which have to be discarded in the counting. 22 Therefore we base our analysis on results using the method of ref. [21].

B Properties of the Scaling Solution and Log Violations
In this Appendix we discuss the properties of the scaling solution in more detail. We emphasize how the logarithmic violations of the naive scaling law affect different observables, such as the number of strings per Hubble volume, the relativistic boost factor of the network, the energy emitted in heavy radial modes and the axion spectrum. The dependence on log(m r /H) of these properties (along with the supporting evidence from the dynamics of single loops studied in Appendix C) points to a consistent picture where the heavy degrees of freedom slowly decouple from the string network in the limit log(m r /H) → ∞.

B.1 The Scaling Parameter
The evolution of scaling parameter ξ provides one of the clearest pieces of evidence of the attractor solution, and of the logarithmic violations of its scaling properties. Both of these features are already manifest for the physical string system in Fig. 1 and are even more evident in the fat string one.
In Fig. 6 we show ξ for the fat string system as a function of log(m r /H) with different initial string densities. For each initial condition we ran multiple simulations to reduce statistical errors. Thanks to larger cosmic time available to reach the same value of log, the results for ξ converge to the attractor at smaller values of the log. The growth of ξ on the attractor solution appears linear over a substantial range of log.
As discussed in the main text, the time-dependence of the scaling parameter is fit well by a universal linear function plus corrections proportional to powers of 1/ log. The latter encode the residual dependence on the initial conditions, and vanish in the large log limit. Including the first two such corrections, we perform a global fit of all of the ξ data (separately for physical and fat strings) with the function in eq. (3), where c 0 and c 1 are universal while c −1,−2 are let differ for each set of initial conditions. We include only points with log > 4.5 and log > 4 for fat and physical strings respectively and weight with the statistical errors. 23 For the latter we rescaled the χ 2 function so as to have one independent contribution for every Hubble e-folding; this is in order to avoid bias from data set with a finer time sampling and decorrelate data from consecutive time shots. The result of the fit is fairly good, with a reduced χ 2 phys ≈ 1.1 and χ 2 fat ≈ 1.5. This indicates that eq. (3) is sufficient to capture the evolution of ξ for the entire broad range of initial conditions considered. By the end of the simulations the fitted values of the parameters are such that the 1/ log corrections are already subleading. 24 This is particularly true for the fat string network simulations that converge to the attractor solution at smaller values of log. As it is clearly noticeable in Fig. 6, for the initial conditions closest to the attractor solution (these correspond to the black data set, which has the smallest values of c −1,−2 ), 1/ log corrections are already negligible for the entire range of log plotted. Indeed, any fitting functions with sizable nonlinearities at late times is highly disfavored.
The coefficient of the linear log term c 1 is particularly important for the extrapolation to the physically relevent regime. The results for this in the fat and physical string systems are 25 Finally, we note that it has previously been shown that the percentage of the total string length in loops with size smaller than Hubble stays constant in time [7]. This means that the logarithmic violation in ξ are reflected in a corresponding increase in the string length contained in small loops as well as long strings. This provides another strong piece of evidence that logarithmic violations are a genuine feature of the scaling solution.
All the simulations in Figures 1 and 6 have 0.5 ≤ m r ∆ ≤ 1, and the curves that reach later times are the average of multiple sets of simulations with different values of m r ∆. For any given initial condition the results from simulations with different m r ∆ give compatible results. This suggests that ξ is already close to the continuum limit for m r ∆ = 1 at least up to log ∼ 8. Similarly, since the higher resolution simulations are stopped at HL = 1.5, when the simulations with poorer resolution have HL 1, the agreement indicates that ξ is close to the infinite volume limit for HL = 1.5.

B.2 String Velocities
Another important quantity characterizing the dynamics of the network is the boost factors of the strings. Indeed, if strings are relativistic with an average boost γ , their energy per unit length is increased by a factor γ . 26 The theoretical expectation for the string tension µ πf 2 a log(m r /H), 23 Since the most overdense set in the physical case reaches scaling only late, we include data from this set at log > 5.5 in the global fit. 24 If data at smaller logs is included higher corrections are needed to get a good fit. Meanwhile, selecting only data at larger logs still leads to a good fit but with greater uncertainties on the coefficients. 25 These are compatible with those reported in our previous analysis [7], in which we studied the network up to log = 6.7. 26 We always refer to the transverse boost, which does not lead to relativistic contraction of the string length. Consequently our definition of γ , which does not have any extra weighting, gives the appropriate modification to the string tension. which holds for strings at rest, is correspondingly modified to µ = γ πf 2 a log(m r /H). Therefore, the value and possible log dependence of γ must be understood so that the energy densities during the scaling regime can be determined correctly.
In Figure 7 we show the time evolution of the network's average γ-factor for fat and physical strings and for different lattice spacings (computed as described in Appendix A.2). Evidently the boost factor is lattice spacing dependent, with γ smaller for coarser lattices. This is not surprising given that the boost factor is measured by the size of the string cores, which might not be resolved when they are relativistically contracted. The continuum extrapolation indicates that, at least for fat strings, γ seems to asymptotically approach a constant mildly relativistic value γ = 1.3 ÷ 1.4. 27 More detailed information about string velocities can be inferred from the γ-factor distribution function. We define ξ γ to be the portion of ξ with boost factor smaller than γ. Then 1 ξ dξγ dγ describes the distribution function of γ-factors in the network; γ = dγ γ 1 ξ dξγ dγ being its first moment. In Figure 8 we plot the velocity distribution function at different times for fat strings (already extrapolated to the continuum limit) and physical strings (for two different lattice spacings). The distribution is strongly peaked at nonrelativistic boost factors at all times with a sharp fall off ∝ γ −6 . This makes the lowest moments of the distribution dominated by low boosts factors and explains why γ appears time independent. On the other hand, the large boost tail of the distribution keeps extending to increasingly large values at later times. Meanwhile, systematic errors from the finite lattice spacing affect the distribution at large γ more at late times. dγ for fat strings (already extrapolated to the continuum limit) at different times. Right: The same function for physical strings for two different lattice spacings (m r ∆ f = 1 has the darker color, m r ∆ f = 1 the lighter). At all times, the distribution is peaked at γ = 1 with a sharp fall off ∝ γ −6 above this, so that γ is nonrelativistic. As the log increases, the high-γ tail grows, suggesting that sub-horizon loops become increasingly relativistic.
This evolution with the log can be understood as follows: γ is dominated by long strings, which are mostly nonrelativistic (due to causality and Hubble friction) and make up the majority of the length at all times (see Section 3.3 of [7]). Therefore the network remains on average only mildly relativistic. Loops with size much smaller than Hubble provide a small, constant, proportion of ξ, but they get more relativistic as the log increases. These contribute to the high-γ tail of the distribution, so this extends to larger γ when the log is bigger. Similarly, kinks and cusps from string recombination also contribute. Such a logarithmic scaling violation is in agreement with the results in Appendix C where we show that single loops with larger initial logs get boosted more as they shrink.
We assume that the behavior identified above persists at large values of the log, and in particular that the mean γ-factor remains approximately constant. In this case the theoretical prediction for the string tension µ πf 2 a log(m r /H) can be extrapolated to the physical scale separation (up to an overall γ ≈ 1.3 ÷ 1.4 constant factor, which we fit at small logs in the next Subsection).

B.3 Effective String Tension and Radial Mode Decoupling
We now show that the string tension calculated in the simulation is in agreement with the theoretical expectation. Moreover, we show that the percentage of energy in radial modes, although non-negligible for small logs, decreases at late times, signaling the decoupling of heavy modes in the limit log → ∞.
The total energy density of the complex scalar field ρ tot ≡ T 00 , where T 00 is the Hamiltonian However, the increasing spread between the m r ∆ = 1 and m r ∆ = 1.5 simulations indicates that this is a lattice effect, and the extrapolated tail at log = 7 would be above that at log = 5.5. density from the Lagrangian in eq. (1), can be split into components as where ρ s is the energy density in strings, ρ a that in axion radiation and ρ r that in radial modes. ρ a is extracted from the kinetic energy density of the axion field 2 1 2ȧ 2 away from string cores, and ρ r from the energy density of the radial field 1 2ṙ 2 + 1 2 |∇r| 2 + V (r) , again away from string cores (see [7] for more details). We then obtain ρ s from the difference ρ s = ρ tot − ρ a − ρ r . The string energy ρ s calculated in this way is expected to match the one predicted from the theoretical expectation for the string tension and the measured values of ξ(t) up to order one coefficients.
We compute the effective tension µ eff = ρ s (t)t 2 /ξ(t) from the definition in eq. (2), using ξ(t) and ρ s (t) from simulations. This is then compared to the theoretically expected form accounting for a non-zero γ-factor, and for the dependence of the average inter-string distance on the string density (via the factor 1/ √ ξ in the log). The coefficient η encodes the string shape, and we chose η = 1/ √ 4π as a reference (we pick this somewhat arbitrarily based on the average distance between strings if they were all parallel, but any roughly similar value would also be reasonable).
In Figure 9, we plot the ratio µ eff /µ th as a function of time for the fat string and the physical systems, where γ in µ th is calculated in simulations as in the previous Section. Different colors represent different lattice spacings, and the blue shaded region shows the effect of varying η in the interval  Figure 10: The fraction of total energy density in simulations that is in radial modes for physical strings (right) and fat strings (left). The results are shown for different lattice spacings, and the continuum extrapolation is also plotted for the fat string system. As log grows this fraction decreases for fat strings (after the continuum extrapolation), which suggests that radial modes decouple and are increasing irrelevant for the string dynamics. The results for physical strings show a similar trend, although we do not perform the continuum extrapolation.
is not unexpected, given that strictly speaking eq. (11) only applies for straight strings and we do not have a reliable way to compute η analytically. Instead, its value is determined by the loop distribution and the shape of the strings. Although γ and ρ s are rather sensitive to lattice spacing effects, µ eff /µ th involves the ratio of the two and seems to have smaller systematic error. The approximately constant behaviour in Figure 9 gives us confidence that eq. (11) can be used to calculate the string tension at large logs.
In Figure 10 we plot the proportion of the total energy that is in radial modes as a function of time for different lattice spacings, i.e. ρ r /ρ tot . We also show the continuum extrapolation in the fat case. 29 The results for fat strings reveal an important feature: As the log increases the fraction of the total energy in radial modes decreases. Lattice spacing effects become increasingly significant, so this behaviour is only seen after the continuum extrapolation. Systematic errors from the lattice spacing also have a significant (and possibly even greater) effect for simulations of physical strings. Even though such simulations have better resolution prior to the final time, lattice spacing effects create a fake increase from log = 6, which is shallower for the data with better resolution. Meanwhile, we will see in Appendix B.4 that the slight difference between the initial values of ρ r /ρ tot for the two resolutions for the physical strings is due to a small difference in the initial conditions.
To sustain the scaling regime, energy density in strings is continuously emitted into axions and radial modes. We denote the emission rates, i.e. the energy released per unit time, by Γ a and Γ r extrapolation m r Δ =1  Figure 11: The ratio r a between the instantaneous energy emission rate from strings to axions and the total emission rate (i.e. to both axions and radial modes), as a function of log for different lattice spacings for fat strings. The continuum extrapolation shows that the network increasingly emits energy to axions (rather than to radial modes) as the log increases. This is consistent with the expectation that heavy modes decouple in the large log limit.
respectively. These can be be computed from where i = a, r, and z a = z r = 4 for fat strings due to the dependence of the radial mode mass on time.
In Figure 11 we plot the proportion of the total energy that is instantaneously emitted from strings that goes into axions, i.e. r a ≡ Γ a / (Γ a + Γ r ), for the fat string system. Similarly to the energy in radial modes, the rate of emission into axions is increasingly sensitive to the lattice resolution as the log grows (with r a larger for finer lattices). The continuum extrapolation suggests that r a increases steadily.
Together, the log dependences identified above indicate that the radial mode plays a decreasing role in the dynamics at late times in simulations. They also suggest that it should decouple in the limit log → ∞. Further, since only axions will get excited in this limit, the details of the particular UV physics that gives rise to the axion field would be unimportant for the dynamics of the strings.

B.4 The Spectrum
As discussed in Section 2, the axion energy density spectrum, and its dependence on log, plays a key role in determining the axion number density when its mass becomes cosmologically relevant. To define the axion spectrum we start from the expression for the axion energy density ρ a = ȧ 2 , where x p = R(t)x are physical coordinates, andã(k) is the Fourier transform ofȧ(x p ). The axion spectrum ∂ρ a /∂|k| is then fixed by requiring d|k| ∂ρ a /∂|k| = ρ a , and is therefore given by where Ω k is the solid angle. In order to exclude strings, the field a(x) needs to be screened. We substituteȧ(x) →ȧ scr (x) ≡ 1 + r(x) fa ȧ(x) in eq. (13), since the factor 1 + r(x) fa = |φ| |φ| r→0 automatically vanishes inside string cores and tends to unity far from strings. 30 In Figure 12 we plot the axion spectrum ∂ρ a /∂k at different times for fat and physical strings. The initial conditions are close to the scaling solution, and correspond to the black lines in Figures 1  and 6. In both cases, the spectrum has a peak at momenta of order 5H ÷ 10H and at smaller k it is power-law suppressed ∝ k 3 . Moreover, the spectrum has a UV-cutoff at momenta k = m r /2, above which it is highly suppressed. The shape of the spectrum remains similar as time passes, modulo the shift in the UV-cutoff.
At the momentum k = m r /2 there is a small peak, which we attribute to the energy exchange between axions and radial modes via parametric resonance. Such an effect can be understood  Figure 13: The energy density spectrum of radial modes for the physical string system as a function of the comoving momentum k com ≡ k m r /H at different times. Solid lines represent results from simulations with m r ∆ f = 1, and dashed lines with m r ∆ f = 1.5. The spectrum is dominated by an IR peak that comes from the initial conditions, with modes at larger comoving momenta generated during the evolution. We also indicate the comoving momentum that corresponds to m r at each of the times plotted.
heuristically. From eq. (8) in flat spacetime, the axion equation of motion is (1 + σ) ∂ µ ∂ µ θ + 2∂ µ σ∂ µ θ = 0, where θ ≡ a/f a and σ ≡ r/(f a / √ 2). 31 In the presence of a spatially homogeneous radial mode σ = σ 0 sin(m r t) with σ 0 < 1, the axion Fourier modes therefore satisfyθ 2 k + k 2 θ k + 2σ 0 m r cos(m r t)θ k = 0, where we kept the first non-vanishing dependence on σ 0 . If σ 0 = 0, this is a parametric resonance equation for the mode k = m r /2. A similar effect also occurs for a nonhomogeneous radial field, as it is in simulations. For fat strings, m r decreases proportionally to the scale factor and the parametric resonance affects a unique comoving momentum at all times, as seen in Figure 12 (left). On the other hand, for physical strings a wide range of comoving momenta are affected, and the resulting oscillations cover almost the entire spectrum. As discussed in Appendix A, this effect is easily triggered if the strings in the initial conditions do not have a core-size equal to m −1 r . In this case, radial modes will be emitted while the string core size is adjusting, and these will produce axions. 32 As the string cores relax the rate of such emission will decrease, and the parametric resonance effect will gradually disappear. This can be seen from the reduction in the amplitude of the oscillations in the final time shots plotted in Figure 12.
We also compute the energy density spectrum of radial modes ∂ρ r /∂k. Since radial modes behave as massive waves, this can be defined similarly to the axion spectrum as dk∂ρ r /∂k ≡ ṙ 2 , so 31 See also eq. (42) with m a = 0. 32 Of course in this case the radial mode will be space-dependent and the discussion above does not strictly apply.
To avoid strings in the determination of the radial spectrum we adopt the same masking technique as for the axion, i.e. in eq. (15) we substituteṙ(x) →ṙ scr (x) ≡ 1 + r(x) fa ṙ(x).
In Figure 13 we plot the spectrum of radial modes ∂ρ r /∂k for physical strings at three different times and for two lattice spacings, m r ∆ f = 1 and 1.5. The spectrum is plotted as a function of the comoving momentum k com ≡ k m r /H and we also indicate the momentum corresponding to the radial mode mass at each time. Figure 13 reveals interesting features of the evolution of radial modes. First, the spectrum is peaked at a fixed comoving momentum, corresponding to k = m r at the time H = m r , and it has a sharp fall off at momenta bigger than m r at any given time. The average momentum is therefore smaller than m r at all times and radial modes are on average nonrelativistic. Second, the height of the peak decreases proportionally to H ∝ R −2 , i.e. as in the free nonrelativistic limit. This shows that this peak is entirely produced at early times and does not receive contributions afterwards. Therefore the slight differences in ρ r /ρ tot at low momenta observed in Figure 10 for the two lattice spacings are only due to the initial conditions. Finally, the lattice spacing effects at high momenta grow at larger logs, producing a fake rise in simulations with the worse resolution. This is in turn related to the fake growth of the ratio ρ r /ρ tot in Figure 10 at late times.

B.4.1 The Instantaneous Emission Spectrum
We now turn to study the spectrum with which axions are instantaneously emitted by the network. To do so, we extract F from its definition F (k/H) = ∂ log(Γ)/∂(k/H) of Section 2.2 and express Γ in terms of the axion spectrum as in eq. (12). This leads to In the last equation A = H/Γ, which is a consequence of the normalization condition ∞ 0 F [x, y] dx = 1, which follows from its definition (see ref. [7] for more details). The time derivative in eq. (16) is calculated from simulation data by taking the difference of spectra with ∆ log = 0.25. 33 In Figure 14 we plot F at different times. We see that F [x, y] has IR and UV cutoffs at x = 5 ÷ 10 and x = y/2 for all y (i.e. at all times). In between it is well approximated with a power law 1/x q . Moreover, similarly to the total spectrum, F has significant fluctuations at the momenta affected by the previously described resonance, which for physical strings encompasses a large range of x (as discussed these do not represent genuine emission from strings, but are an unphysical effect that will disappear in the large log limit). 34 Finally, in Figure 14 a change in the power-law q with the log can be seen by eye. This corresponds to the evolution of q plotted in Figure 2 of the main text. The momentum range [30H, m r /4] over which q is calculated in Figure 2 is highlighted in Figure 14. The increase in q is even clearer in Figure 15 where we plot xF . At the final simulation time the instantaneous emission is not far from reaching q = 1.  17)). An increase in the slope q with log can be seen for both the fat and physical string systems. The highlighted region corresponds to the data points in the range [30H, m r /4], which we use for the fit of q in Section 2. We also indicate m r /(2H) at each time, above which emission is highly suppressed. It is reasonable to expect that at log 8 the slope q will overtake the value q = 1.
We now analyze the possible functional form of q(log) of Figure 2. The result of the linear and the quadratic fits are shown in Figure 2 in dark and light orange respectively. The two fits are both quite good for the fat and also the physical case. However, the result from the quadratic fit is compatible with the linear fit, but with larger errors. This suggests that a linear fit is enough to reproduce the data. For the linear fit q = q 1 log +q 2 , we get q 1 phys = 0.053(5) q 2 phys = 0.51 (7) , q 1 fat = 0.084(2) q 2 fat = 0.28 (2) .
In both the physical and the fat string systems the fit return values of q larger than unity for log 9, which might be accessible with future generation simulations. In particular, for physical strings, the fit gives q phys (log → 70) = 4.1(5).

B.4.2 Systematics
Given the importance of q for the final axion abundance, we now analyze potential systematic errors in its determination in detail.
Lattice spacing effects. For the fat case, we performed multiple sets of simulations with different resolutions (from m r ∆ = 1/1.8 to m r ∆ = 1), and in Figure 16 (left) we plot q as a function of log for each set. We calculate q by fitting the slope of log F (as a function of log x) in each simulation over the momentum range [30H, m r /4], and then averaging over simulations with the same resolution. The result for q is compatible for all lattice spacings. Consequently, in Figure 2 in the main text we report the average of all the sets. The largest values of log, i.e. 7.7 < log < 7.9, can be explored only with the least conservative lattice spacing (m r ∆ = 1), and therefore for those we have no direct comparison. However, given the good agreement for smaller logs, we expect that these values are in the continuum limit as well.
In the physical case we performed two sets of simulations, with final values of m r ∆ f = 1 and 1.5 (at earlier times the resolution is better). The result for q (calculated as before) is shown in Figure 16 (right). Due to the parametric resonance effect, q has unphysical fluctuations at small log. However, for log > 6 the fluctuations are relatively minor and a growth in q is clear for both resolutions. Moreover, the results for the two lattice spacings are mostly compatible with each other. Nevertheless, given the slight difference, in the main text we reported only the results with m r ∆ f = 1, which is the more conservative.
UV and IR cutoffs. In extracting q, the extremes k IR and k UV of the momentum range fitted [k IR , k UV ] need to be sufficiently far from the Hubble peak and the UV cutoff respectively. By construction, q is therefore less prone to lattice spacing and finite volume systematics compared to other quantities, such as energy densities (which are more UV sensitive) and number densities (which are more IR sensitive). For the fat string system, in Figure 17 we plot q at different times as k IR and k UV are varied. This shows that q has already converged to its true value with the choice [30H, m r /4]. varying k IR , with k UV = m r /4 fixed. Right: The best fit q with varying k UV and k IR = 30H fixed. It can be seen that k IR = 30H and k UV = m r /4 are sufficient for q to have convergered to its true value. This analysis has been done for m r ∆ = 1 and the times t 1 , t 2 , t 3 correspond to log = 6.5, 7, 7.5. In the physical case the best fit value of q has a small dependence on the momentum range used. Figure 2 (left) in the main text shows the fit in [30H, m r /4] for m r ∆ = 1. In Figure 18, we show different fits choosing k IR between 30H and 50H and k UV between m r /6 and m r /4, for the two different lattice spacings m r ∆ f = 1 and 1.5. While the numerical value of q at a particular log changes slightly, the choice of the momentum range does not change the trend of q increasing with log. Indeed, in the same plots we also show the linear fit of q as a function of log, and all of these have a positive gradient. We also show the quadratic fits, which are compatible with the linear fits but with much larger uncertainties.
Finite Volume. As mentioned in Appendix A, the simulations have been run until HL = 1.5. In Figure 19 we show F for different choices of HL at log = 6.4 for the fat string system. While for HL = 1.5 the IR part of the spectrum gets modified compared to HL = 2, in the momentum range where q is extracted the two are completely compatible. Although this is shown only for log = 6.4, finite volume effects are not expected to depend strongly on log (and are also expected to be similar for the physical string system). Thus, the choice HL = 1.5 will not introduce a significant systematic error in the calculation of q.

C Log Violations in Single Loop Dynamics
In this Appendix we give additional evidence of logarithmic violations in the dynamics of global string by studying the collapse of circular loops in flat spacetime. Sub-horizon loops make up a small percentage of the total string length during the scaling solution [7]. However, we will see that they can still give useful insights into the properties of the scaling regime, especially in relation to the Nambu-Goto limit.
The dynamics of global string loops can be described by an effective theory in which the fundamental degrees of freedom are the string and the axion radiation [20], with an interaction governed by the Kalb-Ramond action [85]. This theory is valid in the regime log(m r R) 1, where R is the typical loop size, i.e. when the string and the emitted radiation (with frequency ω ∼ 1/R) are not strongly coupled. In particular, this action does not capture the dynamics when strings intersect. As shown in [20], in this theory the coupling of the axion to the string is proportional to 1/ log(m r R). So, in the limit log(m r R) → ∞, the axion radiation decouples from the string. The string loop would therefore behave as in the free (Nambu-Goto) limit, oscillating an infinite number of times. As described in [35], this suggests that for finite but still large log(m r R) the loop might bounce many times before disappearing, emitting radiation with a typical wavelength of the order of the initial loop size. It has been argued in [35] that this will lead to a spectrum with q > 1. However, the previous argument is not definitive because the effective theory breaks down when the loop has shrunk to a small size, and the dynamics when the loop is small are critical in determining whether it bounces.
A complete analysis of the evolution of a string loop, including the bounce, can be carried out by solving the full field equations with the heavy radial mode present. We numerically solved eq. (8) in Minkowski spacetime with initial conditions φ(x) andφ(x) that resemble a static circular loop with initial radius R 0 . Limitations on the gridsize require log(m r R 0 ) 5. In Figure 20 we plot the loop radius R(t) (normalized to its initial value R 0 ) as a function of time for different log(m r R 0 ). We also show the free Nambu-Goto time law, R NG (t) = R 0 cos(t/R 0 ). Figure 20 has a number of interesting features. First, as log(m r R 0 ) increases, R(t) gets closer to the prediction for free strings. 35 As a result, the relativistic boost factor increases with log(m r R 0 ), and tends to infinity in the limit log(m r R 0 ) → ∞. By taking the time derivative of R(t) one can lines correspond to different values of the initial log(m r R 0 ). The dashed red line is the free Nambu-Goto solution R N G (t) = R 0 cos(t/R 0 ). As log(m r R 0 ) grows, the time-law R(t) approaches the Nambu-Goto limit, and the loop tends to bounce more. This is in agreement with the picture proposed in [35]. easily show that, for instance, the boost factor is already of order 10 for log(m r R 0 ) = 5. Second, the loop tends to bounce more for increasing log(m r R 0 ): for example, a loop with log(m r R 0 ) = 5 oscillates producing a loop with log(m r R 0 ) ≈ 4, which subsequently bounces to a loop with log(m r R 0 ) ≈ 2. The larger bounce is likely to be related to the increased boost factor, because a relativistic string loop is less likely to release all its energy at once before disappearing. 36 We also note that correctly evolving strings with a large boost factor requires a fine lattice spacing to resolve the relativistically contracted core and the bounce. For instance, the simulations need to be performed with m r ∆ = 1/20 or smaller for log(m r R 0 ) = 5, otherwise the loop will unphysically collapse as soon as it approaches its center. 37 After the loop disappears, the energy will be released into axions and heavy radial modes. As log(m r R 0 ) increases, we checked that the percentage of the energy emitted in axions gets larger, again pointing to the decoupling of the radial mode in the large log(m r R 0 ) limit. 38 Although not definitive because they are done at small logs, these results support the picture proposed in ref. [35], and agree with the general discussion of the previous Appendix.

D The End of the Scaling Regime
The scaling regime of Section 2 holds in the limit of vanishing axion potential, and it ends once the axion mass becomes cosmologically relevant, which happens when H and m a are of the same order. In this Section we make the previous approximate expectation more precise, and show that the scaling regime is not affected by m a = 0 until H > m a . Indeed we identify H = m a as the point at which the scaling regime begins to break down.
The full Lagrangian of the system is the same as in Section 2 with the addition of a potential for the axion, which we take of the form where φ = r+fa √ 2 e i a fa and m 2 a is the time dependent axion mass introduced in eq. (4). The equation of motion from eq. (18)φ does not depend on f a directly. Instead, the dependence on the two scales f a and m r can be reabsorbed by rescaling the field φ → φf a and the space-time coordinates t → t/m r and x → x/m r . Therefore, up to a trivial field rescaling, the physics is only sensitive to the two ratios m r /H = 2m r t and m r /m a , which we will refer to in the following. We implement eq. (19) numerically in the fat string system, described in Appendix A. As in Section 3, we define H to be the Hubble parameter at the time when H = m a , and correspondingly m a = H (H /H) α/4 and log ≡ log(m r /H ), where m r ≡ m r (t ). The choice of log = log(m r /m a ) determines the scale separation at H = H , sets the time at which the axion mass becomes relevant, and fixes the axion mass in units of m r in eq. (19). Although log ≈ 60÷70 for physical axion masses, we can study the breaking of the scaling regime only for log 6. To do so we analyze when the evolution of the system with finite log starts deviating from the evolution with the same initial conditions but with m a = 0 throughout.
In Figure 21 (left) we show the evolution of the scaling parameter ξ with time for axion masses with α = 7 and different choices of log (solid curves). We also show the evolution for vanishing axion mass (dashed curves). The initial conditions are the same in all of the simulations and are on the scaling solution. Figure 21 (right) shows the relative difference between the curves with and without the axion mass as a function of time. The choice of H /H on the x-axis makes it clear that the critical Hubble at which ξ starts differing by more than one percent is H crit = H for all log . We stopped evolving these simulations at m a /H 10, since at later times systematic errors due to an insufficient hierarchy between m a and m r can become significant (discussed in Appendix F). We have checked however that at later times the effect of a nonzero axion mass is to reduce the string length, until the whole network disappears.
In Figure 22 we show the axion energy spectra at different times (i.e. different H /H) for a system with an axion mass such that log = 5 with α = 7. We again plot the results obtained . The x-axis is normalized with H , so it can be seen that regardless of log , ξ is unaffected by the mass before H = H . Right: The relative difference between the results for ξ with and without a non-zero axion mass, δξ/ξ, as a function of time. We define δξ ≡ (ξ − ξ ma=0 ), and ξ ma=0 is the string length in a simulation with m a = 0 throughout. For all values of log tested the effect of a non-vanishing axion mass is smaller than percent for H > H .
from a system with m a = 0 (dashed lines) for comparison. While H > H there is no significant difference between the two spectra, and at H < H the spectrum starts differing substantially. We note that a nonzero axion mass affects the spectrum earlier than it affects ξ. In particular the IR modes, which contribute the most to the axion number density, are affected first and change soon after H = H . At later times UV modes are also affected. The results obtained are similar for other values of log . We have checked that the behavior seen in Figures 21 and 22 holds for all values 0 ≤ α ≤ 8. As expected, for larger α the mass affects ξ and the spectrum even less while H > H , but the network shrinks within fewer Hubble times after this point.
We therefore conclude that before the critical value H crit = H the string system is not affected by a non-vanishing axion mass. The accuracy of our results is not sufficient to establish whether H crit itself has a small residual dependence on log , and indeed this is plausible considering the log violations discussed in Appendix B. Possible violations are expected to be small since the axion mass changes rapidly, and would therefore not significantly affect the lower bound on the axion abundance calculated in Section 3.
As discussed in Section 4, we do not attempt to calculate the axion number density that is emitted by the string and domain wall network as the latter is destroyed after m a = H. Indeed, we do not expect that a direct calculation in simulations would reproduce the number density in the physically relevant regime. In particular, at the values of log accessible in simulations the strings are still emitting a UV dominated spectrum and have a tension and density ξ that is far from the physically relevant values. Before H = H the non-zero axion mass has a negligible effect on the spectrum, but subsequently this gets modified starting starting from IR momenta.

E Axions through the Nonlinear Regime
In this Appendix we will give further details on the derivation of the analytic prediction for the axion number density after its potential becomes relevant. We also describe the simulations that we performed to confirm its validity and fit its free parameters.

E.1 Derivation of the Analytic Estimate
The initial conditions for the axion equations of motion of the Hamiltonian in eq. (4), correspond to a superposition of waves with energy spectrum ∂ρa ∂k described in detail in Appendix B.4, emitted by strings during the scaling regime prior to H = H . As discussed below, the initial field can be obtained inverting eq. (16) as a function of ξ and F at H = H . In doing so we assume that at large log the energy from strings is emitted purely into axions (as the results in Appendix B.3 suggest). The extrapolation of the effective string tension µ eff is also needed. As shown in Appendix B.3, µ eff in simulations is reproduced well by the theoretical expectation µ th = γ πf 2 a log mr H ηc √ ξ with η c = 1/ √ 4π and γ ≈ 1.3 constant in time, and we assume that this remains approximately true also at large scale separations.
The actual form of F , shown in Fig. 14, has a nontrivial shape. For simplicity here we approximate it with a single power law q and a sharp IR cutoff at x < x 0 , As we will see, the results obtained from this simple form capture the main features of the dynamics. We also considered more complex shapes reproducing the one in Fig. 14 more closely. However, when compared to numerical simulations, the improvement with respect to the simplified approximation in eq. (21) is negligible once compared to the uncertainties induced by the large log extrapolation. By inverting eq. (16), the total energy in axions will be distributed with the spectrum where the primed quantities are computed at the time t , the redshifted momentum is defined as k = kR/R (see eq. (23) of [7] for the explicit derivation). As mentioned in Section 2.2, the emission rate Γ ξµ eff /t 3 8πH 3 f 2 a ξ log(m r /H) is fixed by energy conservation, and we assume a linear logarithmic growth of ξ, as in eq. (3), and q > 1. Using eq. (21), the resulting convoluted spectrum at H = H for momenta k < where k 0 = x 0 H . For k > x 0 √ H m r the spectrum falls faster than 1/k and its precise form is not important since it gives a negligible contribution to the abundance. To a good approximation, we can neglect the effect of γ, η c and the order one factor between µ eff and µ th of Figure 9, so we take µ ≈ πf 2 a log . 40 The terms in the last line of eq. (23) are O(1/ log ) and can be neglected in the large log limit, so in this limit all the dependence of ∂ρa ∂k on q is also suppressed.
Due to the scaling regime, the leading dependence of the spectrum for k > x 0 H is ∂ρa ∂k ∝ 1/k for all q ≥ 1 (i.e. the spectrum obtained after convoluting F is scale invariant). Correspondingly, the energy is distributed equally in logarithmic intervals between the momenta x 0 H and √ H m r .
The logarithmic dependence of ξ and µ on time induces violations of the scale invariance that are proportional to log 2 (k/k 0 ). At least up until H = H , away from strings axions propagate as free waves, and their spectrum can therefore be used to infer the axion field itself via the relations (valid for relativistic waves) where k = | k| andã( k) is the Fourier transform of a(t , x) and V is the volume. From the last equality of eq. (24) it also follows that the average square amplitude is Combining this with eq. (23) we obtain a 2 | = 4ξ µ for log 1 and ξ log 1. Consequently, as mentioned, a 2 | is much larger than f 2 a in this limit. Following the procedure sketched in Section 3, we can now derive the formula for the contribution to the final abundance from the spectrum of eq. (23). In particular, in the limit of large log , the effects from the axion potential can be neglected 41 also after t and the spectrum evolves relativistically as This evolution holds up until the contribution from the axion potential in the Hamiltonian (ρ V = m 2 a (t)f 2 a ) becomes of the same order as the gradient one from IR nonrelativistic modes, i.e. until t = t when the following condition is satisfied: where c V and c m are O(1) coefficients to be determined numerically.
The condition above provides the following implicit equation for t , or equivalently for m a (t )/H : where all the quantities are evaluated at t = t . We note in particular that ρ IR is still much larger than m 2 a f 2 a at H = H because of the enhancement by a factor of ξ µ ∝ log 2 . Introducing the quantity which measures the delay of the nonlinear regime induced by the ξ log enhancement, eq. (28) can be rewritten as 8πξ log log(κ) 1 − 2 log(κ) log + 4 3 The equation can be further simplified by noticing that the first term dominates in the limit log 1, and thus we get Eq. (31) can be solved analytically using the identity a log(bz c ) = z ⇐⇒ z = −acW k (−(acb 1/c ) −1 ) for some k ∈ Z, where W k (z) is the Lambert W -function evaluated on the k-th Riemann sheet and defined by ze z = a ⇐⇒ z = W k (a). The solution is where the choice of the lower branch k = −1 is dictated by the fact that the argument of W is negative (because c m , c V > 0) and the value of W k in eq. (32) must be negative (and large) so that z > 1. In the limit ξ log 1 we can expand W −1 for small negative arguments. Noticing that where the second equality is just the recursion of the first equality and the dots stand for infinitely nested logs, eq. (32) gives where the logarithms are infinitely nested. At t = t the field dynamics is completely nonlinear with most of the spatial gradients of order the axion mass or lower and energy density of order m 2 a (t )f 2 a . As the Universe continues to expand, the energy density and the field value decrease further, and so do nonlinearities. We assume that the transient lasts O(1) Hubble times, so that during this period the total energy is approximately conserved. After the nonlinear transient, the axion field drops below f a , the dynamics is mostly linear and the comoving number density is conserved again. The latter can therefore be derived from the energy density at t as where the O(1) coefficient c n , which we will fit from numerical simulations, parametrizes all matching effects during the nonlinear transient, such as the O(1) effects from the redshift during the transient and the extra contribution from slightly relativistic modes above c m m a (t ).
We finally arrive to an expression for the relative contribution of relic axions from strings during the scaling regime normalized to the reference misalignment value at θ 0 = 1 (i.e. n mis,θ 0 =1 Conservation of the comoving number density implies that Q(t) is constant for t t , so that one can easily compute the number density of axions today, n str a (t 0 ), from the equation above. From the analytic expression in eq. (36) we can notice several important features. First, as a result of the large log approximation, the explicit dependence on q has disappeared (we will show below from the full numerical results how indeed such a dependence is subleading). Second, the dependence on x 0 (the effective IR scale of the spectrum) is only logarithmic. This softens considerably systematic errors from neglecting a possible evolution of x 0 during the scaling regime, and makes manifest the insensitivity of the final result on the details of the shape of the IR part of the spectrum. Third, the unknown parameters c V,m,n enter the final formula as multiplicative factors (or in the logarithmic dependence), so that the functional dependence on ξ log , α and x 0 is a true prediction of the above analytic derivation, which is confirmed by the numerical computations. Notice finally that with α 1, i.e. when the axion potential grows fast and the redshift effects between t and t can be neglected, the formula aboves simplifies further, recovering the simple dependence Q(t ) ∝ (ξ log ) 1/2 , anticipated in the main text. In Appendix E.3 we will see that eq. (36) fits well the numerical results and we will provide numerical fits for the coefficients c V,m,n (see the caption of Fig. 24), which indeed are of order one.

E.2 Setup of the Numerical Simulation
In Section 3 we also studied the evolution of the axion field numerically by solving eq. (20) on a discrete lattice. Unlike our simulations of the scaling regime, eq. (20) only contains the axion and not the radial mode. The absence of the radial mode means that strings are automatically absent, and simulations can reach values of H that are arbitrarily smaller than f a and m r .
The numerical implementation of eq. (20) is very similar to the string system. In particular, it is convenient to use the rescaled field ψ = R(t) a/f a and the conformal time τ , which is defined as In this way eq. (20) simplifies to where ψ and∇ x ψ are derivatives with respect to the dimensionless variables H τ and the comoving coordinate H x, and τ = τ (t ). We solve eq. (38) numerically starting from τ = τ in a box with periodic boundary conditions. Space is discretized on a cubic lattice of comoving side length L c containing up to N 3 x = 3000 3 uniformly distributed grid points. The space-step between grid points in comoving coordinates ∆ c = L c /N x is constant in time. Eq. (38) is discretized following a standard central-difference Leapfrog algorithm for wave-like partial differential equations, and evolved in fixed steps of conformal time ∆ τ , which we fix to ∆ τ = ∆ c /3. 42 The derivatives are expanded to fourth order in the space-step and second order in the time-step.
The physical length of the box is L(t) = L c R(t) and the physical space-step between grid points ∆(t) = L/N x = ∆ c R(t) grows ∝ t 1/2 . Therefore, the number of Hubble patches in the box HL decreases with time, and the axion mass in lattice units m a ∆ rapidly increases. This leads to potential systematic uncertainties in the results, which we analyze in Section E.4. Our final results are obtained with parameter choices that are free from significant uncertainties.
In our analysis we used two different functional forms for the initial spectrum ∂ρa ∂k . The first is that of eq. (23), which is derived from the simplified F consisting of a single power law. For the final results in the main text, we used a more realistic spectrum obtained by inserting a better approximation to the true shape of F (plotted in Figure 14). This comprises four different power laws, chosen to be x 3 , x 1/2 , x −1/2 , x −q joined at the points x = x 0 /4, x 0 , 2x 0 , y, and with F [x > y, y] = 0. The intermediate power laws in this reproduce the broad IR peak in F . For both functional forms we used q = 5, as suggested by the extrapolation of q(log ) of eq. (17) at log = 70. As shown below, provided q 1 its precise value is not important for large ξ log , so this choice does not have any significant effect on the results obtained. 43 42 As shown in ref. [7] this time discretization is sufficient to avoid numerical effects at the per-mille level. 43 The initial spectrum was always generated for log = 70, and different ξ log are obtained by varying ξ (so,

E.3 Further Results from Simulations and Oscillons
To understand the dynamics as the axion mass becomes cosmologically relevant it is useful to study the mean square amplitude a 2 defined in eq. (25), the axion spectrum ∂ρa ∂k and the axion number density In Figure 23 we show the time evolution of the mean square axion field amplitude a 2 and of the comoving number density Q(t) = n str a (t)/n mis,θ 0 =1 a (t). The results are shown for different values of ξ log , and with the simplified initial axion spectrum with IR cutoff x 0 = 10 and with α = 8. As mentioned in Section 3, for sufficiently large ξ log the early evolution of these quantities matches the expectation of a relativistic regime. Calculating n str a makes sense at these times since the whole potential is negligible in the Hamiltonian density (4), which is diagonal in momentum space. After a transient period when nonlinearities dominate the evolution of the system and n str a is not defined (corresponding to the dashed lines in Figure 23), the average amplitude becomes smaller than f a . At such times the majority of the field is in the linear regime, with the exception of objects called oscillons that are produced.
An oscillon is a localized, metastable, time-dependent solution eq. (20), in which the field oscillates with maximum amplitude of order πf a within a region of (decreasing) size m −1 a [45]. It therefore contains an energy density of order 1 2 m 2 a f 2 a π 2 . Although oscillons decay, radiating their energy into axions, they are thought to be very long lived [90,91]. Indeed, we observe that they do not disappear within the range of our simulations. As time passes, however, they occupy an increasingly negligible proportion of space. Consequently it makes sense to calculate the number density of axions by considering only the field far away from oscillons, where the linear approximation to eq. (39) is valid.
We screen oscillons by multiplying the axion field and its time derivative by a window function w(ρ a ) of the local energy density ρ a (x) (defined as in eq. (4)), i.e. a scr (x) ≡ w(ρ a (x))a(x). We choose the window function to be which vanishes for ρ a 1 2 m 2 a f 2 a π 2 and tends to 1 for ρ a m 2 a f 2 a as required. We have checked the results are not sensitive to this particular functional form. The screened average amplitude, spectrum and number density are defined as in eqs. (25) and (39) substituting a(x) → a scr (x) and similarly for the time derivative.
In Figure 23 we plot the mean field amplitude and the axion number density at late times both with and without oscillon screening (respectively, the lower and upper curves between the shaded regions). 44 At earlier times (i.e. during the transient) we do not show the results with oscillons screened since they are not yet clearly defined objected. As expected, with oscillons screened the system reaches a regime in which the comoving number density is conserved, and the average field amplitude decreases as in the limit of nonrelativistic dynamics. 45 Moreover, the effect of oscillons decreases at late times. We have tested that the results at the end of simulations are not sensitive to the minimum energy density that is masked by the window function, provided this is in a reasonable range.
Despite their interesting dynamics, we refrain from analyzing the evolution of the oscillons in detail. Instead we just note a few relevant facts. We estimated the number of oscillons in two ways: (1) by dividing the total volume in oscillons (defined as the points where w(ρ a ) > 1/2) by the volume of one oscillons m −3 a , and (2) by dividing the total energy in oscillons by the energy of one oscillon (m −3 a × 1 2 m 2 a f 2 a π 2 ). The number of oscillons computed in the two ways is compatible and constant in time. This means that, once formed, oscillons do not decay within the simulation time. Additionally, the number of oscillons that form depends on the initial amplitude of the field and increases if the value of ξ log in the initial conditions is larger.
Our analytic estimate for Q has been derived only for the simplified spectrum in eq. (21). However, this differs from the more physical spectrum only in its IR part. As argued before, Q is not very sensitive to this part of the spectrum, so that we can apply the result from eq. (36) to the more physical spectrum. We expect the different shape of the initial spectrum to be reabsorbed in the numerical fit of the coefficients c m , c V , c n .
In Figure 24 and Figure 3 in the main text, we show the comoving number density of axions Q(t f ) at the final simulation time t f as a function of ξ log for α = 4, 6, 8 and x 0 = 5, 10, 30 for the physical initial spectrum. The errors are systematic and come from the uncertainty in the screening of oscillons, which we estimate as the difference between the masked and unmasked number density. The statistical errors, as well as the systematic errors from finite volume and finite UV-cutoff are subdominant. The continuous lines in the same plot represent the analytic estimate in eq. (36) (valid for ξ log 1), where the coefficients c m , c V , c n have been fixed with a global fit of all the data points with ξ log ≥ 100 in Figures 24 and 3. We note that the fit is good, and the dependence on ξ log , α and x 0 in the numerical data is captured well by the analytic result. Equivalent plots for the simplified initial spectrum show a similarly good fit.
Finally, in Figure 25 we show the time evolution of the spectrum ∂ρa ∂k for ξ log = 10 3 with α = 8 and x 0 = 10. As expected, the spectrum evolves as in the free relativistic limit initially, and the nonlinear transient depletes the amplitude of modes k < m a (t ). Oscillons affect the spectrum only at momenta of order of the axion mass at a given time, indicated with an empty dot.
For all of the results presented so far we have fixed q = 5. However, there is some uncertainty in our extrapolation of q. In Figure 26, we show the number density of axions as a function of q before and after the nonlinear regime (at t = t ) for ξ log = 10 3 , x 0 = 10 and α = 8. It is clear that provided q 1, its actual value introduces only a very minor uncertainty on the final axion number density.  Figure 26: The axion number density from strings relative to that from misalignment, as a function of the power law q of the axion emission spectrum during scaling. We have fixed ξ log = 10 3 , x 0 = 10 and α = 8. The results are shown before and after the axion mass becomes cosmologically relevant.

E.4 Systematics
Systematic uncertainties and numerical artifacts in the axion only simulations can arise from various sources. Here we describe the most important of these, and the choices of simulation parameters that were used to obtain our final results. First, we note that the number of Hubble patches in the box and the axion mass in lattice units are The former decreases with time, while the latter increases fast. Therefore the following sources of systematic uncertainty need to be considered.
• The continuum limit corresponds to m a ∆ → 0, and in particular if m a ∆ 1 lattice effects are introduced. All our simulations are stopped when m a ∆ = 1 so that the discretization effects are negligible for practically the whole simulation time.
• The infinite volume limit corresponds to HL → ∞. In Figure 27 (left) we show that the initial number density is free from finite volume effects provided H L ≥ 1, which matches expectations based on the form of F . Although during the subsequent evolution HL < H L , since no strings are present in the simulation, no new IR axion modes are produced and those modes that are present initially still fit in the simulation volume at later times. Therefore, volume effects do not affect the simulation at H < H . Instead, the systematic errors introduced by a finite H L even decrease over the course of a simulation. This is because the nonrelativistic regime tends to wash out the dependence on the IR shape of the spectrum, and this is the most affected by finite volume effects. From eq. (41) the time range H /H of a simulation (before m a ∆ = 1 is reached) is maximized by choosing the smallest viable H L . 46 Thus, all of our results are obtained from simulations that have H L ≥ 1.
• The maximum momentum of modes supported in a simulation is only k max ≡ 2πN x /L, which is far from the UV-cutoff of the (scale invariant part of the) spectrum, i.e. k ∼ √ H f a for the physical parameters. Therefore, given the scale invariance of the convoluted spectrum most of the kinetic energy density of the axion field will not be contained in the grid. However, given the discussion in Section 3, the number density n str a is dominated by momenta of order m a (t ) at t = t . Thus, if the UV-cutoff Λ U V of the grid satisfies Λ U V (H /H ) 1/2 m a (t ) the final number density is expected to be independent of Λ U V .
We tested the dependence of n str a on the UV-cutoff by generating initial conditions from the simplified spectrum, but setting ∂ρa ∂k = 0 for k > Λ U V for different Λ U V . In Figure 27 (left) we plot the axion number density as a function of time during its mass turn on for different Λ U V , relative to that of the largest value tested (the simulations are all identical apart from the value of Λ U V , e.g. they are on the same sized grid). We see that for α = 8, x 0 = 10 and ξ log = 10 3 the dependence of Q on Λ U V is negligible if Λ U V > 10 3 H , 46 A large time range is needed so that the axion field reaches the nonrelativistic regime, and the asymptotic value of n str a can be calculated. We do however note that as ξ log increases so does m a (t ) and therefore a larger UV-cutoff is required. The size of our grids is such that Λ U V > 10 3 H in all our simulations, which is sufficient to obtain accurate results.

F Massive Axions on a String Background
In Section 3 we studied the evolution of the axion field as the axion mass becomes cosmologically relevant considering only the axions emitted at earlier times, i.e. for H > H . In reality, until the string network is completely destroyed the axion field is made of different components: the axion radiation produced up to H by the scaling regime, axion radiation emitted at later times as the string network is destroyed, and strings themselves. Crucially, in the analysis of Section 3 we implicitly assumed that the presence of strings (which actually store an order one fraction of the energy density in their spatial gradients) does not influence the evolution of the axion radiation, or at least it does not weaken it. In this Appendix we show that the number density of axions from the scaling regime is indeed not affected, at least from the axion string network that can be currently simulated.

F.1 The Decoupling Limit
For a simulation to include strings, both the axion and the radial mode must be present as dynamical degrees of freedom, and therefore evolved e.g. following the Lagrangian in eq. (18) Figure 28: A comparison between the comoving axion number density from the homogeneous misalignment with a(t 0 ) = f a using the axion only equations (black) and the full complex scalar field equations with log(m r /m a ) = 7 (blue). The amplitude of the radial mode r/(f a / √ 2) is also shown. Already at m a /m r = 1/3, the radial mode is unphysically displaced from the minimum and the axion number density is not conserved by ≈ 20%.
the same time, the radial mode is increasingly displaced from its minimum, acquiring some of the energy that would otherwise be in the axion field. In other words, light enough degrees of freedom coupled to the axion get excited. This simple analysis shows that any prediction from a simulation where the decoupling limit is not reached will be strongly model dependent, and will not reproduce results in the physically important regime.

F.2 The Effect of Strings
Next we test whether the presence of strings affects the dynamics of the pre-existing radiation, and therefore the number density of axions resulting from the scaling regime. To do so we simulate the fat string system starting from its evolution during the scaling regime, through the time when the axion mass turns on. At H = H we modify the complex scalar field by injecting additional axion radiation, with the spectrum that would be produced at log = 60 ÷ 70 from the scaling regime with q > 1.
The extra radiation must be introduced to account for emission by the scaling regime at large values of the log and to enable a fair comparison with Section 3. In contrast, the radiation component emitted by the string network prior to log ≈ 7, which is directly accessible in simulations, is still UV-dominated. This will therefore not capture the dynamics that we are interested in, and will probably make a negligible contribution to the axion number density compared to misalignment. Of course, with this setup we do not aim to compute the complete axion relic abundance from strings and domain walls. Instead we simply want to confirm that our analysis in Section 3 is not significantly altered by the actual presence of strings and domain walls, and confirm that our analysis does indeed give a lower bound on the relic abundance.
To do so, we solved eq. (19) as in Appendix D starting from H = m r . Then when t = t , in the middle of the evolution, we substituted φ(t , x) → φ(t , x)e iaw(t , x)/fa andφ(t , x) → d dt φ(t, x)e iaw(t, x)/fa t=t .
The field a w (t ) is the additional radiation and is extracted from the kinetic energy spectrum ∂ρa ∂k of the scaling regime at log = 60 ÷ 70 (as in Appendix E.2).
There are a number of potential sources of systematic errors to be taken into account in such simulations, on top of those already present in the analysis of the scaling regime for vanishing axion mass.
• The ratio m r /m a should be large enough at the final time not to introduce the non-decoupling effects discussed in the previous Section. In particular, even if we will use m r ∆ = 1 as in the string network simulations, m a ∆ = 1 at the final time will not be sufficient. Instead we stop the simulations at m r /m a = 10 ( Figure 28 suggests that this is sufficiently large to avoid unphysical energy transfer to the radial mode, although it is not definitive since that plot is for a homogeneous field).
• H L should be large enough so that finite volume effects do not cause the string network to shrink. Instead we want this to occur due to the the axion mass. The resulting constraint on H L is stronger than that discussed in Appendix E.2). For α = 6, we checked that the scaling parameter has dropped by 50% due to the mass at H /H = 4, so choosing H L 4 is sufficiently safe for our purposes.
• The UV cut-off of the injected spectrum (which is scale invariant) should be smaller than m r /2 to satisfy the second requirement for the decoupling limit. We cutoff the injected spectrum for momenta bigger than k UV = 50H (as is described in Appendix E.4). For log = 7 this is sufficiently small compared to m r not to introduce major effects. 50 Accommodating the competing requirements of these sources of systematic errors is challenging. Simulations cannot last as long as the axion-only simulations of Section 3, and we cannot reach the regime where the comoving number density is precisely conserved (after the relativistic period and the nonlinear transient). This is the case even when a spectrum with relatively small ξ log = 200 is injected, for which the nonrelativistic regime is reached relatively early.
Despite these difficulties, results from simulations are still sufficient to show that the presence of strings does not affect the existing radiation, enough for our present work. In Figure 29 we plot (in blue) the evolution of the axion number density when the axion spectrum corresponding to ξ log = 200 is injected into a simulation with strings at the simulation time log = log = 7. 51 50 If substantially larger values of k UV are used the field evolution develops numerical instabilities. 51 The injected spectrum has a UV cutoff at k/H = 50 to maintain the required hierarchy k m r . We choose ξ * log * = 200, somewhat smaller than our central value for the QCD axion, since simulations with larger values require smaller lattice spacing, increasing the computational resources required.
As in Figure 28 we take α = 6. More precisely, we calculate the axion number density from the kinetic component of its energy, i.e. in the first term of eq. (39), multiplied by a factor of two. 52 This gives the correct result at the early and late times, and during the transient the axion number density is not well defined anyway.
For comparison, in Figure 29 we also plot in pink the number density calculated in an axiononly simulation of Section 3 with the same initial spectrum. 53 Additionally, we plot in orange the number density when the axion spectrum is injected into a simulation of the complex scalar field, but with no strings or preexisting radiation present. Finally, we show the number density arising from directly simulated strings when the axion spectrum is not injected.
The effect of the non-decoupling of the radial mode can be seen in these results. The final number density from simulations involving the radial mode is somewhat smaller than that from an axion only simulation, even when there are no strings or radiation present. This happens because we do not have the numerical power to simulate a very large hierarchy between m a and m r at such times. As expected the number density injected is much larger than that emitted by strings, during scaling and also as they are destroyed by the axion mass, at the values of log(m r /H) that can be simulated Most importantly, the final axion number density when the spectrum is injected into a simulation with strings is very close to that when it is injected into a simulation of the complex scalar with no strings. This indicates that the presence of strings, and the dynamics of their annihilation at the end of the scaling regime, does not significantly alter the evolution of the pre-existing axion radiation. Since the number density from the axions in the scaling regime is not depleted by strings with log = 7, there is no reason to expect that this will not also hold for larger log . Indeed is seems implausible that strings could absorb all of the energy in axions previously released by the scaling regime, and then emit this back as high momentum axions, which would be required for the strings to decrease the final axion number density.

G Comparisons
Our final conclusions in Section 4 differ from those obtained by other authors in the literature for various reasons. Here we comment only on those that are relatively close to us in their basic assumptions or techniques used.
• Based on the expected similarity between axion strings at log 1 and Nambu-Goto strings, the authors in [32][33][34][35] assumed values for ξ and q numerically compatible with those inferred by our study, and therefore got an enhanced contribution to the axion abundance from topological defects. Their analysis however did not take into account the effects from nonlinearities induced by the large axion field values. These, as we have shown, crucially  Figure 29: The evolution of the axion number density (extracted as twice the kinetic component of the number density) through the axion mass turn on, for α = 6. We compare simulations in which only the axion is dynamical (pink), starting with the initial field configuration predicted from the scaling regime with ξ * log * = 200, with results when the same axion spectrum is injected into a simulation of the axion and radial mode that has a string background evolved with a potential such that log = 7 (blue) (via eq. (19)). In the latter, the string network is in the scaling regime prior to the axion mass turn on, and is subsequently destroyed. We also plot the evolution of the same axion spectrum injected into a simulation of the axion and radial mode with no strings present (orange), and the result from the scaling regime without the addition of the extra axion field (light blue). The agreement between the results when the axion spectrum is evolved in complex scalar field simulations with and without strings shows that the presence of strings with log = 7 does not have a significant effect on the dynamics of the preexisting radiation.
affect the field evolution during the QCD transition and substantially change the final axion abundance.
• In [24] the authors performed a simulation of the entire axion string/domain-wall system's evolution from the scaling regime until the linear regime after the QCD transition and beyond. The result for the abundance is substantially smaller than the bound in eq. (6), however it is not in contradiction with it. Indeed the range of log(m r /H) that can be directly simulated does not allow large values of ξ log to be reached (and instead ξ log remains a couple of orders of magnitude below the physical one), nor does it allow the spectrum to be seen turning IR dominated (i.e. with q > 1). Consequently, the number of axions produced by strings will be small, thus the smaller abundance measured in such a simulation. In studying the effect of strings and walls in the evolution of the axion field during the nonlinear regime in Appendix F we also performed simulations analogous to those of ref. [24] obtaining compatible results. However from Fig. 29 in Appendix F.2 it is clear that ignoring a proper extrapolation of the parameters could easily lead to the dominant contribution to the abundance being missed and the total contribution being underestimated by more than one order of magnitude.
• As in the previous case, the authors of ref. [92] perform simulations of the entire evolution from the scaling regime to the linear regime including the decay of strings and walls. However, by changing the physics at the string core scale m r they manage to produce an effective string tension which is numerically equivalent to log 70. In doing so they also observe an enhancement of ξ which grows by a factor of a few. At t the system has an effective energy density prefactor ξ log of the same order of magnitude of our extrapolated one. Despite the large energy density, the final axion abundance found is small as if the string contribution was negligible. The result is not in obvious contradiction with our findings and can be understood as follows. Within the range of m r /H that can be simulated we observed that most of the energy of the string network is still dumped into UV modes, which have not yet decoupled and instead influence the IR string dynamics. Despite the higher string tension, the evolution of the string network in ref. [92] is probably still dominated by UV modes (which in such a setup are unphysical) producing a spectrum with q < 1 explaining the suppressed axion abundance. Unfortunately in order to check this interpretation of the disagreement a dedicated study of the axion spectrum in this setup is required, which is currently missing.