text:
# Direct neutrino mass search

## 1 Introduction

Recent experimental results from atmospheric and solar neutrinos give strong evidence that neutrinos oscillate from one flavor state into another. A neutrino of one specific flavor eigenstate $`\nu _\alpha `$ is a non-trivial superposition of different neutrino mass states $`\nu _\mathrm{i}`$ 111Assuming CPT-invariance, here we do not distinguish between masses of neutrinos $`m(\nu _\mathrm{i})`$ and the corresponding antineutrinos $`m(\overline{\nu }_\mathrm{i})`$.. This fact is described by an unitarity mixing matrix $`U_{\alpha i}`$

$$\nu _\alpha =\underset{i}{}U_{\alpha i}\nu _\mathrm{i}.$$
(1)
Future oscillation experiments will determine the elements $`U_{\alpha i}`$ with great precision.

However, $`\nu `$–oscillations experiments do not yield the values of the neutrino masses, they are only sensitive to differences between squared neutrino masses $`\mathrm{\Delta }m_{\mathrm{ij}}^2=|m^2(\nu _\mathrm{i})m^2(\nu _\mathrm{j})|`$ 222Neutrino oscillation effects involving matter effects are in principle able to determine the sign of $`m^2(\nu _\mathrm{i})`$-$`m^2(\nu )`$.. The values $`\mathrm{\Delta }m_{\mathrm{ij}}^2`$ from oscillation experiments only give lower limits on neutrino masses

$$\mathrm{max}(m(\nu _\mathrm{i}),m(\nu _\mathrm{j}))\sqrt{\mathrm{\Delta }m_{\mathrm{ij}}^2}.$$
(2)
On the other hand, if the absolute value of one mass eigenstate $`\nu _\mathrm{i}`$ is known, all other neutrino masses can be reconstructed with the help of the differences of the squared neutrino masses, if the signs of the different $`m^2(\nu _\mathrm{i})m^2(\nu _\mathrm{j})`$ values are known.

The information on the neutrino mass scale is of crucial importance for particle physics as well as for astrophysics and cosmology. The neutrino mass states can be arranged in a hierarchical way like the charged fermions. This would mean that the different neutrino masses are essentially governed by the square root of $`\mathrm{\Delta }m_{\mathrm{ij}}^2`$. On the other hand the neutrino masses could be quasi-degenerate with about the same value – e.g. a few tenth of an eV – and small mass differences between the different states to explain the oscillation signal. Which neutrino mass scenario 333Of course any neutrino mass scenario in between the hierarchical and the quasi-degenerate pattern would be possible. is right is a very important information in order to find the right theory beyond the Standard Model. Big bang theories predict, that a huge abundance of relic neutrinos of all flavors similar to the photons of the cosmic microwave background radiation (CMBR) exists in the universe. An average neutrino mass of 1 eV would contribute to the energy and matter distribution of the universe by 8 % in units of the critical density Therefore any neutrino mass of this order would contribute significantly to the missing dark matter and would influence the evolution of the universe and its structure formation.

Therefore, information on the neutrino masses can also be deduced from cosmology. Independent information on neutrino masses comes from the measurement of the arrival time distribution of neutrinos emitted in a supernova explosion, the present limit from SN1987a is $`m(\nu _\mathrm{e})<23\mathrm{e}V`$  (see section 2 for the remarks on $`m(\nu _\mathrm{e})`$).

Direct information on neutrino masses can be obtained by laboratory experiments using two different approaches: the investigation of the decay kinematics of weak decays and the search for neutrinoless double $`\beta `$ decay. Both methods give complementary information on the neutrino masses $`m(\nu _\mathrm{i})`$.

This paper is organized as follows: Section 2 describes the direct neutrino mass measurements, section 3 concentrates on tritium $`\beta `$ decay experiments. In section 4 the options for future neutrino mass searches is discussed. The planned next generation tritium $`\beta `$ decay experiment KATRIN is presented in section 5. Section 6 gives a short summary.

## 2 Decay Kinematics of Weak Decays

The investigations of the kinematics of weak decays is based on measurements of the charged decay products. Using energy and momentum conservation the missing neutrino mass can be reconstructed from the kinematics of the charged particles. The part of the phase space which is most sensitive to the neutrino mass is the one which corresponds to the emission of a non-relativistic massive neutrino. Therefore decays releasing charged particles with a small free kinetic energy are preferred.

In principle, a kinematical neutrino mass measurement yields information on the different mass eigenstates $`m(\nu _\mathrm{i})`$, since it performs a projection on energy and mass. But usually the different neutrino mass eigenstates cannot be resolved by the experiment. Therefore an average over neutrino mass eigenstates is obtained which is specific for the flavor of the weak decay and hence termed $`m(\nu _\mathrm{e})`$, $`m(\nu _\mu )`$ or $`m(\nu _\tau )`$, respectively444This average value is not a unique quantity but depends also on how the experiment is analyzed, in particular if it is done under the assumption of a single neutrino mass state. Considering however the small differences of squared masses $`\mathrm{\Delta }m_{\mathrm{ij}}^2`$ obtained by neutrino oscillation experiments and comparing them to the experimental resolutions of the present kinematic measurements, this question appears to be of a more academic nature.. This fact will be discussed in more detail for the case of the muon neutrino $`\nu _\mu `$.

### 2.1 $`m(\nu _\mu )`$

The muon neutrino mass $`m(\nu _\mu )`$ has been investigated in the two-body decay of a pion at rest:

$$\pi ^+\mu ^++\nu _\mu \mathrm{o}r\pi ^{}\mu ^{}+\overline{\nu }_\mu $$
(3)
Energy and momentum conservation result in a sharp muon momentum $`p(\mu )`$ from which the mass of the muon neutrino $`m(\nu _\mu )`$ would follow as:

$$m^2(\nu _\mu )=m^2(\pi )+m^2(\mu )2m(\pi )\sqrt{m^2(\mu )+p^2(\mu )}$$
(4)
Eq. (4) only holds, if the muon neutrino $`\nu _\mu `$ is a well-defined mass eigenstate, which does not apply in the case of neutrino mixing. Hence, if the muon momentum $`p(\mu )`$ for pion decay at rest could be measured with sufficient precision one would detect three different values $`p_\mathrm{i}^2(\mu )`$ with relative fraction $`|U_{\mu i}^2|`$ corresponding to the measurements of the corresponding squared mass values $`m^2(\nu _\mathrm{i})`$ contributing to the muon neutrino $`\nu _\mu `$. Up to now no direct neutrino mass measurements has discriminated different neutrino masses or has established a signal for any non-zero neutrino mass.

Therefore on the left part of eq. (4) a mean squared average of mass eigenstates of the muon neutrino is defined as

$$m^2(\nu _\mu )=\underset{i}{}|U_{\mu i}^2|m^2(\nu _\mathrm{i})$$
(5)
To deduce $`m^2(\nu _\mu )`$ from (5) three quantities have to be measured with very high precision: the muon mass $`m(\mu )=105.6583568(52)`$ MeV , the pion mass $`m(\pi )=139.570180(350)`$ bMeV , $`m(\pi )=139.570180(350)`$ MeV, and the muon momentum from pion decay at rest $`p(\mu )=29.791998(110)`$ MeV, which has been derived in a dedicated experiment at Paul Scherrer Institute (Zürich) . Putting these values into equation (4) one obtains :

$$m^2(\nu _\mu )=0.016\pm 0.023\mathrm{M}eV^2$$
(6)
from which an upper limit on the muon neutrino mass can be derived

$$m(\nu _\mu )<190\mathrm{k}eV(90\%C.L.)$$
(7)
### 2.2 $`m(\nu _\tau )`$

The most sensitive information on $`m(\nu _\tau )`$ comes from the investigation of $`\tau `$ pairs produced at electron-positron colliders decaying into multi pions. Due to the large mass of the $`\tau `$ , decays into 5 and 6 pions give the highest sensitivity because they restrict the available phase space of the $`\nu _\tau `$. However, the corresponding branching ratios are rather small.

The quantity looked at is the invariant mass of the multiple pions $`M_\pi `$. Although $`M_\pi `$ does not have a direct physical meaning, the mass of the tau neutrino $`m(\nu _\tau )`$ restricts $`M_\pi `$ due to energy and momentum conservation. In the rest frame of the decaying $`\tau `$ lepton, $`M_\pi ^2`$ is expressed by:

$`M_\pi ^2`$ $`=`$ $`({\displaystyle \underset{\mathrm{j}}{}}E_\mathrm{j}(\pi ),{\displaystyle \underset{\mathrm{j}}{}}\stackrel{}{p_\mathrm{j}}(\pi ))^2=(m(\tau )E(\nu _\tau ),\stackrel{}{p}(\nu _\tau ))^2`$ (8)
$``$ $`\left(m(\tau )m(\nu _\tau )\right)^2`$ (9)
The most sensitive investigation comes from the ALEPH experiment at LEP. Its two-dimensional analysis in the $`M_\pi `$$`_\mathrm{j}E_{\mathrm{j},lab}(\pi )`$ plane restricts $`m(\nu _\tau )`$ :

$$m(\nu _\tau )<18.2\mathrm{M}eV(95\%C.L.)$$
(10)
A further improvement based on data from B-factories can be expected with an estimated sensitivity limit of 3 MeV.

Using the most recent results on atmospheric and solar neutrino oscillation (see chapter 5) the neutrino mixing matrix $`U_{\alpha i}`$ and the squared mass differences $`\mathrm{\Delta }m_{\mathrm{ij}}^2`$ suggest that the averages $`m^2(\nu _\mu )`$ and $`m^2(\nu _\tau )`$ (compare eq. (5)) are rather close due to the strong $`\nu _\mu `$ \- $`\nu _\tau `$ mixing and the very small difference $`\mathrm{\Delta }m_{23}^2`$. Therefore $`m(\nu _\tau )`$ is already constrained by the limit on the muon neutrino mass (7).

### 2.3 $`m(\nu _\mathrm{e})`$

The mass of the electron neutrino is determined by the investigation of the electron energy spectrum ($`\beta `$ spectrum) of a nuclear $`\beta `$ decay . In a $`\beta ^{}`$ decay

$$(\mathrm{Z},\mathrm{A})(\mathrm{Z}+1,\mathrm{A})^++\mathrm{e}^{}+\overline{\nu }_\mathrm{e}$$
(11)
the available energy is shared between the $`\beta `$–electron and the electron antineutrino, because the recoiling nucleus practically receives no kinetic energy due to its much heavier mass. The phase space region of non-relativistic neutrinos, where the highest sensitivity to the neutrino mass is achieved, corresponds to the very upper end of the $`\beta `$ spectrum. To maximize this part, a $`\beta `$ emitter with a very low endpoint energy $`E_0`$ is required. This requirement is fulfilled by 187Re and tritium (T or 3H), which have the two lowest endpoint energies of $`E_0=2.6`$ keV and $`E_0=18.6`$ keV, respectively.

Although tritium has a higher endpoint energy as compared to 187Re its use has several advantages:

* Tritium decays by a super-allowed transition into its mirror nucleus 3He resulting in a half life of 12.3 years, compared to the primordial half life of the forbidden transition of 187Re of $`510^{10}`$ a. The short half life yields a high specific activity and minimizes the inelastic processes of $`\beta `$ electrons within the tritium source.
* Due to the super-allowed decay the transition matrix element does not depend on the electron energy: the $`\beta `$ spectrum is determined entirely by the available phase space.
* Tritium has the simplest atomic shell minimizing the necessary corrections due to the electronic final states or inelastic scattering in the $`\beta `$ source.

These arguments clearly favor tritium for standard setup, which consists of a $`\beta `$ source connected to a $`\beta `$ spectrometer (sometimes called “passive source” setup). The advantage of the lower 187Re endpoint energy can only be exploited if the $`\beta `$ source and the spectrometer are identical (sometimes called “active source” setup), which is realized in the case of cryogenic bolometers for instance.

