Axion dark matter: strings and their cores

Axions constitute a well-motivated dark matter candidate, and if PQ symmetry breaking occurred after inflation, it should be possible to make a clean prediction for the relation between the axion mass and the axion dark matter density. We show that axion (or other global) string networks in 3D have a network density that depends logarithmically on the string separation-to-core ratio. This logarithm would be about 10 times larger in axion cosmology than what we can achieve in numerical simulations. We simulate axion production in the early Universe, finding that, for the separation-to-core ratios we can achieve, the changing density of the network has little impact on the axion production efficiency.


Introduction
The axion [1][2][3][4] is a proposed particle, the angular excitation of a new "Peccei-Quinn" (PQ) field ϕ that would solve the strong CP problem [5][6][7] and which is also a very interesting dark matter candidate [8][9][10], thereby solving two puzzles with one mechanism.In this light, the study of the axion as a dark matter candidate is, in our view, well motivated.
The axion model has one undetermined parameter, the vacuum value of ϕ, f a ; the axion mass m a scales as f −1  a .The value of f a also determines the amount of axion dark matter that would be produced in the early Universe, which means that (for a well-motivated initial condition we will describe) it should be possible to predict the axion mass from the dark matter density.To do so, we need to understand the efficiency of axion production in cosmology.
As we will soon explain, the PQ field's evolution in the early Universe is complicated by the appearance of structures -cosmic strings -that may play a role in determining the production efficiency of axions around the QCD scale T ∼ 1 GeV.Here we will study these structures and their role in axion production via numerical simulations.Ours is not the first study on this problem [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27], but we believe it will also not be the last; we will demonstrate some significant challenges for such studies, one of which we are not yet able to overcome.Axion cosmic string networks are in a family called "global string networks," and Martins and Shellard have argued [27,28] that such networks are sensitive to the sizes of the string cores, which cannot be properly simulated numerically.This implies that numerical simulations will view string networks with very different properties -in particular, a much lower string density -than the ones that would really occur for physical axion field parameter values.We will present strong numerical evidence supporting this claim.In particular, we will demonstrate that the cosmic string density increases as we increase the log of the ratio of the network age to the string core size.This ratio is very large in cosmology, implying that physically relevant string networks may be an order of magnitude denser than those in numerical axion-production simulations.We will also point out another potential pitfall to numerical simulations, which can be overcome but must be monitored.
We use the remainder of the introduction to review the properties of axions relevant to cosmology and the overall picture of axion production in the early Universe.From the point of view of cosmology, the relevant features of the axion are that it is a complex scalar field ϕ with a symmetry-breaking potential 1 , with f a the vacuum expectation value and λ the self-coupling, which we assume is O (1).
The Lagrangian has a U(1) "PQ" symmetry ϕ → ϕe iθ , spontaneously broken when ϕ takes on a vacuum value somewhere on the "vacuum manifold" 2ϕ * ϕ = f2 a .The symmetry is also explicitly broken by an anomalous coupling to QCD and instanton effects, giving rise to an extra contribution to the potential, − L axion,eff.QCD = χ(T ) [1 − cos (arg ϕ)] . (1.2) Here χ(T ) is the temperature-dependent topological susceptibility of QCD, which for our purposes is an input from the theory of QCD.Its vacuum value 2 is χ(0) ≃ (76 MeV) 4 .
The minimum of the axion potential is at √ 2ϕ = σe iθa with σ = f a and θ a = 0. Radial fluctuations σ = f a + s, which we will call saxions, have a mass m 2 s = λf 2 a ; and angular excitations θ a = a/f a , axions, have a mass m 2 a (T ) = χ(T )/f 2 a .If f a ≫ Λ QCD , as required by existing constraints [31][32][33], then there is a large hierarchy between these masses, and only the axion should play a role cosmologically for T < f a .χ(T ) is known to be a strong function of T at high temperature T ≫ Λ QCD , varying roughly as T 7+N f /3 = T 8 [34].Therefore as time progresses and temperature falls cosmologically, the axion goes rapidly from being effectively massless, m a t ≪ 1, to massive, m a t ≫ 1.The dynamics around m a t ∼ π are quite nontrivial; but once m a t ≫ 1, axion fluctuations will be small and the axion number an adiabatic invariant, so the dynamics around m a t ∼ π determine the efficiency of axion production and therefore the amount of axion dark matter.
It seems likely that PQ symmetry is restored near the end of or after inflation.In particular, PQ symmetry is generically restored during high-scale inflation (if the Hubble parameter during inflation is large, H ≥ f a ) [35] 3 .Thermal symmetry restoration after inflation would have occurred if the Universe reheated to a temperature T ≥ f a .If either occurred, then we know the initial conditions for the axion field: θ a would start out uncorrelated at causally-disconnected points, meaning essentially random initial conditions should apply.With the initial conditions known statistically, the axion cosmology model has only one unknown parameter, f a .If we can determine χ(T ) and work out the axion dynamics near m a t ∼ π, then it should be possible to compute the relation between the value of f a and the resulting dark matter density.If we assume that the axion constitutes the dark matter in the Universe, the known dark matter density [38] should give a unique prediction for f a , and with it the axion mass.The axion mass is experimentally measurable, making this scenario testable.And a narrow search window would also be very valuable for the design of the most sensitive type of experiment, resonant microwave cavity detectors [39][40][41][42].
To finish setting the stage, we describe qualitatively how the fields evolve near m a t ∼ π.At early times, m a ≃ 0 and the Lagrangian Eq. (1.1) has a U (1) symmetry.This symmetry is locally broken by the phase choice ϕ = f a e iθa (x) .Random initial conditions give rise to a network of topologically stable cosmic strings [43], which then evolve and untangle, entering a scaling solution in which they maintain a string length per volume of order where γ, L string , and V are the typical Lorentz gamma-factor of a moving string, the physical length of string, and the physical volume under consideration, respectively, and t is the age of the Universe.ξ is an order-1 parameter describing the network evolution.
After the axion mass becomes relevant, the potential has a single global minimum 4 , and there are no true topological structures.However, the strings remain metastable, and there are metastable domain walls associated with the field's phase varying by 2π from one side of the wall to the other; exactly one wall ends on each string.The domain wall tension draws the strings together and accelerates the annihilation of the network [44].
We will argue that the "scaling solution" for the strings is not quite as simple as Eq.(1.3) suggests; the string tension, the parameter ξ, and other properties of the string and string-wall networks should vary logarithmically with the ratio of the inter-string separation t/ √ ξ ∼ t to the string core size ≃ 1/m s .In numerical field-theory simulations this ratio is bounded by the size of the lattice used to solve the field dynamics: m s t will not exceed about 10 3 in 3D simulations or 10 4 in 2D simulations.Physically, we want m s t ∼ f a /H, with f a ∼ 10 11 GeV and H ∼ T 2 /m pl ∼ 10 −19 GeV; so the physically relevant value for the ratio is more like m s t ∼ 10 30 .For a fixed f a , the physical value of m s t gives strings with 10 times higher tension than in simulations, leading to a much denser string network which breaks up more slowly under the action of domain walls.The physical network evolution may be much more efficient at producing axions than recent simulations (including the ones we present), implying a smaller value of f a and heavier value of m a .
The other danger is that even the meta-stability of the domain walls is dependent on the mass ratio m a /m s .We will show below that, for m 2 a /m 2 s > 1/39, the domain walls cease to be even metastable, and the topological structures abruptly collapse.This certainly should not happen cosmologically; but since m a rises rapidly with time, it is a challenge to make a numerical simulation with m s large enough that this condition is maintained until the network breaks up via the expected physics.Therefore, there is some danger that simulations will incorporate unphysical collapse of the string-wall network or must be stopped before the network has finished evolving.
The next section reviews the physics of global cosmic strings and explains the relevance of ln(m s t), using analytical arguments.Next we present numerical evidence that the string network grows denser with increasing ln(m s t).Then we study the network evolution and axion production around m a t ∼ π in more detail.Our study of the string-wall network's collapse actually does not show evidence that the final axion density depends strongly on ln(m s t); but the dynamic range we can study is too narrow to make any brave extrapolations about what happens when the log is increased by another factor of 10.
As an aside, we mention that besides the details of the axion dynamics around m a t ∼ π, there are also uncertainties in the axion production efficiency from our incomplete knowledge of the temperature dependence of the topological susceptibility χ(T ).The hightemperature behavior is only known at first order in perturbation theory [34], so the first unknown corrections are suppressed by O(α s ).Therefore the perturbative treatment may not be very reliable around T ∼ 1 GeV where the interesting dynamics occurs.In this temperature range we only have model calculations [45,46]; lattice calculations [47,48] are currently available only at lower temperatures and/or in the quenched approximation.We will find that the axion production is not too sensitive to the exact value of χ(T ), but it would nevertheless be valuable to have a reliable lattice calculation for T in the range from 0.5 to 1.5 GeV.

