Gravitational Waves Produced by Domain Walls During Inflation

We study the properties of the stochastic gravitational wave background (SGWB) produced by domain walls (DWs) during inflation without forming a network. We numerically simulate the DW production caused by a second-order phase transition and calculate the SGWB spectrum using a $1000\times1000\times1000$ lattice. We show that the SGWB can be observed directly by future terrestrial and spatial gravitational wave detectors and through the B-mode spectrum in CMB. This signal can also explain the common noise process observed by pulsar timing array experiments. With numerical simulations, we derive an empirical formula for the strength and qualitative features of the SGWB spectrum. The details of the SGWB spectrum also contain information about the later evolution of the universe.


INTRODUCTION
The direct discovery of gravitational waves (GWs) produced by black hole binary merger opens up a new era of GW astrophysics [1].GWs can also be produced in the early universe via phase transitions, topological defects, scalar perturbations, and quantum fluctuation during inflation, forming the stochastic GW background.See Ref. [2] for an excellent review and references therein.GWs, once produced, propagate almost freely through the universe, bringing us information about their origin and the evolution of the universe.Terrestrial and spacebased detectors are proposed to search for SGWB [3][4][5][6][7][8][9][10][11][12][13][14][15].SGWB with a frequency around 10 −8 Hz can be detected via radio telescopes using the pulsar timing arrays (PTAs) [11,16,17].The common noise process observed by the PTAs might indicate a signal of SGWB [18][19][20][21].SGWB with longer wavelengths may also leave prints on the cosmic microwave background and be detected via the B-mode spectrum [22][23][24] One important source of SGWB that has been extensively studied in the literature is the cosmic DWs [25][26][27][28][29][30][31][32] (for a review, see [33]).DW-induced GWs as a test for new physics beyond the standard model have attracted a lot of interests in the literature [34][35][36][37][38][39][40][41][42].DWs are 2D topological defects created via the spontaneous breaking of discrete symmetries.Stable domain walls will form a network in the radiation domination (RD) era and overclose the universe [25].Thus an explicit symmetry broken term is usually introduced to ensure the DWs annihilate before the big bang nuclear synthesis [26][27][28].Traditional studies of SGWB produced by DWs focus on production through the DW network.Both qualitative analysis [29] and numerical simulations [30] show that the GW energy density production rate by the DW network is constant during RD.Therefore the GWs produced right before the annihilation of the DWs dominate the contribution to the SGWB energy density today since they are the least redshifted.However, if the DWs annihilate too early, the SGWB signal will be too weak to induce detectable signals in future GW detectors.
In this letter, we point out that if the DWs were produced during inflation, they would generate SGWB during inflation without forming a network, and the strength of the SGWB can reach the sensitivities of the planned GW detectors.In flat space-time or during the RD era, static sources cannot radiate GWs due to energy conservation.However, during inflation, the Hubble expansion rate is significant, and thus energy is no longer conserved.Therefore even comovingly static sources can produce GWs.If the DWs are produced about 60 e-folds before the end of inflation, the SGWB will induce a CMB B-mode spectrum that future CMB observers can detect.If the phase transition happens at about 40 e-folds before the end of inflation, the SGWB the DWs produced will be able to explain the common noise process observed by PTAs [18][19][20][21].

