# I Introduction

## I Introduction

Cosmic strings could be formed in consequence of the cosmological phase transition at the very early universe . They are divided into global and local(gauged) strings according to the property of the broken symmetry. Though they have similar properties, such as they are line objects with the false vacuum energy, intercommute at crossing  and so on, there are also many differences between them. While a local string has two cores comprised of a magnetic flux core and a scalar field core, a global string has only a scalar field core. For local strings, the gradient energy of the scalar field can be canceled out by the gauge field far from the core so that the core is well localized and the vacuum energy of the core is dominant. Therefore, the Nambu-Goto action is adequate to follow the evolution of the local string network except at crossing. On the other hand, there is no gauge field in the global string model so that the total energy of global strings is dominated by not the vacuum energy of the core but the gradient energy of the scalar field, namely, the NG boson field. So, we need consider not only the core but also the associated NG boson field and the coupling between them in contrast with the case of local strings. Thus we must use the Kalb-Ramond action instead of the Nambu-Goto action as an effective action .

The dynamics of the local string network has been examined by use of the Nambu-Goto action. Kibble  first proposed the so-called one scale model, where the behavior of the system can be characterized by only a parameter, namely, the scale length, $`L`$, and an unknown loop-production function. He showed that either the local string network goes into the scaling regime where the scale length, $`L`$, grows with the horizon distance , or, the universe becomes string dominated. In the scaling regime, infinite strings intercommute to produce closed loops so that at any time the number of strings stretching across the horizon distance within each horizon volume is almost constant and produced loops decay through radiating gravitational waves . Bennett  developed the Kibble’s one scale model and showed that unless most of produced loops self-intersect and fragment into smaller loops with the typical length smaller than $`L`$, the reconnection rate is large enough to prevent scaling. Mitchell and Turok  studied the statistical mechanics of the string network in the flat spacetime and found that the equilibrium distribution of the string network is dominated by the smallest loops allowed, which suggested that strings tend to break into very small pieces in the expanding universe. Albrecht and Turok  modeled the network as the hot body radiation where loops are radiated from the long strings and showed that the scaling solution is inevitable. However, the application of the flat spacetime statistical mechanics to the string dynamics in the expanding universe is not necessarily justified. Thus, numerical simulations are unavoidable to decide whether the local string network settles down to the scaling regime or not.

Albrecht and Turok  gave the first numerical investigation for the scaling property. Later three groups  improved numerical codes and all groups concluded that the large scale behavior of the local string network goes into the scaling regime with the scaling parameter $`\xi `$ equal to $`(50\pm 25)`$ , $`(13\pm 2.5)`$ , and $`(16\pm 4)`$  in the radiation-dominated era. Though Albrecht and Turok  found that the loop distribution function also scales, the higher resolution simulations performed by the other two groups  showed that long strings have significant small structures and loops are typically produced at scales much smaller than the horizon distance, which is close to the cut-off scale. In response, Austin et al.  proposed the three scale model with $`\xi `$, the step length $`\overline{\xi }`$, and the scale $`\zeta `$ describing the small scale structure. They found that $`\xi `$ and $`\overline{\xi }`$ grow with the horizon distance but $`\zeta `$ begin growing only when the gravitational back reaction becomes in effect with $`\zeta /\xi 10^4`$. However, Vincent et al  recently claimed from the flat spacetime simulations that instead of loop production due to intercommutation, long strings directly emit massive particles so that the dominant energy loss mechanism of local strings is not gravitational radiation but particle production.

Thus, though there is a consensus that the large scale structure of the local string network obeys the scaling solution, the loop production and the dominant energy loss mechanism are now in dispute. This is mainly because inclusion of gravitational radiation to the numerical simulations is impossible and gravitational radiation is so weak that kinks live for a long time.