For an allowed or super-allowed transition the electron energy spectrum is given by Fermi’s Golden Rule

$`{\displaystyle \frac{dN}{dE}}`$ $`=`$ $`{\displaystyle \frac{G_\mathrm{F}^2}{2\pi ^3\mathrm{}^7}}\mathrm{cos}^2\mathrm{\Theta }_\mathrm{C}|M|^2F(E,Z+1)p(E+m)\epsilon `$ (12)
$`\sqrt{\epsilon ^2m^2(\nu _\mathrm{e})}\mathrm{\Theta }(\epsilon m(\nu _\mathrm{e}))`$
$`=`$ $`AF(E,Z+1)p(E+m)\epsilon `$
$`\sqrt{\epsilon ^2m^2(\nu _\mathrm{e})}\mathrm{\Theta }(\epsilon m(\nu _\mathrm{e}))`$
with the Fermi coupling constant $`G_\mathrm{F}`$, the Cabbibo angle $`\mathrm{\Theta }_\mathrm{C}`$, the nuclear transition matrix element $`M`$ (tritium: $`|M|^2=5.55\mathrm{}^6`$ ), the Fermi function $`F`$ describing the final electromagnetic interaction of the emitted $`\beta `$ electron with the daughter nucleus $`F(E,Z+1)`$, the electron mass, momentum and kinetic energy $`m`$, $`p`$ and $`E`$, and the energy difference $`\epsilon =E_0E`$. The Fermi function is approximately given by

$$F(E,Z+1)=\frac{2\pi \eta }{1exp(2\pi \eta )}$$
(13)
with the Sommerfeld parameter $`\eta =\alpha (Z+1)/\beta `$. Equation (12) only holds for the decay of a bare nucleus. For the more realistic case of an atom or a molecule the possible excitation of the electron shell due to the sudden change of the nuclear charge by one unit has to be taken into account. The atom or molecule will end in a specific state of excitation energy Vj with a probability Wj. The corresponding excitation probabilities can be calculated in the sudden approximation from the overlap of the primary electron wave function $`\mathrm{\Psi }_0`$ with the wave functions of the daughter ion $`\mathrm{\Psi }_{\mathrm{f},j}`$

$$W{}_{\mathrm{j}}{}^{}=|\mathrm{\Psi }_0|\mathrm{\Psi }_{\mathrm{f},j}|^2$$
(14)
Equation (12) is thus modified into a sum of $`\beta `$ spectra of amplitude Wj with different endpoint energies $`E_{0,j}=E_0V_\mathrm{j}`$

$`{\displaystyle \frac{dN}{dE}}`$ $`=`$ $`AF(E,Z+1)p(E+m)`$ (15)
$`{\displaystyle \underset{\mathrm{j}}{}}W{}_{\mathrm{j}}{}^{}\epsilon {}_{\mathrm{j}}{}^{}\sqrt{\epsilon {}_{\mathrm{j}}{}^{}{}_{}{}^{2}m^2(\nu _\mathrm{e})}\mathrm{\Theta }(\epsilon {}_{\mathrm{j}}{}^{}m(\nu _\mathrm{e}))`$
The energy differences $`\epsilon _\mathrm{j}`$ are then defined as $`\epsilon {}_{\mathrm{j}}{}^{}=E_0V{}_{\mathrm{j}}{}^{}E`$.

In case of neutrino mixing the spectrum is a sum of the components of decays into mass eigenstates

$`{\displaystyle \frac{dN}{dE}}`$ $`=`$ $`AF(E,Z+1)p(E+m)`$ (16)
$`{\displaystyle \underset{\mathrm{j}}{}}W{}_{\mathrm{j}}{}^{}\epsilon {}_{\mathrm{j}}{}^{}({\displaystyle \underset{\mathrm{i}}{}}|U_{\mathrm{e}i}|^2\sqrt{\epsilon {}_{\mathrm{j}}{}^{}{}_{}{}^{2}m^2(\nu _\mathrm{i})}\mathrm{\Theta }(\epsilon {}_{\mathrm{j}}{}^{}m(\nu _\mathrm{i})))`$
When this spectrum is convoluted with an experimental resolution function which is much wider than the mass difference $`|m(\nu _\mathrm{i})m(\nu _\mathrm{j})|`$ (which has always been the case so far), the resulting spectrum can be analyzed in terms of a single mean squared electron neutrino mass

$$m^2(\nu _\mathrm{e})=\underset{\mathrm{i}}{}|U_{\mathrm{e}i}|^2m^2(\nu _\mathrm{i})$$
(17)
and eq. (15) applies again.

The square root term of equation (12) shows that the neutrino mass influences the $`\beta `$ spectrum only at the upper end below $`E_0`$ and its relative influence is vanishing as function of $`m^2(\nu _\mathrm{e})/\epsilon ^2`$ (see figure 1) leading far below the endpoint to a small constant offset proportional to $`m^2(\nu _\mathrm{e})`$.

Figure 1 defines the requirements of a direct neutrino mass experiments which investigates a $`\beta `$ spectrum : The task is to resolve the tiny change of the spectral shape due to the neutrino mass in the region just below the endpoint $`E_0`$ , where the count rate is going to vanish. Therefore, high energy resolution is required combined with large source strength and acceptance as well as low background rate.

## 3 Tritium $`\beta `$ Decay Experiments

### 3.1 Overview

The majority of the published direct laboratory results on $`m(\nu _\mathrm{e})`$ originates from the investigation of tritium $`\beta `$ decay, while one single result from 187Re has been reported at conferences 555There are also results from investigations of electron capture  and bound state $`\beta `$ decay , which are about 2 orders of magnitude less stringent on the neutrino mass.. In the early nineties tritium $`\beta `$ decay experiments yielded controversially discussed results: Figure 2 shows the final results of the experiments at Los Alamos National Laboratory and Zürich together with the results from other more recent measurements with magnetic spectrometers at University of Tokio, Lawrence Livermore National Laboratory and at Bejing. The sensitivity on the neutrino mass have improved a lot compared to previous experiments but the values for the observable $`m^2(\nu _\mathrm{e})`$ populated the unphysical negative $`m^2(\nu _\mathrm{e})`$ region. In the case of two experiments significantly negative mass results were obtained. In 1991 and 1994 two new experiments started data taking at Mainz and at Troitsk, which used a new type of electrostatic spectrometer, so-called MAC-E-Filters, which were again superior in energy resolution and luminosity with respect to the previous magnetic spectrometers. However, even their early data were confirming the large negative $`m^2(\nu _\mathrm{e})`$ values of the Los Alamos and Livermore experiments when being analyzed over the last 500 eV of the $`\beta `$ spectrum below the endpoint $`E_0`$. Also a new feature was observed. The large negative values of $`m^2(\nu _\mathrm{e})`$ disappeared when analyzing only small intervals below the endpoint $`E_0`$ (see also figure 6). This effect, which could only be investigated by the high-resolution MAC-E filters, pointed towards an underestimated or missing energy loss process, seemingly to be present in all experiments. The only common feature of the various experiment seemed to be the calculations of the excitation energies $`V_\mathrm{i}`$ of the daughter ions and their probabilities $`W_\mathrm{i}`$. Different theory groups checked these calculations in detail. The expansion was calculated to one order further and new interesting insight into this problem was obtained, but no significant changes were found (see  and references therein).

Then the Mainz group found the origin of the missing energy loss process for its experiment. The Mainz experiment uses as tritium source a film of molecular tritium quench-condensed onto aluminum or graphite substrates. Although the film was prepared as a homogenous thin film with flat surface, detailed studies showed that the film undergoes a temperature activated roughening transition into an inhomogeneous film by formation of microcrystals leading to unexpected large inelastic scattering probabilities.

The Troitsk experiment on the other hand uses a windowless gaseous molecular tritium source, similar to the Los Alamos apparatus. Here, the influence of large angle scattering of electrons magnetically trapped in the tritium source was not considered in the first analysis. The Troitsk group stated, that if this effect is taken into account it gives a correction large enough to make the negative values for $`m^2(\nu _\mathrm{e})`$ disappear.

It is very likely that also for the experiments at Los Alamos and Livermore other experimental effects caused the negative values for $`m^2(\nu _\mathrm{e})`$. Any missed or underestimated experimental correction leads to a negative value for $`m^2(\nu _\mathrm{e})`$. This can be understood by the following consideration: for $`\epsilon m(\nu _\mathrm{e})`$, eq. (12) can be expanded into

$$\frac{dN}{dE}\epsilon ^2m^2(\nu _\mathrm{e})/2$$
(18)
On the other hand the convolution of a $`\beta `$ spectrum (12) with a Gaussian of width $`\sigma `$ leads to

$$\frac{dN}{dE}\epsilon ^2+\sigma ^2$$
(19)
Therefore, in the presence of a missed experimental broadening with Gaussian width $`\sigma `$ one expects a shift of the result on $`m^2(\nu _\mathrm{e})`$ of

$$\mathrm{\Delta }m^2(\nu _\mathrm{e})2\sigma ^2$$
(20)
which gives rise to a negative value of $`m^2(\nu _\mathrm{e})`$.

### 3.2 MAC-E-Filter

The significant improvement in the $`\nu `$–mass sensitivity by the Troitsk and the Mainz experiments are due to MAC-E-Filters. This new type of spectrometer is based on early work by Kruit  and was later redeveloped for the application to the tritium $`\beta `$ decay at Mainz and Troitsk independently . The MAC-E-Filter combines high luminosity at low background and a high energy resolution. Both features are essential to measure the neutrino mass from the endpoint region of a $`\beta `$ decay spectrum. The acronym MAC-E-Filter stands for Magnetic Adiabatic Collimation followed by an Electrostatic Filter.

The main features of the MAC-E-Filter are illustrated in figure 3: two superconducting solenoids are producing a magnetic guiding field. The $`\beta `$ electrons, which are starting from the tritium source in the left solenoid into the forward hemisphere, are guided magnetically on a cyclotron motion along the magnetic field lines into the spectrometer, thus resulting in an accepted solid angle of nearly $`2\pi `$. On their way into the center of the spectrometer the magnetic field $`B`$ drops by several orders of magnitude. Therefore, the magnetic gradient force transforms most of the cyclotron energy $`E_{}`$ into longitudinal motion. This is illustrated in figure 3 at the bottom by a momentum vector. Due to the slowly varying magnetic field the momentum transforms adiabatically keeping the magnetic moment $`\mu `$ constant (equation given in non-relativistic approximation):

$$\mu =\frac{E_{}}{B}=const.$$
(21)
This transformation can be summarized as follows: the $`\beta `$ electrons, isotropically emitted at the source, are transformed into a broad beam of electrons flying almost parallel to the magnetic field lines.

This parallel beam of electrons is energetically analyzed by applying an electrostatic potential made up by a system of cylindrical electrodes. All electrons, which have enough energy to pass the electrostatic barrier are reaccelerated and collimated onto a detector, all other electrons are reflected. Therefore the spectrometer acts as an integrating high-energy pass filter. The relative sharpness of this filter is only given by the ratio of the minimum magnetic field $`B_{\mathrm{min}}`$ in the analyzing plane in the middle and the maximum magnetic field between $`\beta `$ electron source and spectrometer $`B_{\mathrm{max}}`$:

$$\frac{\mathrm{\Delta }E}{E}=\frac{B_{\mathrm{min}}}{B_{\mathrm{max}}}$$
(22)
By scanning the electrostatic retarding potential the $`\beta `$ spectrum can be measured.

Both experiments at Mainz and Troitsk are using similar MAC-E-Filters, which differ slightly in size: The diameter and length of the Mainz(Troitsk) spectrometer are 1 m (1.5 m) and 4 m (7 m). The major differences between the two setups are the tritium sources.

### 3.3 The Troitsk Neutrino Mass Experiment

