Extending the velocity-dependent one-scale model for domain walls

We report on an extensive study of the evolution of domain wall networks in Friedmann-Lema\^{\i}tre-Robertson-Walker universes by means of the largest currently available field-theory simulations. These simulations were done in $4096^3$ boxes and for a range of different fixed expansion rates, as well as for the transition between the radiation and matter eras. A detailed comparison with the velocity-dependent one-scale (VOS) model shows that this cannot accurately reproduce the results of the entire range of simulated regimes if one assumes that the phenomenological energy loss and momentum parameters are constants. We therefore discuss how a more accurate modeling of these parameters can be done, specifically by introducing an additional mechanism of energy loss (scalar radiation, which is particularly relevant for regimes with relatively little damping) and a modified momentum parameter which is a function of velocity (in analogy to what was previously done for cosmic strings). We finally show that this extended model, appropriately calibrated, provides an accurate fit to our simulations.


I. INTRODUCTION
It is generally accepted that phase transitions occurred during the early stages of the evolution of the Universe. Among their possible consequences is the production of topological defects through the Kibble mechanism [1]. Two-dimensional topological defects (domain walls) are tightly constrained, unless they are very light or they decay soon after formation, since otherwise they would dominate the energy density of the universe, in disagreement with observations [2]. On the other hand, one-dimensional objects (cosmic strings) are in principle more benign, although they are also subject to increasingly strong constraints [3]. Nonetheless, cosmic strings could play an important role as a relic of fundamental theories of the early Universe, such as brane inflation scenarios [4,5] or supersymmetric grand unified theories (GUT) [6].
To understand the observational effects of the presence of topological defects, a quantitative understanding of the evolution of their networks is essential. Such quantitative analytic models were first obtained for cosmic strings [7,8], and subsequently for domain walls [9]. Meanwhile, the latter can be more easily simulated numerically at higher spatial resolution and dynamic range. For this reason, in addition to their intrinsic relevance, domain walls also provide a useful testbed for the evolution of cosmic strings and superstrings [10]. Still, the velocitydependent one-scale model (VOS) for domain walls is currently less developed than its cosmic string counterpart. Here we take advantage of recent improvements in hardware and computing power to improve this situation.
Specifically, we build upon the work done in [11,12] and carry out an extensive set of high-resolution field theory simulations of domain wall networks using the PRS algorithm [13]. Compared to this earlier work our simulations are both larger (4096 3 boxes, the largest currently available) and span a more diverse set of conditions, including simulations with fixed expansion rates (radiation era, matter era and 10 other expansion rates) as well as, for the first time for domain walls, series of simulations that accurately span the radiation-matter transition. This extended high-resolution dataset enables us to further calibrate and significantly improve the analytic model, as was previously done for strings [14,15].