On the other hand, the evolution of the global string network has been less studied. Its manifestation, however, is essential to estimation of abundances of relic axions radiated from axionic strings , which may be the cold dark matter. Also, scaling property is indispensable for the scenario where global strings become the seed of the structure formation of the universe and produce the cosmic microwave background anisotropy . Thus, it is important to clarify the dynamics of the global string network, in particular, whether it enters the scaling regime like the local string network.

So far, the result for the local string network has been applied to the global string network because global strings also intercommute with the probability of the order unity . But global strings have the associated Nambu-Goldstone bosons, which lead to long-range forces between strings and become the dominant energy loss mechanism of global strings. Therefore, instead of the Nambu-Goto action we have to use the Kalb-Ramond action , which is comprised of three components, the Nambu-Goto action, the kinetic term of the associated Nambu-Goldstone fields, and the coupling term between them. But the Kalb-Ramond action has difficulty of logarithmic divergence due to self-energy of the string. Dabholkar and Quashnock  solved this difficulty by the similar prescription as given by Dirac in the electromagnetic system , where divergence is renormalized into the electron mass. They gave the renormalized equation of motion for a global string comprised of the free part derived from the Nambu-Goto action with the damping term, which becomes negligible for a circular loop in the cosmological scale where $`\mathrm{ln}(t/\delta )𝒪(100)`$. Then they concluded that the global string network can be well approximated by the motion of the Nambu-Goto action. However, as shown by Battye and Shellard  though the calculation is done in the flat spacetime, the kinks on long string are substantially rounded due to the backreaction of NG boson radiation, which may significantly affect the small scale structure (if at all) of the global network system. Also, the above approach cannot include the long-range force between long strings, which may decrease the energy density of long strings. Thus, the examination of the dynamics of global strings by use of the Kalb-Ramond action is not yet complete.

In the previous paper , we manipulated the equation of motion of the complex scalar field instead of the Kalb-Ramond action. Then, we showed that the large scale behavior of the global string network goes into the scaling regime and $`\xi `$ for the global string network becomes $`𝒪(1)`$, which is significantly smaller than that for the local string network.

In this paper, we give a comprehensive analysis of the evolution of the global string network based on the model adopted in . We show the scaling of the global string network under more general situations and investigate the dependence of $`\xi `$ on boundary conditions and some parameters. Then, the loop distribution function is given in order to decide whether the small scale structure exists like the local string network.

The paper is organized as follows: In the next section, we give the formulation of numerical simulations. In the section 3, we give the method of the identification of string segments and closed loops. Then, the scaling parameter $`\xi `$ and the loop distribution function are given. Finally, we discuss our results and give conclusions.

## II Formulation of numerical simulations

We consider the following Lagrangian density for a complex scalar field $`\mathrm{\Phi }(x)`$,

$$[\mathrm{\Phi }]=g_{\mu \nu }(^\mu \mathrm{\Phi })(^\nu \mathrm{\Phi })^{}V_{\mathrm{eff}}[\mathrm{\Phi },T],$$
(1)   
where $`g_{\mu \nu }`$ is identified with the Robertson-Walker metric and the effective potential $`V_{\mathrm{eff}}[\mathrm{\Phi },T]`$ is given by

$$V_{\mathrm{eff}}[\mathrm{\Phi },T]=\frac{\lambda }{2}(\mathrm{\Phi }\mathrm{\Phi }^{}\eta ^2)^2+\frac{\lambda }{3}T^2\mathrm{\Phi }\mathrm{\Phi }^{}.$$
(2)   
For $`T>T_c=\sqrt{3}\eta `$, the potential $`V_{\mathrm{eff}}[\mathrm{\Phi },T]`$ has a minimum at $`\mathrm{\Phi }=0`$, and the $`U(1)`$ symmetry is restored. On the other hand, new minima $`|\mathrm{\Phi }|_{\mathrm{min}}=\eta \sqrt{1(T/T_c)^2}`$ appear and the symmetry is broken for $`T<T_c`$ (Fig. 1). In this case the phase transition is of second order.