The windowless gaseous tritium source of the Troitsk experiment  is essentially a tube of 5 cm diameter filled with $`\mathrm{T}_2`$ resulting in a column density of $`\rho d10^{17}`$ molecules$`/\mathrm{cm}^2`$. The source is connected to the ultrahigh vacuum of the spectrometer by a series a differential pumping stations.

From their first measurement in 1994 on the Troitsk group reported the observation of a small, but significant anomaly in their experimental spectrum starting a few eV below the $`\beta `$ endpoint $`E_0`$. This anomaly appears as a sharp step of the count rate . Since a MAC-E-Filter is integrating, this step should correspond to a narrow line in the primary spectrum with a relative intensity of about $`10^{10}`$ of the total decay rate. In 1998 the Troitsk group reported that the position of this line oscillates with a frequency of 0.5 years between 5 eV and 15 eV below $`E_0`$ . In 2000 the anomaly did not follow the 0.5 year periodicity anymore but still existed . In total the Troitsk experiment has taken 200 days of tritium data, in almost all of the runs this anomaly has been observed.

The reason for such an anomaly with these features is not clear. Detailed investigations at Troitsk are continuing. In addition synchronous measurements with the Mainz experiment were performed. In 2001 the Troitsk group improved the differential pumping between the gaseous tritium source and the spectrometer, lowered the electric field strength in a critical region and improved the vacuum. The first two runs of 2001 either gave no indication for an anomaly or only showed a small effect with 2.5 mHz amplitude if compared to the previous ones with amplitudes between 2.5 mHz and 13 mHz. These findings also support the assumption that the Troitsk anomaly is due to an still unknown experimental artefact .

Fitting a standard $`\beta `$ spectrum to the data the Troitsk group obtained significantly negative values of $`m^2(\nu _\mathrm{e})`$ of -10 to -20 $`\mathrm{e}V^2`$ (see filled circle in fig. 2). Describing the anomaly phenomenologically by adding a monoenergetic line, free in amplitude and position, to a standard $`\beta `$ spectrum results in values of $`m^2(\nu _\mathrm{e})`$ compatible to zero  (see also open circles in fig. 2). After this correction the average over all runs until 2001 amounts to :

$$m^2(\nu _\mathrm{e})=2.3\pm 2.5\pm 2.0\mathrm{e}V^2$$
(23)
which corresponds - under the assumption that the run-by-run correction by an additional line is correct - to an upper limit  of

$$m(\nu _\mathrm{e})2.2\mathrm{e}V(95\%\mathrm{C}.\mathrm{L}.)$$
(24)
### 3.4 The Mainz Neutrino Mass Experiment

Mainz uses a film of molecular tritium quench-condensed onto a graphite substrate (HOPG). The film has a diameter of 17 mm and a typical thickness of 40 nm, which is measured by laser ellipsometry. The problem of the roughening transitions mentioned above have been investigated by the Mainz group in cooperation with the condensed matter group of P. Leiderer at Konstanz/Germany using conversion electron spectroscopy and scattered light techniques on different hydrogen isotopes . The following results were obtained: The roughening transition follows an Arrhenius-type law. Thus it cannot be avoided but its speed can drastically be slowed down by using lower temperatures. A $`\mathrm{T}_2`$ film at 2 K has a time constant of order 10 y , i.e. much longer than a typical duration of a measurement.

In the years 1995-1997 the Mainz setup was upgraded by a new cryostat providing temperatures of the tritium film below 2 K to avoid the roughening transition. Also a new tilted pair of superconducting solenoids was installed (see figure 4). Consequently $`\beta `$ particles from the source are still guided magnetically into the spectrometer, whereas tritium molecules evaporating from the source are trapped on the bend of the LHe cold tube covered with graphite. This measure eliminated source correlated background and allowed to increase the source strength significantly. The upgrade of the Mainz setup was completed by the application of HF pulses on one of the electrodes in between measurements every 20 s, and a full automation of the apparatus and remote control. The former improvement lowers and stabilizes the background, the latter one allows long–term measurements.

Figure 5 shows the endpoint region of the Mainz 1998 and 1999 data  in comparison with the former Mainz 1994 data . An improvement of the signal-to-background ratio by a factor 10 as well as a significant enhancement of the statistical quality of the data is clearly visible. A fit with $`m^2(\nu _\mathrm{e})`$ fixed to zero perfectly fits the latter data set over these last 15 eV of the $`\beta `$ spectrum. This limits any persistent spectral anomaly in this range to an amplitude below $`10^3`$/s (as against a total flux of $`10^8`$/s entering the spectrometer). A spectral anomaly, like the fluctuating anomaly reported by the Troitsk group , on the other hand, reaches amplitudes up to $`10^2`$/s.

The main systematic uncertainties of the Mainz experiment are originating from the physics and the properties of the quench-condensed tritium film: the inelastic scattering of $`\beta `$ electrons within the tritium film, the excitation of neighbor molecules due to the $`\beta `$ decay, and the self-charging of the tritium film by radioactivity.

These systematic uncertainties were studied in detail by various investigations  and the knowledge of the corresponding corrections could be significantly improved.

Figure 6 shows the fit results of the combined Mainz 1998 and 1999 data on $`m^2(\nu _\mathrm{e})`$ as function of the lower limit of the fit interval. The monotonous trend towards negative values of $`m^2(\nu _\mathrm{e})`$ for larger fit intervals, as observed for the Mainz 1991 and 1994 data , has vanished. This shows that the dewetting of the $`\mathrm{T}_2`$ film from the graphite substrate  indeed was the reason for this behavior. Now this effect is safely suppressed at the much lower temperature of the $`\mathrm{T}_2`$ film. The data do not show any indication for other residual distortions. The energy interval below the endpoint, yielding the smallest combined statistical and systematical uncertainty on the neutrino mass is obtained, corresponds to the last 70 eV below the endpoint $`E_0`$ and gives

$$m^2(\nu _\mathrm{e})=1.6\pm 2.5\pm 2.1\mathrm{e}V^2$$
(25)
which is compatible with a zero neutrino mass. Considering its uncertainties, this value corresponds to an upper limit on the electron neutrino mass of :

$$m(\nu _\mathrm{e})2.2\mathrm{e}V(95\%\mathrm{C}.\mathrm{L}.)$$
(26)
Further data have been taken at Mainz in 2000, but suffering from background problems At the end of 2001 the Mainz group started another 3 month measurement campaign. A very careful maintenance and preparation of the whole setup was done. Especially all parts which needs refreshment from time to time were replaced. These were components related to the tritium source, the vacuum and the high voltage system (e.g. the graphite substrate, out-baking of all vacuum systems and re-activation of the non-evaporable getter pumps, the oil of the high voltage divider, and others). To check the various scanning hysteresis effects, different scanning methods were tested, including a 50 times slower scanning procedure. The summary of all these investigations is: the Mainz experiment has achieved its most stable operation ever in these two runs of 2001 – no scanning direction effect was observed anymore. The background rate was about 13 mHz, lower than before and much more stable without the need of repeating high voltage conditioning during a run (see fig. 5).

The neutrino squared masses obtained from the fit are very stable and compatible with zero within their uncertainties and the previous Mainz results (see figure 7). No indication of a Troitsk-like anomaly or any residual problem in the Mainz data were found. The result for the 2001 data of the last 70 eV of the $`\beta `$ spectrum below the endpoint (Elow=18.5 keV, see fig. 7) on $`m^2(\nu _\mathrm{e})`$ is :

$$m^2(\nu _\mathrm{e})=+0.1\pm 4.2\pm 2.0\mathrm{e}V^2$$
(27)
Combining this value with the one obtained from the data sets Q5–Q8 from 1998 and 1999 over the last 70 eV (25)  gives

$$m^2(\nu _\mathrm{e})=1.2\pm 2.2\pm 2.1\mathrm{e}V^2$$
(28)
which corresponds again to an upper limit  of

$$m(\nu _\mathrm{e})<2.2\mathrm{e}V(95\%\mathrm{C}.\mathrm{L}.)$$
(29)
The inclusion of the high-quality data from 2001 improves the Mainz sensitivity only marginally, showing that the Mainz experiment has reached its sensitivity limit.

In spring 2002 the Mainz group has installed a new electrode system to check new ideas which were deloped for the KATRIN experiment to avoid background and to remove trapped particles . First measurements showed that the new ideas indeed are reducing the background rate by a factor 3.

Indications for a “Troitsk-like” anomaly at Mainz were observed only once in summer 1998. This single coincidence did not appear in previous and later runs . Of special interest are the Mainz 2000 runs, although they were done under less favored conditions. In particular, part of the data have been taken in parallel with Troitsk. Whereas Troitsk again has indications for the anomaly  at Mainz no indication for a Troitsk-like anomaly was found . In summary, the Mainz data clearly state that the anomaly observed at Troitsk is caused by some experimental artefact .

## 4 Future direct neutrino mass searches

As briefly discussed in section 1 the compelling evidence for non-zero neutrino masses from atmospheric and solar neutrino experiments provides squared neutrino mass differences but no absolute neutrino masses. These findings clearly demand for the determination of the absolute neutrino mass scale as one of the most important next steps in neutrino physics since the absolute neutrino mass has strong consequences for astrophysics and cosmology as well as for nuclear and particle physics as discussed in section 1.

There exist different ideas and approaches to determine the absolute neutrino mass with sub-eV sensitivity, which will be briefly discussed in the following:

* Time-of-flight of supernova neutrinos
The spread of the arrival time of supernova neutrinos on earth depends on neutrino energy and mass, thus allowing to extract the neutrino mass by measuring arrival time and energy. A supernova, exploding within our galaxy, would give hundreds to thousand of neutrino events in the current underground neutrino detectors. Although this would exceed the statistics of the only supernova SN1987a observed so far by two orders of magnitude, the systematic uncertainty connected with the not precisely known neutrino emission time spectrum does not allow a sub-eV sensitivity on the neutrino mass.
* Large scale structure
The observation of the structure in the universe at different scales and the angular distribution of the fluctuations of the cosmic microwave background radiation allows to set constraints on the hot dark matter content of the early universe . The relic neutrino density connects this hot dark matter density to the neutrino mass. Although the expected sensitivity on the neutrino mass for the near future is in the sub-eV range, the results derived this way are never model independent. Reversing the arguements, there are strong degeneracies between the different parameters and its therefore very helpful to use information from laboratory neutrino mass experiments to determine the other astrophysics parameters more precisely.
* Neutrinoless double $`\beta `$ decay
The neutrinoless double $`\beta `$ decay is sensitive to the so-called “effective” neutrino mass

$$m_{\mathrm{e}e}=|\underset{\mathrm{i}}{}U_{\mathrm{e}i}^2m(\nu _\mathrm{i})|,$$
(30)
which is a coherent sum over all mass eigenstates contributing to the electron neutrino with fraction $`U_{\mathrm{e}i}`$. The determination of $`m_{\mathrm{e}e}`$ from the measurement of the neutrinoless double $`\beta `$ decay rate is complementary to the direct determination of the mass of the electron neutrino since $`m_{\mathrm{e}e}`$ and $`m(\nu _\mathrm{e})`$ can differ by the following reasons:

1. Double $`\beta `$ decay requires the neutrino to be a Majorana particle.
    2. The values $`U_{\mathrm{e}i}^2`$ in eq. (30) can have complex phases, which could lead to a partial cancellation of the different terms of the sum. Especially that the recent solar neutrino data point to large mixing opens this possibility .
    3. The uncertainty of the nuclear matrix elements of neutrinoless double $`\beta `$ decay still contributes to the uncertainty of $`m_{\mathrm{e}e}`$ by about a factor of 2 .
    4. Non Standard Model processes, others than the exchange of a Majorana neutrino, could enhance the observed neutrinoless double $`\beta `$ decay rate without changing $`m_{\mathrm{e}e}`$.