DW PRODUCTION DURING INFLATION
This work studies the phase transition in a spectator sector during de Sitter inflation.With conformal time τ , the scale factor is a(τ ) = −1/Hτ .We assume the spectator sector is in the Landau-Ginzburg type where σ is the spectator field and serves as the order parameter for the second-order phase transition.We consider two scenarios.In scenario (A), we assume the universe was in RD before inflation.At the beginning of inflation, the temperature of the thermal plasma quickly falls, and thus phase transition happened in the plasma.This scenario has been considered in the context in of grand unification models [43].In this case, we have m 2 eff = yT 2 −m 2 , where T is the temperature in the spectator sector, m is the mass of σ at zero temperature, and y is a constant parameter.In this scenario, phase transition happens at the critical temperature T c = y −1/2 m.In scenario (B), we consider the inflaton-driven phase transition similar to the study for first-order phase transitions in Refs.[44,45].We observe that in large field inflation models, the excursion of the inflaton field can be as large as the Planck scale.Thus, if we assume that the spectator fields are directly coupled to the inflaton field, the evolution of the inflaton field may induce dramatic changes of the masses or couplings in the spectator sector and thus leads to phase transitions during inflation.The effective mass can be written as m 2 eff = yφ 2 − m 2 , where φ is the inflaton field.
In this work, we use lattice simulation to numerically study the evolution of the spectator field σ after the phase transition.The details of the simulation method are presented in the appendix.Here we describe the qualitative features of evolution and the DW formation process.
The production of the topological defects can be described by the Kibble-Zurek mechanism [46,47].Around the critical time τ c , we can use the potential to describe the evolution of σ, where a c ≡ a τ =τc .For scenarios (A) and (B), the Kibble-Zurek scale m KZ can be written as [48] where is the slow-roll parameter, and φ 0 is the homogeneous part of φ.In the simplest slow-roll model, [49], with N c the e-folds between the critical point and the end of inflation.Therefore, compared to (A), m KZ in (B) is suppressed by a factor of N 1/3 c .In this work, for the phase transition to complete during inflation, we assume m is large enough, such that m KZ H. Right after the phase transition, the evolution of the Fourier modes of the σ field is driven by the linearized equation of motion where the primes denote the derivatives to τ , and the temporary frequency ω k (τ ) is given by where σ 2 (τ, x) is the temporary expectation value of σ 2 (τ, x).Around τ c , σ is still located at 0; thus σ 2 (τ ) can be neglected.Therefore ω 2 k is negative for the infrared modes with small enough k 2 , and thus σ k grows exponentially.After a while, when the nonlinear term in the potential becomes important, the exponential growth stops.Our numerical simulation shows that only modes with ka −1 (τ c ) < O(m KZ ) can experience exponential growth.Then we follow the strategy in Ref. [50] to simulate the evolution of the σ field after the phase transition.Due to the exponential growth, the particle numbers for long-wavelength modes are large.Thus we can use the classical Gaussian random field to approximate the quantum modes to set the initial condition for the nonlinear evolution.Then we use a 1000 × 1000 × 1000 lattice to simulate the nonlinear evolution classically.The details of the simulation are presented in the appendix.
Fig. 1 shows the evolution of the temporary rootmean-square of σ for scenarios (A) and (B) together with the temporary vacuum expectation values v temp ≡ λ −1/2 m eff .The curves for σ rms and v temp are normalized to v = λ −1/2 m.One can see that in (A), since the temperature T redshifts as a −1 , v temp approaches v after about one e-fold.However, in (B), since the evolution is driven by φ 0 , v temp will not reach v until the end of inflation.At the beginning of the phase transition, the growth of σ rms lags behind v temp , which causes the oscillation of σ rms around.As a result, σ rms oscillates around v temp in the later oscillation.However, the oscillation is significantly damped after several e-folds due to the Hubble friction, and then the σ field configuration becomes comovingly static.The configurations for σ at different times are shown in Figs. 2 and 3.In both cases, the field configuration becomes static after several e-folds.In Fig. 2, the x and y axes are the comoving coordinates in the unit of the comoving Kibble-Zurek length, m −1 KZ a −1 c .One can see that once the DWs are formed, their comoving density does not change, and the average distance between the walls is about 2πm −1 KZ a −1 c , as predicted, since only modes with ka −1 c < O(m KZ ) can grow exponentially.Figure 3. Evolution of the σ field in scenario (B), the labels of the axes are the same in Fig. 2. The parameters in the potential are the same as in Fig. 1.

