Axions from Strings: the Attractive Solution

We study the system of axion strings that forms in the early Universe if the Peccei-Quinn symmetry is restored after inflation. Using numerical simulations, we establish the existence of an asymptotic solution to which the system is attracted independently of the initial conditions. We study in detail the properties of this solution, including the average number of strings per Hubble patch, the distribution of loops and long strings, the way that different types of radiation are emitted, and the shape of the spectrum of axions produced. We find clear evidence of logarithmic violations of the scaling properties of the attractor solution. We also find that, while most of the axions are emitted with momenta of order Hubble, most of the axion energy density is contained in axions with energy of order the string core scale, at least in the parameter range available in the simulation. While such a spectrum would lead to a negligible number density of relic axions from strings when extrapolated to the physical parameter region, we show that the presence of small logarithmic corrections to the spectrum shape could completely alter such a conclusion. A detailed understanding of the evolution of the axion spectrum is therefore crucial for a reliable estimate of the relic axion abundance from strings.


Introduction
The QCD axion [1][2][3][4][5][6][7] is the simplest and most robust of the known solutions of the Standard Model (SM) Strong CP problem, and it also automatically forms a component of cold dark matter (DM) [8][9][10]. Consequently, a QCD axion that makes up the entire measured DM relic abundance is one of the best motivated scenarios for physics beyond the SM. In addition, numerous experiments aimed at detecting axions are currently running or in development. These will be sensitive to a substantial proportion of the relevant parameter space, and, if an axion is discovered, they could measure its mass and couplings precisely (see e.g. [11,12]).
The dynamics by which QCD axion DM is produced, and its final relic abundance, depends on the cosmological history of the Universe (see e.g. [13,14]). If the Peccei-Quinn (PQ) symmetry that gives rise to the axion was broken after inflation, the axion field had initially random fluctuations over the present day observable Universe. Instead, if the PQ symmetry was broken during inflation and never subsequently restored, the axion field was initially homogeneous. In this case the axion relic abundance is incalculable because it depends on the local value of the axion field after inflation. 1 In this paper we study the class of models in which PQ breaking happens after inflation. This includes theories in which inflation happened at a scale above the axion decay constant f a , and also those with inflation at a lower scale but which were reheated to a temperature above f a . 2 In this case, assuming a standard cosmological history, the relic abundance is calculable in terms of the axion mass due to the random initial conditions. 3 . As a result, there is in principle a unique calculable prediction for the axion mass if it is to make up the complete DM density in such models, which would be extremely valuable for experimental axion searches.
However, in this scenario the mechanisms by which DM axions are produced are complex, and calculating the relic abundance is challenging [18][19][20][21][22]. The random initial axion field after PQ symmetry breaking leads to the formation of axion strings [23][24][25]. These are topologically stable field configurations that wind around the U (1) vacuum manifold of the broken PQ symmetry as a loop in physical space is travelled. Interactions between strings are thought to result in the length of string per Hubble volume, measured in Hubble lengths, remaining approximately of O(1) as the Universe expands [26][27][28][29][30][31]. To maintain such a regime the string network must release energy. This dominantly happens through the production of axions, which form a potentially significant fraction of the total relic abundance [32][33][34][35][36][37][38][39].
The string system persists until the temperature of the Universe drops to around the QCD scale, when the axion mass turns on and becomes cosmologically relevant, leading to the formation of domain walls. The subsequent dynamics depends on the anomaly coefficient between QCD and the PQ symmetry, N W , which is equal to the number of minima that are generated in the axion potential [19,40]. N W = 1 corresponds to the scenario in which the domain walls are automatically unstable and decay, destroying the string network and releasing further DM axions in the process [19,40,41]. If N W > 1 the domain walls are generically stable, and the model is not phenomenologically viable, unless further explicit breaking of the PQ symmetry is introduced [42][43][44][45] (which might reintroduce the strong CP problem) or the Z N W symmetry is gauged [46]. In the following we set N W = 1 so that the PQ breaking scale v PQ = f a , however, since we will focus on the early evolution of the string network, the general case can be recovered from our results by simply replacing f a with v PQ = N W f a . After the QCD crossover the comoving axion number density is conserved, and the coherent axion field behaves as cold DM.
In our present work we consider the string network before the axion mass turns on. An understanding of this stage of its evolution is crucial, both to calculate the relic abundance of axions produced at such times and to set the appropriate initial conditions when analysing the system once the axion mass becomes relevant. The important properties of the network, which we aim to determine, include: the average length of string per Hubble volume, and the way that this is distributed in loops of different lengths; the rate of energy release into axions; and the spectrum of axions emitted. Due to the complex evolution and interactions of the strings, and later the domain walls, an accurate analytic calculation of the axion relic abundance appears hopeless in this scenario, and instead some form of numerical simulation is required. The most direct approach, which we will follow, is to simulate a complete UV theory that gives rise to the axion numerically on a discrete lattice [47] (more details are given in Section 2).
The properties of the axion strings are critically affected by the fact that they come from a global symmetry. As we discuss further in Sections 2 and 3, this means that the energy per unit length of an isolated string is logarithmically divergent. In the early Universe the divergence is cut off by the distance to a string of opposite orientation, or for a small loop its diameter, both of which are typically of order of the Hubble length H −1 [18]. Consequently, the theoretical expectation for the string tension µ, defined as the average string energy per unit length µ = E/L, in this system is roughly µ ≈ πf 2 a log γ where m −1 r is the size of the string core and m r corresponds to the mass scale of the degrees of freedom that UV complete the axion effective theory (in many models m r ∼ f a ). Meanwhile, γ is a numerical factor that to a first approximation is expected to be of order 1, and which is likely to be time dependent due to, for example, the number of strings per Hubble volume changing. Because of the expansion of the Universe, the majority of the DM axions produced by strings are emitted shortly before the network is destroyed, at the time of the QCD crossover. At this point log(m r /H) ≈ 70 is enormous for m r ∼ f a , changing the expected string tension significantly. Further, the large scale separation between m r and H also suppresses the coupling between strings and axions by the same logarithmic factor [34], and is expected to render emission of heavy modes associated to the theory's UV completion irrelevant.
However, the huge scale separation in the physically relevant regime presents an immediate problem in attempting to study the system using numerical simulations. To resolve the dynamics of the strings, a 2 dimensional slice of the lattice perpendicular to a string must contain at least a few grid points inside the string core, and to capture the interactions and dynamics of strings a few Hubble volumes must be simulated (we show this by analysing the systematic errors in simulations in Appendix B). Given the computational power available, the largest grids that can be simulated have N 3 ∼ 1000 3 lattice points. Consequently, the maximum scale separations that can be directly studied correspond to log (m r /H) ≈ 6. In this system, the tension of strings, and their couplings to axion and heavy degrees of freedom, are far from the physically relevant values. Indeed, even if these could be adjusted to the physical values by modifying the UV theory this would not be sufficient. For example, the properties of the string network depend on whether collapsing string loops rebound and oscillate many times before disappearing, and if strings that approach each other recombine. In order to accurately capture such dynamics, processes on all scales between Hubble and the string core size must be resolved.
Making physically relevant predictions about the system at the time of QCD crossover therefore requires that results from simulations are extrapolated over a vast difference in scale separations. What makes such an extrapolation not obviously hopeless is the possible existence of an attractor in the evolution, an understanding of which would allow a controlled extrapolation to be made. A key point of our work is that this is an extremely delicate process. In particular, we stress that a careful analysis of which features of the string network are being assumed to remain constant, or change, between log(m r /H) ≈ 6 and log(m r /H) ≈ 70 is required, and that naive extrapolations can lead to misleading results.
The calculation of the axion relic abundance by strings and domain walls has been the subject of extensive prior investigation. However, there is still substantial disagreement about the qualitative and quantitative dynamics of strings, and predictions of the resulting axion DM density differ dramatically. Previous work has been based on numerical simulations of global strings, and also theoretical analysis and numerical simulations of local strings, which are thought to reproduce the dynamics of global strings in the limit of large scale separation. To enable comparison, we postpone a discussion of the literature until we have presented our results. Instead, here we simply note that there is currently an order of magnitude disagreement about the average numbers of strings per Hubble volume at the time of the QCD crossover, which introduces a similar uncertainty on the axion number density produced; and also an uncertainty on the form of the energy spectrum of the axions produced, which can change the axion number density by almost two orders of magnitude.
Turning to the structure of this paper: In Section 2 we discuss the theory that we numerically simulate, and the technical challenges that this involves. In Section 3 we demonstrate the existence of an attractor solution, which is approached independently of the system's initial conditions, and analyse its properties. In Section 4 we study the energy released by the string network, and the spectrum of axions emitted. There we also discuss the impact on the axion relic abundance, and the process of extrapolating to the physical point in parameter space. In Section 5 we summarise our results and the possibilities for future development. Additional technical details about our simulations may be found in Appendix A, and an extensive analysis of the systematic errors is given in Appendix B. In Appendix C we present a detailed analysis of the distribution of energy into different components, in Appendix D we provide further evidence for the existence of an attractor solution, and in Appendix E we give details of how we fit the parameters of the string network. Finally, in Appendix F we analyse whether the properties of the global strings that we simulate are converging to those of local strings.