The lowest limit of $`m_{\mathrm{e}e}<0.35`$ eV is coming from the Heidelberg-Moscow experiment using an array of semiconductor detectors of enriched 76Ge . Very recently part of the collaboration interpreted the data as a signal for neutrinoless double $`\beta `$ decay , a claim which raised discussions within the community. Not only to check this, but also for the reasons given above double $`\beta `$ decay experiments with much enhanced sensitivity are clearly needed. The proposed double $`\beta `$ decay experiments of the next generation aim for a sensitivity on $`m_{\mathrm{e}e}`$ in the range of 0.1 eV and below .
* Rhenium cryogenic bolometer experiments
Due to the complicated electronic structure of 187Re and its $`\beta `$ decay (compare to section 2.3) the advantage of the 7 times lower endpoint energy $`E_0`$ of 187Re with respect to tritium can only be exploited if the $`\beta `$ spectrometer measures the entire released energy, except that of the neutrino. This situation can be realized by using a cryogenic bolometer as the $`\beta `$ spectrometer, which at the same time contains the $`\beta `$ emitter 187Re (see figure 8).

One disadvantage connected to this method is the fact that one measures always the entire $`\beta `$ spectrum. Even for the case of the very low endpoint energy of 187Re, the relative fraction of events in the last eV below $`E_0`$ is of order $`10^{10}`$ only (compare to figure 1). Considering the long time constant of signal of a cryogenic bolometer (typically several hundred $`\mu `$s) only large arrays of cryogenic bolometers can deliver the signal rate needed.

Up to now two groups are working on 187Re $`\beta `$ decay experiments at Milano  and Genoa . Although cryogenic bolometers with an energy resolution of 5 eV have been produced with other absorbers, this has yet not been reached for rhenium. The two groups are using different ways to produce the crystals. The MANU2 experiment at Genoa succeeded in preparing crystals from metallic rhenium. The group has reported a limit on $`m(\nu _\mathrm{e})`$ of 26 $`\mathrm{e}V`$ . The Genoa group understands their measured spectra well and has seen for the first time the oscillation pattern of the $`\beta `$ environmental fine structure , which describes the interference between the coherent scattered outgoing electron wave function with the crystal similar to the X-ray environmental fine structure (XEFS) known for X-ray absorption edges. The Genoa group expects in the near future a sensitivity on $`m(\nu _\mathrm{e})`$ of 10 $`\mathrm{e}V`$. A significant further improvement could be obtained by improving the energy resolution of the crystals by using new super-conducting transition thermometers. The MiBeta experiment of Milano uses AgReO4 crystals with a typical energy resolution of 35 eV. There is no result on $`m(\nu _\mathrm{e})`$ yet, the expected sensitivity is similar to the one for the Genoa experiment.

To further improve the statistical accuracy the operation of large arrays of micro-calorimeters with better resolution is required. New techniques are explored to enable these improvements. The expected sensitivity on the neutrino mass in the future is in the eV region .
* Next generation tritium $`\beta `$ decay experiment
Summarizing the discussion above clearly means that one or more next generation double $`\beta `$ decay experiments have to be performed due to their very low sensitivity. But considering the complementariness of neutrinoless double $`\beta `$ decay and the direct neutrino mass determination it is also clear that a next generation direct mass search has to be done. None of the alternative direct methods discussed above is able to provide a sub-eV sensitivity in the next decade. Therefore, it is straightforward to explore which sensitivity could be achieved by investigating the tritium $`\beta `$ decay spectrum near its endpoint with the very successful MAC-E-Filter as spectrometer.

Discussions between groups from Mainz, Karlsruhe and Troitsk led to the proposal for a next generation tritium $`\beta `$ decay experiment to be built at Forschungszentrum Karlsruhe/Germany. This idea was strongly supported by the community at the international workshop on “Neutrino masses in the sub-eV range” at Bad Liebenzell/Germany in early 2001 . A strong collaboration including nearly the complete worldwide expertise on tritium $`\beta `$ decay neutrino mass experiments has come together and a letter of intent for the KATRIN experiment (KArlsruhe TRItium Neutrino experiment) has been published . First funding for KATRIN has been obtained. This experiment will be described in the next section in more detail.

## 5 The KATRIN experiment

Figure 9) shows a schematic view of the proposed experimental configuration.

The windowless gaseous tritium source (WGTS) allows the measurement of the endpoint region of the tritium $`\beta `$ decay and consequently the determination of the neutrino mass with a minimum of systematic uncertainties from the tritium source. The WGTS will consist of a 10 m long cylindrical tube of 70 mm diameter, filled with molecular tritium gas of high isotopic purity ($`>95`$ %). The tritium gas will be injected by a capillary at the middle of the tube and then diffuses over a length of 5 m to both end faces of the tube, resulting in a source column density of about $`(\rho d)510^{17}`$ molecules/cm2. This gives a near to maximum count rate at the endpoint666One has to keep in mind, that increasing the source column density further almost does not increase the count rate near the endpoint, since inelastic scattering of electrons on tritium molecules has a threshold of about 12 eV . Any inelastically scattered electron cannot contribute to the count rate in the interesting region of the last 12 eV below the $`\beta `$ endpoint $`E_0`$.. With these values the count rate is increased by factor 40 with respect to the Troitsk experiment.

A quench condensed tritium source (QCTS) following the source concept of the Mainz experiment is considered as a second alternative source, which has complementary systematics.

The electron transport system adiabatically guides $`\beta `$ decay electrons from the tritium sources to the spectrometer while at the same time eliminating any tritium flow towards the spectrometer, which has to be kept practically free of tritium for background reasons. The first part of the transport system consists of a differential pumping section with a tritium reduction of a factor $`10^9`$, the second part of a liquid helium cold cryo-trapping section.

Between the tritium sources and the main spectrometer a pre-spectrometer of MAC-E-Filter type will be inserted, acting as an energy pre-filter to reject all $`\beta `$ electrons except the ones in the region of interest close to the endpoint $`E_0`$. This minimizes the chances of causing background by ionization of residual gas in the main spectrometer. As the designs of the pre- and main spectrometer will be similar, the former will act as a test facility for the larger main spectrometer . The design and construction of the pre-spectrometer has already started. Especially important will be the following tests of : a) the technique to achieve an XUHV of below $`10^{11}`$ mbar, b) the concept of using the vacuum vessel itself as main electrode at high potential, and c) the new electromagnetic design concepts to reduce background.

A key component of the new experiment will be the large electrostatic main spectrometer with a diameter of 7 m and an overall length of about 20 m. This high resolution MAC-E-Filter will allow to scan the tritium $`\beta `$ decay endpoint with increased luminosity at a resolution of 1 eV, which is a factor of 4 better than present MAC-E-Filters at Mainz and Troitsk. The 100 times larger analyzing plane with respect to the Mainz experiment allows the remaining factor 25 to be utilized to increase the source cross section and, correspondingly, the signal rate.

The detector requires high efficiency for electrons at $`E_0=18.6`$ keV and low $`\gamma `$ background. A high energy resolution of $`\mathrm{\Delta }E<600`$ eV for 18.6 keV electrons should suppress background events at different energies. The present concept of the detector is based on a large array of about 1000 silicon drift detectors surrounded by low-level passive shielding and an active veto counter to reduce background. For a later stage of the experiment a segmented bolometer is considered as a possible detector.

The difficulty in measuring a sub eV neutrino mass concerning the statistics, namely the smallness of the interesting region below the endpoint, turns into an advantage with respect to the systematic uncertainties due to energy thresholds for inelastic processes. First simulations with conservative assumptions result in a total 1 $`\sigma `$ uncertainty of $`\mathrm{\Delta }m^2(\nu _\mathrm{e})=0.08`$ eV2 for a measurement time of three years. If no finite neutrino mass is observed, the three year measurement leads to an upper limit of the mass of 0.35 eV at 90% confidence. This sensitivity improves the existing limits by almost one order of magnitude and also demonstrates the discovery potential of KATRIN for an electron neutrino mass in the sub-eV range.

## 6 Summary

Although neutrino oscillation experiments show that neutrinos have non-zero masses, the direct neutrino mass experiments have obtained only upper limits yet. The current tritium $`\beta `$ decay experiments at Mainz and Troitsk are reaching their sensitivity limits. The Mainz upper limit on $`m(\nu _\mathrm{e})`$ is 2.2 eV at 95 % C.L. The Troitsk group gives the same limit under the assumption that an anomalous excess count rate near the endpoint is described correctly. The synchronous measurements at Mainz and Troitsk show that this “Troitsk anomaly” is an experimental artefact.

A mass determination with sub-eV sensitivity is clearly needed to distinguish between hierarchical and degenerate neutrino mass models and to clarify the role of neutrinos in the early universe. The search for the neutrinoless double $`\beta `$ decay is one very important approach. Complementary and equally important is a next generation direct neutrino mass experiment. Discussing the different options shows that this experiment has to be a large tritium $`\beta `$ decay experiment using a MAC-E-Filter. Such an experiment is being prepared by the KATRIN collaboration.

## Acknowledgments

The author would like to thank the various collaborations for providing him kindly the presented informations and results. Especially acknowledged are the valuable remarks by J. Bonn, G. Drexlin, Ch. Kraus and E. Otten. The work of the Mainz and the KATRIN experiments connected to the author is supported by the German Bundesministerium für Bildung und Forschung under contracts 06MZ866I/5 and 05CK2PD1/5.


format: text

text:
# Quantitative evolution of global strings from the Lagrangian viewpoint

## I Introduction

The idea of spontaneous symmetry breaking in high energy physics has profound implications on the vacuum structure of the Universe and our Universe has presumably experienced various phase transitions, which are thermal Kibble  or nonthermal KVY . Their consequences may be traced by topological defects that may have been produced with them. Indeed, many implications of topological defects for cosmology have been investigated VS . Furthermore, recently, defect formations have also been discussed in the context of phase transitions which occur in the laboratory. For example, defect formations in 3He He3  and 4He He4  are studied in detail. Thus the cosmological scenario of defect formation can be tested by experiments in the laboratory Zurek .

Among several types of topological defects, strings hold a unique position in cosmology, because they do not overclose the Universe, unlike magnetic monopoles or domain walls, settling down to a scaling solution in which the typical scale of the network grows in proportion to the horizon scale Kibble ; Kibble2 . The key mechanism to achieve a scaling behavior is intercommutation of infinite strings to dissipate their energy by producing closed loops, which subsequently decay to radiation of relativistic particles or gravitational waves depending on their property. While one may understand the qualitative nature of the scaling solution analytically, more quantitative features such as the number of long strings per horizon volume or the size spectrum of loops produced cannot be obtained unless full numerical analyses are performed.

Although there exist two types of strings, namely, local and global strings depending on the nature of the symmetry breaking, only the former have been investigated extensively for a long time as far as numerical analysis is concerned. By numerically solving the equation of motion of string segments derived from the Nambu-Goto action NG , several groups have confirmed the scaling behavior AT ; AT2 ; BB ; AS , and estimated the scaling parameter $`\xi `$ as $`\xi 10`$ in the radiation dominated universe AT2 ; BB ; AS . Here, $`\xi `$ is defined as

$$\xi =\rho _{\mathrm{}}t^2/\mu $$
(1)
where $`\rho _{\mathrm{}}`$ is the energy density of long strings and $`\mu `$ is the string tension per unit length. Thanks to this feature, local strings generate density fluctuations with a scale-invariant spectrum and their cosmological consequences were investigated extensively some time ago. Recent observations of the cosmic microwave background anisotropy cmb , however, disfavor the cosmic-string scenario of structure formation, and the motivations to investigate local strings as a source of primordial density fluctuations have somewhat diminished, although a hybrid model of structure formation may still be viable, where primordial fluctuations are comprised of adiabatic fluctuations induced by inflation and isocurvature perturbations by topological defects TDCMB .111Note that topological defects can be compatible with cosmic inflation KVY . Note also that local strings may still be important in that they may emit massive particles as sources of ultrahigh energy cosmic rays VHS .