String networks and string cores
We begin with a lightning review of why the strings arising from Eq. (1.1) have logarithmically large string tensions.We assume some familiarity with string defects; a reader who needs some background can look in [49].We will make the classical field approximation throughout, which is an excellent approximation for the IR axion field dynamics since the mean occupancy is ∼ f 2 a /H 2 ∼ 10 60 .

Cosmic strings
Consider a string lying along the z axis in polar coordinates.The field varies as ϕ = f a f (r)e i(φ−φ 0 ) with φ the azimuthal angle and f (r) a function obeying f (r) . The associated energy is (2.1) Far from the string's core, f → 1 and all terms become negligible except for the φ-derivative term, which becomes f 2 a /2r 2 .Therefore the energy density decays as 1/r 2 as we move away from the string core.This is in contrast to "local" strings, where the U(1) symmetry is gauged (local) and a gauge field compensates for the φ-derivative except in the string's core.The local string case has received much more attention in the literature.
The large-distance part of the string tension (energy per length) in Eq. (2.1) is approximately where ℓ is an IR length scale where this description breaks down, for instance the distance to the next string, ℓ ∼ 1/H.This dependence on the distance to the next string indicates that there is a long-range force-per-length (gradient of string energy-per-length) between strings, with a strength of order d(E/L)/dℓ = πf 2 a /ℓ.There are two important facts here.First, there are long-range interactions between strings, mediated by the (nearly) massless axion field.The force between strings of opposite winding sense is attractive, which helps them to find each other and annihilate.Though we have not shown it, the presence of a massless mode also helps accelerating or bent strings to radiate energy more efficiently than for the local string case.Both of these effects help the string network to annihilate more efficiently, leading to a much lower-density string network.
Second, the attractive forces between strings, and the radiation of energy into longwavelength axions, scale with f 2 a in the same way as the string tension, but they are not enhanced by ln(m s ℓ), while the energy-per-length, or string tension, is.Therefore the ratio of tension to radiation/force effects is proportional to this log, but not to f a .To keep track of this difference, we will name this large logarithm κ ≡ ln(m s ℓ) ≃ ln(f a /H).

Scaling in 2D
Let's first see how κ plays a role in string density for 2D networks, where sensitivity to this logarithm is generally accepted [11].The 2D network is described by the 3D theory but enforcing that the field ϕ does not vary along the z axis, ϕ = ϕ(x, y).In the x − y plane, the string becomes a monopole or point-charge, with ϕ(x, y) varying by ±2π as one goes around the charge, depending on whether the string has positive or negative orientation.The two orientations of strings act like two signs of charges.
The analogy to electric charges turns out to be complete [50][51][52][53][54]. Outside of the string cores, we can write ϕ = f a e iθa , and the equation of motion becomes Except in string cores, we can define dual electric and magnetic variables obeying which are the free-space Maxwell equations in 2+1 dimensions.A surface in 2D is a loop, and the electric flux through the surface is 2π times the winding number around the loop, showing that a string is the source for a flux of ±2πf a .The electrical attraction between two strings will be as expected if each has an energy πf 2 a ln(m s r) as we found above 5 .The short-distance part of this energy, πf 2 a ln(m s r 0 ) = πf 2 a κ, should be interpreted as the mass of the charge, m = κπf 2 a .Varying κ (the log of the ratio of separation to core length scales) is varying the charges' mass at fixed charge magnitude.We want to argue that making the charges heavier will make them more non-relativistic, so they move and annihilate less efficiently and have a higher density.
Suppose that at t = 0 there is a very dense starting ensemble of positive and negative charges, which evolve under Hubble drag and their electromagnetic interactions.We can find how the density of charges will evolve by making parametric scaling estimates.The charges find each other under the influence of Coulomb attraction and annihilate off.Suppose that at time t the mean inter-charge separation is ℓ, so the density of charges is 1/ℓ 2 , and let us try to estimate ℓ.On dimensional grounds we expect ℓ ∼ κ −n t, with n an exponent we want to determine.
We see that the density of charges should fall, with O(1) of the charges annihilating in an O(t) amount of time.To do so, the charges have to move a distance ∼ ℓ in a time t, implying that v ∼ ℓ/t.The typical force on a charge is F ∼ 2πf 2 a /ℓ, so the kinetic energy a charge obtains is force times distance, The density of strings is n ∼ ℓ −2 ∼ κ/4t 2 .Therefore we estimate that the number density of strings should increase linearly with the log of the scale ratio κ, and that the velocity should scale as 1/ √ κ, becoming nonrelativistic as the logarithm becomes large.
We therefore have a robust argument that, in 2 space dimensions, the string density will scale with the log of the scale separation κ ≡ ln(m s t) as n ∼ κ/t 2 , with small squared velocities v 2 ∼ 1/κ.

