capbtabboxtable[][]

#  Axion Dark Matter
from Cosmic String Network

Heejoo Kim${}^{\,\dagger}$, Junghyeon Park${}^{\,\dagger}$, and Minho
Son$\,{}^{\dagger}$

(†Department of Physics, Korea Advanced Institute of Science and Technology,
291 Daehak-ro, Yuseong-gu, Daejeon 34141, Republic of Korea
)

###### Abstract

We perform the lattice simulation to estimate the axion dark matter abundance
radiated from the global cosmic strings in the post-inflationary scenario. The
independent numerical confirmation on the recently observed logarithmic growth
in both the number of strings per Hubble patch and the spectral index of the
power law scaling for the axion spectrum is reported. These logarithmic
scalings are checked against two different prescriptions for generating
initial random field configurations, namely fat-string type and thermal phase
transition. We discuss a possible strong correlation between the axion
spectrum and the string evolutions with different initial conditions to
support the insensitivity of scaling behaviors against different initial data
and we provide a qualitative understanding of it. The impact of various
combinations of the power law of the axion spectrum, nonlinearities around the
QCD scale, and average inter-string distances on the axion abundance is
discussed. Additionally, we introduce a new novel string identification
method, based on the tetrahedralization of the space, which guarantees the
connectedness of the strings and provides a convenient way of assigning the
core location. Finally we derive the lower bound on the axion mass.

###### Contents

  1. 1 Introduction
  2. 2 Formation of cosmic strings
  3. 3 Simulation Setup
    1. 3.1 Discretization
    2. 3.2 Initial conditions and Relaxation
  4. 4 String identification on tetrahedron
  5. 5 Scaling regime
    1. 5.1 Scaling of inter-string distances
  6. 6 Cosmological evolution of string network
    1. 6.1 Energy density
    2. 6.2 String tension and boost factor
  7. 7 Axion spectrum
    1. 7.1 Power law scaling
  8. 8 Axion abundance
    1. 8.1 Bound on axion mass
  9. 9 Correlation between strings and axion spectrum
  10. 10 Conclusion
  11. A Equation of motion
  12. B Evolutions on the lattice
    1. B.1 Common setup
    2. B.2 Fat string pre-evolution: $R(t)\propto\sqrt{t}$
    3. B.3 Fat string pre-evolution: $R(t)\propto t$
    4. B.4 From pre-evolution to physical one
    5. B.5 Thermal pre-evolution
    6. B.6 Random initial conditions
    7. B.7 Tetrahedron-based string identification
    8. B.8 Masking
  13. C Dependency on initial conditions and relaxations
    1. C.1 Scaling regime and fitting
    2. C.2 Cosmological evolution of string network
      1. C.2.1 Energy budget
      2. C.2.2 Instantaneous emission
  14. D More on axion spectrum
    1. D.1 Fitting instantaneous emission function
    2. D.2 Analytic understanding of positive correlation
  15. E More on axion abundance
    1. E.1 Estimation for $q>1$
    2. E.2 Estimation for $q=1$

##  1 Introduction

Axion as dark matter has been increasingly attractive, although it was
originally introduced in some extensions of the Standard Model to resolve a
mysterious fine-tuning in quantum chromodynamics (QCD) [1, 2, 3, 4, 5, 6].
While the QCD axion dark matter scenario is a certainly interesting
possibility ([7] for a review), its dominant production mechanism is still
inconclusive. The topological production such as axions radiated from the
global cosmic string network [8] formed during the phase transition of the
Peccei-Quinn (PQ) symmetry breaking [1] after inflation in the early universe
can be predicted in terms of the symmetry breaking scale. However, it entirely
relies on the numerical simulation as the effective field theory (EFT)
description does not work due to the evolution occurring between two vastly
different scales for the PQ symmetry breaking and QCD crossover and the
nonlinearity of the equation of motion. While the precise estimation of the
topological production with respect to the misalignment mechanism has been
subject to a large uncertainty, any progress for the improvement should be
crucial in understanding its origin and role of the QCD axion as the dark
matter. Especially, the outcome will have a large impact on the direct search
for the QCD axion as the ballpark for the natural axion mass will change.

Using the advanced computing resources, some state-of-art simulations have
been performed by several groups on the largest lattice grids of the order
$\mathcal{O}(2^{12\sim 13})$ per a dimension. An intriguing discovery of those
recent simulations was the strong evidence for the scaling regime, or an
attractor solution, exhibiting the logarithmic growth in the number of long
strings per Hubble patch, denoted by $\xi$ [9, 10, 11, 12, 13, 14, 15, 16]. In
this work, we independently confirm the logarithmic growth within the
simulation time coverage. Since the cosmic strings will eventually decay into
axions in late times, the evolution of $\xi$ is an important property that has
to be estimated as precise as possible to determine the axion abundance. If
the logarithmic growth in $\xi$ holds until the late times for the QCD
crossover, its prediction will be very different from the constant scaling
scenario of the order one value [17]. The similar constant scaling behaviors
observed in early simulations on the smaller lattice grids [18, 19, 20, 21,
22, 23] may not be helpful to support either scenario due to the limited
dynamic time range (see related discussions in [24, 25, 26, 27, 28]).

The differential axion energy spectrum is another crucial property for the
accurate estimation of the axion abundance. In the absence of any non-trivial
new scale between the scale of the PQ symmetry breaking and the QCD phase
transition, it is expected to follow the power law fall-off behavior. The
characteristic scaling feature of the axions radiated from the topological
cosmic strings can be captured by the instantaneous emission. However, there
has been a disagreement in the determination of the power law scaling $\sim
k^{-q}$ in the instantaneous axion emission spectrum. While the static lattice
simulation performed in [14] finds a hint on the logarithmic growth of the
spectral index $q$ in time, an independent simulation in [15] by the adaptive
mesh refinement (AMR) 111The AMR method can parametrically extend the
simulation time by locally increasing the lattice resolution with recursively
created finer grids around strings on-the-fly. The AMReX [29, 30], a software
for block-structured AMR, was adopted in [15]. GRChombo [31] is another AMR
based open-source code originally developed for numerical relativity
simulation and it was recently applied to cosmic strings [32]. [33] finds a
consistency with the no-log hypothesis. The former predicts the dominant
topological production over the misalignment mechanism (which favors a heavier
axion mass than the natural ballpark for the QCD axions from misalignments)
whereas the latter predicts a comparable size. Two simulation results in [14,
15] were obtained by different setups, differing by how the spacetime were
discretized, how the initial conditions were generated, and how they were
subsequently evolved over time. This makes the transparent comparison quite
challenging. In this work, we perform independent static lattice simulations
on the grids of $N^{3}=4096^{3}$ by two commonly used approaches in
literature: one using the fat-string type pre-evolution to prepare for initial
conditions for the physical string evolution [9, 14] and the other scheme
where the field evolution occurs in the finite temperature, what we will call
the thermal pre-evolution (for instance, our setup for this approach is
similar to those in [10, 15]). While the fat-string pre-evolution may be
criticized that it is not the same theory as the one for the physical string
evolution, the direct comparison between simulations adopting different
schemes should greatly help establishing the validity of the final
observations or the insensitivity to the aforementioned relaxation schemes.
The result in [10] has not been reported in similar presentations to [9, 14],
making it difficult to be used for the meaningful comparison. The AMR
technique, adopted in [15], allows to cover much longer dynamic time range
than what one has ever imagined to be possible. However, whether it eventually
leads us to the better determination in the region beyond the time coverage of
the static lattice simulation does not seem to be straightforward. The
extending our comprehensive comparison in the context of AMR will be presented
in our companion paper.

The parametric behavior of the axion abundance originated from the scaling
regime can be significantly affected by the nonlinearities around the QCD
scale [14]. In particular, it varies depending on the value of $q$ ($q>1$ vs
$q=1$) and axion field values (large vs small) around the QCD scale
(equivalent to the relevance of the nonlinearities) and the lower momentum
cutoff (dictated by average inter-string distances). As different combinations
match to different observations in literature, it should be highly useful to
work out all relevant cases. Following the similar strategy to [14], we
discuss the impact of the modified momentum cutoff (suggested by the scaling
of the inter-string distances) and axion field values on the axion abundance,
and we extend the discussion to the case with $q=1$ to estimate the relevance
of the nonlinearities around the QCD scale in that situation. This will make
the comparisons between different scenarios transparent.

While there are already various string identification algorithms in
literature, in this work, we propose a new novel string identification
algorithm, based on the tetrahedralization of the space, which guarantees the
connectedness of strings. As the name stands, the algorithm works on the
tetrahedron instead of the cube. It has a few practical advantages that will
be discussed below. Compared to existing techniques [34, 35, 21, 22], its
performance time is comparable with them.

The most severe limiting factor of the simulation is CPU and memory resources.
We have greatly optimized our own independent codes (made from the scratch)
running on OpenMPI to speed up each simulation. Nevertheless, obtaining a
large size of ensemble, performing at the same time all kinds of measurements
such as $\xi$, energy budget, power spectrum, instantaneous emission, string
tension, string velocity and checking out the dependencies of various
prescriptions for each measurement, from the lattice simulation on the grids
of $4096^{3}$ is challenging. Therefore, we extensively used the lattice
simulation on smaller scales such as $1024^{3}$ and $2048^{3}$ as well for
numerous sanity checks and to get guideline results for an educated selection
of the benchmark setups for the simulation on the bigger sized lattice of
$N^{3}=4096^{3}$. Our final result on the axion spectrum, which supports the
soft axions of the order of the Hubble scale around the QCD crossover,
strongly suggests the simulation on the factor of two larger lattice space.

The paper is organized as follows. In Section 2, we introduce the model for
the formation of the global axion strings. In Section 3, brief descriptions on
the discretization of the equation of motion and random field configurations
are given. In Section 4, a new novel string identification algorithm based on
the tetrahedralization is introduced. In Section 5, we present our independent
evidence for the logarithmically growing scaling solution in the number of
strings in the Hubble patch. The scaling of the average inter-string distance
is discussed. In Section 6, the detailed discussion on the energy budget and
string properties such as the string tension and velocity are given. In
Section 7, we provide our fit results on the spectral index of the power law
in the axion spectrum. In Section 8, we discuss the prospect for the axion
abundance for given spectral index from our simulation. In Section 9, a
possible correlation between axion spectrum and string evolutions with
different initial field configurations is discussed. In Appendices, all
simulation details are given. Extra detailed discussions and materials are
given for various subjects. Especially, the parametric behavior of the axion
abundance for several distinctive situations are worked out.

Note added: while this work was being finalized, another work addressing
similar aspects of axions from cosmic strings has appeared [36].

##  2 Formation of cosmic strings

Topological cosmic strings are generated upon the phase transition of the
$U(1)_{PQ}$ symmetry by the complex scalar field $\phi$ in the early Universe
[37, 38, 39]. The PQ symmetry breaking is conveniently parametrized by the
Lagrangian,

$\mathcal{L}=\partial_{\mu}\phi^{*}\partial^{\mu}\phi-\frac{m_{r}^{2}}{2f_{a}^{2}}\left(|\phi|^{2}-\frac{f_{a}^{2}}{2}\right)^{2}\leavevmode\nobreak\
,$ (1)

where the quartic coupling corresponds to
$m_{r}^{2}/(2f_{a}^{2})\equiv\lambda$, often denoted by $\lambda$ in
literature, and $f_{a}$ is the symmetry breaking scale. The dynamical
evolution of the cosmic strings in the expanding universe will be explored,
assuming the radiation domination whose metric is given by

$ds^{2}=dt^{2}-R^{2}(t)d\vec{x}^{2}\leavevmode\nobreak\ ,$ (2)

where $\vec{x}$ denotes the comoving coordinates. The equation of motion from
the Lagrangian in Eq. (1) in the expanding universe is given by

$\ddot{\phi}+3H\dot{\phi}-\frac{1}{R^{2}}\nabla^{2}\phi+\frac{m_{r}^{2}}{f_{a}^{2}}\phi\left(|\phi|^{2}-\frac{f_{a}^{2}}{2}\right)=0\leavevmode\nobreak\
,$ (3)

where dot is the differentiation with respect to the cosmic time $t$ and the
gradient $\nabla$ is evaluated in the comoving coordinate. The equation of
motion in Eq. (3) admits the solitonic solution whose configuration
corresponds to the topological cosmic strings. Cosmic strings formed at the PQ
scale will evolve for a long dynamic time range all the way up to the QCD
crossover where the abundance of the axion dark matter radiated from cosmic
strings are estimated. The dynamics of the cosmic strings are traced through
the evolution of the fields. As they are highly nonlinear in between largely
separated two scales, the EFT approach is not available and it has to proceed
through the numerical simulation on the lattice. The numerical simulation of
the equation of motion in Eq. (3) on the lattice is limited since it involves
three different scaling behaviors. While the lattice spacing scales as
$R\propto\sqrt{t}$ in the radiation-dominated universe, the string core width
$m_{r}^{-1}$ stays the same and the Hubble length $H^{-1}$ scales as $\sim t$
during the evolution. At some point, the simulation runs out of the resolution
within the string core or the correlation length $\sim H^{-1}$ becomes
comparable with the size of the simulation box where the boundary conditions
become relevant.

The field evolution on the lattice in terms of $\phi$ is subject to the
convergence condition, or Courant–Friedrichs–Lewy (CFL) condition, $\Delta
t<R\Delta x$. This condition requires arbitrarily small time step size in the
early stage of the evolution where the scale factor is small. It can be
relaxed by re-expressing the equation of motion in terms of the rescaled field
$\psi(\tau,\,\vec{x})=R\phi(t,\,\vec{x})$ in the comoving spacetime
coordinate. The CFL condition in this situation is relaxed to
$\Delta\tau<\Delta x$. All the simulation details are given in Appendices A
and B.

##  3 Simulation Setup

###  3.1 Discretization

The equation of motion in Eq. (3) is discretized in the conformal frame
according to the leap-frog method whose truncation error is the order of
$\mathcal{O}((\Delta x)^{3})$. The lattice spacing $\Delta x$ is determined by
$m_{r}$ for the given lattice size $N$. It is set to reach the maximum dynamic
time range allowed by the lattice size, namely $\Delta
x=\frac{m_{r}^{-1}}{n_{c}}\left(\frac{2n_{c}n_{H}}{N}\right)^{1/2}$ where
$n_{c}$ is the number of lattice grids inside the core and $n_{H}$ is the
number of Hubble patches in the simulation box at the final time. $n_{c}=1$
and $n_{H}=4^{1/3}$ will be our default choice in the simulation unless
otherwise specified. $n_{H}=4^{1/3}$ ensures at least four Hubble volumes in
the simulation box at the final time. See Appendix B.1 for the detail. The
time interval $\Delta\tau$ is chosen to be $\Delta\tau/\Delta x=1/3$ so that
it satisfies the CFL condition. Our simulation data set is summarized in Table
3.

###  3.2 Initial conditions and Relaxation

The string configurations produced out of initial conditions in the very
beginning of the simulation are highly noisy. Given the limited dynamic time
range of the physical string simulation, a pre-evolution for the purpose of
relaxing the string network to a relatively clean level seems to be a
necessary step. The relaxation step may be justified if the string network
eventually approaches toward the scaling regime [40, 41, 42] and becomes
insensitive to the initial configurations.

The fat-string evolution [43, 44] was adapted in recent works [9, 14] for the
purpose of the relaxation (see Appendices B.2 and B.3 for details). The
outcome of the fat-string pre-evolution feeds into the physical string
evolution as the initial condition. It is called fat-string since the string
core width is forced to scale as $R$ during the pre-evolution so that
$\frac{m_{r}^{-1}}{R}$ stays constant. It is implemented in the simulation by
the quartic coupling scaling as $R^{-2}$ for the fixed $f_{a}$. $R\propto t$
adopted in [14] is the convenient choice for the scale factor during the fat-
string pre-evolution as it ensures the constant Hubble length with respect to
the string core length which we can freely choose. The random field
configurations were generated such that the Fourier amplitude
$\tilde{\phi}(\vec{k})$ below the momentum cutoff $k<k_{\text{max}}=m_{r}$
follows the standard normal distributions. The initial configurations are
rapidly adjusted via the equation of motion in the broken phase. Similarly to
[14], the pre-evolution stops when the requested number of strings per Hubble
length, denoted by $\xi_{0}$, is achieved.