II. THE STANDARD VOS MODEL FOR DOMAIN WALLS
The analytic VOS model for domain wall network evolution was first obtained from arguments on energy conservation in [9]. Later it was shown that the same result can be reached from a microscopic description [16]. We will revisit and clarify this microscopic approach, and further extend it to shed light on the momentum parameter k for the wall network (to be rigorously defined below).
The wall surface M 2 can be parametrized by two parameters, σ 1 and σ 2 . As a result, the wall evolution is described by the vector x µ (σ 1 , σ 2 , τ ), where we identified σ 0 = τ . [17] If the function x µ (σ 1 , σ 2 , τ ) is smooth, it is possible to parametrize the wall surface in such a way arXiv:1602.01322v2 [hep-ph] 17 Feb 2016 that two tangential vectors will be orthogonal Moreover, we can require that the velocity of the wall ∂ τ x µ ≡ẋ µ can be only normal to the tangent surface T M2 (cf. Fig. 1).
To derive the wall equation of motion we start from the worldvolume (Dirac) action, which has the form where σ w is a constant mass per unit area, γ ab = g µν x µ ,a x ν ,b is the induced metric, γ = 1 3! ab cd γ ac γ bd is its determinant, x µ ,a = ∂x µ ∂σ a , ab is the Levi-Civita symbol, and L is the Lagrangian density.
To obtain equations of motion for a domain wall from Eq. (2), it is useful to use the following equality from which one can obtain where g is the determinant of the metric g µν . The energy of the wall in that case is Let us now define the metric g µν as the FLRW metric with conformal time a(τ )dτ = dt where a(τ ) is the scale factor and dl 2 = dx 2 + dy 2 + dz 2 . Then the equation of motion (Eq. 4) can be rewritten aṡ Let us redefine the coordinates σ 1 and σ 2 to s 1 and s 2 in such way that | ∂x i ∂sα | 2 = 1 (α = 1, 2). This means that derivatives will be changed in the following way (no summation over α). In these new coordinates, it is possible to introduce an orthonormal basis (refer to Fig. 1): ξ i α = ∂x i ∂sα , and n i =ẋ i |ẋ i | . Consequently, the zeroth component of Eq. (8) (λ = 0) can be written aṡ The spatial part (λ = i) of Eq. (8) contracted with the vector n i has the form where k i α = ∂ξ i α ∂sα . The scalar products k i α n i project the curvatures corresponding to σ 1 and σ 2 along the normal vector n i . It should be noted that k i α = a Rα u i α , where u i α are unit vectors and R α are the radii of curvature for σ 1 and σ 2 , respectively. Now it is possible to obtain averaged equations, using the same strategy that was used in [7]. One introduces two macroscopic (averaged) quantities, the energy density and the root-mean-squared (rms) velocity and can thus average Eqs. (10)(11), obtaining where t is a physical time, and H = 1 a da dt is the Hubble parameter, and we made the assumption that curvature radii have the same averaged value and are equal to the correlation length: R 1 = R 2 = L. The K 1 and K 2 parameters are curvature/momentum parameters. The component K 1 can be written as suitably averaged over the network, with an analogous definition for K 2 . As a first approximation they may be assumed to be constants, but later on in this work we will address how they may depend on the velocity υ.