Scaling in 3D
The argument in 2D clearly does not translate simply into 3D, since 3D strings move under tension as well as under mutual interactions.Nevertheless we expect strong, though not necessarily linear, κ-dependence in the string network density in 3D as well.The reason is that the long-range interactions of the strings provide rather efficient mechanisms for the strings to radiate energy and to annihilate against each other.This is in contrast to local (gauge) string networks, where there are no long-range interactions and only very inefficient radiation of gravitational waves.As a result, the density of global strings, ξ < 1.2 at achievable κ values, is an order of magnitude smaller than the value ξ ∼ 13 found in numerical lattice field theory [27,55] and Nambu-Goto [24,56,57] evolutions.
However, while the interactions between strings and the radiation of energy from accelerating strings should scale as f 2 a , the tension of strings should scale as f 2 a κ.In the large-κ limit, the tension should dominate the mutual string interactions and radiation.Therefore, as κ is increased global string networks should behave more like Nambu-Goto string networks.For very large values of κ, the network density ξ should be a factor of 10 larger than at currently achievable κ values.Martins and Shellard attempted to model this by amending their 1-scale model for string networks to include the radiation effects associated with global strings [28].Fitting the single parameter describing radiation efficiency in their Eq.(4.2) to the value of ξ = 1.15 at κ = 6 from Fig. 3 below, and solving their equations describing the string density and velocity, we get the results plotted in Fig. 1.According to this model, the string density should vary nearly linearly with κ for κ ∼ 5-10, and should grow by roughly a factor of 5 in going from κ = 6 to κ = 70.
We will compare numerical simulations to these expectations in the next section, by considering string evolution without the "tilt" term, Eq. (1.2).

Implementation
To investigate these issues, we have implemented an MPI-parallelized code to evolve the equations of motion that follow from Eq. (1.1) and Eq.(1.2).We work in a radiation dominated FRW background, so the Hubble parameter6 H = ∂ t a H /a H is related to the time t as H = 1/(2t).We work in comoving coordinates and introduce conformal time with t 0 and a H0 the values of t and a H at some (arbitrary) fixed time.Similarly, we rescale masses by a H so that m phys.dt = m conf.dτ ; in what follows, all masses are m conf. .The temperature scales as T ∝ τ −1 , and the metric is In these units and writing space (i, j) and time (0) components and factors of τ from the metric explicitly, the action, up to an irrelevant multiplicative constant, is We discretize this on a grid that is uniform in these coordinates, with spacing a and temporal spacing a t = a/n t .We give a few more details about our numerical implementation in Appendix A; in particular we explain there how we determine the density of strings and walls, and what we believe is a new technique for establishing the string velocity directly from information in the field variables.Our numerical implementation is rather standard, except for two points.First, like many authors we replace with ϕ r = √ 2 Re ϕ.This form has the advantage of being analytic in field variables, so no trigonometric evaluations are required.It differs from the original form only in string cores, where we expect this term to have little impact; the symmetry-breaking term will always be small, so it is only important because it applies over large regions of space where the field is near its minimum.We assume that χ(τ ) grows as a power of τ in the relevant time regime, so if χ(T ) ∝ T −n then τ 2 χ(τ ) = f 2 a (τ /τ 0 ) n+2 , two powers higher than the temperature dependence of χ(T ).We also scale out an overall factor of f 2 a from the Lagrangian and rescale ϕ → ϕ/f a .
The other nonstandard change we make is to the potential term.Physically, we are interested in large values of m 2 s = λf 2 a .We cannot make the value larger than m 2 s a 2 ∼ 1, since otherwise we encounter numerical artifacts (as discussed in Appendix A); but physically we want it to be as large as possible, since we expect the physical value of m s to be very large.Therefore we remove two powers of τ from this term so that its τ scaling is the same as the gradient terms, so m s a remains fixed in our simulations.In other words, we keep the string core size fixed in lattice units, rather than allowing it to grow smaller in comoving coordinates due to Hubble expansion.This change enlarges the dynamic range where our simulations can see cosmic string networks, making it easier for us to achieve "scaling."

Results for string-only networks
As a first application, we investigate string networks without the symmetry-breaking term, that is, setting χ(τ ) = 0 or m 2 a = 0.In this case, no domain walls form, and the string network evolves until we terminate the simulation at τ ≤ L/2.
We have plotted the length (number) of strings, scaled by τ 2 to account for system expansion, for 3D and 2D simulations in Fig. 2. The left plots are in terms of conformal time measured in units of the saxion mass, τ m s ; they show a rise, at first rapid and then more gradual, in the string density when normalized by a τ 2 factor.To test the hypothesis that this rise represents logarithmic scaling in the separation-to-core ratio κ ≃ ln(τ m s ), we present a semilog plot on the right-hand side of the figure.Indeed, the nearly straightline behavior is striking in both 2D and 3D.In each figure we have shown two or more choices of lattice spacing m s a, and we indicate the 1σ statistical errors with upper and lower thin curves accompanying each best-value thick curve.The figure is also a check that our rather coarse lattice, am s = 1.5, is sufficient to see continuum behavior; we plot curves for am s = 1.5, for am s = 1.0, and (in 3D) for am s = 0.7, which fall on top of each other to within the small statistical errors.
It is common in the cosmic string community to consider also the mean velocity, v 2 , and Lorentz γ-factor of the strings.It is more customary to normalize the network density in terms of the time rather than the conformal time; according to Eq. (3.1) this divides7 the string densities shown in Fig. 2 by 4. In addition, it is customary to track not the length of string network, but the energy in the network divided by the string tension.This is the same as integrating γdl rather than dl in determining the effective length of string in the network.In computing the mean velocity, v 2 , and γ-factor one also uses this normalization, so We plot the normalized network density in 3D and 2D, using this normalization, in Fig. 3.The figures indicate that the γ-weighted string density is also a logarithmic function of τ m s , and that the string density in 3D is over an order of magnitude smaller than the value ≃ 13 for local string networks [56,57].
We also plot the values of v , v 2 , and γ in Fig. 4. The string velocity falls off at large τ m s in 2D as we expected, but in 3D it roughly approaches a constant, which is quite close to the value from Fig. 1.The figure shows that, while m s a = 1.5 was a sufficient value  to obtain continuum-limit string densities (Figure 2 and 14), our algorithm for finding the string velocity (see Subsection A.2) is more sensitive to the size of the string core, and a high-quality determination of the string velocity and especially of its γ-factor by this method requires m s a ≤ 1 in 3D.We also see oscillations in string velocity at small m s τ ; we believe that these are string-core "breathing" modes induced by our initial conditions, which may fool our string-velocity method, since it assumes that the string core takes its unperturbed form.This section's results show that the string density in 2D evolves exactly as we would expect: the string density grows with κ = ln(τ m s ), while the string v 2 falls as κ −1 .In 3D the string network density grows in a way which, for κ values we have available, is roughly consistent with linear behavior.The string velocity in 3D is a very weak function of κ and is near the value predicted by the one-scale model [28].
4 String-wall network evolution