The equation of motion is given by

$$\ddot{\mathrm{\Phi }}(x)+3H\dot{\mathrm{\Phi }}(x)\frac{1}{R(t)^2}^2\mathrm{\Phi }(x)=V_{\mathrm{eff}}^{}[\mathrm{\Phi },T],$$
(3)   
where the prime represents the derivative $`/\mathrm{\Phi }^{}`$ and $`R(t)`$ is the scale factor. The Hubble parameter $`H=\dot{R}(t)/R(t)`$ and the cosmic time $`t`$ are given by

$`H^2={\displaystyle \frac{8\pi }{3m_{\mathrm{pl}}^2}}{\displaystyle \frac{\pi ^2}{30}}g_{}T^4,t={\displaystyle \frac{1}{2H}}={\displaystyle \frac{\xi }{T^2}},`$ (4)   
where $`m_{\mathrm{pl}}=1.2\times 10^{19}`$GeV is the Plank mass, $`g_{}`$ is the total number of degrees of freedom for the relativistic particles, and radiation domination is assumed. We define the dimensionless parameter $`\zeta `$ as

$$\zeta \frac{\xi }{\eta }=\left(\frac{45M_{\mathrm{pl}}^2}{16\pi ^3g_{}\eta ^2}\right)^{1/2}.$$
(5)   
In our simulation, we take $`\zeta =10`$, which corresponds to $`\eta (10^{15}10^{16})`$ GeV with $`g_{}=1000`$ and later investigate $`\zeta `$ dependence on the result. The energy density at each lattice point is written as

$$\rho (x)=\dot{\mathrm{\Phi }}(x)\dot{\mathrm{\Phi }}^{}(x)+\frac{1}{R(t)^2}\mathrm{\Phi }(x)\mathrm{\Phi }^{}(x)+V_{\mathrm{eff}}[\mathrm{\Phi },T].$$
(6)   
We take the initial time $`t_i=t_c/4`$ and the final time $`t_f=75t_i=18.75t_c`$, where $`t_c`$ is the epoch $`T=T_c`$. Since the $`U(1)`$ symmetry is restored at the initial time $`t=t_i`$, we adopt as the initial condition the thermal equilibrium state with the mass squared,

$$m^2=\frac{d^2V_{\mathrm{eff}}[|\mathrm{\Phi }|,T]}{d|\mathrm{\Phi }|^2}|_{|\mathrm{\Phi }|=0},$$
(7)   
which is the inverse curvature of the potential at the origin at $`t=t_i`$. In the thermal equilibrium state, $`\mathrm{\Phi }`$ and $`\dot{\mathrm{\Phi }}`$ are Gaussian distributed with the correlation functions,

$`\beta |\mathrm{\Phi }(𝒙)\mathrm{\Phi }^{}(𝒚)|\beta _{\mathrm{equal}\text{-}\mathrm{time}}`$ $`=`$ $`{\displaystyle \frac{d𝒌}{(2\pi )^3}\frac{1}{2\sqrt{𝒌^2+m^2}}\mathrm{coth}\frac{\beta \sqrt{𝒌^2+m^2}}{2}e^{i𝒌(𝒙𝒚)}},`$ (8)   
$`\beta |\dot{\mathrm{\Phi }}(𝒙)\dot{\mathrm{\Phi }}^{}(𝒚)|\beta _{\mathrm{equal}\text{-}\mathrm{time}}`$ $`=`$ $`{\displaystyle \frac{d𝒌}{(2\pi )^3}\frac{\sqrt{𝒌^2+m^2}}{2}\mathrm{coth}\frac{\beta \sqrt{𝒌^2+m^2}}{2}e^{i𝒌(𝒙𝒚)}}.`$ (10)   
The functions $`\mathrm{\Phi }(𝒙)`$ and $`\dot{\mathrm{\Phi }}(𝒚)`$ are uncorrelated for $`𝒙𝒚`$. We generate these fields for the initial condition in the momentum space, because the corresponding fields $`\stackrel{~}{\mathrm{\Phi }}(𝒌)`$ and $`\stackrel{~}{\dot{\mathrm{\Phi }}}(𝒌)`$ are uncorrelated there. Then these fields are inverse Fourier transformed into the configuration spaces by the FFT formalism.