III. SIMULATIONS AND PRELIMINARY CALIBRATION
An evolving wall network loses energy because of possible intersections and the creation of sphere-like objects that eventually collapse. This energy loss mechanism can be added to Eqs. (14) by analogy to what was originally done by Kibble for cosmic strings [18]. This has the form where c w is a constant which we will call (by analogy to the cosmic strings case) the chopping parameter. Taking into account this energy loss term, we can rewrite Eqs. (14) in terms of the correlation length L = σw ρ , as follows where we further defined k w = K 1 +K 2 as the momentum parameter, which we will initially consider as a constant.
For a FLRW universe expanding as a power law, a ∝ t λ , Eqs. (17) have the asymptotic scaling solution L = t and υ = v 0 , where λ, and v 0 are constants [9]. The two phenomenological parameters c w and k w can then be expressed as It should be noted that the right hand-sides of Eqs. (18)(19) are general expressions for the momentum parameter and energy loss mechanisms that can be measured directly and independently from simulations. For this purpose, we should obtain the asymptotic values of the quantities and v 0 from our simulations.
Building upon the work done in Refs. [11,12], we have carried out field theory simulations of the simplest (single-field) domain wall networks in a FLRW background. The equations of motion for the scalar field φ, adopting the Press, Ryden and Spergel procedure [13], can be written in terms of conformal time τ as Relevant numerical parameters are φ 0 = ±1 for the minima of the potential, while the maximum of the potential is V 0 = π 2 /2W 2 0 (where W 0 = 10 is the initial wall thickness in grid units). All these are similar to the ones used in earlier simulations [11][12][13]. Relative to earlier works our simulations have three key advantages • We used a faster and more memory-efficient version of our earlier WALLS code [11,12], optimized for the Intel Xeon Phi architecture.
• This optimization allows us to increase the box size (and therefore the spatial resolution and dynamic range). Specifically, we ran several series of 4096 3 simulations on the COSMOS supercomputer, thus gaining a factor of 8 in volume and a factor of 2 in dynamic range as compared to Ref. [12]. Each simulation starts with τ i = 1 and is stopped when the horizon becomes half the box size (τ f = 2048), ensuring that the periodic boundary conditions of the simulation boxes do not affect the results. Each such simulation requires 1 Tb of memory and takes about 3.7 hours of wall clock time to run on 512 CPUs.
• We explore a much larger range of fixed expansion rates, including radiation era (λ = 1/2), matter era (λ = 2/3) and 10 other expansion rates, ranging from λ = 1/10 to λ = 19/20. (Additionally we also simulated universes during the transition from radiation to matter, to which we will return below.) For each choice of expansion rate we have carried out 10 simulations with different (random) initial conditions: although each of the 10 choices was made randomly, the same 10 choices were used for each of the simulated expansion rates. (This ensures that any differences can be solely ascribed to the different expansion rates.) Unless otherwise stated, the results presented in what follows correspond to the average of each set of 10 runs. Figure  2 illustrates the results of these constant expansion rate simulations.
Note that our choice of initial conditions will lead to large energy gradients in the early timesteps of the simulation, and the network needs some time (which is proportional to the wall thickness) to wash away these initial conditions. This implies that in many grid points the field will go over the top of the potential to get into the other minimum, transiently leading to a relatively small average velocity (the more so the faster the expansion rate), which is clearly visible in the early timesteps in the bottom panel of Fig. 2-note that τ = 10 is the light-crossing time for walls of the average thickness being simulated. This erasing of initial conditions is done in a quasi-coherent way at the various points in the box, leading to the damped oscillations in the average velocity that are also visible in the bottom panel of Fig. 2 (though in this case they are clearer for the slower expansion rates, corresponding to weaker damping).
In practice, since the simulations are evolved in conformal time, the quantity we measure is the conformal correlation length divided by conformal time, which can be straightforwardly related to the physical time quantities Similarly for the velocity we measure γv (or, more precisely, (γv) 2 ), where γ = 1 √ 1−v 2 is the Lorentz factor [11]. In order to identify accurate asymptotic values, we should find the simulation dynamic range when ξ c /τ and γv have already reached the asymptotic behavior and the simulation box still has enough walls for robust statistics (towards the end of each simulation only a few long walls remain, resulting in comparatively poor statistics). After some tests, we conservatively defined the region τ = 500 − 1500 in which our simulations are generally well behaved for all expansion rates (specifically, they are in scaling solutions without significant fluctuations).
Once this region is specified, the averaged values of ξ c /τ and γv can be obtained. These results are presented in Table I. Together with values of ξ c /τ and γv we also list the scaling exponents ν and µ, quantifying convergence to the attractor scaling solution. These are defined as so for a scaling network these exponents should be numerically consistent with µ = −1 and ν = 0. As expected, one finds that the convergence to the scaling solution is faster for faster expansion rates (corresponding to a larger damping term in the wall equations of motion).
Indeed, the ν diagnostic shows that for the slowest expansion rate we have simulated (λ = 1/10) the network has not converged to the scaling behavior and, as a result, it cannot be used for further analysis. In fact this is also qualitatively clear from a simple visual inspection of Fig. 2.
Using the asymptotic values, one can obtain = ξc τ (1−λ) and the velocity v 0 for each expansion rate. By inserting and v 0 into Eqs. (18)(19) one numerically obtains the momentum and chopping parameters. The values thus obtained for each expansion rate are also listed in Table I. It is noteworthy that, with the exception of the For comparison with previous work [11,12], it is interesting to carry out a joint analysis of the data (except the λ = 1/10 case), and determine the best-fit values for these phenomenological parameters if one imposes that they should have the same constant value for all epochs. The results of this analysis are shown in the top panel of Fig. 3, and the following best-fit parameters and uncertainties were found c w = 0.63 ± 0.36, k w = 0.88 ± 0.51.
It should be noted that in this analysis only the statistical errors were taken into account, because it's not possible to determine the systematic error for these simulations. Such errors would include effects such as the PRS approximation [13] and the choice of the range of conformal times in which to do the fits and numerically measure the scaling parameters. However, it is expected that the systematic error becomes smaller for higher expansion rates (since the code is then more robust [11]), just as the statistical error does. Therefore, while the full errors can be larger than presented here, they are broadly expected to have the correct behavior as a function of the expansion rate: higher expansion rates provide more accurate data. As a result, this set of data can be reliably used to calibrate and extend the VOS model.
If we compare the likelihood contours in the top panel of Fig. 3 with the same plot from Ref. [12] (which only had data from three expansion rates, λ = 1/2, 2/3, 4/5), it is seen that they are statistically consistent, but in our case the error bars are significantly larger (as is the reduced chi-squared for the fit). As a comparison, if we repeat the analysis using only the simulations in the range 0.5 ≤ λ ≤ 0.9 we find c w = 0.48 ± 0.24, k w = 1.12 ± 0.31; the results of this analysis are shown in the bottom panel of Fig. 3. The error bars become smaller and the agreement with Ref. [12] is even better. This implies that assuming c w and k w to be constants is not accurate enough for these simulations, if one aims to model a broad range of expansion rates.
More explicitly this can be shown by plotting the righthand sides of Eqs. (18)(19) in terms of the velocity-cf. Fig. 4. For the first of these (top panel), the right-hand side describes the behavior of the momentum parameter k(v) for different expansion rates, while for the second one (bottom panel) it describes an energy loss function which we will denote F (v). For constant k w and c w , Eqs. (18)(19) would imply a constant value for the first plot and a linear function for the second one. Data from our simulations show that this is not the case. As a final, more straightforward check, Table I also lists the numerically inferred c w and k w for each expansion rate. Hence, the momentum parameter should depend on velocity, and the chopping parameter is not sufficient for describing the energy losses.