Scaling expectations
Consider the axion field in conformal coordinates.As we discussed, the axion mass grows strongly as temperature drops, and therefore as conformal time progresses.We will designate a special time τ 0 , defined as when Roughly speaking, this is when the explicit symmetry breaking and axion mass start becoming physically significant.We will parameterize the topological susceptibility near this time as χ(T ) ∝ T −n .In our numerical work we will take n = 7.When m a (τ )τ ∼ π the dynamics are complex and must be solved numerically.But for n = 7, by τ = 3τ 0 we have m a (τ )τ ∼ 420, which should provide enough field oscillations to force the string-wall network evolution to completion.By this time (dm a /dτ )/m 2 a ≪ 1, ensuring adiabatic evolution of the axion field.In our comoving and conformal coordinates, one then expects adiabatic evolution of the form ε ∝ m a f 2 a /τ 2 and8 n ax ≃ ε/m a ∝ f 2 a /τ 2 .On dimensional grounds, in the radiation era and before the QCD transition when the number of radiation degrees of freedom changes, the axion number should be with K a constant.This constant determines the produced density of axions; since little entropy is produced between GeV temperatures and today, the axion-to-entropy ratio in the modern Universe is where T 0 is the temperature when τ = τ 0 and g * is the effective number of light degrees of freedom at T = T 0 .Therefore, determining K is determining the key input for the axion abundance in the modern universe.
Our goal is to determine this constant by evolving the axionic field, starting with a random space-varying phase at very early time, until the string network is gone and the axion fluctuations are small, and then evaluating which should be independent of the final time τ when we evaluate it, if that time comes after adiabatic behavior sets in.Eq. (4.4) can be evaluated by fast Fourier transform (FFT) methods, or by the method of Appendix B when the FFT is not available.As a baseline for the expected value of K, we can consider the angle-average of the misalignment mechanism; we evolve a spatially uniform initial condition for θ a , leading to with m 2 a = τ n+2 /τ n+4 0 , and then take the average of the resulting axion number over values of the starting angle θ a .For our choice n = 7 we find K = 16.0255.

String-wall network evolution in 2D
We have evolved the lattice axion equations of motion for both 2D and 3D systems, using the coarsest lattice that gives continuum-like string behavior, m s a = 1.5 (see App. A.3).We consider a number of values for the dimensionless ratio m s τ 0 .Depending on one's perspective, a large value for this ratio is either a large value of m s or a large value for the time τ 0 ; we prefer the former, and will express all other time and frequency scales in terms of τ 0 .We make m 2 a ∝ τ n+2 throughout the simulation, that is, we assume that the QCD scale (where m 2 a stops varying) is reached after the simulation ends.We terminate each simulation when m 2 a a 2 = 0.1, so the axion mass is coming on order the lattice spacing.At this time, we evaluate the axion-number content of the field and extract K.
We begin with 2D simulations because they are numerically cheap, so we can achieve quite large values of m s τ 0 .We look first at how the string and wall densities vary with time.Our results, scaling out the expected N str ∝ τ −2 , A wall ∝ τ −1 behavior, appear in Fig. 5.Much of the figure is as expected.The scaled string and wall densities rise with τ below about τ = 1.4τ 0 , due to logarithmic scaling corrections we have already discussed.They are also larger for larger values of m s τ 0 , also as a result of logarithmic corrections to scaling.The wall length starts to fall at about τ = 1.2τ 0 , and the strings begin to move faster around τ = 1.5τ 0 , leading to a brief peak in ξ as the strings' γ-factor rises and then a fall in the string length and extent starting around τ = 1.9τ 0 .Eventually, the extent of both strings and walls falls near zero.For larger values of m s τ 0 , the strings are heavier by a factor of ln(m s τ 0 ), and so they respond to the walls with more inertia.This delays slightly the growth in string velocity and extent, and slightly delays the fall of the wall and string extent.The γ-factors we find at late times, while large, are actually underestimates, as the γ-factor is especially sensitive to the lattice spacing, as we explain in Appendix A.2.
But the figure shows something that might at first be unexpected.For each value of m s τ 0 , there is a point where the wall area abruptly crashes.At the same value, the string number briefly spikes, and then crashes as well.The larger the value of m s τ 0 , the larger the τ 0 value where this occurs.In every case it occurs when m 2 a /m 2 s = 1/(43 ± 3).This collapse of the wall network is a numerical artifact which we will explain in Subsection 4.3.It occurs because the ratio m 2 a /m 2 s is finite in these simulations, and when it becomes large enough, about m 2 a /m 2 s = 1/39, the walls cease to be even metastable and become absolutely unstable.This leads to an abrupt unraveling or collapse of the network.But this collapse has nothing to do with the physics that would occur in any circumstance where m 2 a /m 2 s remains large.We should question the physical relevance of any simulation where this occurs before the string-wall network has largely broken up via the expected physics.We see that this constrains us to consider quite large values of m s τ 0 .
To see this a little better, we plot the axion number determined via Eq.(4.4) (full details in Appendix B) for all but the two largest m s τ 0 values in Fig. 6.The figure clearly shows that larger m s τ 0 values, meaning denser string networks, lead to larger axion number at intermediate times.As the string network breaks up, the axion number falls.When the network abruptly unravels, the axion number drops sharply, and then goes flat after the network is gone.Except for the smallest m s τ 0 , the final value is almost unchanged.Still, the figure leaves us distrustful of overinterpreting simulations where the string network crashes rather than breaking up via natural string-wall network dynamics.Roughly, the m s τ 0 = 900, 600, and 450 simulations look reasonable, but the abrupt drop is rather large for the smaller-m s τ 0 simulations.
As a final way of looking at the network, we plot the energy content as a function of time in Fig. 7.The energy content naturally scales as f 2 a /τ 2 on dimensional grounds and due to Hubble expansion, so we have factored this out.Also, at late times the rising mass causes the energy to rise, ε ∝ m a ; but this behavior only sets in once the network starts responding to m a , at about τ = 1.4τ 0 .Therefore we have also scaled out a factor of (5 + m a τ 0 ), with the 5 chosen so the scaling will turn on at about τ = 1.4τ 0 .To avoid crowding the plots we have not shown error bars.The largest m s τ 0 curves have statistical errors around 2%; the smaller m s τ 0 have smaller statistical errors.In the energy-fractions plot we have shown only the three largest m s τ 0 values.The string energy is estimated as γL str πf 2 a ln(m s r) with r −2 = m 2 a + τ −2 , so the IR cutoff is either the inverse string separation or the axion mass, whichever is larger.The domain wall energy is estimated from domain wall area, without an attempt to include the γ-factor for the walls' motion, thus the wall energy is probably an underestimate.The string energy is an underestimate after γ gets large around τ = 1.8τ 0 , since our method underestimates large string velocities when using such large m s a values.
The figure shows that the energy starts out strongly m s τ 0 dependent, presumably because of the different density and string tension of the string network.By the end of the simulation this m s τ 0 dependence has largely disappeared.Also, the early behavior is dominated by gradient energy, since strings at rest carry almost all their energy in gradients (moving strings have a kinetic energy fraction of v 2 /2).The late behavior is oscillating axions, with energy equipartitioned between kinetic and (potential+gradient). Quadratic fluctuations have 1/2 their energy as kinetic, whereas walls and strings have most energy in gradients and potential; so the extent that the (phase-averaged) kinetic energy is below 1/2 is a reasonable estimate of how much energy is still in strings an walls.Once the energy is mostly in quadratic fluctuations (particles), there is still a shift between gradient and potential energy, due to the increase in m a .The oscillations at late time show that the produced axions are not a phase-random collection, but have some phase coherence.
We should worry about the impact of the abrupt collapse of the string-wall network if the strings and walls still carry significant energy fraction when it happens.The bottom right figure gives a nice criterion for knowing if this happens.In each case shown in the plot, less than 5% of energy in walls when they collapse.This is not true for smaller m s τ 0 values, so the associated simulations cannot be completely trusted.