GW PRODCTION
The Fourier modes of the tensor perturbation, hij (τ, k), satisfy the linearized Einstein equation where T T T ij is the transverse, traceless part of the energymomentum tensor.During inflation, when h(τ, k) leaves the horizon (k|τ | < 1), they will be frozen to a fixed value hf (k).After inflation, when the mode reenters the horizon, hf will serve as the amplitude for later evolution.Then the general form of h where Ẽi 0 and φ are determined by the evolution of the universe before the mode reenters the horizon.The detailed form of Ẽi 0 (k) and φ can be found in [45].Then we calculate the GW energy density spectrum, where G N is the Newton gravity constant, and V is the total comoving volume which will be canceled by the same volume factor in | hf ij (k)| 2 .Using the Green's function method, we can calculate hf [44,45], where K is the integral kernel, Eq. ( 9) indicates that during inflation because the time translation symmetry is badly broken, GWs can even be produced by a stationary source.The plot of K is shown in Fig. 4, from which we can see that the highest peak of K is around kτ ≈ −2, and this is the moment After doing a small kτ expansion in the integral (9), we can see that the integral is finite at τ = 0, and thus the integral is mostly contributed in the region kτ ∼ −2.(B), in which the green curves for ln(a/a c ) = 2.3, 3.3, 4.3 are identical in both scenarios.The coincidences indicate that without DWs, the GW production stops at about ln(a/a c ) = 2.3 because the oscillations of σ are quickly damped by the Hubble expansion, as shown in Fig. 1.This also agrees with the result in [31].From Eq. ( 5), one can also see that without DWs, the GW spectrum would be significantly smaller.Therefore, in the simulation, the DWs are the dominant source for the SGWB.

DETECTABILITY OF THE GW SIGNAL
Today's GW relative abundance can be calculated from Eq. (8), where Ω R and ρ R are today's radiation abundance and energy density.f = k/(2πa today ) gives today's GW fre-quency.The detailed shape of Ω GW (f ) also depends on the universe's evolution when the GW modes reenter the horizon.From the end of inflation to the completion of reheating, the universe may undergo some transition eras, such as matter domination and kination domination [51,52].As shown in the appendix (see also in [45]), the peak value of Ω GW (f ) is not sensitive to the universe's evolution when the modes are outside the horizon.During kination domination, the total energy density of the universe drops as a −6 , whereas the energy density of GW drops as a −4 , and therefore Ω GW (f ) gets a relative enhancement [53,54].On the other hand, if the GW modes reenter the horizon in an intermediate matter dominated era before reheating, the total energy of the universe drops as a −3 .Thus Ω GW , in this case, obtains a relative suppression.From qualitative arguments and numerical simulations in the appendix, we can derive a semi-analytical formula for the peak value of Ω GW , where d w is the physical wall thickness at τ = −2/m KZ .z mp in ( 12) measures the redshift of the universe between the time the peak mode reenters the horizon and the end of the intermediate stage.Future GW detectors can detect SGWB produced during inflation by the DWs.Fig. 6 shows today's SGWB  spectrums and the sensitivities of PTAs and space and ground-based GW telescopes.Here we consider scenario (B), and the parameters for the spectator sector and inflaton model are the same as for the lower panel of Fig. 5.We consider three after-inflation-evolution scenarios, the instantaneous reheating, intermediate MD, and intermediate KD.As shown in Fig. 6, for N c ≈ 40, the SGWB spectrum falls in the region of PTAs, and it is possible to use it to explain the common noise process.For N c ≈ 20 − 25, the SGWB can be detected by the planned space GW detectors, such as LISA, Taiji, and Tianqin.For N c ≈ 10, it can be detected by future terrestrial GW detectors.
If N c ≈ 60, the GW frequency can be about 10 −19 Hz and can be detected via the CMB B-mode spectrum.Here we consider scenario (A), in which the back reaction to the curvature perturbation will not ruin the temperature correlation in the CMB spectrum [43].With reasonable choices of parameters, the future CMB-S4 experiment can detect the signal.Furthermore, we can also see that the shape of the B-mode spectrum induced by the phase transition differs from the quantum-induced one.Therefore, future CMB observations will be able to distinguish the origin of the SGWB source.