Axion Strings and Simulations
We consider a complex scalar field φ taken to have the U (1) PQ invariant Lagrangian in a spatially flat Friedmann-Robertson-Walker background. The metric is ds 2 = dt 2 − R 2 (t)dx 2 , and the Universe is assumed to expand as in radiation domination, so the scale factor R(t) ∝ t 1/2 , and the Hubble parameter H ≡Ṙ/R = 1/(2t). The potential V (φ) leads to φ getting a U (1) PQ breaking vacuum expectation value (VEV) into the radial field r(x), which has a mass m r , and the axion field a(x), which has a period 2πf a . Since we focus on the properties of the system at temperatures above the QCD crossover, the PQ breaking axion potential generated by QCD can be neglected and the axion is massless (at lower temperatures this must be added to eq. (2)). The average Hamiltonian density ρ tot = T 00 of the complex field φ is whereφ = dφ/dt, ∇ is the gradient with respect to the physical spatial coordinates R(t)x, and A ≡ lim V →∞ 1 V V d 3 x A is the spatial average of A. After decomposing φ as in eq. (3), where V (r) = m 2 r 8f 2 a r 2 (r + 2f a ) 2 . The terms on the first line of eq. (5) correspond to the kinetic and potential parts of the axion and radial modes' energies, and the term on the last line is the interaction energy between the two. In the small field limit |r| /f a → 0, the Hamiltonian can be approximated as the sum of that from decoupled axions and radial modes (i.e. by the first line of eq. (5)). However, away from this limit interaction terms between the two fields are not negligible and make the axion-radial system strongly coupled.
Note that the field's equation of motion does not depend on the ratio m r /f a directly. Indeed the dependence on the two scales f a and m r can be reabsorbed by rescaling respectively 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 ratio m r /H = 2m r t. The equations of motion in eq. (6) admit solitonic string-like solutions. As mentioned in the Introduction, these are topologically non-trivial configurations that contain loops in space where the axion field a wraps the fundamental domain [0, 2πf a ] non-trivially. The prototype of such solutions is a static, infinite, string lying along the z-axis. In cylindrical coordinates (ρ, θ, z) this is given by where g is a profile function that satisfies g(ρ) = Cρ+O(ρ 3 ) for ρ → 0 and g(ρ) The string core is defined as the region in which φ is close to the maximum of its potential, i.e. when r/f a ∼ −1, which corresponds to points at a distance less than m −1 r away from the centre of the string ρ = 0. In this part of space the axion-radial mode system is strongly coupled, and all of the terms in eq. (5) contribute to the string energy density. However, for a single string configuration the axion energy density diverges logarithmically for ρ → ∞ due to the angular gradient 1 2 |∇a| 2 . Consequently, the total string energy is dominantly in this component, and is mostly stored away from the string core. In the early Universe this leads to a string tension of the form given in eq. (1).
To analyse the dynamics of the string system in the early Universe, we numerically integrate the equations of motion, eq. (6), in 3 + 1 dimensions. Starting from suitable initial conditions, for example φ random with sufficiently large fluctuations, axion strings automatically form and evolve. In doing so, we are assuming that solutions of classical equations of motion capture the physics of strings and axion radiation. This is justified because strings are themselves intrinsically classical and the relevant part of the axion radiation has large occupation number.
The complex scalar φ is discretised on a lattice with approximately N 3 = 1250 3 grid points, and evolved in fixed steps of conformal time τ ∼ √ t. Our simulation is carried out with periodic boundary conditions, and in comoving coordinates x, so that the comoving distance between grid points remains constant and the physical distance between grid points grows ∼ √ t. Further details of the algorithms used are given in Appendix A. As the system is evolved forward, the number of Hubble lengths contained in the box side decreases ∼ 1/ √ t and the number of lattice points inside a string core also decreases ∼ 1/ √ t, as shown in Figure 1. The maximum accessible scale separation corresponds to an upper bound on the final time that can be simulated. Other possible sources of numerical uncertainty include the time step used in the simulation and the way that the contribution of the string energy is excluded from the calculation of the energy in free axions, and a full analysis is given in Appendix B.
The simulations that we carry out with m r /H ∼ 1 can be interpreted in terms of a theory with m r ≈ f a at a time when H ≈ f a , however this is not a physically relevant part of parameter space. First, for such models the scale separation simulated is far from that at the time of the QCD crossover, when the majority of axions produced by the string network are emitted. Additionally, this regime does not correspond to a system that is realisable even at early times: when H ∼ f a , the temperature of the Universe is T ∼ √ f a M Pl , where M Pl is the Planck mass, and a physical theory will be in the unbroken, PQ symmetric, phase. However, by studying the potential eq. (2) log = logሺ ሻ ≲ 6 − − Figure 1: An illustration of how the size of a string core, shaded red, and a Hubble volume, shaded blue, evolve relative to the lattice points in our simulations, where N is the number of lattice points in a spatial dimension. Requiring that the simulation contains at least a few Hubble volumes and that a string core contains at least ∼ 1 lattice point constrains the maximum scale separation that can be studied.
without including finite temperature effects, we can use the results obtained to extrapolate the properties of the string system to low temperatures, when such corrections are actually negligible. In particular, at the time of the QCD crossover T ∼ GeV and finite temperature corrections to the potential of the radial scalar are irrelevant.
Alternatively, simulations at a scale separation m r /H ∼ 1 can be interpreted as studying a model with m r f a at a time when the Hubble parameter is H f a . Axion theories with a light radial mode are less commonly considered, but in this case the simulations are directly analysing a physically realisable point in parameter space. In particular, taking m r ∼ 10 −18 GeV, the Hubble parameter is that at the time of the QCD crossover (although such a light scalar is excluded by fifth-force experiments). 4 Meanwhile, theories in an intermediate regime, with a radial degree of freedom that has a mass 10 keV m r f a , are not experimentally excluded by evolution of stars and fifth-force experiments. These correspond to scale separations at the time of the QCD crossover that are larger than can be reached with simulations, so that extrapolation is still required, but which are smaller than if m r ≈ f a .
The Lagrangian in eq. (2) includes only one heavy scalar degree of freedom, and it is clearly not the most generic that can arise in UV-complete axion theories. However, when H m r the dynamics of axion strings and radiation are expected to be largely independent of which massive degrees of freedom are included, since these only get excited when strings interact over distances of order m −1 r , e.g. when loops shrink or long strings intersect. The energy of a global string configuration remains logarithmically divergent in more complicated theories, since this comes from the axion angular gradient, which is always present. In principle one could also include interactions with other fields that φ is coupled to -e.g. SM fields -in eq. (2), but away from the string cores the couplings of axions to these are suppressed by powers of f a , and are negligible. Meanwhile the radial mode can have order 1 couplings to SM fields, however it is expected to decouple from the dynamics of the strings and axions at large scale separations. At early times, when the temperature is high, interactions of strings with the visible sector thermal plasma could modify their dynamics [49], however these effects will also be negligible at temperatures around the QCD crossover.
As well as the physical system, the literature has often used a deformed theory in which the mass of the radial mode in eq. (2) is replaced with a time dependent one where m i is the mass of the radial mode at the initial time t i . This is equivalent to a theory with a quartic coupling λ ∼ m 2 r (t)/f 2 a that decreases with time, and is known as the fat string scenario. The size of the string cores increases as m −1 r (t) ∼ t 1/2 , and the number of lattice points inside a string core remains constant throughout a simulation. The maximum scale separation that can be simulated is unchanged compared to the physical potential, however the time taken to reach a particular scale separation (starting from H ∼ m r ) is increased, and simulations can be run for longer before arriving at the upper bound. As a result, energy left over from the initial conditions is redshifted more, and there is more time available for properties of the string system to reach their asymptotic behaviour (there are additional benefits that will be seen in Section 4). Despite these advantages, we stress that by making the potential time dependent the equations of motion of the system are changed by an order 1 amount. Consequently, although the dynamics of axion strings might remain qualitatively similar to those of the physical Lagrangian this is not guaranteed, and the numerical values of the parameters of the scaling solution are not expected to be the same in the two cases. We perform simulations using both techniques and discuss their advantages and disadvantages.
In both the fat string and the physical scenarios, the axions produced in the simulations are massless, and their energy densities redshift as ∼ 1/R (t) 4 . Meanwhile, the radial modes produced are highly non-relativistic. In simulations with the physical Lagrangian their energy density redshifts as ∼ 1/R (t) 3 . In simulations of the fat string Lagrangian the scalar mass decreases with time, and the energy density of these states redshifts as ∼ 1/R (t) 4 , the same as axions.

The Scaling Solution
It has long been claimed that a system of axion strings is driven towards a particular solution, which is independent of its initial conditions [23][24][25]. Indeed, this feature is crucial for making predictions about the properties of the string system at late times, and in particular of the axion relic abundance from strings, that do not depend on the dynamics of the system at early times, which are model dependent.
The existence of such an attractor solution is simple to motivate qualitatively. Strings can lose the energy stored in their length by radiating axions and radial modes. Therefore bends in strings with curvature larger than the Hubble scale tend to straighten, and closed string loops smaller than the horizon are expected to disappear, emitting radiation. Additionally, long strings (or equivalently string loops larger than the horizon) can interact when they enter each other's horizon through a process called recombination: when strings cross they can recombine into a new configuration with a lower tension, and similarly a region of high curvature in a long string can split off forming an isolated loop. The net effect is a reduction of the total string length and the production of smaller loops, string segments with larger curvatures, and radiation. The rate at which such processes occur depends on the density of strings within the horizon. Below some critical density recombination is inefficient. In this case, the number of strings in each Hubble patch increases as the Universe expands and new strings enter the horizon. On the other hand, above a critical density recombination becomes efficient, reducing the number of strings within the horizon. As a result, the density of strings is pushed towards a particular (not necessarily time independent) value. Other statistical properties of the network, such as the distribution of the string density in loops of different length, are expected to converge similarly to a common behaviour.
We define the average number of strings per Hubble patch ξ(t) as where tot (L) is the total length stored in strings in a volume L 3 . Hence the energy density of strings is where, given eq. (1), the effective string tension µ eff (t) is expected to be with 5 µ 0 = πf 2 a . The factor m r /(H √ ξ) is anticipated to capture the main time dependence of µ eff (t) since the logarithm is cut-off by the average distance between strings (∝ t/ √ ξ). The remaining time dependence is encoded in the factor γ(t), which takes into account the non-trivial shape of the strings and is expected either to be constant or to have at most a very mild time dependence. Indeed, it will be a non-trivial check of the string network's properties that the energy density in strings extracted from the simulations is well reproduced by eqs. (10) and (11).
The existence of the scaling law (10) with constant ξ(t) = ξ 0 can easily be understood for local strings. For these the string tension is localised on the core and is constant µ eff (t) = µ. Neglecting the core size the problem only has one scale, H = 1/2t, which completely fixes the scaling law for the energy density ρ s (t) = ξ 0 µ/t 2 . The presence of a single scale suggests that during the scaling regime all the properties of the string configuration should be scale invariant.
On the other hand, for the global case several properties of the strings, including their tension and their coupling to axions [34], depend logarithmically on the core size m r , therefore logarithmic corrections to the scaling law can be expected. To account for these effects we leave an explicit time dependence both in ξ(t) and in γ(t), besides the one contained in m r /H inside µ eff (t).
In the rest of the section we will establish the existence of the (approximate) scaling solution for axion string networks, and study its properties in detail. In particular we will present results from numerical simulations demonstrating the presence of the attractor, its independence from the initial conditions, the behaviour of the parameter ξ(t), and the distribution of loops and long strings during the scaling regime.

The Attractive Solution
The existence of the attractor can be tested by studying if different statistical properties of the string network converge to the same values independently of the initial conditions of the field. As a representative example, here we focus on the evolution of the average number of strings per Hubble patch ξ(t). In Appendix D we present additional results showing that the number density, the total and the instantaneous spectrum of axions emitted from the string network also clearly converge, regardless of the initial conditions. The convergence is particularly evident in the instantaneous emission spectrum, which depends only on the string configuration at a particular moment and has no memory of earlier times.
We set the initial conditions in two different ways. In the first, we just generate sets of random fields. In the second, we construct initial conditions with a fixed number of strings by evolving random configurations until the total string length in the box reaches a required value, and then we reset the clock rescaling Hubble (more details about the procedure can be found in Appendix A). Besides allowing us to start simulations with a predetermined density of strings, the second method produces initial field configurations that have less primordial background radiation (although this would redshift away anyway) and more suitable for a cleaner study of the instantaneous spectrum. In Appendix D we show that the properties and evolution of the string network are independent of the way in which the initial conditions are set.
In Figure 2 we show the evolution of ξ with time in simulations, starting from initial conditions with different numbers of strings using the second method described above. ξ has been computed at different time shots from its definition eq. (9), using the algorithm described in Appendix A. As discussed in Section 2, it takes a longer time to get to the same value of log(m r /H) in the fat string scenario than for the physical theory. As a result the attractor regime is reached at smaller values of the log in the fat string case (here and in the rest of the paper we sometimes use the short-hand notation "log" to refer to log(m r /H)).
In both the fat string and the physical models, the convergence towards a common value of ξ is manifest. In the fat string case, the initial values of ξ span more than three orders of magnitude and, by the end of the simulations, they lead to the same value of ξ to within 10%. For the physical case the convergence is a little slower but it is still clear. In simulations starting at H = m r it is not possible to initially have more than one string per Hubble patch. As a result, to achieve initial conditions with a clear overdensity of strings we started such simulations later, when H < m r . The corresponding data in the Figure 2 have the initial points (the coloured dots) at larger values of log(m r /H). The network of global strings was first studied using field theoretic computer simulations in [37,38], and more recently over a longer time range in [50]. Ref. [51,53] showed evidence for the existence an attractor for the fat string system, respectively in two and three dimensions. Our simulations have a similar time range and constitute an independent check of the convergence of ξ in the fat string system starting from a wide range of initial conditions, and a demonstration of the attractor's existence for the physical case. Further, in Appendix D, we show that other properties of the network, including the spectrum of axions emitted, also converge.

Scaling Violation
Having shown that the attractor solution exists, we now turn to study its properties. One prominent feature is that, although different boundary conditions converge to a common value of ξ, this value does not seem to be constant in time. To see the change more clearly, in Figure 3 we show the plots of Figure 2 on a linear scale. A growth of ξ linear in log(m r /H) is evident both for the fat string and for the physical cases. The fact that simulations with an overdensity of strings first rapidly evolve to smaller values of ξ, converging to the attractor, and then start increasing again is particularly convincing. This strongly indicates that the growth is an intrinsic property of the scaling solution, rather than the sign that the attractor has not yet been reached.
The behaviour shown in Figure 3 is compatible with the asymptotic form In particular, at late times β is subleading, and the value of ξ(t) is dominated by the logarithmic term. From Figure 3 it can be seen that the coefficient α is universal, independent of the initial conditions. Indeed, the derivatives t ξ (t) = ∂ξ/∂[log(m r /H)] of the curves tend towards a common value α more rapidly than ξ itself. We can use the convergence of the slopes to a constant value to select the optimal initial conditions, i.e. those for which the scaling regime is reached the earliest. The corresponding lines are those plotted in solid black in Figure 3, and curves starting from different boundary conditions reach the same constant slope at later times. Considering only simulations that reach the scaling regime (i.e. those that show a sufficiently large region of constant slope) we extract estimates for α α fat = 0.22 (2) , α phys = 0.15 (5) .
Here the errors clearly have no statistical significance, but rather they represent our educated, conservative, guesses about the uncertainty. Plots showing the behaviour of the slope for different boundary conditions, and more details about the fit, can be found in the Appendices D and E. We do not report the values of β, since the uncertainty on these from different initial conditions is larger, and they are irrelevant for the physically interesting values of the log. As mentioned, a logarithmic violation of the scaling behaviour is not completely unexpected, since several properties of the string network, including the string tension and coupling to axions, have a similar dependence. If such behaviour is maintained at later times, as seems plausible, the average number of strings per Hubble patch will grow substantially for values of the log relevant to the QCD axion. For example, if m r ∼ f a eq. (13) would imply ξ = 10.5 ± 3.5 at the physically relevant separation, log(m r /H) ≈ 70. This value is an order of magnitude larger than that found in refs. [37,38,50] from numerical simulations on smaller grids, and the value that is usually assumed in rough estimates of the axion abundance produced by strings. Meanwhile, the extrapolated value of ξ that we obtain is only a factor of 2 larger than that recently obtained in ref. [51], which used the fat string trick and a different UV completion of the core. We stress that although the fat string system shows the same qualitative linear growth with the log as physical strings, the numerical parameters are somewhat different in the two cases. This is not surprising, and in extrapolating to the physical scale separation it is important to study the physical system, not just the fat string case.
The logarithmic enhancement in ξ was first observed and studied in the 2+1 dimensional simulations of ref. [52], but it was missed in the 3+1 dimensional ones of ref. [37,38,50], perhaps due to the use of smaller grids (which limited the time range of their simulations) and the choice of overdense initial conditions. Indeed, the combination of these two factors can produce a fake plateau at the intermediate values of the logs that were analysed (see e.g. the top curve of the right hand plot in Figure 3).
Conversely, the later simulations of the fat string case in ref. [53], made on larger grids, also observe a logarithmic increase, in agreement with our results. Finally, in the recent analysis of ref. [51], which use fat strings and a different UV completion of the core to partly include the effects of a large scale separation, the value of ξ increases with the log. The growth rate is not clear and the authors suggest that part of it may be due to spurious H/m r effects. However, the results of this reference for purely global strings are in agreement with our numerical fit, and correspond to initially overdense networks.
We will resist the temptation to interpret the log-increase of ξ(t) in terms of the reduction of the string-axion coupling or the increase of the string tension. In fact a similar growth also seems to be present for local strings, and can be seen in Figure 3 of ref. [51] and is hidden in Figure 7 of ref. [54]. 6 Neglecting gravity, the only way that local strings can maintain a scaling regime is through the production of heavy modes associated to the core scale. The latter then cannot be neglected and its presence allows for a violation of the scaling law ρ s = ξ 0 µ/t 2 argued before. As we will discuss below, the way that global strings lose energy is not so different from the one above, which might explain the similarity in the way that the scaling is violated.