Domain-wall instability
The most striking feature of the string-wall network simulations we just presented is the sudden collapse in the amount of domain wall, accompanied by a brief spike, and then collapse, in the amount of string.The time value τ /τ 0 when this occurs changes as we vary m s τ 0 .But it always occurs at the same value of m 2 a /m 2 s ≃ 1/43.Here we show that this is not an accident; it happens because the domain walls are only metastable, and at m 2 a /m 2 s = 1/39 they become completely unstable.A domain wall is where the phase θ a changes from 0 to 2π by going around the "valley" of the winebottle potential.Since θ a is a single real parameter, this occurs on a codimension-1 surface.The thickness of the region where θ a is far from 0 or 2π is set by a competition between potential energy, which wants to make the region thin to keep the region with potential energy costs small, and gradient energy, which wants to make the region thick to minimize |∇ϕ| 2 .The thickness of the region is ∼ 1/m a ; when m a τ ≫ 1, this is thin compared to the horizon scale and the walls can be approximated as planar.
In the case m 2 a ≪ m 2 s , we can take √ 2 ϕ = f a e iθa to good approximation.Choosing the θ a -variation along the z axis, the energy to minimize is with boundary conditions θ a → z→−∞ 0, θ a → z→∞ 2π.Here A is the area of wall we consider, and σ = E/A is the surface tension.Extremization gives ) This can be solved by multiplying both sides by ∂ z θ a , and integrating dz to obtain a Virial-type relation showing that potential and gradient terms each represent half the wall's energy.The surface tension is 9 Numerically, we estimate the energy in domain walls as 8m a f 2 a times the area where θ a = π, which we determine as described in Appendix A. Now suppose that m 2 s /m 2 a is large but not enormous.Then ϕ(z) will not strictly lie on the circle f a e iθa , and Eq.(4.6) becomes (introducing φ = ϕ/f a ) 12) The equations of motion are As z is varied, the field φ will follow a curve in the complex φ plane, where ℓ is the affine parameter describing the curve that φ follows; |d φ/dℓ| = 1 and |∂ z φ| = dℓ/dz.The complex equation of motion, Eq. (4.13), can be decomposed into the (r, i) component tangent to the curve and the component normal to the curve.The tangent EOM reads integrating into a Virial relation in the same way as before.The surface tension is again The equation of motion in the field direction normal to the curve is what determines the curve φ(ℓ); it reads This equation is equivalent to asking: what curve φ(ℓ) will produce the lowest surface tension in Eq. (4.16)?By choosing a path with a small value of V , we keep the integrand 9 Some references use the value 9.32 rather than 8.This estimate originates from [58], who find it for the domain wall at T = 0, essentially by incorporating corrections to the strict (1 − cos θa) form of the potential.Such corrections arise at zero temperature because instantons form a correlated liquid; but at high temperatures the dilute instanton gas approximation should be good and we expect such corrections to be tiny.
small; but by choosing a short path, we keep the integration range short.The LHS of Eq. (4.17), ∂ 2 φ(ℓ)/∂ℓ 2 , is the extrinsic curvature of the path.The larger its value, the more we can shorten the curve by rounding it off at this point.This must be balanced against how fast V will rise when rounding off the curve.The equation tells the exact balance between the gain of shortening the curve and the cost of increasing V along the curve.Now V − V0 in the denominator is mostly provided by m 2 a , while the normal derivative in the numerator is almost purely provided by m 2 s .The larger we make m 2 a /m 2 s , the less there is to stop the domain wall from rounding off its path and exploring φ2 r + φ2 i < 1. a /m 2 s path as it appears on the (ϕ r , ϕ i ) potential, illustrating how the curve departs from the "valley" of lowest potential.
We have solved explicitly for the domain wall's shape and surface tension, using V from Eq. (4.12) with different values of m 2 a /m 2 s .We illustrate the results in Fig. 8.The left plot in the figure shows several ϕ(ℓ) curves in the ϕ r , ϕ i plane; the curves shown have m 2 a /m 2 s values spaced in intervals of 0.0032.The dots are the values of φ(ℓ(z)) at a series of z-coordinates, with neighboring dots on a curve separated by ∆z = 1/(4m s ).The wall becomes thinner as the dots become more spread out.The last curve is the last metastable domain wall; a tiny further increase in m 2 a will cause the curve to pull over the top of the potential peak, and the domain wall spontaneously collapses.On the right in the figure, we have shown the path of this domain wall in terms of the potential V (ϕ r , ϕ i ), indicating how the curve "pulls up" partway onto the bump in the potential rather than staying in the "valley" of the potential.
Our study finds that the domain walls lose their metastability when m 2 a /m 2 s reaches about 1/39.For larger m 2 a /m 2 s values, there is no longer any domain wall solution.In simulations we observe the collapse of the wall area to start when the ratio is ≃ 1/43.We believe that this is because fluctuations in the field, hitting the wall, can induce a collapse when m 2 a has not quite reached the value where instability occurs.This will not happen everywhere at once, but only locally where larger fluctuations impact the wall.At these spots the wall will "break"; a loop of string forms the boundary of this break, temporarily raising the total length of string in the simulation.The hole in the wall then rapidly grows and is joined by new breaks, leading to the collapse in the wall area; as the holes in the wall grow and percolate, the amount of string at first rises but then falls essentially to zero as the wall network disappears.This is a good description of both the timing and behavior of the wall collapse we observe.
Physically, small values of m s are experimentally excluded.And it is most natural to expect λ ∼ 1 so m s ∼ f a , which should be orders of magnitude larger than m a .So while this physics is the correct dynamics of a string-wall network with m 2 a /m 2 s ∼ 1/40, it does not describe the evolution of physical interest.Therefore we cannot rely on the results of any simulation in which a significant amount of energy still resides in the string-wall network when this collapse occurs.This criterion places a limit on what values of m s τ 0 give reliable answers, pushing us towards fine lattice spacings and towards the largest values of m s a that still give continuum behavior.

3D simulations
Our experience with 2D simulations tells us that we must consider the largest possible values of m s τ 0 to avoid the unphysical collapse of the string network.This means we need to choose the largest m s a we can, subject to the constraint m s a ≤ 1.5, found in Sec.A.3 to ensure continuum behavior.We should also choose the largest aτ 0 value we can, subject to the constraints that τ /τ 0 must get large enough for the network evolution to complete, with L ≥ 2τ to ensure no finite-volume errors.Unfortunately, our limited numerical resources limit us to boxes of 1600 3 or smaller, constraining us to consider aτ 0 = 300 or m s τ 0 = 450 but not larger.We have not considered m s τ 0 smaller than 225, since the 2D simulations showed that the wall network breaks up too early in such simulations to learn anything of value.So we have less dynamic range than in the 2D simulations.
The 3D simulations show surprisingly similar behavior to those in 2D.Rather than repeating all plots we presented in 2D, we plot just the string length and speed, the wall area, and the energy ratios in Fig. 9.The most significant difference from 2D is that, while the strings and walls start to decay at about the same time, the network decays more quickly, with very little string left by τ = 2.57τ 0 , when the network collapse occurs for m s τ 0 = 450.In particular, the collapse of the walls, in the plot of wall area and of energy ratios, is almost invisible for this largest m s τ 0 value.Therefore, while this final simulation does not provide the right string tension due to the missing core tension, it at least presents a case where the network breaks up via a physical mechanism rather than the loss of wall stability.
For further comparison between the 2D and 3D cases, we have plotted the energy and the energy fractions for each dimensionality at a common value of m s τ 0 = 450, in Fig. 10.The energy density starts out somewhat higher in the 3D case, corresponding to the larger string density obtained in 3D; but the later stages of the evolution are strikingly similar, except that the 3D string/wall network decays before the unphysical collapse, while the 2D network still carries some energy when the collapse occurs.We also plot the final axion-production efficiency K as a function of m s τ 0 , for both 2D and 3D, in Fig. 11.Besides the very smallest values of m s τ 0 , it is a very weak function.Surprisingly, K is about half of the angle-averaged misalignment-mechanism estimate, and it appears to be a weakly declining function of m s τ 0 .Note however that it is very dangerous to extrapolate based on this result, since the physical value of m s τ 0 ∼ 10 30 is about 27 orders of magnitude larger.