On the other hand, global strings are much better motivated in the context of axion cosmology. They are formed as a consequence of the breaking of the Peccei-Quinn U(1) symmetry PQ ; VE , which was introduced to solve the strong $`CP`$ problem in quantum chromodynamics. These global strings radiate axions as associated Nambu-Goldstone (NG) bosons Davis , which are one of the most promising candidates for cold dark matter. Despite their importance, the cosmological evolution of global strings has been less studied and the results of local strings have often been borrowed even though there is a decisive difference between them. For local strings, the gradient energy of scalar fields is canceled out by gauge fields far from the core. The string core is well localized and the false vacuum energy of the core dominates the system. Hence the Nambu-Goto action is suitable as an effective action in order to study the cosmological evolution of local strings except at crossing NG . On the other hand, for global strings, there are no gauge fields to cancel the gradient energy of the NG scalar field, which dominates over the false vacuum energy of the core. The effective action appropriate for global strings is not the Nambu-Goto action but the Kalb-Ramond action, which incorporates NG bosons and their couplings with the core KR . Also, due to the gradient energy of the NG scalar field, a long-range force works between global strings. Thus, the behavior of the two types of string is expected to be different and it is nontrivial whether global strings relax to a scaling regime. In fact, for example, the behavior of global monopoles is quite different from that of local monopoles due to long-range forces. While the former may be useful in cosmology monopole ; Yamaguchi , the latter causes a disaster Preskill  unless diluted by inflation. Furthermore, it has been shown in the literature that local and global strings behave differently in two dimensional space YB ; 2dim .

There have been several attempts to investigate the cosmological evolution of the global string network by use of the Kalb-Ramond action BS ; sharp . However, the Kalb-Ramond action is too complicated to be dealt with numerically. It has difficulty because of logarithmic divergence due to the self-energy of the string. In such a situation the authors and Kawasaki made the first numerical investigations of the cosmological evolution of global strings without resort to the Kalb-Ramond action. We instead solved the equations of motion for scalar fields forming strings in three dimensional Eulerian meshes YKY . We found that the global string network would also go into a scaling regime but the scaling parameter $`\xi `$ was found to be of the order of unity YKY ; YYK , which is significantly smaller than the case of local strings ($`\xi 10`$).

Recently, however, this quantitative difference was questioned and it was claimed that the smallness of the scaling parameter of global strings might be a fake due to the small dynamic range of the numerical simulations MS ; MSM . The authors of MS ; MSM  claimed that in our simulations global strings lost their energy by excessive direct emission of NG bosons. Based on the speculation that such direct emission of NG bosons from long strings would be negligible on cosmological scales, they reached a conjecture that both types of strings behave quantitatively in the same way in the cosmological context, namely, $`\xi 10`$, with the only difference between them being the energy loss mechanism of closed loops.

In order to clarify the evolution of global strings on cosmological scales, we should first study which is the dominant energy loss mechanism of long strings in numerical simulations, loop production or direct emission of NG bosons. We note that the master equation for the energy density of long strings, $`\rho _{\mathrm{}}`$, can be expressed as

$$\frac{d\rho _{\mathrm{}}}{dt}=2H(1+v^2)\rho _{\mathrm{}}\mathrm{\Gamma }_{\mathrm{loop}}\rho _{\mathrm{}}\mathrm{\Gamma }_{\mathrm{NG}}\rho _{\mathrm{}},$$
(2)
where the second and the third terms on the right-hand side represent energy loss due to loop formation and direct emission of NG bosons or axions, respectively, and $`v^2`$ denotes the average square velocity of string segments. For our purpose, we need to calculate both the loop production rate $`\mathrm{\Gamma }_{\mathrm{loop}}`$ and the emission rate of NG bosons $`\mathrm{\Gamma }_{\mathrm{NG}}`$, and compare their magnitudes in simulations.

If the system relaxes to the scaling regime, which will be confirmed to be the case shortly, the string network is described by the so-called one-scale model222Although improvements of the simplest one-scale model have been proposed by taking into account effects of small-scale structures ACK  and time evolution of the velocity MS ; MS2 , the original one-scale model is sufficient here because small-scale structures are smeared in the case of global strings and the velocity remains constant in the scaling regime, as will be shown later.  with the characteristic scale $`L\sqrt{\mu /\rho _{\mathrm{}}}`$, which grows with the horizon scale $`Lt`$. If we introduce the loop production coefficient $`c`$ and the emission coefficient $`\kappa `$ of the NG bosons by

$`\mathrm{\Gamma }_{\mathrm{loop}}\rho _{\mathrm{}}=c{\displaystyle \frac{\rho _{\mathrm{}}}{L}},\mathrm{\Gamma }_{\mathrm{NG}}\rho _{\mathrm{}}=\kappa {\displaystyle \frac{\rho _{\mathrm{}}}{L}},`$ (3)
these parameters remain constant in the one-scale picture and are related to the scaling parameter $`\xi `$ as

$$\xi =\left(\frac{1v^2}{c+\kappa }\right)^2.$$
(4)
Indeed, if $`\kappa `$ incorrectly turned out to be much larger than $`c`$, we would find a smaller value of $`\xi `$ than it should really be. Hence, it is essential to evaluate the parameters $`c`$ and $`\kappa `$ in the scaling regime.

In our previous simulations YKY , however, it was impossible to calculate these quantities for two reasons. First, in the previous work a lattice was identified as part of a string based on the value of the potential energy there, which had the microscopic problem that we occasionally found disconnected string pieces, although the overall features were traced reasonably well. Second, it was impossible to monitor intercommutation of strings for lack of dynamical information, namely, the velocity of each string segment. These problems are overcome by our new identification scheme and Lagrangian analysis of the evolution of strings YY .

The rest of the paper is organized as follows. In Sec. II we present the method of our new procedure to follow the Lagrangian evolution of global strings in Eulerian simulations, that is, new methods of identification of strings, measurement of string velocity, and estimation of the intercommutation rate. In Sec. III the results are described and applied to cosmological situations. Finally, Sec. IV is devoted to a discussion and conclusion.

## II numerical simulations

### II.1 Formulation

We consider the following Lagrangian density for two-component real scalar fields $`\varphi _a(x)`$ ($`a=1,2`$) which can accommodate global strings:

$$[\varphi _a]=\frac{1}{2}g_{\mu \nu }(^\mu \varphi _a)(^\nu \varphi _a)V[\varphi _a,T],$$
(5)
in the spatially flat Robertson-Walker spacetime,

$$ds^2=g_{\mu \nu }dx^\mu dx^\nu =dt^2R^2(t)d𝒙^2,$$
(6)
with $`R(t)`$ being the scale factor. We adopt the following potential at finite temperature $`T`$:

$`V[\varphi _a,T]`$ $`=`$ $`{\displaystyle \frac{\lambda }{4}}(\varphi ^2\sigma ^2)^2+{\displaystyle \frac{\lambda }{6}}T^2\varphi ^2`$ (7)
$`=`$ $`{\displaystyle \frac{\lambda }{4}}(\varphi ^2\eta ^2)^2+{\displaystyle \frac{\lambda }{4}}(\sigma ^4\eta ^4).`$ (8)
Here $`\varphi ^2\varphi _1^2+\varphi _2^2`$ and $`\eta ^2\sigma ^2T^2/3=\sigma ^2(1T^2/T_c^2)`$ with $`T_c\sqrt{3}\sigma `$ being the critical temperature. Below the critical temperature, global U(1) symmetry is broken to form global strings.

The equations of motion for the scalar fields are given by

$$\ddot{\varphi _a}(x)+3H(t)\dot{\varphi _a}(x)\frac{1}{R(t)^2}^2\varphi _a(x)+\frac{V}{\varphi _a}=0,$$
(9)
where an overdot represents the time derivative. The discretization of these differential equations is given in Appendix A. Our numerical calculations are based on the staggered leapfrog method with second order accuracy both in time and in space. In the radiation dominated Universe, the Hubble parameter $`H(t)=\dot{R}(t)/R(t)`$ and cosmic time $`t`$ are given by

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

$$\zeta \frac{ϵ}{\sigma }=\left(\frac{45M_{Pl}^2}{16\pi ^3g_{}\sigma ^2}\right)^{1/2},$$
(11)
and take $`\zeta =10`$ and the self-coupling $`\lambda =0.08`$ for definiteness, but these particular choices do not affect the results. We start the numerical simulation at the temperature $`T_i=2T_c`$ corresponding to $`t_i=t_c/4`$ and adopt as an initial condition the thermal equilibrium state with a mass equal to the inverse curvature of the potential at that time. In this state $`\varphi _a`$ and $`\dot{\varphi }_a`$ are Gaussian distributed with the correlation functions

$`\beta |\varphi _a(𝒙)\varphi _b(𝒚)|\beta _{\mathrm{equal}\mathrm{time}}`$ $`=`$ $`{\displaystyle \frac{d^3k}{(2\pi )^3}\frac{1}{2\sqrt{𝒌^2+m^2}}\mathrm{coth}\frac{\beta \sqrt{𝒌^2+m^2}}{2}e^{i𝒌(𝒙𝒚)}\delta _{ab}},`$ (12)
$`\beta |\dot{\varphi }_a(𝒙)\dot{\varphi }_b(𝒚)|\beta _{\mathrm{equal}\mathrm{time}}`$ $`=`$ $`{\displaystyle \frac{d^3k}{(2\pi )^3}\frac{\sqrt{𝒌^2+m^2}}{2}\mathrm{coth}\frac{\beta \sqrt{𝒌^2+m^2}}{2}e^{i𝒌(𝒙𝒚)}\delta _{ab}},`$ (13)
where $`m^2=V^{\prime \prime }[\varphi _a,T_i]`$ and $`\beta =T_i^1`$. $`\varphi _a(𝒙)`$ and $`\dot{\varphi }_a(𝒚)`$ are uncorrelated for $`𝒙𝒚`$. We generate the initial configuration in the momentum space, where the scalar fields are uncorrelated. Then they are transformed into the position space by fast Fourier transformation.

We perform simulations in five different sets of lattice sizes and spacings as shown in Table 1 to investigate their effects on the results. In all cases, the time step is taken as $`\mathrm{\Delta }t=0.01t_i`$, and the periodic boundary condition is adopted. In the case (a), the box size is nearly equal to the horizon length $`H^1`$ and the lattice spacing to the typical core width $`\delta 1.0/(\sqrt{2\lambda }\sigma )`$ of a string at the final time $`t_f=200t_i`$. The other cases have equal or larger simulation volumes.

In this type of numerical calculation, it is a nontrivial task to identify string cores from the data of scalar fields because a point with $`\varphi _a=0`$ corresponding to a string core is not necessarily situated at a lattice point. So, in the next subsection, we give our new method of identification of string cores.

### II.2 Identification method of strings

In previous work YKY ; YYK , whether a lattice point was a part of string core was judged based on the potential energy density there. That is, a lattice point was identified to be at the core of a string if the potential energy density there was found to be larger than that of a static cylindrically symmetric solution of a global string at $`r=\mathrm{\Delta }x_{\mathrm{phys}}/\sqrt{2}`$ off center, where $`\mathrm{\Delta }x_{\mathrm{phys}}`$ is the physical length of the lattice spacing at each time.

Although this method worked fairly well to evaluate the scaling parameter, it was inadequate to fully identify strings particularly for small loops which do not resemble the static cylindrically symmetric solution due to the curvature. As a result we occasionally found disconnected string segments that should not exist in reality. Furthermore, it was impossible to find a more correct position of the string core in a box beyond the lattice spacing, which is extremely important to evaluate the length and velocity of the string correctly.

Here we develop a new method of string identification based on the fact that strings lie on the intersection of two surfaces $`\varphi _1(𝒙,t)=0`$ and $`\varphi _2(𝒙,t)=0`$, so that if a string penetrates a (sufficiently small) plaquette, $`\varphi _a`$ has a different sign at one or two corners of the plaquette from the rest for each $`a`$. First, we classify the relative phase of the scalar fields into three groups as shown in Fig. 1 (left), that is,