Another approach is to evolve the equation of motion from the epoch right
before the phase transition where the symmetry was restored (see Appendix B.5
for the detail). It is realized by adding the temperature dependent term to
the potential [10, 15] (and accordingly modifying the equation of motion),
$\Delta\mathcal{L}=-\frac{m_{r}^{2}}{6f_{a}^{2}}T^{2}|\phi|^{2}$ where
$T^{4}=f_{a}^{4}\cdot\frac{H^{2}}{m_{r}^{2}}\cdot 4\zeta^{2}$ with
$\zeta^{2}=\frac{45}{2\pi^{2}g_{*}}\frac{M_{p}^{2}m_{r}^{2}}{f_{a}^{4}}$ and
the Hubble parameter $H^{2}=\frac{\pi^{2}}{90}g_{*}\frac{T^{4}}{M_{p}^{2}}$.
The correlation length over the core size at the critical temperature $T_{c}$
($=\sqrt{3}f_{a}$) is given by
$\left.\frac{m_{r}}{H}\right|_{T_{c}}=\frac{2}{3}\zeta$. The choice
$\zeta=\frac{3}{\sqrt{2}}=2.12123$ corresponds to the case with
$H|_{T_{c}}=f_{a}$. While $\zeta\sim\mathcal{O}(1)$ ensures the order one size
of the correlation length in the unit of the core length at the phase
transition where cosmic strings are formed, the natural size of it is roughly
$\zeta\sim 10^{7}$ for $g_{*}\sim\mathcal{O}(10^{2})$, the order one size of
the quartic coupling, and $f_{a}\sim\mathcal{O}(10^{10}\ \text{GeV})$ (see Eq.
(43) for the explicit expression). The order one size of $\zeta$ may
correspond to the large symmetry breaking scale of $f_{a}\sim 10^{17}$ GeV, a
large number of relativistic degrees of freedom of $g_{*}\sim 10^{14}$, or the
combination of them for the fixed order one quartic coupling. This might be
considered as a drawback of this approach. The random initial conditions for
field configurations were generated, assuming the Gaussian random field
configuration following the thermal distributions as in [10] in the symmetric
phase. We impose the cut on the momentum $k<k_{\text{max}}=10m_{r}$ at the
initial time $\tau_{i}=0.1\tau_{c}$ ($\tau_{c}$ as the time at the phase
transition). Since $k_{\text{max}}\propto\tau^{-1}$, $k|_{\tau_{i}}\leq
10\,m_{r}$ implies $k|_{\tau_{c}}\leq m_{r}$ with the choice of
$\tau_{i}=0.1\tau_{c}$, and this looks consistent with the cut imposed in the
aforementioned fat-string type relaxation. Unlike the approach exploiting the
fat-string technique, the generated initial configurations smoothly transit to
those of the broken phase via the equation of motion with the thermal
potential naturally realizing the phase transition. More details regarding
initial conditions are given in Appendix B.6.

The higher modes above $m_{r}$ scale can be considered as a short distance or
ultraviolet (UV) physics smaller than the string core size. Having the
momentum cutoff imposed, the simulation results become insensitive to $\Delta
x$ or equivalently the lattice size $N$ (for instance, see Fig. 20).
Otherwise, the higher momentum modes set by a smaller $\Delta x$ for a larger
$N$ kick in, the spectrum becomes more noisy, and it takes longer time to be
relaxed. The impact of UV modes above $m_{r}$ on initial conditions in the
properties of the string network will be discussed in the Appendix C.

##  4 String identification on tetrahedron

A string is basically identified by examining the pattern of phases of the
complex scalar fields. While there are various options for string
identification algorithm, we yet introduce a new novel string identification
algorithm which we find numerically very efficient and guarantees the
connectedness of strings. It starts with the tetrahedralization, sub-dividing
the lattice space into a series of tetrahedrons.

Figure 1: Illustration of a cosmic string configuration in the tetrahedron
(left) and in a cube consisting of six plaquettes (right).

Each tetrahedron, illustrated in the left panel of Fig. 1, has four triangular
faces, and each triangular face has three vertices which we can label as
$v_{i=1,2,3}$. Let us define
$\theta_{123}=\theta_{12}+\theta_{23}+\theta_{31}$ where
$\theta_{ij}=\arg(\phi_{j}\phi^{*}_{i})$ which measures the phase difference
of fields at two adjacent vertices. A string core on the triangular face is
declared for the configuration with $\theta_{123}=2\pi$ ($-2\pi$). This last
statement for the string declaration can be replaced with more numerically
efficient and equivalent statement, as was explained in Appendix B.7. With the
assignment of $\omega_{ij}=\Re\phi_{i}\Im\phi_{j}-\Im\phi_{i}\Re\phi_{j}$ and
$\omega_{123}\equiv\omega_{12}+\omega_{23}+\omega_{31}$, $\theta_{123}=2\pi$
($-2\pi$) can be shown to be equivalent to the case where all $\omega_{ij}$’s
are non-negative (non-positive) while $\omega_{123}$ is positive (negative):

$\begin{split}\theta_{123}=2\pi&\leftrightarrow\omega_{ij}\geq
0\quad\text{and}\quad\omega_{123}>0\leavevmode\nobreak\ ,\\\\[3.0pt]
\theta_{123}=-2\pi&\leftrightarrow\omega_{ij}\leq
0\quad\text{and}\quad\omega_{123}<0\leavevmode\nobreak\ .\end{split}$ (4)

The location of the string core inside the triangle $r_{\text{core}}$ is
assigned via the linear interpolation, or
$r_{\text{core}}=(\omega_{23}v_{1}+\omega_{31}v_{2}+\omega_{12}v_{3})/\omega_{123}$
such that the field value vanishes at the core. This may be considered as a
leading approximation of more complete expansion in higher orders of vertices.
In this prescription, the string core is strictly inside the triangle. This
approach is well-defined as each tetrahedron includes either none or a pair of
faces on which string cores are declared. This property also proves that
strings are always connected without any discontinuity. See [35, 21, 22] for
other types of string identifications that ensure the connectedness of
strings. In a typical string identification, each cubic cell is associated
with one volume and three faces (plaquette) 222For $n$ cells in the simulation
box, there are total $6^{n}$ faces and each face is shared by two adjacent
cells.. Upon tetrahedralization of the volume out of existing lattice points
as vertices, each cubic cell contains 6 tetrahedrons and it is associated only
with 12 triangular faces (only 4 times more). Although there are multiple
choices of assigning 6 tetrahedrons, it will be interesting if this method can
improve the spatial resolution out of only existing vertices. While our
approach has more faces to examine, our criteria based on the signs of
$\omega_{ij}$ for the string core declaration and the prescription for the
string core location save computational times a lot. Overall, our approach for
the string identification takes a similar amount of the computation time
compared to others in literature.

##  5 Scaling regime

A pair of string cores found on a tetrahedron forms a line segment, and
connecting them forms the cosmic string. In our simulation, the total string
length in the simulation box is obtained by summing over all line segments.
The average number of strings per Hubble patch is measured by

$\xi=\frac{\ell_{\text{tot}}(L)t^{2}}{L^{3}}\leavevmode\nobreak\ ,$ (5)

where $L$ is the size of the simulation box and $\ell_{\text{tot}}(L)$ is the
total string length inside the box. The maximum dynamic time range is limited
by $\log\frac{m_{r}}{H}\leq\log\frac{N}{n_{c}n_{H}}$. Early simulations have
shown the strong hints that the cosmic string network approaches the scaling
regime where the reduction of $\xi$ due to the radiation is balanced with the
enhancing $\xi$ by more strings entering into the Hubble volume over time.
Once the string network is believed to flow into the scaling regime at some
point, its property will have to be insensitive to the initial configurations
and the relaxation types. However, apparently the dynamic time range limited
by the lattice size is not long enough to reach the asymptote that we expect
it to match the scaling solution, and choosing the benchmark parameter sets
for the further analyses is rather subjective. We performed a large number of
simulations on the lattice of $N^{3}=1024^{3}$ for various initial conditions
and two different relaxation schemes explained in Section 3.2, and outcomes
are illustrated in Fig. 2.

Figure 2: The evolutions of the averaged $\xi$ in the logarithmic scales from
simulations on the lattice of $N^{3}=1024^{3}$ with various initial conditions
(distinguished by different curves within the same panel) made using fat-
string pre-evolutions (left) and thermal pre-evolutions (right).

While the curves of $\xi$ in the logarithmic time $\log\frac{m_{r}}{H}$ in
Fig. 2 clearly show the strong evidence for the scaling regime (also looks
consistent with the result in [9, 14]), they demonstrate that the residual
dependency on initial conditions and relaxation types is not settled down
within the time coverage. It takes different amounts of time to be relaxed to
a relatively clean level, especially, the $\xi$ evolution started in the
finite temperature, in the right panel of Fig. 2, has relatively short time
range for a reliable fit 333We empirically observe that the curves with this
type of the relaxation requires a larger ensemble size for the smooth curves
than the case with the fat-string pre-evolution..

Figure 3: The evolutions of the averaged $\xi$ in the logarithmic scales from
simulations on the lattice of $N^{3}=4096^{3}$ for two relaxations whose
initial conditions were prepared using fat-string pre-evolutions (left) and
thermal pre-evolutions (right). The width of each curve represents the
statistical errors.

For the physical simulation whose relaxed initial conditions were prepared by
the fat-string pre-evolution, following the suggestion in [14], we similarly
choose $\xi_{0}=0.2$ at the initial time of the physical string simulation as
our benchmark curve for the bigger-sized lattice simulation. This curve with
$\xi_{0}=0.2$ can be considered in between over-generated and under-generated
string networks, and thus it is hoped to be close to the scaling solution. For
simulations starting right before the epoch of the phase transition in finite
temperature, the curve corresponding to $\zeta=2.12$ is chosen as our
benchmark. While looking at curves in the right panel of Fig. 2, the one with
$\zeta=1.5$ may be considered to be more consistent with the selection
criteria for the benchmark curve in terms of over- and under-produced string
networks, $\zeta=2.12$ can provide with more consistent comparison with the
results in literature. The $\xi$ curves for two benchmark points from the
simulation on the lattice of $N^{3}=4096^{3}$ are illustrated in Fig. 3. While
it clearly shows the strong dependency on the relaxation methods for
$\log\frac{m_{r}}{H}<5$ (also supported by the results in Fig. 2), $\xi$ from
two benchmarks exhibit similar scaling behaviors for $\log\frac{m_{r}}{H}>5$.

If the long-term scaling behavior is assumed to be characterized by the single
linear logarithmic term to which all trajectories with different initial
conditions are attracted over time, deviations from the linear term in early
times may be parametrized in terms of various mixtures of $\log^{-n}$ ($n\geq
1$) terms. We fit the $\xi$ curves with the ansatz including up to $\log^{-2}$
[14],

$\xi=\frac{d_{2}}{\log^{2}\frac{m_{r}}{H}}+\frac{d_{1}}{\log\frac{m_{r}}{H}}+c_{0}+c_{1}\log\frac{m_{r}}{H}\leavevmode\nobreak\
.$ (6)

We refer this to global log-scaling hypothesis. Along this line of reasoning,
all curves are fitted together by taking $c_{0}$ and $c_{1}$ as the global
parameters while taking $d_{1}$ and $d_{2}$ as local ones which are specific
to each curve. Variations of $c_{0}$ and $c_{1}$ when adding $\log^{-3}$ as
well to Eq. (6) or when changing the fit interval (two choices, $\log\geq 4.5$
and $\log\geq 5$ are selected) may be taken as the part of the systematic
error. The detail of the fit procedure is presented in Appendix C.1. While we
find that fitted curves are sensitive to both fit ansatz and range used for
the fit, our results from the fit over the range $\log>5$ are presented in
Table 1 444Enlarging the interval to $\log>4.5$ gives $c_{1}=0.26$, and
including $\log^{-3}$ term, while keeping the same interval of $\log>5$, leads
to $c_{1}=0.23$ for the simulation with the fat-string pre-evolution. As the
scaling regime occurs at a bit later times in the simulation using the thermal
pre-evolution, a similar exercise may not be meaningful..

Pre-evolution type | Fit with global log-scaling hypothesis | Interval for fit
---|---|---
Fat-string |  | $\xi\sim-0.81+0.21\log\frac{m_{r}}{H}$ | $\log\frac{m_{r}}{H}=[5.0,\,-]$
Thermal |  | $\xi\sim-1.15+0.26\log\frac{m_{r}}{H}$ | $\log\frac{m_{r}}{H}=[5.0,\,-]$

Table 1: The fitted curve taking the global log-hypothesis from two benchmark
simulations differing by the relaxation schemes. The combined data from
simulations on the grids of $N^{3}=1024^{3}$ and $4096^{3}$ (results in Figs.
2 and 3) were used in the fit.

###  5.1 Scaling of inter-string distances

The average inter-string distance can be measured by
$d_{\text{str.}}=(L^{3}/\ell_{\text{tot}}(L))^{1/2}$. Using the definition for
$\xi$ in Eq. (5), the inter-string distance is expected to scale like
$d_{\text{str.}}=1/(2H\sqrt{\xi})$. Taking the log-scaling hypothesis, the
simulation data is equivalently fitted with $d_{\text{str.}}=t/\sqrt{\xi}$
where $\xi$ is given by Eq. (6), and $d_{\text{str.}}\propto t/\sqrt{\log t}$
is expected at late times. On contrary, the scenario favoring the constant
scaling of $\xi$ at late times is realized as the linear scaling in time for
the inter-string distances, or $d_{\text{str.}}\propto x_{*}t$ with the
constant $x_{*}$. If the long-term behavior of $d_{\text{str.}}$ is assumed to
be attracted to a global linear term in time $t$, possible deviations from the
linear term in early times may be accounted by extra time-dependent terms. We
refer this to global constant-scaling hypothesis as it implies that
$\xi=1/(4H^{2}d^{2}_{\text{str.}})$ asymptotes to a constant. The evolution of
$d_{\text{str.}}$ in the $m_{r}^{-1}$ unit is illustrated in Fig. 4. As is
evident in Fig. 4, the curves of $d_{\text{str.}}$ from different initial
conditions do not look possible for them to be attracted toward the single
linear term in $t$ at late times. While fitting by taking all parameters as
local ones for each curve can be carried out, it looks challenging to
distinguish between the log-scaling and constant-scaling hypotheses since the
evolution of $d_{\text{str.}}$ at late times is dominated by the $t$ term.
Instead, the quantity $2Hd_{\text{str.}}$, measured by evaluating
$2H(L^{3}/\ell_{\text{tot}}(L))^{1/2}$, may be better choice to isolate the
scaling behavior of the slope of $d_{\text{str.}}$ itself, or $1/\sqrt{\xi}$.
Its evolution is illustrated in Fig. 5. Although $m_{r}d_{\text{str.}}$ at
late times looks like following almost straight lines in time $t$, the
asymptote behavior of $2Hd_{\text{str.}}$ in Fig. 5 does not support the
constant scaling behavior (as it should be obvious from our observed
logarithmic growth behavior of $\xi$).

Figure 4: The evolutions of the averaged $m_{r}d_{\text{str.}}$ from
simulations on the lattice of $N^{3}=1024^{3}$ and $4096^{3}$. Initial
conditions were prepared using fat-string pre-evolutions (left) and thermal
pre-evolutions (right). The width of each curve represents the statistical
errors.

Figure 5: The evolutions of the averaged $2Hd_{\text{str.}}$ from simulations
on the lattice of $N^{3}=1024^{3}$ and $4096^{3}$ for relaxations. Initial
conditions were prepared using fat-string pre-evolutions (left) and thermal
pre-evolutions (right). The width of each curve represents the statistical
errors.

##  6 Cosmological evolution of string network

###  6.1 Energy density

The cosmological evolutions of the energy densities of individual components,
namely axions, radial modes, and strings, demonstrate the interplay among
those over time, and provide the dynamic picture of how the string network
approaches the scaling regime. The total Hamiltonian density of the PQ field
$\phi$ is given by
$T_{00}={\dot{\phi}}^{*}{\dot{\phi}}+\frac{1}{R^{2}}\nabla\phi^{*}\cdot\nabla\phi+V$
where $\nabla$ is the gradient with respect to the comoving coordinate. It can
be expressed explicitly in terms of the radial mode $r$, axion mode $a$, and
their interactions,

$T_{00}=\frac{1}{2}\left(1+\frac{r}{f_{a}}\right)^{2}\left(\dot{a}^{2}+\frac{1}{R^{2}}\left|\nabla
a\right|^{2}\right)+\frac{1}{2}\left(\dot{r}^{2}+\frac{1}{R^{2}}\left|\nabla
r\right|^{2}\right)+V\leavevmode\nobreak\ .$ (7)

The total energy density is the spatial average of the total Hamiltonian
$\rho_{\text{tot}}=\langle T_{00}\rangle$ over the physical coordinates. The
energy densities of axions and radial modes are measured by evaluating

$\begin{split}\rho_{a}&=2\times\langle\frac{1}{2}\dot{a}^{2}\rangle\leavevmode\nobreak\
,\\\\[3.0pt]
\rho_{r}&=\langle\frac{1}{2}\left(\dot{r}^{2}+\frac{1}{R^{2}}\left|\nabla
r\right|^{2}\right)+V\rangle\leavevmode\nobreak\ ,\end{split}$ (8)

where fields are properly masked according to the prescription described in
Appendix B.8 to screen off the effect from string cores. Assuming that the
total energy density of the PQ field is conserved, the energy density of
strings is extracted by subtracting $\rho_{a}$ and $\rho_{r}$ in Eq. (8) from
the total energy density, $\rho_{s}=\rho_{\text{tot}}-\rho_{a}-\rho_{r}$. The
gradient of axions in Eq. (7) mainly contributes to the string tension and it
is taken as part of the string energy density.

Figure 6: Energy budgets of axions, radial modes, and strings from the
simulations with $N^{3}=4096^{3}$ (top) and also with smaller lattice size
(bottom) whose initial data was relaxed by the fat-string pre-evolution (left)
and thermal pre-evolution (right).