Discussion and Conclusions
The axion abundance is determined by dynamics occurring around the conformal time τ 0 where τ m(τ ) = 1; we have seen that the axion number is essentially fixed by τ = 3τ 0 .Entering this time period, there is a network of axionic cosmic strings, which evolves and breaks up due to the appearance of the axionic mass.
We have shown that in both 2D and in 3D, the initial density of the cosmic string network is a strong, roughly linear, function of the log of the horizon-to-core ratio ln(m s τ 0 ).In simulations we can make ln(m s τ 0 ) as large as ln(450) = 6 in 3D and ln(1800) = 7.5 in 2D; but in nature we expect it to be ln(∼ 10 30 ) ∼ 69.Therefore we can only study axion production (so far) for cases with a much more dilute starting string network than we expect in the cosmological axion context.
We also showed that, for too small a value of m s τ 0 , the network dynamics suffer an additional numerical artifact; the axionic walls and strings become absolutely unstable and the network suddenly collapses.This occurs whenever the ratio of angular-to-radial masses exceeds m 2 a /m 2 s = 1/39.This effect compels us to work at large (numerically expensive) values of m s τ 0 ; but we can achieve values large enough that the instability only occurs after almost all strings and walls have disappeared.
We do not believe that this problem influenced the results of [17] because their m 2 a /m 2 s ratio stopped increasing partway through their simulations and then remained fixed, preventing it from growing too large.On the other hand, for the axion to play a part in dark matter we expect τ 0 to occur at temperature T 0 ∼ 1.5 GeV, so this flattening-off behavior is unphysical, raising some questions about their results.Perhaps the most significant, and to us surprising, result of our simulations is that, even over a range of κ = ln(m s τ 0 ) where the starting string density varies by a factor of 2, the final axion number is unchanged at the 20% level, and in fact trends slightly down with increasing string network density.So at least within the range of string network densities we have been able to study, the string density appears to play little role in establishing the final axion number density.
Suppose this is the case, so the produced axion number is well described by Eq. ( 4.3) with K = 8.What do we learn about the axion decay-constant f a and mass m a ?Combining Eq. (4.3), the Friedmann equation for a hot gas, an estimate of the hot topological susceptibility from Wantz and Shellard [46], with Λ ≡ 400 MeV, n ≃ 6.68, and α W S = 1.68 × 10 −7 , and the relevant cosmological information from Planck [38] along with an estimate g * = 64 based on a plasma of (eµτ, ν 123 ,γ,g,udsc) with the QCD degrees of freedom contributing 80% of the free-particle value to account for strong interaction corrections at this temperature [59], we find This value gives m a = (76 MeV) 2 /f a = 18 µeV.The transition temperature is T 0 ∼ 1.5 GeV.In our simulations, most of the nontrivial 3D network evolution took place between τ = 1.4τ 0 and τ = 2.6τ 0 , corresponding to T = 1100-580 MeV.This is the range where we need to understand the topological susceptibility better -though f a is only sensitive to a rescaling of the topological susceptibility through the 0.079 power, so even a factor of 10 error in the estimate of Eq. (5.2) (as suggested by recent work [48]) makes a modest 20% shift in f a and m a .A significant change in our estimate for K would have a larger effect (though a very small effect on the relevant T 0 ).Over the range we studied, there was very little change to K in going from 2D to 3D and in varying the string network density by about a factor of 2. We have made no attempt to separate which axions arise from misalignment, which from strings, and which from walls; indeed it is not clear to us that doing so is either well defined or terribly useful.But it appears that there is not a large component strictly proportional to the density of strings.Nevertheless, we find it rather brave to extrapolate that the same value K = 8 should apply when the string network is 5 to 10 times denser than in the simulations we considered.
We believe that it is well motivated to look for a way of simulating axion production from global networks with much larger core tension.It is clear that the enormous factor we need cannot be achieved simply by shrinking the lattice spacing.Rather, a new strategy is needed.We see the need for a method to excise and treat explicitly the physics of the string core, by adding degrees of freedom to describe its tension and inertia.We are close to presenting such a technique for 2 dimensional networks, taking advantage of the dual electromagnetic description.Details will appear elsewhere.action is Here the time derivative term stretches between times τ n−1 and τ n , so we replace τ 2 in its coefficient with the product of the starting and ending time.The coefficient a τ,n on this term accounts for the dτ running from τ n−1 to τ n .For the space terms, which occur at time τ n , we replace τ 2 with τ 2 n , and we treat the amount of dτ contributing to the term to be half the interval before the term appears plus half the interval after the term appears, hence the factor (a τ,n + a τ,n+1 )/2.We have also implemented a next-neighbor improved gradient term, but most of our results are based on the above gradient term.
For uniform temporal spacing, variation with respect to ϕ i gives One often writes this in terms of the conjugate momentum The factor τ n−1 /τ n+1 , which arises from the overall τ 2 factor in the action, accounts for Hubble drag.The middle term is the gradient term, the final term is the radial potential term, and the ϕ r equation is the same but with the addition of a linear m 2 a term.
In practice we use a very fine temporal spacing at small τ , a τ = a/40 for τ < a and a τ < τ /20 thereafter, to avoid errors when the Hubble drag is large.This is probably unnecessary.At larger times we go over to a fixed value for a τ /a.
As initial values, we first set the field to be of unit magnitude and independent random phase (ϕ r (x) = cos θ x and ϕ i (x) = sin θ x , with each θ x chosen uniformly from [0, 2π)) and then apply smearing to a coherence length ℓ smear .We will check for dependence on parameters such as m s a, ℓ smear , and a t /a momentarily.

