Understanding Domain Wall Network Evolution

We study the cosmological evolution of domain wall networks in two and three spatial dimensions in the radiation and matter eras using a large number of high-resolution field theory simulations with a large dynamical range. We investigate the dependence of the uncertainty in key parameters characterising the evolution of the network on the size, dynamical range and number of spatial dimensions of the simulations and show that the analytic prediction compares well with the simulation results. We find that there is ample evidence from the simulations of a slow approach of domain wall networks towards a linear scaling solution. However, while at early times the uncertainty in the value of the scaling exponent is small enough for deviations from the scaling solution to be measured, at late times the error bars are much larger and no strong deviations from the scaling solution are found.


Introduction
Topological defects are generic in nature and may be formed whenever a phase transition occurs. In an expanding universe cooling down from a very hot ini-tial state it is to be expected that topological defects may provide a unique window onto the physics of the early universe offering perhaps the best hope of a clear observable link between cosmology and particle physics [1,2]. Most cosmological studies of cosmic defects have focused on cosmic strings due to their interesting properties and strong motivation from fundamental physics (see for example [3,4,5] and references therein). Although standard cosmic strings are now ruled out as the sole contribution to the large scale structure of the universe [6,7] they may still be the dominant source of perturbations on small cosmological scales and may give rise to a number of interesting cosmological consequences [8,9,10]. However, at present there are only two observations for which cosmic strings seem to offer the most natural explanation [11,12].
Domain wall scenarios have attracted less attention since heavy domain walls in a linear scaling regime rapidly dominate the energy density of the universe. Moreover, domain walls which are light enough to satisfy current CMB constraints have a negligible direct contribution to structure formation. However, in this case a number of interesting consequences are possible such as a contribution to the dark energy [13,14] (if domain walls are frozen in comoving coordinates) and a small but measurable contribution to the CMB anisotropies at large angular scales [15]. Domain walls may also separate regions in the universe with different values of the cosmological parameters and/or fundamental constants of nature [16,17].
In this paper we perform a quantitative study of the cosmological evolution of domain wall networks. We investigate the dependence of the uncertainty in key parameters characterising the evolution of the network on the size, dynamical range and number of spatial dimensions of the simulations using a simple analytic model. We then compare our analytic predictions with the results of a large set of high-resolution simulations of domain walls in two and three spatial dimensions [18], using the standard Press-Ryden-Spergel (PRS) algorithm [19] (see also [20,21,22,23]), and discuss the evidence from the simulations of a slow approach towards a linear scaling regime. Previous studies of domain wall network evolution [19,20,21,22,23] having a smaller number of simulations with smaller size and dynamical range than the present one have found some hints for deviations from a scale-invariant evolution. It is therefore crucial to investigate if these are only transitory or if there is a more fundamental reason for such deviations.
The present article is a follow-up of [18]. There, we concentrated on the overall (global) dynamical features of the simulations. On the other hand, having a large dynamic range means that a more localised analysis is also possible, and in particular local exponents can be calculated with relatively small errors. In the present paper we explore this possiblity, and also make use of the large number of simulations to discuss some analytic ways to estimate statistical errors.