IV. EXTENDING THE VOS MODEL
We now describe how to extend the analytic VOS model, by more accurately modeling the momentum parameter and the energy loss term.

A. Momentum parameter
The momentum parameter can be estimated in an analogous way to what was done for cosmic strings in Ref. [8]. As we saw, the momentum parameter in our VOS wall model is given by k w = K 1 + K 2 (refer to Sec. II). The component K 1 was previously defined in Eq. (15), suitably averaged over the network, and a similar definition applies to K 2 . We now need to estimate this scalar product in terms of the velocity υ. As can be seen in Fig. 1, there is an orthonormal basis ξ i 1 , ξ i 2 , n i , and therefore we can decompose the vector note that vector ξ i 1 is orthogonal to u i 1 . Therefore u i 1 n i = A and u i 1 u 1i = A 2 + B 2 = 1. When the expansion rate is slow, the velocity squared tends to some maximal value 1/q and perturbations on the wall surface increase. Since the wall surface is highly perturbed in that regime, the averaged value of u i 1 n i goes to zero (A → 0). In the opposite limit when the rate of expansion is fast, the velocity squared tends to zero, and perturbations on the wall surface are very small. As a result, the scalar product u i 1 n i goes to some value k 0 /2. The same considerations apply for K 2 .
Hence, k(v) should reach some value k 0 when the ve- , as numerically determined from the right-hand side of Eqs. (18)(19). The red line in the energy loss plot is a linear function of the rms velocity cwv fitted for high λ (hence low velocity). The blue lines are from the extended analytic model, using phenomenological forms of the momentum parameter (Eq. 25) and energy loss due to scalar radiation (Eq. 27) with the following best-fit parameters d = 0.28, r = 1.30, β = 1.69, k0 = 1.73 and q = 4.27, discussed in the text.
locity is zero and tend to zero when the velocity squared is 1/q. In that case k(v) can be written similarly to the momentum parameter of the string network [8] where β, k 0 and q are unknown parameters. At this point there is one difference between the string and wall cases: there are no non-trivial analytic solutions for walls (like the helicoidal solution for strings) that can be used to infer exact values of q, k 0 and β, as it was done for strings. Consequently, it is only possible to impose physical restrictions on these parameters. The constant k 0 characterizes the maximum value of the momentum parameter: it is positive, but cannot be bigger than 2. The parameter 1/q is an averaged maximal velocity for the wall network. Similarly to what was done for strings in [18], using the general expression for the n-dimensional topological defect dynamics [16], it can be shown that the maximal possible velocity is v 2 max = n n+1 . For walls this is v 2 w = 2/3, as expected, but this result requires a set of assumptions that need not be satisfied. In that case the maximal averaged velocity of the network can be smaller (but not larger). As a result we have Other than these general physical constraints, these parameters must be calibrated numerically. Fortunately, the resolution of our simulations is high enough to enable this calibration, as we will show below.