$`\theta \mathrm{arccos}{\displaystyle \frac{\varphi _2}{\sqrt{\varphi _1^2+\varphi _2^2}}}+2\pi \left[1\mathrm{\Theta }(\varphi _2)\right],`$ (14)
$`(\mathrm{i})`$ $`0\theta <{\displaystyle \frac{\pi }{2}},`$
$`(\mathrm{ii})`$ $`{\displaystyle \frac{\pi }{2}}\theta <{\displaystyle \frac{3\pi }{2}},`$ (15)
$`(\mathrm{iii})`$ $`{\displaystyle \frac{3\pi }{2}}\theta <2\pi ,`$
where the arccosine should take the principal value and $`\mathrm{\Theta }(\varphi _2)`$ is the step function.

Then, we judge whether a string penetrates a plaquette by monitoring the phase rotation around it just as in the Vachaspati-Vilenkin algorithm VV . If we assigned an equal range $`2\pi /3`$ of the relative phase to all regions (i)-(iii) as in the original Vachaspati-Vilenkin algorithm and judged the presence of a string from the phase rotation, we would occasionally identify a plaquette as containing a string even if $`\varphi _1`$ or $`\varphi _2`$ takes the same sign at its four corners. This is the main reason why all regions (i)-(iii) do not have equal ranges of the relative phase $`\theta `$ in our scheme, in which this case is avoided and $`\varphi _a`$ takes a different sign for each $`a`$ at one or two corners of each plaquette along which a circular phase rotation is observed. Then, using the values of $`\varphi _a`$ at its four corners we can draw two lines corresponding to $`\varphi _1=0`$ and $`\varphi _2=0`$ and their intersection is identified with the point where a string penetrates the plaquette, as shown in Fig. 2.

In fact, the intersection point could be found outside the plaquette. In this case, the nearest point on the edge of the plaquette from the intersection is identified as the position of the string as shown in Fig. 3. In our simulations we encountered such a case very rarely and it did not cause any serious problems.

By joining these points, the strings are completely connected and a more accurate total length of the global string can be obtained. One may wonder if our biased classification of the relative phase may cause some artificial effects on the result of numerical simulations. So we have tried another classification of the relative phase as shown in Fig. 1 (right), that is,

$`(\mathrm{i})`$ $`0\theta <{\displaystyle \frac{\pi }{2}},`$
$`(\mathrm{ii})`$ $`{\displaystyle \frac{\pi }{2}}\theta \pi ,`$ (16)
$`(\mathrm{iii})`$ $`\pi \theta <2\pi ,`$
and confirmed that the results do not depend on the classification scheme for the relative phase.

At each time step, we can evaluate the total length of the global string by the method shown above. Since strings typically move at a speed close to unity, we should multiply the length of each string segment by $`\mu =\gamma \mu _s`$ to calculate the total energy density of strings. Here $`\gamma `$ is the Lorentz factor and $`\mu _s2\pi \eta ^2\mathrm{ln}\left[t/(\delta \xi ^{1/2})\right]`$ is the line density of a static string MSM . So it is important to establish a method to calculate the velocity and Lorentz factor of the string segments, which will be presented shortly. We can then evaluate the scaling parameter $`\xi `$ from Eq. (1).

### II.3 Velocity of strings

Here we describe our method to evaluate the velocity of strings, which is a nontrivial task in Eulerian calculations of scalar field configurations.333The measurement of the velocity was also attempted in Ref. MSM . In their method, the velocity is obtained by comparing the positions of strings at different times. Although some efforts have been made to reduce errors, their method may cause some systematic errors because the time separation cannot be shorter than that sufficient to ensure a string can move to another lattice. On the other hand, our method is much superior; it is based on first principles and the velocity of each string segment is given by a local quantity.

First we expand the scalar fields $`\varphi _a(𝒙,t_0+\delta t)`$ around $`\varphi _a(𝒙_0,t_0)`$ up to first order,

$`\varphi _a(𝒙,t_0+\delta t)\varphi _a(𝒙_0,t_0)+\varphi _a(𝒙_0,t_0)(𝒙𝒙_0)+\dot{\varphi }_a(𝒙_0,t_0)\delta t(a=1,2).`$ (17)
Suppose that a string core exists at a point $`𝒙_0`$ at time $`t_0`$ and moves to a point $`𝒙`$ at $`t=t_0+\delta t`$, that is, $`\varphi _a(𝒙_0,t_0)=0`$ and $`\varphi _a(𝒙,t_0+\delta t)=0`$ for each $`a`$. Then, from Eq. (17) we find that the loci of the string core at $`t=t_0+\delta t`$ lie on the intersection of two planes,

$$\varphi _a(𝒙_0,t_0)(𝒙𝒙_0)+\dot{\varphi }_a(𝒙_0,t_0)\delta t=0,$$
(18)
with $`a=1,2`$. Since motion tangential to a string is a gauge mode, we should evaluate the velocity normal to it. Suppose that the line normal to the string segment at $`(𝒙_0,t_0)`$ reaches across the above intersection line at $`𝒙=𝒙_l(𝒙_0,t_0,\delta t)`$. Then we can easily obtain the velocity of this string segment as

$`𝒗(𝒙_0,t_0)`$ $`=`$ $`\underset{\delta t0}{lim}{\displaystyle \frac{𝒙_l(𝒙_0,t_0,\delta t)𝒙_0}{\delta t}}`$ (19)
$`=`$ $`{\displaystyle \frac{\dot{\varphi }_1\varphi _2\dot{\varphi }_2\varphi _1}{|\varphi _1\times \varphi _2|}}|_{𝒙_\mathrm{𝟎},t_0}.`$
We calculate the velocity at each point where the strings cross plaquettes. For this purpose we need to evaluate $`\dot{\varphi }_a`$ and $`\varphi _a`$ at arbitrary points on a plaquette, as shown explicitly in Appendix B, where $`\dot{\varphi }_a`$ and $`\varphi _a`$ are given by the quantities on lattice points up to the second order. Collecting the value of Eq. (19) at each intersection of strings with plaquettes, we can obtain the average of the velocity, the velocity squared, and the Lorentz factor.

### II.4 Intercommutation of strings

Now we discuss the energy loss rate from infinite strings to loops by calculating the intercommutation rate $`c`$. In simulations of cosmic strings based on the Nambu-Goto action, one has to assign the intercommutation probability of intersecting strings by hand due to the lack of microscopic information. Since the strings are identified from the scalar field configuration in our simulation, we can unambiguously calculate how strings intersect with each other, namely, whether they simply pass through or intercommute upon collision. We can therefore identify new loops by monitoring the lengths of all long strings and loops at each time step and comparing the data with those at the previous time step.

### II.5 NG boson emission

It is very difficult to directly evaluate the NG boson emission rate because separation of the emitted axions and NG phase associated with the string core is nontrivial. Since loops also emit these particles, it is a formidable task to identify axions radiated from long strings only. However, we have already described how we can obtain $`\rho _{\mathrm{}}`$, $`v^2`$, and $`\mathrm{\Gamma }_{\mathrm{loop}}`$ from simulation data at each time step, and $`d\rho _{\mathrm{}}/dt`$ can also be calculated from $`\rho _{\mathrm{}}`$ at two adjacent time steps. So all the quantities in the master equation (2) can be calculated from our simulation data except for the last term, which can be found from the equation itself. In particular, in the scaling regime, when the one-scale description suffices, we can find $`\kappa `$ from

$$\kappa =\frac{1v^2}{\sqrt{\xi }}c.$$
(20)
## III Results

### III.1 Scaling parameter

In Fig. 4, the time evolution of $`\xi `$ is depicted for all cases \[(a)-(e)\]. The filled squares represent the time development of $`\xi `$ for the new identification method. For comparison, the results of our previous identification method based on the potential energy density YKY ; YYK , which slightly overestimates $`\xi `$ by a factor of $`1.2`$, are shown by blank circles. The results of the original Vachaspati-Vilenkin identification scheme are indicated by blank squares, obtained by counting the number of boxes through which a string passes as was done in Ref. MSM . This method overestimates $`\xi `$ by a factor of $`1.4`$.

We find that differences in the simulation settings except the box size normalized by the final horizon scale do not affect the overall behavior significantly. $`\xi `$ becomes smaller in the smallest box size with $`=H^1`$ at the final time. This is mainly because under the periodic boundary conditions, a string feels an attractive force from the boundary and tends to disappear. Indeed, in the case (a) this artificial effect starts to operate before the system really relaxes to the scaling solution, so $`\xi `$ continues decreasing without a plateau. On the other hand, in the case (b), the boundary effect becomes operative after $`t130`$ when $`\xi `$ starts to decrease from the scaling value. The other three cases, whose box size is larger than $`2H^1`$, are free from the boundary effects and the scaling parameter $`\xi `$ remains constant. Hence we conclude that after the relaxation period the global string network goes into the scaling regime. From Table 2, which shows average values of $`\xi `$ after some relaxation period ($`t>80`$), $`\xi `$ is found to be $`\xi 0.80`$ in the scaling regime.

### III.2 String velocity

In Figs. 57, the time evolutions of the average velocity $`v`$, the average square velocity $`v^2`$, and the Lorentz factor $`\gamma `$ are shown. After some relaxation period, they become fairly constant in all cases. The average values are given in Table 2. We find that the average values become larger in the smallest box size. This is mainly because under the periodic boundary conditions, a string feels an attractive force from the boundary and tends to be accelerated. However, such an effect becomes negligible for a box size larger than $`2H^1`$ as observed in Table 2. Thus we find the average velocity $`v0.60`$ and the average square velocity $`v^20.50v^2`$, in the scaling regime. On the other hand, the average Lorentz factor has a large scatter in time, although the long-time average is fairly constant with $`\overline{\gamma }1.8`$. Since string segments moving with a speed close to the light velocity have extremely large Lorentz factors and push up the average dramatically, the fluctuation in the number of such string segments results in this large scatter. This is the reason the average Lorentz factor $`\gamma =1.8\pm 0.2`$ is much larger than $`1/\sqrt{1v^2}`$ and $`1/\sqrt{1v^2}`$. Thus the energy per unit length of a string is enhanced by a factor $`\overline{\gamma }1.8`$ compared with a static string.

### III.3 Energy dissipation coefficients

As explained in Sec. II.4 we can obtain the intercommutation rate $`c`$ without assigning its probability by hand. The result is given in Table 2, where we find that $`c`$ is larger in the smallest box size. This can be understood in the same way, namely, under the periodic boundary conditions, a string feels an attractive force from the boundary and tends to intercommute more often. Since such an effect is unimportant for a box size larger than $`2H^1`$ as discussed above, $`c`$ is found to be $`c=0.40\pm 0.04`$.

Taking account of the relation (20), the emission parameter $`\kappa `$ is given in Table 2. In the simulations with the smallest box, strings intercommute more often due to the boundary effect, which suppresses NG boson emission. But again such an effect is inoperative for cases with box size larger than $`2H^1`$. So we conclude that the NG boson emission rate is $`\kappa =0.16\pm 0.04`$.

## IV Discussion and conclusion

In this paper we have investigated the cosmological evolution of the global string network in detail. In our numerical simulations, the equations of motion of the two-component real scalar field are solved on Eulerian meshes. In order to follow the time evolution of global strings, we have developed a new identification method which enables us to find a more correct position of the string core in a box beyond the lattice spacing. Furthermore, we have given a detailed explanation of the method to extract Lagrangian quantities, such as velocity and intercommutation, characterizing the evolution of global strings. The NG boson emission rate is obtained from the master equation (2). Thus the quantitative nature of the cosmological evolution of the global string network is elucidated without setting the intercommutation probability of two intersecting strings by hand. Specifically, we find the scaling parameter characterizing the energy density $`\xi 0.80`$, the peculiar velocity $`v0.60`$, the velocity squared $`v^20.50`$, the Lorentz factor $`\gamma =1.8\pm 0.2`$, the formation rate of loops $`c=0.40\pm 0.04`$, and the emission rate of Nambu-Goldstone bosons $`\kappa =0.16\pm 0.04`$ from our results.