A.2 Strings and string velocities
We want to know the density of the string network, which requires identifying where the strings are.To do this we define the plaquettes that are pierced by a string, and we count these plaquettes, applying a statistical correction which we now explain.
Each plaquette has four corners, and we say that a string goes through a plaquette if the tetragon in ϕ space, whose corners are the ϕ values at the corners of the plaquette, encloses the point ϕ = 0, see Figure 12.The direction or sense of the string is determined by whether the origin is enclosed in a clockwise or a counterclockwise sense.
ϕ Re Imϕ Figure 12.Four points around a plaquette map to four points on the complex ϕ plane.If the resulting tetragon encloses the origin, we identify the plaquette as being pierced by a string.
Our algorithm for determining if a string pierces a plaquette is as follows.For the tetragon in Fig. 12 to contain a string, the real axis must be crossed twice with the same handedness (clockwise or counterclockwise); that is, the number of strings piercing the plaquette is half the signed sum of real axis crossings, where the sign is determined by whether the crossing is clockwise or counterclockwise.The real axis is crossed between points x 1 and x 2 with values ϕ 1 = ϕ(x 1 ) and ϕ 2 = ϕ(x 2 ) if Im (ϕ 1 ) Im (ϕ 2 ) < 0. The axis crossing is clockwise if Im (ϕ 1 ϕ * 2 ) > 0 and is counterclockwise if Im (ϕ 1 ϕ * 2 ) < 0. The winding number is the sum of + 1 2 (− 1 2 ) for each clockwise (counterclockwise) axis crossing, as we consider each pair of corners going clockwise around the plaquette.
We also identify which links have a domain wall go through them by finding links where Im (ϕ 1 ) Im (ϕ * 2 ) < 0 (the real axis is crossed) with Im (ϕ 1 ϕ * 2 ) Im (ϕ 1 − ϕ 2 ) < 0 (it is the negative half-axis which is crossed).This occurs an odd number of times for each plaquette containing a string, and an even number of times for any other plaquette.Each algorithm (string and wall) requires only multiplication and comparison; no divisions or trigonometric functions are needed, providing good numerical efficiency.We can sample the full lattice for strings every lattice unit of time (∆τ = a) while taking less numerical effort than the field updates.We believe that our plaquette identification is equivalent to that of Hiramatsu et al [16], but by avoiding explicit reference to angles we avoid the trigonometric evaluations needed there.
We make no attempt to connect plaquettes pierced by strings to identify the strings and their lengths.Instead we rely on counting plaquettes with string and links with wall.In 2 dimensions this works for strings, but in 3d (and in both 2d and 3d for walls) there are normalization issues.If the unit normal along a string is (l x , l y , l z ), then a length L of string will pierce L|l x | yz-plaquettes, L|l y | xz-plaquettes, and L|l z | xy-plaquettes; for different directions the total number of plaquettes is between L and L √ 3. We are only trying to determine the string density statistically, and all directions are statistically equally likely, so we can find the right average string density statistically by counting string-pierced plaquettes and dividing by the angle-averaged number of plaquettes per unit length of string, which is count-factor = where Ω x,y,z are the x, y, z-components of the unit vector in the φ, θ direction.The same overcounting factor applies for domain walls in 3D, since the number of links piercing the wall depends on the wall-normal in the same way as plaquettes depend on the string-normal.
In 2D the domain walls are overcounted by with f (r) solving the equation of motion The equation of motion near 0 enforces that The moving string solution at t = 0 is found by replacing x with γx; and the time derivative for the moving string is At lowest order in m 2 r 2 we find (A.9) The next-order term in Eq. (A.8) can be used to engineer a correction for fields close to but not at the center of the string: which we found by establishing the leading-order small-distance behavior of ϕ * ϕ and |ϕ * ∂ t ϕ| and finding a combination that would cancel the subleading contributions in ∂ t ϕ * ∂ t ϕ.
We apply this approach by sampling over points that are near the string's center, so that the second-order, ∼ m 2 s r 2 terms are small, and higher (uncomputed) corrections should be negligible.We choose for our sample of points near the string core the set of points on corners of plaquettes that are pierced by a string; points on more than one string-pierced plaquette are counted once per plaquette.For each plaquette we find v 2 /(1 − v 2 ) explicitly by averaging Eq. (A.10) over the four plaquette corners; we then compute the local value of v , v 2 , and γ .We average (1, v, v 2 , γ) over all pierced plaquettes, weighting with weight γ in order to follow the literature convention that the string network should be weighted by γdl (energy content), not dl (length).The numerical overhead is small because the computation need only be performed on sites that are identified as corners of a string-pierced plaquette.
The method assumes that points with distance ∼ a away from the string core are still "relatively near" the string's core, which is an expansion in (cm s a) 2 /16 < 1.Note that, in 3D, the average distance-squared of a point on a pierced-plaquette to the nearest point on the piercing string is x 2 = a 2 /2, so the sampled points are surprisingly close to the strings, and the higher-order corrections are expected to be small.
Another weakness of the approach is that it is based on the structure of a straight string; it will commit errors for bent, accelerating strings.However, in this case the string's length and velocity are anyways not uniquely defined; it is not clear to us that the method works any worse than other approaches.Similarly, it could be confused by any radial fluctuations (breathing modes) in f (r); but these modes are heavy and we do not expect them to occur with large amplitude except perhaps at first due to initial conditions.Also note that, since the algorithm determines γ 2 v 2 in terms of a manifestly positive expression, it never concludes that a string is moving faster than light speed, something which could in principle occur for some string-velocity approaches.
We mention one final challenge for our method.We have assumed that the string's shape correctly reflects the Lorentz contraction of the string's structure.On a sufficiently coarse lattice, Lorentz invariance is not respected in the string core and the Lorentz contraction will not be as strong as it should be in the continuum.This means that the method will under-estimate the γ factor of very fast strings on coarse lattices.If the main goal is a high-precision determination of string velocities and γ-factors, one should use a smaller value of m s a than we have used in most of our studies.
To test this last point, we have repeated the simulation of a m s τ 0 = 450 2D string network using m s a = 0.75 and m s a = 1.0 lattices, to compare with our results in the main text using m s a = 1.5.Most properties (energy, string length, wall extent, final axion number) agree very well; the final axion number is the same within 1% statistical errors.But once the strings start to move very fast late in the simulation, the finer lattice observes a significantly larger γ-factor for the strings, as shown in Fig. 13.
We think a similar velocity-finding approach could be applied to field-theory simulations of local strings (Abelian Higgs model simulations), but since the string's structure m s a = 0.75, 1.0, 1.5.Our string velocity method finds a higher γ-factor on the finer lattice.(The crash in the γ-factor at τ = 2.57τ 0 coincides with the collapse of the network by direct wall instabilities and should not be taken seriously.)then depends on two parameters (the Higgs mass and the gauge boson mass), the application is more complicated, with the constant c above replaced by some m H /m A dependent value; finding the NLO corrections would also be significantly more complicated.