SUMMARY AND DISCUSSIONS
In this work, we numerically study the SGWB produced by DWs during inflation.We show that due to the severe violation of the time translation symmetry during inflation, even the cosmologically static DW configurations, without forming a network, can create detectable SGWB for future GW detectors.We also show the SGWB produced in scenarios can be used to explain the common red noise observed by NanoGrav.
To emphasize the DW-induced SGWB during inflation, we assume that the DWs decay sufficiently fast such that the DW network disappears before dominating the universe, and the GWs produced by the DW network are negligible.
Scenario (A) discussed in this work may also be realized in warm inflation models, where the decay of the inflaton field supports the thermal plasma [55,56].In scenario (B), the interaction between the φ and σ leads to a back reaction to curvature perturbation, and the Lagrangian can be written as yσ 2 φ 0 δφ.This term may induce secondary GWs once the induced curvature perturbations reenter the horizon.This term may also lead to the production of primordial back holes.We leave these interesting phenomena to future studies.

Appendix
In the Appendix, we present the details of our numerical simulation of the formation of the DWs and calculation of the SGWB.Additionally, we derive a semi-analytical formula for the peak value of the SGWB spectrum.Thirdly, we present the detailed calculation of the influence on the SGWB spectrum induced by the evolution of the universe after inflation.

DETAILS OF THE NUMERICAL SIMULATION
Here we present the detailed simulation of the evolution of the spectator field σ during the second-order phase transition.
We use comoving coordinate with the conformal time, and the background de Sitter metric can be written as where a(τ ) = −1/Hτ is the scale factor.The Lagrangian of the spectator σ can be written as where, the primes indicate derivatives with respect to τ .In this work, we assume the potential is in the Landau-Ginzburg form, We discuss two scenarios.In scenario (A), the phase transition is triggered by evolution of the temperature at the beginning of inflation, and m 2 eff = yT 2 − m 2 .In scenario (B), the phase transition is triggered by the evolution of the inflaton field, and m 2 eff = yφ 2 − m 2 .Around the critical time τ c , the potential can be written as Tachyonic growth of the long wavelength modes Now we consider the quantum behavior of σ at the beginning of the phase transition.Around τ c , the σ field still locates at σ = 0, and thus the interaction term can be neglected.As a result, m 2 eff < 0, and the long wavelength σ modes grow tachyonically.Here we generalize the discussions in [50] to dS space.
The Hamiltonian governs the evolution of σ in the linear region is where is the canonical momentum.Then the Fourier modes of σ and φ satisfy the Hamiltonian equations of motion To quantize the theory, we impose the equal-time canonical commutation relation to σ We expand σ by ladder operators a k , a † −k , The solution can be expressed as where the mode function f depends only on k since the space is homogeneous and isotropic and satisfies the equation of motion We require the ladder operators to satisfy the commutation relation This leads to the normalization condition for the mode function The time scale of the phase transition is determined by m KZ , which is much larger than H, and thus when setting the initial condition of f at τ c , we can neglect the Hubble expansion.Therefore just as in [50], we choose (S14) as the initial condition for (S11).
The modes become classical when the anticommutation of σk and πk is significantly larger than their commutation, where the expectation value is over the temporary vacuum state at τ c .This is equivalent to require where the function F is defined as A plot of F (k) is shown in Fig. S2 for different values of k and H.In this plot, to compare with Minkowski space (H = 0) result in [50], we use the physical time t in Fig. S2, defined as dt = adτ .One can see that F grows exponentially after the phase transition for the modes with k < m KZ as long as H is significantly smaller than m KZ .Then the canonical phase space distribution of the long-wavelength modes can be described by the Gaussian ensemble given by the Wigner function [57][58][59]  After τ c , due to the tachyonic growth, the nonlinear term λσ 4 becomes more and more important, and then when we consider the evolution of σ(k), we must include the back reaction from the long-wavelength modes.Thus the linearized equation of motion for σ(k) can be written as with where the expectation value of σ 2 (τ, x) includes all the contributions from the long-wavelength modes that have grown exponentially.Fig. S3 shows the evolution of ω 2 for different values of k with H = 0.1m KZ and λ = 0.1.Note that σ(k) can grow tachyonically only when ω 2 < 0. Thus one can see that only modes with k < O(m KZ ) can grow tachyonically.This determines that the size of the topological structure is about 2π/m KZ as shown Figs. 2  and 3.
We should note that the time span for ω 2 < 0 is not sensitive to the value of λ, since after the phase transition, σ 2 (x) grows as e [mKZ(t−tc)] 3 .Follow [50], we choose to use m KZ (t − t c ) = 2 as the matching point between tachyonic growth and the nonlinear classical evolution.In detail, we solve numerically Eq. (S11) with the initial condition (S14).Then at m KZ (t − t c ) = 2, we randomly generate σ(k) and π(k) according to the distribution function (S18).Then from σ and π we calculate the initial configuration of σ(x) and π(x) for the classical lattice simulation of the nonlinear evolution.