In our simulations, due to the limit of the dynamic range of the lattices, the logarithm of the ratio of the Hubble radius to the string width, which appears in the expression for the effective line density of global strings, took $`\mathrm{ln}(t/\delta )5`$, while in the cosmological setting it can be as large as $`𝒪(100)`$. The authors of MS ; MSM  argue that $`\kappa `$ is inversely proportional to $`\mathrm{ln}(t/\delta )`$, although $`c`$ and $`v^2`$ do not have such dependence. Based on this speculation, they claimed that we obtained a smaller value of the scaling parameter $`\xi `$ in our previous work YKY  than it should really be, because $`\kappa `$ was incorrectly large there. They also argued that in the cosmological situation loop production is the dominant mechanism of energy dissipation of long strings and that $`\xi `$ would take the same value as for local strings.

Now that we have all the values of relevant quantities from our simulation data, we can disprove their claim. Indeed, we find that $`\kappa `$ is smaller than $`c`$ and hence direct emission of NG bosons does not dominate the energy loss mechanism even in our setting of numerical simulations with a relatively small dynamic range, $`\mathrm{ln}(t/\delta )`$. If we extrapolate the value of $`\kappa `$ with the scaling $`\kappa 1/\mathrm{ln}(t/\delta )`$, we find $`\kappa 0.01`$ in the cosmological situation. Even then the scaling parameter (4) does not increase much, yielding $`\xi =1.6\pm 0.3`$, which is still much smaller than that of local strings.

Thus we conclude that there is a quantitative difference between the cosmological evolution of global strings and that of local strings. It is based on the qualitative difference that global strings have a long-range force and intercommute more often with larger $`v^2`$ than do local strings. Note that it is by no means surprising that global defects behave differently from local defects in cosmology. In the case of monopoles, as discussed in the Introduction, global monopoles evolve in a strikingly different manner than do magnetic monopoles. In the case of strings, both local and global defects relax to a scaling solution and their difference is more subtle. So it is not until our numerical analysis of the evolution of the latter from the Lagrangian point of view is performed that their difference is fully elucidated.

Finally, we use our results to obtain a constraint on the symmetry breaking scale $`\eta f_a`$ of the Peccei-Quinn U(1) symmetry. From $`\xi =1.6\pm 0.3`$ and $`\gamma =1.8\pm 0.2`$, we find $`f_a(0.161.2)\times 10^{12}`$ GeV for the normalized Hubble parameter $`h=0.7`$ YKY . We also note that the Lagrangian method developed in this paper is directly applicable to other species of extended objects, for instance, topological defects such as local strings and nontopological solitons such as Q balls as well.

## Note added in proof

After we submitted the original manuscript, we became aware of an analytic estimate of $`\xi `$ which reports that it is $`1`$ even in the case of the local string network. We are grateful to M. Hindmarsh for informing us of his result Hindmarsh .

## Acknowledgments

M.Y. is grateful to Robert Brandenberger for his hospitality at Brown University, where the final part of the work was done. J.Y. would like to thank Toru Tsuribe for useful comments on numerical analysis. This work was partially supported by JSPS Grants-in-Aid for Scientific Research, No. 12-08555 (M.Y.) and No. 13640285 (J.Y.).

## Appendix A Discretization of differential equations

The equation of motion of the scalar fields is given by

$$\ddot{\varphi _a}(x)+3H(t)\dot{\varphi _a}(x)\frac{1}{R(t)^2}^2\varphi _a(x)+\frac{V}{\varphi _a}=0.$$
(21)
In order to discretize the above equation, it is reduced to first-order differential equations,

$`\dot{\varphi }_a`$ $``$ $`\pi _a`$
$`\dot{\pi }_a`$ $`=`$ $`3H(t)\pi _a{\displaystyle \frac{1}{R(t)^2}}^2\varphi _a(x){\displaystyle \frac{V}{\varphi _a}}.`$
Expanding $`\varphi _a(t,𝒙)`$ and $`\varphi _a(t+\mathrm{\Delta }t,𝒙)`$ around the intermediate time step $`t+\frac{1}{2}\mathrm{\Delta }t`$, we find

$$\varphi _a(t+\mathrm{\Delta }t,𝒙)\varphi _a(t,𝒙)=\mathrm{\Delta }t\dot{\varphi }_a(t+\frac{1}{2}\mathrm{\Delta }t,𝒙)+𝒪\left((\mathrm{\Delta }t)^3\right),$$
(23)
up to the third order. Thus $`\pi _a`$ at the intermediate time step is given by

$$\pi _a(t+\frac{1}{2}\mathrm{\Delta }t,𝒙)=\frac{\varphi _a(t+\mathrm{\Delta }t,𝒙)\varphi _a(t,𝒙)}{\mathrm{\Delta }t}+𝒪\left((\mathrm{\Delta }t)^2\right).$$
(24)
In the same way, $`\dot{\pi }_a`$ and $`\pi _a`$ at the time step $`t`$ are represented by the quantities at the two adjacent intermediate steps $`t\pm \frac{1}{2}\mathrm{\Delta }t`$,

$`\dot{\pi }_a(t,𝒙)`$ $`=`$ $`{\displaystyle \frac{\pi _a(t+\frac{1}{2}\mathrm{\Delta }t,𝒙)\pi _a(t\frac{1}{2}\mathrm{\Delta }t,𝒙)}{\mathrm{\Delta }t}}+𝒪\left((\mathrm{\Delta }t)^2\right),`$
$`\pi _a(t,𝒙)`$ $`=`$ $`{\displaystyle \frac{\pi _a(t+\frac{1}{2}\mathrm{\Delta }t,𝒙)+\pi _a(t\frac{1}{2}\mathrm{\Delta }t,𝒙)}{2}}+𝒪\left((\mathrm{\Delta }t)^2\right).`$ (25)
The second-order derivatives are approximated up to the second order in $`\mathrm{\Delta }`$ as

$`{\displaystyle \frac{^2\varphi _a(t,𝒙)}{x^2}}`$ $`=`$ $`{\displaystyle \frac{\varphi _a(t,x+\mathrm{\Delta },y,z)2\varphi _a(t,x,y,z)+\varphi _a(t,x\mathrm{\Delta },y,z)}{\mathrm{\Delta }^2}}+𝒪(\mathrm{\Delta }^2),`$
$`{\displaystyle \frac{^2\varphi _a(t,𝒙)}{y^2}}`$ $`=`$ $`{\displaystyle \frac{\varphi _a(t,x,y+\mathrm{\Delta },z)2\varphi _a(t,x,y,z)+\varphi _a(t,x,y\mathrm{\Delta },z)}{\mathrm{\Delta }^2}}+𝒪(\mathrm{\Delta }^2),`$
$`{\displaystyle \frac{^2\varphi _a(t,𝒙)}{z^2}}`$ $`=`$ $`{\displaystyle \frac{\varphi _a(t,x,y,z+\mathrm{\Delta })2\varphi _a(t,x,y,z)+\varphi _a(t,x,y,z\mathrm{\Delta })}{\mathrm{\Delta }^2}}+𝒪(\mathrm{\Delta }^2).`$ (26)
Thus the fundamental equations are discretized up to the second order in both space and time,

$`\ddot{\varphi _a}(x)+3H(t)\dot{\varphi _a}(x){\displaystyle \frac{1}{R(t)^2}}^2\varphi _a(x)+{\displaystyle \frac{V[\varphi _a(x)]}{\varphi _a}}`$
$`={\displaystyle \frac{\dot{\varphi }_a(t+\frac{1}{2}\mathrm{\Delta }t,𝒙)\dot{\varphi }_a(t\frac{1}{2}\mathrm{\Delta }t,𝒙)}{\mathrm{\Delta }t}}+3H(t){\displaystyle \frac{\dot{\varphi }_a(t+\frac{1}{2}\mathrm{\Delta }t,𝒙)+\dot{\varphi }_a(t\frac{1}{2}\mathrm{\Delta }t,𝒙)}{2}}`$
$`{\displaystyle \frac{1}{R(t)^2}}{\displaystyle \underset{l=x,y,z}{}}\left\{{\displaystyle \frac{\varphi _a(t,𝒙+𝚫_l)2\varphi _a(t,𝒙)+\varphi _a(t,𝒙𝚫_l)}{\mathrm{\Delta }^2}}\right\}+{\displaystyle \frac{V[\varphi _a(t,𝒙)]}{\varphi _a}}+𝒪((\mathrm{\Delta }t)^2,\mathrm{\Delta }^2)`$
$`=0,`$ (27)
where $`\mathrm{\Delta }_x=(\mathrm{\Delta },0,0),\mathrm{\Delta }_y=(0,\mathrm{\Delta },0)`$, and $`\mathrm{\Delta }_z=(0,0,\mathrm{\Delta })`$.

In summary, in our numerical simulations

$`\dot{\varphi }_a(t+{\displaystyle \frac{1}{2}}\mathrm{\Delta }t,𝒙)`$ $`=`$ $`{\displaystyle \frac{1}{1+3H(t)\mathrm{\Delta }t/2}}[(1{\displaystyle \frac{3H(t)\mathrm{\Delta }t}{2}})\dot{\varphi }_a(t{\displaystyle \frac{1}{2}}\mathrm{\Delta }t,𝒙)`$
$`+\mathrm{\Delta }t({\displaystyle \frac{1}{R(t)^2}}{\displaystyle \underset{l=x,y,z}{}}\left\{{\displaystyle \frac{\varphi _a(t,𝒙+𝚫_l)2\varphi _a(t,𝒙)+\varphi _a(t,𝒙𝚫_l)}{\mathrm{\Delta }^2}}\right\}{\displaystyle \frac{V[\varphi _a(t,𝒙)]}{\varphi _a}})],`$
$`\varphi _a(t+\mathrm{\Delta }t,𝒙)`$ $`=`$ $`\varphi _a(t,𝒙)+\mathrm{\Delta }t\dot{\varphi }_a(t+{\displaystyle \frac{1}{2}}\mathrm{\Delta }t,𝒙).`$ (28)
The value of $`\dot{\varphi }_a`$ at the time step $`t`$, which is required to calculate the string velocity, is evaluated from Eq. (25).

## Appendix B Quantities at an arbitrary point on a plaquette

In order to evaluate the velocity of a string correctly, we need to evaluate quantities at an arbitrary point on a plaquette. That is, $`\dot{\varphi }_a`$ and $`\varphi _a`$ within a plaquette should be expressed by their values at its four corners. As an example, let us consider a plaquette parallel to the $`z`$ plane with four corners $`(x,y,z)`$, $`(x+\mathrm{\Delta },y,z)`$, $`(x,y+\mathrm{\Delta },z)`$, and $`(x+\mathrm{\Delta },y+\mathrm{\Delta },z)`$ and express $`\dot{\varphi }_a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z)`$ and $`\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z)`$, which we denote collectively by $`\mathrm{\Pi }(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })`$ below, in terms of their values at these four points. Here $`0\alpha 1`$ and $`0\beta 1`$.

We express $`\mathrm{\Pi }`$ at the four corners by an expansion around $`(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })`$ as