Domain Wall Network Evolution
We study the evolution of a domain wall network in a flat homogeneous and isotropic Friedmann-Robertson-Walker (FRW) universe. We consider a scalar field φ with the Lagrangian density and we will take V (φ) to be the generic φ 4 potential with two degenerate minima given by which obviously admits domain wall solutions. Following the procedure described in ref. [19] we modified the equations of motion in such a way that the co-moving thickness of the domain walls is fixed in co-moving coordinates allowing us to resolve the domain walls throughout the full dynamical range of the simulations. With this modification implemented the equations of motion for the field φ become: where a is the scale factor, η is the conformal time and α and β are constants. We take β = 0 in order to have constant co-moving thickness and α = 3 to ensure that the momentum conservation law of the wall evolution in an expanding universe is maintained [19]. Equation (3) is then integrated using a standard finite-difference scheme.
We have verified that the PRS alghoritm gives the correct results in some special cases such as the dynamics of a plane wall or the collapse of a spherical or cilindrical domain wall. We have also verified that it appears to have a small impact on the large-scale dynamics of domain wall networks and does not seem to affect the quantities we want to measure for the purpose of testing scaling properties provided a minimum acceptable tickness is used. However, it is only possible to test the performance of the PRS alghoritm over a narrow window since the 'true' equation of motions for the domain walls rapidly make the wall thickness smaller than the grid size.
In addition to these simple tests with domain walls, the PRS algorithm has been much more extensively used and tested in the case of cosmic strings (see for example ref. [24]). From those one can infer that with the use of the PRS algorithm some quantitative features of the networks will indeed differ (for example, the distribution of small-scale features on the networks), but the broad features will be largely unchanged. An example of the latter is the existence of an attractor scaling solution, and how fast it is reached starting from some given initial configuration. Since this is the issue we are studying here, we believe that the PRS algorithm is adequate for our purposes.
The ratio between the kinetic and potential energy of the domain walls is approximately given by where A is the total co-moving area of the domain walls determined using the algorithm described in ref. [19,22,18] and we are measuring length in units of the grid spacing ∆x (so that ∆x = 1). This quantity is related to the root-mean squared velocity of the domain walls which should be conserved in a linear scaling regime. We assume the initial value of φ to be a random variable between −φ 0 and +φ 0 and the initial value of ∂φ/∂η to be zero. See [22,18,24] for further discussion of these and other issues.

Analytic Modelling
We consider a simple model for the evolution of the uncertainty in key parameters characterising the evolution of a domain wall network which we then test against domain wall network simulations. Let us define the comoving correlation length of the network as ξ ≡ V /A. Consider a cubic grid in N D dimensions where each cube of comoving volume V ξ = ξ N D has N D faces of its own. Note that the number of faces of a cube in N D dimensions is 2N D but each one of them belongs to two adjacent cubes so that in practise each cube has N D faces of its own. We also assume that there is on average one domain wall of comoving area ξ N D −1 per cube of comoving volume V ξ occupying one of the faces of the cube. Hence, the probability that a face of such a cube is occupied by a domain wall is p = 1/N D so that the variance in the number of domain walls per cube is σ 2 If we now consider a cube of comoving volume V = L N D it will contain N = V /V ξ cubes of comoving volume V ξ . The standard deviation, σ X , of X = n/n (where n denotes the number of domain walls of comoving area ξ N D −1 andn its average value) on a given volume, V , is proportional to √ N /N and consequently This can also be calculated for any other variable, Y , characterising the network if Y is proportional to X. In this case Here the random variable Y may represent the total area A, or the ratio between the kinetic and potential energy parametrised by F . Although in the first case it is a fair assumption to draw a direct relation between the number of domain walls and the total area, in the second case we expect the relation to be less direct due to the dispersion in the domain wall velocities and to the existence of an important fraction of the kinetic energy which is not associated with the domain walls themselves but with their decay products.
The scaling exponent may be calculated from the value of R ≡ ηA/V at two different values of the conformal time η 1 and η 2 as so that assuming that ∆ 1 and ∆ 2 are small (compared withR 1 andR 2 respectivelly) and uncorrelated. Here ∆ = (R −R)/R whereR(η) denotes the average value of R(η) over a given number, N S , of simulations. We also may also estimate σ Y for several key parameters, Y = R, F, λ, describing the evolution of the network directly from the simulations as Note that the average value of Y for a given sample of simulations,Ȳ , is the best estimator of the real average of Y with and error of σ Y / √ N S . Also s Y is the best estimator of σ Y and s Y ∼ σ Y if N S is large enough.  Fig. 1. The evolution of R ≡ ηA/V and F for one hundred 4096 2 2D matter era simulations. We also plot the average evolution of these parameters and the expected 1-sigma interval around the mean (using the analytic estimates for σ R and σ F described in Sec. III).