Energy densities of individual components for two different simulation setups
are illustrated in Fig. 6. In the region before the string network enters the
scaling regime, the separation between different components should be
ambiguous, and their physical meanings become more transparent at late times.
Contrary to similar scaling behaviors in $\xi$ curves in Fig. 3 for
$\log\frac{m_{r}}{H}>5$ for both types of pre-evolutions, the apparent
evolution of the energy density for each component look dramatically
different. Nevertheless, all energy densities for axions, strings, and radial
modes near the end of simulation amazingly agree well. Since two simulations
in Fig. 6 differ by the way the initial field configurations were generated
and subsequent relaxation schemes, the distinctive behaviors until late times
may indicate a source of a large systematic errors or an ambiguity in
selecting the would-be attractor solution. Unlike the energy budget in [9]
where the different prescription for the fat-string pre-evolution was used
555The simulation in [9] implements the fat-string pre-evolution with the
scale factor scaling as $R\propto\sqrt{t}$ whereas, in our prescription, the
scale factor $R\propto t$ as in [14] is used during the pre-evolution., the
contributions from strings and axions in our simulation cross around $\log\sim
5.5$, as is seen in the left panels of Fig. 6, and the axion-to-string energy
density ratio seems to follow the expectation in the scaling regime,
$\rho_{a}/\rho_{s}\sim\log\frac{m_{r}}{H}$. The contribution from radial modes
is supposed to decrease at late times whereas $\rho_{r}$ in Fig. 6 shows the
arising behavior around $\log\frac{m_{r}}{H}\sim 6$ in both benchmark
simulations. The similar observation was made in [14] except that the overall
rate is smaller, for instance, $\rho_{r}\sim 20$% in the scaling regime in
Fig. 6 versus $\sim 14$% in [14]. The latter [14] compares better with our
simulation with $\xi_{0}\sim 0.1$ as is indicated in the left of Fig. 21 in
Appendix C.2.1 where the radial modes appear taking roughly 14% of the total
energy density. However, the crossing point between strings and axions for the
case with $\xi_{0}=0.1$ is delayed to more late times, around
$\log\frac{m_{r}}{H}\sim 6.5$. The bottom-left panel of Fig. 6 overlays three
simulation results only differing by the lattice size which equivalently can
be taken to extract the dependency on the lattice spacing in the string core
length unit $m_{r}\Delta$ 666Since the bigger lattice size $N$ is associated
with the smaller lattice spacing $\Delta$, increasing $N$ is equivalent to
taking the continuum limit from the smaller-sized-lattice point of view at a
fixed time., and it strongly indicates that the arising behavior after around
$\log\frac{m_{r}}{H}\sim 6$ could be an artifact due to the finite size of the
lattice spacing. A similar conclusion from the simulation using the thermal
pre-evolution is rather difficult due to undamped noises as is evident in the
bottom-right panel of Fig. 6.

In the right panels of Fig. 6, the energy budget from simulations using
thermal pre-evolution is presented in a similar manner. While energy densities
from the simulation with $N^{3}=4096^{3}$ was coarsely recorded starting only
from $\log\frac{m_{r}}{H}\sim 2.3$ to save the computational time, more dense
time samplings starting from earlier times were selected to record the data in
the simulations with smaller lattice sizes. While the bottom-right panel of
Fig. 6 shows apparently different evolution patterns of axions and strings
compared to those using fat-string pre-evolution, the energy budgets of two
simulations using the fat-string and thermal pre-evolutions amazingly agree
well at the final time, $\log\frac{m_{r}}{H}\sim 7.8$, as was mentioned above.
This agreement may justify the choice of two benchmark simulations on the
lattice with $N^{3}=4096^{3}$. The analogous dependency on the lattice spacing
$m_{r}\Delta$ does not have a clear pattern, for instance, the slightly
decreasing $\rho_{a}$ (and slightly increasing $\rho_{s}$) near their final
times in simulations with $N^{3}=1024^{3}$, $2048^{3}$ is flipped when the
dynamic time range is extended with the bigger lattice size. The overall
relaxation in the simulation using the thermal pre-evolution seems to take
longer dynamic time and is subject to the larger uncertainty than the case
using the fat-string pre-evolution.

###  6.2 String tension and boost factor

Figure 7: The cosmological evolution of the string tension with respect to the
theory expectation from two benchmark simulations whose initial field
configurations were prepared by fat-string pre-evolution (left) and thermal
pre-evolution (right).

The total energy stored in strings $E_{s}$ is given by
$E_{s}=\mu_{\text{eff}}\,\ell_{\text{tot}}(L)$ where $\mu_{\text{eff}}$ is an
effective tension and $\ell_{\text{tot}}(L)$ is the total string length in the
simulation box of size $L$. The effective tension can be extracted from
measurements of $\xi$ and the string energy density $\rho_{s}=E_{s}/L^{3}$.
That is, using the relation in Eq. (5), $\mu_{\text{eff}}$ can be expressed as

$\mu_{\text{eff}}=\frac{E_{s}}{\ell_{\text{tot}}(L)}=\rho_{s}\frac{t^{2}}{\xi}\leavevmode\nobreak\
,$ (9)

or the energy density of strings is given by

$\rho_{s}=\mu_{\text{eff}}\frac{\xi}{t^{2}}\leavevmode\nobreak\ .$ (10)

The string tension of the global strings in Eq. (9) is logarithmically
divergent. Its divergence is cut off by the UV cutoff, $r_{\text{UV}}$,
roughly string core length $m_{r}^{-1}$, and the inter-string distance, which
is roughly $(L^{3}/\ell_{\text{tot}}(L))^{1/2}$, and it leads to the form that
suits well for the straight string [34],

$\mu^{\text{th}}_{\text{eff}}=\langle\gamma\rangle\mu_{0}\log\frac{\left(L^{3}/\ell_{\text{tot}}(L)\right)^{1/2}}{r_{\text{UV}}}=\langle\gamma\rangle\mu_{0}\log\frac{1}{2r_{\text{UV}}H\xi^{1/2}}\leavevmode\nobreak\
,$ (11)

where $\langle\gamma\rangle$ is the averaged (over the length) boost factor of
string velocity, and $\mu_{0}$ denotes the proportionality factor. On
contrary, the logarithmic divergence is absent in the local strings [24, 41,
42, 45, 44] due to the compensating contribution from the gauge fields. The
logarithmic dependence on distances in global strings indicates the long-range
force mediated by massless axions between strings. However, as the coupling of
large wavelength axions with strings scales as $1/\log(m_{r}/H)$ [46], the
string network at late times is expected to become similar to the situation of
the local strings. By comparing the measured effective string tension in Eq.
(9) against the theoretical expectation in Eq. (11), the non-trivial scaling
property of the string energy density can be checked. The effective tension
$\mu_{\text{eff}}$ with respect to $\mu_{\text{eff}}^{\text{th}}$ is
illustrated in Fig. 7 where $\mu_{\text{eff}}^{\text{th}}$ is given by Eq.
(11) with $\langle\gamma\rangle\mu_{0}=\pi f_{a}^{2}$ and
$r_{\text{UV}}=m_{r}^{-1}(4\pi)^{1/2}/2$ [14]. As is evident in Fig. 7, the
effective string tension-to-theoretical expectation ratio is relaxed to the
roughly constant value at later times, and it implies that the expression in
Eq. (11) is intriguingly applicable to more generic configurations.
Approaching the scaling regime looks more dramatic in the simulation with the
thermal pre-evolution as is seen in the right panel of Fig. 7.

Figure 8: The probabilistic distribution $P(\gamma)$ of the boost factor
$\gamma$ from two benchmark simulations on the grid with $N^{3}=4096^{3}$.

The time variation of $\mu^{\text{th}}_{\text{eff}}$ in Eq. (11) is given by

$\dot{\mu^{\text{th}}_{\text{eff}}}=\frac{\dot{\langle\gamma\rangle}}{\langle\gamma\rangle}\mu^{\text{th}}_{\text{eff}}+\langle\gamma\rangle\mu_{0}\left(\frac{1}{t}-\frac{\dot{\xi}}{2\xi}\right)\leavevmode\nobreak\
.$ (12)

If the averaged boost factor $\langle\gamma\rangle$ approaches the scaling
regime, the prefactor $\langle\gamma\rangle\mu_{0}$ in Eq. (11) is expected to
behave like

$\langle\gamma\rangle\mu_{0}=\dot{\mu^{\text{th}}_{\text{eff}}}\left(\frac{1}{t}-\frac{\dot{\xi}}{2\xi}\right)^{-1}\leavevmode\nobreak\
,$ (13)

which can be explicitly checked from lattice simulations. The average boost
factor is obtained by evaluating

Figure 9: The averaged boost factor $\langle\gamma\rangle$ from two benchmark
simulations on grids with $N^{3}=2048^{3}$ and $4096^{3}$.

$\langle\gamma\rangle=\int_{1}^{\infty}\gamma
P(\gamma)d\gamma=\int_{1}^{\infty}\gamma\
\frac{1}{\xi}\frac{d\xi_{\gamma}}{d\gamma}d\gamma\leavevmode\nobreak\ ,$ (14)

where $\xi_{\gamma}$ is the fraction of $\xi$ with the boost factor below
$\gamma$. The Lorentz boost factor $\gamma$ of the string core is computed
from field configurations through the relation
$\gamma^{2}-1=\beta^{-2}|\dot{\phi}|^{2}$ evaluated at every string core.
Recall that, in our tetrahedron-based string identification method, a pair of
string core locations are assigned via the linear interpolation from field
values at three vertices on each triangular face. Similarly for $\dot{\phi}$
at the string core. The corresponding tetrahedron includes string line segment
(if exists), and the string boost factor $\gamma$ is evaluated at the mid-
point of the string line segment. The probability distribution $P(\gamma)$ of
the boost factor is illustrated in Fig. 8. $P(\gamma)$ roughly scales like
$\sim\gamma^{-6}$ as was also observed in [14]. Finally, the time variation of
the averaged boost factor is illustrated in Fig. 9. To save the computational
time, the boost factors were computed only at the coarsely selected sampling
points. In the simulation with the lattice size of $N^{3}=4096^{3}$, only time
slices after $\log\sim 5.5$ were picked. The earlier times were covered by the
simulation on the smaller lattice of $N^{3}=2048^{3}$ with the small overlap
with those from the former. The overall qualitative behavior and averaged
boost factor $\langle\gamma\rangle=1.2\div 1.3$ agrees well with the result in
[14] as is seen in the left panel of Fig. 9. On contrary, the early time
behavior appears pronouncedly different for the cosmic strings in the
simulation using the thermal pre-evolution although both become similar at
late times. The gray dashed lines in both panels in Fig. 9 represent the
expectation from the ansatz $P(\gamma)=5\gamma^{-6}$. Two curves from
simulations differing by the lattice sizes, at an instant time, can be
effectively thought of as two cases differing by $m_{r}\Delta$, for instance,
$m_{r}\Delta$ of the red curve at a fixed time within $\log=[5.5,\,7]$ will be
larger than that of the blue curve which is closer to the continuum limit.

##  7 Axion spectrum

Figure 10: Fraction of the transfer rate from strings to axions,
$r_{a}=\Gamma_{a}/(\Gamma_{a}+\Gamma_{r})$, from two benchmark simulations
using the fat-string ($\xi_{0}=0.2$) and thermal pre-evolutions
($\zeta=2.12$).

The strings lose their energy into axions and radial modes, denoted by
$\Gamma$. Individual transfer rates into axions and radial modes are
represented by $\Gamma_{a}$ and $\Gamma_{r}$, respectively, and
$\Gamma=\Gamma_{a}+\Gamma_{r}$. Similarly to [9], our simulation confirms that
strings mostly decay into axions with small fraction to radial modes, namely
$\Gamma\sim\Gamma_{a}$, as is illustrated in Fig. 10. The fluctuating
evolution of $r_{a}$ from the simulation using the thermal pre-evolution (red
solid line in Fig. 10) is due to the oscillatory behavior in the energy budget
as was illustrated in the right panels in Fig. 6. The evolution of the
relativistic axion energy density follows the equation

$\frac{1}{R^{4}(t)}\frac{\partial}{\partial
t}\left(R^{4}(t)\rho_{a}(t)\right)=\Gamma_{a}+\cdots$ (15)

where $\cdots$ corresponds to the effect from the interaction between axions
and radial modes. Similarly for the radial energy density except that it
scales like matter. Assuming that $\Gamma_{a}$ eventually dominates at late
times, $\Gamma\sim\Gamma_{a}$, in the evolution of $\rho_{a}$, the axion
energy density can be expressed as the integration of the transfer rate,

$\rho_{a}(t)=\int^{t}dt^{\prime}\left(\frac{R(t^{\prime})}{R(t)}\right)^{4}\,\Gamma(t^{\prime})\leavevmode\nobreak\
.$ (16)

The inclusive axion energy density does not tell about the characteristic
hardness or softness of the axions radiated from strings, and it is more
informative to look into the differential distribution of the axion energy
density in momentum $\partial\rho_{a}/\partial k$, defined by $\rho_{a}=\int
dk\,(\partial\rho_{a}/\partial k)$. In the absence of no new scales roughly
between $m_{r}$ and $H$, the existence of the scaling regime indicates the
power law fall-off behavior, $\sim k^{-q}$. The characteristic feature of the
axion spectrum falls into roughly three categories depending on the value of
$q$ [9]. When $q>1$ [47, 48, 26, 27], the spectrum is more pronouncedly peaked
near the order of Hubble which implies that the spectrum is IR-dominated and
axions radiated from strings are soft. On contrary, when $q<1$, the axion
spectrum is relatively suppressed near the Hubble scale and more pronounced
near the string core scale. Axions in this category are supposed to be hard
and the spectrum is UV-dominated. The situation with $q=1$ [49, 50] is in
between the above two cases.

Figure 11: The spectral axion energy density from simulations with
$N^{3}=4096^{3}$ for two benchmark points. Dashed lines represent $m_{r}/(2H)$
for three choices of time snapshots.

The distributions of the differential axion energy density, divided by
$Hf_{a}^{2}$ to cancel the overall time dependence, in two different
simulation setups, differing by the relaxation scheme, on the lattice with
$N^{3}=4096^{3}$ are illustrated in Fig. 11 at three instant times, $\log=6$
(purple), 7 (blue), 7.5 (red), respectively. The positions of $m_{r}/(2H)$,
where the parametric resonances are expected, are also presented by dashed
lines with the same color coding. The shape of $\partial\rho_{a}/\partial k$
at each instant time in Fig. 11 nicely captures the power law fall-off scaling
behavior between two scales $m_{r}$ and $H$, and the distributions rapidly
drop off at the momentum below the Hubble scale $H$ and above the string core
scale $m_{r}$ as is expected. Specifically, from the curves in Fig. 11, the
peaks in IR region are located at $k/H\sim 5\div 10$ and the drop-offs in UV
occur around $m_{r}/(2H)$.

From the transfer rate from strings to axions in the differential form
$\Gamma(t)=\int dk\,(\partial\Gamma(t)/\partial k)$, assuming
$\Gamma\sim\Gamma_{a}$, the normalized instantaneous emission function $F$ can
be defined as

$1=\int\frac{dk}{H(t)}\frac{H(t)}{\Gamma(t)}\frac{\partial\Gamma}{\partial
k}(k,t)\equiv\int\,dxF[x,y]\leavevmode\nobreak\ ,$ (17)

where $x=k/H$ and $y=m_{r}/H$. The instantaneous emission $F$ is nothing but
the normalized differential transfer rate measured in $k/H$, and it
characterizes the distribution of the axion spectrum in momentum, contributing
to the axion energy density in Eq. (16). With the choice of $k/H$, the peak
location at order of the Hubble remains constant. In the simulation, $F$ is
measured from the axion spectral energy density $\partial\rho_{a}/\partial k$,

$F\left[\frac{k}{H},\,\frac{m_{r}}{H}\right]=\frac{1}{\Gamma_{a}/H}\frac{1}{R^{3}}\frac{\partial}{\partial
t}\left(R^{3}\frac{\partial\rho_{a}}{\partial k}\right)\leavevmode\nobreak\ ,$
(18)

by taking two time slices of the data, for instance, separated by
$\log\frac{m_{r}}{H}=0.25$ as our default choice. The time intervals used for
the instantaneous emission spectrum differ in literature, and its dependency
on the time interval is discussed in Appendix D.1 in detail.

Figure 12: The spontaneous emission $F$ from simulations on the lattice of
$N^{3}=4096^{3}$ for $\xi_{0}=0.2$ (left) and $\zeta=2.12$ (right) at three
instant times snapshots, $\log\frac{m_{r}}{H}=$ 6, 7, 7.5. Dashed vertical
bars with the same color coding correspond to $m_{r}/(2H)$ at each time
instant.

The instantaneous emission $F$ from our simulations on the lattice of
$N^{3}=4096^{3}$ for two simulation setups are illustrated in Fig. 12 for
three logarithmic time slices, $\log=6$, 7, 7.5. While the overall shapes in
both panels of Fig. 12 are similar, the emission spectra using the thermal
pre-evolution are more noisy although they become cleaner at late times.