Hereafter we measure the scalar field in units of $`t_i^1`$, $`t`$ and $`x`$ in units of $`t_i`$, and the energy density in units of $`t_i^4`$. The equation of motion and the total energy density are given by

$`\ddot{\mathrm{\Phi }}(x)+{\displaystyle \frac{3}{2t}}\dot{\mathrm{\Phi }}(x){\displaystyle \frac{1}{t}}^2\mathrm{\Phi }(x)=\left(|\mathrm{\Phi }|^2+{\displaystyle \frac{\zeta ^2}{36t}}{\displaystyle \frac{\zeta ^2}{144}}\right)\mathrm{\Phi }^{},`$ (12)   
$`\rho (x)=\dot{\mathrm{\Phi }}(x)\dot{\mathrm{\Phi }}^{}(x)+{\displaystyle \frac{1}{t}}\mathrm{\Phi }(x)\mathrm{\Phi }^{}(x)+{\displaystyle \frac{1}{2}}\left(|\mathrm{\Phi }|^2{\displaystyle \frac{\zeta ^2}{144}}\right)^2+{\displaystyle \frac{\zeta ^2}{36t}}|\mathrm{\Phi }|^2,`$ (13)   
where $`\lambda `$ is set to unity for brevity. The scale factor $`R(t)`$ is normalized as $`R(1)=1`$.

We perform the simulations in four different sets of lattice sizes and spacings (See TABLE I.): (1) $`128^3`$ lattices with the physical lattice spacing $`\delta x_{\mathrm{phys}}=2\sqrt{3}t_iR(t)/25`$. (2) $`64^3`$ lattices with $`\delta x_{\mathrm{phys}}=4\sqrt{3}t_iR(t)/25`$. (3) $`256^3`$ lattices with $`\delta x_{\mathrm{phys}}=\sqrt{3}t_iR(t)/25`$. (4) $`256^3`$ lattices with $`\delta x_{\mathrm{phys}}=2\sqrt{3}t_iR(t)/25`$. In all cases, the time step is taken as $`\delta t=0.01t_i`$. In the case (1), the box size is nearly equal to the horizon volume $`(H^1)^3`$ and the lattice spacing to a typical width $`\delta 1.0/(\sqrt{2}\eta )`$ of a string at the final time $`t_f`$. Furthermore, in order to investigate the dependence of $`\zeta `$, we arrange the case with $`\zeta =5`$, (7) $`128^3`$ lattices with the physical lattice spacing $`\delta x_{\mathrm{phys}}=2\sqrt{6}t_iR(t)/25`$. In this case we follow the time development of the system until the final time, $`t_f=150t_i=37.5t_c`$, when the box size is nearly equal to the horizon volume $`(H^1)^3`$ and the lattice spacing to a typical width of a string. For each case, we simulate the system from 10 ((2), (3), (4), and (7)) or 300 ((1)) different thermal initial conditions. Since the simulation box is larger than the horizon volume even at the final time of the simulation, we adopt the periodic boundary condition. But, under the periodic boundary conditions, there exists no infinite string so that it is possible that string completely disappears in the simulation box. Therefore, in order to verify the dependence of boundary conditions, we also simulate the cases in TABLE II under the reflective boundary condition where $`^2\mathrm{\Phi }(x)`$ on the boundary points disappears <sup>*</sup><sup>*</sup>*Note that this condition is different from the open boundary condition where $`\mathrm{\Phi }(x)`$ on the boundary points disappears. Under the open boundary condition, the string feels attractive forces from the boundary so that the number of the strings tends to decrease than that in the real universe as under the periodic boundary condition..

