{'fraction_non_alphanumeric': 0.03903045369794904, 'fraction_numerical': 0.015537600994406464, 'mean_word_length': 4.438485804416404, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability of domain wall network with initial inflationary fluctuations, and its implications for cosmic birefringence

We study the formation and evolution of domain walls with initial inflationary fluctuations by numerical lattice calculations that, for the first time, correctly take into account correlations on superhorizon scales. We find that, contrary to the widely-held claim {over the past few tens of years}, the domain wall network exhibits remarkable stability even when the initial distribution is largely biased toward one of the minima. This is due to the fact that the domain wall network retains information about initial conditions on superhorizon scales, and that the scaling solution is not a local attractor in this sense. Our finding immediately implies that such domain walls will have a significant impact on cosmology, including the production of gravitational waves, baryogenesis, and dark matter from domain walls. Applying this result to the axion-like particle domain wall, we show that it not only explains the isotropic cosmic birefringence suggested by the recent analysis, but also predicts anisotropic cosmic birefringence that is nearly scale-invariant on large scales and can be probed by future CMB observations.

Introduction. -Topological defects are formed by spontaneous symmetry breaking and have been widely studied. [1-3] (see also e.g. [4] for a review). While the stability of individual topological defects is guaranteed by the conservation of topological charge, the evolution of the large number of topological defects formed in the early universe is not obvious because it depends on their evolution and interactions. Such topological defects and their traces may be explored experimentally in the present universe, or in some cases they may pose serious cosmological problems.
We focus on the cosmological evolution of the simplest class of topological defects, domain walls (DWs). DWs are field configurations connecting distinct degenerate vacua and are formed from spontaneous breaking of discrete symmetries. The simplest DWs are the ones associated with the Z 2 symmetric potential, Here φ is the scalar field, λ(> 0) is the quartic coupling, and two degenerate vacua are located at φ min = ±v. Once such DWs separated by the two vacua are formed in the early universe, it is known that the subsequent evolution converges to the so-called scaling solution, a situation where there is on average about one DW per the Hubble horizon [5][6][7][8][9][10][11][12].
The energy density of DWs that obey the scaling solution decreases more slowly than radiation or matter, and thus easily dominates the universe, making it significantly inhomogeneous. Even if it were subdominant, it could be inconsistent with current observations by making the cosmic microwave background (CMB) spectrum largely anisotropic. This is the notorious cosmological DW problem [1,13]. This problem can be avoided if the DW tension is small enough, or if DWs are unstable. A well-known way to destabilize DWs is to introduce a symmetry-breaking bias in the potential and/or the initial distribution of the scalar field [3,[14][15][16][17][18][19][20][21][22][23][24][25].
Our focus is on the case in which the initial distribution is biased. Such a biased initial distribution can be realized when fluctuations are generated by non-thermal mechanisms, for example, fluctuations generated during inflation [26][27][28]. It was argued in Refs. [17][18][19][20][21][22] that since the DW is a nonlinear object and the field distribution is highly non-Gaussian, the DW network with initial inflationary fluctuations can be understood based on percolation theory.
Under this assumption, previous numerical calculations showed that DWs decay quickly even when the bias in the initial distributions is very small [18][19][20][21][22][23][24]. In this Letter we focus on the DWs generated by inflationary fluctuations and show that, contrary to this widely held conclusion, the DWs are very stable against biased initial distributions.
The purpose of this Letter is two-fold. First, we study the formation and evolution of DWs with the initial inflationary fluctuations by numerical lattice simulations, correctly taking into account correlations on superhorizon scales. We then show that the resultant DW networks are very stable against biased initial fluctuations. This stability is due to the DW structure on large scales, which reflects the initial fluctuations at scales beyond the Hubble horizon.
Second, as a phenomenological application of this result, we consider axion-like particle (ALP) DW with inflationary fluctuations. In Ref. [29] two of the present authors (FT and WY) proposed the ALP DW to explain a hint for the isotropic cosmic birefringence (CB) [30,31], β = 0.36 ± 0.11 deg, where β denotes the isotropic rotation angle of the CMB polarization. Interestingly, the ALP DWs naturally account for the closeness of the rotation angle to the fine structure constant, α = 1/137 rad 0.42 deg, over a very wide range of the ALP mass and decay constant. Furthermore, they predict a characteristic anisotropic CB determined by the DW distribu-tion on the last scattering surface (LSS), but the case with the initial correlation on superhorizon scales has not been studied. In this Letter, we investigate for the first time the anisotropic CB power spectrum predicted by ALP DWs with inflationary fluctuations and show that it can be tested in future CMB observations such as Lite-BIRD [32], Simons Observatory [33], CMB-S4-like [34], and PICO [35].
Here we mention some relevant past literature. The axion DW with inflationary fluctuations was first considered in [26,27] in the context of the QCD axion, and since then it has been discussed in various context [29,[36][37][38][39][40][41][42][43]. In the symmetry breaking scenario during inflation, superhorizon-scale fluctuations are generated over a certain range of scales, and their spatial distribution was studied in Ref. [28]. See also Refs. [29,40,43] for discussions on the relation between the superhorizon fluctuations and the DW formation. However, to the best of our knowledge, detailed numerical lattice simulations on the formation and evolution of DWs that properly incorporate long-wavelength correlations have not been performed.
Setup. -Let us consider a scalar field φ with the quartic potential (1) as a simple example. The scalar mass at the potential minimum is m φ = √ 2m 0 . A planar DW perpendicular to the z-axis is given by φ DW = v tanh (zm φ /2). Thus the DW width is of m −1 φ , and the tension σ is 4 √ λv 3 /3. As long as we focus on the macroscopic behavior of DW network at scales much larger than the DW width, the following discussion is also applicable to other potential forms, such as the sine-Gordon potential, which will be discussed later.
The equation of motion for the minimally coupled φ in the flat Friedmann-Lemaître-Robertson-Walker universe with three spatial dimensions is where t is the cosmic time, x is the comoving coordinate, a(t) is the scale factor, H = (da/dt)/a is the Hubble parameter, and ∂ 2 i is the Laplace operator. Here we neglect the gravitational self-interaction of DWs, which is justified if the DW energy density is subdominant.
Clarification of initial fluctuations. -The initial conditions of the scalar field are of essential importance for the formation and subsequent evolution of DW, which is one of the findings of this Letter. We assume that the scalar field initially fluctuates around φ 0 , where φ 0 is the field value averaged over the spatial region of our interest. We can bias the initial fluctuations by making φ 0 nonzero. For simplicity, let us first assume φ 0 = 0 and consider fluctuations around the origin.
The scale-dependence of the fluctuations is encoded in the power spectrum P (k) defined by where is the Fourier transform, the angle brackets mean an ensemble average, and we introduced the spatial dimensions d = 2 or 3 for later use. Here, the rotational invariance is assumed to use k = |k| for the power spectrum. Then, the variance of fluctuations in real space can be expressed as where P(k) is the reduced power spectrum defined by P(k) = k 2 P (k)/2π for d = 2 and P(k) = k 3 P (k)/2π 2 for d = 3. We mainly consider inflationary fluctuations; if the scalar is sufficiently light, i.e. m φ H inf , it acquires Gaussian fluctuations with almost scale-invariant power spectrum P(k) (H inf /2π) 2 during inflation, where H inf is the inflationary Hubble parameter.
The spatial average φ (= φ 0 ) is generically non-zero. In this case the initial distribution is biased, and it is obtained simply by shifting the above distribution by φ 0 . To quantify the bias, we define the bias parameter as For the inflationary fluctuations, we replace the denominator with P(k) to avoid the logarithmic dependence on the lattice box size. For the numerical lattice calculations, the initial conditions for the inflationary fluctuations are prepared as follows. First, we generate scalar fluctuations in the momentum space so that they satisfy the scale-invariant power spectrum. Then, we set the initial scalar field fluctuations in the real space by the inverse Fourier transform. On the other hand, for the Gaussian white noise, as in the case of thermal fluctuations, we randomly generate a scalar field value at each spatial point, using the Gaussian probability distribution.
Scaling solution. -For the white noise or thermal initial conditions, the DW network quickly converges to the scaling solution, where there is on average about one DW per the Hubble horizon [5][6][7][8][9][10][11][12], given by with L DW being the average DW length (area) per Hubble horizon for d = 2(3). We will see that, while the evolution of the DW length is given by Eq. (6) for inflationary fluctuations, the spatial DW distribution is significantly different from that for the white noise initial condition.
Stability of the DWs against the initial condition bias. -We solve numerically the equation of motion (2) on a lattice in the comoving coordinate space with two different initial conditions: inflationary fluctuations and Gaussian white noise. In the case of the inflationary fluctuations, in particular, the lattice box size should be sufficiently large because of the correlations on large scales. To this end, we focus on DWs with their motion restricted in 2D space. We have performed numerical simulations in 3D space to confirm that similar results are obtained both quantitatively and qualitatively. We have set the box size to contain 200 2 Hubble horizons at the initial conformal time τ i = 1/m 0 on a lattice with a grid size of 4096 × 4096 and a periodic boundary condition. Here we assume the radiation dominated universe, and define the conformal time τ = 2t/m 0 for convenience. Then, we have solved the equation of motion (2) using the second order leapfrog method with the two initial conditions, varying the bias b d of the initial distributions.
The snapshots of the scalar field at τ = 10/m 0 in the case of b d = 0 are shown in Fig. 1, where the DWs separate the two vacua, φ = ±v. Compared to the white noise case, there are very large regions in the case of inflationary fluctuations. Such a region much larger than the Hubble horizon is considered to be stable, at least until the horizon grows to its size, due to the causality. Thus, in the case of scale-invariant fluctuations, the DW network is expected to be very stable because such a large region exists at any scale.
To see the stability of DWs with biased initial fluctuations, we show in Fig. 2 the time evolution of the averaged DW length per Hubble horizon. In each panel the top line corresponds to b d = 0 which satisfies Eq. (6) 1 , while the lower lines are for more biased cases. In the case of white noise, the DWs collapse in a very short timescale even with a very small bias b d 1, whereas in the case of inflationary fluctuations, DWs are very long-lived even with a very large bias b d = O(1). The reason for the stability is the correlation of the scalar fluctuations at large scales. This leads to the characteristic void structure of domains, which cannot collapse until the Hubble horizon becomes comparable to its size. This is the central result of this Letter, which differs significantly from the claims in the literature [17][18][19][20][21][22].
Lastly, we show in Fig. 3 the reduced power spectrum 1 The slight increase of L DW /H −1 is mainly due to the finite box effect, which becomes milder as the lattice size increases.  Here k H is the comoving wavenumber such that its corresponding wavelength is equal to the Hubble horizon, k H /a = 2πH. The red and blue lines correspond to the initial conditions of inflationary fluctuations and white noise [43], respectively. Each light-colored band represents the 1σ statistical uncertainties based on the 50 times simulations. One can see that the power spectrum is peaked around k = k H and there is no correlation at large scales in the case of white noise initial condition.
The suppression at small scales k k H is due to the DW dynamics. In fact, these features can be well explained by a simple model in which infinitely thin DWs have a uniform probability distribution on an arbitrary 1D line in space [29,43]. In the case of inflationary fluctuations, on the other hand, the power spectrum is nearly scaleinvariant at large scales, and suppressed at small scales as in the white noise case. The slight tilt of the power spectrum at large scales can be understood as follows. First, the initial power spectrum is scale-invariant and grows as the scalar field rolls down to either potential minimum. Then, the enhancement factor for a mode k  [44] and the sensitivity reach of LiteBIRD for the scale-invariant spectrum, respectively. They are adopted from Ref. [45] (see also Refs. [32][33][34][35]).
is inversely proportional to the typical initial value of the coarse-grained field,φ R , averaged over the comoving scale R = 2π/k, andφ R grows as R increases; where R 1 > R 2 , i.e. k 1 < k 2 , is assumed in the last inequality. We confirmed that this nicely explains the slightly blue-tilted spectrum at large scales. An application of the DW stability to the cosmic birefringence from ALP DW. -The dynamics of the DW network is not sensitive to the detailed shape of the Z 2 symmetric potential, and in particular, the above argument can be applied to the sine-Gordon potential, where m φ and f φ are the mass and decay constant of φ. In the following φ is identified with an ALP with the above potential. The sine-Gordon potential is the simplest of the potentials with discrete shift symmetry, φ → φ + 2πf φ , and there are mulitple degenerate vacua. We assume that the two vacua at φ min = ±πf φ are populated by the initial fluctuations, for simplicity. The ALP is coupled to photons as where c γ is the anomaly coefficient of O(1), α em the finestructure constant, and F µν andF µν the photon field strength and its dual, respectively. Suppose that during the propagation of the CMB photon from the LSS to us, the ALP field value changes in time and space slowly enough compared to the frequency of the photon. Then the polarization plane rotates in proportion to the change in the ALP field [46][47][48]. If the ALP DWs have existed from the time of recombination to the present, the rotation angle is [29] with Ω being the spherical polar coordinates. Here φ today and φ LSS is the vacuum at the solar system and at the LSS, respectively. The parenthesis in Eq. (10) is approximately equal to 1/2 when averaged over the entire sky 2 , and it predicts the isotropic CB, β 0.21c γ deg [29], independent of m φ and f φ . This is consistent with the recent hint, β = 0.36 ± 0.11 deg [30,31].
Unlike models based on spatially uniform ALP [30,[49][50][51][52][53][54][55], the ALP DWs predict anisotropic CB that reflects the DW distribution on the LSS. We can map the reduced power spectrum obtained above to the LSS by using the relations, where d L is the comoving distance to the LSS, and C Φ is the angular power spectrum of ∆Φ(Ω) [43,56]. The mapped and C Φ are shown on the upper and right axes of the panel in Fig. 3, where β = 0.36 deg is assumed. Note that the anisotropic CB is correlated with the isotropic one as C Φ ∝ β 2 . We also show the future reaches and present bound for the white noise and scale-invariant power spectra, including the BICEP/Keck bound on the scale-invariant spectrum, ( + 1)C Φ /(2π) < 0.014 deg 2 (95 % CL) [44]. While the prediction shown in the figure barely satisfies the current limit, one can relax the tension by taking slightly smaller β or including the effects of cosmic variance and/or bias, while explaining the isotropic CB [30,31]. We have checked that, even for b d = 1, the reduced power spectrum gets suppressed up to by 20% at large scales and is stable over time. Note also that the lowest multipoles < 20 are most important for constraining the scale-invariant CB. Thus, the predicted anisotropic CB can be confirmed or refuted in near-future experiments including LiteBIRD.
Conclusions.-We have shown that DW network is very stable against the bias of initial inflationary distribution. This is at odds with the understanding of the past few tens of years. The unexpected stability is due to correlations on large scales generaterd by inflation, which lead to the formation of DWs of superhorizon size.
Once formed, such DWs are very stable due to causality in the present set-up. Thus, various studies on the production of the gravitational waves, baryon asymmetry, dark matter, and CB from DWs must be reconsidered. We studied the last one in this Letter, but there is still room for further research on the rest.

Time evolution of the DW network
In the main text we have considered the characterization of the DW network through the total DW length per Hubble horizon, L DW /H −1 . There are other interesting measures in the system, like the total positive or negative area. We define the ratio r as where S + (S − ) is the total area where field values of φ are positive (negative). Here r is defined so that it is bounded as r ∈ [0, 2], making it easy to study in a logarithmic plot. It is straightforward to define the ratio for the 3D case. In a d dimensional space, L DW /H −(d−1) and r encode information on the d − 1 dimensional hypersurface and d dimensional volume, respectively. In principle, these are independent measures, all their correlations stemming from the DW dynamics. We show in Fig. 4 the time evolution of the ratio r for the inflationary fluctuations (top) and Gaussian white noise (bottom). As with the DW length, the ratio decreases only slowly even for large biases in the case of inflationary fluctuations. Also, as we can see in Fig. 4, after m 0 τ = 5 when the DW network is fully formed, the evolution of r can be well described by an exponential decay over the time range shown. However, unless b d is very large, the DW lifetime in the case of inflationary fluctuations is so long that it is not possible to determine numerically whether r eventually approaches 0 and the DW network decays completely, or whether r asymptotes to a certain value.
In Fig. 5 we show snapshots of the DW network for various biases in the case of inflationary fluctuations. For b d = 0.15 -a very large bias in the case of white noise DW networks-we can see that even at later times there are still "infinitely" large DWs within the box. Once the bias becomes large enough, as is in the case of b d = 1.25, then no "infinitely" large DWs are found at later times. This suggests that for the inflationary fluctuations, if the fluctuation at the largest scale of interest is greater than the bias, DWs separating domains of the largest size are likely formed. Such DWs will be stable or very long-lived due to the causality. Thus, there might be a critical value of the bias parameter b d,crit , above which the DWs are finite and completely decay in a finite time. On the other hand, for biases smaller than the critical value the DW length per horizon as well as the ratio r may asymptote to some constant values. This kind of arguments will help us to understand the stability of the DW network originated from scale-dependent superhorizon fluctuations. For instance, if the initial fluctuations are not scale-invariant, but blue-tilted, there may be an upper bound on the size of the DWs, which could be many orders of magnitudes larger than the Hubble horizon at the formation time. In this case, the largest DW size roughly determines the lifetime of the whole DW network.

A DW model of subhorizon fluctuation with bias
In the main text, we focused on the inflationary fluctuations, while the Gaussian white noise is considered only for comparison. On the other hand, the formation and evolution of DWs based on such biased white noise initial fluctuations have been investigated in the past literature [18][19][20][22][23][24]. Therefore, it would be useful to provide an explicit UV model that actually realizes such initial condition.
One way is to introduce another axion, a , with a mixing with an axion a that constitutes the DW in the low energy. If the global U(1) symmetry for a is spontaneously broken after inflation, axionic strings are produced. On the other hand, the U(1) symmetry for a is assumed to be already spontaneously broken during inflation. If a combinations of two U(1)'s for a and a is explicitly broken, those strings annihilate when the corresponding DWs are generated. Then, we are left with almost white noise initial condition for the remaining lighter degrees of freedom [29].
As another example, let us build a model with ap- proximate Z 2 symmetry, such that the adjacent vacua are physically identical and therefore strictly degenerate in energy. We introduce an axion, a, which enjoys the discrete shift symmetry, with an approximate (dark) charge conjugation symmetry, C dark . The interacting Lagrangian is given as, Here Φ is a Peccei-Quinn (PQ) field for the axion, and Ψ L,R are PQ fermions, while the fermion mass M is the explicit breaking of the PQ symmetry. Under the C dark symmetry, the fields transform asΨ L ↔ Ψ R , Φ ↔ Φ * , and therefore, both y and M must be real. Now let us consider the low energy effective theory for the axion, a, after the PQ spontaneous symmetry breaking. The Nambu-Goldstone boson, a, resides in the phase as Φ = fa √ 2 exp [ia/f a ], and it obtains the potential of The axion effective (thermal) potenital is given by and M 2 eff ≡ yfa √ 2 e ia/fa + M 2 ⊃ √ 2(yf a M ) cos(a/f a ). Here we assume that Ψ is thermalized at temperature T . We obtain a thermal mass contribution at T 2 M 2 eff as Since both yf a and M are real thanks to the C dark symmetry, the thermal mass at high enough T can bring the axion to a = 0 if yM > 0. We may also have another source for the axion potential from e.g. non-perturbative effects of non-Abelian gauge symmetry, in which case the total potential is V (a) = Λ 4 cos(a/f a ) + δV eff .
As long as the explicit breaking term satisfies C dark , the potential always has the symmetric enhanced point at a = 0 since a → −a under C dark . Thus a = 0 is either the potential maximum or minimum. Since the axion is at the hilltop due to the thermal mass at a high enough temperature, the initial condition for the domain wall formation is set by thermal effects. At the low temperature, the symmetry enhancement point may turn out to be a local maximum, and thus the domain wall is formed at H ∼ m a with m a being the axion mass in the vacuum. If we introduce a tiny breaking the C dark symmetry, e.g., by adding a phase, θ = 0 for the coupling, y → exp (iθ)y. Then the thermal potential, in general, does not necessarily bring the axion exactly at the origin, and the misalignment will induce the bias of the axion distribution. On the other hand, in spite of such bias, we still have degenerate vacua due to the periodicity of the potential, i.e., the potential does not change with a → a + 2πf a .