Approach to scaling in axion string networks
Mark Hindmarsh, Joanes Lizarraga, Asier Lopez-Eiguren, Jon Urrestilla
HHIP-2021-7/TH
Approach to scaling in axion string networks
Mark Hindmarsh,
1, 2, ∗ Joanes Lizarraga, † Asier Lopez-Eiguren,
1, 4, ‡ and Jon Urrestilla § Department of Physics and Helsinki Institute of Physics, PL 64, FI-00014 University of Helsinki, Finland Department of Physics and Astronomy, University of Sussex, Falmer, Brighton BN1 9QH, U.K. Department of Physics, University of the Basque Country UPV/EHU, 48080 Bilbao, Spain Institute of Cosmology, Department of Physics and Astronomy, Tufts University, Medford, MA 02155, USA (Dated: February 16, 2021)We study the approach to scaling in axion string networks in the radiation era, through measuringthe root-mean-square velocity v as well as the scaled mean string separation x . We find good evidencefor a fixed point in the phase-space analysis in the variables ( x, v ), providing a strong indicationthat standard scaling is taking place. We show that the approach to scaling can be well describedby a two parameter velocity-one-scale (VOS) model, and show that the values of the parameters areinsensitive to the initial state of the network. The string length has also been commonly expressedin terms of a dimensionless string length density ζ , proportional to the number of Hubble lengthsof string per Hubble volume. In simulations with initial conditions far from the fixed point ζ is stillevolving after half a light-crossing time, which has been interpreted in the literature as a long-termlogarithmic growth. We show that all our simulations, even those starting far from the fixed point,are accounted for by a VOS model with an asymptote of ζ ∗ = 1 . ± .
09 (calculated from the stringlength in the cosmic rest frame) and v ∗ = 0 . ± . I. INTRODUCTION
The axion [1, 2] is a hypothetical pseudoscalar parti-cle predicted in the Peccei-Quinn (PQ) mechanism [3, 4].The PQ mechanism extends the Standard Model (SM) ofparticle physics to solve the so-called strong-CP problemof Quantum Chromodynamics by adding an extra spon-taneously broken global U(1) symmetry, which is anoma-lous. If the symmetry is broken at a high scale [5–8], theparticle is very weakly interacting and long-lived, andbecomes a well-motivated dark matter candidate [9–11].The PQ symmetry can be spontaneously broken be-fore or after inflation, leading to very different axionicdark matter scenarios. If the symmetry is broken beforeinflation, the Universe is filled by a homogeneous axionfield, which produces zero-momentum axions through thevacuum realignment mechanism [9–11]. On the otherhand, in the post-inflationary PQ symmetry breaking,the phase transition happens in the standard cosmol-ogy, producing axionic cosmic strings [12, 13]. These de-fects are a variety of global cosmic string [14, 15], whichlive until the QCD confinement transition, when theyform hybrid string-wall composites and are annihilated[12, 16, 17]. The strings formed in this last scenario arethe ones studied in this paper.Axion strings release energy mainly into pseudo-Goldstone radiation, both during their evolution as wellas in the string-wall system collapse. This radiation con-stitutes an initially degenerate gas of axions with a non-thermal distribution. The complicated non-linear dy-namics at the QCD transition imprints density fluctu- ∗ mark.hindmarsh@helsinki.fi † [email protected] ‡ asier.lopez [email protected] § [email protected] ations which provide the seed for the axion miniclusterformation through early gravitational collapse [18–21].These miniclusters could be detected by their distinctivesmall-scale lensing signals [18, 21–23].The evolution of axion strings and axion radiation isgoverned by the classical field equations of the underlyingscalar field theory, and due to their non-linearities, latticesimulations are required to go beyond order-of-magnitudeestimates. In recent years, several groups have studiedthe evolution of axion strings and the production of ax-ions using lattice simulations [24–42].An accurate calculation of the total number densityof axions produced is essential for an accurate calcula-tion of the axion energy density, which if matched to thedark matter density today, gives a prediction for the ax-ion mass. While the axion number density is not verysensitive to the string density at the QCD transition (see[43] for a discussion), high accuracy in the mass estimateis required for resonant cavity searches, which are cur-rently targeting axion masses appropriate for productionby vacuum realignment [44].The prediction of the axion density depends on havingan accurate description of axion string evolution. As thestring evolution takes place from the PQ transition ataround 10 GeV, to the QCD transition at 100 MeV (afactor 10 ), the results from numerical simulations mustbe extrapolated.It is important to have a physical basis for the extrapo-lation. Such a physical model is the one-scale model [45]and its velocity-dependent improvement [46–48], whichwe discuss in more detail below. It predicts that thestring network should approach a scaling solution, wherethe mean string separation grows in proportion to cosmictime t , and the RMS velocity of the strings is constant.By scaling we mean that at distances much larger thanthe string width, network length scales such as the meanstring separation are proportional to the cosmic time t . a r X i v : . [ a s t r o - ph . C O ] F e b In the standard scaling picture the dynamical evolutionof the network is independent of the string width andtension. The justification is based on approximating thestring dynamics by the Nambu-Goto equations of motion,in which the string tension drops out and the string widthplays no role.The general picture of network evolution is that stringsare initially in a dense tangle of loops and infinite strings,with mean separation set by the correlation length of thefield and the cooling rate [49, 50]. They decay by thecollapse of loops of string initially present, and choppedoff from the infinite strings. The mean separation grows,until it becomes of order t , the cosmic time. This isapproximately as fast as causality allows. By this timethe temperature has dropped many orders of magnitude,and friction with the cosmic plasma is negligible. Stringsevolve essentially in vacuum, in the background providedby the the rest of the matter in the universe.In a previous paper [40], we presented results on thescaling dynamics of axion strings. We showed that thenetwork evolution was consistent with standard scaling,and obtained an asymptotic value for the dimensionlesslength density of axion string ζ , which is proportional tothe mean number of Hubble lengths of string per Hubblevolume, ζ ∞ = 1 . ± .
20. This is consistent with, andimproves in accuracy, estimates in earlier works [24–30,32]. Equivalently, the mean string separation ξ is alwaysabout half a horizon length.In [40] we also showed how the presence of the scalingsolution can be disguised, either by the choice of vari-able to study, or by the choice of initial conditions. Inthis light, claims of a slow or logarithmic growth in thedimensionless length density [35–39, 41, 48] are to beinterpreted as an approach to scaling from initially lowvalues of ζ that it is not completed before the simulationends.It is important to note that the only significant dif-ference between the simulations of the different groupsis the method for preparing the initial conditions of thefield, which determines the initial string separation, andthat there are no significant differences in the subsequentevolution of the string network. All but one group report ζ (cid:46) l φ in order to cover a range of initial string sep-arations, which tests the sensitivity of the system to theinitial conditions. The evolution of the system has beencarried out using both the true physical field equations,and also using the Press-Ryden-Spergel (PRS) method[51] to simulate strings with constant comoving width.Simulating (a priori unphysical) strings with constant co-moving width allows for a longer period of scaling, thusgiving insight into the long-term behaviour of a systemof strings. We extend the study in [40] by analysing the root-mean-square velocity v of the networks alongside themean string separation in units of cosmic time, x = ξ/t .We demonstrate that the evolution of the simulationsat later stages of the simulation can be well describedby a two-parameter velocity-dependent one-scale (VOS)model [46–48] where all simulations tend asymptoticallyto a common point in the phase space ( x ∗ , v ∗ ), a fixedpoint of the VOS dynamical system. The dynamical sys-tems analysis predicts that the approach to the fixedpoint is governed by a pair of complex exponents withnegative real parts, a stable spiral. We find good quan-titative accord with the prediction near the fixed point,where the model is supposed to be a good description.Due to this good accord, we obtain a more precise es-timate of the scaling density of strings than in our previ-ous analysis [40]. We find that the physical and constantcomoving width systems have fixed points which are con-sistent with each other.Further away from the fixed point, the qualitativeagreement is good. The VOS model predicts that ini-tially overdense ( x < x ∗ ) networks will accelerate, andevolve towards an underdense ( x > x ∗ ) network as theenergy-loss mechanism (production and decay of stringloops) overcompensates. The approach to scaling fromthe underdense side is a common feature of simulations.The model also predicts that very underdense networkstake a long time to reach scaling, often longer than half-box crossing time, consistent with the very underdensesimulations of Refs. [35, 41]. II. MODEL AND NETWORK PARAMETERSA. Field dynamics
The simplest axion models include a singlet scalar fieldwith a U(1) symmetry, Φ, with action S = (cid:90) d x √− g (cid:16) ∂ µ Φ ∂ µ Φ − λ (Φ − η ) (cid:17) , (1)where we have written the field as a two-component vec-tor, and the U(1) symmetry is realised as a rotation onthe vector.In a FLRW metric, and when the field is coupled to athermal bath of weakly-coupled particles, the equationsof motion take the formΦ (cid:48)(cid:48) + 2 a (cid:48) a Φ (cid:48) − ∇ Φ = − a λ (Φ − η ( T ))Φ , (2)where a is the scale factor, a prime denotes differentiationwith respect to conformal time τ , and in the radiation era a ∝ τ . The free energy of the system is minimised at thefield magnitude η ( T ), where η ( T ) = d ( T − T ), T c (cid:39) η is the critical temperature of the PQ phase transition,and d is a constant computable in perturbation theory.For T (cid:29) T c , it is energetically favourable for the field tofluctuate around Φ = 0. Well below the critical temper-ature it is energetically favourable for the magnitude ofthe field to take the value η , with a massless pseudoscalarfluctuation mode (the axion) and a scalar mode of mass m sca = √ λη . During the phase transition, the direc-tion in field space is chosen at random in uncorrelatedregions of the universe, with the result that the field isforced to stay zero along lines [49]. These lines form thecores of the axion strings [13]. The size of the core isapproximately w = m − . B. Network parameters from field averages
The subsequent evolution of the string network can betracked by the string length (cid:96) and the RMS velocity v ofthe strings.A couple of estimators for (cid:96) are possible. We define thewinding length (cid:96) w as the number of plaquettes piercedby strings multiplied by the physical lattice spacing aδx ,corrected by factor of 2 / E = E π + E D + E V (3)with the functions E π = 12 (cid:90) d x Π W (Φ) , (4) E D = 12 (cid:90) d x ( ∇ Φ) W (Φ) , (5) E V = (cid:90) d xV (Φ) W (Φ) , (6)where Π = ( ∂ t Φ) and V (Φ) = λ (Φ − η ) . The function W (Φ) is a local function of the fields which is stronglypeaked near Φ = 0, and zero for | Φ | = η , so that itpicks out strings. We call the three functions definedabove the weighted kinetic, gradient and potential energyrespectively.Suppose that all the energy in the volume V is in theform of global strings, centered on the line X ( σ, t ). Thecoordinate σ is chosen so that | X (cid:48) | = (1 − ˙ X ) , wherethe prime represents the derivative with respect to σ , andthe dot the derivative with respect to t . We denote thetotal rest-frame length of string (cid:96) r = (cid:90) dσ. (7) Writing local rest frame space coordinates x s , and fieldsmeasured in the local rest frame with the subscript s, thefields of a piece of string moving with orthogonal velocity˙ X are Π( x , t ) = γ ˙ X · ∇ Φ s ( x s ) , (8) ∇ Φ( x , t ) = γ ˆ v (ˆ v · ∇ Φ s ( x s )) + ∇ ⊥ Φ( x , t ) , (9)where ˆ v is a unit vector in the direction of ˙ X , γ =1 / (cid:112) − ˙ X is the boost factor, and ∇ ⊥ i Φ( x , t ) = ( δ ij − ˆ v i ˆ v j ) ∇ j Φ( x , t ) . (10)Choosing the local rest frame so that the string is orientedin the z s direction, a string moving with velocity ˙ X hasscalar kinetic energy E π = 14 (cid:90) dx s dy s ( ∇ Φ s ) W (Φ s ) (cid:90) dσ ˙ X , (11)gradient energy E D = 14 (cid:90) dx s dy s ( ∇ Φ s ) W (Φ s ) (cid:90) dσ (cid:18) γ (cid:19) , (12)and potential energy E V = (cid:90) dx s dy s V (Φ s ) W (Φ s ) (cid:90) dσ γ . (13)The total energy is therefore E = µ (1 − f V v ) (cid:96) r , (14)where µ = (cid:90) dx s dy s (cid:20)
12 ( ∇ Φ s ) W (Φ s ) + V (Φ s ) W (Φ s ) (cid:21) (15)is the W -weighted mass per unit length of a static string,with f V the fraction contributed by the potential energydensity, and we have defined an RMS velocity v through v = 1 (cid:96) r (cid:90) dσ ˙ X . (16)A convenient choice for the weight function is W = V (Φ) , (17)for which µ = 0 . η and f V = 0 . Besides the energy, we can also calculate the La-grangian L = E π − E D − E V , (18)(19) These numbers are obtained from a code implementing a relax-ation method on the discretised radial energy functional, with800 lattice points and lattice spacing ηdr = 0 .
01. The conver-gence criterion was that the change in energy in an update shouldbe less than 10 − η . finding L = − µ (1 − v ) (cid:96) r . (20)Combining the total energy E and the Lagrangian esti-mators, an estimate for both the rest-frame length (cid:96) r andthe mean square velocity can be obtained (cid:96) r = E + f V Lµ (1 − f V ) , (21) v L = E + LE + f V L , (22)where the subscript L denotes the use of the Lagrangianto obtain the estimate. An alternative way of estimatingthe string velocity comes from the pressure, p V = E π − E D − E V , (23)which depends on the rest frame length and RMS velocityas p V = 13 µ(cid:96) r (cid:2) (2 v − − f V (2 − v ) (cid:3) . (24)It is then straightforward to derive another mean squarevelocity estimator v ω = 1 + 3 ω + 2 f V f V (1 + 3 ω ) , (25)where ω = p V /E is the equation of state parameter ofthe strings. A third estimate for the string velocity canbe constructed from the ratio of the kinetic to gradientenergies [53], R s = E π E D , (26)which can be rearranged to give v = 2 R s R s . (27)Given that we only have three independent underlyingquantities E π , E D , and E V , only three independent es-timators can be derived from them: one length, and twovelocity estimators.The winding length is not derived from the weightedenergies, and so is an independent length estimator. Asit is the ordinary Euclidean length of the curve traced bythe string, it can be represented as (cid:96) w = (cid:90) dσ | X (cid:48) | = (cid:96) r (cid:10) γ − (cid:11) . (28)Note that the average of γ − is not in general equal to(1 − v ) / .In a cosmological simulations one can express thestring length in terms of Hubble lengths per Hubble vol-ume, or ζ = (cid:96)t V . (29) When investigating scaling in string networks, it is moretransparent to parametrise the string density by themean string separation, which is obtained from measuresof the string length via ξ = (cid:114) V (cid:96) . (30)In this work we will use two length estimators, whichwill define two different mean string separation estima-tors: when the length estimator used is the rest-frameestimator (cid:96) r , we will define ξ r ; and when the length es-timator used is the winding length estimator (cid:96) w , we willdefine ξ w .The above estimators were derived for a Minkowskispace-time. In an expanding background, one can viewthe space-time coordinates as representing comoving po-sition and conformal time, from which physical lengthsfollow by multiplication by the scale factor a . III. SIMULATIONS AND SCALINGOBSERVABLE RESULTS
We solve a discretised version of the equations of mo-tion (2) on a cubic lattice with periodic boundary condi-tions, evolving the system in conformal time t The resultswe present in this section are extracted from the same setof simulations analysed in [40], where lattices with 4096sites per dimension were used with spatial resolution of δxη = 0 . δτ η = 0 .
1. Inaddition, we use a set of simulations with a larger initialcorrelation length, but otherwise identical. In the follow-ing lines we will only summarise the procedure and referthe reader to [32, 40] for more detailed descriptions onthe method.The field configuration is initiated at conformal time τ start by setting the scalar field canonical momentum tozero and the scalar field to be a Gaussian random fieldwith power spectrum P Φ ( k ) = A (cid:2) kl φ ) (cid:3) − , were A is chosen so that (cid:104) Φ (cid:105) = η and l φ is the field correlationlength in comoving coordinates. We use different valuesof l φ in order to cover a range of string separations inthe initial conditions. In order to allow the strings toform, and to remove the energy excess in the field fluc-tuations around the string configurations, we evolve thisconfiguration with a diffusion equation with unit diffu-sion constant until conformal time τ diff . We then applythe second-order time evolution equation (2).Similarly to our previous paper, we extract data fromsimulations with both fixed comoving string width andfixed physical string width. We promote the scalar self-coupling constant to be a time dependent parameter λ = λ /a − s ) following the PRS method [51]. This makesthe comoving string width decrease with conformal timeas: w ( τ ) = w a s ( τ ) . (31)The physical equation of motion, where the physicalstring width remains constant at w = 1 / √ λ η , andthe comoving width decreases with time, corresponds to s = 1. With s = 0 the comoving width is constant at w and the physical string width increases in time.For the s = 1 case, it is difficult to avoid the stringwidth being larger than the Hubble length at early times,which also means that the relaxation of the field to itsequilibrium value is longer than a Hubble time. In orderto speed up the string formation, we arrange the time-dependence of the coupling so that strings are formedand diffused with a constant comoving width, equal totheir final comoving width. At the end of the diffusionperiod, the string width is much smaller than its physicalvalue w . The string width is then allowed to grow bysetting s = − τ cg , which is when the string corehas expanded to its correct physical width w . Afterconformal time τ cg , the physical evolution with s = 1starts. We call this procedure core growth. Simulationsend at conformal time τ end .Table I contains all simulation parameter choices thathave been considered in the procedures described above.Four simulations with different random number seedswere carried out at each parameter choice, for a to-tal of 28 runs. The data are analysed in cosmic time t = ( τ /τ end ) τ end / Model s = 1 s = 0 l φ η (5,10,20,40) (5,10,20) τ start η
50 50 τ diff η
70 70 s cg -1 – τ cg η τ end η Figure 1 shows the comparison of the evolution of themean string separation for ξ r and ξ w presented in theprevious section (30) for a single run with l φ η = 5. Theupper panel is for simulations with s = 1 and the lowerpanel for s = 0. Their growth is consistent with a linearasymptote, as extensively studied in [40]. This is the ex-pectation from the standard picture of scaling in axionstring networks [12, 45, 46]. Note that in [40], the wind-ing length estimator ξ w was used to establish the asymp-totic linear growth; here we establish that the rest-frameestimator also grows linearly, as expected.The ratio of the winding length estimator (cid:96) w to the restframe estimator (cid:96) r is plotted in Fig. 2 (solid line), whichaccording to Eq. (28) is an estimate of (cid:10) γ − (cid:11) , where γ isthe Lorentz factor of the string. For comparison, we plot(1 − v ) / , using the scalar field estimator (27), whosetime-dependence in the simulations is discussed later inthis section. As pointed out in the previous section, the tη ξ η ξ w ξ r ξ w (1 − v ) / tη ξ η ξ w ξ r ξ w (1 − v ) / FIG. 1. Comparison of the mean string separation definedfrom the winding length estimator ξ w (solid black) and therest-frame estimator ξ r (solid blue) as presented after Eq. (30),from a single simulations with correlation length l φ η = 5 and s = 1 (upper panel) and s = 0 (lower panel). The dashedblack line corresponds to the winding estimator ξ w modifiedvia Eq. (28). two quantities are not necessarily equal, but empiricallywe observe that they are close by the end of the simula-tion.The closeness of (1 − v ) / to (cid:10) γ − (cid:11) is also observedin Fig. 1, where we show as a dashed line the windingestimator multiplied by (1 − v ) / . We choose (cid:96) r as thelength estimator for the rest of this work, which is bettersuited to the dynamical modelling we carry out. It canbe related to the winding length through the factor ofapproximately 0 . σ standard deviations. Uncertain-ties are calculated by propagating the fluctuations in theweighted energies (11), (12) and (13). We use black for l φ η = 5, red for l φ η = 10, blue for l φ η = 20, and green tη . . . . . . ‘ w / ‘ r s = 1 ‘ φ η = 5 s = 1 ‘ φ η = 10 s = 1 ‘ φ η = 20 s = 1 ‘ φ η = 40 tη . . . . . . ‘ w / ‘ r s = 0 ‘ φ η = 5 s = 0 ‘ φ η = 10 s = 0 ‘ φ η = 20 FIG. 2. Ratios of the winding length estimator (cid:96) w to rest-frame estimator (cid:96) r (solid line) for different correlation lengths.The lines correspond to a single simulations. In dashed, weplot (1 − v ) / . for l φ η = 40 (only in simulations with s = 1). The endof the core growth period is shown as a vertical greendashed line. These figures extend the results of Fig. 1 inRef. [40], which shows the winding length estimator onlyfor s = 1, and a subset of the initial correlation lengths.We now turn to the velocity estimators. To establishtheir consistency, we plot all three for the same run inFig. 4, with s = 1 in the top panel and s = 0 in the bot-tom panel. As mentioned in the previous section, onlytwo are independent, but the fact that all three are soclose gives confidence that they are indeed estimatinga global translational velocity of a string-like solution,rather than field fluctuations in regions where Φ is closeto zero, which is a potential contaminant of velocity es-timators.Uncertainties in velocities are calculated by propagat-ing the fluctuations in the weighted energies (11), (12)and (13). We find that the largest fluctuations are in theweighted potential energy E V . We therefore choose the tη ξ r η s = 1 ‘ φ η = 5 s = 1 ‘ φ η = 10 s = 1 ‘ φ η = 20 s = 1 ‘ φ η = 40 tη ξ r η s = 0 ‘ φ η = 5 s = 0 ‘ φ η = 10 s = 0 ‘ φ η = 20 FIG. 3. Mean string separation ξ r from simulations with s = 1 (top panel) and s = 0 (bottom panel) with initial fieldcorrelation lengths l φ η = 5 (black), l φ η = 10 (red), l φ η = 20(blue) and l φ η = 40 (green - only for s = 1). The solid linerepresents the mean over realisations of ξ at each time, withthe shaded regions showing the 1- σ variation. The verticalgreen line corresponds to the end of the core growth period,after which the system is evolved with the physical equationsof motion in the s = 1 case. estimator v s derived from kinetic and gradient energiesonly (27) as the mean square string velocity estimator,and show the means and uncertainties in Fig. 5.We see that after an initial period of acceleration, thereis a decreasing trend, approaching what appears to be aconstant value at the end of the simulations. The max-imum velocity is larger for the fields with smaller cor-relation length in the initial conditions, as is consistentwith string-like behaviour, where acceleration is propor-tional to curvature. For the case of strings with constantphysical width ( s = 1), the RMS velocity is approxi-mately constant during the core growth phase, and thenapproaches an asymptote more slowly than the s = 0simulations.Figures 3 and 5 show that, independently of the initial tη . . . . . . v v s v ω v L tη . . . . . . v v s v ω v L FIG. 4. Comparison of velocity estimators presented inSec. II B for a simulation with correlation length l φ η = 5 and s = 1 (upper panel) and s = 0 (lower panel). The valuesof the scalar field velocity estimator v s (27), the equation ofstate velocity estimator v ω (25) and the Lagrangian-derivedvelocity estimator v L (23) are shown in black, red and blue,respectively. field correlation length, all simulations are compatible, i.e. all of them give separation and velocity data whichare within 1 σ of each other. Moreover, the behaviourof both estimators ( ξ and v ) qualitatively agrees withthe standard scaling, showing a tendency towards lineargrowth in ξ and a constant RMS velocity.There is a departure from standard scaling in the ear-lier phases of the simulations, which needs to be under-stood in order to improve the estimates of the asymptoticbehaviour of v and ξ , or more precisely the asymptoticvalues of the scaled mean string separation, x = ξ/t. (32)In our previous paper [40], which studied scaling using ξ w only, we observed that the average slope ∆ ξ/ ∆ t con-verged more quickly to a constant than ξ/t . We therefore tη . . . . . . v s s = 1 ‘ φ η = 5 s = 1 ‘ φ η = 10 s = 1 ‘ φ η = 20 s = 1 ‘ φ η = 40 tη . . . . . . v s s = 0 ‘ φ η = 5 s = 0 ‘ φ η = 10 s = 0 ‘ φ η = 20 FIG. 5. Velocities from scalar field estimator v s (27) fromsimulations with s = 1 (top panel) and s = 0 (bottom panel)with initial field correlation lengths l φ η = 5 (black), l φ η = 10(red), l φ η = 20 (blue) and l φ η = 40 (green - only for s = 1).The solid line represents the mean over realisations at eachtime, with the shaded regions showing the 1- σ variation. Thevertical green line corresponds to the end of the core growthperiod. used the slope of the curve of ξ against t as our estimatorof the asymptotic value x ∗ . The linear fit can have a sig-nificant constant term, which we parameterised in termsof the intercept with the time axis, the time offset, t .The value of t has no physical importance, and insteadparametrises an effect of the initial conditions.In this paper we make use of RMS velocity data, whichgives extra information about the approach to scaling,and avoids the need for t . We will see that in doingso we improve the accuracy of the estimate of x ∗ , whileremaining consistent with our previous estimate.In Fig. 6 we show the evolution of the network in thephase space ( x, v ), where the scaled mean string sepa-ration is measured with the rest frame string length x r and the RMS velocity is measured using the scalar fieldenergies v s . . . . . . . . x r . . . . . v s . . . . . . . x r . . . . . v s FIG. 6. Phase space plot for v s and x r for s = 1 (top panel)and s = 0 (bottom panel), and the same colour scheme corre-spond to different initial correlation lengths, as the previousfigure. Larger dots are plotted every 20 cosmic time units,starting at cosmic time tη = 4. As time increases, all curvesspiral in towards an apparent fixed point, analysed in SectionV. The phase space representation shows clearly the dif-ferent regimes in which the network evolves. At the endof the diffusion period the strings are accelerated undertheir curvature, and the RMS velocity increases rapidlywhile the inter-string distance remains nearly constant.For s = 1 simulations, the diffusive evolution is followedby the core growth period, which is part of the prepa-ration of the initial conditions. During the core growthperiod, velocities remain approximately constant. Thescaled mean string separation, however, changes, and itchanges differently for different initial correlation lengths:For correlation lengths l φ η = 5 ,
10, the scaled meanstring separation grows, whereas for correlation lengths l φ η = 20 ,
40 it decreases. Finally, when the physicalequations of motion (2) are being solved, the systemstarts to spiral towards an apparent common fixed pointfor all simulations. Estimating the position of the fixed point and hence the asymptotic values of x and v is thesubject of the next two sections.It is interesting to note the qualitative difference inthe velocity evolution between the core growth era andthe physical equations of motion. In the core growth erathe velocity remains constant after the initial accelera-tion: this constant depends on the initial conditions. Itis only after the physical evolution sets in that the veloc-ity starts evolving towards its asymptote. Note that thecore growth era corresponds to evolution with the fixedscale hierarchy m sca /H explored in Ref. [39]. We willdiscuss this observation in the final section.The initial conditions of the field and the time at whichthey are set determine the simulation’s starting point inthe phase space. Figure 6 shows that varying the initialfield correlation length one can choose whether to starton the left hand side or on the right hand side of thehypothetical fixed points, corresponding to strings beingeither above or below their scaling density. As mentionedbefore, the time offset t depends on the the initial condi-tion, and for approaches to the fixed point from the rightcorresponds to t < IV. PHASE SPACE ANALYSIS WITH THE VOSMODEL
In this section we model our results as a dynami-cal system. The model best adapted to a network ofstrings is the velocity-dependent one-scale (VOS) model[45, 46, 48]. This class of models assumes a statisticaldistribution of string configurations and velocities whichhas a universal form, parametrised by the string separa-tion ξ (or equivalently the length (cid:96) in a volume V ) andRMS velocity v .When applied to Nambu-Goto strings, the VOS modeldescribes “long” strings only, that is, either infinitestrings, or string loops with total length greater thansome threshold of order ξ . In our simulations, stringlengths and velocities are measured over the whole stringnetwork, including loops. As a string network’s totallength is dominated by strings winding around the sim-ulation’s periodic box, the distinction should not be im-portant for a first approximation.The movement of the long strings results causes a seg-ment of string of length ξ to encounter others at a rate oforder v/ξ . The encounter causes the string to reconnect,producing a loop. If this loop is smaller than ξ it radiatesand shrinks without further encounters, apart from self-intersections. The net result is the loss of energy from thestring network into axions and massive scalar radiation.The string motion is a balance between the accelera-tion caused by the curvature, and the Hubble damping.There can also be damping due to the preferential lossof energy from fast-moving segments of string, which en-counter others more rapidly: this effect is neglected in thesimplest models. There is also direct energy loss fromlong string in the form of radiation, which we will notdistinguish from energy loss via loops in our modelling.The equations of motion for a system of Nambu-Gotostrings, and the assumptions above lead to the followingdynamical system, dξdt = Hξ (1 + v ) + 12 cv, (33) dvdt = (1 − v ) (cid:18) kξ − Hv (cid:19) , (34)where H is the Hubble parameter. The model has twophenomenological parameters, k and c . The parameter c describes the efficiency of the energy loss mechanism,while the parameter k describes the correlation betweenthe string curvature vector and the velocity. For exactlyNambu-Goto strings, k is a function of velocity. However,as a first approximation, and because the strings we arestudying are not exactly Nambu-Goto, we will take k tobe a constant.Using the dimensionless mean string separation vari-able x (32), and taking, H = 1 / t as appropriate for aradiation-dominated universe, t ˙ x = 12 x (cid:0) v − (cid:1) + c v (35) t ˙ v = (cid:0) − v (cid:1) (cid:18) kx − v (cid:19) (36)This dynamical system has a fixed point in the relevantregion 0 ≤ x , 0 ≤ v < x ∗ = (cid:112) k ( c + k ) v ∗ = (cid:115) k ( c + k ) . (37)From here, one can express the parameters c and k interms of the fixed point values, k = x ∗ v ∗ , c = x ∗ v ∗ (1 − v ∗ ) . (38)Small perturbations ( δx, δv ) evolve to the fixed point ac-cording to t ddt (cid:18) δxδv (cid:19) = M ∗ (cid:18) δxδv (cid:19) , (39)where M ∗ = (cid:18) ( v ∗ − x ∗ v ∗ (1 + v ∗ ) v ∗ x ∗ ( v ∗ − v ∗ − (cid:19) (40)The eigenvalues of the matrix σ ± are σ ± = − κ ± (cid:114) κ − κ where κ = 1 − v ∗ . Since 0 < κ < σ ± are complex, with negative real part. Therefore, the fixedpoint is a stable spiral.We plot flows in the phase diagram predicted by theVOS model in Fig. 7 for the global best fit ( x ∗ , v ∗ ), alongwith the mean values of selected ( x, v ) from the simu-lations. The stable spiral form is clearly visible in thestreamlines. In the next section we explain the fittingprocedure from which the global best fit ( x ∗ , v ∗ ) was ob-tained. V. FITS AND ASYMPTOTIC BEHAVIOUR
In this section we measure the degree at which theevolution dictated by the VOS model presented in theprevious section, Eqs. (33) and (34), is compatible withour simulations, by performing a fitting analysis. Wecompute the χ value for each set of runs with a given( s, l φ ) as χ = (cid:88) i ( O i − E i ) σ i , (41)where O i is the observed value, E i the expected valueon the basis of the model, and σ i the uncertainty in theobserved value. In our case, the observed values are thetime series data ( x, v ) recorded from our simulations. Weuse a bootstrapping method to create the time series for aspecific l φ . For each case we have four different runs, outof which we create the bootstrapped time series by choos-ing randomly a value at specific time step i . This proce-dure is performed four times so that four different boot-strapped time series are created. The observed valuesand their uncertainties are then obtained by averagingand computing the standard deviation from those boot-strapped realisations. The expected value is the valuepredicted by the VOS model, as described below. Theset of observations is taken in the time range [ t fit , t end ],where the start of the fitting period is t fit η = 171 . τ fit η = 600). Note that the l φ η = 40 case ispresent only for s = 1.We explore the two dimensional parameter space usinga grid of size 100 × . < k < . . < c <
1. These are set by preliminary analysisof a wider parameter space. s l φ η k c x ∗ v ∗ . ± .
006 0 . ± .
010 0 . ± .
009 0 . ± . . ± .
003 0 . ± .
017 0 . ± .
007 0 . ± . . ± .
010 0 . ± .
010 0 . ± .
010 0 . ± . . ± .
013 0 . ± .
032 0 . ± .
012 0 . ± . . ± .
008 0 . ± .
016 0 . ± .
013 0 . ± . . ± .
011 0 . ± .
032 0 . ± .
021 0 . ± . . ± .
004 0 . ± .
020 0 . ± .
006 0 . ± . . ± .
005 0 . ± .
015 0 . ± .
006 0 . ± . . ± .
027 0 . ± .
039 0 . ± .
037 0 . ± . TABLE II. Inferred best-fit values of model parameters c and k , and asymptotic values of x and v for each correlation lengthin s = 0 and s = 1. These values are obtained by fitting a setof 20 bootstrap realisations for each correlation length, usingdata with t > . s = 0 and s = 1 are the standarddeviations of the mean values for each correlation length. The best-fit values of c and k , and asymptotic scaledmean string separation and velocity ( x ∗ and v ∗ ) for each0( s, l φ ) can be found in Table II, where fits were takenwith t fit η = 171 .
4. The final mean values for s = 0 , l φ , and thequoted uncertainties correspond to the resulting standarddeviations. The errors obtained by quadrature combina-tion of bootstrap errors were systematically smaller thanthe standard deviations. Mean values for the parametersfor other values of t fit η = 171 . . . . . . . . x r . . . . . v s . . . . . . . x r . . . . . v s FIG. 7. Phase plane for s = 1 (top), s = 0 (bottom),with stream lines of the best-fit values of the VOS modelparameters shown in Table III, with t fit η = 50. The colourscheme is the same as in previous figures. The larger markerscorrespond to the same points as in Fig. 6, with empty circlesdenoting points with t fit η < Figure 7 shows the evolution of the simulations along-side the streamlines of the VOS dynamical system cal-culated using the inferred global mean values ( x ∗ , v ∗ ),obtained by fitting with t fit η = 50, the earliest time fromwhich we start fitting (see Table III). It can be seen thatafter an initial relaxation period, the simulation data fol-low the spiral-like evolution towards the fixed point. Akey feature is that initial conditions with x < x ∗ tend toflow to states with v > v ∗ , and then around to x > x ∗ , corresponding to string networks less dense than scaling.For a more quantitative comparison between our simu-lations and the VOS model, we show in Figure 8 the rel-ative difference between the simulation time series dataand the VOS best-fit model for each ( s, l φ ), where the ini-tial conditions for the integration of the VOS equationsare set at t fit . Shaded regions correspond to the uncer-tainties propagated from simulation estimators. It canbe observed that the mean relative difference always liesbelow 5% level of deviation, with zero deviation alwayswithin the errors. s t fit η k c x ∗ v ∗ . ± .
027 0 . ± .
008 0 . ± .
030 0 . ± . . ± .
018 0 . ± .
018 0 . ± .
022 0 . ± . . ± .
017 0 . ± .
034 0 . ± .
020 0 . ± . . ± .
012 0 . ± .
041 0 . ± .
013 0 . ± . . ± .
011 0 . ± .
053 0 . ± .
026 0 . ± . . ± .
015 0 . ± .
077 0 . ± .
025 0 . ± . . ± .
036 0 . ± .
115 0 . ± .
029 0 . ± . . ± .
057 0 . ± .
013 0 . ± .
062 0 . ± . . ± .
043 0 . ± .
015 0 . ± .
045 0 . ± . . ± .
031 0 . ± .
021 0 . ± .
036 0 . ± . . ± .
030 0 . ± .
039 0 . ± .
037 0 . ± . . ± .
028 0 . ± .
046 0 . ± .
037 0 . ± . . ± .
028 0 . ± .
083 0 . ± .
032 0 . ± . . ± .
060 0 . ± .
102 0 . ± .
050 0 . ± . TABLE III. Global mean values of model parameters c and k , and asymptotic values of x and v , with different start timesfor the fit t fit . Table III contains the global mean values obtained ap-plying the same bootstrap procedure as in the previousanalysis. The last four fit start times are those used forthe linear fitting procedure carried out in Ref. [40]. A re-markably good agreement is obtained in the global meanswhen comparing early and late fits. Earlier fits have asmaller scatter in c , but a larger scatter in k . This canbe related to the spread of velocities in the initial condi-tions The spread is minimised for the intermediate times t fit = 150 ,
171 and 233. We quote results from t fit = 171,which is also the earliest fitting time from our previouspaper [40], meaning that a direct comparison of the meth-ods can be made.The simplest VOS model therefore gives a good quanti-tative description of the joint evolution of the string sep-aration and RMS velocity. We have experimented withfitting for additional parameters q and d (with β = 1 and r = 1) in the VOS model presented in Ref. [54]. but ourpreliminary analysis shows that the preferred values forthese additional parameters are compatible with zero.It is also interesting to study the network evolutionin terms of the length density parameter ζ (29), and by1
200 250 300 350 400 450 500 tη − . − . . . . ( x V O S − x ) / x V O S
200 250 300 350 400 450 500 tη − . − . . . . ( v V O S − v ) / v V O S FIG. 8. Relative difference of the VOS prediction and thesimulation data of the dimensionless string separation x andrms velocity v . The shaded bands represent errors propagatedfrom simulation’s estimators. plotting against the logarithm of time one can emphasisethe earlier times when the network is further away fromscaling. Fig. 9 shows our s = 1 rest-frame length dataplotted this way, along with the best-fit VOS models foreach correlation length, and their extrapolation to largervalues of time. The asymptotic ζ r , ∗ obtained from theoverall mean values of fit parameters in Table II is alsodepicted, for which we obtain ζ r , ∗ = 1 . ± .
11. Thecentral value is shown as solid purple line and its corre-sponding errors in shaded purple bands. Note that allsimulations approach the asymptotic ζ from below, andare still slowly increasing at the end of the simulation,but within 20% of its asymptotic value. The increaseis most noticeable for simulations which start very un-derdense, and therefore have further to evolve to reachscaling.We also include the value of the asymptotic lengthdensity parameter reported in our previous paper [40]in dashed black with its corresponding uncertainty ingrey. In [40] the universe-frame string length (cid:96) w wasused as the measure of the string length, and the valueof x ∗ estimated by linear fits to ξ w against t , fromwhich we obtained an estimate of the asymptotic value ζ ∗ = 1 . ± .
20. In this work we use the rest-framelength estimator, the corresponding string density pa-rameter is (28), for which ζ r , ∗ = ζ ∗ (1 − v ∗ ) − / , and gives ζ r , ∗ = 1 . ± .
25. As the figure shows, the agreementwith our previous result is very good. log( tη ) . . . . . . . . . ζ r s = 1 ‘ φ η = 5 s = 1 ‘ φ η = 10 s = 1 ‘ φ η = 20 s = 1 ‘ φ η = 40 log( tη ) . . . . . . v s s = 1 ‘ φ η = 5 s = 1 ‘ φ η = 10 s = 1 ‘ φ η = 20 s = 1 ‘ φ η = 40 FIG. 9. String network evolution expressed in terms of therest-frame length density parameter ζ and RMS velocity v ,plotted against log( tη ), with VOS models (dotted line) thatcorrespond to best fit values for each l φ shown in Table II.The prediction of the VOS model is plotted from t fit on. Thehorizontal dashed black line and grey band in the top panelshow the mean and uncertainly obtained from our previousanalysis [40], translated from the universe-frame length usedin that paper by multiplying by (1 − v ∗ ) − / . The solid purpleline and shaded purple bands are the mean and uncertainlyobtained from the analysis in this work. The vertical greendashed line marks the end of the core growth phase, and thefit start time t fit is indicated by the vertical grey dashed line. VI. DISCUSSION AND CONCLUSIONS
In this paper we have studied axion string networksin the radiation era by measuring the root-mean-squarevelocity of the strings v and the mean string separation ξ . The strings are modelled by a scalar field with tworeal components with a spontaneously broken O(2) sym-metry, simulated in periodic cubic lattices. Performinga phase space analysis in the variables x = ξ/t and v ,we find good evidence for the existence of a fixed point,which shows that the system reaches a scaling regime.2These prompted us to continue the analysis in the frame-work of the velocity-dependent one-scale (VOS) model.The VOS model assumes a statistical distribution ofstring positions and velocities which can be adequatelydescribed by two parameters mentioned above. By as-suming that the strings follow approximately Nambu-Goto trajectories, that they reconnect with a fixed highprobability when they cross, and that the loops so formedannihilate quickly, the VOS model reduces the networkevolution to a simple dynamical system, a pair of first or-der non-linear ordinary differential equations. The equa-tions have a fixed point ( x ∗ , v ∗ ), which describes a scalingnetwork, that is, one whose mean separation increaseslinearly with time, and with constant RMS velocity.We have fitted the results of a set of numerical simula-tions in the radiation era to a two-parameter VOS model,with initially random fields with several different initialcorrelations lengths l φ , using both the true physical fieldequations and the PRS approximation. In terms of thecore growth parameter s , the comoving string width be-haves as w = w a − s . In this paper we have used s = 1,which corresponds to the true physical case, and s = 0which corresponds to a string with constant comovingwidth.We find that the two-parameter VOS model gives agood qualitative and quantitative description of the net-work evolution, with parameters given in Table II.Qualitatively, the initial acceleration of the string net-work results in a RMS velocity which overshoots the fixedpoint v ∗ , as the Hubble length in our initial conditionsis larger than the string separation, meaning that thedynamical system is underdamped. The higher velocityresults in more rapid loop formation, and hence an in-crease in the mean string separation. This decreases theacceleration, and hence the RMS velocity. The net resultis a curved approach to the fixed point in the ( x, v ) plane,clearly visible in Fig. 7.In assessing the quantitative success of the two-parameter VOS model, we observe that the residuals tothe fits points in Fig. 8 are consistent with zero, and thatthe fixed points given in Table II for differential initialcorrelation lengths are remarkably similar. The fluctua-tions between the fixed point estimates are slightly largerthan the bootstrap fitting errors would predict, whichsuggests that the model could be tuned slightly, or thatthe fitting errors have been underestimated. A prelim-inary investigation shows that the more complex modelof Ref. [54] does not improve the fit. A more thoroughexploration of VOS models and a more accurate estimateof the fixed point could be obtained with a wider rangeof initial correlation lengths and initial times.Translating the values of Table II to the universe-framelength density parameter ζ , estimated in our previouspaper by linear fitting, we find ζ ∗ = 1 . ± .
09 ( s = 1) , (42) ζ ∗ = 1 . ± .
04 ( s = 0) . (43)These values are consistent with our previous determina- tion, with an improved accuracy arising from the joint fitwith the velocity data in the context of the VOS model.An important consequence of the description in termsof a dynamical system is that the approach to the fixedpoint is determined by a pair of complex exponents σ ± ,whose real part is − (1 − v ∗ ) (cid:39) − .
47. Hence, even whenclose to the fixed point, the approach can be rather slow.If the initial string separation ξ i is chosen far awayfrom its scaling value x ∗ t i , it may not get within 1 σ ofits scaling value (as determined by the VOS model) bythe end of the simulation, which has to be chosen as L/ L in systems like this one with degreesof freedom propagating at the speed of light.This is particularly noticeable for initial conditionswhich are very underdense, i.e. with x i (cid:29) x ∗ . Whenthe length density parameter ζ is plotted against the log-arithm of cosmic time log( tη ), one sees a slow drift uptowards the fixed point value. Other groups have alsonoticed this feature of underdense initial conditions [35–39, 41]. As we have explained elsewhere [40] this doesnot signal a breakdown of the standard scaling picture.Our analysis in the framework of the VOS model showsthat slow approaches to the fixed point from values of ζ less than its fixed point value are to be expected, andindeed, nearly all simulations to date have final values of ζ less than our estimated fixed point.If the slow upward drift in ζ were an asymptotic featureof axion string networks, one would expect to see finalvalues of ζ significantly above the fixed-point value inthe largest simulations. However, the maximum valueof ζ obtained in the most recent (and therefore largest)simulations are nearly all below value of ζ ∗ computed inthis work. These values are (all of them measured inthe universe frame): ζ (cid:39) . ζ (cid:39) . ζ (cid:39) . ζ (cid:39) . ζ (cid:39) ζ . They also suggest that abetter method will render their results comparable withthe ones in [41], and therefore, also with other groups.In summary, our data and fits already show that thestraightforward and physically motivated picture pro-vided by the simplest VOS model provides a good de-scription of the evolution of the string network consistentwith the standard scaling picture, and an asymptoticallyconstant dimensionless length density and RMS velocity.This gives confidence that our results can be extrapolatedover the many orders of magnitude required for predic-tions of the axion number density.Finally, we comment on a suggestion that simulationswith growing comoving string width shed light on theasymptotic behaviour of axion string networks [39]. If thestring width grows in proportion to the horizon ( w ∝ τ ),the field equations have a scale symmetry, and it canbe argued that a fixed point must exist. We use thisgrowth in width in the core growth phase of our s = 13simulations, where it can can be seen that this phasecharacterised by a constant RMS velocity, but not a con-sistent one. A smaller initial correlation length gives alarger RMS velocity in the core growth phase (see thetop panels in Figs. 9 and 5). This is understandable inthat the initial acceleration is proportional to the curva-ture. However, as the mean string separation increases,the RMS velocity of the strings does not decrease, aswould be expected from the decrease in the average ac-celeration. While it seems that x is evolving towards avalue x ∗ (cid:39) .
8, there is little sign of a definite value ofthe RMS velocity. It seems therefore that networks inthe core growth phase can be used to estimate the fixedpoint in x , but not the RMS velocity. ACKNOWLEDGMENTS
We are grateful to Daniel Cutting and Daniel G.Figueroa for comments on the draft manuscript. MH (ORCID ID 0000-0002-9307-437X) acknowledges sup-port from the Science and Technology Facilities Council(grant number ST/L000504/1) and the Academy of Fin-land (grant number 333609). JL (ORCID ID 0000-0002-1198-3191) and JU (ORCID ID 0000-0002-4221-2859) ac-knowledge support from Eusko Jaurlaritza (IT-979-16)and PGC2018-094626-B-C21 (MCIU/AEl/FEDER,UE).ALE (ORCID ID 0000-0002-1696-3579) is supported bythe National Science Foundation grant PHY-1820872.ALE is grateful to the Early Universe Cosmology groupof the University of the Basque Country for their gener-ous hospitality and useful discussions. This research wassupported by the Munich Institute for Astro- and Par-ticle Physics (MIAPP) which is funded by the DeutscheForschungsgemeinschaft (DFG, German Research Foun-dation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. This work has been possible thanks tothe computational resources on the STFC DiRAC HPCfacility obtained under the dp116 project. Our simula-tions also made use of facilities at the i2Basque academicnetwork and CSC Finland. [1] S. Weinberg, Phys. Rev. Lett. , 223 (1978).[2] F. Wilczek, Phys. Rev. Lett. , 279 (1978).[3] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. , 1440(1977).[4] R. D. Peccei and H. R. Quinn, Phys. Rev. D16 , 1791(1977).[5] J. E. Kim, Phys. Rev. Lett. , 103 (1979).[6] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov,Nucl. Phys. B166 , 493 (1980).[7] A. R. Zhitnitsky, Sov. J. Nucl. Phys. , 260 (1980), [Yad.Fiz.31,497(1980)].[8] M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. , 199 (1981).[9] J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B120 , 127 (1983).[10] L. F. Abbott and P. Sikivie, Phys. Lett.
B120 , 133(1983).[11] M. Dine and W. Fischler, Phys. Lett.
B120 , 137 (1983).[12] A. Vilenkin and A. Everett, Phys. Rev. Lett. , 1867(1982).[13] R. L. Davis, Phys. Lett. B180 , 225 (1986).[14] M. Hindmarsh and T. Kibble, Rept.Prog.Phys. , 477(1995), arXiv:hep-ph/9411342 [hep-ph].[15] A. Vilenkin and E. P. S. Shellard, Cosmic Strings andOther Topological Defects (Cambridge University Press,2000).[16] P. Sikivie, Phys. Rev. Lett. , 1156 (1982).[17] H. Georgi and M. B. Wise, Phys. Lett. , 123 (1982).[18] C. Hogan and M. Rees, Phys. Lett. B , 228 (1988).[19] E. W. Kolb and I. I. Tkachev, Phys. Rev. Lett. , 3051(1993), arXiv:hep-ph/9303313.[20] E. W. Kolb and I. I. Tkachev, Phys. Rev. D , 769(1994), arXiv:astro-ph/9403011.[21] E. W. Kolb and I. I. Tkachev, Astrophys. J. Lett. ,L25 (1996), arXiv:astro-ph/9510043.[22] M. Fairbairn, D. J. E. Marsh, and J. Quevillon, Phys. Rev. Lett. , 021101 (2017), arXiv:1701.04787 [astro-ph.CO].[23] M. Fairbairn, D. J. E. Marsh, J. Quevillon,and S. Rozier, Phys. Rev. D , 083502 (2018),arXiv:1707.03310 [astro-ph.CO].[24] M. Yamaguchi, M. Kawasaki, and J. Yokoyama, Phys.Rev. Lett. , 4578 (1999), arXiv:hep-ph/9811311 [hep-ph].[25] M. Yamaguchi, M. Kawasaki, and J. Yokoyama, in Proceedings, 7th International Symposium on Particles,Strings and Cosmology (PASCOS 99): Lake Tahoe, Cal-ifornia, December 10-16, 1999 (1999) pp. 291–294.[26] M. Yamaguchi, M. Kawasaki, and J. Yokoyama, in
Pro-ceedings, 3rd International Workshop on The identifica-tion of dark matter (IDM 2000): York, UK, September18-22, 2000 (2000) pp. 297–304.[27] M. Yamaguchi and J. Yokoyama, Phys. Rev.
D67 ,103514 (2003), arXiv:hep-ph/0210343 [hep-ph].[28] T. Hiramatsu, M. Kawasaki, T. Sekiguchi, M. Yam-aguchi, and J. Yokoyama, Phys. Rev.
D83 , 123531(2011), arXiv:1012.5502 [hep-ph].[29] T. Hiramatsu, M. Kawasaki, K. Saikawa, andT. Sekiguchi, Phys. Rev.
D85 , 105020 (2012), [Erratum:Phys. Rev.D86,089902(2012)], arXiv:1202.5851 [hep-ph].[30] M. Kawasaki, K. Saikawa, and T. Sekiguchi, Phys. Rev.
D91 , 065014 (2015), arXiv:1412.0789 [hep-ph].[31] L. Fleury and G. D. Moore, JCAP , 004 (2016),arXiv:1509.00026 [hep-ph].[32] A. Lopez-Eiguren, J. Lizarraga, M. Hindmarsh, andJ. Urrestilla, JCAP , 026 (2017), arXiv:1705.04154[astro-ph.CO].[33] V. B. Klaer and G. D. Moore, JCAP , 043 (2017),arXiv:1707.05566 [hep-ph].[34] V. B. Klaer and G. D. Moore, JCAP , 049 (2017),arXiv:1708.07521 [hep-ph].[35] M. Gorghetto, E. Hardy, and G. Villadoro, JHEP ,
151 (2018), arXiv:1806.04677 [hep-ph].[36] M. Kawasaki, T. Sekiguchi, M. Yamaguchi,and J. Yokoyama, PTEP , 091E01 (2018),arXiv:1806.05566 [hep-ph].[37] A. Vaquero, J. Redondo, and J. Stadler, JCAP ,012 (2019), arXiv:1809.09241 [astro-ph.CO].[38] M. Buschmann, J. W. Foster, and B. R. Safdi, Phys.Rev. Lett. , 161103 (2020), arXiv:1906.00967 [astro-ph.CO].[39] V. B. Klaer and G. D. Moore, JCAP , 021 (2020),arXiv:1912.08058 [hep-ph].[40] M. Hindmarsh, J. Lizarraga, A. Lopez-Eiguren, andJ. Urrestilla, Phys. Rev. Lett. , 021301 (2020),arXiv:1908.03522 [astro-ph.CO].[41] M. Gorghetto, E. Hardy, and G. Villadoro, (2020),arXiv:2007.04990 [hep-ph].[42] M. Gorghetto, E. Hardy, and H. Nicolaescu, (2021),arXiv:2101.11007 [hep-ph].[43] M. Dine, N. Fernandez, A. Ghalsasi, and H. H. Patel,(2020), arXiv:2012.13065 [hep-ph].[44] T. Braine et al. (ADMX), Phys. Rev. Lett. , 101303(2020), arXiv:1910.08638 [hep-ex].[45] T. W. B. Kibble, Phase transitions in the very early Univesre. proceedings, International Workshop, Biele-feld, F.R. Germany, June 4-8, 1984 , Nucl. Phys.
B252 ,227 (1985), [Erratum: Nucl. Phys.B261,750(1985)].[46] C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev.
D54 , 2535 (1996), arXiv:hep-ph/9602271 [hep-ph].[47] C. Martins and E. Shellard, Phys.Rev.
D65 , 043514(2002), arXiv:hep-ph/0003298 [hep-ph].[48] C. J. A. P. Martins, Phys. Lett.
B788 , 147 (2019),arXiv:1811.12678 [astro-ph.CO].[49] T. Kibble, J.Phys. A9 , 1387 (1976).[50] W. Zurek, Phys. Rept. , 177 (1996), arXiv:cond-mat/9607135.[51] W. H. Press, B. S. Ryden, and D. N. Spergel, Astrophys.J. , 590 (1989).[52] T. Vachaspati and A. Vilenkin, Phys. Rev. D30 , 2036(1984).[53] M. Hindmarsh, J. Lizarraga, J. Urrestilla, D. Dave-rio, and M. Kunz, Phys. Rev.
D96 , 023525 (2017),arXiv:1703.06696 [astro-ph.CO].[54] J. R. C. C. C. Correia and C. J. A. P. Martins, Phys. Rev.D100