$`\mathrm{\Pi }(x,y)`$ $`=`$ $`\mathrm{\Pi }(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })\alpha \mathrm{\Delta }{\displaystyle \frac{\mathrm{\Pi }}{x}}\beta \mathrm{\Delta }{\displaystyle \frac{\mathrm{\Pi }}{y}}`$
$`+{\displaystyle \frac{1}{2}}\alpha ^2\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{x^2}}+\alpha \beta \mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{xy}}+{\displaystyle \frac{1}{2}}\beta ^2\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{y^2}}+𝒪(\mathrm{\Delta }^3),`$
$`\mathrm{\Pi }(x+\mathrm{\Delta },y)`$ $`=`$ $`\mathrm{\Pi }(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })+(1\alpha )\mathrm{\Delta }{\displaystyle \frac{\mathrm{\Pi }}{x}}\beta \mathrm{\Delta }{\displaystyle \frac{\mathrm{\Pi }}{y}}`$
$`+{\displaystyle \frac{1}{2}}(1\alpha )^2\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{x^2}}(1\alpha )\beta \mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{xy}}+{\displaystyle \frac{1}{2}}\beta ^2\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{y^2}}+𝒪(\mathrm{\Delta }^3),`$
$`\mathrm{\Pi }(x,y+\mathrm{\Delta })`$ $`=`$ $`\mathrm{\Pi }(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })\alpha \mathrm{\Delta }{\displaystyle \frac{\mathrm{\Pi }}{x}}+(1\beta )\mathrm{\Delta }{\displaystyle \frac{\mathrm{\Pi }}{y}}`$
$`+{\displaystyle \frac{1}{2}}\alpha ^2\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{x^2}}\alpha (1\beta )\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{xy}}+{\displaystyle \frac{1}{2}}(1\beta )^2\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{y^2}}+𝒪(\mathrm{\Delta }^3),`$
$`\mathrm{\Pi }(x+\mathrm{\Delta },y+\mathrm{\Delta })`$ $`=`$ $`\mathrm{\Pi }(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })+(1\alpha )\mathrm{\Delta }{\displaystyle \frac{\mathrm{\Pi }}{x}}+(1\beta )\mathrm{\Delta }{\displaystyle \frac{\mathrm{\Pi }}{y}}`$ (29)
$`+{\displaystyle \frac{1}{2}}(1\alpha )^2\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{x^2}}+(1\alpha )(1\beta )\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{xy}}+{\displaystyle \frac{1}{2}}(1\beta )^2\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{y^2}}+𝒪(\mathrm{\Delta }^3).`$
Making an appropriate combination, $`\mathrm{\Pi }(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })`$ can be expressed by its values at the four vertices in the plaquette,

$`(1\alpha )(1\beta )\mathrm{\Pi }(x,y)+\alpha (1\beta )\mathrm{\Pi }(x+\mathrm{\Delta },y)+(1\alpha )\beta \mathrm{\Pi }(x,y+\mathrm{\Delta })+\alpha \beta \mathrm{\Pi }(x+\mathrm{\Delta },y+\mathrm{\Delta })`$
$`=\mathrm{\Pi }(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })+{\displaystyle \frac{1}{2}}\alpha (1\alpha )\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{x^2}}+{\displaystyle \frac{1}{2}}\beta (1\beta )\mathrm{\Delta }^2{\displaystyle \frac{^2\mathrm{\Pi }}{y^2}}+𝒪(\mathrm{\Delta }^3)`$ (30)
$`=\mathrm{\Pi }(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })+𝒪(\mathrm{\Delta }^2).`$ (31)
In particular, replacing $`\mathrm{\Pi }(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })`$ by $`\dot{\varphi }_a(x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta })`$, it can be expressed by the quantities on the lattice points up to second order,

$`\dot{\varphi }_a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z)`$ $`=`$ $`(1\alpha )(1\beta )\dot{\varphi }_a(t,x,y,z)+\alpha (1\beta )\dot{\varphi }_a(t,x+\mathrm{\Delta },y,z)`$ (32)
$`+(1\alpha )\beta \dot{\varphi }_a(t,x,y+\mathrm{\Delta },z)+\alpha \beta \dot{\varphi }_a(t,x+\mathrm{\Delta },y+\mathrm{\Delta },z)+𝒪(\mathrm{\Delta }^2).`$
In the same way, inserting $`\varphi _a/x`$ into $`\mathrm{\Pi }`$, we find

$`{\displaystyle \frac{\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z)}{x}}`$ $`=`$ $`(1\alpha )(1\beta ){\displaystyle \frac{\varphi _a(t,x,y,z)}{x}}+\alpha (1\beta ){\displaystyle \frac{\varphi _a(t,x+\mathrm{\Delta },y,z)}{x}}`$ (34)
$`+(1\alpha )\beta {\displaystyle \frac{\varphi _a(t,x,y+\mathrm{\Delta },z)}{x}}+\alpha \beta {\displaystyle \frac{\varphi _a(t,x+\mathrm{\Delta },y+\mathrm{\Delta },z)}{x}}+𝒪(\mathrm{\Delta }^2)`$
$`=`$ $`{\displaystyle \frac{1}{2\mathrm{\Delta }}}[(1\alpha )(1\beta )\{\varphi _a(t,x+\mathrm{\Delta },y,z)\varphi _a(t,x\mathrm{\Delta },y,z)\}`$
$`+\alpha (1\beta )\left\{\varphi _a(t,x+2\mathrm{\Delta },y,z)\varphi _a(t,x,y,z)\right\}`$
$`+(1\alpha )\beta \left\{\varphi _a(t,x+\mathrm{\Delta },y+\mathrm{\Delta },z)\varphi _a(t,x\mathrm{\Delta },y+\mathrm{\Delta },z)\right\}`$
$`+\alpha \beta \{\varphi _a(t,x+2\mathrm{\Delta },y+\mathrm{\Delta },z)\varphi _a(t,x,y+\mathrm{\Delta },z)\}]+𝒪(\mathrm{\Delta }^2).`$
Here we have used the relations

$`{\displaystyle \frac{\varphi _a(x+\mathrm{\Delta },y,z)\varphi _a(x\mathrm{\Delta },y,z)}{2\mathrm{\Delta }}}`$ $`=`$ $`{\displaystyle \frac{\varphi _a(x,y,z)}{x}}+𝒪(\mathrm{\Delta }^2),`$
$`{\displaystyle \frac{\varphi _a(x+2\mathrm{\Delta },y,z)\varphi _a(x,y,z)}{2\mathrm{\Delta }}}`$ $`=`$ $`{\displaystyle \frac{\varphi _a(x+\mathrm{\Delta },y,z)}{x}}+𝒪(\mathrm{\Delta }^2).`$ (35)
Interchanging $`x`$ for $`y`$ and $`\alpha `$ for $`\beta `$, $`\varphi _a/y`$ is also given by

$`{\displaystyle \frac{\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z)}{y}}`$ $`=`$ $`(1\alpha )(1\beta ){\displaystyle \frac{\varphi _a(t,x,y,z)}{y}}+\alpha (1\beta ){\displaystyle \frac{\varphi _a(t,x+\mathrm{\Delta },y,z)}{y}}`$ (37)
$`+(1\alpha )\beta {\displaystyle \frac{\varphi _a(t,x,y+\mathrm{\Delta },z)}{y}}+\alpha \beta {\displaystyle \frac{\varphi _a(t,x+\mathrm{\Delta },y+\mathrm{\Delta },z)}{y}}+𝒪(\mathrm{\Delta }^2)`$
$`=`$ $`{\displaystyle \frac{1}{2\mathrm{\Delta }}}[(1\alpha )(1\beta )\{\varphi _a(t,x,y+\mathrm{\Delta },z)\varphi _a(t,x,y\mathrm{\Delta },z)\}`$
$`+\alpha (1\beta )\left\{\varphi _a(t,x+\mathrm{\Delta },y+\mathrm{\Delta },z)\varphi _a(t,x+\mathrm{\Delta },y\mathrm{\Delta },z)\right\}`$
$`+(1\alpha )\beta \left\{\varphi _a(t,x,y+2\mathrm{\Delta },z)\varphi _a(t,x,y,z)\right\}`$
$`+\alpha \beta \{\varphi _a(t,x+\mathrm{\Delta },y+2\mathrm{\Delta },z)\varphi _a(t,x+\mathrm{\Delta },y,z)\}]+𝒪(\mathrm{\Delta }^2).`$
Using the relations obtained above and Eq. (30), $`\varphi _a/z`$ is calculated as

$`{\displaystyle \frac{\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z)}{z}}`$ $`=`$ $`{\displaystyle \frac{\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z+\mathrm{\Delta })\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z\mathrm{\Delta })}{2\mathrm{\Delta }}}+𝒪(\mathrm{\Delta }^2)`$ (38)
$`={\displaystyle \frac{1}{2\mathrm{\Delta }}}[(1\alpha )(1\beta )\{\varphi _a(t,x,y,z+\mathrm{\Delta })\varphi _a(t,x,y,z\mathrm{\Delta })\}`$
$`+\alpha (1\beta )\left\{\varphi _a(t,x+\mathrm{\Delta },y,z+\mathrm{\Delta })\varphi _a(t,x+\mathrm{\Delta },y,z\mathrm{\Delta })\right\}`$
$`+(1\alpha )\beta \left\{\varphi _a(t,x,y+\mathrm{\Delta },z+\mathrm{\Delta })\varphi _a(t,x,y+\mathrm{\Delta },z\mathrm{\Delta })\right\}`$
$`+\alpha \beta \{\varphi _a(t,x+\mathrm{\Delta },y+\mathrm{\Delta },z+\mathrm{\Delta })\varphi _a(t,x+\mathrm{\Delta },y+\mathrm{\Delta },z\mathrm{\Delta })\}]`$
$`{\displaystyle \frac{1}{4}}\alpha (1\alpha )\mathrm{\Delta }\left\{{\displaystyle \frac{^2\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z+\mathrm{\Delta })}{x^2}}{\displaystyle \frac{^2\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z\mathrm{\Delta })}{x^2}}\right\}`$
$`{\displaystyle \frac{1}{4}}\beta (1\beta )\mathrm{\Delta }\left\{{\displaystyle \frac{^2\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z+\mathrm{\Delta })}{y^2}}{\displaystyle \frac{^2\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z\mathrm{\Delta })}{y^2}}\right\}`$
$`+𝒪(\mathrm{\Delta }^2).`$
Although these approximations appear of first-order accuracy, using the relations,

$`{\displaystyle \frac{^2\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z+\mathrm{\Delta })}{x^2}}{\displaystyle \frac{^2\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z\mathrm{\Delta })}{x^2}}`$
$`=2\mathrm{\Delta }{\displaystyle \frac{^3\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z)}{x^2z}}+𝒪(\mathrm{\Delta }^3),`$
$`{\displaystyle \frac{^2\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z+\mathrm{\Delta })}{y^2}}{\displaystyle \frac{^2\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z\mathrm{\Delta })}{y^2}}`$
$`=2\mathrm{\Delta }{\displaystyle \frac{^3\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z)}{y^2z}}+𝒪(\mathrm{\Delta }^3),`$ (39)
we find they have second-order accuracy. As a result, up to the second order, we find

$`{\displaystyle \frac{\varphi _a(t,x+\alpha \mathrm{\Delta },y+\beta \mathrm{\Delta },z)}{z}}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{\Delta }}}[(1\alpha )(1\beta )\{\varphi _a(t,x,y,z+\mathrm{\Delta })\varphi _a(t,x,y,z\mathrm{\Delta })\}`$ (40)
$`+\alpha (1\beta )\left\{\varphi _a(t,x+\mathrm{\Delta },y,z+\mathrm{\Delta })\varphi _a(t,x+\mathrm{\Delta },y,z\mathrm{\Delta })\right\}`$
$`+(1\alpha )\beta \left\{\varphi _a(t,x,y+\mathrm{\Delta },z+\mathrm{\Delta })\varphi _a(t,x,y+\mathrm{\Delta },z\mathrm{\Delta })\right\}`$
$`+\alpha \beta \{\varphi _a(t,x+\mathrm{\Delta },y+\mathrm{\Delta },z+\mathrm{\Delta })\varphi _a(t,x+\mathrm{\Delta },y+\mathrm{\Delta },z\mathrm{\Delta })\}]`$
$`+𝒪(\mathrm{\Delta }^2).`$