Long vs Short: the scale-invariant distribution of loops
In order to further characterise the attractor solution, and to better understand its properties, we also study how the total string length per Hubble patch ξ is distributed over different loop sizes. Hℓ/π ξ ℓ ξ ∞ fat Figure 4: The fraction of the total string length ξ /ξ ∞ that is contained in loops smaller than for different time shots.
of loop length, then the quantity represents the contribution to ξ from loops of size smaller than and, in particular, ξ ∞ = ξ.
We have performed a large number of simulations using the fat string trick to acquire enough statistics to study ξ as a function of time and loop size . The initial conditions were fixed as for the thick line in Figure 3 left, so that the system started close to the scaling solution.
The results for the ratio ξ l /ξ ∞ , which gives the percentage of string length contained in loops of size smaller than l, are plotted in Figure 4, and reveal several features of the string network. Most of the string length, more than 80%, is contained in loops much larger than Hubble, of order of the full box size (in fact it seems that most of the string length is contained in a single loop wrapped around the simulation box multiple times). This leads to the abrupt increase in the right part of the plot, which saturates ξ to its asymptotic value ξ. Less than 20% of the total contribution to ξ is contained in loops of size of order H −1 and smaller, which results in the appearance of the plateau with ξ /ξ ∞ ≈ 0.2 at H 1. 7 On the left of the plot the UV cut-off corresponding to the smallest possible loops, of order the core size, is also visible. As the Universe expands the physical size of the simulation box in units of Hubble shrinks, and as a result ξ saturates its asymptotic value ξ at smaller and smaller values of H . At the same time, the value of m r /H grows so the UV cut-off moves to the left. The lines in Figure 4 corresponding to different times approximately overlap for values of sufficiently far from the UV and IR cutoffs. Therefore, since ξ ∞ grows logarithmically with time, the corresponding growth in ξ is homogeneous in . This signals that the logarithmic increase of ξ is equally distributed over all scales, and that the ratio between long and short strings stays constant in time (see also Figure 17 in Appendix B). The fact that ξ /ξ ∞ remains constant in time for 10π/H also shows that the number of loops of a particular length per Hubble patch does not change, apart from this logarithmic increase. As loops shrink and disappear (or recombine with other strings) they are replaced at the same rate by larger loops themselves shrinking, or by new loops being produced from interactions of long strings, which is an indication that the attractor solution has been reached. When, at the final times, the Hubble scale becomes of order of the box size there is no longer a sharp distinction between long and short strings.
Another feature of the loop distribution is evident from the plot: for loop lengths smaller than Hubble ξ ∝ (the dashed line in the plot), so dξ /d ∝ dn /d = const. This means that the number of loops in each logarithmic interval of length is constant, over almost two orders of magnitude. Equivalently, the 10% of the full string length contained in sub-Hubble loops is equally distributed over all loop sizes on a linear scale. This approximate power law seems to become a better fit as time progresses, suggesting that it is an intrinsic property of the attractor solution, and further confirming that the attractor regime has been reached within the time range of the simulation. 8 Our analysis suggests that, in the infinite volume limit, the distribution of string length in the attractor solution is of the form depicted in Figure 5. Roughly 80% of the string length per Hubble patch is contained in long strings (infinite string loops), while the remaining 20% is distributed in loops ranging from the core to the Hubble size, with equal numbers of loops in each decade of length. The total string length grows logarithmically according to eq. (12) with the relative 4:1 ratio fixed.
In the literature the parameter ξ is sometimes defined restricting to long strings only. However, the loop distribution that we observe implies that the two definitions only differ by 20%, and more importantly the factor of proportionality is constant in time.
Since we do not directly use the quantitative behaviour of ξ in our subsequent estimates of the final axion abundance, we have only performed this analysis in the fat string case. The similarity in the behaviour of all of the other properties studied suggests that the picture for the physical case should be qualitatively similar.

The Spectrum
We now discuss how the string network radiates energy, and study in detail how the energy in the system is split among the different components (strings, axions, radial modes), and the way that the energy in axions is distributed in modes of different frequencies. We then analyse the evolution of the axion number density, and the extrapolations of this to the physically relevant point. However, before turning to the results of our simulations, we first establish the physically relevant quantities and describe how these affect the axion relic abundance.
The conservation of energy implies that, in order to maintain the scaling regime, the string network must constantly lose energy into radiation. This is because, in the absence of interactions, the number of strings per Hubble patch would increase fast as more strings re-enter the horizon. To keep ξ approximately constant, the excess string length must be destroyed, emitting energy.
The rate at which the system of interacting strings releases energy can be calculated by comparing the energy density in the scaling regime (parameterised by eq. (10)), to that of a "free" network of strings. By free we mean that long strings remain essentially at fixed comoving coordinates, so that ξ(t) ∝ R 2 (t) ∝ t. The energy density of such a system is where µ(t) = πf 2 a log(m r d(t)) is the tension of free strings and d(t) ∝ 1/R(t) ∝ 1/ √ t parametrises the average distance between strings.
We consider a network of free strings that has the same string configuration as the interacting system with energy ρ s (t) given by eq. (10) at a time t 0 . The energy of such a system is with µ(t 0 ) = µ eff (t 0 ), so that this matches that of the interacting network at t = t 0 . The string energy density of the interacting system in the scaling regime decreases faster than ρ free s , and the difference corresponds to energy that is released. The rate Γ at which the interacting network emits energy into radiation is therefore given by The last limit holds at late times when the log is sufficiently large, as is the case in the physically relevant regime. 9 The rate of energy loss from strings can be split into two contributions Γ = Γ a + Γ r , corresponding to the rate of energy transfer to axions and radial modes respectively. The corresponding continuity equations for the axion and radial energy densities are theṅ where z is a factor ranging from 3 (for non-relativistic radial modes) to 4 (for relativistic radial modes, and in the fat string case), and the dots represent subleading contributions from energy transfer via axion-radial mode interactions. At sufficiently late times in the scaling system's evolution, most of the energy released by the string network is expected to go into axions, since radial modes are heavy relative to Hubble and harder to excite. The energy density in axions is therefore mostly fixed by conservation of energy. In contrast, the number of axions produced depends on the energy spectrum with which they are emitted by the string network. 10 Indeed it is the nature of the axion spectrum that is the source of the largest uncertainty in the relic abundance of axions produced from strings, and this has been the subject of disagreement for many years [21,35,39,49]. Before reporting the arguments underlying the different possibilities, we first review how the axion number density depends on the properties of the spectrum.
Since strings typically have curvature of order Hubble, the natural expectation -always assumed in the literature, but never confirmed in simulations -is that the spectrum of axions emitted at each instant is peaked at momenta of order the Hubble parameter at that time. Meanwhile production of modes with momentum below Hubble or above the string core scale is expected to be strongly suppressed. Between these scales an approximate power law is usually assumed, which determines the hardness of the spectrum. If the spectrum is soft, meaning that it is sharply peaked in the IR scale (around Hubble in this case), a relatively large number of axions will be released to account for the total energy lost by strings. If the spectrum is harder, with a larger UV tail, fewer axions will be produced, although each will be more energetic. The expectation that the attractor solution is approximately scale invariant corresponds to a prediction that the location of the spectrum's peak relative to Hubble, and the power law fall off, are constant up to possible logarithmic corrections.
From eq. (18), if we neglect the energy emitted in radial modes (and axion-radial interactions), which as we will see is a small fraction, we have that and therefore the energy density in axions at a time t, when they are still massless is that is, integral of the energy emitted at each previous instant, appropriately redshifted (and we omitted the initial time in the integral since it is dominated by late times).
We also introduce the differential energy transfer rate which depends only on the axion momenta k, the time, and the core size m r . It is convenient to further split this up as where the function F (x, y) fully characterises the shape of the spectrum (through the variable x), and its time dependence (through the variable y). Combining eqs. (20), (21) and (22) we get a formula for the axion spectral energy density where primed quantities are computed at the time t , the redshifted momentum is defined as k = kR/R , and we have left a possible time dependence in the core mass scale m r to include the fat string case. Eq. (23) is just the time integral of the instantaneous spectra appropriately redshifted, and the change in power of the redshift factors compared to eq. (20) is due to the extra power of k in the differential spectrum. The total number density of axions is therefore given by To see how the number density depends on the shape of the spectrum we consider an analytic form that reproduces the theoretical expectation: i.e. a single power law 1/k q with an IR cutoff at k = x 0 H and a UV one at k = yH ≈ m r (the extra factors are required to have the right normalisation). Substituting this into eq. (24) and taking the large time limit, which corresponds to keeping the leading log contributions, we get Given the large size of the log, the last factor strongly depends on whether the power q is larger, equal, or less than 1 (see Figure 6). In fact we can rewrite the expression above as in which the last factor varies from O(1) for q > 1 to O(10 −2 ) for q = 1, and is exponentially small (in terms of the log) for q < 1.
The dependence of the axion number density on q can be easily understood in terms of our previous qualitative discussion. For q > 1 the spectrum is soft, most of the energy is emitted with momenta of order Hubble, and the final number density is of order the total energy density contained in strings, H 2 ξµ eff , divided by the average axion momentum O(H). For q = 1 energy is equally distributed in logarithmic intervals of momentum. Therefore, although most of the axions are still emitted with momenta of order H, the total number of axions emitted is smaller by a factor of the log. For q < 1 the spectrum is UV dominated, and the majority of the energy is distributed to axions with large momentum so that the axion number density is power suppressed by the UV scale.
The different behaviour of the number density for different choices of q can be linked to the change in the average momentum of the axions in the spectrum. If we define the inverse average momentum as the number and energy densities are related via n a = k −1 ρ a . Depending on whether q is larger or smaller than unity the average momentum is parametrically of order H or m r respectively. The huge ratio of scales m r /H in the physically relevant part of parameter space results in an enormous range of possible values of n a . It is therefore clear that understanding the spectrum is of paramount importance if results obtained at the values of the log accessible in the simulations are to be extrapolated to the physical values. In particular, even a small change in the behaviour of the spectrum could change the extrapolated value of the relic abundance by many orders of magnitude.
We can now identify the main source of disagreement in the literature. Refs. [21,22,35,36] claim that at late times, when the scale separation is large, the coupling of strings to axions is small and the rate that axions are emitted is suppressed, and as a result the dynamics of axion strings are close to those of local strings. If this is the case, the expectation based on the Nambu-Goto effective theory is that loops will oscillate many time before emitting their energy, producing a spectrum that is sharply peaked at small frequencies, of order their initial size, i.e. Hubble. Consequently, they predict that q > 1, and that the number density of axions produced by strings in the scaling regime will dominate over the contribution from misalignment, here taken as a reference value n mis a = θ 2 0 Hf 2 a , with θ 2 0 ≈ 5. 11 Setting µ eff = πf 2 a log following the theoretical expectation in eq. (11), log ≈ 70 and x 0 ≈ 10 we get n q>1 a n mis which ranges from 30 to 300 depending if ξ is taken to be close to 1 or 10. Conversely, refs. [49,56] claim that string loops do not oscillate, but instead efficiently shrink emitting all of their energy at once and producing a spectrum with q = 1. In this case the number density from strings is suppressed n q=1 a n mis and can even be subleading with respect to that from misalignment if ξ is taken to be 1. 11 In this expression for n mis a , H is the Hubble parameter when the axion mass becomes cosmologically relevant, which is around the time of the QCD crossover although the exact value depends on f a .
In the rest of this section we present a detailed analysis of the spectrum obtained from simulations with the aim of understanding which of these possibilities is more likely. We therefore postpone further discussion to the end of the section, when we compare our findings with the existing literature.
To analyse the spectrum emitted by the scaling solution we fixed the initial conditions in simulations to be as close to the attractor as possible, corresponding to the black curves in Figure 3. This isolates the radiation emitted in the scaling regime as much as possible, and reduces contamination from pre-existing radiation. In Appendix D we show that, starting with different initial conditions, the spectrum and number density converge to those of the scaling solution, so that the results we obtain do not depend on this convenient choice.