Using the second order leap-frog method and the Crank-Nicholson scheme, the discretized equation of motion reads

$`\dot{\mathrm{\Phi }}_{𝒊,n+1/2}`$ $`=`$ $`{\displaystyle \frac{1}{1+\frac{3\delta t}{4t}}}\left[\left(\mathrm{\hspace{0.17em}1}{\displaystyle \frac{3\delta t}{4t}}\right)\dot{\mathrm{\Phi }}_{𝒊,n1/2}+{\displaystyle \frac{\delta t}{t}}^2\mathrm{\Phi }_{𝒊,n}\delta t\left\{|\mathrm{\Phi }_{𝒊,n}|^2+{\displaystyle \frac{\zeta ^2}{36t}}{\displaystyle \frac{\zeta ^2}{144}}\right\}\mathrm{\Phi }_{𝒊,n}^2\right],`$ (14)   
$`\mathrm{\Phi }_{𝒊,n+1}`$ $`=`$ $`\mathrm{\Phi }_{𝒊,n}+\delta t\dot{\mathrm{\Phi }}_{𝒊,n+1/2},`$ (15)   
$`^2\mathrm{\Phi }_{𝒊,n}`$ $``$ $`{\displaystyle \underset{s=x,y,z}{}}{\displaystyle \frac{\mathrm{\Phi }_{i_s+1_s,n}2\mathrm{\Phi }_{i_s,n}+\mathrm{\Phi }_{i_s1_s,n}}{(\delta x)^2}},`$ (16)   
where $`𝒊`$ represents spatial index and $`n`$ temporal one.

## III Results

In order to judge whether the global string network relaxes into the scaling regime, we give time development of $`\xi `$, which is defined as

$$\rho =\xi \mu /t^2,$$
(17)   
where $`\mu 2\pi \eta ^2\mathrm{ln}(t/(\delta \xi ^{1/2}))`$ is the average energy per unit length of global strings.

### A String identification

Before obtaining $`\xi `$ and the loop distribution function, we must identify the string segment. Since spacetime is discretized in our simulations, a point with $`\mathrm{\Phi }=0`$ corresponding to a string core is not necessarily situated at a lattice point. In the worst case, a point with $`\mathrm{\Phi }=0`$ lies at the center of a plaquette. Therefore, we use a static cylindrically-symmetric solution, which is obtained by solving the equation

$$\frac{d^2f}{dr^2}+\frac{1}{r}\frac{df}{dr}\frac{f}{r^2}V_{\mathrm{eff}}^{}[f,T]=0,$$
(18)   
with $`\mathrm{\Phi }(r,\theta )f(r)e^{i\theta }`$ and the winding number $`n=1`$. The boundary conditions are given by

$`f(r)`$ $``$ $`|\mathrm{\Phi }|_{\mathrm{min}},(r\mathrm{}),`$ (19)   
$`f(0)`$ $`=`$ $`0.`$ (20)   
We require that a lattice is identified with a part of a string core if the potential energy density there is larger than that corresponding to the field value of a static cylindrically-symmetric solution at $`r=\delta x_{\mathrm{phys}}/\sqrt{2}`$. Then only one lattice within a section of a straight string core is identified with a string segment except the case where a point with $`\mathrm{\Phi }=0`$ lies at the center of a plaquette. Of course, the real string is not exactly straight but bent and more complex. But, as seen in Fig. 2, our identification is good and only one lattice within a section of a string core is identified with a part of a string core. In order to solve the above equation of motion, we have to use the standard shooting technique, which is the repeated process so that it is bothersome to follow the above procedure for each time step. Instead, we obtain the solutions every $`500\delta t`$ and make the fitting formula. For intermediate time steps, we judge whether a lattice point belongs to the string segment by comparing the potential energy obtained from the simulation and that from the fitting formula. Thus, by counting the number of the lattices identified with a part of a string core, we can evaluate the total length of strings within the horizon volume, $`(2t)^3`$, from which the energy density can be obtained.