Results and discussion
A total of several thousand matter and radiation era simulations in two and three spatial dimensions were run for various box sizes and dynamical ranges. Here, we compare some of these numerical results with analytic expectations and discuss the main results. In Fig. 1 we plot the evolution of R ≡ ηA/V and F for one hundred 2D matter era simulations. We also plot the average evolution of these parameters and the expected 1-sigma interval around the mean (using the analytic estimates for σ R and σ F described in the previous section). Fig. 2 shows a similar plot but now for one hundred 4096 2 2D radiation era simulations. Figs. 3 and 4 display analogous results for one hundred 256 3 3D matter and radiation era simulations respectivelly.
To zeroth order the results do not seem to be very sensitive on the number of  spatial dimensions, N D . We clearly see that although the network is not scaling during most of the simulation it seems to approach a scaling solution very slowly at late times. Also, the main difference between matter and radiation era simulations lies in the larger value of the typical velocity of the domain walls in the latter case although we also find a considerable difference in the evolution of the matter and radiation era runs at early times when the initial conditions are still important.
As expected we see that in the case of the domain wall area the analytic estimate is a good approximation (σ R is only slightly overestimated -see Figs. 7 and 8) while in the case of the parameter F the uncertainty is underestimated due to not taking into account fact the uncertainty in the velocity of the domain walls as well as the part of the kinetic energy which is not associated with the domain walls themselves but with the particle radiation emitted by them.
Our analytic estimates for the variance of the scaling exponent, λ, were based on the assumption that the statistical properties of the domain wall network  at two different times η 1 and η 2 were uncorrelated. Here we test for the validity of that assumption for different sizes of the conformal time interval η 2 − η 1 . In Fig. 5 we compare the numerical estimates of the standard deviation of the scaling parameter in an ensemble of N S simulations, s λ / √ N S , with s λ computed either directly from the simulation (inner bars) or calculating s ∆ from the simulations and using eqn. (8) to compute s λ assuming no correlation (outer bars) for η 1 = η * = 8 −N D /2 L and various values of η 2 = η for one hundred 4096 2 2D and one hundred 256 3 3D matter era simulations (left and right plots respectivelly). We clearly see that when we increase η 2 − η 1 our assumption of no correlation becomes increasingly more justified and that  Fig. 7. Evolution of the scaling exponent λ as a function of the conformal time η for one hundred 4096 2 2D and one hundred 256 3 3D matter era simulations (left and right plots respectivelly). The binnings have constant dynamical range (with η 2 /η 1 = 2) and the error bars are the standard deviation of the scaling parameter in an ensemble of N S simulations, s λ / √ N S , calculated either using our analytic estimate for σ λ (outer bars) or by calculating s ∆ from the simulations and using eqn. (8) to compute s λ assuming no correlation (inner bars). All the error bars (except the last one in each plot) were artificially enlarged by a factor of (η/η f ) −N D /2 where η f is the conformal time corresponding to the last error bar. there is no significant dependence of the correlation time on the number of spatial dimensions N D . Fig. 6 displays analogous results but now for radiation era simulations. However, no significant differences from the matter era runs are found.
In Fig. 7 we plot the evolution of the scaling exponent λ as a function of the conformal time η for one hundred 4096 2 2D and one hundred 256 3 3D matter era simulations (left and right plots respectivelly). Here, we chose binnings of constant dynamical range (with η 2 /η 1 = 2) and the error bars are the standard deviation of the scaling parameter in an ensemble of N S simulations, s λ / √ N S , calculated either using our analytic estimate for σ λ (outer bars) or by calculating s ∆ from the simulations and using eqn (8) to compute s λ assuming no correlation (inner bars). Note that the error bars in Fig. 7 were artificially enlarged by a factor of of (η/η f ) −N D /2 where η f is the conformal time corresponding to the last error bar on the right so that at early times the real error bars are much smaller than at late times. We clearly see that our simple model successfully predicts the evolution of the uncertainties in the value of the scaling exponent, λ, although the performance of the method may slightly depend on the number of spatial dimensions, N D . We find that, as expected, at early times there are strong deviations from the linear scaling solution but the networks seem to slowly approach a linear scaling solution at late times. Fig. 8 displays analogous results but now for radiation era simulations but no significant differences (as far as the approach to scaling is concerned) were found.
Our results are consistent with those of previous studies [19,20,21,22,23] in finding a slow approach to linear scaling as well as strong deviations from a linear scaling solution at early times. However, we do not find any evidence for strong deviations from a linear scaling solution at late times. We note that this is not in disagreement with the results obtained by other authors but it is a consequence of the larger number, size and dynamical range of the simulations analysed in our paper.