###  7.1 Power law scaling

The power law fall-off behavior of $F\propto k^{-q}$ between $m_{r}$ and $H$
is the most crucial part of the cosmic string network for the plausibility of
the axion as the dark matter candidate. The precise determination of the
spectral index $q$ will allow us to estimate the relic abundance of axions. To
avoid contaminations from the UV and Infrared (IR) regions, the fit needs to
be done inside an appropriate interval away from the string core and Hubble
scales, or parametrized as $k=[x_{\text{IR}}H,\,m_{r}/x_{\text{UV}}]$. Since
there is no obvious choice for the interval for fitting $F$ and it rather
depends on the situation, we fit the power law fall-off region for various
choices of $x_{\text{UV}}$, $x_{\text{IR}}$ to screen off UV and IR
contaminations. It is instructive to estimate the number of sampling points
within the interval for given $x_{\text{IR}}$, $x_{\text{UV}}$, and the
lattice size $N$. The physical momentum on the lattice is
$\vec{k}=\frac{2\pi}{RN\Delta x}\vec{n}_{k}$ where the index runs up to
$\frac{N}{2}$ 777It should be distinguished from the physical momentum $k$
truncated at $m_{r}$ when generating initial conditions for scalar fields. The
maximum momentum set by the lattice spacing is
$k_{\text{max}}=\frac{2\pi}{R(2\Delta x)}$., and in the Hubble unit,
$x=k/H=\frac{2^{3/2}\pi}{N(\Delta xm_{r})}\sqrt{y}\,n_{k}$. The number of
sampling points used for the fit within the interval $x_{\text{IR}}\leq x\leq
y/x_{\text{UV}}$ for the static lattice simulation, using $\Delta x$ in Eq.
(29), is given by

$N_{\text{fit}}=\frac{\text{Interval of
}x}{x\text{-spacing}}=\frac{1}{2^{3/2}\pi}\left(\frac{\sqrt{y}}{x_{\text{UV}}}-\frac{x_{\text{IR}}}{\sqrt{y}}\right)\left(\frac{2n_{H}N}{n_{c}}\right)^{1/2}\leavevmode\nobreak\
.$ (19)

Figure 13: The number of sampling points inside $x_{\text{IR}}\leq
x\leq\frac{y}{x_{\text{UV}}}$ for fitting the power law for the lattice size
of $2048^{3}$ and $4096^{3}$. $n_{c}=1$, $n_{H}=4^{1/3}$ were used. Curves are
plotted from 5 to the maximal time coverage $\sim\log\frac{N}{n_{c}n_{H}}$.

As is illustrated in Fig. 13 for various intervals, the number of sampling
points increases with the increasing time. It is because the spacing of $x$
scales as $\sim\sqrt{y}$ with the time whereas the upper end of the interval
scales as $\sim y$ for the fixed lower end value. At later times, both higher
and lower $k$ modes previously staying outside the interval enter into the
interval. In terms of $k$ scaling as $\sim R^{-1}$, the lower end of the
interval $\sim H$ decreases faster than $R^{-1}$ and more lower $k$ modes
enter into the interval. More higher modes from above the $m_{r}$ scale get
redshifted below $m_{r}$ and part of them may also enter into the interval. As
is evident in Fig. 13, the lattice size of $N^{3}=2048^{3}$ seems to be enough
for the loose interval, $15\leq x\leq\frac{y}{4}$, adopted in [14], and not
enough for the tight interval, $40\leq x\leq\frac{y}{16}$ adopted in [15]
888When using AMR, assuming the locally refined region around string cores are
completely masked for simplicity, the lattice spacing which is relevant for
the power spectrum analysis will be similar to the static case except for the
modified lattice spacing in Eq. (32). However, as the simulation time extends
over beyond the coverage for the static simulation, more momentum modes enter
into the interval of $x$, increasing $N_{\text{fit}}$. One limitation in this
situation with AMR will be that the upper end of the interval should be
smaller than the maximally allowed $x$ set by $\Delta x$, or
$\frac{y}{x_{\text{UV}}}\leq x_{\text{max}}$ where
$x_{\text{max}}=\frac{\pi}{HR\Delta
x}=\pi\sqrt{\frac{n_{c}N}{2^{k}n_{H}}\,y}$, and it gives the upper limit on
$\log_{i}\leq\log[x_{\text{UV}}^{2}\pi^{2}(n_{c}/2^{k})\frac{N}{n_{H}}]$. For
the simulations with $N^{3}=2048^{3}$, $n_{c}=4$, $k=4$, and $n_{H}=4^{1/3}$,
the limit on the time coverage is estimated to be $\log_{i}\leq 10.8$ for
$x_{\text{UV}}=4$ and $\log_{i}\leq 13.6$ for $x_{\text{UV}}=16$..

The modes above $m_{r}$ at the phase transition correspond to the short
distance physics above the PQ symmetry breaking scale. At later times, those
UV modes such as radial modes are redshifted and can enter the interval.
Requiring no contamination from those redshifted UV modes sets the upper bound
on the $\log_{i}$, or $\frac{y}{x_{\text{UV}}}\leq\frac{m_{r}}{HR}$ which
leads to $\log_{i}\leq 2\log[\sqrt{2}\,x_{\text{UV}}]$. $\log_{i}\leq 6.2$
(3.47) for $x_{\text{UV}}=16$ (4). There will be an unavoidable contamination
unless the radial modes get decayed away before reaching the scaling regime.
Matching this upper limit to the maximal time coverage set by the lattice size
leads to $x_{\text{UV}}=\sqrt{N/(2n_{c}n_{H})}$. Assuming $n_{c}=1$ and
$n_{H}=4^{1/3}$, $x_{\text{UV}}$ needs to be set to $x_{\text{UV}}=25.4$ for
$N^{3}=2048^{3}$ and $x_{\text{UV}}=35.9$ for $N^{3}=4096^{3}$. Applying this
tight UV cutoff will make extracting the power law fall-off behavior
challenging due to the limited sampling points. However, upon generating the
initial conditions for scalar fields, only momentum modes below $m_{r}$ were
populated. Although higher momentum modes are generated via the equation of
motion, they do not look as harm as those originally generated at the
generation of the random field configurations.

The results from our simulations following the fit prescription explained in
Appendix D.1 are illustrated in Figs. 14 and 15 for our default choice of
$\Delta\log=0.25$ (see Figs. 25 and 26 for different choices of $\Delta\log$).
We observe the similar logarithmic growth to those in [14] except that our
result shows a bit steeper slope than that in [14], for instance, the value of
$q$ reaches unity within the dynamic time range of the simulation for a
similar choice of the interval. It will be interesting to see how $q$ further
evolves beyond the time coverage of our current simulation, and this justifies
improving the lattice resolution by the factor of two, namely
$N^{3}=8192^{3}$. In Fig. 14, the UV cutoff is varied while IR cutoff
$x_{\text{IR}}$ is fixed to 15. In Fig. 15, the IR cutoff is varied while UV
cutoff $x_{\text{UV}}$ is fixed to 10. The error bars in all curves represent
the statistical ones over 100 independent simulation runs. Interestingly, the
patterns of central values at late times look similar irrespective of types of
pre-evolutions. We suspect that the fluctuations of the central values are due
to effects from the residual red-shifted high-frequency modes. Based on
distributions in Figs. 14 and 15, the interval $k=[15H,\,m_{r}/10]$ looks
reasonable choice to get the fitted spectral index with a relatively small
fluctuations.

Figure 14: The spectral index $q$ of the instantaneous emission $F$ fitted in
the interval with the varying $x_{\text{UV}}$ while $x_{\text{IR}}=15$. The
error bars are statistical.

Figure 15: The spectral index $q$ of the instantaneous emission $F$ fitted in
the interval with the varying $x_{\text{IR}}$ while $x_{\text{UV}}=10$. The
error bars are statistical.

The fit result of the spectral index $q$, assuming the log-hypothesis, in the
interval $\log=[6.5,\,7.8]$ for two benchmark simulations on lattice with
$N^{3}=4096^{3}$ differing by relaxation types is presented in Table 2.

Pre-evolution type | Fit with $(x_{\text{IR}},\,x_{\text{UV}})=(15,\,10)$ | Interval for fit
---|---|---
Fat-string |  | $q\sim(-0.349\pm 0.116)+(0.174\pm 0.016)\log\frac{m_{r}}{H}$ | $\log\frac{m_{r}}{H}=[6.5,\,7.8]$
Thermal |  | $q\sim(-0.460\pm 0.079)+(0.186\pm 0.011)\log\frac{m_{r}}{H}$ | $\log\frac{m_{r}}{H}=[6.5,\,7.8]$

Table 2: The fit result taking the log-hypothesis from two benchmark
simulations on $N^{3}=4096^{3}$ differing by the relaxation schemes. The
presented errors are merely statistical ones within the particular choices of
cuts on $x_{\text{IR}}$, $x_{\text{UV}}$ and the fit intervals.

The outcome of fits with different assumptions on the form of the Likelihood
function, UV/IR cutoffs on the momentum, and time interval $\Delta\log$ are
summarized in Figs. 27 and 28. The discrepancy, compared to the result in
Table 2, may be considered as part of systematic errors. Most importantly,
despite all variant situations, our conclusion on the spectral index
predicting $q\sim O(10)$ at $\log{\frac{m_{r}}{H}}\sim 70$ does not change.

##  8 Axion abundance

The axion string network evolves until the time of the QCD crossover, denoted
by $t_{\star}$ (roughly $\log_{\star}\equiv\log\frac{m_{r}}{H(t_{\star})}\sim
70$), around which the axion mass $m_{a}$ becomes comparable with the Hubble,
or $m_{a}(t_{\star})\sim H_{\star}$ and the axion potential becomes relevant.
The impact of the logarithmically growing spectral index $q$ of the
instantaneous emission $F$, obtained from fitting the data in Section 7.1, on
the axion abundance is dramatic. Despite the UV-dominated spectrum within the
simulation time range, extrapolating the spectral index $q$ until the time of
the QCD crossover indicates the IR-dominated spectrum of axions radiated from
strings at those late times. The system of strings and domain walls are formed
around this time [51, 52]. Especially, when the number of domain walls are
one, $N_{W}=1$, the strings-and-walls system decay into axions [53] (see [54,
13, 16] for simulations). Similarly to [14], our estimate will not include
axions from the collapse of strings-and-walls system. We will focus only on
the axions radiated strings in the scaling regime and its evolution afterwards
(later times than $t_{\star}$) through the nonlinear regime in presence of the
axion potential. Therefore, it sets the conservative lower bound on the axion
abundance (thus the upper bound on the axion decay constant).

The axion energy density from strings in the scaling regime is still bigger
than the axion potential at $t=t_{\star}$ due to the dominant energy density
stored in gradient terms. It continues evolving until the time, denoted by
$t_{\ell}$, and it becomes comparable to the axion potential where axions
transit to nonrelativistic ones. If the average axion field values are smaller
than the axion decay constant, $\langle a^{2}\rangle^{1/2}\ll f_{a}$, the
nonlinearities due to the axion potential in the axion equations of motion can
be ignored and the axion abundance will be obtained simply by evaluating $\int
dkk^{-1}(\partial\rho_{a}/\partial k)$. Whereas, in the opposite situation,
$\langle a^{2}\rangle^{1/2}\gtrsim f_{a}$, the axion abundance at $t_{\ell}$
is affected by the nonlinearities and the estimation of the axion abundance is
more tricky. The parametric behavior of $n^{\text{str}}_{a}(t_{\ell})$, axion
abundance from cosmic strings, in terms of string properties takes different
forms depending on the spectral index $q$, the value of ${\langle
a^{2}\rangle}^{1/2}/f_{a}$, and the IR cutoff of the momentum. The detailed
derivation of all distinctive cases are postponed to Appendix E. Here, we will
discuss only relevant case to our situation.

The axion spectrum at the time $t=t_{\star}$ can be obtained by evaluating

$\frac{\partial\rho_{a}}{\partial
k}=\int^{t}dt^{\prime}\frac{\Gamma(t^{\prime})}{H(t^{\prime})}\left(\frac{R(t^{\prime})}{R(t)}\right)^{3}F\left[\frac{k^{\prime}}{H(t^{\prime})},\,\frac{m_{r}}{H(t^{\prime})}\right]\leavevmode\nobreak\
,$ (20)

where the form of $F$ for $q>1$ is given in Appendix E and the transfer rate
at late times is approximated by $\Gamma\sim\xi\mu_{\text{eff}}/t^{3}\sim 8\pi
H^{3}f_{a}^{2}\xi\log\frac{m_{r}}{H}$. When $q>1$, the contribution to the
axion abundance at $t=t_{\star}$ from higher momenta is power suppressed as
$\partial\rho_{a}/\partial k$ decreases faster than $\sim 1/k^{q}$, and only
low momenta roughly $k<\sqrt{m_{r}H_{\star}}$ which leads to the scaling
$(\partial\rho_{a}/\partial k)|_{t_{\star}}\propto k^{-1}$ dominantly
contributes. It has been analytically approximated in [14] (and supported by
numerical simulation) to estimate the axion abundance around the time
$t=t_{\ell}$, taking into account the nonlinearities. It was shown that the
topological production of axions with respect to the contribution from the
misalignment mechanism with the angle of unity is expected to be scaling like
$(n_{a}^{\text{str},q>1}/n_{a}^{\text{mis},\theta_{0}=1})|_{t_{\ell}}\propto(\xi_{\star}\log_{\star})^{1/2+\cdots}$
where $\xi_{\star}$ is the number of strings per Hubble patch at the time of
the QCD crossover, $\log_{\star}\sim 70$ and $\cdots$ denotes the dependency
on the axion mass evolution. The exact value of $q$ is not relevant as long as
it is well above unity at the time of the QCD phase transition. Assuming that
$q$ continuously increases logarithmically, the ballpark of the axion dark
matter mass will be similar to roughly what was estimated in [14].

As was discussed in Section 5.1, the average inter-string distance is roughly
$(H\sqrt{\xi})^{-1}$ instead of $H^{-1}$. One may worry that axions whose
momentum softer than $\sim H\sqrt{\xi}$ are suppressed and it may cause the
parametrically reduced axion abundance. In Appendix E, we have shown that it
can cause the non-negligible change of the overall size of the axion abundance
at late times due to the slightly modified parametric dependence on string
properties. It is also demonstrated in Fig. 16 for a set of selected
parameters that will be further discussed in Section 8.1. However, the
modified IR cutoff can significantly change parametric behavior if the small
average axion field values $\langle a^{2}\rangle^{1/2}\ll f_{a}$ are assumed.
It is also true for the case with $q=1$ similarly to [15] whose spectral index
was claimed to be consistent with no-log hypothesis. The parametric behaviors
of all distinctive cases are explicitly worked out in Appendix E.

###  8.1 Bound on axion mass

Figure 16: The analytic prediction of the axion abundance from strings with
respect to the misalignment with $\theta_{0}=1$ as a function of $\xi_{\star}$
for varying $x_{0}$ while $\log_{\star}=64$ and $\alpha=8$ are fixed. The
solid lines correspond to those obtained with the IR momentum cutoff of
$k_{\text{IR}}=x_{0}H\sqrt{\xi_{\star}}$ whereas the dashed ones are those
with $x_{0}H$. $c_{n}=1.35$, $c_{V}=0.13$, $c_{m}=2.08$ and
${c_{n}}^{\prime}=2.81$ were taken from [14]. The right panel is zoomed-in
version of the left one, and the color coding is the same as those in the left
panel.

The number of strings per the Hubble patch around time $\log_{\star}=64$
similarly to [14] is estimated from our simulation result, and they are
$\xi_{\star}=13.4$ (fat-string pre-evolution) and $\xi_{\star}=16.6$ (thermal
pre-evolution). Adopting fit values $c_{n}=1.35$, $c_{V}=0.13$, $c_{m}=2.08$
and ${c_{n}}^{\prime}=2.81$ from [14] and setting $\alpha=8$, the axion
abundance with respect to the misalignment mechanism with the unit angle,
given in Eq. (75) along with Eq. (76) in Appendix E, can be estimated to give
$(n_{a}^{\text{str},q>1}/n_{a}^{\text{mis},\theta_{0}=1})|_{t_{\ell}}\sim 119$
(fat-string pre-evolution) and
$(n_{a}^{\text{str},q>1}/n_{a}^{\text{mis},\theta_{0}=1})|_{t_{\ell}}\sim 134$
(thermal pre-evolution) with the choice of $x_{0}=10$ as an illustration.
While the analytic predictions as a function of $\xi_{\star}$ for the
different IR cutoff $k_{\text{IR}}$ (with varying $x_{0}$) is illustrated in
Fig. 16, more general expressions are given in Appendix E. The lower bound on
the axion mass can be derived by requiring that the axion dark matter
abundance does not exceed the currently observed value. Taking the reference
value $\Omega_{a}^{\text{mis}}h^{2}\approx
0.12\left(28\,\mu\text{eV}/m_{a}\right)^{7/6}$ in the post-inflationary
scenario as in [14] where we assume $\langle\theta_{0}^{2}\rangle\sim 5$ for
translating it to $\Omega_{a}^{\text{mis},\theta_{0}=1}h^{2}$ 999The precise
estimate of the contribution from the misalignment [55, 7, 56] amounts to
shifting the bound.,