### B Scaling property

Time development of $`\xi `$ in the cases from (1) to (4) is described in Fig. 3. We find that after some relaxation period, $`\xi `$ becomes a constant irrespective of time with (1)$`0.99\pm 0.09`$, (2)$`0.97\pm 0.04`$, (3)$`0.90\pm 0.06`$, and (4)$`0.90\pm 0.05`$. They all are consistent within the standard deviation. Thus we can conclude that a global string network relaxes into scaling regime in the radiation domination. We also show time development of $`\xi `$ in the case (7) in Fig. 4. Then $`\xi `$ asymptotically becomes a constant, $`0.88\pm 0.07`$, which is also consistent with the above all cases with $`\zeta =10`$ within the standard deviation. Hence we can also conclude that $`\zeta `$ does not change the essential result. Note that the standard deviation in the case (2) is much smaller than the other cases because the box in the case (2) includes more horizon volumes at each time. Also, $`\xi `$ seems to oscillate at the early epoch because the homogeneous mode of the field oscillates in the radial direction of the potential, which rapidly decays. In fact, the period of the oscillation coincides with $`2\pi `$ times the inverse mass at the potential minimum.

Fig. 5 represents the results under the reflective boundary condition, where $`\xi `$ becomes a constant irrespective of time with (1)$`1.77\pm 0.03`$, (2)$`1.57\pm 0.04`$, (3)$`2.00\pm 0.05`$, (4)$`1.30\pm 0.02`$, (5)$`1.18\pm 0.03`$, and (6)$`1.03\pm 0.03`$. Though all results are consistent with that under the periodic boundary condition within the factor two, the former tends to become larger than the latter. This is because under the reflective boundary condition, strings are repulsed by the boundary and the string near the boundary intercommutes less often than that near the center of the simulation box because the partner to intercommute only lies in the inner direction of the boundary. The results in large box simulations, (4)-(6), are converging so that it is safe to say that if we take the box size larger than $`2^3`$ times the horizon volume, the reflective boundary condition has no effect on the results.

### C Loop distribution

We also investigate the loop distribution. Since in our simulations it is judged by the potential energy whether a lattice point is a part of the string, the string is not necessarily continuous. Therefore, we identify a closed loop as follows; First we select a lattice point belonging to the string segment. Then we connect it with the nearest neighbor among lattice points belonging to the string segment. We proceed this process one after another until the connection returns to the starting lattice point.

Kibble’s one scale model predicts the loop distribution as

$$n_l(t)=\frac{\nu }{t^{\frac{3}{2}}(l+\kappa t)^{\frac{5}{2}}},$$
(21)   
where $`\nu `$ is a constant, $`l`$ is the length of a loop, and the log-dependence of $`\mu `$ is neglected. Different from a local string, the dominant energy loss mechanism for global strings is radiation of the associated Nambu-Goldstone field . We define the radiation power, $`P`$, as $`P=\kappa \mu `$ where $`\kappa `$ is a constant. An example of decay of a closed loop is depicted in Fig. 6.

We determine whether the loop distributions in the simulation coincides with the above predicted loop distribution function. The loop distributions in the case (1) at different times ($`t=45,55,65,75t_i`$) are described in Fig. 7. Since long strings are rare, we cut the length of loops into bins with the width 5$`\times \delta x`$. Also, we divide 300 realizations into 6 groups comprised of 50 realizations and we summed the number of loops over 50 realizations for each groups. The dot represents the number of loops averaged over 6 groups and the dash line represents the standard deviation. They can be simultaneously fitted with the above formula if one takes $`\kappa 0.535`$ and $`\nu 0.0865`$. Fittings for $`\kappa `$ and $`\nu `$ are also given in Fig. 8. Thus, the loop production function as well as the large scale behavior of the string scales together for the global string network.