Energy Budget
In analysing the distribution of energy in the scaling solution, and the rate at which axions are produced by strings, it is useful to organise the total energy density stored in the complex scalar field into three components, namely Here ρ tot = T 00 is the total energy density of the scalar field as given by the average Hamiltonian density eq. (4); ρ s is the contribution contained in strings; ρ a is the energy density in axion particles; and ρ r is that in radial modes. At early times axions, strings, and radial modes are strongly coupled to each other so this separation is ill-defined, but at later times the individual components decouple and the separation becomes meaningful. In the scaling regime, the theoretical expectation is that when the energy density in strings is parameterised as in eq. (11), γ(t) will be of order 1 and vary only slowly with time.
The way that we actually compute the various components in eq. (31) is as follows: The axion energy density is calculated from the spatial average ofȧ 2 away from the core of the strings (close to the string cores, the motion of strings gives a significant contribution toȧ 2 ). By using 2 1 2ȧ 2 rather than 1 2ȧ 2 + 1 2 (∇a) 2 we avoid the part of the energy that corresponds to the strings' tension, which is mostly contained in (∇a) 2 . We have checked that our results for ρ a redshift as expected (i.e. as relativistic matter) and that they are robust against different types of string-core masking (more details about our procedure for screening strings and the consistency tests can be found in Appendix B). The radial energy density is computed by averaging the part of the Hamiltonian density that involves only the radial mode, i.e. 1 2ṙ 2 + 1 2 (∇r) 2 + V (r) , again away from the string cores. Finally ρ s is simply extracted from the difference ρ s = ρ tot − ρ a − ρ r , which avoids doublecounting energy contributions in the other components. The string energy density defined in this way includes the energy density stored in the axion-radial interactions, corresponding to the terms (r/f a +r 2 /2f 2 a )(∂a) 2 in the Hamiltonian (second line of eq. (5)), only part of which (that in regions close to string cores) genuinely contributes to the string energy. The remainder corresponds to interaction energy between axion and radial modes. Such interactions could in principle trigger parametric resonance and have a substantial effect, however we have checked that they only give negligible oscillating corrections to the energy densities in our simulations. 12 We can compare ρ s extracted in this way to the prediction obtained using the theoretical expectation for the string tension, based on the typical separation between strings, and the measured values of ξ(t). In particular, we compute the effective tension µ eff = ρ s (t)t 2 /ξ(t) from the definition eq. (10), using ξ(t) and ρ s (t) from the simulation. This can be compared to the theoretically expected form which is obtained by replacing γ(t) in eq. (11), by a constant γ c = 1/ √ 4π as a reference (we choose 1/ √ 4π somewhat arbitrarily based on the average distance between strings if they were all parallel, but any roughly similar value would also be theoretically reasonable).
In Figure 7, we plot the ratio µ eff /µ th as a function of time, for the fat string and the physical cases. The closeness of µ th and µ eff over the entire time range that the system is in the scaling regime is highly non-trivial. These quantities could have differed by orders of magnitude, or had different time dependences, but instead the naive theoretical prediction reproduces the results from the simulation to within 20% throughout. As well as showing that our method of computing ρ s is meaningful, it is a strong sign that eq. (10) with µ eff replaced by the theoretical expectation eq. (32) correctly captures the dynamics of the string system. This includes the logarithmic growth of both ξ (t) and the string tension due to the increasing scale separation, as well as the variation of ρ s compared to ρ a and ρ r with time.
The small difference between µ eff and µ th is not worrisome. First of all µ th does not include relativistic effects. These are expected to be small on the basis that the energy density is dominated by long strings, whose motion is damped by the Hubble expansion; indeed such effects are expected to increase the ratio µ eff /µ th which is instead smaller than (and close to) unity. Second, we do not have a reliable way to compute γ analytically since its value is determined by the loop distribution and the shape of the strings, which can also depend on time. Further, the parametrisation µ = µ 0 log (m r γ/H √ ξ) with an IR cutoff ∼ m r /H √ ξ and constant γ applies only to long strings (for loops with radius smaller than Hubble, which make up less than 10% of ξ, the tension is expected to be cutoff at smaller values). A simple modification of γ can fix the ratio µ eff /µ th = 1 at all times, however for the moment we are not interested in such a detailed understanding of µ eff , and we content ourselves with the degree of agreement obtained in Figure 7.
Turning to consider the energy in axions and radial modes, in Figure 8 we show the proportion of the total energy that is in the three components of eq. (31) as a function of time, in the fat string and the physical cases. We plot this only for times corresponding to the range log(m r /H) = 3.5÷6 for the fat string case, and log(m r /H) = 4 ÷ 6 for the physical one. Data at later time may not be safe from finite volume effects (see Appendix B) while data at earlier times is not representative of the scaling regime since there is not yet a sufficient hierarchy between Hubble and the core scale to consider these decoupled (as will be clear when we analyse the axion spectrum, shortly).
First, we note that in the fat string case all three components redshift as 1/t 2 , up to logarithmic  Figure 7: The closeness of µ eff = t 2 ρ s (t)/ (ξ(t)) to the theoretical prediction µ th , defined in eq. (32), (plotted in terms of the ratio of these two quantities) is a non-trivial check that our procedure to extract the string energy density ρ s is reliable. More importantly, it shows that the relation eq. (10) can be used to predict the energy in the string network for a given string density and time, by replacing µ eff by µ th , and that the theoretically predicted logarithmic growth in the string tension is seen in simulations.
corrections. Therefore the time dependence in Figure 8 is only a result of the increase in energy in strings (due to the factor ξ(t)µ eff (t) growing) and energy transfer from strings to axions and radial modes. Indeed, using eqs. (17) and (20) these two effects combined predict that ρ a /ρ s ∝ log(m r /H), which explains the relative change in ρ a and ρ s . The proportion of energy in radial modes appears to be approximately constant in time, accounting for around 13% of the energy budget.
The situation is less clear for the physical case. The radial mode now redshifts as nonrelativistic matter 1/t 3/2 , slower than the other components. If we use again eqs. (18) and (17), and assume constant energy transfer rates Γ i , then ρ r /ρ a ∝ t 1/2 /(ξ(t)µ eff (t) log(m r /H)). Therefore, if the rates remain constant until late times, energy in radial modes will eventually dominate that in axions. However, over the time range plotted in Figure 8 this expression predicts that ρ r /ρ a remains approximately constant, which matches our results from simulations. The fact that the string energy density is not decreasing with respect to that in axions, as happens in the fat string case, may be an indication that larger scale separations are required for the asymptotic behaviour Γ ∝ ξ(t)µ eff (t)/t 3 in eq. (17) to be reached in the physical case, or due to the shorter time that simulations can be run compared to the fat string scenario.
Perhaps more revealing are the rates at which energy is transferred into axions and radial modes. In Figure 9 we show the time dependence of the ratio r a = Γ a /Γ in the simulations, which we compute from eq. (18) by taking derivatives of the energy densities calculated as above (setting z = 3 in the physical case, z = 4 in the fat string one, and neglecting the dots, i.e. including possible energy transfer between axions and radial modes in Γ a ). It seems that in both cases the fraction of energy that is transferred to axions is roughly constant, and the value in the fat string scenario is compatible with the approximately constant proportion of the total energy that is in radial modes, plotted in Figure 8. As discussed, in the physical case a constant energy transfer will eventually lead to radial modes dominating the total energy density, but is compatible with the results in Figure 8 for the time range that we can simulate.
Similarly to the fraction of the total energy in strings ρ s /ρ tot , the constant value of r a in the physical case might only be a transient effect. In fact we expect that the rate of energy transfer from radial modes to axions will increase if the abundance of radial modes grows, so that at sufficiently late times the dominant net effect of energy loss from the string network should be the production of axions. Because of this, any interpretation of the fractional rate r a computed from simulations in the physical case should be taken with a grain of salt.
To summarise, in both the fat string and the physical cases we find that most of the energy radiated by strings in the scaling regime goes into axions, and these give the largest contribution to the energy density at late times. Meanwhile, the string energy density is well reproduced by the ansatz from the scaling solution, with the theoretically expected string tension and the values of ξ(t) measured in simulations. There is a non-negligible component of the energy density in radial modes, and an approximately constant proportion of the energy emitted by strings goes into such states. This surprising result clashes with the expectation that heavy modes decouple from the evolution of macroscopic soft objects, and that their production is suppressed. We will see that this phenomenon has a close analogue in the rate at which high momentum axions are emitted.

Axion Spectrum
We finally arrive at our analysis of the spectrum of axions emitted, the shape of which has a dramatic effect on the axion relic density, as previously discussed.  Figure 9: The fraction of the instantaneous energy density emitted by strings that is converted into axions (as opposed to radial modes) r a = Γ a /Γ as a function of time, for the fat string and the physical simulations. At late times more than 70% (80%) of the energy released from the string network in the scaling regime goes into axions in the physical (fat string) case. The rest of the energy goes into the production of massive radial mode states.
In Figure 10 we show the differential energy density ∂ρ a /∂k at different time shots, for both the fat string and the physical cases. The spectrum is such that dk ∂ρ a /∂k = ρ a , and it has been computed from the Fourier transform ofȧ, with the strings appropriately masked out, as explained in Appendix C. We plot the spectrum as a function of k/H to highlight that it remains peaked around momenta of order Hubble as this decreases, and we divide it by the factor Hf 2 a to remove the main time dependence (see eq. (20)).
The spectra have a number of relevant features. We start with the fat string case, for which the core scale m r evolves in exactly the same way as physical momenta redshift. As a result, radiation produced by the cores with k ∼ m r remains at k ∼ m r at all later times, and contributes to the same part of the differential spectrum regardless of when it was originally produced. Similarly, since the simulations begin with H ∼ m r , the typical momentum of pre-existing radiation from the initial conditions is k ∼ m r and states from these early times will subsequently remain at ∼ m r . Therefore the part of the spectrum sufficiently below m r is entirely from genuine radiation produced by strings, and is independent of the physics of their cores. Figure 10 clearly shows that the spectrum is peaked in the IR at around the Hubble scale, in particular at k/H ∼ 5÷10 for all values of the log larger than 3. There is another small peak in the UV, which is exactly at k = m r /2, the typical frequency of parametric resonance. These modes could be produced by the string core itself, or by the conversion of non-relativistic radial modes into axions. The position of this peak at a particular time identifies the part of the spectrum that is sensitive to core scale (i.e. UV) physics. For times corresponding to log(m r /H) 3 the peaks in the UV and the IR cannot be distinguished, and it is only for log(m r /H) 3.5 that there starts to be enough of a hierarchy to justify considering the IR dynamics as decoupled from the core scale physics. For this reason we do not consider smaller values of the log in our analyses  Figure 10: The axion energy density spectrum ∂ρ a /∂k as a function of the ratio of the physical momentum k to the Hubble parameter. The results are shown for different times (i.e. different values of log(m r /H)) for a string system in the scaling regime, in the fat string (left) and the physical (right) scenarios. In both cases the spectrum is dominated by a broad peak at around k/H = 10, and emission at lower momenta is suppressed. As a reference, for each time shot we also show the value of k corresponding to the core (k = m r /2), above which there is very little emission, as expected.
involving the spectrum.
At late times the size of the simulated box, L(t), starts cutting off the low momentum part of the spectrum, as can be seen from the interruption of the spectrum at the minimum non-trivial k = 2π/L(t) in the final shots. Subsequently this would start altering the IR peak, so we only analyse the spectrum until log(m r /H) = 6 when this effect is still harmless (a more detailed study of the various systematics can be found in Appendix B).
As expected, the spectrum is power suppressed for momenta smaller than the IR peak or larger than the UV peak -long wavelength modes are inhibited by the horizon, while the high energy ones are suppressed by decoupling. In particular, the spectrum seems to fall as ∼ k 3 in the far IR, and as ∼ 1/k 2 in the far UV. The region of physical interest is in between the two peaks. The spectrum reaches a stable form towards the end of the simulation, and in the last Hubble e-folding (i.e. for logs in the interval 5 to 6) its shape remains very similar, modulo the shift of the UV peak. However, even though ∂ρ a /∂k is largest in the IR, its area is dominantly in the UV, since the slope between the two peaks is less steep than 1/k. This can be seen more clearly by plotting ∂ρ a /∂ log k on a log scale (shown in Figure 27 in Appendix E). We conclude that although most of the axions produced by the evolution of the string network at late times are soft, the majority of the energy density is contained in UV axions with energy of order the inverse core size.
The situation for the physical case seems similar, although the uncertainties are larger. In fact it is in studying the spectrum that the advantages of the fat string trick are most pronounced. In the physical case, axions produced at early times with momentum of order the core-scale redshift with respect to m r , which is constant. This means that, although the distance between the IR peak and m r /2 at late times is the same as in the fat string case (for the same value of log(m r /H)), the spectrum is contaminated by core-scale radiation produced earlier and redshifted down to k ∼ Hm r /2. This effectively halves the region of the spectrum that is free from UV dependent contributions. Indeed the spectrum shows that the IR and the UV dynamics are not decoupled before logs of order 4.5 or 5, and for this reason we do not consider quantities that rely on the spectrum at times corresponding to log < 4.5. Similarly, the UV peak is now replaced by a broader feature, caused by the convolution of core-scale radiation emitted at different times and redshifted by different amounts.
Because of these disadvantages, extracting the shape of the spectrum between the peaks is more challenging in the physical case. The rest of the spectrum shows similar features to the fat string case, with a strong suppression of modes below k ∼ (5 ÷ 10) H and above k = m r /2 (for the UV modes, this time with a stronger suppression, ∼ 1/k 3 instead of ∼ 1/k 2 ). As for the fat string case, the IR peak seems to dominate the spectrum, especially at late times, while the area is still dominated by the UV region. Therefore, most of the axions produced by the evolution of the string network at late times are again soft, but most of the energy density is in UV axions. However, we stress that the results in the physical case are less clear and should be interpreted with caution.
To the best of our knowledge, the only other serious attempt to extract the axion spectrum in the scaling regime was in ref. [50], based on results of simulations carried out on a somewhat smaller grid than ours. In that paper the authors observe an exponential suppression of the spectrum at large momenta, which they interpreted as indicating a strongly IR peaked distribution. However, the range of momenta that shows such behaviour seems to lie above the scale of the core, m r /2. Unfortunately, the region at smaller momenta has been very poorly binned, so that little information on the actual behaviour of the spectrum in this region is available.

