Long-term dynamics of cosmological axion strings

We present results of new field-theoretic simulation of cosmological axion strings, which are eight times longer than previous ones. We have upgraded our simulation of physical strings in Hiramatsu et al. (2011) in terms of the number of grids as well as the suite of analysis methods. These improvements enable us to monitor a variety of quantities characterizing the dynamics of the physical string network for the longest term ever. Our extended simulations have revealed that global strings do not evolve according to the scaling solution but its scaling parameter, or the number of long strings per horizon, increases logarithmically in time. In addition, we have also found that the scaling parameter shows nontrivial dependence on the breaking scale of the Peccei-Quinn symmetry.


Introduction
The axion is one of the best motivated particles beyond the Standard Model [1,2]. Its existence is inevitably derived from the Peccei-Quinn (PQ) mechanism [3], which dynamically solves the strong CP problem in quantum chromodynamics (QCD). Concurrently, it can account for the cold dark matter (CDM) in the observed Universe and offer rich phenomenology in the cosmological context (see e.g. [4,5] for recent review). There are a few major production mechanisms of the axion CDM in the early Universe, including the misalignment and topological defects. A number of experiments are going on aiming at the detection of the CDM axion [6] (for review we refer to [7,8,9,10]).
Formation of topological defects associated to the broken global PQ U(1) symmetry takes place if the symmetry breaking occurs after the observable Universe exits horizon during inflation. At first, axion strings form following the PQ phase transition. As the cosmological network of the strings evolves, its energy is released in the form of axion radiation. Later on, the QCD phase transition takes place and domain wall (DW) stretches out between the axion strings. If the number of DWs attached to each string is unity, the network of string-DW is unstable, so that it disappears soon after the QCD phase transition leaving additional contribution to the axion CDM [11]. Otherwise, DWs are stable and dominate the Universe, which spoils successful subsequent cosmology (see also [12]).
In this paper, we will focus on the network of axion strings in between the two phase transitions. This subject has already been studied by many different people for decades [13,14,15,16,17,18,19,20,21,22,23,24,25,26]. However, the current picture of the string dynamics could not be developed straightforwardly. Notably, there had been long-standing controversies concerning, for instance, the number of strings per horizon and energy spectrum of radiated axions [13,14,15,16,17,18,19,20,21,24]. Those controversies were settled principally with the aid of field theoretic simulation of the strings based on the first principles [21,22,23,25,26,27,28]. Nonetheless, as we shall reveal its novel dynamics in this Letter, axion strings seem still concealing its nature from us. We need to pursue better understanding of dynamics of global strings in order for accurate prediction of the axion CDM abundance.
The primary purpose of this paper is to update the result of our previous simulation [27]. In that paper, we performed field theoretic simulation of cosmological axion strings and estimated the relic abundance of axion radiated from the strings. Our simulation is improved mainly in two respects. Exploiting massively parallel computation on computer clusters, the number of grids is increased from 512 3 to 4096 3 , which enhances the simulation time by a factor of eight. This allows us to examine long-term behavior of the string network. We also improved our analysis method. In [27], we introduced a novel string identification method and statistical reconstruction of the axion energy spectrum. In the present analysis, we incorporated estimation of velocity of strings and loop identification of strings in the suite of analysis. This allows us to monitor a variety of quantities characterizing the string network and examine its dynamics in detail.
This paper is organized as follows. In the next section, we will first describe the Lagrangian of the PQ scalar we simulate as well as the essence of our numerical simulation and analysis. Then in Section 3 we will present relevant results from our simulation. In particular, we pay particular attention to the long-term behavior of the string network, that is found in simulation of physical strings for the first time. The final section will be devoted to discussion.
2 Setup and analysis methods