B. Energy loss mechanisms
The modification of the momentum parameter described above is not sufficient to account for the mismatch between the simulation data and the analytic prediction for the energy losses. We should also improve the modeling of the latter for a better description of the wall network evolution. In addition to the chopping mechanism, another significant contribution to energy losses is expected to be from scalar radiation. Moreover, one may expect it to be proportionally more important (compared to the chopping mechanism) for slower expansion rates.
Energy loss due to scalar radiation was considered in Ref. [19]. It was shown that the uniformly moving wall does not radiate. Only perturbations on the wall surface produce scalar radiation. We have already estimated the level of perturbations in the momentum parameter expression k(v). The maximal value k 0 corresponds to the minimal rms velocity and hence to minimal perturbations on the wall surface. Conversely the case when the momentum parameter is zero corresponds to a maximal rms velocity and a maximally perturbed surface. It looks reasonable to anticipate that the amount of radiation is proportional to the surface perturbations. As a result, we can introduce a modified analytic description of energy losses where d and r are constants. In the maximally perturbed (slow expansion) limit v 2 → 1/q this behaves as and we expect the scalar radiation term to be the dominant one. Conversely in the uniform surface (fast expan-sion) limit we have and in this case we expect the chopping term to be more important, and possibly dominate (depending on the values of the free parameters to be calibrated numerically). Putting together these extensions, the VOS model equations (Eq.17) can finally be rewritten as where k(v) is defined by Eq. (25).

C. Calibrating the extended model
One can easily confirm that the extended VOS model given by Eqs. (30) possess the same scaling behavior as the original one, given by Eqs. (17). In the extended model we have in principle 6 undefined parameters that should be determined from numerical simulation data. By using bootstrapping techniques one finds that the chopping parameter is negligibly small in comparison with the contribution from scalar radiation and may be neglected as a first approximation (specifically, we find c w = 0.00 ± 0.01), while the other five parameters have the following values which are fully consistent with the ones obtained for the full range of expansion rates. While this is not entirely surprising (since the fit is dominated by the high expansion rates, for which the statistical uncertainties of the parameters measured in the simulations are smaller) it is supporting evidence for the fact that the model can accurately describe all expansion rates. We note that in this restricted range the chopping parameter is still negligible (c w = 0.00 ± 0.03), indicating that much smaller velocities (and therefore faster expansion rates than λ = 19/20) would be needed to make this term comparable to the scalar radiation one. Numerically exploring this ultra-fast expansion regime is an interesting but computationally challenging task which we leave for subsequent work.
Alternatively, by solving Eqs. (30) in the scaling regime for different expansion rates λ with the best-fit parameters determined above we can compare the model predictions with the numerically determined quantities. Figure 5 displays this comparison (both for the full and restricted ranges of the expansion rate λ), confirming that the extended model accurately reproduces the simulation data. Moreover, it can be concluded that the main energy loss mechanism for the wall network is generically scalar radiation, with the chopping term only becoming important for fast expansion rates.