A.3 Numerical tests
We want to find the largest values of m s a and a t /a that are compatible with a continuum interpretation, and we want to check for sensitivity to initial conditions such as ℓ smear .We will test these parameters using m a = 0 or string-only simulations, because there are then fewer scales to consider.
We start with m s a.If we choose this too large, the string core is 1 lattice site across, which is too small for the lattice to resolve properly.In this case, the UV edge of the integral in Eq. (2.2) becomes sensitive to the exact location of the string relative to the lattice, leading to a string energy varying periodically with period a, rather than being translation independent.This makes it possible for a string to "stick" in the most energetically favorable lattice location, which will interfere with string evolution and impede the annihilation of the network.To test for this problem, we made a series of evolutions with different values of m s a, but identical other properties as measured in terms of m s .These can be interpreted as varying the lattice spacing.
We show how the string length (or, in 2D, the number of strings) vary with time for several values of m s a, in both 2 and 3 dimensions, in Fig. 14.In both cases, the results are consistent within 2% for all m s a ≤ 1.5, but deviate for larger m s a, especially at later times.We did a similar study with the next-neighbor improved action (not shown), which showed continuum behavior slightly sooner, at m s a = 1.8.The reduced number of lattice points this allows, (1.8/1.5)d+1 in d space dimensions, does not quite make up for the factor of 2 in numerical cost for that algorithm, so we have generally stuck with nearest-neighbor interactions in the remainder of our work.
Note that the result in Fig. 14 clearly shows that the string length rises as ln(m s t), rather than approaching a flat value, both in 2D and in 3D.Therefore we already see the logarithmic corrections to scaling in this figure.If we re-plot the figure in terms of τ /a, rather than m s τ , the lines do not fall on top of each other; the relevant physical length scale for comparison is m s , which sets the string core size, not the lattice spacing.
Next we check for independence on initial conditions, by varying the length scale ℓ smear over which the initial conditions are smeared.The result, shown in Fig. 15, indicates that different initial conditions rather quickly converge to the same string network density, which again scales logarithmically with the system age as measured in m s units.This means the choice for this smearing length is not important.We typically choose 2.1/m s in this work.
We should also check what temporal spacing is sufficient to approximate continuous time.There is a Courant condition beyond which the update algorithm becomes unstable: a t /a = 2/ 4d + m 2 s with d = 2, 3 the number of space dimensions.We choose a/a t integer, so this limits the available values to a/a t = 2, 3, . ... The axion is very light, and at late (interesting) times we expect that the physics is primarily in long wave lengths except near string cores; so one might expect less severe dependence on this ratio than for some lattice  In each case we compare with an extrapolation to zero spacing, based on a t /a = 1/6 and 1/8 and assuming (a t /a) 2 scaling of errors, which is the expected scaling form.
problems.This turns out to be the case; as we show10 in Fig. 16, a temporal spacing of a t /a = 1/2 actually only leads to percent differences in string length and energy, compared to small temporal spacing.Nevertheless, out of paranoia we have generally used a t /a = 1/6 or 1/8 in other measurements, which should keep a t errors below the per-mille level.Finally, although there are strong theoretical arguments that the volume cannot affect the network's statistical properties so long as τ < L/2 (in conformal coordinates), we check this nevertheless in Fig. 17.The figure shows that τ must grow to nearly twice L/2 before significant finite-volume corrections occur; but they then rapidly become severe.All results in the rest of the paper consistently use volumes large enough that evolutions stop at τ < L/2.
We have generally re-checked the 2D-only results of this section with 3D simulations, which give consistent results but with poorer statistics because the added numerical costs of treating the 3D system make it harder to achieve high statistics.

B Counting axions
At the end of a simulation we have ϕ(x), from which we want to extract an axion number.To do so we first convert to the axion field amplitude θ a , and then extract the axion number stored in the field.The axion field amplitude θ a is determined as processing the data; they drop in size when the energy becomes one digit shorter so one more digit is written out.
Since the linear term, Eq. (3.3), shifts the absolute minimum of the potential slightly away from ϕ r = f a , we make an addative shift to ϕ r so that its minimum is at f a before applying these rules.
We then assume that this axion field obeys the quadratic Hamiltonian (writing θ = ∂ τ θ) 2 The Hamiltonian is diagonal in Fourier space (V is the volume of space) 3) The particle number associated with a mode with oscillation frequency ω = k 2 + m 2 a is the energy over the frequency, so the particle number density is Evaluating this in Fourier space is straightforward and is efficient when the FFT is available.If the box size is not a power of two or if the data is divided over processors in an inconvenient way, we can use Laplace methods instead.It is easy to write an algorithm to evolve θ, θ in dissipative "time" τ , with boundary conditions that θ(τ = 0) is the original value of θ a : ∂ τ θ a (x, τ ) = (∇ 2 − m 2 a )θ a (x, τ ) , θ a (x, τ = 0) = θ a (x) , (B.5) and the same evolution for θa (in terms of τ evolution, θa is considered an independent field).The reason to do so is that we can find θ a (x, τ ) by position-space methods, but we know that the evolution in Fourier space will be By evolving in τ and numerically implementing the first integral, which can all be done in coordinate space, we get the desired second integral, which correctly counts axion number.
We have compared this method to the FFT method where both are applicable and have confirmed that they give the same answer up to controllable numerical (τ -spacing and large-τ extrapolation) issues, which are easily held below 1% with smaller numerical cost than the time evolution needed to find the final θ a configuration.However, this method is not practical for making frequent determinations of n ax over the course of a simulation.

Figure 1 .
Figure 1.Estimated κ dependence of the string density ξ and velocity v according to the one-scale model of Martins and Shellard.

Figure 3 .
Figure 3. String network density in 3D (left) and 2D (right), varying with the log of conformal time.Different curves refer to different lattice spacings m s a.

Figure 4 .
Figure 4. String mean velocity, v 2 , and γ-factor in 3D (left) and 2D (right), varying with the log of conformal time.

Figure 5 .
Figure 5. String and wall density in 2D simulations, for a number of choices of m s τ 0 .Top left: string number.Top right: wall length.Bottom left: string ξ-factor.Bottom right: velocity, v 2 , and γ factor of strings.

Figure 7 .
Figure 7. Energy density, scaled by τ 2 /(f 2 a (5 + m a τ )), for 2D simulations.Top left: total energy.Top right: potential.Bottom left: gradient energy.Bottom right: energy fractions and estimated energy in strings and walls.

Figure 8 .
Figure 8. Left: domain wall solution, as seen in φ-space, for several values of m 2 a /m 2 s .Each curve is the path φ follows through field-space for a given m 2 a /m 2 s value.Neighboring dots are points that are separated by 1/(4m s ) in coordinate space.Right: the largest-m 2a /m 2 s path as it appears on the (ϕ r , ϕ i ) potential, illustrating how the curve departs from the "valley" of lowest potential.

Figure 9 .
Figure 9. String length, wall area, string velocity, and energy ratios as a function of time in 3 dimensions for 3 values of m s τ 0 .

Figure 10 .
Figure 10.Comparing 2D and 3D simulations at the same m s τ 0 value: Scaled energy (left) and energy fractions (right).

Figure 11 .
Figure 11.Axion production efficiency K as a function of m s τ 0 , in 2D and in 3D

Figure 13 .
Figure 13.String velocity for 2D, m s τ 0 = 450 simulations with three values of lattice spacing: m s a = 0.75, 1.0, 1.5.Our string velocity method finds a higher γ-factor on the finer lattice.(The crash in the γ-factor at τ = 2.57τ 0 coincides with the collapse of the network by direct wall instabilities and should not be taken seriously.)

Figure 14 .
Figure 14.Dependence of string density on the lattice spacing (in units of the saxion mass, am s ).Left: string length in 3 dimensions.Right: number of strings in 2 dimensions.Lines represent 1σ statistical range based on averaging several samples.

Figure 15 .
Figure 15.Dependence of string density, in 2D, on the initial correlation length of the random initial conditions.

Figure 16 .
Figure 16.Temporal spacing dependence of string density (left) and energy (right).In each case we compare with an extrapolation to zero spacing, based on a t /a = 1/6 and 1/8 and assuming (a t /a) 2 scaling of errors, which is the expected scaling form.