Instantaneous Emission
In addition to the overall spectrum, the shape of the instantaneous axion spectrum, i.e. the function F (x, y) in eq. (22), is crucial for understanding the properties of the emitted axions, and in particular for inferring the evolution of the axion number density at later times.
We compute F from the spectrum by inverting eq. (23), namely where the factor A = H/Γ is fixed by requiring that F is normalised to 1. To evaluate the derivative numerically we took the difference of R 3 ∂ρ a /∂k between two subsequent time shots (separated by ∆log = 1/4). The results are shown in Figure 11. Since interactions with radial modes induce small oscillations in time with frequency ∼ m r of the axion energy density (see Appendix C.1) the procedure to extract F is subject to fluctuations at frequencies near the core, as evident in the plots (as explained in Appendix C.1 this effect is more pronounced for physical strings than for fat ones).
For the fat string case we plot F at three time shots, i.e. at three values of the log. At all of these times F has an IR peak at momentum around the Hubble scale (and suppressed emission at lower momenta), and a UV peak at the scale of the string cores. In particular, the IR peak is located around k/H ≈ 5 and extends to k/H ≈ 10, while the UV peak is at k ∼ m r /2. As a result, the ratio between the maximum and the minimum momenta between which the instantaneous spectrum can show a power law behaviour is bounded by k max /k min m r /40H N/120, where N is the number of points in the spatial grid of the simulations (assuming a lattice spacing a = m −1 r and that the simulations are stopped when the box size L = 3/H, and also requiring k/H > 20 to be safely away from the IR peak). For our simulations this means that k max /k min ∼ 10.
Although the separation of the UV and IR peaks is not yet large for log = 4, at values above 5 an intermediate momentum range can be clearly identified in which the differential spectrum shows a definite power law behaviour. In this region the instantaneous spectrum of axion emitted is compatible with a behaviour 1/k q with q ≈ 0.7 ÷ 0.8. More details about how we extract q, and further plots, can be found in Appendix E. Consistent with our analysis of the convoluted spectrum, this value means that most of the energy released by the string network goes into high energy axions, although most of the axions are soft.
We also note the constant form of F at different times (over a range in which H changes by 2 orders of magnitude) is another demonstration of the scaling behaviour of the system. In particular, the IR and UV peaks remain their expected positions, and the intermediate power law with q ≈ 0.75(5) is constant with time to within the uncertainty. Although there is no noticeable change in q in the interval of logs between 5 and 6, the uncertainty is relatively large. Consequently, we cannot exclude the possibility that, as for ξ, a logarithmic scaling violation occurs. This could lead to, for example, a behaviour q ∼ 0.7 + log(m r /H) for some small constant . In Appendix D we also show that F is independent of the string systems initial conditions, confirming that we are indeed analysing the properties of the attractor solution.
As expected, the analysis of the physical case is more difficult. UV modes pollute the convoluted spectrum well below the scale m r , so useful results only start appearing at late times. Because of this, in Figure 11 we only show F for the last time shot at log = 6. Although there are still large fluctuations due to higher frequency modes, an approximate power law can be recognised at this time, with a value of q that is again smaller than 1 and which appears to be similar to in the fat string scenario.
The form of the instantaneous emission found above clashes with either of the theoretical expectations discussed in eqs. (29) and (30), which predicted q ≥ 1. Similarly to the radial modes, axions with momentum of order the string core scale have not decoupled from the dynamics of strings at the final times log(m r /H) 6, despite the relatively large hierarchy m r /H ∼ 500. To investigate this surprising result further, in Appendix F we study the dynamics of the collapse of a single axion string loop. In particular we analyse whether this is converging to the prediction for a Nambu-Goto string in the limit that its initial radius is large compared to its core size. We find that the radius of the axion string loop as a function of time during its initial collapse does indeed get closer to the cosine prediction (and is extremely close to the prediction based on the effective theory of global strings coupled to radiation in the limit of large scale separation [34]). At the scale separations that can be simulated (log(m r /R 0 ) 5 where R 0 is the initial loop radius) the loop does not rebound significantly, as opposed to the case where log(m r R 0 ) 1 in which the loop is expected to rebound many times according to arguments based on the use of the Nambu-Goto string action [31]. On the other hand, whether the loop bounces also depends on its dynamics at small radius, a regime in which the Nambu-Goto approximation is not valid. Since it is at such times that energetic modes are more efficiently produced, this result is consistent with our observation that the spectrum also deviates from the expectation based on a loop oscillating many times.

Number density
We are now ready to consider the axion number density. We evaluate this in simulations using the differential spectrum ∂ρ a /∂k at different time shots, and eq. (24). The results for n a , normalised with the factor Hf 2 a , as a function of time are shown in Figure 12 for the fat string and the physical scenarios. In both cases n a /H increases logarithmically over the time range plotted, although the values and the slopes are different. Therefore, similarly to the fit of ξ(t), our first conclusion is that at the quantitative level the fat string system differs from the physical theory, although they seem to have the same qualitative features.
As we have stressed, the values of the log that can be analysed with numerical simulations are nowhere near large enough to reach the physically relevant region of parameter space. Consequently, n a must be extrapolated to extract predictions, and the way that this is carried out has a dramatic impact on the results obtained. Such an extrapolation is a viable possibility because we have shown the existence of an attractor solution, and we have found that the energy of the string network is accurately reproduced by eq. (10) with the theoretical prediction for µ eff , eq. (32). Additionally, given our results for ξ(t) it is plausible that the fit in eqs. (12) and (13) can be extended all the way to the physical scale separation. The remaining component needed in order to predict the axion relic abundance is an understanding of the form of the instantaneous emission spectrum at late times.
As usual we will discuss the fat string case before moving to the physical one. From Figure 12  it is tempting to extrapolate n a / (Hf 2 a ) linearly in the log, however this procedure is too naive. From eq. (27), a linear behaviour is expected only for very large logs and even then only if the power law of the instantaneous emission spectrum is q = 1 (and ξ(t) also increases linearly with the log). This is not the case in simulations since, as we saw in the previous section, q ∼ 0.75 and the logs are not large enough to neglect 1/ log corrections in eqs. (17) and (24). We conclude that the linear behaviour observed is most likely a transient effect. Indeed, should the power q remain constant below 1 the number density would start decreasing exponentially with the log at late times (from eq. (27)). Meanwhile, if q grows to be significantly larger than 1 the number density will increase as the log squared.
The different behaviours at large values of the log are shown in Figure 13, assuming that ξ continues to increase logarithmically. It is clear that a naive linear extrapolation leads to significantly different results compared to if q is assumed to remain constant at 0.75. Since we are not able to reliably study the behaviour of q with time, we cannot exclude the possibility that a small scaling violation results in q growing to 1 or even larger values, and we also show the axion number densities at large values of the log in these cases. It can be seen that if such a change does occur, the large log behaviour of n a is completely different to both the linear extrapolation and the q = 0.75 possibilities, leading to much larger number densities. The extrapolations in Figure 13 have been obtained from eq. (24) using a form for F fitted from the results shown in Figure 11 with q modified between the IR and UV peaks, although the results are not very sensitive to the exact shape of F away from the power law region (more details on the form of F used are given in Appendix E).
For the physical case there is even more uncertainty in the value of q, and consequently on the late time number density. In Figure 13  1 it grows with log squared (again assuming that ξ continues to increase logarithmically), while if q < 1 it is exponentially suppressed in terms of the log.
It follows that if the spectrum remains UV dominated with small q (i.e. with q 0.85) for large values of the log, the number of axions produced from strings during the scaling regime is negligible. However, if q increases at larger log, the axions produced by strings can easily match or dominate those from misalignment. Indeed, even a very mild log dependence of the power q, for example of the form q = 0.7 + 0.01 log(m r /H), too small to be excluded by our simulations, would be enough to make the final abundance at log = 70 more than 2 orders of magnitude larger than the misalignment contribution.
As discussed in Section 2, the results of simulations only depend on the ratio m r /H, and can be interpreted either in terms of a theory with m r ∼ f a at an early time when H ∼ f a , or a theory with a much smaller m r f a at a later time H f a . Fixing the Hubble scale to its value when the axion mass becomes relevant H ≈ m a ∼ Λ 2 QCD /M Pl , the values of n a calculated in our simulations can therefore be used to extract the physically relevant axion number density for theories with extremely small values of m r (as would occur in a UV completion with a complex scalar field that had a tiny quartic coupling). However, constraints from observations of the evolution of stars and fifth force experiments require that m r keV in viable theories, corresponding to rather large values of the log (∼ 23 for the reference value of f a = 10 11 GeV), which is still beyond the reach of simulations. Due to this phenomenological requirement, we only begin the extrapolations in Figure 13 at this scale separation. Meanwhile, the typical axion models with m r ∼ f a correspond to logs ≈ 70 at the physically relevant time as usual, at the right hand boundary of the plot.
The main message from this section is that an understanding of the late time evolution of the properties of the string network is of paramount importance for a correct extrapolation of the axion abundance to the physically relevant parameter range. If only small-log data are available, very precise computations of the spectrum and its time dependence will be necessary, or if theoretical arguments that q → ∞ in the limit of large scale separation are believed, a prediction for the relic abundance can already be made. In contrast, a naive extrapolation of the axion number density (such as the linear-log one plotted above) could easily give results that are off by several orders of magnitude.
Our final results for the axion number density calculated in simulations are not in disagreement with those obtained by other groups. Despite this, the conclusions that have been reached in the previous literature about the abundance at the physically relevant scale separations are incompatible with each other, primarily due to the different ansatz used in extrapolating. In particular: • Refs. [21,22,35,36] assume that q > 1 based on arguments related to the evolution and emissivity of effective Nambu-Goto strings, which are expected to reproduce the dynamics of global strings at large logs. Assuming ξ = O(10) they conclude that strings should produce a large number of axions. Our result would only be compatible with this estimate if the power q of the axion spectrum increases at large scale separation, however we are not aware of any reliable evidence that this is indeed the case beyond the Nambu-Goto approximation. In fact, as mentioned, energetic modes are expected to be efficiently produced from collapsing loops when their radius is comparable to the string thickness (or by long strings that are similarly close together), and in this regime the Nambu-Goto effective description is not valid.
• Refs. [49,56] on the contrary assume that q = 1, supported by computations of the axion spectrum produced by collapsing loops at values of the log 5. Using ξ = O(1) they conclude that the abundance of axions produced by strings is not larger than that from misalignment. While the value of ξ assumed is compatible with those measured in our simulations, the time evolution of ξ that we observe suggests that larger values need to be used in the physical regime. Additionally, the study of the spectrum from collapsing loops is not sufficient to infer the spectrum emitted by the full string network in the scaling regime. Instead, this is the result of a combination of spectra produced by loops and long strings, and therefore also depends strongly on the distribution and the evolution of the latter. While the value of q observed in our simulation is not far from 1, we have no evidence to support an expectation of an asymptotic value of q = 1 at large logs.
• Refs. [50] performed similar simulations of global strings to those that we have carried out. They find a value of ξ ≈ 1, in the same approximate range as that obtained in our simulations, but they did not observe any linear increase with the log (probably for the reasons discussed in Section 3.2). For the reasons discussed in Section 4.2, they inferred a very IR dominated spectrum and conclude that strings would efficiently produce a large axion number density, even though their ξ is lower than our extrapolated value. This example shows how a different interpretation of the spectrum can lead to a completely different conclusion, despite the fact that the spectrum measured in the simulation for small logs might be qualitatively in agreement. At small logs, a very detailed study of the spectrum is mandatory to reliably extrapolate the number density to the physical point.
• Finally in refs. [51,57] a different UV completion of the string cores was introduced in order to simulate fat strings with a large effective tension µ eff ∼ 70µ 0 , despite the actual core size remaining bounded by the lattice spacing size log(m r /H) ∼ 6. While this trick is claimed to remove the need for a large extrapolation connected to the string tension and the effective axion-string coupling, the one associated to the decoupling of the radial modes of the string core remains. The results in ref. [51,57] confirm the non-trivial dependence of ξ on the logs, although the actual functional dependence is not clear from the analysis. The authors did not present explicit results for the spectrum. The final number density of axions observed, which also includes the contribution from the domain walls and the destruction of the string network owing to the axion mass, turns out to be smaller than that from misalignment. This result seems to be compatible with an underproduction of axions from strings, associated to a spectrum with q < 1, in agreement with our results if we assume that q remains below 1 at larger logs.
However, note that in this case the spectrum (and the way that strings sustain the scaling regime) is dominated by core scale physics. This makes the extrapolation of the string core size from log(m r /H) ∼ 6 to log(m r /H) ∼ 70, as subtle as in the normal case. In fact we cannot exclude the possibility that at larger logs (i.e. relatively thinner strings) the production of energetic modes gets suppressed and q increases. The extrapolation of log(m r /H) to the physically relevant values made in ref. [57] is therefore just as sensitive to the uncertainties discussed as it is in conventional simulations. A high precision study of the spectrum is therefore required to identify the correct extrapolation.