Nonlinear classical evolution
The evolution of classical systems can be described by the Hamiltonian equations, It's convenient to use the dimensionless variables (S25) Then we have For the case studied in this work, the Hamiltonian can be divided into two parts, with (S30) We use the symplectic integrator method [60,61] to numerically calculate the evolution of Σ. Define Here the infinitesimal evolutions A and B are generated by the dimensionless kinetic energy potential energies K /m and V /m.An nth order symplectic integrator is defined as a symplectic approximation of the evolution operator with a difference of order O(∆T n+1 ) e ciA∆T e diB∆T +O(∆T n+1 ), (S33) where k = 8 for sixth order accuracy used in this work.In our simulation we use the Solution A of 6th order symplectic integrator proposed in [61], and the values for the coefficients are On lattice, derivative is replaced by difference.We use the neutral differential defined as The Laplacian is defined as Both ∆ x and ∆ 2 x are at O(δx 4 ) accuracy.∆ x f (n) can be written as a general form and in our case we have With D(n) we can define the effective momentum Then we define the transverse projective operator as [62,63] to suppress the unphysical modes in lattice calculation of the GWs.Specifically, in our choice of D(n) in (S42) we have Here we present a derivation to Eq. ( 12) in the main text and discuss the IR and UV behavior of the SGWB spectrum, and compare the analytical result with numerical simulations.
Due to the shape of the kernel K, as shown in Fig. 4, the dominant contribution of T T T ij (τ , k) to hf ij (k) happens at kτ ≈ −2.To estimate the peak value of Ω GW we only keep this dominant contribution.Thus we have Here • • • means taking average over the direction of k.
The space average value of T T T (τ, x) 2 can be written as As discussed in the main text, the largest correlation of σ field after the phase transition is determined by m KZ , and the typical size of independent regions in the universe is about 2πm −1 KZ .Thus, to calculate the spatial integral in (S47), we can first define a window function W (a) for each region.The center position of each region locates at x (a) .Then the energy-momentum tensor can be decomposed as Then we can define ) .