$\begin{split}\Omega_{a}^{\text{str}}h^{2}\approx\frac{n_{a}^{\text{str},q>1}}{n_{a}^{\text{mis},\theta_{0}=1}}\Omega_{a}^{\text{mis},\theta_{0}=1}h^{2}\leq
0.12\rightarrow&\left\\{\begin{array}[]{lll}m_{a}\gtrsim
420\leavevmode\nobreak\ (\mu\text{eV})&:&\text{fat-string pre-
evolution}\\\\[15.0pt] m_{a}\gtrsim 470\leavevmode\nobreak\
(\mu\text{eV})&:&\text{thermal pre-evolution}\leavevmode\nobreak\
.\end{array}\right.\end{split}$ (21)

Using $m_{a}\approx 5.7\,\mu\text{eV}(10^{12}\,\text{GeV}/(f_{a}/N_{W}))$
[56], $f_{a}/N_{W}\leq 1.3\times 10^{10}\text{GeV}$ (fat-string pre-evolution)
and $f_{a}/N_{W}\leq 1.1\times 10^{10}\text{GeV}$ (thermal pre-evolution) for
QCD axions are derived.

##  9 Correlation between strings and axion spectrum

Figure 17: The correlation between $q$ and $d\log(\xi
H)/(d\log\frac{m_{r}}{H})$ over 100 ensemble elements from two benchmark
simulations using the fat-string type pre-evolution (left) and thermal pre-
evolution (right) on the lattice with $N^{3}=4096^{3}$.

While our computing resources were limited only to two benchmarks, it would be
interesting to explore possible observables or relations that can hint for the
behavior for the benchmark points with different initial conditions. In this
section, we investigate one such a correlation between $q$ and time variation
of $\xi H$. While $\xi$ in the scaling regime is the consequence of two
competing factors, namely decays of strings and more strings entering into the
Hubble patch, the quantity $\xi H$ will count only the effect from string
decays. Each benchmark scenario for $\xi_{0}=0.2$ and $\zeta=2.12$ consists of
100 independent simulations, and individual elements are spread around the
central values in $\xi$ and $q$ distributions. Over all samples, the
normalized differential change of $\xi H$ in logarithmic time, $d\log(\xi
H)/(d\log\frac{m_{r}}{H})$, and spectral index $q$ are computed at some time
slices for various choices of $\Delta\log$. An illustrative example taken at
$\log\sim 7.5$ is presented in Fig. 17 for the choice of $\Delta\log=0.5$, and
it reveals a strong correlation. The correlation is slightly less pronounced,
when $\Delta\log=0.25$ is used, and more pronounced with $\Delta\log=0.75$ in
that scattered data points are slightly less and more densely aligned along
the correlation axes for the cases with $\Delta\log=0.25$ and 0.75,
respectively. While the correlation is not stable in early times, it becomes
robust at late times, for instance, the orientation of the correlation axis
stays same for $\log\gtrsim 7$. A higher decay rate of strings is correlated
with a smaller $q$ according to Fig. 17. This looks consistent with the
analytic study in [32, 26] that the faster string decay rates are associated
with the emission of harder axions and vice versa. The detailed discussion is
postponed to Appendix D.2.

Figure 18: The $\xi H$ evolution from simulations using the fat-string type
pre-evolution (left) and thermal pre-evolution (right) on the lattice with
$N^{3}=1024^{3}$ for various initial conditions.

Assuming that the similar correlation holds for different averaged $\xi$
curves from different initial conditions of $\xi_{0}$ or $\zeta$, the
correlation in Fig. 17 can imply that the spectral index will have similar
values in the vicinity of the benchmark points. For instance, the
distributions of $\xi H$ from simulations on the smaller lattice with
$N^{3}=1024^{3}$ are illustrated in Fig. 18. Although the maximal dynamic time
range in Fig. 18 is smaller than $\log\sim 7$, it looks strong enough to
support that the values of $d\log(\xi H)/(d\log\frac{m_{r}}{H})$ will become
even closer at late times in the vicinity of the benchmark curves and thus the
spectral index will not be significantly different from those from our
benchmark simulations.

##  10 Conclusion

We have studied the dynamics of the global cosmic strings, originated from the
PQ symmetry breaking, in the post-inflationary scenario where the definite
prediction on the axion dark matter abundance can be made. We performed
independent numerical checks on recently discovered various logarithmic
scalings in properties of the cosmic string network. Especially, we provided
the strong evidences for two logarithmically growing scalings in the number of
strings in the Hubble patch and the spectral index of the power law of the
axion spectrum against two differently prepared initial field configurations.
Having established the strong evidence of $q>1$ around the QCD scale,
extrapolated from the simulation time, we pointed out a possible strong
correlation between the axion spectrum and the string evolution with different
initial conditions as a way of supporting the robust behavior of $q$ against
nearby string evolution trajectories with different initial data. We showed
that this strong correlation agrees with the theoretical prediction from a
simple sting vibration model. The IR dominated spectrum of axions radiated
from cosmic strings in the scaling regime may indicate the insensitivity to UV
physics of cosmic string modeling. It will be interesting to see if this
property can lead us an extra handle for better understanding the nature of
cosmic strings.

We surveyed various parametric behaviors of the axion abundance including the
impact from the value of the spectral index, nonlinearities around the QCD
scale, and average inter-string distances. This exercise should greatly help
in extracting the lower bound on the axion mass accurately from the simulation
result and comparing with different scenarios. Finally we estimated the lower
bound on the axion mass based on our simulation.

## Acknowledgments

We thank Masahide Yamaguchi for useful discussions. We thank Marco Gorghetto
for discussions about their previous work. We thank the department of physics
at KAIST for supporting with the computing resources in the cluster center.
HK, JP, and MS were supported by National Research Foundation of Korea under
Grant Number NRF-2021R1A2C1095430. The simulation in this work was supported
by KISTI National Supercomputing Center under Project Number
KSC-2023-CRE-0078.

##  Appendix A Equation of motion

The equation of motion from the Lagrangian in Eq. (1) in the expanding
universe with the scale factor of $R$ is given by

$\ddot{\phi}+3H\dot{\phi}-\frac{1}{R^{2}}\nabla^{2}\phi+\frac{m_{r}^{2}}{f_{a}^{2}}\phi\left(|\phi|^{2}-\frac{f_{a}^{2}}{2}\right)=0\leavevmode\nobreak\
,$ (22)

where dot is the differentiation with respect to the cosmic time $t$ and the
gradient $\nabla$ is evaluated in the comoving coordinate. To simply our
implementation in coding and numerical operation, we rescale field and
spacetime coordinates to make them dimensionless in such a way that
$m_{r}^{2}/f_{a}^{2}$ appears:

$\phi\rightarrow f_{a}\,\phi\leavevmode\nobreak\ ,\quad t\rightarrow
m_{r}^{-1}t\leavevmode\nobreak\ ,\quad\vec{x}\rightarrow
m_{r}^{-1}\vec{x}\leavevmode\nobreak\ ,$ (23)

under which $H\rightarrow m_{r}H$. After the rescaling, the equation of motion
becomes

$\ddot{\phi}+3H\dot{\phi}-\frac{1}{R^{2}}\nabla^{2}\phi+\phi\left(|\phi|^{2}-\frac{1}{2}\right)=0\leavevmode\nobreak\
,$ (24)

where now dot is the differentiation with respect to the dimensionless cosmic
time $t$ and the gradient $\nabla$ is evaluated in the dimensionless comoving
coordinate. Similarly $H$ is dimensionless.

Upon the discretization, the field evolution in terms of $\phi$ is subject to
the CFL condition, $\Delta t\lesssim R\Delta x$. This condition requires
arbitrarily small time step size in the early stage of the evolution where the
scale factor is small, and it is rather inefficient from the simulation point
of view. Since we want more sampling points at later times, we re-express the
equation of motion in terms of the rescaled field
$\psi(\tau,\,\vec{x})=R\phi(t,\,\vec{x})$ in the comoving spacetime
coordinate,

$\psi^{\prime\prime}-\nabla^{2}\psi+\psi\left(|\psi|^{2}-\frac{R^{2}}{2}-\frac{R^{\prime\prime}}{R}\right)=0\leavevmode\nobreak\
,$ (25)

where the $H$-dependent term disappeared and it makes the simulation simpler.
The CFL condition in this situation is relaxed to $\Delta\tau\lesssim\Delta
x$. While the last term $R^{\prime\prime}/R$ vanishes in the radiation
dominated era, it may not vanish in the universe with a generic scale factor.

In the approach using the relaxation during the thermal evolution, we add the
thermal mass term to the original Lagrangian,

$\Delta\mathcal{L}=-\frac{m_{r}^{2}}{6f_{a}^{2}}T^{2}|\phi|^{2}\leavevmode\nobreak\
.$ (26)

After the rescaling in Eq. (23) along with $T\rightarrow f_{a}T$ and the
replacement $\psi(\tau,\,\vec{x})=R\phi(t,\,\vec{x})$, the equation of motion
for $\phi$ in the radiation dominated era becomes

$\psi^{\prime\prime}-\nabla^{2}\psi+\psi\left(|\psi|^{2}-\frac{R^{2}}{2}+\frac{1}{6}R^{2}T^{2}\right)=0\leavevmode\nobreak\
,$ (27)

where $T^{2}=2\zeta H$ (the definition of $\zeta$ will be given in Appendix
B.5).

##  Appendix B Evolutions on the lattice

The approach toward the scaling solution can be boosted by the relaxation
which cleans up the noisy short-distance string structure. The pre-evolution
runs over the time interval $\tau=[\tau_{pi},\,\tau_{pf}]$. The outcome of the
pre-evolution is taken as the initial condition at $\tau=\tau_{i}$ for the
actual physical string evolution. In our numerical simulation, three different
time values, $\tau_{pi}$, $\tau_{pf}$, and $\tau_{i}$ are chosen to satisfy
certain properties. In this section, we describe the detail of our simulation
setup regarding two different types of the pre-evolutions.

In what follows, we present all the expressions in terms of the dimensionful
spacetime to avoid any confusion. To make it easier to convert the
dimensionful spacetime to the dimensionless one in the $m_{r}$ unit for the
convenience in the simulation, we adopt $R(t)=\sqrt{t/t_{0}}$ with
$t_{0}=1/m_{r}$ for the radiation dominated era. In terms of the comoving time
$\tau$, $R(\tau)=\frac{\tau}{2}m_{r}$ as $t=(\frac{\tau}{2})^{2}m_{r}$. For
the situation with the scaling of $R(t)\propto t$, we adopt $R(t)=t/t_{0}$
again with $t_{0}=1/m_{r}$. In terms of $\tau$, $R(\tau)=e^{m_{r}\tau}$
through $t=\frac{1}{m_{r}}e^{m_{r}\tau}$. One can pair $t$ and/or $x$ with
$m_{r}$ to convert it to the dimensionless one in the $m_{r}$ unit.

###  B.1 Common setup

We first explain how we set $\Delta x$ which allows the longest simulation
time in $m_{r}/H$. For given $N$, one can require the minimum number of the
lattice points within the core at the final time, denoted by $n_{c}$, and the
minimum number of the Hubble length within the simulation box $L$ at the final
time, denoted by $n_{H}$:

$\frac{m_{r}^{-1}}{R\Delta x}\geq n_{c}\leavevmode\nobreak\
,\quad\frac{L}{H^{-1}}=\frac{NR\Delta x}{H^{-1}}\geq n_{H}\leavevmode\nobreak\
\quad\rightarrow\quad\frac{m_{r}}{H}\leq\frac{N}{n_{c}n_{H}}\leavevmode\nobreak\
,$ (28)

where the maximum time coverage in $m_{r}/H$ is reached when two inequalities
in Eq. (28) for $n_{c}$ and $n_{H}$ are saturated at the same time. The
lattice spacing $\Delta x$ is set to guarantee the maximal $m_{r}/H$, and it
is given by

$\Delta
x=\left.\frac{m_{r}^{-1}}{n_{c}R(t)}\right|_{t=(1/2H)_{\text{max}}}=\frac{m_{r}^{-1}}{n_{c}}\left(\frac{2n_{c}n_{H}}{N}\right)^{1/2}\leavevmode\nobreak\
.$ (29)

The step size of the conformal time $\Delta\tau$ is chosen to
$\Delta\tau=r_{\text{CFL}}\,\Delta x$ where $r_{\text{CFL}}=1/3\sim 1/10$
depending on the simulation approach.

It is useful to compare this setup with that in the simulation assisted by the
AMR. When turning on the AMR which zooms in a local region recursively by the
factor of $2^{k}$, one can repeat the similar steps as in Eq. (28) with the
replacement of either $n_{c}\rightarrow n_{c}/2^{k}$ or $N\rightarrow 2^{k}N$
(up to the implication to $\Delta x$) to estimate the maximal time coverage.
In either choice, the maximum dynamic time range in $\log\frac{m_{r}}{H}$ is
extended to (by $k\ln 2$)

$\log\frac{m_{r}}{H}\leq\log\frac{N\cdot
2^{k}}{n_{c}n_{H}}\leavevmode\nobreak\ .$ (30)

Since $N\rightarrow 2^{k}N$ implies that the local resolution effectively
increases by the factor of $2^{k}$, $\Delta x$ in Eq. (28) corresponds to the
lattice spacing of the finest grids in the comoving coordinate,

$\left.\Delta x\right|_{\text{Finest girds by
}2^{k}}=\frac{m_{r}^{-1}}{n_{c}}\left(\frac{2n_{c}n_{H}}{2^{k}N}\right)^{1/2}\leavevmode\nobreak\
.$ (31)

On contrary, $n_{c}\rightarrow n_{c}/2^{k}$ in Eq. (28) implies that the
minimum number of grids within the core is relaxed by the factor of $2^{k}$
since the local region around the string core will be resolved later up to by
$2^{k}$. In this case, $\Delta x$ in Eq. (28) corresponds to the lattice
spacing of the coarsest grids in the comoving coordinate,

$\left.\Delta x\right|_{\text{Coarest
girds}}=\frac{m_{r}^{-1}}{n_{c}}\left(\frac{2n_{c}n_{H}}{N}\right)^{1/2}\times
2^{k/2}\leavevmode\nobreak\ ,$ (32)

which is $2^{k}\times\left.\Delta x\right|_{\text{Finest girds by }2^{k}}$ and
it should be compared to the lattice spacing of the static simulation for the
region away from strings. A peculiar feature of the scaling in Eq. (32) is
that the lattice spacing away from strings in AMR is larger than that of the
static case by the factor of $2^{k/2}$ for fixed $n_{c}$ and $N$. This can be
compensated with a larger $n_{c}$.

###  B.2 Fat string pre-evolution: $R(t)\propto\sqrt{t}$

The first scheme introduced in [9] is that, during the fat-string pre-
evolution, the core radius $m_{r,\text{fat}}^{-1}$ scales as
$R\propto\sqrt{t}$ while the VEV of $\phi$ is fixed to be $f_{a}$, or
equivalently, the string is made to be fat such that its size $r_{\text{fat}}$
in the comoving coordinate stays the same during the evolution. The quartic
coupling scales as $1/R^{2}$ 101010The decreasing quartic coupling in this
scheme indicates that the field becomes less and less interacting with the
evolution, and it might be the reason for why the simulation reach the
attractor solution faster than the case for the physical string simulation.,

$\frac{m_{r,\text{fat}}^{2}}{2f^{2}_{a}}=\frac{1}{2f_{a}^{2}R^{2}\left(m_{r,\text{fat}}^{-1}/R\right)^{2}}=\frac{1}{2f_{a}^{2}R^{2}r_{\text{fat}}^{2}}\propto\frac{1}{R^{2}}\leavevmode\nobreak\
.$ (33)

We choose the string core size during the fat string pre-evolution to match
the value at the initial time $\tau_{i}$ of the physical evolution, or

$r_{\text{fat}}(\tau_{i})=\left.\frac{m_{r,\text{fat}}^{-1}}{R}\right|_{\tau=\tau_{pi}}=\frac{m_{r}^{-1}}{R(\tau_{i})}=m_{r}^{-1}\times\frac{2}{m_{r}\tau_{i}}\leavevmode\nobreak\
.$ (34)

The number of lattice points within the string core at the initial time
$\tau_{i}$, denoted by $N_{r}$, of our physical string simulation is given by

$\begin{split}N_{r}=\left.\frac{m_{r,\text{fat}}^{-1}}{\Delta
xR}\right|_{\tau=\tau_{pi}}=\frac{m_{r}^{-1}}{\Delta
x}\times\frac{2}{m_{r}\tau_{i}}\leavevmode\nobreak\ .\end{split}$ (35)

This allows us to start the physical evolution from the same initial time
$\tau_{i}$ for the given $\Delta x$ and $N_{r}$. We choose the initial time of
the pre-evolution $\tau_{pi}$ such that the number of lattice within the
Hubble length at $\tau_{pi}$ is roughly a similar order of the one within a
wavelength of the string core, $2\pi m_{r}^{-1}$, at the time $\tau_{i}$ up to
an arbitrary constant $\alpha$,

$\begin{split}\underbrace{\frac{\tau_{pi}}{\Delta
x}=\left.\frac{H^{-1}}{\Delta
xR}\right|_{\tau_{pi}}}_{{\begin{subarray}{c}\text{Number of lattice
within}\\\ \text{the Hubble
length}\end{subarray}}}=\alpha\times\underbrace{\left.\frac{m_{r}^{-1}}{\Delta
xR}\right|_{\tau_{i}}}_{\begin{subarray}{c}\text{Number of lattice}\\\
\text{within the string core}\end{subarray}}\leavevmode\nobreak\ ,\end{split}$
(36)

where the second term is equal to $\alpha N_{r}$ which stays the same over the
entire fat string evolution. The role of the condition in Eq. (36) is to
achieve roughly the reasonable size of the correlation length of the cosmic
string at $\tau_{i}$. The ambiguity in the identification in Eq. (36) is taken
into account by the constant $\alpha$ that can be empirically determined. The
setting in Eq. (36) ensures that
$\frac{m_{r,\text{fat}}}{H}|_{\tau_{pi}}\sim\alpha$ which is taken to be order
one. The ending time of the fat-string evolution $\tau_{pf}$ is determined by
demanding the requested value for the number of strings per Hubble patch
$\xi_{0}$, namely, $\xi(\tau_{pf})=\xi_{0}$111111While the fat string
evolution can be designed to be terminated at $\tau_{pf}$ when it reaches the
same requested correlation length, the initial value $\xi_{0}$ will not be
same in this case.. The value of $\tau_{pf}$ varies depending on $\xi_{0}$,
and the correlation length at $\tau_{pf}$ (which will be the property of the
initial condition for the physical string evolution) is determined posteriori.
A drawback is that the correlation length with respect to the core radius can
significantly grow until $\tau_{pf}$ due to the scaling,
$\frac{m_{r,\text{fat}}}{H}|_{\tau_{pf}}\propto\tau$. While the value of
$\tau_{i}$ for the physical string evolution can be freely chosen irrespective
of $\tau_{pf}$ ($\tau_{i}$ is not necessarily larger $\tau_{pf}$), the
implicit correlation length at $\tau_{i}$ is genuinely determined by the
configuration of fields at $\tau_{pf}$ instead of being determined in terms of
$H^{-1}(\tau_{i})$. While the condition in Eq. (36) looks a bit ad-hoc, it
will be replaced by more consistent one in the second scheme.

###  B.3 Fat string pre-evolution: $R(t)\propto t$

The second scheme that we take as a default option for the pre-evolution was
introduced in [14]. In this scheme, during the fat-string evolution, the core
radius $m_{r}^{-1}$ scales as $R_{\text{fat}}\propto t\sim
H_{\text{fat}}^{-1}$ while the VEV of $\phi$ is fixed to be $f_{a}$. Similarly
to the option 1 in Section B.2, we match the core radii in the comoving
coordinate in the transition from the pre-evolution to the physical one,

$r_{\text{fat}}(\tau_{i})=\left.\frac{m_{r,\text{fat}}^{-1}}{R_{\text{fat}}}\right|_{\tau=\tau_{pi}}=\frac{m_{r}^{-1}}{R(\tau_{i})}=m_{r}^{-1}\times\frac{2}{m_{r}\tau_{i}}\leavevmode\nobreak\
.$ (37)

The number of lattice points within the string core at the initial time
$\tau_{i}$ will be the same as that during the pre-evolution. Since
$H_{\text{fat}}^{-1}/R_{\text{fat}}$ stays constant during the pre-evolution,
we can maintain the exactly same correlation length in the comoving coordinate
at $\tau_{i}$ as the one at $\tau_{pi}$ by demanding

$\begin{split}\left.\frac{H_{\text{fat}}^{-1}}{\Delta
xR_{\text{fat}}}\right|_{\tau_{pi}}=\left.\frac{H^{-1}}{\Delta
xR}\right|_{\tau_{i}}=\frac{\tau_{i}}{\Delta x}\leavevmode\nobreak\
,\end{split}$ (38)

which sets the comoving Hubble length $H^{-1}/R$ during the pre-evolution to
$\tau_{i}$. Combining two relations in Eqs. (37) and (38), it is
straightforward to see that the correlation length over the core size stays
the same, that is,

$\frac{m_{r,\text{fat}}(\tau)}{H_{\text{fat}}(\tau)}=\frac{R_{\text{fat}}(\tau)}{m^{-1}_{r,\text{fat}}(\tau)}\frac{H^{-1}_{\text{fat}}(\tau)}{R_{\text{fat}}(\tau)}=\left.\frac{R}{m_{r}^{-1}}\right|_{\tau_{i}}\left.\frac{H^{-1}}{R}\right|_{\tau_{i}}=\frac{m_{r}}{H(\tau_{i})}\leavevmode\nobreak\
,$ (39)

where the initial correlation length over the core size is given by
$\frac{m_{r}}{H(\tau_{i})}=\frac{1}{2}m_{r}^{2}\tau_{i}^{2}$. The initial time
for the physical string evolution $\tau_{i}$ is fixed by Eq. (38) 121212The
freedom to choose any value of $\tau_{i}$ (or the freedom to choose any
desired value of $\frac{m_{r}}{H(\tau_{i})}$) effectively amounts to
introducing $\alpha$ in $R(t)=\alpha(t/t_{0})$, or $\tau_{i}=\frac{1}{\alpha
m_{r}}$ from Eq. (38). in this scheme which accordingly determines $N_{r}$
through Eq. (37). Similarly to the option 1 in Section B.2, the final time of
the pre-evolution is determined by requesting a certain number of strings per
Hubble patch $\xi_{0}$, namely, $\xi(\tau_{pf})=\xi_{0}$. While the initial
time $\tau_{pi}$ of the pre-evolution is apparently a free parameter, we set
it to $\tau_{pi}=0$.

###  B.4 From pre-evolution to physical one

Since field equations in two regions in options 2 and 3 are qualitatively
different and not physically equivalent, there may not be a unique way of
gluing two regions. In this work, we join two regions by two relations,

$\begin{split}\left.(R^{-1}\psi)\right|_{\tau=\tau_{pf}}&=\left.(R^{-1}\psi)\right|_{\tau=\tau_{i}}\leavevmode\nobreak\
,\\\\[5.0pt]
\left.\left(-\frac{R^{\prime}}{R^{2}}\psi+\frac{1}{R}\psi^{\prime}\right)\right|_{\tau=\tau_{pf}}&=\left.\left(-\frac{R^{\prime}}{R^{2}}\psi+\frac{1}{R}\psi^{\prime}\right)\right|_{\tau=\tau_{i}}\leavevmode\nobreak\
,\end{split}$ (40)

where the second relation is the continuity of $(R^{-1}\psi)^{\prime}$, or the
velocity of $R^{-1}\psi$ with respect to the conformal time and it ensures
that the energy density is smoothly connected between two regions 131313While
the velocity $\dot{\phi}$ is discontinuous in this scheme, we have numerically
checked that the impact of the discontinuity on the scaling solution is
negligible.. We empirically find that maintaining the continuity of the energy
density in the transition is strongly favored not to artificially introduce
the discrete enhancement of the radial modes in the transition from the pre-
evolution to the physical one.

###  B.5 Thermal pre-evolution

We continue presenting the relations in terms of the dimensionful parameters.
The expressions in terms of the dimensionless quantities can be
straightforwardly obtained as the previous discussion along with the rescaling
$T$ by $f_{a}$. The critical temperature $T_{c}$ is given by
$T_{c}=\sqrt{3}f_{a}$. Using the relation for the Hubble parameter
$H^{2}=\frac{\pi^{2}}{90}g_{*}\frac{T^{4}}{M_{p}^{2}}$, the temperature can be
expressed as

$T^{4}=f_{a}^{4}\times\frac{H^{2}}{m_{r}^{2}}\times
4\,\zeta^{2}\leavevmode\nobreak\
\quad\text{with}\quad\zeta^{2}=\frac{45}{2\pi^{2}g_{*}}\frac{M_{p}^{2}m_{r}^{2}}{f_{a}^{4}}\leavevmode\nobreak\
,$ (41)

from which the correlation length over the core size at the critical
temperature $T_{c}$ is given by

$\left.\frac{m_{r}}{H}\right|_{T_{c}}=\frac{2}{3}\zeta\leavevmode\nobreak\ .$
(42)

While $\zeta\sim\mathcal{O}(1)$ ensures the order one size of
$\frac{m_{r}}{H}$ at the phase transition where cosmic strings are formed, the
natural size of it is given by

$\zeta=2\times
10^{7}\left(\frac{100}{g_{*}}\right)^{1/2}\left(\frac{\lambda}{1}\right)^{1/2}\left(\frac{M_{p}}{1.\times
10^{18}\text{GeV}}\right)\left(\frac{1.\times
10^{10}\text{GeV}}{f_{a}}\right)\leavevmode\nobreak\ .$ (43)

Eq. (43) implies that the order one size of $\zeta$ corresponds to the
symmetry breaking scale of $f_{a}\sim 10^{17}$ GeV, a huge number of
relativistic degrees of freedom of $g_{*}\sim 10^{14}$, or an extremely small
quartic coupling 141414This looks like a similar unnatural situation in the
fat-string type pre-evolution which relies on the quartic coupling scaling
away as $R^{-2}$.  $\lambda=\frac{m_{r}^{2}}{2f_{a}^{2}}\sim 10^{-15}$ (or the
combination of three factors). This might be considered to be a drawback of
this option. The initial time of the simulation $\tau_{i}$ can be chosen to be
any value earlier than the phase transition time $\tau_{c}$ which is given by,
using Eq. (42),

$\tau_{c}^{2}=m_{r}^{-2}\times\frac{4}{3}\zeta\leavevmode\nobreak\ .$ (44)

In our numerical simulation, we fix the initial time to $\tau_{i}=0.1\tau_{c}$
with the varying $\zeta$. The value of $\zeta$ can be fixed to have a specific
Hubble length in the $m_{r}^{-1}$ unit at the critical temperature as was
indicated in Eq. (42), for instance, $\zeta=\frac{3}{\sqrt{2}}\sim 2.12132$
for $\left.\frac{m_{r}}{H}\right|_{\tau_{c}}=\sqrt{2}$, and it implies
$\left.\frac{m_{r}}{H}\right|_{\tau_{i}}=(\frac{\tau_{i}}{\tau_{c}})^{2}\left.\frac{m_{r}}{H}\right|_{\tau_{c}}$
at the initial time $\tau_{i}$.

###  B.6 Random initial conditions

In fat string type relaxations explained in Sections B.2 and B.3, a random
field configuration $\phi$ (and its velocity $\dot{\phi}$) is generated with
the Fourier modes with the wavenumber truncated at $m_{r}$, or
$|\vec{k}|_{\text{max}}=\frac{2\pi}{N\Delta
xR(\tau_{i})}|\vec{n}_{k}|_{\text{max}}=m_{r}$. The field can be discretized
as $\phi(\vec{x}=\Delta x\,\vec{n}_{x})=\sum_{\vec{n}_{k}}\frac{1}{(N\Delta
xR)^{3}}\tilde{\phi}\left(\frac{2\pi}{N\Delta
xR}\vec{n}_{k}\right)\exp{\frac{2\pi i\vec{n}_{k}\cdot\vec{n}_{x}}{N}}$ where
the Fourier mode $\tilde{\phi}(\vec{k})$ is randomly generated below the
cutoff scale of the momentum as
$\tilde{\phi}_{0}\,\frac{\tt{Rand}_{1}+i\tt{Rand}_{2}}{\sqrt{2}}$ where
$\tt{Rand}_{1,2}$ are random numbers following the standard normal
distribution such that $\langle{\tt Rand}_{i}(\vec{n}_{k_{1}})\,{\tt
Rand}_{j}(\vec{n}_{k_{2}})\rangle=\delta_{ij}\delta_{\vec{n}_{k_{1}}\vec{n}_{k_{2}}}(N\Delta
xR)^{3}$. The overall size $\tilde{\phi}_{0}$ is chosen such that
$\langle|\phi(\vec{n}_{x})|^{2}\rangle=f_{a}^{2}\times\frac{1}{2}$:

$\langle\phi^{*}(\vec{n}_{x_{1}})\phi(\vec{n}_{x_{2}})\rangle=|\tilde{\phi}_{0}|^{2}\sum_{\vec{n}_{k}}\frac{1}{(N\Delta
xR)^{3}}\exp^{-\frac{2\pi
i\vec{n}_{k}\cdot(\vec{n}_{x_{1}}-\vec{n}_{x_{2}})}{N}}\leavevmode\nobreak\ ,$
(45)

where the sum is truncated at $|\vec{n}_{k}|_{\text{max}}$. The randomly
generated Fourier mode $\tilde{\phi}$ is given by

$\tilde{\phi}=\frac{f_{a}}{\sqrt{2}}\frac{(N\Delta
xR)^{3/2}}{\sqrt{N_{k}}}\frac{{\tt{Rand}_{1}}(\vec{n}_{k})+i{\tt{Rand}_{2}}(\vec{n}_{k})}{\sqrt{2}}\theta(|\vec{n}_{k}|_{\text{max}}-|\vec{n}_{k}|)\leavevmode\nobreak\
,$ (46)

where $N_{k}$ is the total number of lattice points within the sphere in the
momentum space, which is defined as $N_{k}\equiv\frac{1}{(\Delta k)^{3}}\int
d^{3}\vec{k}\,\theta(|\vec{n}_{k}|_{\text{max}}-|\vec{n}_{k}|)$ in terms of
the notation in the continuum limit, and $\theta$ is a step function. The
random velocity distribution for $\phi^{\prime}$ is similarly generated with
the Fourier mode further weighted by $RE$, compared to $\tilde{\phi}$ in Eq.
(46), where $E(\vec{k}=2\pi\vec{n}_{k}/(N\Delta
xR))=\sqrt{m_{r}^{2}+\vec{k}^{2}}$.

In the thermal relaxation scheme in Section B.5, the initial conditions for
the field and velocity are generated assuming the Gaussian random field
configuration following the thermal distributions as in [10]. Similarly to the
fat string type relaxation, a random field configuration and its velocity will
be generated with an appropriate cut on the wavenumber. For the periodic
boundary condition, the lowest wavenumber in the $m_{r}$ unit at the initial
time is given by
$\left.\frac{k_{\text{min}}}{m_{r}}\right|_{t_{i}}=\left.\frac{2\pi}{m_{r}RN\Delta
x}\right|_{t_{i}}=\frac{10\sqrt{3}}{\sqrt{\zeta}}\frac{2\pi}{\left(2n_{H}N/n_{c}\right)^{1/2}}$.
The cutoff on the physical momentum will be imposed as the multiple of
$m_{r}$. It will correspond to including only $|\vec{n}_{k}|_{\text{max}}$
number of lowest modes, or $k_{\text{max}}=\frac{2\pi}{RN\Delta
x}|\vec{n}_{k}|_{\text{max}}=\kappa\,m_{r}$ at $\tau=\tau_{i}$ from which
$|\vec{n}_{k}|_{\text{max}}=\kappa\left.(k_{\text{min}}/m_{r})^{-1}\right|_{\tau_{i}}$.
For $\zeta=2.12$ and $N=1024/2048/4096$ along with our default setting
$n_{H}=4^{1/3}$ and $n_{c}=1$,
$\left.(k_{\text{min}}/m_{r})^{-1}\right|_{\tau_{i}}\sim 0.763/1.079/1.526$.
Therefore, the cut $k_{\text{max}}|_{\tau_{i}}\leq 10\,m_{r}$ is equivalent to
including only $|\vec{n}_{k}|\lesssim$ 8/11/15 lowest modes, respectively. As
$k_{\text{max}}\propto\tau^{-1}$, $k_{\text{max}}|_{\tau_{i}}\leq 10\,m_{r}$
will imply $k_{\text{max}}|_{\tau_{c}}\leq m_{r}$ with the choice of
$\tau_{i}=0.1\tau_{c}$ (note that $m_{r}$ denotes the constant string core
scale).

###  B.7 Tetrahedron-based string identification

We subdivide the entire lattice box into a series of tetrahedrons as in the
left panel of Fig. 1, and decide whether a string penetrates a face or not by
examining the pattern of the phases in three vertices. Once the criterion for
the existence of a string is met, we simply determine the string core by the
linear interpolation (although it is not guaranteed to be the truth position).
In this approach, we estimate total string length in the simulation box by
summing over all line segments connecting string cores. On the other hand,
since the position of the string core over the plaquette (see right panel of
Fig. 1) through which a string penetrates is likely randomly distributed from
the statistics point of view, the total string length in the simulation box
can be estimated by just counting the number of plaquettes and multiply it by
a weighting factor [34]. Using the simulation on the grid with
$N^{3}=1024^{3}$, we have compared our algorithm with the string
identification in [34], and the discrepancy $\Delta\xi/\xi$ over the entire
simulation time was found to be smaller than a few %. Since we want to
visualize the string formation and take snapshots at some selected times, we
take our approach for the estimation of the total string length. The number of
long strings per Hubble patch in our simulation is calculated as

$\xi=\frac{(L_{\text{Lattice unit}}\Delta xR)t^{2}}{(N\Delta
xR)^{3}}=\frac{L_{\text{Lattice unit}}t^{2}}{N^{3}\Delta
x^{2}R^{2}}\leavevmode\nobreak\ .$ (47)

where $L_{\text{Lattice unit}}$ is the total length of strings measured in a
lattice unit.

In our prescription using the tetrahedron, the surface of the tetrahedron in
the left panel of Fig. 1 is the triangle with three vertices, call them
$v_{i=1,2,3}$. Given the complex field values on three vertices, one can
define $\theta_{ij}=\arg(\phi_{j}/\phi_{i})=\arg(\phi_{j}\phi_{i}^{*})$. A
string core on the triangular surface is declared if
$\theta_{123}\equiv\theta_{12}+\theta_{23}+\theta_{31}$ has the value of $\pm
2\pi$. Note that $\theta_{123}$ can have only $0,\,\pm 2\pi$. In practice, the
computation of the arguments for all complex fields are time-consuming. To
reduce the computational time, we replace the above criterion with more
efficient equivalent step. $\theta_{123}=2\pi$ ($-2\pi$) is equivalent to the
situation where all $\theta_{ij}$’s are non-negative (non-positive) while not
all of them are zero. The sign of $\theta_{ij}$ is equivalently replaced with
that of $\omega_{ij}=\Re\phi_{i}\Im\phi_{j}-\Im\phi_{i}\Re\phi_{j}$.
Therefore, $\theta_{123}=2\pi$ ($-2\pi$) is equivalent to the case where all
$\omega_{ij}$’s are non-negative (non-positive) while $\omega_{123}$ is
positive (negative) ($\omega_{123}\equiv\omega_{12}+\omega_{23}+\omega_{31}$).
Once the string core is identified, the location of the string core inside the
triangle where the field value vanishes is assigned via the linear
interpolation, or
$(\omega_{23}v_{1}+\omega_{31}v_{2}+\omega_{12}v_{3})/\omega_{123}$ where
$v_{i}$ is the position of the $i$-th vertex on the triangle. It can be easily
proven that the string core is strictly inside the triangle. This approach is
well-defined as each tetrahedron includes either a pair of two faces on which
string cores are found or none, and thus strings are guaranteed to be
connected without any discontinuity. We find that our tetrahedron-based string
identification algorithm is quite efficient in the numerical simulation.

The connectedness of the string in our tetrahedron approach of the string
identification can be rigorously proven as follows. One can introduce an
arbitrary complex field as a linear function of coordinates,
$\phi_{\text{inter}}(x,\,y,\,z)=\beta_{0}+\beta_{1}x+\beta_{2}y+\beta_{3}z$
with four complex constants which we fix them by matching the interpolating
field $\phi_{\text{inter}}$ to field values at four vertices, labeled as
$v_{i=0,1,2,3}$, of the tetrahedron, or $\phi_{\text{inter}}(v_{i})=\phi_{i}$.
That is, $\phi_{\text{inter}}(x,\,y,\,z)$ interpolates between vertices of the
tetrahedron in three-dimensional space. The linear interpolation of the field
on the triangular face of the tetrahedron is equivalent to restricting
$\phi_{\text{inter}}$ onto the triangle. In this picture, the corresponding
triangular face will be declared to have a string core if and only if the
triangular face and the locus of zero, $\phi_{\text{inter}}(x,\,y,\,z)=0$,
intersect. $\phi_{\text{inter}}(x,\,y,\,z)=0$ represents the straight line as
it is collection of two constraints (for real and complex components of
$\beta$’s) in three-dimensional space. A straight line can penetrate only a
pair of triangular faces of the tetrahedron if it intersects unless the line
passes through a vertex or an edge of the tetrahedron. Therefore, there must
be none or exactly two triangular faces which are penetrated by a cosmic
string. This prescription will work well for long strings compared to the size
of the tetrahedron.

###  B.8 Masking

The axion energy density is approximated by
$2\times\langle\frac{1}{2}\dot{a}^{2}\rangle$ where $\langle\cdots\rangle$
denotes (with the abuse of the notation) the spatial average over the physical
space (denoted by $\vec{x}_{p}=R\,\vec{x}$),

$\rho_{a}=\frac{1}{L^{3}}\int
d^{3}\vec{x}_{p}\,\dot{a}^{2}(\vec{x}_{p})=\frac{1}{L^{3}}\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\,|\tilde{\dot{a}}(\vec{k})|^{2}\leavevmode\nobreak\
,$ (48)

from which the axion energy spectrum for the momentum $k$ is given by

$\frac{\partial\rho_{a}}{\partial
k}=\frac{k^{2}}{2\pi^{2}L^{3}}\int\frac{d\Omega_{k}(\hat{k})}{4\pi}\left|\tilde{\dot{a}}(\vec{k}=k\hat{k})\right|^{2}\leavevmode\nobreak\
.$ (49)

Importantly, $\dot{a}(\vec{x}_{p})$ in Eq. (48) should be the axion field
velocity with the string effect removed. To this end, in presence of strings,
the masking function needs to be applied to the unmasked axion velocity field
to screen off the effect from string cores and extract the pure axion energy
density,

$\tilde{\dot{a}}_{\text{masked}}(\vec{k})=\int
d^{3}\vec{x}_{p}\,W(\vec{x}_{p})\,\dot{a}_{\text{unmasked}}(\vec{x}_{p})e^{-i\vec{k}\cdot\vec{x}_{p}}\leavevmode\nobreak\
,$ (50)

where $W(\vec{x}_{p})$ is the masking function. We adopt the natural masking
$W(\vec{x}_{p})=(1+r(\vec{x}_{p})/f_{a})$ for the smooth screening the string
core effect as a default. Another choice for the masking adopted in literature
is the hard-masking which sets $W=0$ for the region with distance from string
$r<d_{s}m_{r}^{-1}$ (for some positive $d_{s}$) and $W=1$ outside the region.
The masked Fourier transform $\tilde{\dot{a}}_{\text{masked}}(\vec{k})$ is
related to the Fourier transform $\tilde{\dot{a}}(\vec{k})$ of the axion
velocity field $\dot{a}(\vec{x}_{p})$, now separated from the string effect,
by the convolution,

$\tilde{\dot{a}}_{\text{masked}}(\vec{k})=\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\,\tilde{W}(\vec{k}-\vec{k}^{\prime})\,\tilde{\dot{a}}(\vec{k^{\prime}})\leavevmode\nobreak\
,$ (51)

where $\tilde{W}(\vec{k})$ is the Fourier transform of the masking function
$W(\vec{x}_{p})$ itself. In our simulation,
$\mu_{a,\text{masked}}(k)=\int\frac{d\Omega_{k}}{4\pi}|\tilde{\dot{a}}_{\text{masked}}(\vec{k})|^{2}$
and similarly
$\mu_{W}(k)=\int\frac{d\Omega_{k}}{4\pi}|\tilde{W}(\vec{k})|^{2}$ are
numerically computed. The value of
$\mu_{a}(k)=\int\frac{d\Omega_{k}}{4\pi}|\tilde{\dot{a}}(\vec{k})|^{2}$ for
the axion velocity with the string core effect screened off can be obtained by
inverting the following relation,

$\mu_{a,\text{masked}}(k)=\frac{1}{4\pi^{2}kL^{3}}\int
dk^{\prime}k^{\prime}\left(\int^{k+k^{\prime}}_{|k-k^{\prime}|}dk^{\prime\prime}k^{\prime\prime}\mu_{W}(k^{\prime\prime})\right)\mu_{a}(k)\leavevmode\nobreak\
.$ (52)

This relation is derived under the assumption that each $\vec{k}$ modes of
$\dot{a}$ are independent and has random amplitude from mean zero Gaussian
ensemble whose variance only depends on the size of the momentum $k$. Here,
$\mu_{a,\text{masked}}$ and $\mu_{a}(k)$ need to be the average over the
ensemble to make the relation mathematically exact. Using the relation to get
$\mu_{a}$ from the measured $\mu_{a,\text{masked}}$ is secretly relying on the
law of large number applied to $\mu_{a,\text{masked}}(\vec{k})$ for
$\vec{k}$’s with large number of different directions and the same magnitude.
However, once space is discretized, and so is the momentum, there are only a
few number of different momentum $\vec{k}$ for small magnitude $k$, and this
may result in larger uncertainty at low momentum region while estimating the
unmasked spectrum. Using $\mu_{a}(k)$ extracted from Eq. (52), the axion
energy spectrum for the momentum $k$ with the string core effect screened off
is given by

$\frac{\partial\rho_{a}}{\partial
k}=\frac{k^{2}}{2\pi^{2}L^{3}}\mu_{a}(k)\leavevmode\nobreak\ .$ (53)

Similar masking is applied to the radial modes.

##  Appendix C Dependency on initial conditions and relaxations

###  C.1 Scaling regime and fitting

The underlying assumption behind the pre-evolution is that the string network
eventually approaches the scaling solution. Any residual influence from the
different pre-evolutions or/and different initial conditions should be faded
away as the string network enters the scaling regime. Any justification for it
within the consistent numerical setup should be highly beneficial. To this
end, we extend our discussion in Section 5 to include more comprehensive and
detailed comparisons between results obtained by two different pre-evolution
schemes for various initial data. We do this task mainly with simulations on
the lattice of $N^{3}=1024^{3}$ to accumulate as many simulation runs as
possible for various initial conditions within a reasonable time. As was
discussed in Section 5 and illustrated in Fig. 2, this task helped us in
choosing our benchmark points to be run on a larger lattice space for the
study at much later dynamic time.

The left panel of Fig. 19 is similar to the right panel of Fig. 2 except that
the cutoffs on the momentum were not imposed. That is, the initial
configurations with $k$ modes between maximally allowed value (set by the
lattice spacing and lattice size) and $10m_{r}$ at $\tau_{i}=0.1\tau_{c}$ were
included. While the $\xi$ curves in the left panel of Fig. 19 illustrate the
approach toward the scaling regime, the $\xi$ values are overally larger than
those with momentum cutoffs as is evident in the right panel of the same Fig.
19 where the dependence on the $k$ cutoffs is demonstrated for two selected
$\zeta$ values.

Figure 19: The impact of higher modes above $k_{\text{max}}=10m_{r}$ on $\xi$
for physical string simulations with initial data generated with the thermal
pre-evolution. Note that $k_{\text{max}}=10m_{r}$ at $\tau_{i}=0.1\tau_{c}$ is
equivalent to $k_{\text{max}}=m_{r}$ at $\tau_{c}$.

Since the maximally allowed $k$ modes varies with the lattice size, $\xi$
curves obtained from simulations with different lattice sizes do not overlap,
and the simulation on the larger lattice space populates more strings and
takes longer time to be relaxed. The cutoffs on the momentum removes this
ambiguity as is illustrated in Fig. 20 for two benchmark simulations:
$\xi_{0}=0.2$ using fat-string pre-evolution (left panel) and $\zeta=2.12$
(with $\tau_{i}=0.1\tau_{c}$) using thermal pre-evolution (right panel). Since
the same momentum modes are initially randomly generated irrespective of the
lattice size, the $\xi$ curves exactly on top of each other 151515The slight
mismatch between $N^{3}=1024^{3}$ and $N^{3}=2048^{3}$ in the right panel of
Fig. 20 is due to the round-off error, for instance, the cutoff $k>10m_{r}$
includes $k$ modes within the radius of roughly $\sim 8$ and $\sim 11$ for
$N^{3}=1024^{3}$ and $N^{3}=2048^{3}$, respectively, (see Appendix B.6) and
rounding-off causes an error up to 10%. We have checked via independent
simulations that increasing the cutoff, allowing more $k$ modes, makes the
effect due to the round-off error negligible. A similar statement applies to
the case with $N^{3}=4096^{3}$.. A good agreement between different lattice
sizes in Fig. 20 indicates that the impact of the lattice spacing is not
relevant as long as the string core includes at least one grid at the final
time, or $n_{c}=1$. At the final time of the simulation with $N^{3}=1024^{3}$
where $n_{c}^{1024}=1$, the number of grids within the string core for the
simulation with $N^{3}=2048^{3}$ is $n_{c}^{2048}|_{t_{f}^{1024}\sim
6.3}=\sqrt{2}$.

Figure 20: The independence of $\xi$ curves on the lattice spacing $\Delta x$
(and thus lattice size $N$ as well) for simulations with the momentum cutoffs
imposed.

To quantify the quality of the scaling solution, we fit the curves to the form
in Eq. (6). Given an ensemble of $N_{\text{ens}}$ independent runs for $\xi$
curves, we select the $N_{\text{fit}}$ number of sampling points in the
logarithmic time $\log\frac{m_{r}}{H}$ for fitting the curve. The Gaussian
likelihood function is constructed as

$L(c_{0},c_{1},d_{1})=\frac{1}{(2\pi)^{N_{\text{fit}}/2}\sqrt{\det(S)}}e^{-\frac{1}{2}\sum_{i,j=1}^{N_{\text{fit}}}\left(\bar{\xi}_{i}-\xi_{i}^{\text{sc}}(c_{0},\,c_{1},\,d_{1},\,d_{2})\right)\left(S^{-1}\right)_{ij}\left(\bar{\xi}_{j}-\xi_{j}^{\text{sc}}(c_{0},\,c_{1},\,d_{1},\,d_{2})\right)}\leavevmode\nobreak\
,$ (54)

where $\bar{\xi}_{i}$ is the mean of $\xi$ data at the time
$\log_{i}=\log\frac{m_{r}}{H}|_{i}$ and the ansatz for the scaling solution is
given by
$\xi_{i}^{\text{sc}}(c_{0},\,c_{1},\,d_{1},\,d_{2})=c_{0}+c_{1}\log_{i}+\frac{d_{1}}{\log_{i}}+\frac{d_{2}}{\log^{2}_{i}}$.
The matrix $S$ in Eq. (54) is the $N_{\text{fit}}\times N_{\text{fit}}$
covariance matrix of the sampling points along the curve whose element is
given by
$S_{ij}=\frac{1}{N_{\text{ens}}-1}\sum_{k=1}^{N_{\text{ens}}}\left(\xi_{i,k}-\bar{\xi}_{i}\right)\left(\xi_{j,k}-\bar{\xi}_{j}\right)$.
The uncertainties $\sigma_{c_{0}}$, $\sigma_{c_{1}}$ and $\sigma_{d_{1}}$ for
the $c_{0}$, $c_{1}$ and $d_{1}$ are estimated to be
$\frac{1}{\sigma_{X}^{2}}=-\frac{\partial^{2}}{\partial X^{2}}\log L$ for
$X=\\{c_{0},\,c_{1},\,d_{1}\\}$ at the maximizing point. One may ignore the
correlation among $\xi_{i}$ at different times for a large enough separation
$\Delta\log$ between two consecutive $\log_{i}$, or

$L(c_{0},c_{1},d_{1})=\prod_{i=1}^{N_{\text{fit}}}\frac{1}{\sqrt{2\pi}\sigma_{i}}e^{-\frac{1}{2\sigma_{i}^{2}}\left(\bar{\xi}_{i}-\xi_{i}^{\text{sc}}(c_{0},c_{1},d_{1},d_{2})\right)^{2}}\leavevmode\nobreak\
,$ (55)

where
$\sigma_{i}^{2}=\frac{1}{N_{\text{ens}}-1}\sum_{j=1}^{N_{\text{ens}}}\left(\xi_{i,j}-\bar{\xi}_{i}\right)^{2}$.
Our fitting result is presented in Table 3. We find that the fitting results
are sensitive to the fit ansatz and the interval.

Fat string pre-evolution | Thermal pre-evolution
---|---
$\xi_{0}$ | Lattice size ($N^{3}$) | Cuts | $N_{\text{runs}}$ | $\zeta$ | Lattice size ($N^{3}$) | Cuts at $\tau_{i}=0.1\tau_{c}$ | $N_{\text{runs}}$
0.1 | $1024^{3}$ | $k\leq m_{r}$ | 160 | 0.5 | $1024^{3}$ | $k\leq 10m_{r}$ | 130
0.2 | $1024^{3}$ | $k\leq m_{r}$ | 160 | 1.5 | $1024^{3}$ | $k\leq 10m_{r}$ | 160
0.3 | $1024^{3}$ | $k\leq m_{r}$ | 160 | 1.5 | $1024^{3}$ | no cuts | 130
0.4 | $1024^{3}$ | $k\leq m_{r}$ | 160 | 2.12 | $1024^{3}$ | $k\leq 10m_{r}$ | 160
0.6 | $1024^{3}$ | $k\leq m_{r}$ | 100 | 2.12 | $1024^{3}$ | no cuts | 130
0.9 | $1024^{3}$ | $k\leq m_{r}$ | 100 | 3.0 | $1024^{3}$ | $k\leq 10m_{r}$ | 130
0.2 | $2048^{3}$ | $k\leq m_{r}$ | 100 | 2.12 | $2048^{3}$ | $k\leq 10m_{r}$ | 100
0.2 | $4096^{3}$ | $k\leq m_{r}$ | 100 | 2.12 | $4096^{3}$ | $k\leq 10m_{r}$ | 100

Table 3: Summary of our data set used for the analyses. Note that the cut
$k\leq 10m_{r}$ at $\tau_{i}=0.1\tau_{c}$ is equivalent to $k\leq m_{r}$ at
$\tau_{c}$

Based on our fitting result, we select two benchmark points for further
simulations on the larger lattice. For the simulation with the fat-string type
relaxation, $\xi_{0}=0.2$ is chosen. For the one with the thermal relaxation,
$\zeta=2.12$ and $\tau_{i}=0.1\tau_{c}=0.1414$ is chosen.

###  C.2 Cosmological evolution of string network

####  C.2.1 Energy budget

Figure 21: Energy budget for axion, radial, and string modes from simulations
with different lattice sizes: $N^{3}=1024^{3}$ (solid), $N^{3}=2048^{3}$
(dashed), and $N^{3}=4096^{3}$ (dotted), which equivalently represents the
dependence on $m_{r}\Delta$.

In addition to illustrations in Section 6.1, Fig. 21 shows the dependence of
the energy budgets on the initial conditions for two benchmark simulations
using our large ensemble data from simulations on the lattice of
$N^{3}=1024^{3}$. As is evident in Fig. 21, the fraction of the string energy
density, $\rho_{s}/\rho_{\text{tot}}$
($\rho_{s}=\rho_{\text{tot}}-\rho_{a}-\rho_{r}$), tends to have a larger value
for a bigger size of $\xi_{0}$ or $\zeta$ in early times. However, the
hierarchy switches over the time, namely $\rho_{s}/\rho_{\text{tot}}$
eventually tends to have a smaller value for a bigger $\xi_{0}$ or $\zeta$ at
late times. The transition time for switching the hierarchy, or
$\log\frac{m_{r}}{H}\sim 2.5$ (fat-string pre-evolution) and
$\log\frac{m_{r}}{H}\sim 4.5$ (thermal pre-evolution), looks consistent with
the evolution of $\xi$ in Fig. 2. While the decreased fractional string energy
density, when increasing $\xi_{0}$, goes democratically into axions and radial
modes for the case with the fat-string pre-evolution (see left panel of Fig.
21), the situation, when increasing $\zeta$, in the case with the thermal pre-
evolution (see right panel of Fig. 21) looks more difficult to interpret due
to a large migration of radial modes to axions (see also more noisy spectrum
in $F$ for a smaller $\zeta$ in the right panel of Fig. 23). More detailed
study on the distinctive behaviors of energy budgets originated from different
relaxations will be interesting.

####  C.2.2 Instantaneous emission

We continue investigating the dependency of our simulation result on the
initial configurations and relaxation schemes with the large ensemble data
from simulations on the lattice of $1024^{3}$. The differential axion energy
spectrum $\partial\rho_{a}/\partial k$ and instantaneous emission $F$ for the
curves with various initial values of $\xi_{0}$ and $\zeta$ are illustrated in
Figs. 22 and 23, respectively. Overall, the instantaneous emission obtained by
the simulation with the thermal relaxation shows more noisy spectrum. While
the larger (smaller) value of $\xi$ populates more (less) pronounced high
momentum modes which cause more (less) oscillatory instantaneous emission
spectrum in the left panel of Fig. 23, the larger value of $\zeta$ appears to
populate more broad range of momentum modes (see the right panel of Figs. 22)
and the strength of the oscillation amplitudes in the corresponding
instantaneous emissions for both $\zeta=2.12$ and 1.5 look similar (although
the case with $\zeta=1.5$ is more noisy) as is evident in the right panel of
Fig. 23. The axion emission spectra in Fig. 23 look more oscillatory than the
case in [9] where the average over samples with finite-ranged initial times
was considered. It may be reasonable to assume the stochastic distribution of
the Hubble length during the phase transition, instead of being fixed to one
specific value, and consider its influence on the instantaneous emission
spectrum. It can be implemented by sampling physical string simulations over
finite-ranged initial time and average over them to derive the emission
spectrum as was done in [9]. We also performed the simulation taking into
account the stochastic distribution of the initial Hubble length, and we
confirmed that distributions become significantly smoothened. While this
effect is visibly pronounced in early times where the higher modes were not
relaxed enough, it gets diminished in later times beyond the time coverage of
$1024^{3}$, for instance, it becomes less relevant in the simulation with the
lattice size of $N^{3}=4096^{3}$.

Figure 22: The differential axion energy spectrum rescaled by $Hf_{a}^{2}$ at
the time slice $\log\frac{m_{r}}{H}=6$ for various initial conditions in two
different simulation setups: fat-string type pre-evolution (left) and thermal
pre-evolution (right).

Figure 23: The instantaneous emission $F$ at the time slice
$\log\frac{m_{r}}{H}=6$ for various initial conditions in two different
simulation setups: fat-string type pre-evolution (left) and thermal pre-
evolution (right).

Figure 24: The differential axion energy spectrum (left) and instantaneous
emission $F$ (right) at the time slice $\log\frac{m_{r}}{H}=6$ for two
situations with and without the cuts on UV modes above $m_{r}$ in the
simulation using the thermal pre-evolution.

Finally we demonstrate in Fig. 24 the impact on the differential axion energy
spectrum and the instantaneous emission from UV modes shorter than string core
length. The momentum cutoff set by the lattice spacing is higher than the
string core scale $m_{r}$ and generating those UV modes above the $m_{r}$
scale significantly distorts the spectrum as is seen in Fig. 24. Although the
initial field configurations above $m_{r}$ are not generated at the start of
the simulation, higher modes above $m_{r}$ show up later through the equations
of motion and their effects are much smaller compared to those from directly
generated UV modes from initial random fields generation.

##  Appendix D More on axion spectrum

###  D.1 Fitting instantaneous emission function

We explain how the instantaneous emission was fitted using the ansatz for the
power law fall-off profile. The spectral index $q$ for each individual
simulation is fit with the ansatz for $F(x,y)$ assuming the Gaussian
distribution with the mean $ax^{-q}$ and the standard deviation
$\sigma(x,B,p,C)=Bx^{-p}+C$ (similarly to [15]) where $x=\frac{k}{H}$ is the
momentum in the Hubble unit. The fit is performed within the finite range
$[x_{\text{IR}},\,\frac{y}{x_{\text{UV}}}]$ at a given time
$y\equiv\frac{m_{r}}{H}$, and the fit parameters $a$, $q$, $B$, $p$ and $C$
are estimated by maximizing the likelihood function,

$L(a,q,B,p,C)=\prod_{j=1}^{N_{\text{fit}}}\frac{1}{\sqrt{2\pi}\sigma(x_{j},B,p,C)}e^{-\frac{1}{2\sigma(x_{j},B,p,C)^{2}}\left(F_{j}-ax_{j}^{-q}\right)^{2}}$
(56)

where $\\{x_{j},F_{j}\\}$ are the $N_{\text{fit}}$ number of data points
within the fit interval $x_{\text{IR}}\leq x\leq y/x_{\text{UV}}$. Our ansatz
in Eq. (56) assumes no correlation among the measurements of instantaneous
power spectrum at different momenta to be consistent with string core masking
procedure in which the mode independence is assumed. The results from
individual fits are averaged to yield the spectral index for the corresponding
ensemble. They are illustrated in Figs. 14, 15, 25 and 26. The fit results of
the spectral index in Figs. 25 and 26 differs from those in Figs. 14 and 15 in
Section 7.1 only by the $\Delta\log$. With the increasing time interval from
the default value, $\Delta\log=0.25$, to $\Delta\log=0.5$ and 0.75 makes the
distributions smoother. As was illustrated in Figs. 14 and 15,
$x_{\text{UV}}\sim 10$ looks good enough to get the stable result, and the
distribution is less sensitive to $x_{\text{IR}}$.

While those in Figs. 14, 15, 25 and 26 are our fiducial fit results, the fits
of $q$ with no $Bx^{-p}$ term in $\sigma$ are carried out as well by setting
$B=0$ by hand. It leads to $30\%$ deviation in slope $q_{1}$ compared to our
fiducial result. Nevertheless, the conclusion from linear growth model fit
expecting $q\sim O(10)$ at $\log{\frac{m_{r}}{H}}\sim 70$ is not altered.

The final fit results with two different forms of the standard deviations as
well as varying UV, IR cutoffs and the time interval $\Delta\log$ are
summarized in Figs. 27 and 28.

Figure 25: The spectral index $q$ of the instantaneous emission $F$ fitted in
the interval with the varying $x_{\text{UV}}$ while $x_{\text{IR}}=15$. The
error bars are statistical.

Figure 26: The spectral index $q$ of the instantaneous emission $F$ fitted in
the interval with the varying $x_{\text{IR}}$ while $x_{\text{UV}}=10$. The
error bars are statistical. Figure 27: Our fit results on the spectral index
with log-hypothesis ($q=q_{0}+q_{1}\log\frac{m_{r}}{H}$), assuming the
standard deviation of the form $\sigma=Bx^{-p}+C$, for two types of
relaxations (fat-string and thermal pre-evolutions) while varying
$x_{\text{IR}}$ or $x_{\text{UV}}$. See Appendix D.1 for the detail. Figure
28: Our fit results on the spectral index with log-hypothesis
($q=q_{0}+q_{1}\log\frac{m_{r}}{H}$), assuming the constant standard deviation
for two types of relaxations (fat-string and thermal pre-evolutions) while
varying $x_{\text{IR}}$ or $x_{\text{UV}}$. See Appendix D.1 for the detail.

###  D.2 Analytic understanding of positive correlation

Figure 29: The correlation between the differential rate of the comoving
number of strings per Hubble with respect to the logarithmic time change,
$d\log(\xi H)/(d\log\frac{m_{r}}{H})$, and the spectral index $q$ at various
instant times from two benchmark simulations. The straight lines are major
axes of the ellipses and each ellipse denotes 1-sigma contour of the Gaussian
distribution over 100 individual runs.

Fig. 29 shows the correlation between the spectral index $q$ and $d\log(\xi
H)/(d\log\frac{m_{r}}{H})$ (the differential rate of the comoving number of
strings per Hubble with respect to time) for multiple instant time slices to
illustrate that the orientation of the correlation axes get saturated at late
times while the spectral index $q$ increases with the logarithmic time. The
ellipse represents the 1$\sigma$ contour around the central values, assuming
the Gaussian distribution of the 100 ensemble samples at each instant time. In
both simulations using the fat-string and thermal pre-evolutions, the
correlation axes, which is taken from the major axis of the ellipse, similarly
stay constant at later times for $\log>7$. This correlation can hint us on how
the spectral index might have shifted when a different initial condition
following a different $\xi$ evolution curve was chosen.

The positive correlation between $q$ and $d\log(\xi H)/(d\log\frac{m_{r}}{H})$
can be qualitatively understood by the analytic study on the axion radiation
power from sinusoidal string oscillations with small-amplitudes [32, 26]. The
time averaged radiation power per unit length by long strings aligned in
$z$-direction is given by the second harmonic at leading-order,

$\frac{\partial P}{\partial z}\simeq\frac{\beta
f_{a}^{2}}{8\lambda_{\text{str}}}\varepsilon^{4}\leavevmode\nobreak\ ,$ (57)

in terms of the wavelength of the string $\lambda_{\text{str}}$, the relative
amplitude $\varepsilon\equiv 2\pi A/\lambda_{\text{str}}$ which is the ratio
of the amplitude $A$ to its wavelength $\lambda_{\text{str}}$, and a constant
$\beta$ encoding the string geometry, for instance, $\beta=\pi^{3}/4$ for the
sinusoidal oscillation. The relative amplitude $\varepsilon$ is expected to
not vary much with respect to $\lambda_{\text{str}}$ for being a normalized
amplitude. Assuming no energy transfer to radial modes, the radiation power
can be related to the decay rate of the string energy per unit length,
$\partial E_{\text{str}}/\partial z$, via the energy conservation,

$\frac{\partial}{\partial t}\frac{\partial E_{\text{str}}}{\partial
z}=-\frac{\partial P}{\partial z}\leavevmode\nobreak\ .$ (58)

With the relation of the string energy to the string length assuming the
constant string tension $\mu_{\text{eff}}$, $\frac{\partial
E_{\text{str}}}{\partial z}=\mu_{\text{eff}}\frac{\partial
L_{\text{str}}}{\partial z}$, the back-reaction leads to shortening of string
length $L_{\text{str}}$,

$\frac{1}{L_{\text{str}}}\frac{\partial L_{\text{str}}}{\partial
t}=-\frac{1}{\partial E_{\text{str}}/\partial z}\frac{\partial P}{\partial
z}\simeq-\frac{\beta
f_{a}^{2}\varepsilon^{4}}{8\mu_{\text{eff}}\lambda_{\text{str}}}\leavevmode\nobreak\
.$ (59)

As the energy is mostly radiated into axions of wavelength
$\lambda_{\text{str}}/2$ [32], the relation in Eq. (59) indicates more rapid
decay rate of the string length for higher-frequency (or shorter wavelength)
fluctuations, which are associated with smaller spectral indices, and vice
versa.

##  Appendix E More on axion abundance

The parametric behavior of the axion abundance crucially depends on the size
of the spectral index $q$ (whether $q>1$ or $q=1$), the average axion field
values in the $f_{a}$ unit (whether $\langle a^{2}\rangle^{1/2}/f_{a}\ll 1$ or
$\langle a^{2}\rangle^{1/2}/f_{a}\gtrsim 1$), and the lower boundary of the
momentum integral (whether $k_{\text{IR}}=H$ or $H\sqrt{\xi}$). The nature of
the evolution after the scaling regime ends until the axions become
nonrelativistic crucially relies on the value of $\langle
a^{2}\rangle^{1/2}/f_{a}$. The average inter-string distance is order of
$1/(H\sqrt{\xi})$ as opposed to $H^{-1}$. The discrepancy between literature
is basically originated from different set of those factors. Therefore,
sorting out all possible parametric behaviors should be highly beneficial for
the clearer comparison. In this section, we follows the similar computation to
[14] not only filling more details, but also extending to other cases.

The axion spectrum is given by

$\frac{\partial\rho_{a}}{\partial
k}=\int^{t}dt^{\prime}\frac{\Gamma(t^{\prime})}{H(t^{\prime})}\left(\frac{R(t^{\prime})}{R(t)}\right)^{3}F\left[\frac{k^{\prime}}{H(t^{\prime})},\,\frac{m_{r}}{H(t^{\prime})}\right]\leavevmode\nobreak\
,$ (60)

where $k^{\prime}=(R(t)/R(t^{\prime}))k$ and the normalized instantaneous
emission in the interval $x=[x_{\text{IR}},\,y/x_{\text{UV}}]$ is given by,
depending on $q$,

$\begin{split}F\left[x,\,y\right]=&\left\\{\begin{array}[]{lll}\displaystyle\frac{q-1}{x_{\text{IR}}^{1-q}-\left(\frac{y}{x_{\text{UV}}}\right)^{1-q}}\,\frac{1}{x^{q}}&\text{for}&q>1\\\\[25.0pt]
\hskip
10.11775pt\displaystyle\frac{1}{\log\frac{y}{x_{\text{UV}}x_{\text{IR}}}}\,\frac{1}{x^{q}}&\text{for}&q=1\leavevmode\nobreak\
.\end{array}\right.\end{split}$ (61)

The transfer rate from strings to axions $\Gamma$ is approximated at late
times, $\Gamma\sim\xi\mu_{\text{eff}}/t^{3}\sim 8\pi
H^{3}f_{a}^{2}\xi\log\frac{m_{r}}{H}$ [9] where
$\langle\gamma\rangle\mu_{0}=\pi f_{a}^{2}$ in Eq. (11) was used in the second
relation.

###  E.1 Estimation for $q>1$

Following similar steps in [14], we show that modifying IR momentum cutoff can
induce non-negligible change in the overall size of the axion abundance due to
the slightly changed parametric behavior. Using the approximation $\xi\sim
c_{1}\log\frac{m_{r}}{H}$ at late times, the axion spectrum in Eq. (60) can be
rewritten as

$\frac{\partial\rho_{a}}{\partial k}(k,\,t)=\frac{4\pi
f_{a}^{2}c_{1}m_{r}}{y^{3/2}}\int^{y}_{y_{0}}dy^{\prime}\frac{(\log
y^{\prime})^{2}}{{y^{\prime}}^{1/2}}F\left[\left(\frac{y^{\prime}}{y}\right)^{1/2}x,\,y^{\prime}\right]\leavevmode\nobreak\
,$ (62)

where $x=\frac{k}{H}$, $y=\frac{m_{r}}{H(t)}$, and the integration over
$t^{\prime}$ in the interval $[t_{0},\,t]$ was converted in terms of
$y^{\prime}=\frac{m_{r}}{H(t^{\prime})}$ with the interval $[y_{0},\,y]$. We
take into account the modified IR cutoff by setting
$x_{\text{IR}}=x_{0}\sqrt{\xi}$. $x_{\text{UV}}=1$ is chosen similarly to [14]
and the result is insensitive to this choice. The instantaneous emission in
Eq. (61) for $q>1$ can be approximated as