V. THE RADIATION-MATTER TRANSITION
Thus far we tested the model against numerical simulations with a fixed expansion rate λ. As an additional test, we have carried out analogous field theory simulations of the radiation-matter transition. In this case the scale factor has the following exact analytic expression where τ * = τ eq /( √ 2 − 1) and the parameters a eq and τ eq are constants denoting the scale factor and conformal time at the epoch of equal radiation and matter densities. For illustration purposes we can also calculate an 'effective' expansion rate during the transition as expected this interpolates between the radiation and matter era values.
In this case we ran various sets of simulations with the same parameters and (random) initial conditions that were described in Sect. III, except that the scale factor obeys Eq. (41). The requirement of sufficient resolution implies that there is not enough memory available for a single simulation to span the entire transition epoch; instead, various sets of runs were carried out starting at (30) with the best-fit parameters described in the text, compared to the data (with statistical error bars) from the numerical simulations for different expansion rates. The solid blue line corresponds to the best-fit parameters for the full range of expansion rates considered while the red dashed one corresponds to the best-fit parameters for the restricted range.
various different conformal times relative to the transition epoch, equally spaced in the logarithm of τ i /τ eq . Figure 6 (to be compared to Fig. 2) summarizes the results of these simulations. Note that the two black solid lines correspond to the radiation (λ = 1/2) and matter (λ = 2/3) simulations already discussed in Sect. III. This is an important test of our code: it shows that simulations evolving sufficiently early and sufficiently late in the transition behave exactly like radiation and matter era simulations-as they must. This figure also makes it visually clear that although the 'early' and 'late' simulations reach scaling (since they are effectively evolving with a constant or quasi-constant expansion rate) this is not case for the ones evolving during the transition itself: in that case the effective expansion rate is changing and the network is constantly trying to adapt (as fast as allowed by causality) to these changing conditions. This is clear in the cyan, green and yellow lines in the plots.
By inserting the scale factor expression (Eq. 41) with the corresponding constants a eq and τ * in the system of Eqs. (30), we can now compare the dynamics of the extended analytic model and the simulations. This comparison is summarized in Fig. 7, where the results of simulations (solid color lines) and the extended analytic model (dashed black line) are compared. It is seen that the analytic model provides an excellent description of the radiation-matter transition. Figure 7. Evolution of the conformal correlation length divided by conformal time ξc/τ (top panel) and of (γv) 2 (bottom panel) during the radiation-matter transition, plotted as a function of the natural logarithm of the scale factor (relative to aeq). The simulations are denoted by solid color lines (each line being an average of 10 simulations with random initial conditions) while the prediction of the extended analytic model with the best-fit parameters discussed in the text is shown by the black dashed lines. The plot only includes the dynamic range 20 ≤ τ ≤ 1500 of each set of simulations; the earlier part (which is dominated by the initial conditions in the box rather than converging to the attractor solution) and the latter part (which has comparatively poor statistics, as discussed in Sect. III) have been omitted for clarity.

VI. CONCLUSIONS
In this paper we revisited the evolution of domain wall networks in expanding FLRW universes. We took advantage of recent progress in computing power and hardware to carry out the largest and most extensive set of field theory simulations of domain walls, using the PRS algorithm. We have simulated expanding universes with 12 different fixed expansion rates, as well as sets of simulations which together span the entire radiation-matter transition.
Our simulations allowed us to significantly improve the analytic description of wall network evolution, based on the quantitative VOS model. We have explicitly shown that a constant momentum parameter k w and chopping parameter c w cannot fully reproduce the simulations for different expansion rates. As a result of this mismatch, we used phenomenological arguments to introduce an extended model, given by Eqs. (30). In this model the momentum parameter is described by a velocity-dependent function k(v) (Eq. 25), and there is a generalized energy loss function F (v) (Eq. 27), which in addition to chopping losses also includes scalar radiation of walls. We did not address the issue of possible losses to gravitational radiation, which is left for subsequent work.
Fitting the phenomenological parameters to the simulations, we found that energy losses due to creation of sphere-like objects are typically subdominant in comparison with scalar radiation, except in the case of fast expansion rates. We have confirmed that the extended analytic model can describe both the fixed expansion rate cases and the transition from the radiation to the matterdominated era. The latter one is an important test of the model, since the network is not scaling during the transition (while the model parameters were calibrated from fixed expansion rate data in the scaling regime).
In the future it would be interesting to extend this analysis to the case of cosmic strings. In that case the chopping term is known to be more important, but signi-ficant and not fully understood differences exist between the results of Goto-Nambu and field theory simulations. As already shown in Ref. [14], the VOS model allows a direct comparison of the results of both types of simulations, and the recent progress in computing power should permit a clarification of this issue.