(S49)
Therefore, we have Then the righthand side of (S47) can be written as The interference contributed from different patches is suppressed.Therefore it can be further simplified as In comoving coordinates the size of each independent region is about 2πm −1 KZ a −1 c .Therefore the above expression can be simplified as where A 3 is a numerical factor and | T T T (A) ij (τ, k)| 2 can be seen as the space and angular averaged value of We know that the modes with k ∼ m KZ a c dominates the contribution to the SGWB and as we will see later the p integral in (S54) is dominated by the region p ∼ m m KZ .Therefore, as a qualitative estimate, we can neglect the k in σ.Then the angular average of T T T ij (τ, k) T T T ij (τ, k) contributed by a single region can be written as A typical DW configuration is shown in Fig. S4, where we only present the x and z coordinates of the space.We can decompose the comoving momenta p and q in Eq. (S56) into the direction that is in perpendicular to the surface of the wall (the z direction as in Fig. S4), and the directions that in parallel to the wall surface.For a generic DW, if we neglect its thickness it can be seen as a Heaviside theta function.In our case, if we only consider the field configuration along the z direction, it can be written as The Fourier transformation of the above configuration is Therefore, we expect that, for a general configuration, σ(p) has the form in the region that m KZ a c < p < d w a(τ ), where d w is the physical wall thickness.As an example, in the Landau-Ginzburg type model, as discussed in this work, the DW configuration can be approximated as The Fourier transformation is Then we know that for p z < m, σ(p z ) drops as p −1 z , and for p z > m, σ(p z ) drops exponentially.
In the parallel directions, since the oscillations on the surface of the wall are fast diluted by the expansion of the universe, the configuration changes slowly.Therefore we don't expect σ to have significant support for p > m KZ a −1 c .Now, we can decompose the p and q integrals in Eq. (S56) into d 2 p dp ⊥ d 2 q dq ⊥ , and then it becomes where the p ⊥ and q ⊥ integrals are from m KZ to d −1 w a(τ ).From (S61) one can see that the p ⊥ and q ⊥ integrals are dominated in the region p ⊥ ∼ d −1 w a(τ ) and q ⊥ ∼ d −1 w a(τ ), whereas the supports of the p and q integrals are around m KZ a c .Thus, when estimate the value of (S61), we can use the approximation p , q p ⊥ , q ⊥ .Another thing to notice is that σ must be proportional to v temp , the temporary vev of the σ field.Then (S61) can be estimated as where [m KZ a −1 c ] 4 in the denominator is from dimensional analysis and A 4 collects all the numerical factors when calculating the integral.Together with (S46), (S53), (S62), Eq. ( 8) in the main text, and a(τ = −2/k) = k/(2H) , we arrive at where A 1 collects all the numerical factors.Now, let's focus on the case that reheating happens fast and finish within one e-fold.Then, as shown in [45], where a end is the scale factor at the end of inflation.Thus we have, in the RD era, In the instantaneous reheating case, neglecting the change of the number of relativistic degrees of freedom, the radiation energy density can be written as Thus, we have In Landau-Ginzburg type model, as shown in Eq. ( 1) we have Therefore we have During inflation, we have Then we have This result assumes k ma c and thus we can neglect the k dependence in σ when calculating | T T T ij | 2 .However, we see in Fig. S5, | T T T ij | 2 drops significantly when k > m KZ a c .Therefore the peak position of Ω GW is around k = m KZ a c , and we arrive at the formula where A collects all the numerical factors.For scenario (A), as shown in Fig. 1 in the main text, v temp arrives at mλ −1/2 within about one e-fold after the critical time.Thus, to estimate Ω GW , we can just use the parameters in the Lagrangian.However, for scenario (B), v temp changes slowly with φ 0 after the phase transition.Therefore, to have an empirical formula for Ω peak GW for scenario (B), we take τ = −2/(m KZ a c ) to evaluate d w and ∆ρ.
A by-product of (S71) is that it predicts that on the left side of the peak, Ω GW increases as k 1 until it reaches the peak.For k m KZ a c , from causality argument, we have Ω GW ∼ k 3 [64].Thus Eq. (S71) predicts that there is a transition of Ω GW from k 3 to k 1 before reaching the peak frequency.This behavior is clearly shown by the numerical simulation presented in Fig. S6.The k 1 behavior before reaching the peak can be seen as a nontrivial check to the numerical simulation.Here we assume that the modes reenter the horizon during RD era.
In the far IR region, ΩGW increases as k 3 , and then transits into k 1 when approaching the peak.Then it decays as k −4 in the UV region.
The analytical formula (S72) also shows that Ω peak GW is proportional to H 2 .This is also confirmed by the numerical simulations.In Fig. S7, we show the Ω peak GW as a function H/m.We can see clearly that the qualitative behavior shown in (S72) is correct.