As an another evidence for the scaling of the loop production, we consider the NG boson spectrum radiated from strings. If loops are formed at scales much smaller than the horizon distance, there should be significant power of radiated NG bosons for modes corresponding to the scale at which loops are formed. For the purpose, we first represent the complex field $`\mathrm{\Phi }(t,𝒙)`$ in terms of the radial mode $`\chi (t,𝒙)`$ and the NG boson field $`\alpha (t,𝒙)`$ as

$$\mathrm{\Phi }(t,𝒙)=\left[\eta +\frac{\chi (t,𝒙)}{\sqrt{2}}\right]\mathrm{exp}\left(\frac{i\alpha (t,𝒙)}{\sqrt{2}\eta }\right).$$
(22)   
Then the kinetic energy density of NG bosons is given by

$`{\displaystyle \frac{1}{2}}\dot{\alpha }(t,𝒙)^2`$ $`=`$ $`{\displaystyle \frac{\eta ^2}{|\mathrm{\Phi }(t,𝒙)|^4}}`$ (24)   
$`\times \left[\mathrm{Im}\mathrm{\Phi }(t,𝒙)\mathrm{Re}\dot{\mathrm{\Phi }}(t,𝒙)+\mathrm{Re}\mathrm{\Phi }(t,𝒙)\mathrm{Im}\dot{\mathrm{\Phi }}(t,𝒙)\right]^2.`$   
One may wonder if the power spectrum obtained by the Fourier transformation of the above kinetic energy density is what we want. But it includes the contribution from NG bosons formed at the symmetry breaking as well as that radiated from strings. Also, the decomposition of the field into the radial and phase modes is inadequate for the lattice point near the string segment.

Then, we evaluate the average energy density of NG bosons radiated in the period between $`t_1`$ and $`t_2`$, $`\overline{\rho }[t_1,t_2]`$. For the purpose, we subtract the redshifted kinetic energy at $`t_1`$ from the kinetic energy at $`t_2`$ since emitted NG bosons damp like radiation. Thus, $`\overline{\rho }[t_1,t_2]`$ is given by