Model
We adopt the following Lagrangian density for the PQ complex scalar Φ as in [26]: where the effective potential V eff [Φ; T ] is dependent on the temperature T : with v being the vacuum expectation value of Φ at T = 0. The QCD axion ϕ is identified as the phase of Φ. At low energy, Φ can be written as Φ = fa 2v is the decay constant of the axion. The mass of PQ field is also temperature dependent and given by m 2 (T ) = λ( T 2 3 − 2v 2 ). Given the potential, the PQ phase transition occurs at the critical temperature T c = √ 6v. In what follows, quantities measured at T = T c are subscripted with c.
We assume a flat Universe with its line element given by We normalize the scale factor at the PQ phase transition (i.e. a c = 1). We assume the radiation domination and the Hubble expansion rate H is given by where g * and M * are the relativistic degrees of freedom and the reduced Planck mass, respectively. Note that t ∝ a 2 and the conformal time τ is given by t/2a in the radiation domination. In what follows, we denote a derivative with respect to t and τ with a dot and a prime (i.e. ∂f ∂t =ḟ and ∂f ∂τ = f ), respectively. Following [21], we adopt a parameter Granting that the physical width of the axion string d is given by d ( √ λv) −1 , ζ parameterizes the ratio of the Hubble scale to the string width at the time of the PQ phase transition. In the following analysis, we fix g * = 1000 #1 and λ = 1 but vary v among v/M * = 0.005, 0.002 and 0.001, which respectively give ζ = 9.5, 23.9 and 47.8.