Figure 2 .
Figure 2. Evolution of the σ field in scenario (A) for ln a/ac = 1.3, 2.3, 3.3 and 4.3, where the x and y axes show the comoving coordinates in the unit of a −1 c m −1 KZ , and the z axis shows σ/v.The parameters in the potential are the same as in Fig. 1.

Figure 4 .
Figure 4.The kernel K as a function of η ≡ kτ .

Fig. 5
shows the spectrum of | hf ij | 2 /V accumulated from τ c to ln(a/a c ) = 1.3, 2.3, 3.3, 4.3 for both scenarios (A) and (B).One can see that the GW productions for both scenarios stop at about four e-folds after the phase transition.Th As comparisons, we also did simulations without DW formations.We use the same parameters in the simulations without DW formation to generate the initial configuration σ(x) induced by the exponential growth.Then we replace σ(x) with |σ(x)| for the initial condition for the nonlinear growth.This way, we eliminate all the DWs while keeping the initial fluctuations of the σ field.The green curves in Fig.5show the accumulated | hij | 2 /V without DWs for both scenarios (A) and

Figure 5 .
Figure 5.The spectrum of | hf ij | 2 for scenarios (A) (up) and (B) (down).The choices of parameters in the potential are the same as in Fig. 1.The spectrum is shown accumulated up to ln(a/ac) = 1.3, 2.3, 3.3 and 4.3.As a comparison, the spectrum without DW production is also shown by the green curves, where for both (A) and (B), the ln(a/ac) = 2.3, 3.3, 4.3 curves are identical, which indicates that the GW production ceases completely after ln(a/ac) = 2.3.

Figure 6 .
Figure 6.Spectrums of today's SGWB energy density distribution for scenario (B), together with the sensitivity curves of future GW detectors and the region favored by the NanoGrav observation.

2 ]Figure 7 .
Figure 7. B-mode power spectrums induced by the SGWB in scenario (A) for different choices of phase transition time and potential energy in the spectator sector.As a comparison, the B-mode spectrum generated via tensor mode quantum fluctuations with tensor-scalar ratio r = 0.003 is also shown.The orange dots and triangles are for the data from BICEP/KECK array.

Figure S3 .
Figure S3.Evolution of the mode function f (k) as a functions of the physical time t for different values of k and H.

Figure S4 .
Figure S4.The illustration of DW configuration.
Fig. S5 shows the | T T T ij | 2 at different time after the phase transition.Once can see that the distributions are flat for small k, and drops as k −3 ∼ k −4 when k > m KZ a c .

5 V - 1 Figure
Figure S5.| T T T ij | 2 at different time as a function of k.

1 .
Figure S6.Qualitative behavior of the GW spectrum.Here we assume that the modes reenter the horizon during RD era.In the far IR region, ΩGW increases as k 3 , and then transits into k 1 when approaching the peak.Then it decays as k −4 in the UV region.

2 ν[ 2 x 2 xA 2 |S| 2 +ΩFigure S8 .
Figure S8.The distortion function D(k) as functions of x for MD intermediate stage (up), and KD intermediate stage (down), respectively. ) kz -1/2 Figure S1.Evolution of the mode function f (k) as a functions of the physical time t for different values of k and H. KZ t)| Figure S2.F (k) as functions of physical time t for different values of H.
After inflation, the universe might undergo an intermediate stage before going into the RD era.The GW modes may reenter the horizon in the intermediate stage, and