$`\overline{\rho }[t_1,t_2]`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle d^3𝒙\rho [t_1,t_2;𝒙]}`$ (25)   
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle d^3𝒙\left[\frac{1}{2}\dot{\alpha }(t_2,𝒙)^2\frac{1}{2}\dot{\alpha }(t_1,𝒙)^2\left(\frac{t_1}{t_2}\right)^2\right]}`$ (26)   
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \frac{d^3𝒌}{(2\pi )^3}\left[\frac{1}{2}|\dot{\alpha }_𝒌(t_2)|^2\frac{1}{2}|\dot{\alpha }_𝒌(t_1)|^2\left(\frac{t_1}{t_2}\right)^2\right]}`$ (27)   
$``$ $`{\displaystyle \frac{d^3𝒌}{(2\pi )^3}\stackrel{~}{\rho }_𝒌[t_1,t_2]}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dk}{2\pi ^2}}\rho _k[t_1,t_2],`$ (28)   
where $`V`$ is the simulation volume and $`\dot{\alpha }(t,𝒙)=\frac{d^3𝒌}{(2\pi )^3}\dot{\alpha }_𝒌(t)\mathrm{exp}(i𝒌𝒙)`$.

Furthermore, to avoid contamination of string cores to the spectrum of emitted axions, we divide the simulation box into 8 cells and stock the field data of a cell if there are no string cores in that cell between $`t_1`$ and $`t_2`$. Over all such cells, we average power spectra of kinetic energy of axions obtained through Fourier transformation. We follow the above procedure between $`t_1=65t_i`$ and $`t_2=75t_i`$ for the case (3) under the periodic boundary condition, which is the highest resolution simulation, but with the zero temperature potential $`V_{\mathrm{eff}}[\mathrm{\Phi },T=0]`$ after $`t=20t_i`$ because the decomposition of the field becomes well-defined. One may suspect that this spectrum is dominated by the kinetic energy of NG bosons associated with global strings (string NG bosons) rather than that of free NG bosons radiated from global strings. But this is incorrect for the following reason; The total energy of string NG bosons is almost as much as that of free NG bosons. However, the contribution to the energy of string NG bosons is dominated by the gradient energy because the kinetic energy of string NG bosons is much smaller that their gradient energy by the factor $`v^2`$ ($`v`$ is the velocity of the string core and $`v1`$ except just before the disappearance of the loop). On the other hand, the gradient and the kinetic energy of free NG bosons is the same. Therefore the kinetic energy of free NG bosons is much larger than that of string NG bosons, that is, our spectrum is dominated by the former. Also, even if it contributed, the kinetic energy of string NG bosons would decay in proportion to $`t^2`$, so that its contribution has been removed from the spectrum in our method as done in Eq. (28).

The result is depicted in Fig. 9, which has already been given in  in the different context. The spectrum of emitted NG bosons is highly peaked at the horizon scale. Thus, loops are formed not at much smaller scale than the horizon distance but around the horizon scale.

## IV Discussions and conclusions

In this paper we gave a comprehensive investigation on the evolution of the global string network in the radiation dominated universe by use of the numerical simulation based on the complex scalar field model which spontaneously breaks the U(1) symmetry.

In order to decide whether the large scale behavior of the global string network goes into the scaling regime, we followed time development of the scaling parameter $`\xi `$ which characterizes the average energy density of global strings. We found that $`\xi `$ is almost constant irrespective of cosmic time under both the periodic and the reflective boundary condition. All the results are consistent within the factor two. But $`\xi `$ under the reflective boundary condition tend to become larger than that under the periodic boundary condition. This can be understood as follows; Under the periodic boundary condition, there are no infinite strings and strings with two boundary points on the opposite planes always intercommute with the partner so that $`\xi `$ tends to become small. On the other hand, under the reflective boundary condition, strings are repulsed by the boundary. Furthermore the string near the boundary intercommutes less often than that near the center of the simulation box because the partner to intercommute only lies in the inner direction of the boundary. Thus, $`\xi `$ tends to become large. Therefore, $`\xi `$ in the real world should lie between that under the periodic boundary condition and that under the reflective boundary condition. Considering all the cases under the periodic boundary condition and the cases (4)-(6) of the large box simulations under the reflective boundary condition, it is safe to say that $`\xi (0.91.3)`$.

You should note that $`\xi `$ is much smaller than that of a local string, which is of the order of ten. This is mainly because global strings can intercommute more often than local strings since an attractive force proportional to the inverse separation works between a global string and a global anti-string.

We have also investigated the loop distribution. It can be well fitted to that predicted by the one scale model if we take $`\nu 0.0865`$ and $`\kappa 0.535`$. Thus, the loop production grows with the horizon distance and we did not observe the small scale structure. This is because the Nambu-Goldstone(NG) boson radiation from strings is so efficient that the small scale structures on strings are damped out. The damping scale is typically $`\kappa t`$, which is near the horizon distance in contrast with $`G\mu tt`$ (G : the gravitational constant) if at all for the local string network. Bennett  showed that unless produced loops with the length of the horizon distance often self-intersect and fragment into smaller loops with the typical length smaller than the horizon distance, the reconnection rate is large enough to prevent scaling. In the global string network, parent loops do not fragment into smaller loops with the typical length smaller than the horizon distance, instead they rapidly shrink through radiating NG bosons so that the reconnection rate of parent large loops becomes small and scaling should be kept.

### Acknowledgments

The author is grateful to Jun’ichi Yokoyama and Masahiro Kawasaki for useful discussions. This work was partially supported by the Japanese Grant-in-Aid for Scientific Research from the Monbusho, Nos. 10-04558.