Conclusions
In this paper we have studied the dynamics of the global strings that form in QCD axion models when the U (1) PQ symmetry is broken after inflation. Using numerical simulations, we have shown that the string network approaches an attractor solution that is independent of its initial conditions, and which is approximately scale invariant. We have also seen that this solution has a number of interesting properties: • The string length per Hubble patch ξ(t), defined by eq. (9), increases logarithmically with the ratio of the Hubble parameter and the string core size, with best fit parameters given in eqs. (12) and (13).
• At any time more than 80% of the string length is contained in long strings and the rest is in loops of size Hubble and smaller, which follow a scale invariant distribution. While such loops shrink and disappear they are replaced, at the same rate, by loops produced from longer strings.
• The energy density of the string network is determined by eq. (10) with -nontrivially -an effective string tension µ eff that is close to the theoretical expectation eq. (32). This means that to a good approximation, the energy density of the string network in the scaling regime at a particular time can be determined solely from the density of strings ξ(t).
• Over the scale separations that can be simulated, a substantial, and approximately constant, proportion of the energy emitted by the string network goes into heavy radial modes. Contrary to expectations, the radial modes have therefore not decoupled from the string dynamics by the end of the simulations, despite the relatively large scale separations m r /H ∼ 500 reached at these times.
• The instantaneous spectrum of energy emitted into axions, plotted in Figure 11 in terms of F defined in eq. (22), has the theoretically expected shape with a peak at momentum around the Hubble scale. At higher momenta it follows a power law until the string core scale above which emission is suppressed. The slope of the power is such that the majority of energy emitted by strings goes into high momentum axions for the scale separations that we can study (although a greater number of low momentum axions are produced).
• In simulations with physical strings, the instantaneous emission spectrum of the scaling solution can only be evaluated close to the end of the simulation because the residual energy from the initial conditions has less time to redshift. The result obtained matches the conclusions drawn from our analysis of the overall spectrum. The fat string model ameliorates these problems, and for this system we see that F remains similar at different times. The slope of the power law shows no sign of changing as the scale separation increases, although the uncertainty is substantial.
Since the scale separations that can be simulated are far from the physically relevant values, extrapolation over a vast distance is necessary if predictions relevant to the relic abundance of the QCD axion are to be made. Given our results, it is plausible that ξ continues to grow logarithmically (although we cannot exclude the possibility that it eventually saturates). Such an increase would enhance the number of axions emitted by a factor ∼ 10. Having tested that the energy density in strings is well matched by eq. (11) with the theoretical expectation for the string tension, the remaining ingredient necessary for a precise prediction of the relic abundance is the instantaneous axion emission spectrum at large scale separations, in particular the power q in F . While q < 1 at the scale separations that we can study, there is space within the current uncertainty for a small logarithmic correction that would result in it increasing to ≥ 1 in the physically relevant part of parameter space. On the other hand, if q < 1 persists at large scale separations then the number of axions produced by strings is presumably sensitive to the UV completion of the theory, as a significant proportion of the energy is going into exciting heavy modes, although in this case the number density is suppressed relative to that from misalignment.
To determine whether q has a scale dependence it is desirable to carry out simulations at larger scale separations. Indeed, even with a relatively modest increase in range, evidence for a change in q might be seen, which would allow an extrapolation. Improvement could come from different directions. One possibility is to simply perform simulations on larger high performance computing clusters, after implementing more efficient parallelisation than we have in our present work. A more involved possibility would be to develop an adaptive mesh algorithm, in which the scalar field is simulated on a grid with finer mesh spacing in regions of space where this is required, close to the cores of strings. The mesh must be updated as the strings move, which requires a more sophisticated algorithm, however this approach has proved beneficial in simulations of, for example, astrophysical black hole mergers [58]. It may also be interesting to further study the dynamics of individual strings, or small numbers of strings, which could give an indication of the behaviour of the full network. A qualitatively different approach would be to develop a numerical simulation in which the strings themselves were treated as the fundamental degrees of freedom, with parameterised dynamics and interactions. Such an effective theory style simulation could allow much larger scale separations to be studied. In [59] this idea was used to numerically simulate "strings" in 2 dimensions, although the extension to three dimensions is challenging, and a careful analysis would be necessary to ensure that the full dynamics of the underlying theory is capture.
A full calculation of the relic abundance must also include the axions produced when axion mass turns on, at which time domain walls form and the string network is destroyed. After the axion mass first becomes cosmologically relevant, when m a (T ) ∼ H, it continues to increase fast, and quickly H (T ) m a (T ) f a . As a result, this system has three widely separated scales, and is even harder to numerically simulate than the string network alone. Additionally, if N W > 1 it is currently unknown if the explicit breaking permitted such that the axion still solves the strong CP problem allows the domain walls to decay fast enough to avoid axions overclosing the Universe. We leave a study of the dynamics at these times to future work. time τ , which is defined as The time-dependent mass in the fat string scenario, eq. (8), is m r (τ ) = ( is the time at which m r = m i . In this way, the equations of motion simplify to where ψ and ∇ x ψ are derivatives with respect to the dimensionless time and distance variables m r τ and m r x respectively (or m i τ and m i x for the fat string case). In the physical scenario u(τ ) = 1, while in the fat string case u(τ ) = τ i /τ . We then solve eq. (35) numerically on a cubic lattice with periodic boundary conditions. 13 Space is discretised in a box of comoving side length L c containing N 3 = 1250 3 uniformly distributed grid points, where the upper value of N is limited by our memory budget. Consequently the space-step between grid points in comoving coordinates is a c = L c /N , which is constant in time. It is again convenient to work in terms of the dimensionless comoving space-step m r a c in the physical case, and m i a c for the fat string system.
The physical length of the box is L(t) = L c R(t) and the physical space-step between grid points a(t) = L/N = a c R(t) grows ∼ t 1/2 . In the physical string scenario, the string core size m −1 r is constant, and therefore the number of grid points in a core (m r a(t)) −1 ∼ 1/t 1/2 decreases. Meanwhile, in the fat string scenario the string core size increases ∼ t 1/2 , and as a result the number of grid points in a string core (m r (t)a(t)) −1 = (m i a c ) −1 is constant. When considering systematic errors from the space steps in the fat string system we use the notation m r (t)a(t) for the size of the space step due to its direct physical interpretation (although this is entirely equivalent to m i a c ).
The equations of motion are discretised following a standard central-difference Leapfrog algorithm for wave-like PDEs (see e.g. [60]). The system is evolved in fixed steps of conformal time a τ , and we work in terms of the dimensionless time-step m r a τ in the physical case, and m i a τ in the fat string case. The derivatives are expanded to fourth order in the space-step and second order in the time-step. 14 In Appendix B we extensively study the systematics from the discretisation of space and time, as well as from finite volume effects.
As discussed in Section 3, we set the initial conditions in two ways, the second of which is used to produce a cleaner initial configuration with a fixed number of strings per Hubble volume. for | k| ≤ k max and zero for | k| > k max for a fixed choice of k max ∈ [0, m r ] in each set of simulations. Larger values of k max lead to initial fluctuations with smaller wavelength, of order 2π/k max , and more initial strings. As a result, the initial string density is controlled by the parameter k max /m r . Even though this method produces configurations in which the axion field winds the fundamental domain [0, 2πf a ] nontrivially, it does not lead to a clean string configuration at the initial time, because |φ| does not typically resemble the profile function of string-like solutions. Instead, the system takes some time to relax to a string network solution and in doing so releases a large quantity of energy, which produces extra contamination to the axion spectrum measured at later times. To solve this issue, we employ the method below that also allows us to construct axion field configurations with a predetermined string length, or equivalently with predetermined initial ξ.
(b) Fixed string number: This approach simply consists of evolving a field configuration produced by method (a) with k max = m r until the required string length inside the box is obtained, and then using that configuration as initial condition for the actual simulation, with a different initial value of the Hubble parameter. Since this involves resetting the Hubble parameter, the strings produced do not in general have the right core-size, but we have checked that they quickly relax to their expected thickness long before the scaling regime is reached.

A.2 String Identification and String Length
To identify strings and measure their length, we first identify grid points that are likely to be close to a string core. We have carried this out in two ways, and have verified that the results obtained are extremely similar. In our main approach, we flag points such that, as a loop that surrounds it is travelled, there is at least one change in the axion field ∆a between consecutive lattice points encountered that satisfy |∆a| /f a > π/2. In particular, the loop is taken to be a square of side length 2 grid points. In order to capture strings with all possible orientations, at each point in the grid we consider loops in three orthogonal planes, and we flag a point if a loop in any of these satisfies the condition. We have checked that the results obtained are extremely similar for any reasonable choice of the threshold value for the change between adjacent grid points.
Having identified points close to the string core, we then combine these into strings. In particular, we cluster together flagged points that are adjacent in the x-y plane into an individual string point located at the mean of the flagged points. If such a cluster has a non-zero overlap with a cluster at the next level up or down in the z direction, these are connected into a string segment. As expected, the reconstructed strings form loops, and the length of each individual loop (which is required to analyse the loop distribution) as well as the total string length is recorded.
We have validated our string identification algorithm by comparing to the results obtained following the procedure adopted in [50]. In this, a lattice plaquette is identified as containing a string if the minimum axion field range that includes the field values on the four surrounding vertices satisfies |∆a|/f a > π. It can be seen that this leads to the correct results for the prototype string solution, and indeed for any string solution for which the 2π field change incurred in traversing an enclosing path is distributed sufficiently homogeneously around the loop.

B Analysis of Systematic Errors
To check that our results are free from numerical artifacts, we study how the key observables depend on the unphysical parameters in our simulations. Of particular importance are: (1) the lattice spacing a, (2) the time spacing a τ , (3) the number of Hubble patches HL contained in the box, and (4) the way that the string cores are screened when evaluating the axion spectrum and energy. In this section we study the first three of these, postponing (4) until Section C. The results that we present here are all obtained from the fat string scenario, however we have tested that the conclusions we reach are also valid for the physical case.

B.1 Lattice Spacing
We study the dependence of ξ(t) and of the axion spectrum on the discretisation parameter m r a to find the largest value compatible with the continuum limit m r a = 0. As mentioned, m r a is constant and is equal to the number of grid points per string core at all times in fat string simulations (in contrast, for the physical system the number of points per string core decreases with time, and this source of systematic errors only becomes relevant towards the end of simulations). If m r a 1, the number of gridpoints per string core is much smaller than 1 and the evolution of strings is not correctly captured. Conversely, making the space-step unnecessarily small would restrict the maximum scale separations that could be analysed.
We have performed sets of simulations averaging over 20 samples, keeping all parameters fixed except m r a, with identical initial conditions.
In Figure 14 we plot the continuum extrapolation for ξ and its dependence on log(m r /H) for the different values of the space-step, normalised to those from the finest lattice spacing tested (m r a = 0.67). For all m r a ≤ 1.33, ξ(t) differs from the most precise result by less than 1%. However, for larger values of m r a, it is systematically larger, especially at later times. In the main text, for our results of ξ(t) and the loop distribution in the fat string scenario, we used the rather conservative choice m r a = 1.33. In the physical scenario the number of grid points inside a string core decreases as the simulation progresses. We use parameters such that, at the final time, we match the resolution in the fat string case m r a = 1.33 (which means that at early times there are many more grid points inside a string core).
In Figure 15 we plot the axion spectrum at the time log(m r /H) = 6, for different space steps, again normalised to the results with m r a = 0.67. This is slightly more sensitive to the value of m r a than ξ(t) is, and we see that discretisation effects increase the production of UV states, reducing the energy emitted in the form of low momentum axions. In particular, only simulations with m r a ≤ 1 seem to have converged to within few percent of the continuum limit, while m r a = 1.33 and 1.6 result in a spectrum that is systematically smaller in the IR. In the main text we have chosen m r a = 1 in our analysis of the spectrum in the fat string system; this value is sufficiently small that we are confident that the UV dominated spectrum obtained is not an artifact of the finite space-step. Similarly, in the physical case we use parameters such that m r a = 1 when the final time shot used in calculating the spectrum is taken.

B.2 Time Spacing
We study the dependence of our results on the time-step similarly. Given the form of eq. (35), the general theory of numerical solutions of PDEs tells us that the relevant quantity to which m i a τ a τ /a c =1/5 a τ /a c =1/4 a τ /a c =1/3 a τ /a c =1/2.75 a τ /a c =1/2. should be compared is the comoving space-step m i a c [60]. 15 Consequently, we perform sets of simulations that differ only in the value of a τ /a c , averaging over 10 samples. In Figure B.2 we plot the results for ξ(t) and the spectrum. For all values of a τ /a c , ξ(t) is within 0.5% of the continuum limit, but for a τ /a c ≤ 1/3 the difference is less that 0.1%. The spectrum is also affected at less than percent level for a τ /a c < 1/3, and we use a τ /a c = 1/3 for our simulations in the main text in both the fat string and the physical scenarios.

B.3 Finite Volume
Finite volume effects are a particularly delicate issue. These are controlled by the parameter HL, which counts the number of Hubble lengths per box side length, with the physical limit for a spatially flat Universe corresponding to (HL) −1 → 0.
Since HL decreases over the course of a simulation, we want to determine the smallest value that it can take at the final time, H f L f , such that the results obtained match those in the physical limit. Although we use periodic boundary conditions, if H f L f > 2 every point is causally disconnected from itself from the beginning to the end of the simulations. Nevertheless, even choosing H f L f > 2, finite volume effects might still affect some observables such as the loop distribution (see Section 3.3).
We performed a set of simulations keeping the values of all the parameters except for H f L f fixed (all ending at a final scale separation log(m r /H f ) = 6). In Figure 17 we show the convergence of ξ to the infinite volume limit for HL > 2 for different values of log(m r /H) and as a function of log(m r /H) (where we indicate with vertical lines the times at which HL = 2 for each simulation). Finite volume effects do not affect ξ while HL > 2, and they remain small even slightly later. However, for HL 1 they result in a dramatic change, since the the whole network is in causal contact and starts to be destroyed.
We also study how the loop distribution is affected by finite volume effects. In Figure 18 we show the ratio ξ /ξ ∞ , defined in Section 3.3, for different values of , for a set of simulations with H f L f = 2 and log(m r /H f ) = 6.7. Vertical lines correspond to the times at which HL = 5, 4, 3. As mentioned in Section 3.3, the constant value of ξ /ξ ∞ in time is a strong indication that the system is in the scaling regime. However, the finite box size limits the maximum loop radius that can be contained. As a result, finite volume effects lead to ξ /ξ ∞ growing at late times once larger loops are no longer possible, and the larger is the earlier this occurs. For our analysis of ξ and the loop distribution we have chosen H f L f = 2 for both the fat string and the physical systems. This is rather conservative for ξ, but slightly sub-optimal for the study of the distribution of relatively long loops with lH/π 1/2 (although earlier time shots are also plotted in Figure 4, so the effect of the finite volume is clear).
In Figure 19 (left) we plot the axion spectrum at log(m r /H f ) = 6 for simulations with different values of H f L f . The peak at the Hubble scale is well reproduced for H f L f ≥ 2, even if its position is dangerously close to the IR cut-off (indeed, both are of order Hubble). On the other hand, for H f L f 1, not only is the Hubble peak not present, but there is also a significant overproduction of UV modes, related to the shrinking of a significant fraction of the string network, which towards the end of the simulation is made up of loops with radius smaller than the Hubble distance.
In Figure 19 (right) we plot the axion number density n a as a function of log(m r /H) for simulations with different values of H f L f ; with vertical lines corresponding to times at which HL = 3 for each simulation. n a has been normalised to the values obtained in the simulation with H f L f = 4, which is the least affected by finite volume effects. The number density seems to be slightly more sensitive to the finite box size than ξ. While HL > 3 the effect on n a is less than one percent, for HL < 3 finite volume effects result in a systematic underestimation of a few percent. This is reasonable, since n a is strongly dependent on the IR of the axion spectrum, which is the part most sensitive to finite volume effects. For this reason in the main text we chose H f L f = 3 when studying both the spectrum and the number density, for the fat string and also the physical system. We have also studied the dependence on H f L f of the effective string tension µ eff . This depends on the distribution of string length in loops of different shapes and sizes, since e.g. the logarithmic divergence in the string tension is expected to be cut off at a smaller scale for small loops than long strings. Consequently, µ eff could in principle be sensitive to finite volume effects. However, we do not observe any change in the string tension for H f L f ≥ 2.

C.1 Components of the energy
In this section we discuss in more detail the way that the total energy density of the system ρ tot splits into the contributions defined in eq. (5), and demonstrate that our method of screening the strings when extracting the energy in free axions and radial modes is consistent and does not introduce systematic errors. We also study the energy contributions in the absence of strings, and the effect of the interaction term between axion and radial modes. Finally, we describe the algorithm that we use to calculate the axion spectrum.
As summarised in Section 4.1, the total energy density is written as the sum ρ tot = ρ s + ρ a + ρ r . The energy density in free axions ρ a has been calculated as the spatial average ρ a = ȧ 2 far away from string cores, and in these regions axion and radial modes are to a good approximation decoupled at sufficiently late times. Similarly, the radial energy density is computed from ρ r = 1 2ṙ 2 + 1 2 ( ∇r) 2 + V (r) , again away from the string cores. We define the parameter d s such that the minimum distance from the strings' centre of the regions that are included in the average is d s m −1 r , so that, at all times, d s is the screening distance in units of the core size for both the fat string and physical scenarios.
There is actually an alternative approach to calculating the energy density of axions, which does not require screening. This consists of including the interaction terms in the axion energy density, i.e. evaluating (1 + r/f a ) 2ȧ2 over the whole space. On the string centre the factor (1 + r/f a ) 2 vanishes, and the contribution from the core regions is thus suppressed. Calculated in this way, the axion energy density will include the contribution from the axion-radial interaction energy. When we discuss the energy contributions to ρ tot , we will show that energy densities in axions calculated in both ways, i.e. with ȧ 2 screened and with (1 + r/f a ) 2ȧ2 , are in close agreement.
In Figure 20 we plot the ratio of different components of the energy density to ρ tot as a function of m r t/2π for a single simulation of the fat string system carried out on a small grid, with a spacestep m r a = 1. The goal here is not to understand how energies behave in the scaling regime (which has not been reached in the simulation analysed here), but rather to check that our approach to calculating the energy in axion and radial modes is consistent. All of the energy contributions plotted are calculated as the spatial average over regions a distance of at least d s = 1 away from string centres, except for 1 2 (1 + r/f a ) 2ȧ2 , which is averaged over the whole space. The simulation has been run until the box only contains 1/8 of a Hubble volume, i.e. H f L f = 1/2. As a result the entire system is in causal contact, and all the strings are destroyed at about m r t/2π ≈ 19, before the end of the simulation.
Many key features concerning the way that the total energy is split up can be understood from  Figure 20: Different contributions to the total energy density ρ tot , each one normalised to ρ tot , as a function of time for fat string system. All contributions are calculated as the spatial average at a distance d s = 1 from the strings' center, except for 1 2 (1 + r/f a ) 2ȧ2 , which is calculated over the whole space. The string network begins to be destroyed by the finite box size at m r t/2π ∼ 19, corresponding to the purple line going to 1. Figure 20. First, the kinetic part 1 2 ȧ 2 of the axion energy outside the string cores is systematically smaller than the gradient part 1 2 ( ∇a) 2 outside the string cores while strings are present. This is expected since the former only gets contributions from free axions (in the approximation that the motions of strings does not have a significant effect onȧ far from string cores), while the latter contains energy from both axion waves and a large fraction of the tension of strings. However, once the string network starts shrinking 1 2 ( ∇a) 2 decreases, and both 1 2 ȧ 2 and the radial energy increase. This indicates that as the strings are destroyed the energy stored in their tension is transferred to free axion and radial modes.
As expected for free classical waves, after the network shrinks the axion energy 1 2ȧ 2 + 1 2 ( ∇a) 2 is repeatedly interchanged between its kinetic and gradient parts, which on average are equal. Similarly, the total energy in radial modes 1 2ṙ 2 + 1 2 ( ∇r) 2 + V (r) behaves as expected (for free states with a mass decreasing adiabatically in time, due to the fat string trick), with 1 2ṙ 2 and the sum 1 2 ( ∇r) 2 + V (r) giving equal contributions. After the string network disappears, the ratios of axion and radial energy to the total energy stay approximately constant. This is a sign that axions and radial models (as well as the total energy) redshift at the same rate, i.e. as massless radiation, or equivalently as massive radiation with a mass decreasing with time as ∼ 1/R(t).  Figure 21: The spectrum calculated with method (1), described in the text, for different screening distances d s , normalised to the spectrum calculated with method (2).
The two ways of computing the axion energy density: evaluating ȧ 2 outside strings or (1 + r/f a ) 2ȧ2 over the whole space, are compatible within few percent while strings are present. The small discrepancy arises because the factor (1 + r/f a ) is not exactly a step function at the edge of the string cores. As we will see in the calculation of the spectrum, the difference is stored in UV axion modes (as expected) and therefore do not affect the calculated axion number density. Once the strings disappear the two measurements of the energy density give almost identical results since the difference is due to the non-vanishing interaction energy 1 2 (r 2 /f 2 a + 2r/f a )ȧ 2 , which goes to zero as the Universe expands.
The sum 1 2ȧ 2 + 1 2 ( ∇a) 2 + 1 2ṙ 2 + 1 2 ( ∇r) 2 + V (r) converges to the total energy ρ tot only once all the strings have been destroyed. Prior to this, a significant proportion of the energy is stored in string cores, and this is largest at small values of the log. The difference corresponds to part of the string energy, with the remainder coming from the difference 1 2 ( ∇a) 2 − 1 2ȧ 2 away from the string cores. The small mismatch between 1 2ȧ 2 + 1 2 ( ∇a) 2 + 1 2ṙ 2 + 1 2 ( ∇r) 2 + V (r) and ρ tot in the absence of strings is due to the non-vanishing axion-radial interaction energy, which however decreases as the Universe expands.
We also note that the energy densities in axions and radial modes have small oscillations with frequency O(m r ). The amplitude of the oscillations decreases with time due to redshifting, and as a result for a fixed value of the log they are larger in the physical case compared to the fat string scenario. Calculating the instantaneous axion emission spectrum involves taking the difference between the axion spectra (appropriately redshifted) at two times, which are separated by more than m −1 r . Consequently, the result obtained for momenta around m r is sensitive to the phase of the oscillations at the two times. This is the origin of the fluctuations of the instantaneous spectrum in Figure 11, which as expected are more significant in the physical case. To reduce this effect in the physical case we averaged over simulations starting at slightly different times, with relative differences of order 2π/m r .

C.2 The Axion Spectrum
We now describe the details of the calculation of the differential axion spectrum ∂ρ a /∂k. In the absence of strings, the axion energy density ρ a = ȧ 2 is where x p = R(t)x are physical coordinates, andã(k) is the Fourier transform ofȧ(x p ). The axion spectrum ∂ρ a /∂|k| is defined by d|k| ∂ρ a /∂|k| = ρ a , and is therefore given by In the following, as in the main text, we use the notation ∂ρ a /∂k ≡ ∂ρ a /∂|k|. In performing the integral dΩ k in eq. (37), which converts the three-dimensional spectrum to the one-dimensional spectrum, the three-dimensional momenta k ≡ 2π n L with | n| ≤ N 2 and | n| ∈ m − 1 2 , m + 1 2 , m ∈ N, are grouped into the same one-dimensional momentum bin labeled 2πm L . Momenta with N 2 < | n| < √ 3 N 2 have not been included since they lie far above the string core scale and do not matter for either the spectrum or the number density.
In the presence of strings, eqs. (36) and (37) are no longer valid because of the contribution of the string cores toȧ(x p ) and so toã(k). We adopted two methods to mask strings out of the calculation of ∂ρ a /∂k: 1. The first is the Pseudo Power Spectrum Estimator (PPSE) introduced in the analysis of cosmic microwave background data [61], and first used in the context of cosmic strings in ref. [50]. This method involves calculating the spectrum ∂ρ a /∂k after removing the regions of space that are less than d s m −1 r from the strings' center. Briefly, it is based on eq. (37) withã(k) taken to be the Fourier transform of θ(x p )ȧ(x p ), where θ(x p ) = 1 for x p a distance of more than d s m −1 r from the strings' center, and θ(x p ) = 0 otherwise. ∂ρ a /∂k is then appropriately corrected to account for the bias introduced by the window function θ(x p ) (more details may be found in [61]).
2. The second approach takes advantage of the automatic screening of string cores provided by the factor (1 + r/f a ), and consists of using eq (37) withã(k) taken to be the Fourier transform of (1 + r(x p )/f a )ȧ(x p ). This is similar to the previous method, except that the string cores do not have to be identified, and the degree of masking varies smoothly as the core is approached. We have tested that the computation of the axion spectrum ∂ρ a /∂k is independent of the method used, and in particular we studied the dependence of the first method on the screening distance d s . In Figure 21 we plot the ratio of spectra at log = 6 computed using the two methods for different values of the screening parameter d s , including the case d s = 0, where the strings are not screened. We see that the difference between the spectrum from method 1 with d s = 1 and the one from method 2 is 5% at k = m r /2 and rapidly decreases for smaller momenta. This is expected since, given the form of the string profile function for the radial mode, method 2 is roughly equivalent to masking distances up to m −1 r . Meanwhile, increasing d s suppresses the spectrum at momenta of order 2π/(d s m −1 r ), but leaves IR momenta unchanged. Again this is not surprising, since masking distances of order d s m −1 r only loses information about momenta larger than 2π/(d s m −1 r ). The spectrum for d s = 0, which is affected by the presence of strings, is much more UV dominated than the others, but at momenta lower than m r /2 still does not differ by more than 10%.
In the main text we calculated ρ a and ∂ρ a /∂k using method 1 with d s = 1. The resulting systematic uncertainty from the masking algorithm at k ∼ m r /4 is less than 1%, while for k ∼ m r /2 it can be estimated to be of order 5 ÷ 10%. As explained in Appendix D, we extract the power law of the spectrum between the IR and the UV peaks, q, by considering momenta smaller than m r /6. As a result the systematic uncertainty on q introduced by the masking procedure is negligible. The fact that ρ a and ∂ρ a /∂k are stable under the change of algorithm, and that ρ a redshifts as radiation, gives us further confidence that, even in presence of strings, our results only include free axion modes, as required.

D Convergence from Different Initial Conditions
In Section 3.1 we demonstrated the existence of the attractor solution by considering the convergence of the string density ξ(t), starting from initial conditions with fixed values of ξ. Here we present further evidence for the attractor solution. First, rather than fixing initial conditions with a predetermined initial string density by the method (b) described in Section 3.1 and A.1, we use random initial conditions given by the method (a) varying the maximum coefficients k max /m r . Second, we show how other key properties of the string network, in particular the axion spectrum and axion number density, also converge to the same late time values starting with different string densities.
In Figure 22 we plot ξ(t) starting from random field initial conditions for different choices of k max /m r at the initial time H i = m r . Although the results obtained are slightly less regular than those starting with fixed numbers of strings, shown in Figure 3, the convergence to the attractor, and the logarithmic increase, is clear.
To analyse the convergence of other properties of the string network, we present results starting from fixed string number densities, which lead to the values of ξ(t) plotted in Figure 23 (left). We have also confirmed that the quantities that we study converge to the same attractor solution starting from random initial conditions. In Figure 23 (right) we show the results for the axion number density obtained from the same set of simulations as ξ(t) in the left panel.
Notice that while at log = 3 the values for ξ and n a are spread respectively by a factor of 3 and 15, at log = 6 the spread reduces to around 10%.
The convergence of the axion spectrum ∂ρ a /∂k to the attractor, shown in Figure 24 (left), is also revealing. At log(m r /H) = 3, the axion spectrum is highly suppressed in simulations with a lower initial string density than in the others. However, by log(m r /H) = 6 the spectra are very similar. This is especially the case for the modes around Hubble, which have nearly the same amplitude for all initial conditions. More dramatically, the instantaneous emission shown in Figure 24 (right) converges extremely fast at late times. This is because, unlike the overall spectrum, it only depends on the properties of the network at a fixed time instant, which are practically the same for all initial conditions by the end of the simulations. We see that underdense networks tend to have more instantaneous emission in the IR at early times, which is expected since at these times the typical distance between strings is larger in this case. Taken all together, our results provide convincing evidence that the properties of the string scaling solution at late enough times are sufficiently insensitive to the initial conditions chosen.

E Extraction of the Scaling Parameters
In this appendix we describe how we extract the slope α in the fit of ξ(t), eq. (13), and the power law q of the instantaneous spectrum, and estimate the associated uncertainties. We also describe how we obtain the extrapolations of the number density in Figure 13.
Given the convergence shown in Figure 3, initial conditions that result in a constant value of dξ/d log(m r /H) ≡ dξ/d log for the longest time correspond to string networks that are the closest to the scaling regime. Therefore, we estimate α by calculating dξ/d log for initial conditions with different fixed initial string densities, and then restrict to those such that dξ/d log is approximately constant at late times. 16 The constant common value that dξ/d log reaches in such simulations is a good estimate of α, and the spread indicates the uncertainty.
In Figure 25 we plot dξ/d log for different initial string densities, for both the fat string and  Figure 25: The derivative dξ/d log(m r /H) for different initial conditions. The black curves are approximately constant for the largest amount of time and therefore the closest to scaling, and their spread is used to estimate the error on α. Blue dashed lines represent the estimate of α, and are plotted at times when the system is reasonably close to scaling. Error bands represent statistical errors in the average over different samples. the physical case. Slopes of different simulations tend to converge asymptotically at to a common constant value. This is evident in the fat case, for which the constant approached is dξ/d log = 0.22 ± 0.02, where the error has been estimated from the spread of the results from simulations that have constant dξ/d log for log(m r /H) 4 (plotted in black). In the physical case, the scaling regime is reached at larger values of the log and dξ/d log is changing over most of the simulated range. As a result, the uncertainty on dξ/d log = 0.15 ± 0.05, estimated as the spread of dξ/d log over the simulations for which dξ/d log is constant for log(m r /H) 5, is larger.
We now turn to the extraction of the power law of the spectrum q. For this purpose, we consider simulations starting from the initial conditions that lead to string networks that are the closest to scaling (corresponding to results for ξ(t) plotted in black in Figure 3). 17 As shown in Figure 11, for log(m r /H) 5 the instantaneous emission spectrum, parametrised by F (k/H, m r /H), has an approximate power law behaviour 1/k q for momenta that are large enough with respect to the peak at around the Hubble scale (at k/H ∼ 5 ÷ 10) and are small compared to the UV cutoff at k = m r /2.
The fact that q is less than 1 can be more clearly seen from Figure 26, in which we plot xF (x, y). In these plots q < 1 corresponds to the increase between the IR and UV peaks, which is evident both in the fat string and physical cases. Moreover, most of the area under the curves in Figure 26 is at UV momenta k ∼ m r /2, which shows that the energy density is dominantly emitted in the form of very high momentum axions, although the axion number density is dominated by   Figure 27: The axion energy density spectrum ∂ρ a /∂ log k for the fat string (left) and the physical (right) cases (the data is the same as in Figure 10 but represented on a log scale). low momentum states since q > 0. In Figure 27 we plot the total energy spectrum of Figure 10 on a log scale, i.e. ∂ρ a /∂ log k. The positive gradients show that the energy density in axions is dominated by UV modes at all times, for both the fat string and the physical case.
Due to the challenges in the physical case, discussed in Section 4.2.1, we only attempt a detailed analysis of q in the fat string scenario. To do this, we consider F at late times, log(m r /H) 5, and focus on the region with momenta a factor of 3 larger than the Hubble peak and a factor of 3 smaller than the core peak, i.e. 30H k m r /6, which is sufficiently uncontaminated by the cutoffs. The power q is then given by −d log F (x, y)/d log(x). 18 In Figure 28  The black sections corresponds to momenta 30H < k < m r /6, which are sufficiently far from the IR and UV peaks to be uncontaminated by effects at these scales. The roughly constant values in these regions correspond to the value of −q in the approximate power law, and the blue dashed lines represent our estimated range for this.
The sections of the curves plotted in black indicate momenta in the range 30H < k < m r /6, which are safe from contamination from the IR and UV peaks.
d log F/d log(x) crosses zero at about k ∼ 5H, signalling the position of the Hubble scale peak, and it is smaller than −1 at large momenta, meaning that after the UV cutoff the instantaneous spectrum falls off steeply. Moreover, in the intermediate region of interest −d log F/d log(x) changes slowly as a function of momentum and across different time shots, ranging from about 0.7 to 0.8. We note that in this range −d log F/d log(x) (and therefore q) shows a tendency to increase at later times. However, given the present uncertainties this apparent change is not significant enough to draw any conclusion, in particular we are not able to assess whether it is due to a residual contamination from the nearby UV peak or it constitutes a genuine feature of the spectrum (similar analysis with larger grids and more statistics could be beneficial to understand this point). In the regime 5 log(m r /H) 6, we therefore estimate q = 0.75 ± 0.05, where the uncertainly mostly comes from d log F/d log(x) not being constant over the momentum range of interest.
In order to extrapolate the axion number density with different power laws, we used an analytic form of F (x, y) that closely matches the three approximate power laws visible in Figure 11, where θ(x) is the step function. In this expression x 1 and x 2 are the positions of the IR and UV cutoffs in units of Hubble, q 1 , q 2 are the powers that suppress the spectrum in the IR and in the UV respectively, and q is the power between the two cutoffs. The normalisation N is required so that F (x, y)dx = 1. A reasonable fit for the 5 free parameters is x 1 ∼ 3, x 2 = y/2 and q = 0.7 ÷ 0.8, q 1 ∼ 3, q 2 ∼ 2 in the fat string case, and similarly in the physical case except that q 2 ∼ 3. The extrapolations in Figure 13 are carried out by keeping all parameters fixed except for q. What matters most for the final extrapolation is the parameter q, which regulates the hardness of the spectrum, while the remaining parameters only lead to changes of order 1.

F Comparison with EFT Estimates
As mentioned in Section 4, the dynamics of global string loops can be equivalently described by an effective theory where the fundamental degrees of freedom are the string and the axion radiation [34], with an interaction governed by the Kalb-Ramond action [62]. This theory is valid in the regime where the string and the emitted radiation (with frequency ω ∼ 1/R) are not strongly coupled, which corresponds to log(m r R) 1, where R is the loop radius (as shown in [34], the effective coupling of the emitted axion radiation to the string is proportional to 1/ log(m r R)). As a result, for loops with a large hierarchy between the radius and the core size, the emitted axion radiation should approximate the one predicted by the Nambu-Goto effective theory. Indeed, using such a theory it was shown [34] that a circular loop starting at rest with log(m r R 0 ) = 100 follows the cosine time-law for the Nambu-Goto strings with percent precision, at least for values of the loop radius such that log(m r R) 1, where the theory is applicable. We show now that the evolution of circular loops provided by the solutions of the field equations matches the one predicted by the effective theory of strings, for the values of log(m r R 0 ) reachable in our field theoretic simulations. We solve eq. (6) in Minkowski space, i.e. with a time independent scale factor, H = 0, and initial conditions φ(x) andφ(x) that approximately resemble a static circular loop with initial radius R 0 . 19 In Fig. 29 we plot the time-law for the loop radius R(t) normalised to the initial radius R 0 for log(m r R 0 ) = 4 and 5. We also plot the prediction given by the effective theory for log(m r R 0 ) = 5 and the free Nambu-Goto time law, R NG (t) = R 0 cos(t/R 0 ). The result of the simulation for log(m r R 0 ) = 5 matches very well the EFT prediction where this is valid. Moreover, as log(m r R 0 ) increases, the circular loop time law gets closer to the free Nambu-Goto prediction, indeed indicating that at large log(m r R 0 ) global string dynamics converge to that of free Nambu-Goto strings.
Although the convergence R(t) → R NG (t) is good when the loop radius is sufficiently larger than the core size m −1 r , i.e. for log(m r R) 1, there is a substantial deviation when the string becomes strongly coupled, log(m r R) 2, which prevents the loop from bouncing many times. This however is not in conflict with the EFT prediction, which correctly reproduces the right time law before the loops collapses even for log(m r R 0 ) as small as 5.
Given the limitations in performing direct field theoretic simulations at much larger values of log(m r R 0 ) we cannot test whether the EFT expectation that loops will bounce many times (thus emitting an IR dominated spectrum) is correctly reproduced or whether core effects when the loops shrink keep inhibiting the rebounce.