Employing the publicly available osmoattice code we conduct numerical simulations of a domain wall network and the resulting gravitational waves GWs in a radiationdominated Universe in the ZSUB2SUBsymmetric scalar field model. In particular the domain wall evolution is investigated in detail both before and after reaching the scaling regime using the combination of numerical and theoretical methods. We demonstrate that the total area of closed walls is negligible compared to that of a single long wall stretching throughout the simulation box. Therefore the closed walls are unlikely to have a significant impact on the overall network evolution. This is in contrast with the case of cosmic strings where formation of loops is crucial for maintaining the system in the scaling regime. To obtain the GW spectrum we develop a technique that separates physical effects from numerical artefacts arising due to finite box size and nonzero lattice spacing. Our results on the GW spectrum agree well with refs. 2930 which use different codes. Notably we observe a peak at the Hubble scale an exponential falloff at scales shorter than the wall width and a plateaubump at intermediate scales. We also study sensitivity of obtained results on the choice of initial conditions. We find that different types of initial conditions lead to qualitatively similar domain wall evolution in the scaling regime but with important variations translating into different intensities of GWs.
['Cosmic strings', 'domain walls', 'monopoles', 'gravitational waves / theory', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Phenomenology', 'High Energy Physics - Theory']
Revisiting evolution of domain walls and their gravitational radiation with CosmoLattice
{'Revisiting evolution of domain walls and their gravitational radiation with CosmoLattice': "I. Dankovsky a,b , E. Babichev c , D. Gorbunov b,d , S. Ramazanov e , A. Vikman f  \na Faculty of Physics, Lomonosov Moscow State University, Leninskiye Gory 1-2, 119991 Moscow, Russia  \nb Institute for Nuclear Research of the Russian Academy of Sciences, 117312 Moscow, Russia c Universit'e Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay  \nd Moscow Institute of Physics and Technology, 141700 Dolgoprudny, Russia  \ne Institute for Theoretical and Mathematical Physics, MSU, Moscow 119991, Russia  \nf CEICO, FZU-Institute of Physics of the Czech Academy of Sciences,  \nNa Slovance 1999/2, 182 00 Prague 8, Czech Republic", 'Abstract': 'Employing the publicly available C osmo L attice code, we conduct numerical simulations of a domain wall network and the resulting gravitational waves (GWs) in a radiation-dominated Universe in the Z 2 -symmetric scalar field model. In particular, the domain wall evolution is investigated in detail both before and after reaching the scaling regime, using the combination of numerical and theoretical methods. We demonstrate that the total area of closed walls is negligible compared to that of a single long wall stretching throughout the simulation box. Therefore, the closed walls are unlikely to have a significant impact on the overall network evolution. This is in contrast with the case of cosmic strings, where formation of loops is crucial for maintaining the system in the scaling regime. To obtain the GW spectrum, we develop a technique that separates physical effects from numerical artefacts arising due to finite box size and non-zero lattice spacing. Our results on the GW spectrum agree well with Refs. [29, 30], which use different codes. Notably, we observe a peak at the Hubble scale, an exponential falloff at scales shorter than the wall width, and a plateau/bump at intermediate scales. We also study sensitivity of obtained results on the choice of initial conditions. We find that different types of initial conditions lead to qualitatively similar domain wall evolution in the scaling regime, but with important variations translating into different intensities of GWs.', '1 Introduction': 'The first evidence of stochastic gravitational wave (GW) background has been recently found by several collaborations including NANOGrav [1], EPTA 1 with InPTA [2], PPTA [3], and Chinese PTA [4]. While supermassive black hole binaries (SMBHBs) are often thought to be the most likely source of the signal detected [6, 7], the latter can be also explained by various primordial mechanisms [5, 7]. Moreover, the explanation invoking SMBHBs may require highly nontrivial astrophysical mechanisms, see, e.g., Refs. [8, 9, 10]. In light of this and future searches for stochastic GW backgrounds with LISA [11] and the Einstein Telescope [12], it is timely to investigate the fine structure of the GW spectrum generated by a specific primordial source. This detailed analysis could enable us to distinguish it from the GW spectrum of astrophysical objects. In this work we focus on GWs from evolving cosmic domain walls [13].  \nDomain walls arise in field theories with scalars which transform under discrete symmetries spontaneously broken by non-zero vacuum expectation values. For concreteness, we focus on the simplest model of Z 2 -symmetric scalar field ϕ :  \nS = ∫ d 4 x √ -g [ 1 2 g µν ∂ µ ϕ∂ ν ϕ -λ 4 ( ϕ 2 -η 2 ) 2 ] . (1)  \nHere η is the field ϕ expectation value spontaneously breaking Z 2 -symmetry. Domain walls are often viewed as problematic in cosmology because, while the Universe expands, their energy density grows relative to the energy density of the surrounding matter. Thus, unless their tension is chosen to be very small [13, 14], domain walls would promptly become the dominant driving source for cosmological evolution, invalidating the standard expansion history. This problem is cured, e.g., by adding terms, which explicitly break Z 2 -symmetry [13, 15], introducing biased initial conditions [16], or by promoting η to a dynamical quantity [17, 18, 19, 20]. In this work, however, we do not discuss the domain wall problem assuming that it is resolved in one or another way, so that the domain wall network gets destroyed at some moment of time in the early Universe without any contradiction with cosmological observations.  \nEvolution of domain walls has been first studied numerically in Ref. [21]. To make the problem tractable, a non-physical assumption was made that the wall width grows linearly with the scale factor. It was concluded in Ref. [21] that the domain wall network enters a scaling regime with one or a few long domain walls per horizon volume. The scaling regime was confirmed in later simulations [22, 23], which, contrary to [21], kept the wall width constant. Furthermore, Ref. [22] estimated the energy density and spectrum (frequency distribution) of GWs emitted by domain walls. The latter also serve as a source of particle  \nemission studied in the context of axion string-wall networks [24, 25, 26, 27]. Most recently, numerical simulations of domain walls have been focused on proper incorporation of the potential bias [28, 29, 30, 31], which triggers formation of primordial black holes through the collapse of false vacuum bubbles [32, 33, 34]. Note that besides numerical approaches, (semi)-analytical modelling of domain walls is also being actively developed [35, 36, 37, 38, 39, 40, 41].  \nThe main goal of this work is to perform an independent study of domain wall evolution and the resulting GW production at radiation domination, with the help of the publicly available code C osmo L attice [42, 43, 44]. Since the latter can be used by a broad audience, it is an important task to analyse, whether the results obtained with C osmo L attice and different codes [22, 29, 30], as well as analytical results of Sec. 2, agree between each other. Most simulations are performed assuming either vacuum initial conditions (with a cutoff at large momenta) or thermal ones defined in Eqs. (36) and (37), respectively.  \nThe structure of the domain wall network revealed with C osmo L attice is as follows. We observe in Sec. 4 that the network mainly consists of a single wall stretching throughout the box, which simulates the early Universe comprising multiple Hubble patches. At the same time, closed walls constitute only a very small fraction of the total network, see also Ref. [45]. We demonstrate this explicitly for the first time in the literature using histograms, see Fig. 5 in Sec. 4. Neglecting closed walls, there is on average one long domain wall per Hubble volume in the scaling regime, - in accordance with the previous studies [21, 22]. We also provide a quantitative analysis of the domain wall scaling. We find that the scaling is reached when the wall width becomes ∼ 1 / 20 fraction of the horizon radius 2 . This fact is remarkably independent of the choice of initial conditions, even though approaching the scaling solution may occur in drastically different ways. On the other hand, the domain wall area is slightly sensitive to the choice of initial conditions. The similar observation has been made in Ref. [31], when comparing thermal and inflationary initial conditions. We prove that this difference is independent of lattice spacing, box size, and UV properties of the initial spectra, and hence it is likely to have a physical origin.  \nPrevious simulations of GWs from domain walls using C osmo L attice [46, 47] indicated strong departures from the results of Ref. [22]. In particular, Refs. [46, 47] noticed appearance of a pronounced peak in the ultraviolet (UV) part of the spectrum. We also observe the UV peak as the result of our simulations, but we prove that it is an artefact of the non-zero lattice spacing. We propose a technique to attain trustworthy predictions of the GW spectrum by comparing GW spectra obtained with different space resolutions. For this purpose we use N × N × N lattices with different grid numbers N = 512, N = 1024, and N = 2048.  \nOur main results for the GW differential spectra, defined through the energy density of  \nFigure 1: Spectrum of GWs emitted by the domain wall network at radiation domination starting with vacuum (top panel) and thermal (bottom panel) initial conditions defined in Eqs. (36) and (37), respectively. Conformal momenta and conformal times are in units of √ λη and 1 √ λη , respectively. The sharp upper cutoff at k cut = 1 is applied in the case of vacuum initial conditions. The expectation value η is set at η = 6 · 10 16 GeV. Rescaling to arbitrary η is achieved by multiplying the spectra by ( η/ 6 · 10 16 GeV) 4 . Simulations have been performed on a lattice with a grid number N = 2048. The positions of diamonds correspond to the comoving Hubble scale k = 2 πHa at the time associated with the corresponding curves, while stars show the inverse domain wall width 1 /δ w , i.e., k = 2 πa/δ w .  \n<!-- image -->  \nGWs ρ gw ( τ ) and total energy density ρ ( τ ) as  \nΩ gw ( τ ) ≡ dρ gw ( τ ) ρ ( τ ) d log k ,  \nare summarised in Fig. 1 at various conformal times τ and conformal momenta k for vacuum (with a cutoff at large momenta) and thermal initial conditions. At this stage we do not switch on the mechanism of the domain wall destruction (which is needed for the healthy cosmology) and consequently ignore its possible impact on the GW production. The overall growth of the spectrum reflects constancy of GWs produced by domain walls in the scaling regime, which we prove analytically in Sec. 2.2. The infrared (IR) tail of the GW spectrum as well as the behaviour in the vicinity of the main (scaling) peak are in a good agreement with those of Refs. [22, 29, 30]. Nevertheless, the spectral index on the right of the peak is slightly different compared to Refs. [22, 29, 30], as we demonstrate in Sec. 5. At intermediate scales, we observe a plateau for the thermal initial conditions and a bump for the vacuum initial conditions. Such features have been absent in the earlier simulations [22], but can be seen in later analyses of Refs. [29, 30]. Note that presence of the bump in the case of vacuum initial conditions appears to be the consequence of the sharp cutoff. Indeed, artificially applying the sharp cutoff to thermal initial conditions also leads to the bump feature, cf. Fig. 9. Our simulations reveal other differences of GW spectra depending on the choice of initial conditions. In particular, at any time after 3 τ = 5 the heights of GW peaks on the top plot of Fig. 1 exceed those on the bottom plot. That tendency can be traced back to the aforementioned difference of domain wall areas for the two types of initial conditions, see also Sec. 4.  \nThe outline of the paper is as follows. In Sec. 2, we provide a theoretical background for domain wall evolution in the scaling regime and the associated GW emission. In Sec. 3, we set up the strategy for the optimal analysis of the domain wall network on a lattice. In particular, we establish the optimal parameters to master the natural limitations such as a finite box size and a non-zero lattice spacing. We present results of simulations of the domain wall evolution and GW production in Secs. 4 and 5, respectively. We conclude in Sec. 6 with a summary and discussion of the future prospects.', '2 Scaling regime: analytic approach': "In this section we discuss analytic approach to domain wall network evolution and the resulting GWs. On scales much longer than the width, domain walls evolve towards a self-similar and scale-free ('scaling') dynamical attractor regime described in Sec. 2.1. In Sec. 2.2, we calculate the energy density of GWs from domain walls, and show that it is constant in the  \nscaling regime. Unlike Ref. [22], we demonstrate this analytically without referring to the Einstein quadrupole formula. We also briefly review, how the spectral shape of the IR part of GW spectrum can be deduced solely from generic considerations of causality.", '2.1 Domain walls in the scaling regime': "Let us consider some stochastic quantity s ( x ) describing the domain wall network and having zero expectation value, i.e., ⟨ s ( x ) ⟩ = 0, where we average over many realisations of the network with the same microscopic physics and cosmological conditions. This could be, e.g., the expectation value of the scalar field ϕ comprising domain walls. Another example will be given in Sec. 2.2. We characterise statistical properties of the field s ( x ) by the normalised 2-point function:  \n⟨ s ( x , τ ) s ( y , τ ' ) ⟩ √ ⟨ s 2 ( x , τ ) ⟩⟨ s 2 ( y , τ ' ) ⟩ ≡ F ( x , y , τ, τ ' ) . (2)  \nAccounting for the statistical homogeneity and isotropy, one can write  \nF ( x , y , τ, τ ' ) = F ( | x -y | , τ, τ ' ) . (3)  \nIn terms of the Fourier transformed quantity,  \ns ( k ) = 1 (2 π ) 3 ∫ d x e -i kx · s ( x ) , (4)  \nthe correlation function respecting these properties has the form:  \nF ( k , q , τ, τ ' ) ≡ 1 (2 π ) 6 ∫ d x e -i kx ∫ d y e -i qy · F ( x , y , τ, τ ' ) = ⟨ s ( k , τ ) s ( q , τ ' ) ⟩ √ ⟨ s 2 ( τ ) ⟩ · ⟨ s 2 ( τ ' ) ⟩ = δ ( k + q ) P ( k, τ, τ ' ) , (5)  \nwhere the power spectrum P ( k, τ, τ ' ) is defined as  \nP ( k, τ, τ ' ) = 1 (2 π ) 3 ∫ d z e -i kz F ( | z | , τ, τ ' ) . (6)  \nWe also defined z ≡ x -y . In Eq. (5), we used that ⟨ s 2 ( x , τ ) ⟩ = ⟨ s 2 ( τ ) ⟩ and ⟨ s 2 ( y , τ ' ) ⟩ = ⟨ s 2 ( τ ' ) ⟩ , which follow from statistical homogeneity.  \nFrom Eq. (6), it is straightforward to deduce the behaviour of the function P ( k, τ, τ ' ) in the IR limit k → 0. For this purpose, one observes that the r.h.s. of Eq. (6) remains finite, because the integral there is converging at spatial infinity | x -y | → ∞ by causality. Indeed, correlations of random variables s ( x ) and s ( y ) must be zero, if they are placed in  \ncausally disconnected patches. At equal times, it implies that the separation | x -y | exceeds the horizon. This should also hold true in the limit k → 0. Hence,  \nP ( k, τ, τ ' ) | k → 0 → K ( τ, τ ' ) , (7)  \nwhere K ( τ, τ ' ) < ∞ is a finite function of times τ and τ ' . This will be important when discussing properties of the IR tail of GW spectrum.  \nThe non-linear differential equations may exhibit attractor solutions, so that the asymptotic solution does not depend on the initial conditions, at least on large scales. Moreover, these attractor solutions may reveal self-similarity through dynamic scaling. Assuming such a case for the domain wall network, we restore the time and space dependence of the correlation function on dimensional grounds as  \n⟨ s ( x , τ ) s ( y , τ ) ⟩ ⟨ s 2 ( τ ) ⟩ = F ( | x -y | τ ) . (8)  \nAny dependence on the model dimensionful parameters, that is the domain wall width, is irrelevant once we are interested in the long range dynamics of the domain wall network. Namely, the domain wall width plays the similar role as the UV cutoff scale in effective field theories.  \nNumerical simulations confirm that the domain wall network settles to the scaling regime some time after its formation [21, 22]. One manifestation of the scaling regime is that the number of domain walls in the horizon volume is fixed. Assuming that there is one (a few) domain wall(s) stretching throughout the horizon, we can estimate the energy density of the network as  \nρ dw ∼ ( σ dw · H -2 ) /H -3 ∼ σ dw H , (9)  \nwhere σ dw is the domain wall tension defined as  \nσ dw = 2 √ 2 λη 3 3 . (10)  \nThe domain wall energy density is proportional to the Hubble parameter, while the energy density of any dominant component obeying the Friedmann equation is proportional to the squared Hubble parameter. Hence, if the Hubble parameter decreases with time, domain walls eventually dominate the energy density of the Universe, unless the tension is chosen to be very small or there is a mechanism of domain wall annihilation.  \nFurthermore, the assumption of the scaling regime imposes important constraints on the form of the correlation function P ( k, τ, τ ' ) at equal times, which we write as  \nP ( k, τ, τ ) = τ 3 · P ( kτ ) , (11)  \nwhere we introduced the dimensionless reduced power spectrum P ( kτ ). The expression (11) follows from Eq. (6), where we set τ = τ ' and assume the form of the 2-point correlation  \nfunction F as in Eq. (8). The reduced power spectrum P ( kτ ) is related to the function F by  \nP ( kτ ) = 1 (2 π ) 3 ∫ d u e -i u · k τ F ( u ) , (12)  \nwhere u ≡ z /τ . Note that the expression (11) could be written on dimensionful grounds and from the scaling considerations that there is the unique mass parameter defining the large scale evolution of the domain wall network - Hubble rate, or the inverse horizon size. For the same reasons and recalling the symmetry τ ↔ τ ' , we can immediately write the generic expression for the unequal time correlation function:  \nP ( k, τ, τ ' ) = ( ττ ' ) 3 / 2 · P ( kτ, kτ ' ) . (13)  \nThis consideration has important implications for GWs generated by domain walls, as we show in due course.", '2.2 Energy density of gravitational waves in the scaling regime': "GWs sourced by the domain wall network in the scaling regime have a constant energy density. Previously, this conclusion has been obtained from rough estimates based on the Einstein quadrupole formula [22]. Here we provide a more rigorous analytical derivation of the GW energy density. In what follows, it is assumed that the generation of GWs takes place at radiation domination.  \nEquation of motion describing production of GWs is given by  \n( ∂ 2 ∂τ 2 +2 a ' a ∂ ∂τ -∂ 2 ∂ x 2 ) h ij = 16 πGa 2 ( τ )Π ij ( x , τ ) , (14)  \n(the metric signature is (+ ---)) where G is the Newton's constant, h ij is the transverse traceless (TT) part of the metric perturbation, and  \nΠ ij ( x , τ ) = 1 a 2 ( τ ) Λ ij,kl T kl ( x , τ ) (15)  \nis the TT part of the stress energy tensor; Λ ij,kl is the so called Lambda tensor projecting any symmetric tensor into TT form.  \nWe again switch to Fourier components and perform splitting into polarisations:  \nh ij ( k , τ ) = ∑ A e A ij ( k ) h A ( k , τ ) , (16)  \nand  \nΠ ij ( k , τ ) = ∑ A e A ij ( k )Π A ( k , τ ) , (17)  \nwhere the pair of polarisation tensors e A ij ( k ) obeys the normalisation condition  \ne A ij ( k ) e A ' ij ( k ) = 2 δ AA ' . (18)  \nInitial conditions for Eq. (14) are set at the time τ i , when the domain wall network is formed: h ij = 0 and h ' ij = 0. With these initial conditions, solution of Eq. (14) in the radiationdominated Universe reads  \nh A ( k , τ ) = 16 πG a ( τ ) k ∫ τ τ i dτ ' a 3 ( τ ' ) sin k ( τ -τ ' )Π A ( k , τ ' ) . (19)  \nThe energy density of GWs is given by  \nρ gw ( x , τ ) = 1 32 πGa 2 ( τ ) ⟨ ∂h ij ( x , τ ) ∂τ ∂h ij ( x , τ ) ∂τ ⟩ , (20)  \nwhere averaging is with respect to multiple realisations of a source. In Eq. (20), we expand the perturbation h ij ( x , τ ) into Fourier modes h ij ( k , τ ), then use Eq. (16) and substitute the solution (19). We get  \nρ gw ( τ ) = 8 πG a 4 ( τ ) ∫ d k d q e i ( k + q ) x e A ij ( k ) e A ' ij ( q ) × (21) × ∫ τ τ i ∫ τ τ i dτ ' dτ '' a 3 ( τ ' ) a 3 ( τ '' ) ⟨ Π A ( k , τ ' )Π A ' ( q , τ '' ) ⟩× × [ sin k ( τ -τ ' ) sin q ( τ -τ '' ) k 2 τ 2 -sin( kτ + qτ -kτ ' -qτ '' ) kτ +cos k ( τ -τ ' ) cos q ( τ -τ '' ) ] .  \nFrom the discussion in Sec. 2.1, where we identify s ( k , τ ) ≡ Π A ( k , τ ), one can write for the correlation function on the r.h.s.:  \n⟨ Π A ( k , τ ' )Π A ' ( q , τ '' ) ⟩ = δ AA ' δ ( k + q ) √ ⟨ Π 2 ij ( τ ' ) ⟩ √ ⟨ Π 2 ij ( τ '' ) ⟩ · P ( k, τ ' , τ '' ) , (22)  \nwhere we used Eqs. (2) and (5).  \nFrom Eqs. (21) and (22), the spectral shape of GWs is given by  \ndρ gw d ln k = 128 π 2 Gk 3 a 4 ( τ ) · ∫ τ τ i ∫ τ τ i dτ ' dτ '' a 3 ( τ ' ) a 3 ( τ '' ) √ ⟨ Π 2 ij ( τ ' ) ⟩ √ ⟨ Π 2 ij ( τ '' ) ⟩ · P ( k, τ ' , τ '' ) × × [ sin k ( τ -τ ' ) sin k ( τ -τ '' ) k 2 τ 2 -sin k (2 τ -τ ' -τ '' ) kτ +cos k ( τ -τ ' ) cos k ( τ -τ '' ) ] . (23)  \nThe oscillating term in the square brackets here is independent of k in the limit k → 0. The same is true for the function P ( k, τ ' , τ '' ), as it has been shown earlier, see Eq. (7). We end up with a well-known fact that for kτ ≪ 1, the spectral shape of GWs is fixed to be [48, 49, 50]  \ndρ gw d ln k ∝ k 3 . (24)  \nLet us stress that this behaviour is independent of the assumption of the scaling regime.  \nNow let us use our assumption of the scaling regime (13). Furthermore, we assume that the quantity √ ⟨ Π 2 ij ( τ ) ⟩ evolves with time as the domain wall energy density ρ dw ( τ ) ∼ σ dw H ( τ ), i.e., ⟨ Π 2 ij ( τ ) ⟩ = Cσ 2 dw H 2 ( τ ), and C is constant in time. This can be inferred from the fact that ρ dw = T µν u µ u ν ∼ T ij v i v j = 1 a 2 T ij ( av i )( av j ). From Π ij ∼ 1 a 2 T ij , see Eq. (15), and using that domain walls move with relativistic velocities, i.e., av i ∼ 1, we conclude that Π ij ∼ ρ dw . Using the above assumptions, we obtain from Eq. (21):  \nρ gw ( x , τ ) = 128 π 2 CGσ 2 dw ∫ ∞ 0 dζ ζ 5 ∫ ζ ζ i ∫ ζ ζ i dζ ' dζ '' ( ζ ' · ζ '' ) 5 / 2 · P ( ζ ' , ζ '' ) × × [ sin( ζ -ζ ' ) sin( ζ -ζ '' ) ζ 2 -sin(2 ζ -ζ ' -ζ '' ) ζ +cos( ζ -ζ ' ) cos( ζ -ζ '' ) ] , (25)  \nwhere we defined ζ ≡ kτ and, consequently, ζ ' ≡ kτ ' and ζ '' ≡ kτ '' . It is straightforward to see that the integral on the r.h.s. converges. It depends on time only through the integration limits ζ i = ζτ i /τ . Since ζ i tends to zero in the limit τ ≫ τ i , thus one can set ζ i = 0 and the integral in this limit does not depend on time:  \nρ gw ( τ ) = const τ ≫ τ i . (26)  \nIn view of the lattice simulations, let us note that the integration with respect to ζ in practice does not start from zero, because the comoving momentum k is typically bounded from below by the box size, k ≥ k low > 0. Physically, this means that we neglect GWs with wavelengths longer than the box. Since the comoving box size and hence k low , are constant, the lower integration limit ξ low ≡ k low τ is growing with time. This may induce a slight time dependence, which vanishes as the box size goes to infinity.  \nNote that the (approximate) constancy of the GW energy density explicit in the form (25) holds only as soon as the GW background is supported by the domain wall network. After dissolution of the domain wall network, one naturally recovers the standard radiation-like behaviour ρ gw ( τ ) ∝ 1 /a 4 ( τ ). The spectrum of GWs in terms of the conformal momentum remains intact in this regime.", '3 Setting up lattice simulations': 'Lattice simulations are performed in the cubic box described by N 3 lattice points and a constant comoving volume V = L 3 = const, where L is the comoving box size. The scalar field ϕ evolving in the box is equipped with periodic boundary conditions. This box mimics a portion of the early Universe expanding with the scale factor a (hence, its physical volume grows as V ph ∝ a 3 ) and initially comprising many causally disconnected Hubble patches.  \nSince V ph is growing slower than the Hubble volume ∝ H -3 , the latter eventually exceeds the former. This is the situation to be avoided, because domain walls tend to escape from the box in that case (recall that there is one domain wall in the horizon volume in the scaling regime). Together with a finite lattice spacing L/N , this may pose serious obstacles for the reliable interpretation of our simulations. In Sec. 3.3, we discuss the technique that allows us to circumvent these limitations and to optimise the parameters of simulations in order to achieve trustworthy results. Before doing that we set up the system to make it suitable for lattice simulations by switching to dimensionless variables and fixing initial conditions in Secs. 3.1 and 3.2.', '3.1 Field and coordinate redefinitions': "In what follows, we switch to dimensionless variables defined as  \n˜ τ = √ λη · τ , ˜ x i = √ λη · x i , ˜ k = k √ λη , ˜ ϕ = ϕ η . (27)  \nIn terms of dimensionless conformal time, the cosmological scale factor takes the form:  \na (˜ τ ) = √ λη H i · ˜ τ ˜ τ 2 i , (28)  \nwhere H i is the initial Hubble parameter. It is convenient to set the initial scale factor a (˜ τ i ) to unity in what follows:  \na (˜ τ i ) = 1 . (29)  \nNow, omitting tilde notations, we can write the equation of motion for the field ϕ following from the action (1):  \nϕ '' + 2 a ' a ϕ ' -∂ 2 i ϕ + λη 2 τ 2 H 2 i τ 4 i · ϕ · ( ϕ 2 -1) = 0 , (30)  \nwhere prime refers to the derivative with respect to the dimensionless conformal time. We are free to start the simulations at the time, when the Hubble rate is given by  \nH i = √ λη . (31)  \nThis corresponds to the choice of the initial condition for the Hubble rate. The latter condition is equivalent to setting the initial time of simulations to unity,  \nτ i = 1 , (32)  \nwhere we took into account Eqs. (28) and (29). Using Eqs. (31) and (32), one can simplify the equation of motion for the scalar field:  \nϕ '' + 2 τ ϕ ' -∂ 2 i ϕ + τ 2 · ϕ · ( ϕ 2 -1) = 0 . (33)  \nThe condition (31) is physically justified: while one has H ≫ √ λη , where √ λη is of the order of the field mass, ϕ is almost constant being stuck by the Hubble friction. The field configuration starts to evolve only once the Hubble rate drops down to the field ϕ mass, that explains our choice (31).  \nRemarkably, evolution of the system is independent of the coupling constant λ , and thus λ can be chosen arbitrary. The same concerns the initial Hubble rate H i as well as the expectation value η , as soon as the condition (31) is fulfilled. In particular, we set η = 6 · 10 16 GeV, when performing lattice simulations. Such a large value is rather unconventional in the setup, where domain walls are supposed to be formed during radiation domination. This choice is, however, purely technical, since our results can be easily adjusted to any value of η .", '3.2 Initial conditions for the scalar field': 'One of the goals of this work is to study the effect of initial conditions on evolution of domain wall network. Generically, one writes  \n⟨ ϕ ( k ) ϕ ( q ) ⟩ = A ( k ) δ ( k + q ) (34)  \nand  \n⟨ ˙ ϕ ( k ) ˙ ϕ ( q ) ⟩ = B ( k ) δ ( k + q ) . (35)  \nAs for the first set, we choose vacuum initial conditions with a cutoff k cut applied at high momenta:  \n√  \nA ( k ) = Θ( k -k cut ) 16 π 3 k , B ( k ) = k 16 π 3 · Θ( k -k cut ) . (36)  \nIn what follows, we will mostly assume k cut = λη , or simply k cut = 1 in dimensionless units. Strictly speaking, the initial conditions (36) correspond to the massless scalar field in the Minkowski spacetime. We do not assume that these correlators correctly describe quantum vacuum of the scalar field in our case, rather for the sake of comparison we follow Ref. [22], which was motivated by Ref. [51].  \nFor comparison, for the second set we choose thermal initial conditions:  \nA ( k ) = 1 8 π 3 · k ( e k T -1 ) , B ( k ) = k 8 π 3 · ( e k T -1 ) , (37)  \nwhere T is the temperature of ϕ -particles. We choose T to be equal to the Universe temperature for simplicity. Such initial conditions are somewhat less natural in our setup, because one should equip the field ϕ with a thermal mass in this case 4 . This is, however, irrelevant  \nfor our purposes, because here we aim to check, if the system really approaches the scaling regime, - an attractor solution, which should be reached independently of the starting point. Hence, we are free to pick arbitrary (but clearly distinct) initial conditions.', '3.3 Optimal parameters for lattice simulations': "There are two major challenges when running simulations in the presence of domain walls: i) the wall width δ w = √ 2 /λ · 1 /η must be kept wider than the lattice spacing; ii) one must ensure that the horizon volume does not exceed the simulation box volume. Optimally these two conditions should be violated at about the same time τ f , and then we stop running the code.  \nThe first of the conditions above can be rewritten as  \nδ w = √ 2 λ 1 η = 2 L i N · a ( τ f ) a ( τ i ) , (38)  \nwhere L i is the physical box size at the moment τ i . The factor ' 2 ' on the r.h.s. is introduced such that we are restricting to the cases, when the domain wall width is at least twice larger than the lattice spacing. The second condition reads  \nH -1 ( τ f ) = L i α · a ( τ f ) a ( τ i ) . (39)  \nHere we introduced the 'optimisation' parameter α , which is assumed to be close to unity. Varying this parameter, we will be able to control the level of contamination of the spectrum due to IR and UV obstructions of the problem. Here we explicitly assume that the time τ f is the same as in Eqs. (38) and (39). By virtue of Eq. (38), one can write Eq. (31) as  \nH -1 i = 1 √ λη = √ 2 · L i N · a ( τ f ) a ( τ i ) . (40)  \nDividing Eq. (39) by Eq. (40) and using that the Hubble parameter scales with the conformal time as H ∝ τ -2 at the radiation domination, one gets  \nτ f τ i = 1 2 1 / 4 √ N α . (41)  \nThen Eq. (40) can be rewritten as  \nH -1 i = 1 √ λη = 2 1 / 4 a ( τ i ) L √ αN , (42)  \nwhich is the desired constraint on the set of model parameters. Here we have rewritten L i through L i = a ( τ i ) L ; recall that L is the constant comoving box size. In what follows, for the sake of concreteness we will make the following choices of the parameter α :  \nα = 2 π (basic) , α = 8 π (IR) , α = π 18 (UV) . (43)  \nThe choice α = 2 π is the basic one giving rise to the GW spectrum with the IR and UV parts equally respected. On the other hand, with the choice α = 8 π , one focuses on the IR part of the spectrum somewhat neglecting the UV part, and vice versa in the case α = π/ 18.  \nNote that the condition (41) limits the final time of simulations τ f . In what follows, when discussing GW spectra, we will primarily refer to the lattice with N = 2048, in which case one has τ f /τ i ≃ 15 for α = 2 π .", '4 Results: evolution of the domain wall network': 'Evolution of the domain wall network in the case of vacuum initial conditions is shown in Fig. 2. Thermal initial conditions lead to qualitatively similar pictures, and thus they are not shown. One observes that the wall area inside the simulation box decreases, and the overall network is mainly comprised of a long single wall and a bunch of small closed walls at late times. We compare the domain wall energy density ρ dw and particles ϕ energy density ρ ϕ in Fig. 3. Note that evolution of ρ dw can be fitted by ρ dw ∝ 1 /a 2 ( τ ), which is in accordance with the estimate of ρ dw in the scaling regime ρ dw ∼ σ dw H . Furthermore, the energy density of ϕ -particles ρ ϕ ∝ 1 /a 3 ( τ ) seen in Fig. 3 is in accordance with Ref. [24] for the case of axion models characterised by the domain wall number 5 N dw = 2. Such a behaviour of the particles ϕ energy density is not quite trivial: it suggests that production of particles ϕ is efficient only at early times, while at late times the source of particle production is switched off (see, however, a comment in the end of this section). In particular, this means that closed walls are not efficiently formed during the network evolution, so that scalar radiation resulting from their collapse is negligible (at least in the late time regime). In what follows, we give an additional confirmation of this statement. The characteristic time when the energy density of ϕ -particles starts tracking the behaviour ρ ϕ ∝ 1 /a 3 ( τ ) is τ sc ≃ 4 τ i , which is about the time, when the domain wall network reaches the scaling regime, as we show below.  \nLet us study evolution of the domain wall network more quantitatively. The energy density of the network is given by  \nρ dw ≃ σ dw A a ( t ) V , (44)  \nwhere A is the comoving domain wall area, and V is the constant comoving volume of the simulation box. It is convenient to define the dimensionless scaling parameter ξ [22] also referred to as the area parameter in what follows:  \nξ ≡ At a ( t ) V . (45)  \nFigure 2: Snapshots of domain wall evolution in the case of vacuum initial conditions at different conformal times τ in units of 1 √ λη . Simulations have been performed starting from vacuum initial conditions on a lattice with the grid number N = 512. The visible dot-like structures are small size domain walls.  \n<!-- image -->  \nFigure 3: The energy densities of domain walls and ϕ -particles are shown with orange and black lines, respectively, in units of λη 4 for the cases of vacuum (left) and thermal (right) initial conditions. The total energy density of the system is shown with the blue solid line. Conformal time τ is in units of 1 √ λη . Simulations have been performed on a lattice with the grid number N = 512.  \n<!-- image -->  \nNote that the parameter ξ is constant in the scaling regime. This corresponds to the area A shrinking as 1 /τ in the simulation box V . We used that a ( τ ) ∝ τ and t = ∫ a ( τ ) dτ . Smoothing of the domain wall area seen in Fig. 2 is a reflection of this fact. Conversely, departures of ξ from a constant would signal violation of scaling. Numerically, one obtains the area using the formula [21]:  \nA = ∆ A · ∑ links δ · |∇ ϕ | | ϕ ,x | + | ϕ ,y | + | ϕ ,z | . (46)  \nThe summation here is over the pairs of grid points (links). The area ∆ A = (∆ x ) 2 , where ∆ x is the comoving area between the grid points. The quantity δ takes zero value, if the scalar ϕ has same signs at different grid points, and equals unity, if the scalar ϕ flips its sign.  \nWe show the results of simulations for the parameter ξ in Fig. 4. For this purpose we make use of lattices with the grid numbers N = 512 and N = 1024 assuming vacuum and thermal initial conditions, with and without a cutoff imposed. As the results in cases N = 512 and N = 1024 are in a very good agreement between each other, we deem them trustworthy. The area parameter ξ indeed takes a constant value at sufficiently late times, and hence the domain wall network enters the scaling regime. Interestingly, the onset of scaling occurs universally at the conformal time τ sc ≃ 5 in dimensionless units, or returning to real units, at  \nτ sc ≃ 5 τ i ≃ 5 √ λη . (47)  \nThis can be interpreted in terms of domain wall width δ w relative to the horizon size H -1 ,  \nFigure 4: The area parameter ξ inferred in Eq. (44) is obtained from numerical simulations performed on lattices with the grid numbers N = 512 and N = 1024 starting from vacuum and thermal initial conditions with and without cutoffs at high momenta. Conformal time τ and conformal momentum k are in units of 1 √ λη and √ λη , respectively. The parameter ξ taking a constant value reflects that the domain wall network enters the scaling regime. Expectation values and error bars are obtained from 10 simulations run with different base seed values.  \n<!-- image -->  \ni.e., the scaling starts, when δ w constitutes about ∼ 5% of the Hubble radius H -1 :  \nδ w ≃ 0 . 06 H -1 ( τ sc ) . (48)  \nAt the same time, one observes in Fig. 4 that settling to scaling occurs in a drastically different way in the cases with and without a cutoff. While in the former case, the parameter ξ approaches a constant from below, in the latter case it approaches from above after a period of initial rapid growth. This can be attributed to the fact that initial scalar fluctuations are larger in the no-cutoff case leading to formation of a more curved domain wall initially.  \nNote that the constant values reached by the parameter ξ are clearly distinct in the two cases with thermal and vacuum initial conditions. A similar observation has been made earlier in Ref. [31]. Namely, in the vacuum case, we obtain by fitting to the data  \nξ = 0 . 85 ± 0 . 04 (vacuum) , (49)  \nwhile in the thermal case we get  \nξ = 0 . 63 ± 0 . 04 (thermal) . (50)  \n= 3  \nFigure 5: Histograms showing the distribution of the quantity n · A versus A , where A is the comoving domain wall area in units of 1 λη 2 , and n is the number of domain walls with the area A . Distributions are considered at different conformal times τ in units of 1 √ λη . Simulations have been performed starting from vacuum initial conditions on a lattice with the grid number N = 512.  \n<!-- image -->  \n(The error bars correspond to 1 σ CL.) This result is somewhat surprising, because the definition of scaling itself implies that the domain wall evolution should be independent of initial conditions. The discrepancy is unlikely to be due to differences in the UV parts of the initial scalar spectra. Indeed, analogous spectra with cutoffs imposed lead to the same value of ξ , as it is clearly seen in Fig. 4. We also checked that the discrepancy cannot be attributed to the dependence on the finite box size L . Indeed, the change of L in the case of thermal initial conditions does not lead to considerable changes in ξ , leaving it lower than in the case of vacuum initial conditions. This suggests that the discrepancy is likely to have a physical origin.  \nAs it follows from Fig. 5, the total domain wall area is strongly dominated by a single  \nwall stretching throughout the simulation box, while the area of closed walls is negligible. To the best of our knowledge, it is the first time that this statement is clearly demonstrated with histograms. As a consequence, closed domain walls cannot serve as an efficient source of ϕ -particles. This agrees with the observation from Fig. 3 and our discussion in the beginning of this section that the number density of ϕ -particles is conserved in a comoving volume. At the same time, this appears to be in contrast with the case of cosmic strings, where settling to the scaling regime occurs through formation of loops (analogous to closed walls) evaporating by emission of particles or GWs [53, 54, 55, 56].  \nTo summarise, we confirm a picture of domain wall network that is essentially comprised of a single long wall when the scaling regime is reached [21, 22]. In this regime the selfintersection and formation of small closed structures are rare, as it is indicated in Fig. 5. The physical area of domain walls captured inside the horizon volume exhibits the growth ∝ H -2 independently of initial conditions, - in accordance with the scaling law. At the same time, values of the area parameter shown in Eqs. (49) and (50) are slightly different, suggesting that the scaling solution is sensitive to the choice of initial conditions.  \nIn the end of this section, let us mention one caveat in our discussions. Maintaining the scaling regime generically requires an energy loss of domain walls, which should occur via particle or GW emission. On the other hand, Fig. 3 suggests that efficient particle production takes place only before the scaling is reached. Furthermore, GWs discussed in the next section carry a tiny fraction of the network energy density. With the current resolution, however, it is difficult to estimate significance of the problem. Indeed, it can be of numerical origin, since the code misses wavelengths longer than the box size or particles with de Broglie wavelength shorter than the lattice spacing. We postpone a more detailed study of particle emission from domain walls with C osmo L attice for the future.', '5 Results: gravitational waves': 'In this section we take a closer look at GW production as a result of formation and evolution of a network of domain walls discussed above in Sec. 4. Using the code inbuilt in C osmo L attice we obtain the total energy density and the spectrum of GWs, depending on initial conditions and parameters of simulations. In particular, we check whether the results from analytical approach for the scaling regime presented in Sec. 2.2 agree with the outcome of numerical simulations. We identify characteristic properties of different parts of GW spectra and discuss possible new features. We also pay a particular attention to distinguishing physical effects in the GW spectra from numerical artefacts associated with the finite size of the grid.', '5.1 Energy density of gravitational waves': 'Figure 6: The energy density of GWs in units of λη 4 emitted by the domain wall network is obtained from numerical simulations on lattices with the grid numbers N = 1024 and N = 2048 starting from vacuum and thermal initial conditions with and without cutoffs. Conformal time τ is in units of 1 √ λη . The expectation value η is set at η = 6 · 10 16 GeV. Rescaling to arbitrary η is achieved by multiplying the energy density ρ gw by ( η/ 6 · 10 16 GeV) 2 .  \n<!-- image -->  \nWe show results for the energy density of GWs ρ gw in Fig. 6. One observes that ρ gw takes a constant value, which reflects the fact that the domain wall network enters scaling regime, in accordance with our discussion in Sec. 2. For the curves corresponding to vacuum initial conditions with the cutoff k cut = √ λη and for thermal initial conditions, the energy density ρ gw at large τ takes the values ρ gw ≈ 8 . 1 × 10 -6 λη 4 and ρ gw ≈ 5 . 8 × 10 -6 λη 4 , respectively. Making the latter estimates, we primarily focused on curves in Fig. 6 corresponding to N = 2048. Recalling dependence on model constants and the reduced Planck mass M Pl = 1 / √ 8 πG ≈ 2 . 44 · 10 18 GeV, i.e., ρ dw ∝ σ 2 dw /M 2 Pl (see Eq. (25)), and that we have set η = 6 × 10 16 GeV, one gets  \nρ gw ≈ 0 . 015 · σ 2 dw M 2 Pl (vacuum) , (51)  \nand  \nρ gw ≈ 0 . 011 · σ 2 dw M 2 Pl (thermal) . (52)  \nWe see that the GW energy density in the case of thermal initial conditions is slightly lower compared to the case of vacuum initial conditions. We have already observed the similar tendency in Sec. 4 for the area parameter ξ . Furthermore, as in the case of the area parameter, the energy density of GWs approaches a constant from below, if the initial conditions involve a cutoff, and from above, if the cutoff is absent.  \nFor comparison with the results of Ref. [22], let us quantify the energy density of GWs in terms of the dimensionless quantity:  \nϵ gw = ρ gw Gξ 2 σ 2 dw . (53)  \nWe infer from our numerical simulations:  \nϵ gw ≈ 0 . 52 ± 0 . 05 (vacuum) , (54)  \nand  \nϵ gw ≈ 0 . 70 ± 0 . 09 (thermal) . (55)  \nNote that these values are an order of magnitude lower than those in Ref. [22], see Fig. 6 there. The origin of this discrepancy is unclear to us, taking into account that other results regarding the domain wall evolution, in particular the value (49) for the area parameter, and GW spectra discussed below are in a rather good agreement with Ref. [22].', '5.2 Gravitational wave spectrum': "Results of numerical simulations of the GW spectrum for different initial conditions and different lattice resolution are shown in Figs. 7 and 8. Besides the scaling peak at k ∼ 1 in dimensionless units, one can clearly see the peak in the UV part of the spectrum (see top and bottom left panels). Notice, however, that the height and location of the UV peak depend on the lattice resolution: the value of conformal momentum k associated with the UV peak is twice smaller for N = 1024 as compared to that for N = 2048. Moreover, the UV peak is absent for the UV-optimised lattice (i.e., when the resolution is increased in the UV by the price of losing predictivity in the IR, see the discussion in Sec. 3.3) shown on the bottom right panels of Figs. 7 and 8. This clearly demonstrates that the UV peak is an artefact of numerical simulations due to a finite lattice spacing and it has no physical origin.  \nThe non-physical UV peak is not localised, but it propagates towards lower frequencies of the GW spectrum due to the non-linear dynamics. This may strongly compromise interpretation of the GW spectrum emerging from numerical simulations. To solve this problem, one observes that the level of contamination depends on lattice resolution. Therefore, parts of the spectrum, which do not change when we modify the lattice resolution, are considered to be trustworthy. Hence, we can rely on the results obtained on lattices with N = 1024  \nk  \nFigure 7: Comparison of GW spectra from the domain wall network obtained by making use of numerical simulations performed on lattices with different grid numbers N = 1024 (dashed lines) and N = 2048 (solid lines) starting from vacuum initial conditions with the sharp cutoff set at k cut = 1. Conformal momenta and conformal times are given in terms of √ λη and 1 √ λη , respectively. The expectation value η is set at η = 6 · 10 16 GeV. Rescaling to arbitrary η is achieved by multiplying the spectra by ( η/ 6 · 10 16 GeV) 4 . The spectrum shown on the top panel is calculated for the optimal parameters given by Eq. (42) with the optimisation parameter α = 2 π . The spectra on the bottom left (right) panel refer to optimisation parameter α = 8 π ( α = π/ 18) allowing us to probe deeper the IR (UV) part of the spectrum. Diamonds correspond to the comoving Hubble scale k = 2 πHa , while stars indicate the inverse domain wall width 1 /δ w , i.e., k = 2 πa/δ w .  \n<!-- image -->  \nk  \n3  \n4  \n3  \nFigure 8: Same as in Fig. 7, but with thermal initial conditions.  \n<!-- image -->  \nand N = 2048, provided only that they overlap. Our method to obtain the GW spectrum is then the following. For the vacuum initial conditions we reconstruct central, IR and UV parts of the GW spectrum using different values (43) of the optimisation parameter α entering Eq. (42); the results of reconstruction are shown in the top, bottom left, and bottom right panels of Fig. 7, correspondingly. Then we combine these three plots at wavenumbers k ≈ 0 . 3 and k ≈ 10 in dimensionless units to obtain the full GW spectrum, which does not exhibit the pathological peak in the UV. We repeat the procedure for the thermal initial conditions, combining the central, IR, and UV parts of the GW spectrum shown in the top, bottom left and bottom right panels of Fig. 8, correspondingly, now at wavenumbers k ≈ 0 . 1 and k ≈ 7. Note that such a procedure is legitimate, since the central, IR, and UV parts of the GW spectrum follow from solving the same equation of motion with the same initial conditions in dimensionless units. This is how our main results demonstrated in Fig. 1 have been obtained.  \nFrom Figs. 7 and 8, we estimate the peak frequency of GWs in both cases of vacuum and thermal initial conditions to be  \nF peak = k peak 2 πa ≃ 0 . 7 H . (56)  \nNamely, the peak frequency flows as the Hubble rate, and this flow reflects the domain wall scaling. From the observational point of view, however, only the last instants of GW emission corresponding to the times τ ≃ τ dec are relevant. Here the subscript ' dec ' stands for the moment of time, when the domain wall network collapses. Recall that such a collapse is generically required to avoid overclosing of the Universe. We assume that this moment of time occurs during radiation domination. Taking into account redshifting, the present peak frequency reads:  \nf peak = F peak ( τ dec ) · a dec a 0 ≃ 0 . 7 H dec · a dec a 0 ≃ 7 . 5 nHz ( T dec 100 MeV ) · ( g ∗ ( T dec ) 10 ) 1 / 6 , (57)  \nwhere a 0 is the present day scale factor of the Universe. At the last step in Eq. (57), we have used the normalisation to the decay temperature, which is mostly interesting from the viewpoint of PTAs. Interestingly, the inferred peak frequency f peak is slightly smaller than the one obtained in Ref. [22], where F peak has been found to be even closer to the Hubble rate. This difference is interesting, but we deem it not to be very significant.  \nFrom Figs. 7 and 8 we can also infer the spectral energy density of GWs at the peak:  \nΩ gw,peak ( τ ) ≈ 7 . 7 × 10 -10 · ( H i H ( τ ) ) 2 · ( η 6 · 10 16 GeV ) 4 (vacuum) , (58)  \nand  \nΩ gw,peak ( τ ) ≈ 4 . 6 × 10 -10 · ( H i H ( τ ) ) 2 · ( η 6 · 10 16 GeV ) 4 (thermal) . (59)  \nThese expressions are true during radiation domination prior to the domain wall collapse. Again for the purpose of comparison with Ref. [22], we can express Ω gw,peak as  \nΩ gw,peak ( τ ) = 8 π ˜ ϵ gw ξ 2 G 2 σ 2 dw 3 H 2 ( τ ) . (60)  \nHere following Ref. [22] we introduce the dimensionless quantity ˜ ϵ gw . Comparing Eqs. (58) and (59) with Eq. (60), one gets  \n˜ ϵ gw ≈ 0 . 25 ± 0 . 03 . (61)  \nWe obtain approximately the same value ˜ ϵ gw for both vacuum and thermal initial conditions. Note that the value (61) is in excellent agreement with results of high resolution simulations  \nFigure 9: Same as in Fig. 7 (top panel), but for thermal initial conditions with the sharp cutoff set at k cut = 1.  \n<!-- image -->  \nof Ref. [29], see Eq. (40) there. Furthermore, integrating the spectra over ln k , one obtains the energy density ρ gw and the scaling parameter ϵ gw defined in Eq. (53) consistent with the values in Sec. 5.1. On the other hand, Ref. [30] shows higher GW peak power in the case of thermal initial conditions compared to vacuum ones, contrary to our findings, see Eqs. (58) and (59). We do not know the source of this contradiction; let us note, however, that the value of Ω gw,peak in our case appears to be correlated with the domain wall area, or the scaling parameter ξ . Namely, larger wall area leads to larger GW signal, consistently with naive expectations.  \nRedshifting the gravitons at matter-dominated and dark energy dominated epochs, one can write for the GW relic abundance:  \nΩ gw,peak h 2 0 ≈ 1 . 0 × 10 -10 · ( 100 MeV T dec ) 4 · σ 2 dw (100 TeV) 6 · ( 10 g ∗ ( T dec ) ) 4 / 3 (vacuum) , (62)  \nand  \nΩ gw,peak h 2 0 ≈ 0 . 6 × 10 -10 · ( 100 MeV T dec ) 4 · σ 2 dw (100 TeV) 6 · ( 10 g ∗ ( T dec ) ) 4 / 3 (thermal) . (63)  \nHere we again adopt the normalisation parameters, which are most suitable from the viewpoint of PTAs.  \nThe IR part of the GW spectrum behaves in accordance with causality arguments, i.e.,  \nΩ gw h 2 0 ∝ k n k ≪ k peak , (64)  \nwhere the exponent n ≃ 3, cf. Eq. (24). It is important to stress though that our simulations reveal small, but significant, departures from the exponent n = 3 suggested by causality. For example, by fitting to the numerical results we obtain n = 2 . 66 ± 0 . 06 for vacuum initial conditions and n = 2 . 69 ± 0 . 11 for thermal initial conditions, respectively, at the moment τ = 9. We attribute this difference to the presence of the low limit for momenta k in the simulations due to a finite box size, while the case n = 3 is suggested to be recovered in the limit k → 0. Indeed, at the time τ = 5, one approaches closer to n = 3 getting from the fit to simulations n = 2 . 74 ± 0 . 04 and n = 2 . 87 ± 0 . 06 for vacuum and thermal initial conditions, correspondingly.  \nTo the right of the peak, the slope of GWs taken at the time τ = 13 behaves as (according to the fit to our numerical results)  \nΩ gw h 2 0 ∝ k -1 . 53 ± 0 . 04 k peak < k < k ∗ (vacuum) , (65)  \nand  \nΩ gw h 2 0 ∝ k -1 . 28 ± 0 . 01 k peak < k < k ∗ (thermal) . (66)  \nThis decrease is sharper than in Refs. [22, 29] giving -1 for the exponent, but milder than that of Ref. [30] giving -1 . 7. In Eqs. (65) and (66), we have introduced a new scale k ∗ , at which the falloff of the spectrum is followed by a (near) plateau-type of behaviour, i.e.,  \nΩ gw h 2 0 ≃ const k ∗ < k < 2 πa/δ w . (67)  \nNote that the plateau region is not perfectly flat, especially in the case of vacuum initial conditions, when there is a bump at intermediate scales. The same observation has been made in Ref. [30]. It appears, however, that the bump is a consequence of the sharp cutoff on the vacuum spectrum. Indeed, the bump also arises in the case of thermal initial conditions with the cutoff applied, as it follows from Fig. 9. On the other hand, we see no clear relation between bump locations on the spectra and cutoff scales; therefore the mechanism of a bump creation is not obvious to us. Finally, there is the exponential falloff of the spectrum for the momenta exceeding the inverse domain wall width/mass of scalar, i.e., at k > 2 πa/δ w , which has a clear physical origin.  \nLet us comment on the possible origin of the scale k ∗ . As it can be seen in Figs. 7 and 8, the location of k ∗ approximately stays constant with time, and it is not sensitive to the choice of initial conditions. That is, the quantity k ∗ is likely to be related to the Hubble rate and the inverse wall width by k ∗ ≃ 2 π √ H/δ w · a . This, nevertheless, remains a conjecture, which may be worth investigating in more details in the future. If this hypothesis is correct, then  \nthe plateau (bump) is moving deep into the UV region relative to the main peak at large times. Its width in terms of the ratio of the largest and smallest momenta bounding the plateau is growing with time as √ t/δ w . Consequently, the longer the domain wall lifespan is, the better the prospects of observing the plateau in future GW experiments.", '6 Discussion': 'In this work we revisited evolution of the domain wall network and the generated GW spectrum using the C osmo L attice code. Since C osmo L attice is publicly available, we believe that the results of this work will pave the way for a broader audience aiming at the analysis of domain walls. The most non-trivial part of such analysis is separating physical effects from artefacts triggered by a limited grid number N and a finite box size. In particular, the artefact caused by finite lattice spacing is manifested as a peak in the UV part of the spectrum. Due to non-linearity of equations, the contamination from the UV peak propagates towards lower momenta with time. Consequently, the longer the time of simulations, the smaller part of the spectrum is available for physical interpretation. We developed a technique allowing to discard nonphysical artefacts and extract trustworthy information from simulations.  \nAfter excluding non-physical artefacts we find that our results summarised in Fig. 1 are in a good agreement with those of Refs. [29, 30]. It is important to mention that unlike Ref. [29] we did not invoke a potential bias or any other mechanism of domain wall destruction in our simulations. This suggests that the potential bias has little effect on the shape of GW spectrum. However, this statement should be taken with a grain of salt, as some properties of GWs are not yet well established. For example, compared to Refs. [22, 29, 30], we get slightly different values for the power-law exponent on the r.h.s. of the (scaling) peak, see Eqs. (65) and (66). Notably, our values are in between -1 and -1 . 7 that were obtained in Refs. [22, 29] and Ref. [30], respectively.  \nThere are some interesting features in the domain wall evolution and GW spectrum revealed after excluding the artefacts, which deserve a closer look. In particular, the area parameter ξ defined in Eq. (4) takes different values in the cases of thermal and vacuum initial conditions, although it is believed to have same values due to the scaling regime, cf. Ref. [31]. We have checked that the difference is unlikely to be attributed to the lattice spacing or the box size. Hence, it may have a physical origin, which is yet to be understood. Furthermore, the spectral energy density of GWs at peak also turns out to be different, which can be related to the aforementioned differences in the parameter ξ evolution. Finally, the origin of the plateau region (bumpy region) at intermediate scales in the GW spectrum is unclear. It appears that the region starts at the comoving scale k ∗ ≃ 2 π √ H/δ w · a , but the physical origin of such a scale is lacking.  \nThe presence of the plateau can have interesting consequences for GW searches. In  \nparticular, assuming a suitable formation time and a lifespan of the domain wall network, one can have a sizeable GW signal from domain walls in the form of the plateau in the LISA frequency range, which is also predicted from cosmic strings and inflationary models. This increases the level of degeneracy between different GW sources, which can be broken, e.g., by GW observations in the nHz frequency range covered by PTAs. Further discussion of potential phenomenological applications of the GW spectra shown in Fig. 1 is out of the scope of this work, and we plan to investigate them in more details elsewhere.  \nIn the future, we are also planning to extend our analysis to include domain walls described by a time-varying tension, with a particular focus on melting domain walls [18, 19, 20]. The principal advantage of this scenario in comparison to the standard one is that the decreasing domain wall tension solves the problem of overclosing the Universe [18, 19, 20]. Furthermore, the time dependence of the wall tension follows from well motivated particle physics scenarios [18]. In this situation, one expects a drastically different GW spectrum compared to that of this work. Importantly, the predicted power-law shape of the GW spectrum Ω gw ∝ k 2 [19, 20, 57] is favoured by the recent PTA observations. It is therefore crucial to perform numerical simulations, similar to those presented in this paper, to support this prediction of GW spectrum in the case of melting domain walls.', 'Acknowledgments': "Numerical simulations were performed on the computer cluster of the Theoretical Division of INR RAS and on the 'Phoebe' computer cluster at CEICO/FZU. We are indebted to Josef Dvoˇr'aˇcek for his assistance in using the latter cluster and to Daniel Figueroa for useful discussions. I. D. and D. G. acknowledge the support of the scientific program of the National Center for Physics and Mathematics, section 5 'Particle Physics and Cosmology', stage 2023-2025. E. B. acknowledges the support of ANR grant StronG (ANR-22-CE31-0015-01). A. V. was supported by the Czech Science Foundation (GA ˇ CR), project 24-13079S.", 'References': "- [1] G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 951 (2023) no.1, L8; [arXiv:2306.16213 [astro-ph.HE]].\n- [2] J. Antoniadis et al. [EPTA and InPTA:], Astron. Astrophys. 678 (2023), A50; [arXiv:2306.16214 [astro-ph.HE]].\n- [3] D. J. Reardon, A. Zic, R. M. Shannon, G. B. Hobbs, M. Bailes, V. Di Marco, A. Kapur, A. F. Rogers, E. Thrane and J. Askew, et al. Astrophys. J. Lett. 951 (2023) no.1, L6; [arXiv:2306.16215 [astro-ph.HE]].  \n- [4] H. Xu, S. Chen, Y. Guo, J. Jiang, B. Wang, J. Xu, Z. Xue, R. N. Caballero, J. Yuan and Y. Xu, et al. Res. Astron. Astrophys. 23 (2023) no.7, 075024; [arXiv:2306.16216 [astro-ph.HE]].\n- [5] A. Afzal et al. [NANOGrav], Astrophys. J. Lett. 951 (2023) no.1, L11; [arXiv:2306.16219 [astro-ph.HE]].\n- [6] G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 952 (2023) no.2, L37; [arXiv:2306.16220 [astro-ph.HE]].\n- [7] J. Antoniadis et al. [EPTA], [arXiv:2306.16227 [astro-ph.CO]].\n- [8] A. Sesana, Class. Quant. Grav. 30 (2013), 224014; [arXiv:1307.2600 [astro-ph.CO]].\n- [9] J. Ellis, M. Fairbairn, G. Franciolini, G. Hutsi, A. Iovino, M. Lewicki, M. Raidal, J. Urrutia, V. Vaskonen and H. Veermae, Phys. Rev. D 109 (2024) no.2, 023522 [arXiv:2308.08546 [astro-ph.CO]].\n- [10] G. Sato-Polito, M. Zaldarriaga and E. Quataert, [arXiv:2312.06756 [astro-ph.CO]].\n- [11] P. Amaro-Seoane et al. [LISA], [arXiv:1702.00786 [astro-ph.IM]].\n- [12] S. Hild et al. Class. Quant. Grav. 28 (2011), 094013; [arXiv:1012.0908 [gr-qc]].\n- [13] Y. B. Zeldovich, I. Y. Kobzarev and L. B. Okun, Zh. Eksp. Teor. Fiz. 67 (1974), 3-11\n- [14] A. Lazanu, C. J. A. P. Martins and E. P. S. Shellard, Phys. Lett. B 747 (2015), 426-432; [arXiv:1505.03673 [astro-ph.CO]].\n- [15] G. B. Gelmini, M. Gleiser and E. W. Kolb, Phys. Rev. D 39 (1989), 1558;\n- [16] D. Coulson, Z. Lalak and B. A. Ovrut, Phys. Rev. D 53 (1996), 4237-4246;\n- [17] A. Vilenkin, Phys. Rev. D 23 (1981), 852-857;\n- [18] S. Ramazanov, E. Babichev, D. Gorbunov and A. Vikman, Phys. Rev. D 105 (2022) no.6, 063530; [arXiv:2104.13722 [hep-ph]].\n- [19] E. Babichev, D. Gorbunov, S. Ramazanov and A. Vikman, JCAP 04 (2022) no.04, 028; [arXiv:2112.12608 [hep-ph]].\n- [20] E. Babichev, D. Gorbunov, S. Ramazanov, R. Samanta and A. Vikman, Phys. Rev. D 108 (2023) no.12, 123529; [arXiv:2307.04582 [hep-ph]].\n- [21] W. H. Press, B. S. Ryden and D. N. Spergel, Astrophys. J. 347 (1989), 590-604.  \n- [22] T. Hiramatsu, M. Kawasaki and K. Saikawa, JCAP 02 (2014), 031; [arXiv:1309.5001 [astro-ph.CO]].\n- [23] K. Saikawa, Universe 3 (2017) no.2, 40; [arXiv:1703.02576 [hep-ph]].\n- [24] T. Hiramatsu, M. Kawasaki, K. Saikawa and T. Sekiguchi, JCAP 01 (2013), 001; [arXiv:1207.3166 [hep-ph]].\n- [25] T. Hiramatsu, M. Kawasaki, K. Saikawa and T. Sekiguchi, Phys. Rev. D 85 (2012), 105020 [erratum: Phys. Rev. D 86 (2012), 089902]; [arXiv:1202.5851 [hep-ph]].\n- [26] M. Kawasaki, K. Saikawa and T. Sekiguchi, Phys. Rev. D 91 (2015) no.6, 065014; [arXiv:1412.0789 [hep-ph]].\n- [27] C. F. Chang and Y. Cui, [arXiv:2309.15920 [hep-ph]].\n- [28] T. Krajewski, J. H. Kwapisz, Z. Lalak and M. Lewicki, Phys. Rev. D 104 (2021) no.12, 123522; [arXiv:2103.03225 [astro-ph.CO]].\n- [29] N. Kitajima, J. Lee, K. Murai, F. Takahashi and W. Yin, [arXiv:2306.17146 [hep-ph]].\n- [30] R. Z. Ferreira, S. Gasparotto, T. Hiramatsu, I. Obata and O. Pujolas, [arXiv:2312.14104 [hep-ph]].\n- [31] N. Kitajima, J. Lee, F. Takahashi and W. Yin, [arXiv:2311.14590 [hep-ph]].\n- [32] R. Z. Ferreira, A. Notari, O. Pujol'as and F. Rompineve, [arXiv:2401.14331 [astroph.CO]].\n- [33] D. I. Dunsky and M. Kongsore, [arXiv:2402.03426 [hep-ph]].\n- [34] Y. Gouttenoire and E. Vitagliano, Phys. Rev. D 109 (2024) no.12, 123507; [arXiv:2311.07670 [hep-ph]].\n- [35] M. Hindmarsh, Phys. Rev. Lett. 77 (1996), 4495-4498; [arXiv:hep-ph/9605332 [hep-ph]].\n- [36] P. P. Avelino, C. J. A. P. Martins and J. C. R. E. Oliveira, Phys. Rev. D 72 (2005), 083506; [arXiv:hep-ph/0507272 [hep-ph]].\n- [37] C. J. A. P. Martins, I. Y. Rybak, A. Avgoustidis and E. P. S. Shellard, Phys. Rev. D 93 (2016) no.4, 043534; [arXiv:1602.01322 [hep-ph]].\n- [38] C. J. A. P. Martins, I. Y. Rybak, A. Avgoustidis and E. P. S. Shellard, Phys. Rev. D 94 (2016) no.11, 116017 [erratum: Phys. Rev. D 95 (2017) no.3, 039902]; [arXiv:1612.08863 [hep-ph]].  \n- [39] P. P. Avelino, JCAP 04 (2020), 012; [arXiv:2001.06318 [astro-ph.CO]].\n- [40] O. Pujolas and G. Zahariade, Phys. Rev. D 107 (2023) no.12, 123527; [arXiv:2212.11204 [hep-th]].\n- [41] D. Gruber, L. Sousa and P. P. Avelino, [arXiv:2403.09816 [gr-qc]].\n- [42] D. G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, JCAP 04 (2021), 035; [arXiv:2006.15122 [astro-ph.CO]].\n- [43] D. G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, Comput. Phys. Commun. 283 (2023), 108586; [arXiv:2102.01031 [astro-ph.CO]].\n- [44] D. G. Figueroa, A. Florio and F. Torrenti, [arXiv:2312.15056 [astro-ph.CO]].\n- [45] T. Garagounis and M. Hindmarsh, Phys. Rev. D 68 (2003), 103506; [arXiv:hepph/0212359 [hep-ph]].\n- [46] Y. Li, L. Bian and Y. Jia, [arXiv:2304.05220 [hep-ph]].\n- [47] Y. Li, L. Bian, R. G. Cai and J. Shu, [arXiv:2311.02011 [astro-ph.CO]].\n- [48] C. Caprini, R. Durrer, T. Konstandin and G. Servant, Phys. Rev. D 79 (2009), 083519; [arXiv:0901.1661 [astro-ph.CO]].\n- [49] R. G. Cai, S. Pi and M. Sasaki, Phys. Rev. D 102 (2020) no.8, 083528; [arXiv:1909.13728 [astro-ph.CO]].\n- [50] A. Hook, G. Marques-Tavares and D. Racco, JHEP 02 (2021), 117; [arXiv:2010.03568 [hep-ph]].\n- [51] G. N. Felder, L. Kofman and A. D. Linde, Phys. Rev. D 64 (2001), 123517 [arXiv:hepth/0106179 [hep-th]].\n- [52] S. Blasi, A. Mariotti, A. Rase, A. Sevrin and K. Turbang, JCAP 04 (2023), 008; [arXiv:2210.14246 [hep-ph]].\n- [53] T. Vachaspati and A. Vilenkin, Phys. Rev. D 31 (1985), 3052;\n- [54] A. Vilenkin and E. P. S. Shellard, Cambridge University Press, 1994.\n- [55] M. B. Hindmarsh and T. W. B. Kibble, Rept. Prog. Phys. 58 (1995), 477-562; [arXiv:hep-ph/9411342 [hep-ph]].  \n- [56] J. J. Blanco-Pillado, K. D. Olum and B. Shlaer, Phys. Rev. D 89 (2014) no.2, 023512; [arXiv:1309.6637 [astro-ph.CO]].\n- [57] S. Ramazanov, JCAP 02 (2024), 027; [arXiv:2310.19148 [gr-qc]]."}