Simulation
We start the simulation immediately before the PQ phase transition, and initial condition of the real and imaginary parts of the PQ field Φ = (φ 1 + iφ 2 )/ √ 2 is generated assuming Gaussian random field with thermal distribution: where a, b = 1, 2 and E k = k 2 a 2 + m 2 . We integrate the classical equation of motion of Φ using the leapfrog method. For detailed description of how the equation is discretized on the lattice, we refer readers to the appendices of [26]. The number of grids is N g = 4096 3 , which is 512 times as large as in the previous simulations [27]. In order to suppress boundary effects, the size of the comoving simulation box L is set to twice as large as the horizon scale at the final time (i.e.
where the subscript f indicates the final time. The final time t f is determined by t f /d N 1/3 g /4 10 3 so that the lattice spacing does not exceed the string width at t f .

Analysis
We perform a suite of analysis in order to extract a variety of dynamical quantities associated to the string network. One of the essential techniques in the analysis is identification of axion strings from the discretized data of the PQ field on grids. We adopt the same string identification method as in [27], which tells us not only the existence of a string inside a given cell, but also its location in the cell if it does exist.
Furthermore, this time we introduce a novel method to identify string loops from the string points obtained from the string identification method. This is realized by grouping those string points by the Friends-of-Friends algorithm. We found that setting the physical linking length to (d 2 t) 1/3 works in practice. Once a distinct loop is identified, #1 The choice of g * is not relevant in our analysis, where we focus on slowly time-varying quantities such as the string parameter ξ and the ratio of radiation momentum to the Hubble rate . We chose the value of g * so that the direct comparison with the previous studies [21,27] can be straightforward.
we can compute its circumference length. By accumulating them we estimate the string parameter ξ, which is defined by where µ is the tension of the axion string. The physical meaning of ξ is the average number of strings per horizon volume. We also implant an estimation method of string speed following [26]. On a point in the vicinity of a string, the string velocity can be estimated by This is not exactly the same estimation method as in [26], but coincides when Φ vanishes on the point where we are trying to estimate the string velocity. Sometimes we found the magnitude of v exceeds unity. This failure of estimation is caused where a string has large curvature and/or strings are colliding. However, the fraction of the failure is at worst five percent and does not significantly bias the estimate of the rms of v that will be presented in the next section.
On the other hand, the energy spectrum of the axion radiation is computed according to the pseudo-power spectrum estimator [27]. The method consists mainly of two processes. First, we configure a window function which masks cells in proximity to string cores and compute the energy spectrum of axion convolved with the window function. Then we reconstruct the energy spectrum statistically by deconvolving the mode-mixing caused by the window function. This procedure can remove the contamination from the string cores, which dominates over the contribution genuinely from the axion radiation at high wave numbers.
More detailed description of our simulation and analysis will be presented in a separate paper [34].  [27] is also shown (red filled square), whose dynamic range is only as large as t f /d 40. #2 As can be read from the figure, ξ grows in time and the time scale of the growth can be asymptotically characterized in ln(t/d). #3 Although the magnitude of ξ does not agree among v/M * = 0.005, 0.002 #2 We note that our previous choice of t f /d = 40 is conservative compared to that expected from N g (magenta arrow). This makes the actual enhancement of the dynamic range more than eight and actually close to twenty. #3 We note that when we fit ξ with a power function of log(t/d), that is, [log(t/d)] α at log(t/d) ≥ 300, the best-fit power indices α are 0.64 ± 0.02, 0.98 ± 0.03 and 0.45 ± 0.02 for v/M * = 0.005, 0.002, and 0.001, respectively. . For reference, we plotted the result of simulation corresponding to our previous study [27] (red filled square). Error bar shows the standard deviation among realizations. Magenta arrow shows the dynamic range that should have been available in the previous simulation [27]. and 0.001, the asymptotic logarithmic growth is common. Such logarithmic growth has not been reported in the previous studies of physical strings including ours [27] due to the limitation of the dynamic range. Those simulations with such small dynamic range could trace the evolution only up to t/d O(10) where such long-term dynamics has been buried in statistical uncertainties. The larger number of grids in our simulation, which enables longer duration of simulation, is essential in finding this novel property.

Results
On the other hand, ξ shown in Figure 1 incorporates both infinite strings and loops. Thus, one may wonder the logarithmic growth occurs only among loops or infinite strings. To answer this question, in Figure 2 we plot the string parameter only contributed from loops with length less than t (left) and πt (right). From the figure, one can see that the contribution of the loops in the total length is around 10 percent. Thus, the dominant contribution in ξ comes from the infinite strings. In addition, the figure shows that the string parameter of loops does not apparently grow in time. Therefore, we conclude that the logarithmic growth seen in Figure 1 originates from the infinite strings.
Velocity of strings also exhibits some degree of discrepancy from the previous studies. The root mean square of string velocity v 2 is plotted in Figure 3. We found 0.6 v 2 0.7 around the initial stage, which is consistent with [26]. However, as time advances, it decreases gradually and reads around 0.5 at the later time of the simulation. This is substantially smaller compared to the velocity estimated from simulation of local-string based on the Nambu-Goto action [29,30] but consistent with field-theoretic simulation of local strings [31]. At this moment, we are unable to tell if the velocity has settled down already within the simulation time or continues to decrease subsequently.
In Figure 4 we plot the evolution of the differential energy spectrum of the axion radiation, d(a 4 ρ rad ) d ln k , for v/M * = 0.005. This demonstrates two aspects of the axion radiation  Figure 2: Evolution of string parameter contributed only from the loops with physical circumference length less than t (left) and πt (right). Color and symbol are the same as in Figure 1.  from the axion strings. One is that energy spectrum has a peak at low momenta and decays quickly towards higher ones. The other is that the spectral peaks moves towards lower momenta as time advances.
To quantify the evolution of the spectral peaks, let us define where k −1 is the mean inverse momentum of radiated axion at t defined as with a −4 dkP rad being the energy density of the axion radiation. Roughly speaking, represents the mean wave number of radiated axion in units of the Hubble expansion rate. Figure 5 shows the evolution of . We confirm lies between two and four, which is largely consistent with [21,27]. This supports the claim that radiated axions have typical momenta comparable to the Hubble expansion rate [13,14,16,19,20]. Due to the limitation of the number of bins, which is taken to 75 in our analysis, we cannot resolve the spectral peaks any longer when the mean wave number becomes comparable to the bin size. This is why we can trace only up to t/d ∼ 100.
On the other hand, the wavelength of radiated axions are expected to be related to the correlation length of the strings, which is proportional to 1/ √ ξ. This implies should be proportional to √ ξ. Nonetheless, at this moment we cannot confirm such evolution of due to the lack of late-time measurements of .

Discussion
Among the results we presented in the previous section, the growth of string parameter would have the most profound significance. While field theoretic simulation of physical strings cannot trace t/d (N 1/3 g ), log(t/d) can be as large as 70 at the time of the QCD phase transition given 10 9 GeV v 10 11 GeV, where the lower and upper bounds come from astrophysical consideration and cosmological axion abundance, respectively [32,33,10]. This means even the merely logarithmic growth can enhance ξ by a factor of seven or larger #4 .
The number of axions radiated from strings is proportional to the number of axion strings and inversely proportional to the typical momentum of radiated axions. The former is proportional to ξ. On the other hand, the latter is expected to be proportional to √ ξ, though we couldn't confirm this in our present simulation. Assuming that the axion momentum grows in proportion to √ ξ, we expect the axion abundance is proportional √ ξ. Our results indicates the abundance of axion can be enhanced by a factor of two or three. This should give a large impact not only on the constraints on the decay constant f a from overabundance of the axion CDM but also detection experiments, which target the mass of axion or f a ∝ v with which the axion accounts for the observed CDM density.
Logarithmic growth in the string parameter is reported in [35,36,37]. However, the relation between our observation and theirs is not at all clear. This is because there is a crucial difference between the two simulations. Our strings are physical (i.e. constant string width), while the strings in [35,36,37] are so-called fat strings (i.e. string width growing in time). As discussed in [17], the radiation backreaction, which should be crucial in the dynamics of global strings [13,38], is suppressed by the logarithm of t/d (see also [30]). Since d increases in fat strings, the extent of suppression of the radiation backreaction may deviate from that in physical strings. Therefore, results of fat strings cannot be easily #4 Fitting ξ with a power-law dependence on log(t/d) at log(t/d) ≥ 300, the enhancement factor can be 4.6 ± 0.68, 8.8 ± 1.6 and 3.8 ± 0.4 for v/M * = 0.005, 0.002, and 0.001, respectively. identified with physical ones.
What is more puzzling is the nontrivial dependence of ξ on the PQ breaking scale v. As seen in Figure 1, we found the magnitude of ξ depends on v. Although previously reported values of ξ have also been varying among literature [21,30,26,27], the discrepancy has not been paid much attention due to statistical errors. Thanks to larger box size available at given t, the discrepancy in ξ is manifested in our analysis with enough significance.
Such dependence apparently conflicts with the prediction of one-scale model [39,40,41]. Moreover, the dependence of ξ on v is not monotonic. As can be seen from Figure 1, the magnitude of ξ fluctuates as v/M * is decreased from 0.005 to 0.001. This may indicate existence of competing effects, which explicitly depend on the PQ breaking scale.
At the moment, the physical origin of the logarithmic growth of the string parameter ξ and its dependence on v is still cloaked. As noted above, the efficiency of the energy loss of strings through radiation is suppressed by log(t/d). We speculate this can at least partially account for the logarithmic growth of ξ. Nonetheless, more rigorous investigation is required before we conclude this is the dominant cause. In view of the impact, scrutiny of the origin cannot be stressed too much. We pursue more detailed analysis with extensive parameter range and physical interpretation in the future publication [34].
We have also computed the spectrum of the axion radiated from the axion strings. The result is largely in agreement with the previous studies [21,27]. We have confirmed the spectral peak stays around a few time larger than the Hubble expansion rate.
Note added: After we have finished writing the substantial part of the present manuscript we became aware of a paper by Gorghetto, Hardy, and Villadoro, who performed very similar numerical analysis [42].