Phase Transitions in an Expanding Universe: Stochastic Gravitational Waves in Standard and Non-Standard Histories
PP REPARED FOR SUBMISSION TO
JCAP
Phase Transitions in an ExpandingUniverse: Stochastic GravitationalWaves in Standard andNon-Standard Histories
Huai-Ke Guo a Kuver Sinha a Daniel Vagie a Graham White b a Department of Physics and Astronomy, University of Oklahoma, Norman, OK 73019, USA b TRIUMF Theory Group, 4004 Wesbrook Mall, Vancouver, B.C. V6T2A3, CanadaE-mail: [email protected], [email protected], [email protected], [email protected]
Abstract.
We undertake a careful analysis of stochastic gravitational wave production from cos-mological phase transitions in an expanding universe, studying both a standard radiation as well asa matter dominated history. We analyze in detail the dynamics of the phase transition, including thefalse vacuum fraction, bubble lifetime distribution, bubble number density, mean bubble separation,etc., for an expanding universe. We also study the full set of differential equations governing theevolution of plasma and the scalar field during the phase transition and generalize results obtained inMinkowski spacetime. In particular, we generalize the sound shell model to the expanding universeand determine the velocity field power spectrum. This ultimately provides an accurate calculationof the gravitational wave spectrum seen today for the dominant source of sound waves. For the am-plitude of the gravitational wave spectrum visible today, we find a suppression factor arising fromthe finite lifetime of the sound waves and compare with the commonly used result in the literature,which corresponds to the asymptotic value of our suppression factor. We point out that the asymptoticvalue is only applicable for a very long lifetime of the sound waves, which is highly unlikely due tothe onset of shocks, turbulence and other damping processes. We also point out that features of thegravitational wave spectral form may hold out the tantalizing possibility of distinguishing betweendifferent expansion histories using phase transitions. a r X i v : . [ h e p - ph ] J u l ontents β and Mean Bubble Separation ( R ∗ ) 17 Primordial stochastic gravitational waves from first order cosmological phase transitions have becomea new cosmic frontier to probe particle physics beyond the standard model [1–6]. Alongside exten-sive studies on the theory side, direct searches for stochastic gravitational waves at LIGO and Virgohave also been performed using their O1 and O2 data sets [7, 8]. Perhaps even more significantly,many space-based detectors have been proposed, such as the Laser Interferometer Space Antenna(LISA) [9], Big Bang Observer (BBO), DECi-hertz Interferometer Gravitational wave Observatory– 1 –DECIGO) [10], Taiji [11], and Tianqin [12]. They will come online within the next decade or so andcan probe lower frequencies coming from an electroweak scale phase transition [13–28]. Precise calculations of the gravitational wave power spectrum are required to have any hope ofinferring parameters of the underlying particle physics model. There have been significant advancesin this direction in recent years. In particular, it is now generally accepted that the dominant source forgravitational wave production in a thermal plasma is the sound waves [48], although a more preciseunderstanding of the onset of the turbulence is still needed to settle this issue. For the acousticproduction of gravitational waves, many large scale numerical simulations have been performed [49,50], with the result that standard spectral formulae are now available for general use. These resultshave also been understood reasonably well for relatively weak transitions, through the theoreticalmodeling of the hydrodynamics [51] and with the recently proposed sound shell model [52, 53].The first major goal of this paper is to undertake a careful analysis of the gravitational wavepower spectrum in a generic expanding universe. This is necessary, since the standard result forthe spectrum is obtained in Minkowski spacetime where the effect of the expansion of the universeis neglected. The spectrum depends on the lifetime of the source (the sound waves); it has beenshown in [49], based on simple rescaling properties of the fluid, that the effective lifetime of thesource is a Hubble time, such that H ∗ τ sw ≈ , where H ∗ is the Hubble parameter at the phasetransition and τ sw the lifetime of the source. We note, however, that many approximations were usedin arriving at this result, some of which may bear deeper scrutiny. We present a comprehensive andvery careful analysis of the spectrum, clarifying subtle issues when the calculation is generalized fromMinkowski spacetime to an expanding universe, and ultimately providing an accurate spectrum in astandard radiation dominated universe and in other expansion scenarios. We also perform a detailedcalculation of the nucleation and growth of bubbles in an expanding background, including trackingthe shrinking volume available for new bubbles to nucleate in as well as the total area of uncollidedwalls. Both are needed for an accurate understanding of how the volume fraction and mean bubbleseparation evolve throughout the phase transition. We then derive and solve the equations governingthe evolution of the fluid velocity field in an expanding Universe and then proceed to a derivation ofthe spectrum for different expansion scenarios.The second major goal of this paper is encapsulated in the title: after having calculated thegravitational wave spectrum in an expanding universe, we want to explore the extent to which thephase transition can distinguish between different expansion histories . In other words, we would liketo interrogate how well a phase transition can serve as a cosmic witness . This is important, sincegrowing evidence suggests that the standard assumption of radiation domination prior to Big BangNucleosynthesis may be too naive [54, 55]. An early matter dominated era, for example, is motivatedby the cosmological moduli problem [56–59], hints from dark matter searches [60–68], and perhapseven baryogenesis [69]. Another possibility of a non-standard expansion history is kination, whichwe do not cover in this paper but can be explored by our methods [70–79]. We note that gravitationalwaves have been previously employed to investigate early universe cosmology [80–85].Our goal is to provide a general theoretical framework to calculate the gravitational wave spec-trum in different cosmic expansion histories. This includes scrutiny for changes in different aspects.The dynamics of the phase transition in an expanding universe is studied in Sec. 3, the velocity fieldpower spectrum is calculated in Sec. 4 and the gravitational wave spectrum in Sec. 5. The mainfindings of the first two aspects are as follows.1. The mean bubble separation R ∗ is related to β through a generalized relation for the exponential Note that they are also poised to probe hidden sector transitions [29–43] and transitions from multi-step GUT breaking[44–47] – 2 –ucleation (Eq. 3.50): R ∗ ( t ) = a ( t ) a ( t f ) (8 π ) / v w β ( v w ) , (1.1)where t f is the time when the false vacuum fraction is /e , at which β ( v w ) is evaluated, and β ( v w ) can vary by ∼ for different v w . This relation is also confirmed by numericalcalculations and is accurate up to an uncertainty of . If one uses the conformal version of R ∗ and β , then they satisfy the same relation as in Minkowski spacetime (see Eq. 3.46).2. We derived the bubble lifetime distribution in a generic expanding universe in Eq. 3.29, andthe conformal lifetime η lt rather than ordinary lifetime t lt should be used. It coincides with thedistribution e − ˜ T found in Minkowski spacetime [53] for exponential nucleation.3. We derived the full set of differential equations in an expanding universe for the fluid and orderparameter field model as used in numerical simulations. We find that in the bubble expansionphase the full field equations do not admit rescalings of the quantities that would reduce theexpressions to their counterparts in Minkowski spacetime; this rescaling does, however, doeswork in the bag equation of state model. This implies the velocity profile maintains the sameform when appropriate rescalings and variable substitutions are used.4. We generalized the sound shell model to an expanding universe and calculated the velocityfield power spectrum [52, 53].For the gravitational wave energy density spectrum, the main results are:1. The peak amplitude of the gravitational wave spectrum visible today has the form (see Eq. 5.45) h Ω GW = 8 . × − (cid:18) g s ( T e ) (cid:19) / Γ ¯ U f (cid:20) H s β ( v w ) (cid:21) v w × Υ . (1.2)Here Γ ∼ / is the adiabatic index, g s ( T e ) is the relativistic degrees of freedom for entropy at T e when the gravitational wave production ends, ¯ U f is the root mean square fluid velocity (seeFig. 18), v w is the wall velocity, H s is the Hubble rate when the source becomes active, and Υ is the suppression factor arising from the finite lifetime, τ sw , of the sound waves. For radiationdomination, it is given by Υ = 1 − √ τ sw H s , (1.3)where the standard spectrum generally used corresponds to the asymptotic value Υ = 1 when τ sw H s → ∞ . However the onset of non-linear shocks and turbulence which can disrupt thesound wave source occurs at around τ sw H s ∼ H s R ∗ / ¯ U f . This means the asymptotic valuewill not be reached and there is a suppression to the standard spectrum. In Fig. 1 we compareour result with the suppression factor recently proposed in [23] (see also [86, 87]). Similarly,the spectrum for matter domination has also been derived in our work and a similar suppressionfactor Υ is observed, which has an asymptotic value of / .2. We find a change to the spectral form, depending upon whether the phase transition occurs dur-ing a period of matter or radiation domination. The change in the form is not leading order, dueto the fact that the velocity profiles remain largely unchanged and that the autocorrelation timeof the source is much smaller than the duration of the transition. This is in contrast to gravita-tional waves generated from cosmic strings [85]. Even then, the modification of the spectrum– 3 – .0 0.2 0.4 0.6 0.8 1.0 1.2 1.40.00.20.40.60.81.0 Figure 1 . The suppression factor (blue solid line) as a function of the lifetime of the dominant source, thesound waves, in unit of the Hubble time at t s , the time when the source becomes active. The black dashed linedenotes Min[ τ sw H s , . presents an enticing possibility that the gravitational waves formed during a phase transitioncan bear witness to an early matter dominated era. We leave a further detailed exploration ofthe change of the spectral form for future work.The remainder of this paper is organized as follows. We firstly lay out the theoretical frameworkfor the stochastic gravitational wave calculation in the next Sec. 2 and study the details of the phasetransition dynamics in an expanding universe in Sec. 3. After that, we summarize the full set offluid equations applicable in an expanding universe and study the velocity profile and velocity powerspectrum using the sound shell model in Sec. 4. We then analytically calculate the gravitationalwaves from sound waves in both radiation dominated and matter dominated scenarios in Sec. 5. Wesummarize our results in Sec. 6. In this section, we set up the framework for calculating the stochastic gravitational waves in the pres-ence of a source, which also serves to define our notation. The power spectrum of the gravitationalwaves, as will be discussed, depends on the unequal time correlator of the source. Therefore thiscorrelator is of central importance in this work and is discussed in the second subsection.
The gravitational wave is the transverse traceless part of the perturbed metric. Neglecting the non-relevant scalar and vector perturbations, the metric is defined in the FLRW universe as: ds = − dt + a ( t ) ( δ ij + h ij ( x )) d x , (2.1)– 4 –here h ij is only the transverse traceless part of the perturbed × metric matrix (see, e.g., [88]for a detailed discussion). It is convenient, most often, to work in Fourier space, with the followingconvention: h ij ( t, x ) = (cid:90) d q (2 π ) e i q · x h ij ( t, q ) , (2.2)where q is the comoving wavenumber, in accordance with the comoving coordinate x . The physicalcoordinate is a x and the physical wavenumber is q /a . The Fourier component h ij ( t, q ) is thus ofdimension − .Gravitational waves are sourced by the similarly defined transverse traceless part of the per-turbed energy momentum tensor of the matter content, defined by [88] T ij = a π Tij + · · · , (2.3)where “ · · · ” denotes the neglected non-relevant parts. Its Fourier transform is defined by π Tij ( t, x ) = (cid:90) d q (2 π ) qe i q · x π Tij ( t, q ) . (2.4)Since π Tij is of dimension 4, the dimension of its Fourier component π Tij ( t, q ) is 1. The Einsteinequation leads to a master equation governing the time evolution of each Fourier component of thegravitational waves, which is decoupled from the scalar and vector perturbations, h (cid:48)(cid:48) ij ( t, q ) + 2 a (cid:48) a h (cid:48) ij ( t, q ) + q h ij ( t, q ) = 16 πGa π Tij ( t, q ) . (2.5)Here (cid:48) ≡ ∂/∂η , with η being the conformal time. Derivatives with respect to the coordinate time willbe denoted by a dot. The gravitational wave energy density, as denoted by ρ GW here, is defined as ρ GW ( t ) = 132 πG (cid:104) ˙ h ij ( t, x ) ˙ h ij ( t, x ) (cid:105) , (2.6)with the angle brackets, (cid:104)· · · (cid:105) , denoting both the spatial and ensemble average. Due to the over-all spatial homogeneity of the universe, we can define the power spectrum of the derivative of thegravitational wave amplitude as: (cid:104) ˙ h ij ( t, q ) ˙ h ij ( t, q ) (cid:105) = (2 π ) δ ( q + q ) P ˙ h ( q , t ) . (2.7)Then the gravitational wave energy density follows ρ GW ( t ) = 132 πG π (cid:90) dq q P ˙ h ( t, q ) , (2.8)and the gravitational wave energy density spectrum: dρ GW ( t ) d ln q = 164 π G q P ˙ h ( t, q ) . (2.9)It is conventional to use the dimensionless energy density fraction of the gravitational waves Ω GW ( t ) = ρ GW ( t ) /ρ c ( t ) where ρ c is the critical energy density at time t . The corresponding dimensionless ver-sion of the spectrum is P GW ( t, q ) ≡ d Ω GW ( t ) d ln q = 124 π H q P ˙ h ( t, q ) = 124 π H a q P h (cid:48) ( t, q ) , (2.10) P GW is also denoted as Ω GW ( t, q ) . – 5 –here in the last step P h (cid:48) ( t, q ) is defined by replacing ˙ h with h (cid:48) in Eq. 2.7.We thus need to solve for h ij ( η, q ) by solving Eq. 2.5 together with equations governing theevolution of the source. We will follow the conventional approach by neglecting the back-reactionof the metric on the source and calculate the stress tensor with a modelling of the phase transitionprocess. Once π Tij ( t, q ) is provided in this way, then h ij ( t, q ) can be solved from Eq. 2.5 with Green’sfunction and with the following boundary conditions G (˜ η (cid:54) ˜ η ) = 0 , ∂G (˜ η, ˜ η ) ∂ ˜ η | ˜ η =˜ η +0 = 1 , (2.11)where ˜ η = qη , which is a dimensionless quantity and ˜ η is the time when the phase transition starts.With the Green’s function, the solution of the inhomogeneous Eq. 2.5 is given by h ij ( t, q ) = 16 πG (cid:90) ˜ η ˜ η d ˜ η (cid:48) G (˜ η, ˜ η (cid:48) ) a ( η (cid:48) ) π Tij ( η (cid:48) , q ) q , (2.12)and its derivative with respect to the conformal time follows simply: h (cid:48) ij ( η, q ) = 16 πG (cid:90) ˜ η ˜ η d ˜ η (cid:48) ∂G (˜ η, ˜ η (cid:48) ) ∂ ˜ η a ( η (cid:48) ) π Tij ( η (cid:48) , q ) q . (2.13)Then we can calculate the 2-point correlation function: (cid:104) h (cid:48) ij ( η, q ) h (cid:48) ij ( η, q ) (cid:105) = (16 πG ) (cid:90) ˜ η ˜ η d ˜ η (cid:90) ˜ η ˜ η d ˜ η ∂G (˜ η, ˜ η ) ∂ ˜ η ∂G (˜ η, ˜ η ) ∂ ˜ η × a ( η ) a ( η ) q (cid:104) π Tij ( η , q ) π Tij ( η , q ) (cid:105) . (2.14)Supposing that the gravitational wave generation finishes at ˜ η f , the upper limits for the integrals in theexpression above will be ˜ η f . Subsequently, the energy density of the gravitational waves for modesinside the horizon will be simply diluted as /a . We thus see that at the core of the gravitationalwave energy density spectrum calculation is the unequal time correlator (UETC) of π Tij . It can beparametrized in the following way due to the overall spatial homogeneity of the universe [49] (cid:104) π Tij ( η , q ) π Tij ( η , q ) (cid:105) = Π ( q , η , η )(2 π ) δ ( q + q ) . (2.15)It is obvious that the dimension of Π ( k, η , η ) is 5. Let us first write down the energy momentum tensor of the matter content in the universe. Herewe keep the dominant contribution from the fluid and assume the fluid velocities are non-relativisticfollowing Ref. [53], then T ij = a (cid:2) pδ ij + ( p + e ) γ v i v j (cid:3) ,T i = a (cid:2) − ( p + e ) γ v i (cid:3) ,T = γ ( e + pv ) , (2.16)where e is the energy density, p is the pressure and the velocity is defined w.r.t the conformal time v i = dx i /dη . Then, comparing with Eq. 2.3 and neglecting the non-relevant parts, we have π ij = ( p + e ) γ v i v j , (2.17)– 6 –ere the scale factor dependent ( p + e ) , takes its homogeneous value (defined with a bar) to leadingorder ¯ e + ¯ p ≡ ¯ ω and scales as /a . γ is the Lorentz factor. The calculation of the correlator of π Tij parallels that in Minkowski spacetime: (cid:104) π Tij ( η , k ) π Tij ( η , q ) (cid:105) = Λ ij,kl (ˆ k ) 1(2 π ) (cid:90) d x (cid:90) d y e − i k · x e − i q · y (cid:104) π Tkl ( η , x ) π Tij ( η , y ) (cid:105) , = Λ ij,kl (ˆ k )¯ ω π ) (cid:90) d x (cid:90) d y e − i k · x e − i q · y (cid:104) v k ( η , x ) v l ( η , x ) v i ( η , y ) v j ( η , y ) (cid:105) , = ¯ ω Λ ij,kl (ˆ k ) 1(2 π ) (cid:90) d q (cid:90) d q (cid:104) ˜ v k q ( η )˜ v l ∗ q − k ( η )˜ v i q ( η )˜ v j ∗ q − q ( η ) (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) ≡ X klij . (2.18)Here Λ ij,kl is the standard projection operator and Λ ij,kl (ˆ k ) = P ik (ˆ k ) P jl (ˆ k ) − P ij (ˆ k ) P kl (ˆ k ) with P ij (ˆ k ) = δ ij − ˆ k k ˆ k j . ˜ v i q is the Fourier transform of the velocity field v i ( x ) . Due to the natureof the first order phase transition process and according to the central limit theorem, ˜ v i q ( η ) followsthe Gaussian distribution to a good approximation. Also as in Ref. [53], we neglect the rotationalcomponent, then the two point correlator can be defined in the following way: (cid:104) ˜ v i q ( η )˜ v j ∗ k ( η ) (cid:105) = δ ( q − k )ˆ q i ˆ k j G ( q, η , η ) . (2.19)Any higher order correlator can be reduced to the two point correlator. Defining ˜ q ≡ q − k and ˜ q ≡ q − q , then X klij = (cid:104) ˜ v k q ( η )˜ v l ∗ ˜ q ( η ) (cid:105)(cid:104) ˜ v i q ( η )˜ v j ∗ ˜ q ( η ) (cid:105) + (cid:104) ˜ v k q ( η )˜ v i q ( η ) (cid:105)(cid:104) ˜ v l ∗ ˜ q ( η )˜ v j ∗ ˜ q ( η ) (cid:105) + (cid:104) ˜ v k q ( η )˜ v j ∗ ˜ q ( η ) (cid:105)(cid:104) ˜ v l ∗ ˜ q ( η )˜ v i q ( η ) (cid:105) . (2.20)The first term contributes trivially to k = 0 and, collecting all other contributions, we have (cid:104) π Tij ( η , k ) π Tij ( η , q ) (cid:105) = δ ( k + q )¯ ω π ) (cid:90) d q G ( q , η , η ) G (˜ q , η , η )(1 − µ ) q ˜ q . (2.21)Comparing with Eq. 2.15, it follows that Π ( k, η − η ) = ¯ ω (cid:90) d q (2 π ) G ( q, η , η ) G (˜ q, η , η ) q ˜ q (1 − µ ) , (2.22)where ˜ q = | q − k | and µ = ˆ q · ˆ k . Here Π depends on η − η rather than on η and η separately.This is because the source is largely stationary.We will later see that the fluid equations maintain the same form as in the Minkowski spacetimeonce properly rescaled quantities and previously defined v i ( x ) are used (see also Ref. [49]). Inparticular it means that we can define a rescaled stress energy tensor ( ˜ π Tij ) for the fluid: π Tij ( q , η ) = a s a ( η ) ˜ π Tij ( q , η ) , (2.23)– 7 –here a s is a reference scale factor when the source becomes active. Similarly we can define arescaled and dimensionless two point correlator ˜Π following Ref. [49] by Π ( q, t , t ) ≡ a s a ( η ) a ( η ) (cid:2) (¯˜ e + ¯˜ p ) ¯ U f (cid:3) L f ˜Π ( qL f , qη , kη ) , (2.24)where ¯˜ e and ¯˜ p are the rescaled average energy density and pressure, which correspond to the quantitiesmeasured at t s . The quantity ¯ U f describes the magnitude of the fluid velocity and is dimensionless.The correlator, Π , on the left hand side of the equation has dimension 5. Therefore, the additionallength factor L f is inserted here to make ˜Π dimensionless. Since this length scale is free from theeffect of the expanding universe, it is a comoving length scale. It is found from numerical simu-lations [49, 50] that the typical scale in the gravitational wave production is the (comoving) meanbubble separation R ∗ c . So we will choose L f = R ∗ c .The calculation of the UETC requires us to scrutinize the entire process of the phase transitionand the gravitational wave production. This task can be separated into two parts. The first part is astudy of the bulk parameters characterizing the process of the phase transition, which we will performin the next section. The second part is understanding the evolution of the source, which we go on toperform in Sec. 4. In this section, we study the changes to the dynamics of the phase transition in an expanding uni-verse. This includes parameters characterizing the behavior of the bubble formation, expansion andpercolation: the bubble nucleation rate, the fraction of the false vacuum and the unbroken area ofthe walls at a certain time. These will eventually be incorporated in the calculation of the veloc-ity power spectrum in the sound shell model. Another set of important quantities characterize thestatistics of the bubbles ever formed: the bubble lifetime distribution, as well as the bubble numberdensity. These are also needed in the velocity power spectrum calculation. Moreover, the timing ofsome important steps in the phase transition are also included, like the nucleation temperature andthe percolation temperature. Other changes to the parameters entering the gravitational wave powerspectrum calculation are also included, with β/H a representative example. We now proceed to adetailed discussion of these quantities.
The first and most basic ingredient in the analysis of a first order cosmological phase transition is thenucleation rate of the bubbles in the meta-stable vacuum at finite temperature [89, 90]. The numberof bubbles nucleated per time per physical volume is given by the following formula: p = p exp (cid:20) − S ,b ( T ) T (cid:21) . (3.1)Here S is the Euclidean action of the underlying scalar field (cid:126)φ that minimizes the solution S ( (cid:126)φ, T ) = 4 π (cid:90) drr (cid:32) d(cid:126)φ ( r ) dr (cid:33) + V ( (cid:126)φ, T ) , (3.2)with the following bounce boundary conditions: d(cid:126)φ ( r ) dr (cid:12)(cid:12)(cid:12) r =0 = 0 , (cid:126)φ ( r = ∞ ) = (cid:126)φ out , (3.3)– 8 – Figure 2 . The representative profile of S ( T ) /T for the example used in Sec. 3. See Appendix. A for detailson how to reproduce this. where (cid:126)φ out are the components of the vev outside the bubble. For the pre-factor, we see that p ∝ T on dimensional grounds, while its precise determination requires integrating out fluctuations aroundthe bounce solution (see e.g., [91, 92] for detailed calculations or [93] for a pedagogical introduction).The function S ( T ) /T generally starts from infinity at T c and drops sharply as temperaturedecreases, with a typical profile shown in Fig. 2. Bubbles will be nucleated within a short range oftime, say at t ∗ , when this rate changes slowly, which admits the following Taylor expansion: p ( t ) = p exp [ − S ∗ + β ( t − t ∗ )] , (3.4)where S ∗ ≡ S ( T ∗ ) /T ∗ , β ≡ d ln p ( t ) /dt | t = t ∗ . More explicitly, we have S T = S T (cid:12)(cid:12)(cid:12)(cid:12) t ∗ + d ( S /T ) dT dTdt (cid:12)(cid:12)(cid:12)(cid:12) t = t ∗ (cid:124) (cid:123)(cid:122) (cid:125) ≡− β ( t − t ∗ ) , (3.6)and thus βH ∗ = − H ∗ dTdt d ( S /T ) dT (cid:12)(cid:12)(cid:12)(cid:12) t = t ∗ . (3.7) If there exists a barrier at zero temperature, then S ( T ) /T will reach a minimum, say at t ∗ . The rate can be expandedaround the minimum: p ( t ) = p exp (cid:20) − S ∗ − β ( t − t ∗ ) (cid:21) , (3.5)with β ≡ S (cid:48)(cid:48) ( t ∗ ) and the first derivative vanishing. The bubble nucleation will happen mostly around t ∗ , making it looklike an instantaneous nucleation [52]. – 9 –e will later see how t ∗ should be chosen. For now, we provide a generic expression for β duringan expanding universe, which needs the relation between t and T . Suppose the universe is expandingas a = c a t n and the radiation sector is expanding adiabatically such that entropy s R is conserved percomoving volume for the radiation sector: s R ( T ) a = const . (3.8)Here s R ∝ T , giving then T ∝ /a ∝ t − n . This is the case for a radiation dominated universe,and for a matter dominated universe where the non-relativistic matter does not inject entropy to theradiation sector. However when the matter decays into radiation, entropy injection into the radiationsector gives a different dependence T ∝ a − / [94]. Generically, we can assume T ∝ a − γ , (3.9)which then leads to T = c T t − nγ , with c T being another constant. We thus have dTdt = − c T nγ t − nγ − . (3.10)Moreover H = ˙ a/a = n/t . Then H dTdt = − c T γ t − nγ = − γ T. (3.11)Therefore β/H ∗ reduces to the following form βH ∗ = γ T d ( S /T ) dT (cid:12)(cid:12)(cid:12)(cid:12) t = t ∗ . (3.12)It is obvious from this result that β/H ∗ does not depend on n , i.e., it does not depend on how thescale factor evolves with time but rather on how T decreases with the scale factor through γ . For boththe standard radiation dominated universe and an early matter dominated universe wherein the matteris decoupled from the radiation, γ = 1 . For the matter dominated universe wherein the matter decaysinto radiation, γ = 3 / , which gives a smaller β/H ∗ . The false vacuum fraction g ( t c , t ) at t > t c can be obtained following the derivation in Ref. [95] g ( t c , t ) = exp (cid:20) − π (cid:90) tt c dt (cid:48) p ( t (cid:48) ) a ( t (cid:48) ) r ( t (cid:48) , t ) (cid:21) ≡ exp [ − I ( t )] . (3.13)Here I ( t ) corresponds to the volume of nucleated bubbles per comoving volume, double counting theoverlapped space between bubbles and virtual bubbles within others. r ( t (cid:48) , t ) is the comoving radiusof the bubble nucleated at t (cid:48) and measured at t , r ( t (cid:48) , t ) = (cid:90) tt (cid:48) dt (cid:48)(cid:48) v w a ( t (cid:48)(cid:48) ) . (3.14)For Minkowski spacetime, r ( t (cid:48) , t ) = v w ( t − t (cid:48) ) . For a FLRW spacetime r ( t (cid:48) , t ) = v w ( η (cid:48) − η ) , whichtakes the same form as the Minkowski spacetime, irrespective of the detailed expansion behavior,when conformal time is used. In obtaining the above results, a constant bubble wall velocity v w has– 10 – - - - Figure 3 . The false vacuum fraction as defined Eq. 3.13 for different fractions of matter energy density at T c ( κ M = 0 , . , defined in Eq. 3.18) and for several bubble wall velocities ( v w = 0 . , . , . ). The case of κ M = 0 corresponds to a radiation dominated universe and κ M = 0 . for matter domination. The horizontalline at g = 0 . is roughly the time when the bubbles percolate. been assumed and the initial size of the bubble has been neglected. This is justified as the initial sizeis very small.Eq. 3.13 can be recast in a form that is convenient for calculations, in terms of the temperature.Suppose that the scale factor at the time of the critical temperature is a c and that the scale factor at alater time is related to it by aa c ≡ (cid:18) T c T (cid:19) /γ . (3.15)The comoving bubble radius can be conveniently expressed with an integral over temperature: r ( T (cid:48) , T ) = v w a c (cid:90) T (cid:48) T dT (cid:48)(cid:48) T (cid:48)(cid:48) γH ( T (cid:48)(cid:48) ) (cid:18) T c T (cid:48)(cid:48) (cid:19) − /γ . (3.16)Accordingly I ( T ) can be written as I ( T ) = 4 π (cid:90) T c T dT (cid:48) T (cid:48) γH ( T (cid:48) ) ¯ p T (cid:48) exp (cid:20) − S ( T (cid:48) ) T (cid:48) (cid:21) (cid:18) T c T (cid:48) (cid:19) /γ [ a c r ( T (cid:48) , T )] . (3.17)– 11 –ere the ¯ p is defined by p = ¯ p T and we choose ¯ p = 1 in the examples of analysis as is usuallydone in the literature. A different choice of ¯ p would, of course, affect the resulting false vacuumfraction and thus the relevant temperatures defined [96]. Since the focus here is on the changes dueto different expansion histories, a fixed choice of ¯ p serves our purpose well. For the Hubble rate, weneed to be more precise with regard to the matter content. We consider a universe consisting of bothradiation and non-relativistic matter and define κ M to be the fraction of the total energy density at T c that is non-relativistic matter: κ M = ρ Matter ρ Total (cid:12)(cid:12)(cid:12)(cid:12) T = T c . (3.18)We also neglect the vacuum energy for these examples, though it certainly exists during a phasetransition. H = H ( T c ) (cid:114) κ M y + 1 − κ M y , (3.19)where y = a/a ( T c ) . We show in Fig. 3 the false vacuum fraction during the phase transition, fora purely radiation dominated universe with κ M = 0 and a matter dominated one with κ M = 0 . ,and for three choices of bubble wall velocities v w = 0 . , . , . . For both choices of κ M , it is clearfrom these figures that increasing v w speeds up the process of phase transition. From κ M = 0 to κ M = 0 . , a larger energy density and thus a larger Hubble rate is obtained, which decreases thefunction r ( T (cid:48) , T ) and I ( I ) and thus slows down the drop of g ( T c , T ) .One often encounters the percolation temperature, which is defined such that the fraction in truevacuum is about of the total volume [86], i.e., when g ( T c , T p ) ≈ . , or I ( T p ) ≈ . , (3.20)and corresponds to the intersection points of the horizontal line with the curves in Fig. 3. Sincedifferent choices of v w and κ M lead to different g ( T c , T ) , the corresponding values of T p are alsodifferent. With the false vacuum fraction in Eq. 3.13, the unbroken bubble wall area during the phase transitioncan be derived [53] and will be used in the derivation of the bubble lifetime distribution. Consider acomoving volume of size V c and a sub-volume occupied by false vacuum V c, False . Then the comovingunbroken bubble wall area A c ( t ) at t satisfies the following relation: dg ( t , t ) = dV c, False V c = −A c ( t ) v w dta ( t ) = −A c ( t ) v w dη. (3.21)Then A c is given by A c ( t ) = − v w dg ( t , t ) dη = a ( η ) H ( T ) γTv w dg ( T c , T ) dT . (3.22)One can also define the proper area per proper volume A = Proper AreaProper Volume = a × Comoving Area a × Comoving Volume = 1 a A c . (3.23)– 12 – - - - - - - - - Figure 4 . The dimensionless comoving uncollided bubble wall area as defined in Eq. 3.22 and Eq. 3.41 fordifferent values of κ M (defined in Eq. 3.18) and v w . Since A c ( t ) and A are the area per volume, they are of dimension , and can be presented in unitsof m − or GeV. A more meaningful representation can be obtained by comparing it with the typicalscale at the corresponding temperature. One such quantity is β c , to be defined later, which is thecomoving version of the β parameter and is related to the mean bubble separation (also to be definedlater). We show A c /β c in Fig. 4 for different choices κ M and v w , similar to what are used in Fig. 3.We can see the area first increases as more bubbles are formed and expanding. It decreases as bubblescollide with each other and the remaining false vacuum volume is shrinking to zero. The differentbehaviors when changing v w and the amount of non-relativistic matter contents coincide with whatwe observe in Fig. 3. The bubble lifetime distribution describes the distribution of bubble lifetime for all the bubbles everformed and destroyed during the entire process of the phase transition. This can be obtained with thehelp of the unbroken bubble wall area derived earlier, by generalizing the result derived in Ref. [53]to the expanding universe. We start by considering the number of bubbles that are created at t (cid:48) andare destroyed with comoving radius r . Here a bubble is defined as destroyed when approximatelyhalf of its volume is occupied by the expanding true vacuum space. These bubbles are therefore ata comoving distance of r at t (cid:48) from the part of the unbroken bubble wall, assuming constant anduniversal bubble wall velocity v w . The time t when this set of bubbles is destroyed is connected with– 13 – igure 5 . Illustration for the calculation of the bubble lifetime distribution. At t (cid:48) , there is a central blue blobcomposed of two already collided bubbles depicting a region of true vacuum space which is expanding intothe surrounding false vacuum space, and also a small red nucleus denoting a bubble starting to form. At thistime, the comoving distance between the red dot and the nearest blue boundary is r . At t fc , the walls of theblue blob and the fledged red bubble advance to the place denoted by blue and red dashed circles respectively,where they make the first contact. At t , they reach the place denoted by the solid blue and red circles, wherehalf of the red bubble is devoured by the blue one, and the red bubble is defined to be destroyed with a finalradius r . t (cid:48) and r by r = (cid:90) tt (cid:48) v w dt (cid:48)(cid:48) a ( t (cid:48)(cid:48) ) . (3.24)Since only two quantities out of ( r, t (cid:48) , t ) are independent, we denote A c ( t ( t (cid:48) , r )) as A c ( t (cid:48) , r ) anddefine the number of bubbles per comoving volume as n b,c . We then have (see an illustration andmore details in Fig. 5): d n b,c = p ( t (cid:48) ) (cid:2) a ( t (cid:48) ) A c ( t (cid:48) , r ) dr (cid:3) dt (cid:48) . (3.25)This implies that: d (cid:18) dn b,c dr (cid:19) ≡ dn b,c ( r ) = p ( t (cid:48) ) (cid:2) a ( t (cid:48) ) A c ( t (cid:48) , r ) (cid:3) dt (cid:48) . (3.26)Now, for fixed r , we consider all the bubbles ever formed before a time t f : n b,c ( r ) | t f t c = (cid:90) t f t c dt (cid:48) p ( t (cid:48) ) (cid:2) a ( t (cid:48) ) A c ( t (cid:48) , r ) (cid:3) , (3.27)and n b,c ( r ) = 0 at t c for all r . Consider a time when all bubbles have disappeared, when t f is largeenough. Now n b,c ( r ) | t f becomes a constant ˜ n b,c ( r ) . We can then relate r with the lifetime of the– 14 –ubbles. For the bubble nucleated at t (cid:48) and destroyed at t , we have r = (cid:90) tt (cid:48) dt (cid:48)(cid:48) v w dt (cid:48)(cid:48) a ( t (cid:48)(cid:48) ) = v w η lt , (3.28)where η lt is the conformal lifetime of the bubble. Thus, r has the same relation with the confor-mal lifetime as its relation with t lt in Minkowski spacetime. We can therefore proceed to derive aconformal lifetime distribution for all bubbles ever formed and destroyed: ˜ n b,c ( η lt ) ≡ dn b,c dη lt = v w ˜ n b,c ( r ) = v w (cid:90) t f t c dt (cid:48) p ( t (cid:48) ) a ( t (cid:48) ) A c ( t (cid:48) , v w η lt ) . (3.29)Remember A c ( t (cid:48) , v w η lt ) = A c ( t ( t (cid:48) , v w η lt )) and it is evaluated at t , which should be determinedthrough Eq. 3.28 given t (cid:48) and η lt . To present a numerically convenient representation of the aboveresult, we convert coordinate time t to conformal time η and then to temperature. For the bubbleformed at t (cid:48) , the corresponding conformal time is related to temperature by η (cid:48) − η c = (cid:90) t (cid:48) t c dt (cid:48)(cid:48) a ( t (cid:48)(cid:48) ) = 1 a c (cid:90) T c T (cid:48) dT (cid:48)(cid:48) T (cid:48)(cid:48) γH ( T (cid:48)(cid:48) ) (cid:18) T c T (cid:48)(cid:48) (cid:19) − /γ ≡ ∆ η ( T (cid:48) , T c ) . (3.30)Then for the bubble with conformal lifetime η lt , the conformal time for its destruction is given by η lt + ( η (cid:48) − η c ) , with the corresponding temperature T determined through η lt + ( η (cid:48) − η c ) = ∆ η ( T, T c ) . (3.31)This temperature, or time, is what should be used in A c , rather than T (cid:48) . With the relation between T and T (cid:48) found, it is then straightforward to do the integral in Eq. 3.29, which requires only converting t (cid:48) to temperature. The evolution of the bubble number density per proper volume n b = N b /V is governed by thefollowing equation d [ n b a ( t )] dt = p ( t ) g ( t c , t ) a ( t ) , (3.32)which can be integrated to give (noting that n b ( t c ) = 0 ): n b ( t ) = 1 a ( t ) (cid:90) tt c dt (cid:48) p ( t (cid:48) ) g ( t c , t (cid:48) ) a ( t (cid:48) ) . (3.33)This does not include the decrease of bubble number due to collisions and n b thus includes all thebubbles ever formed. The result for n b ( t ) can be similarly transformed into a function of temperature. n b ( T ) = (cid:18) TT c (cid:19) /γ (cid:90) T c T dT (cid:48) T (cid:48) γH ( T (cid:48) ) ¯ p T (cid:48) exp (cid:20) − S ( T (cid:48) ) T (cid:48) (cid:21) g ( T c , T (cid:48) ) (cid:18) T c T (cid:48) (cid:19) /γ . (3.34)We show n b in units of m − in the left panel of Fig. 6 and the total bubble number per Hubble volume n b /H ( T ) in the right panel. We can see that the bubble number density increases for a delayed false– 15 – Figure 6 . The number of bubbles (see Eq. 3.33) per m (left) and per Hubble volume(right) as a function oftemperature for difference fractions of non-relativistic matter content at the critical temperature κ M (definedin Eq. 3.18) and for different bubble wall velocities v w . vacuum fraction, which is consistent with physical intuition. From n b , we can define the mean bubbleseparation, R ∗ , as R ∗ ( t ) = (cid:20) V ( t ) N b ( t ) (cid:21) / = (cid:20) n b ( t ) (cid:21) / . (3.35)This is shown in Fig. 7. For both n b and R ∗ , it appears they both reach an asymptotic value after thebubbles have disappeared when the curves in these figures become flat. This is misleading as afterthe time the bubbles have disappeared, n b will be diluted as /a and accordingly R ∗ increases as a . The flat curves in the figures are simply due to the very tiny change of temperature plotted. Fromnumerical simulations [49, 50], it is found that the peak frequency of the gravitational wave spectrumis related to R ∗ . Therefore any change on R ∗ will translate into a shift of the peak frequency ofthe gravitational waves. Since R ∗ is of particular importance, it is convenient to use the comovingversion of it R ∗ c = ( V c /N b ) / , which will reach an asymptotic value after the bubble disappearance.From the right panel of Fig. 6, we can easily read off the nucleation temperature T n , which isdefined such that at this temperature there is about one bubble within a Hubble volume [97]. Note T n obtained this way differs slightly from the usually used, and a bit crude, criterion: (cid:90) t n t c dt p ( t ) H ( t ) = 1 , (3.36)which for radiation dominated universe where H = 8 πGρ/ and ρ = π g ∗ ( T ) T translates intothe condition: (cid:90) T c T n dTT (cid:18) π g ∗ (cid:19) (cid:16) m Pl T (cid:17) exp (cid:20) − S ( T ) T (cid:21) = 1 . (3.37)– 16 – - - - - - - - Figure 7 . Mean bubble separation R ∗ (defined in Eq. 3.35) for difference fractions of the non-relativisticmatter content at the critical temperature κ M and for different bubble wall velocities v w . The left panel is inunit of meter and the right in unit of Hubble radius. Here m Pl is the Planck mass. A further simplification says that T n is determined by S ( T n ) /T n =140 [97]. These determined T n differs slightly from the more accurate result obtained by solvingdirectly for n b with Eq. 3.32. β and Mean Bubble Separation ( R ∗ ) It was found from numerical simulations that the peak of the gravitational wave power spectrum islocated at kR ∗ ∼ [50], where R ∗ is the mean bubble separation defined earlier. However thestandard spectrum people generally use is expressed in terms of β (see, e.g., [1, 2]). So the relationbetween β and R ∗ is needed. It can be derived analytically under reasonable assumptions as wasshown in Ref. [53], which says R ∗ = (8 π ) / β ( v w ) v w . (3.38)Here we emphasize that β varies when v w is changed. The question is then will this relation still holdin an expanding universe. We will answer this question by giving a detailed derivation here, whichparallels and generalizes the derivation in Ref. [53].We rewrite Eq. 3.32 in terms of the conformal time (we still use the same function labels though t is replaced by η ) d ( n b,c ) dη = p ( η ) g ( η c , η ) a ( η ) , (3.39)where n b,c = n b a and is the comoving bubble number density. Here the false vacuum fraction g decreases sharply when its exponent I ( T ) becomes of order . Since p ( η ) increases exponentially,there is a peak for the r.h.s in above equation, at which time the bubbles are mostly nucleated. As g decreases much more sharply than p increases, the rate p only increases slowly during this time– 17 –uration and it can be Taylor expanded at around this time. This time can be conveniently chosen tobe η which satisfies I ( η ) = 1 . Then similarly to Eq. 3.4, we define a Taylor expansion but w.r.t theconformal time: p ( η ) = p ( η ) exp [ − S + β c ( η − η )] , (3.40)where we have neglected the very slow change of p ( η ) and defined a comoving version of η : β c = d ln pdη (cid:12)(cid:12)(cid:12)(cid:12) η = η . (3.41)Now lets see how n b,c in Eq. 3.39 can be solved in terms of β c . To do it, lets firstly see how g or itsexponent I can be expressed in terms of β c . From Eq. 3.13, we can write I in the following way: I ( η ) = 4 π (cid:90) ηη c dη (cid:48) a ( η (cid:48) ) p ( η (cid:48) ) r ( η (cid:48) , η ) = 4 π v w (cid:90) ηη c dη (cid:48) p ( η ) e − S + β c ( η (cid:48) − η ) ( η − η (cid:48) ) = 8 π v w β c p ( η ) e − S + β c ( η − η ) . (3.42)Now define a time η f such that I ( η f ) = 1 , (3.43)then at a later time much simpler expressions can be obtained: I ( η ) = e β c ( η − η f ) , g ( η c , η ) = e − I ( η ) . (3.44)As I ( η ) depends on the bubble wall velocity v w , the resulting t f and more importantly β c is a functionof v w . Plugging above expressions of g ( η c , η ) , p ( η ) into Eq. 3.39, and integrating over η , we have n b,c = 1 β c p ( η ) e − S + β c ( η f − η ) = β c ( v w )8 πv w . (3.45)Here the second equality comes from the relation in Eq. 3.43. As noted in Ref. [53], the best choiceof t is t f so that the Taylor expansion of p ( η ) converges more quickly. This result gives the relationbetween the comoving mean bubble separation R ∗ c and η c : R ∗ c = (8 π ) / v w β c ( v w ) . (3.46)We can also write all results in terms of physical quantities. From Eq. 3.41 and enforcing t = t f ,we have β c = a ( η f ) β = a ( η f ) (cid:20) γH ( T ) T d ( S /T ) dT (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) T = T f , (3.47)and thus n b,c = a ( η f ) β πv w . (3.48)– 18 – .0 0.2 0.4 0.6 0.8 1.0024681012 Figure 8 . The left panel shows the mean bubble separation R ∗ immediately after all the bubbles havedisappeared versus bubble wall velocity v w . The right panel shows β ( v w ) calculated using Eq. 3.7 at t f , ascompared with that calculated using R ∗ using Eq. 3.50. The dotted line shows these differ by roughly . Note n b,c becomes a constant number as N b reaches its maximum and the comoving volume is fixed.The physical number density after all the bubbles have vanished will be diluted by the expansion.Suppose we consider the physical number density n b at time η , then n b ( η ) = (cid:18) a ( η f ) a ( η ) (cid:19) β πv w . (3.49)The corresponding physical mean bubble separation would be R ∗ ( η ) = a ( η ) a ( η f ) (8 π ) / v w β ( v w ) . (3.50)Therefore the relation between R ∗ and β is similar to that derived in Minkowski spacetime and needsonly additional attention on the scale factors. If one uses R ∗ c and β c , then the relation is exactly thesame as in Minkowski spacetime. We emphasize again that β and β c are functions of v w . To see this,we plot R ∗ at a time immediately after all the bubbles have disappeared, as a function of v w , in theleft panel of Fig. 8. For each v w , we find the corresponding β ( v w ) as implied in above equation andcompare with β ( v w ) directly calculated using Eq. 3.47. This comparison is shown in the right paneland the two different determinations differ by at most , where the uncertainty can be attributed tothe approximations made. The dominant source of gravitational wave production is the sound waves in a perturbed plasmadue to the advancing bubble walls and their interaction with the surrounding fluid. In the soundshell model [52, 53], the total velocity field is modelled as a linear superposition of the individualcontribution from each bubble. The first step is then to understand the velocity profile of the fluidaround a single bubble. This topic has been extensively studied several decades ago and is reviewed– 19 –ith a complete treatment in Ref. [51]. However the analysis is set in Minkowski spacetime and itis not clear whether it needs changes in an expanding universe. Ref. [98] studied the velocity profilein an expanding universe and found that there is a significant change to the velocity profile and areduction of energy fraction going into the kinetic energy of the sound waves. But we will see inthis section the velocity profile actually remains unchanged. We first review the full set of fluid andfield equations and then analyze the fluid velocity profile around a single bubble. Armed with thisinformation, we then find the total velocity field from a population of bubbles in the sound shellmodel and calculate the velocity field power spectrum.
Numerical simulations that are performed to provide the widely adopted gravitational wave formulaeare based on the fluid-order parameter field model [49, 99, 100] in Minkowski spacetime. Here wegeneralize the full set of equations used in the simulations to the FLRW universe. Our purpose isto understand whether simulations can be done in Minkowski spacetime and then generalized to anexpanding universe by simple rescalings of the physical quantities. This is an important question asit is computationally very expensive to do a numerical simulation.The universe consists of: (1) the underlying scalar field(s) responsible for the phase transition;(2) the relativistic plasma whose constituent particles can interact with the scalar field(s); (3) magneticfield produced from the phase transition; (4) other sectors which do not directly interact with eitherthe scalar field, the plasma or the magnetic field, though they do interact gravitationally. We willneglect (3) by focusing on the dominant source for gravitational wave production, and only consider(4) through its effect on the expansion. Given our cosmological context, the total energy momentumtensor for (1) and (2) is given by [49] T µν = ∂ µ φ∂ ν φ − g µν ∂ µ φ∂ µ φ + ( e + p ) U µ U ν + g µν p, (4.1)where U µ = γ (1 , v /a ) with γ = 1 / √ − v and v = d x /dη . The energy and momentum densityare given by e = a B T + V ( φ, T ) − T ∂V∂T ,p = 13 a B T − V ( φ, T ) , (4.2)where a B = g ∗ π / and g ∗ is the relativistic degrees of freedom. It is certainly conserved, i.e., T µν ; µ = 0 , and it is usually split into two parts by adding and subtracting a friction term δ ν [99]: T µν ; µ | field = ( ∂ φ ) ∂ ν φ + 1 √ g ( ∂ µ √ g )( ∂ µ φ )( ∂ ν φ ) − ∂V∂φ ∂ ν φ = δ ν ,T µν ; µ | fluid = ∂ µ [( e + p ) U µ U ν ] + (cid:20) √ g ( ∂ µ √ g ) g νλ + Γ νµλ (cid:21) ( e + p ) U µ U λ + g µν ∂ µ p + ∂V∂φ ∂ ν φ = − δ ν . (4.3)Note here the appearance of ∂ µ g and Γ νµλ as we are using a generic metric. The friction term δ ν ismodelled by δ ν = ηU µ ∂ µ φ∂ ν φ . For high temperatures it can be chosen as η = ˜ ηφ /T [50], whichworks well in that case [101] but may lead to numerical singularities for small temperature. The The subscript “;” denotes covariant derivative. – 20 –umerical simulations on sound waves adopted a constant value for the lower temperature case [102].The exact set of equations can also be derived from field theory [96, 103].In an FLRW universe, the field energy momentum conservation leads to a scalar equation: − ¨ φ + 1 a (cid:79) φ − ∂V∂φ − aa ˙ φ = ηγ ( ˙ φ + 1 a v · (cid:79) φ ) , (4.4)which is just the Klein-Gordon equation for the scalar field when the friction term is absent, i.e., when η = 0 . The vector part of the fluid energy-momentum conservation gives: ˙ Z i + 1 a (cid:79) · ( v Z i ) + 5 ˙ aa Z i + 1 a ∂ i p + 1 a ∂V∂φ ∂ i φ = − a ηγ ( ˙ φ + 1 a v · (cid:79) φ ) ∂ i φ, (4.5)where Z i ≡ γ ( e + p ) U i = γ ( e + p ) v i /a . The parallel projection along U ν for the fluid gives anotherscalar equation: ˙ E + p [ ˙ γ + 1 a (cid:79) · ( γ v )] + 1 a (cid:79) · ( E v ) − γ ∂V∂φ ( ˙ φ + 1 a v · (cid:79) φ ) + 3 ˙ aa γ ( e + p )= ηγ ( ˙ φ + 1 a v · (cid:79) φ ) , (4.6)where E ≡ eγ . While the above equations form a complete set, the velocity profile is usually derivedfrom a different scalar equation, the perpendicular projection for the fluid along the direction ¯ U ν ,which is defined by ¯ U µ U µ = 0 , ¯ U µ ¯ U µ = 1 , (4.7)and takes the explicit form ¯ U µ = γ ( v, ˆ v i /a ) . This gives (cid:20) ˙ aa v + γ (cid:18) ˙ v + 12 a ˆ v · (cid:79) v (cid:19)(cid:21) ( e + p ) + v ˙ p + 1 a ˆ v · (cid:79) p + ∂V∂φ ( v ˙ φ + 1 a ˆ v · (cid:79) φ )= − ηγ ( v ˙ φ + 1 a ˆ v · (cid:79) φ )( ˙ φ + 1 a v · (cid:79) φ ) . (4.8)These equations are direct generalizations of those in Ref. [49] to an FLRW universe. It is notpossible, however, to express the above equations in a form used in Minkowski spacetime and theproblem lies with the scalar field. Despite this, the effect on the bubble and fluid motions should beminor, since the bubble collision process is fast compared with the long duration of the ensuing soundwaves.The process of the phase transition can thus be divided into two stages. The first stage is thebubble collision and disappearance of the symmetric phase, and the second is the propagation ofsound waves. The difference between them is that the first stage takes a much shorter time, while thesecond is long-lasting. This is indeed what is observed from numerical simulations and should welljustify simply neglecting the change of the scale factor during the first stage [49]. In this sense, thenumerical simulations as performed in Ref. [49, 50] still give a faithful account of the first step for anexpanding universe. However we will see in the next subsection that the analytical modelling of thisfirst stage still admits simple rescaling properties and takes the same form as in Minkowski.During the second stage gravitational waves are dominantly produced due to the long-lastingnature of the sound waves. Therefore the change of the scale factor cannot be ignored. The questionis: can we still only perform numerical simulations in Minkowski spacetime. Fortunately, during thisstage, the scalar field plays no dynamical role and we can consider only the fluid. The corresponding– 21 –quations can indeed be reduced to the Minkowski form. This is achieved by using the conformaltime, neglecting the scalar field as well as the friction terms and using p = e/ for the plasma. ThenEq. 4.5, Eq. 4.6 and Eq. 4.8 reduce to (again, (cid:48) ≡ ∂/∂η ): ( a S i ) (cid:48) + (cid:79) · ( a S i v ) + ∂ i ( a p ) = 0 , ( a eγ ) (cid:48) + [ γ (cid:48) + (cid:79) · ( γ v )]( a p ) + (cid:79) · ( a eγ v ) = 0 ,γ ( v (cid:48) + 12 ˆ v · (cid:79) v )[ a ( e + p )] + v ( a p ) (cid:48) + ˆ v · (cid:79) ( a p ) = 0 , (4.9)where S i = aZ i = γ ( e + p ) v i . The Minkowski counterpart of these equations can be obtained bysetting a = 1 . This suggests that we can define rescaled quantities ˜ e = a e/a s and ˜ p = a p/a s ,where a s is the scale factor when the source becomes active. They are free from the dilution dueto the expansion, and that the equations governing ˜ e , ˜ p and v take exactly the same form as theirMinkowski counterparts, as long as the time t is interpreted as the conformal time η . We will seehow these rescaled quantities can be used to derive the modified gravitational wave spectrum in latersections.We note here that these equations were derived earlier in Ref. [104, 105] when also consideringelectromagnetism and it was shown that the above rescaling works not only for the purely fluid systembut also for a system containing both fluid and electromagnetism. Including electromagnetism willadd additional terms to the right hand side of the above equations. Solving the velocity profile for a single expanding bubble depends on analyzing the behavior of thesystem consisting of both the fluid and the scalar field. This is usually done in the so called bagequation of state model, as summarized in Ref. [51]. The energy momentum tensor for the fluid plusscalar field system is assumed to take the following form (“ + ” for outside the bubble and “ − ” forinside): T µν ± = p ± g µν + ( p ± + ρ ± ) U µ U ν , (4.10)with the bag equation of state: p + = 13 a + T − (cid:15), e + = a + T + (cid:15).p − = 13 a − T − , e − = a − T − , (4.11)where (cid:15) is the vacuum energy difference between the false and true vacua. One can also find theenthalpy ω = e + p . Here v , T and thus e, p, ω all vary from the bubble center to the region faroutside the bubble where there is no perturbation. The task is to solve for these fields at regions bothinside and outside the bubbles and smoothly match these two sets of solutions through the junctionconditions across the bubble wall. In this region, we drop all terms related to φ including the vacuum energy from (cid:15) , and we also applythe relation p = e/ . The resulting equations are already given in Eq. 4.9 and the equations areexactly the same as the Minkowski counterpart when the rescaled quantities are used. Now, assuming Of course, we are assuming a constant value of the speed of sound, i.e., c s = 1 / √ . Without doing so, the equationscannot be put into the form in Eq. 4.9. We also dropped any spatial variation of the scalar field and its time variationfollowing the conventional analysis. – 22 – spherically symmetric profile and denoting the comoving bubble radius with r and the conformaltime elapsed since its nucleation as ∆ η , the solution should be a self-similar one which dependssolely on the ratio ξ ≡ r/ ∆ η . Then we can obtain the same equations as in Minkowski spacetime: ( ξ − v ) ∂ ξ ˜ e = ˜ w (cid:20) vξ + γ (1 − ξv ) ∂ ξ v (cid:21) , (1 − vξ ) ∂ ξ ˜ p = ˜ wγ ( ξ − v ) ∂ ξ v, (4.12)which can then be combined to give an equation for the velocity field: vξ = γ (1 − vξ ) (cid:20) µ c s − (cid:21) ∂ ξ v. (4.13)Here µ ( ξ, v ) = ( ξ − v ) / (1 − ξv ) , which is the Lorentz boost transformation. This equation can bedirectly solved given a boundary condition at the wall, to be specified later. Outside the bubble, the presence of the constant vacuum energy term (cid:15) seemingly does not allow usto reach Eq. 4.9 for two possible reasons: (1) we can not apply p = e/ since p = − e for vacuumenergy; (2) (cid:15) does not scale like radiation with the behavior /a and the rescaled quantity a e stillcontains the expansion effect. Let us look more closely at the equations. The parallel projection inEq. 4.6, when the friction and scalar gradient terms are neglected, becomes (cid:20) ( γe ) (cid:48) + 3 a (cid:48) a γ ( e + p ) (cid:21) + p [ γ (cid:48) + (cid:79) · ( γ v )] + (cid:79) · ( γe v ) = 0 . (4.14)Correspondingly, the perpendicular projection in Eq. 4.8 reduces to (cid:20) a (cid:48) a v ( e + p ) + vp (cid:48) (cid:21) + γ ( v (cid:48) + 12 ˆ v · (cid:79) v )( e + p ) + ˆ v · (cid:79) p = 0 . (4.15)In the absence of the vacuum energy inside e and p , both of above equations can be put into the formin Eq. 4.9, by combining the terms in [ · · · ] and using e = 3 p . The resulting equations for the rescaledquantities are the same as in Minkowski spacetime. The presence of (cid:15) makes this impossible. InRef. [98], the self-similar velocity profile is assumed anyway. But the existence of an explicit timedependence from a (cid:48) makes it impossible to solve, except in corners of the parameter space where itvanishes numerically. It is also in doubt if there exists a self-similar solution at all for these equationsand we refrain from going in that direction.Despite this dilemma, we can still cast above equations in the form 4.9 under the assumptionthat (cid:15) is a constant of time during this very short period of time. Then the first equation can bereorganized in the following way: (cid:20) γ (cid:48) + (cid:79) · ( γ v ) + 3 a (cid:48) a γ (cid:21) ( e + p ) + γe (cid:48) + γ v · (cid:79) e = 0 . (4.16)Then (cid:15) cancels out in ( e + p ) and drops out in e (cid:48) , and of course also in (cid:79) e . So above e and p can include only the fluid part. Then one can put it back into the previous form 4.14 and definethe rescaled quantities: ˜ e , ˜ p , which obey exactly the same equation as in the Minkowski spacetime.Therefore we obtain the second equation in Eq. 4.9 and the first in Eq. 4.12. Similarly for Eq. 4.15, (cid:15) drops out in all terms and one can safely define the rescaled quantities, and obtain the third equationin Eq. 4.9 and the second in Eq. 4.12. Combining these two equations again gives the same equationfor the velocity field as in Eq. 4.13. – 23 – .2.3 Matching at Bubble Wall The equation 4.13 for both regions needs the junction conditions at the wall to connect them. Theyare derived by integrating the conservation of energy momentum tensor across the bubble wall, whichgives in the wall frame (note + , − denote quantities at positions immediately outside and inside thewall) T rη + = T rη − , (4.17) T rr + = T rr − , (4.18)where v − and v + are both at wall frame. These two equations imply ( e + + p + ) v + γ = ( e − + p − ) v − γ − , (4.19) ( e + + p + ) v γ + p + = ( e − + p − ) v − γ − + p − . (4.20)Here both e ± and p ± are the ordinary energy density and pressure and include the vacuum energy (cid:15) . The reason is while they can be neglected away from the bubble wall due to the vanishing spa-tial gradient, they jump across the bubble wall and give non-negligible contributions to the aboveequations. The junction equations can be solved by making the change of variables v = tanh( φ ) and γ = cosh ( φ ) which, after simplifying, will yield two linear equations in cosh ( φ + ) and cosh ( φ − ) .The solution will give v + = (cid:115) ( p − − p + )( e − + p + )( e − − e + )( e + + p − ) ,v − = (cid:115) ( p + − p − )( e + + p − )( e + − e − )( e − + p + ) . (4.21)The product and ratio of v + and v − can further be found, v + v − = p + − p − e + − e − , v + v − = e − + p + e + + p − . (4.22)Plugging e ± , p ± as specified by the bag equation of state in Eq. 4.11 leads to v + v − = 1 − (1 − α + ) σ − α + ) σ , (4.23) v + v − = 3 + (1 − α + ) σ α + ) σ , (4.24)where α + and σ are defined by α + = (cid:15)a + T (cid:12)(cid:12)(cid:12)(cid:12) wall , σ = a + T a − T − (cid:12)(cid:12)(cid:12)(cid:12) wall . (4.25) α + characterizes the amount of vacuum energy released from the phase transition normalized by thetotal radiation energy density immediately outside the bubble (as denoted by the subscript “wall”).It is not the α usually used in phase transition analyses. Rather, its value should be solved from the Also we follow the conventional procedure by neglecting the time dependence of the various quantities. – 24 – .0 0.2 0.4 0.6 0.8 1.00.00.10.20.30.40.5 Deflagration Hybrid Detonation
Figure 9 . Representative velocity profiles surrounding the bubble walls. requirement that far from the bubble where the plasma is not perturbed (denote by ∞ ), the corre-sponding α + at ∞ matches α . The two equations in Eq. 4.24 can be solved for both r and v + to givetwo branch solutions for the velocity in the symmetric phase, v + = 11 + α + (cid:18) v − v − (cid:19) ± (cid:115)(cid:18) v − v − (cid:19) + α + 23 α + − . (4.26)Up to this point, the results for the velocity profile are exactly the same as in Minkowski space-time, but with the understanding that the time t is replaced by the conformal time η , v = d x /dη and ( e, p ) are replaced by (˜ e, ˜ p ) . We will not go into the details of the physics of above results butonly summarize the main features of the velocity profile relevant for this study and refer the reader toRef. [51] for a more detailed analysis.The fluid admits three modes of motion: deflagration, detonation and supersonic deflagration(also called hybrid) [106], with representative velocity profiles shown in Fig. 9. For deflagration, thevelocity inside the bubble vanishes and is only non-zero outside. Detonation is the opposite, withnon-zero velocity inside the bubble. Supersonic deflagration has non-zero velocity both inside andoutside the bubble. Therefore for deflagration, v − = v w which should be used in Eq. 4.26 to find v + , choosing a value of α + . This v + is Lorentz transformed to the plasma static frame to find v ( v w ) immediately outside the wall, which is then used as the boundary condition to solve for v ( ξ ) outsidethe wall. It might not consistently drop to zero, in which case a shock front is encountered and shouldbe determined. Beyond the shock v ( ξ ) = 0 . This gives a complete profile, but not yet the correctone, since a specific value of α + is used in above determination of the profile. This value needs to betuned such that α + = α far outside the bubble. For detonation, v + = v w and v − can be determinedfrom Eq. 4.26 with α + = α as outside the bubble the plasma is not perturbed. Then one can Lorentztransform v − to v ( v w ) immediately inside the wall and use it as a boundary condition to determinethe full profile. No inconsistency or shock front will be encountered in this case. For supersonicdeflagration, the condition v − = c s is the boundary condition used. Shock front can exist in this caseand should be treated similarly. We refer the reader for more details in Ref. [51].– 25 – .3 Velocity Field in the Sound Shell Model With the velocity profile surrounding a single bubble determined, we can now find the total velocityfield, as needed in Eq. 2.19. As we have already seen, in an expanding universe the equations ofmotion of the fluid are exactly the same as those in non-expanding Minkowski spacetime. This meansthat the equation of motion for the sound waves remain the same as its Minkowski counterpart, aslong as we replace t by η and interpret the velocity as obtained by differentiation with respect to theconformal time. So the procedure parallels that in Ref. [53].Lets start with the contribution from one bubble. Before it collides with another bubble, thevelocity profile is governed by equations given in previous sections. After the collision, the frictionvanishes and the velocity field starts freely propagating and becomes sound waves, with the speed ofsound c s . So we need to match the velocity profile surrounding this bubble with the velocity field atthe time when the friction vanishes. Before collision, we can Fourier decompose the velocity field as v i ( η < η fc , x ) = 12 (cid:90) d q (2 π ) (cid:2) ˜ v i q ( η ) e i q · x + ˜ v i ∗ q ( η ) e − i q · x (cid:3) , (4.27)with x being the comoving coordinate and q the comoving wavenumber. After collision, the velocityfield freely propagates as sound waves and admits the following decomposition: v i ( η, x ) = (cid:90) d q (2 π ) (cid:2) v i q e − iωη + i q · x + v i ∗ q e iωη − i q · x (cid:3) , (4.28)where ω = qc s . Since the plasma consists of relativistic particles, c s = 1 / √ . Here v i q is independentof η , different from ˜ v i q ( η ) .The task is then to find the contribution to v i q from ˜ v i q ( η ) at η fc . Since the equation governingthe sound waves is of second order, we need the following initial conditions: ˜ v i q ( η ) and ˜ v i (cid:48) q ( η ) at η fc .While one can obtain ˜ v i q ( η ) directly from the velocity profile in the previous section, one subtletyappears here for ˜ v i (cid:48) q ( η ) . As demonstrated in Ref. [53], the equation governing ˜ v i (cid:48) q ( η ) before thecollision relies on a force term from the scalar field, which disappears once the collision occurs. Sothe value ˜ v i (cid:48) q ( η ) calculated with this force (as was previously used in Ref. [52]) is different from thecorresponding value without it. It is the latter one that should enter the initial conditions for the soundwaves. In this case, rather than calculating ˜ v i (cid:48) q ( η ) from the velocity profile ˜ v i q ( η ) , we need to calculateit directly from the energy fluctuations: λ ( x ) = ˜ e ( x ) − ¯˜ e ¯˜ ω , (4.29)where a bar denotes averaged quantity and tilde denotes rescaled quantity. Similarly its Fouriercomponent ˜ λ q can be defined in analogy to Eq. 4.27, The equation for sound waves then follows: ˜ λ (cid:48) q + iq j ˜ v j q = 0 , ˜ v j (cid:48) q + c s iq j ˜ λ q = 0 . (4.30)Therefore ˜ v j (cid:48) q = − c s iq j ˜ λ q , and one needs to calculate v i ( η, x ) and λ ( η, x ) from the self-similarvelocity profile for one bubble. In coordinate space, the velocity profile for the n-th bubble can bewritten as v ( n ) ( η, x ) = ˆ R ( x ) v ( ξ ) , (4.31)– 26 –
10 20 30 40 - - - Figure 10 . The real (blue dotted), imaginary (red dashed) parts and absolute value (magenta solid) of A ( z ) (defined below Eq. 4.35) for v w = 0 . and α = 0 . . where R ( x ) ≡ x − x ( n ) , ξ ≡ | R ( n ) | /T ( n ) and T ( n ) ( η ) ≡ η − η ( n ) , with x ( n ) and η ( n ) the coordinateof the bubble center and the conformal time when the bubble is nucleated. Similarly for λ , as itis a scalar field, we can define λ ( η, x ) ≡ λ ( ξ ) . With the profile specified in coordinate space, thecorresponding Fourier coefficients can be obtained straightforwardly ˜ v j ( n ) q ( η fc ) = e − i q · x ( n ) ( T ( n ) ) i ˆ z j f (cid:48) ( z ) | η = η fc , ˜ λ ( n ) q ( η fc ) = e − i q · x ( n ) ( T ( n ) ) l ( z ) | η = η fc , (4.32)with z ≡ q T ( n ) and the two functions f ( z ) and l ( z ) given by f ( z ) = 4 πz (cid:90) ∞ dξ v ( ξ ) sin( zξ ) ,l ( z ) = 4 πz (cid:90) ∞ dξ ξ λ ( ξ ) sin( zξ ) . (4.33)Then the n-th bubble’s contribution to the Fourier coefficient of the sound waves is v j ( n ) q = 12 (cid:104) ˜ v j ( n ) q ( η fc ) + c s ˆ q j ˜ λ ( n ) q ( η fc ) (cid:105) e iωη fc , (4.34)and after using the explicit expression of the bubble profile, v j ( n ) q = i ˆ z j ( T ( n ) fc ) e iωη fc − i q · x ( n ) A ( z fc ) , (4.35)where A ( z fc ) = [ f (cid:48) ( z fc ) − ic s l ( z fc )] / , with an example shown in Fig. 10. Thus we have calculatedthe contribution to v i q from one bubble that is nucleated randomly. The randomness of this bubble is– 27 –eflected in its formation time, location, collision time and its radius. Since the radius at collision isfixed once its formation and collision times are given, we have three independent random variables.The velocity field after all bubbles have disappeared, can be assumed to be the linear additionof the contributions from all bubbles, which is the essence of the sound shell model [52, 53]. Supposethe total number of bubbles nucleated within a Hubble volume with comoving size V c is N b . Thenthe velocity field can be assumed, according to the sound shell model, to be given by v i q = N b (cid:88) n =1 v i ( n ) q . (4.36) As these N b bubbles are just one realization of the phase transition, the resulting v i q has a randomnature with it and follows a Gaussian distribution to a good approximation according to the centrallimit theorem . Randomness of this kind can be removed by doing an ensemble average of theproduct: (cid:104) v i q v j ∗ q (cid:105) , which is all needed for a Gaussian distribution. Now let us see how this is achieved.The N b bubbles can be separated into groups with the bubbles within each group sharing acommon formation and collision time. Then the only variable that is random across the bubbles ofone group, e.g., group g with N g bubbles, is the spatial locations of the bubbles when they form. Nowconsider group g . Its contribution to the correlator is (cid:104) v i q v j ∗ q (cid:105) g = ˆ q i ˆ q j [ T ( g ) fc ] A ( z ( n ) fc ) A ( z ( m ) fc ) ∗ e i ( ω − ω ) η ( g ) fc (cid:104) N g (cid:88) m,n =1 e i q · x ( m ) − iq · x ( n ) (cid:105) . (4.37)Here the order of the ensemble average and the summation can be switched. Since the ensembleaverage of each of these N g terms gives the same result and oscillatory cross terms vanish, we have (cid:104) N g (cid:88) m,n =1 e i q · x ( m ) − iq · x ( n ) (cid:105) = N g δ mn (cid:104) e i q · x ( m ) − iq · x ( n ) (cid:105) = N g V c (cid:90) d x ( ∗ ) e i ( q − q ) · x ( ∗ ) = N g V c (2 π ) δ ( q − q ) . (4.38)The constraint q = q removes the η ( g ) fc dependence, leading to a result solely dependent on theconformal lifetime of the bubble T ( g ) fc ≡ η lt but not their absolute formation or destruction time: (cid:104) v i q v j ∗ q (cid:105) g = ˆ q i ˆ q j η lt | A ( qη lt ) | N g V c (2 π ) δ ( q − q ) . (4.39)This result means that we can combine groups with the same η lt , and of course, different formationtime, by solely enlarging the value of N g . In the following we will simply stick to the group label“ g ”, though its definition is changed and now includes all bubbles with the same η lt . Restricting to asufficiently small region centered at η lt , the number N g is an still an infinitesimally small fraction of N b and can be written as N g = N b P ( η lt ) dη lt , (4.40) If there is a sufficiently large population of bubbles within this single volume, the summation of these contributionscan also remove the randomness, equivalent to an ensemble average. – 28 – - - - - Figure 11 . The dimensionless bubble lifetime distribution ν ( β c η ) defined in Eq. 4.45 and more explicitly inEq. 4.47. All previously used choices of κ, v w give the same blue line. The gray dashed line is the analyticallyderived result e − βt lt in Ref. [53]. where P ( η lt ) is the probability density for bubbles to have conformal lifetime in the range ( η lt , η lt + dη lt ) , thus with dimension and normalized by (cid:90) dη lt P ( η lt ) = 1 . (4.41)Adding the contributions from all the groups and noting that cross terms vanish due to the oscillatorybehavior, we have (cid:104) v i q v j ∗ q (cid:105) = ˆ q i ˆ q j (2 π ) δ ( q − q ) (cid:90) dη lt (cid:20) P ( η lt ) N b V c (cid:21) η lt | A ( qη lt ) | . (4.42)One can now identify the quantity in the square bracket as the conformal lifetime distribution definedin Eq. 3.29: P ( η lt ) N b V c = ˜ n b,c ( η lt ) . (4.43)Since P ( η lt ) is of dimension , it is convenient to define a dimensionless version of it ν : P ( η lt ) ≡ β c ν ( β c η lt ) . (4.44)Then ˜ n b,c ( η lt ) = β c R ∗ c ν ( β c η lt ) , (4.45)where R ∗ c is the asymptotic comoving mean bubble separation. Then we have– 29 – − kR ∗ c − − − − − − − P v ( k R ∗ c ) k k − Figure 12 . Representative velocity power spectrum calculated in the sound shell model for a weak phasetransition with α = 0 . and v w = 0 . . The bubbles are assumed to nucleate exponentially. The low andhigh frequency regimes follow the k and k − power law fits respectively (black solid lines). (cid:104) v i q v j ∗ q (cid:105) = ˆ q i ˆ q j (2 π ) δ ( q − q ) 1 R ∗ c β c (cid:90) d ˜ T ˜ T ν ( ˜ T ) | A ( q ˜ Tβ c ) | , (cid:124) (cid:123)(cid:122) (cid:125) ≡ P v ( q ) (4.46)with here ˜ T = β c η lt , and we have defined the spectral density P v ( q ) for the plane wave amplitude v i q .Lets write down the explicit expression for ν ( ˜ T ) . From Eq. 4.45 and 3.29, we have ν ( ˜ T ) = v w R ∗ c (cid:90) t f t c dt (cid:48) p ( t (cid:48) ) a ( t (cid:48) ) A c ( t (cid:48) , v w ˜ T /β c ) β c , (4.47)which can be directly used for numerical calculations once t (cid:48) is transformed to T (cid:48) as was done inprevious sections. The numerically calculated distribution for the examples we have been using isshown in Fig. 11. For all choices of κ, v w , the distributions are almost indistinguishable, shownas the blue curve, and it coincides with the gray dashed curve which denotes the distribution e − ˜ T ,derived analytically in Ref. [53]. With ν ( ˜ T ) obtained, the spectral density P v ( ˜ T ) can be calculatedstraightforwardly from its definition in Eq. 4.46.To calculate the velocity power spectrum, we need to evaluate the correlator (cid:104) ˜ v i q ( η )˜ v j ∗ k ( η ) (cid:105)| = δ ( q − k )ˆ q i ˆ k j G ( q, η , η ) , (4.48)and it can be shown that G ( q, η , η ) = 2 P v ( q ) cos[ ω ( η − η )] . (4.49)– 30 –lugging it into Eq. 2.21 or 2.22 gives the stress energy correlator. Also the velocity field powerspectrum P v follows naturally, P v = q π [2 P v ( q )]= 164 π v w ( qR ∗ c ) (cid:90) d ˜ T ˜ T ν ( ˜ T ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A (cid:32) ( qR ∗ c ) ˜ T (8 π ) / v w (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (4.50)and we have used β c R ∗ c = (8 π ) / v w . It is obvious to see that P v is dimensionless, as it is con-structed with purely dimensionless quantities. A representative profile for the velocity power spec-trum is shown in Fig. 12 assuming an exponential bubble nucleation rate, and more details about itsproperties can be found in Ref. [52]. We can now go back to Eq. 2.14 and collect all the pieces to calculate the gravitational power spec-trum. It only remains to calculate the Green’s function, and it requires to specify an expansionscenario. We will as usual focus on the RD and MD scenarios as examples, but the method here isapplicable to any expansion history.
First, we choose a parameter to measure the time of the cosmic history. It can either be the actualtime t , the conformal time η , the redshift z or the scale factor a . To present a result independent ofthe origin of the time coordinate, we choose the dimensionless scale factor ratio y ≡ a/a s , givingthen d/dt = ˙ a/a s d/dy . Here a s is the time when the source, the sound waves, becomes active, sothat y starts from . The Friedmann equation gives the relation between y and the conformal time y = κ M a s H s ) ( η − η s ) + a s H s ( η − η s ) + 1 . (5.1)It is obvious that when η = η s , we have y = 1 . Also it does not matter how the origin of theconformal time is chosen as it only depends on ∆ η ≡ η − η s . For RD, where κ M ∼ , we have y = a s H s ( η − η s ) + 1 . For MD, where κ M ≈ , y = [ a s H s ( η − η s ) + 1] . In the literature, itis usually approximated that a ∝ η deep inside the radiation era or a ∝ η deep inside the matterera. However we remain agnostic about when the phase transition happens and do not require it tostart deep inside the radiation or matter era. Also the duration of the phase transition is very smallcompared with the conformal time, which makes such approximation quite crude. But our choiceusing y is free from above limitations and offers a more accurate description of phase transitionprocess.With y , the Hubble rate, when assuming the existence of both matter and radiation components,takes the following form H = H s (cid:114) κ M y + 1 − κ M y , (5.2)where κ M is the matter fraction of the total energy density at t s . Note this κ M is defined differentlyfrom that in Eq. 3.19, which is defined at T c . If the lifetime of the sound waves is sufficiently long,we can neglect this difference. – 31 –witching from the conformal time η to y in Eq. 2.5, the Einstein equation becomes : ( κ M y + 1 − κ M ) d h q dy + (cid:20) κ M + 2(1 − κ M ) y (cid:21) dh q dy + (cid:101)(cid:101) q h q = 16 πGa ( y ) π Tq ( y )( a s H s ) . (5.3)Here (cid:101)(cid:101) q ≡ q/ ( a s H s ) , and characterizes the number of wavelengths contained within a Hubble radiusat t s . The Green’s function can be found by solving the homogeneous version of this equation,together with a slightly modified boundary conditions compared with Eq. 2.11: G ( y (cid:54) y ) = 0 , ∂G ( y, y ) ∂y | ˜ η =˜ y +0 = 1 κ M y + 1 − κ M . (5.4)The solution to the homogeneous equation is a linear combination of the hypergeometric functionand Bessel functions. For the case of radiation domination κ M (cid:28) and matter domination κ M ≈ ,the solutions take simpler forms that can be expressed in terms of elementary functions. For RD, theequation becomes simpler when expressed using the parameter ˜ y , defined by ˜ y = y (cid:101)(cid:101) q = q ( η − η s ) + (cid:101)(cid:101) q = ∆˜ η + (cid:101)(cid:101) q. (5.5)Then the Einstein equation becomes d h q d ˜ y + 2˜ y dh q d ˜ y + h q = 16 πGa ( y ) π Tq ( y ) q . (5.6)The corresponding Green’s function can be easily solved: G (˜ y, ˜ y ) = ˜ y sin(˜ y − ˜ y )˜ y . (5.7)For MD, the wave equation can be similarly simplified with ˜ y = y (cid:101)(cid:101) q = (cid:20)
12 ∆˜ η + (cid:101)(cid:101) q (cid:21) . (5.8)Note this definition is different from that in the radiation dominated case. Then the Einstein equationbecomes ˜ y d h q d ˜ y + 52 dh q d ˜ y + h q = 16 πGa (˜ y ) π Tq (˜ y ) q . (5.9)The homogeneous equation for h q can be transformed into the Bessel equation for a different variable Z ( λ ) defined by h q = ( λ/ − / Z ( λ ) with λ = 2 √ ˜ y : λ Z (cid:48)(cid:48) ( λ ) + λZ (cid:48) ( λ ) + (cid:34) λ − (cid:18) (cid:19) (cid:35) Z ( λ ) = 0 . (5.10)The two independent solutions are the first and second kind Bessel functions both with order / ,which can all be expressed in elementary functions. Upon using the boundary conditions, the Green’s We are using a simplified notation for h and π T – 32 –unction is found to be : G (˜ y, ˜ y ) = ( λλ + 1) sin( λ − λ ) − ( λ − λ ) cos( λ − λ ) λ / . (5.12)Finally in both cases, the gravitational wave amplitude is given by h ij (˜ y, q ) = (cid:90) ˜ y ˜ y s d ˜ y (cid:48) G (˜ y, ˜ y (cid:48) ) 16 πGa (˜ y (cid:48) ) π Tij (˜ y (cid:48) , q ) q . (5.13) The spectral density for h (cid:48) , when using ˜ y and the dimensionless stress energy tensor correlator ˜Π defined in Eq. 2.24, becomes P h (cid:48) = [16 πG (¯˜ (cid:15) + ¯˜ p ) ¯ U f ] L f (cid:90) ˜ y ˜ y s d ˜ y (cid:90) ˜ y ˜ y s d ˜ y (cid:18) ∂ ˜ y∂ ˜ η (cid:19) ∂G (˜ y, ˜ y ) ∂ ˜ y ∂G (˜ y, ˜ y ) ∂ ˜ y × a s a (˜ y ) a (˜ y ) ˜Π ( kL f , kη , kη ) k . (5.14)From the explicit form of the Green’s functions derived earlier, we can see P h (cid:48) has the correct behav-ior ∝ /a (˜ y ) for the mode deep inside the horizon . The dimensionless source correlator can beobtained from Eq. 2.22, 2.24, 4.49: ˜Π ( kR ∗ c , β c | η − η | ) = π U f (cid:90) d ˜ q P v (˜ q ) P v (˜¯ q ) (1 − µ ) ˜ q ˜¯ q × cos (cid:20) c s ˜ q β c ( η − η ) β c R ∗ c (cid:21) cos (cid:20) c s ˜¯ q β c ( η − η ) β c R ∗ c (cid:21) . (5.15)Here ˜ q = qR ∗ c , a dimensionless quantity, and we use L f = R ∗ c . In Fig. 13, we show this auto-correlator of the source as a function of β c | η − η | . We can see the correlation is quickly lost as β c | η − η | becomes larger than O (1) . Since the source correlator depends only on η − η , wecan change the integration variables from ˜ y , to a quantity proportional to ( η − η ) and anotherindependent linear combination. For RD and MD, the relation between ( η − η ) and y , is givenfrom Eq. 5.5, 5.8: β c ( η − η ) β c R ∗ c = 1 R ∗ c a s H s (cid:26) y − y √ y − √ y ) , (5.16)where the upper row applies to RD and lower one to MD. Then for RD, we can make the followingchange of variables: (cid:26) y y ⇒ (cid:26) y − y ≡ y − , y + y ≡ y + . (5.17) Alternatively, one can express above Green’s functions using the conformal time. The corresponding Green’s functionsare defined to be zero for η (cid:54) η and for η > η . G (˜ η, ˜ η ) = (cid:40) ˜ η ˜ η sin(˜ η − ˜ η ) , RD ˜ η ˜ η [(˜ η − ˜ η ) cos(˜ η − ˜ η ) + (˜ η ˜ η + 1) sin(˜ η − ˜ η )] . MD (5.11)We note that there is a typo in the Green’s function for the matter dominated universe given in Ref. [80], where instead of (˜ η − ˜ η ) cos(˜ η − ˜ η ) , they have − (˜ η − ˜ η ) cos(˜ η − ˜ η ) . For modes deep inside the horizon, ˜ y (cid:29) and ˜ y (cid:29) . Then both Green’s functions take a universal form a a sin(˜ η − ˜ η ) . This implies that h (cid:48) ∝ /a , P h (cid:48) ∝ /a and P GW ∝ /a , behaving like radiation which is true for masslessgravitons. – 33 – β c | η − η |− . . . . . . . . . ˜ Π ( k R ∗ c = , β c | η − η | ) Figure 13 . Autocorrelation of the source for kR ∗ c = 10 , calculated with the explicit expression in Eq. 5.15. The integration range is − y − (cid:54) y + (cid:54) y + y − when − y (cid:54) y − (cid:54) , and y − (cid:54) y + (cid:54) y − y − when (cid:54) y − (cid:54) y − . Similarly for MD, we can perform the following transformations: (cid:26) y y ⇒ (cid:26) λ − λ ≡ y − , λ + λ ≡ y + , (5.18)where λ i = 2 √ y i and the Jacobian is √ y y . The range of integration is y − (cid:54) y + (cid:54) √ y − y − when (cid:54) y − (cid:54) √ y − and − y − (cid:54) y + (cid:54) √ y + y − when − √ y ) (cid:54) y − (cid:54) .It turns out the relation y − (cid:28) y + generally holds, barring special parameter space. This canbe seen from Eq. 5.16 by noting that β c R ∗ c = (8 π ) / v w ≈ v w < , R ∗ c a s H s ∼ O (10 − ) fromFig. 7, and thus y − ∼ O (10 − ) /v w × β c ( η − η ) . Except for extremely small v w , which giveshighly suppressed gravitational waves, we have y − (cid:28) . On the contrary, y + ∼ O (1) . Then we have y − (cid:28) y + . This means in the integration over y + , we can keep the leading order in y − .Now lets look in more detail at the integrand. For RD and MD, the factor containing Green’sfunction can be written as ∂G (˜ y, ˜ y ) ∂ ˜ y ∂G (˜ y, ˜ y ) ∂ ˜ y = y (cid:104) c R ˜ y + c R − y + · · · ] (cid:105) y (cid:104) c M ˜ y + c M − y + · · · ] (cid:105) ≡ (cid:40) y y (cid:41) × G (˜ y, ˜ y , ˜ y ) . (5.19)Then P GW ( y, kR ∗ c ) = [16 πG (¯˜ (cid:15) + ¯˜ p ) ¯ U f ] π H H s y ( kR ∗ c ) × (cid:90) dy − ˜Π ( kR ∗ c , β c | η − η | ) (cid:34)(cid:90) dy + G (˜ y, ˜ y , ˜ y )˜˜ k (cid:40) y − y − y − / y − / (cid:41)(cid:35) . (5.20)In the square bracket, y , are understood to be functions of y ± (note that ˜ y is defined differently formatter and radiation cases). The reason we associate a factor of ˜˜ k − with G is that G ∝ ˜˜ k to a goodapproximation. For both RD and MD, the integral over y + leads to a result in the following form: (cid:20)(cid:90) dy + · · · (cid:21) = 12 Υ( y ) cos (cid:16) ˜˜ ky − (cid:17) . (5.21)– 34 – .0 1.5 2.0 2.5 3.0 - - - Figure 14 . The integrand of y + integration, with y = 3 . Left is RD and right is MD. The blue is the dominantnon-oscillatory part, the magenta dashed is the oscillatory part( kR ∗ c chosen to be . ) and the dark green isthe total contribution. The profile in a wide range of y is shown in Fig. 15. We can see Υ of RD is slightly larger than MD.For both cases, Υ approaches an asymptotic value: for RD and / for MD, irrespective of howlong the source lasts. This is due to the dilution of the source over time, which makes the contributionfrom later time increasingly suppressed. To have a better understanding of the behavior of Υ( y ) , letssee how they can be obtained in a simpler analytical way.First for RD, neglecting terms suppressed by ( R ∗ c a s H s ) or y − , the dominant contributions tothe integrand of the power spectrum are G RD = 12 y (cid:110) cos (cid:104) ˜˜ ky − (cid:105) + cos (cid:104) k ( y − y + ) (cid:105)(cid:111) + · · · . (5.22)The second term is y − independent and is a highly oscillatory function of y + , which averages to zeroduring the integration over y + (see Fig. 14 for the non-oscillatory and oscillatory contributions). Onthe other hand, the first term, a function of y − , when integrated, gives the dominant contribution: Υ RD = 1 − y . (5.23)For y (cid:29) , it approaches an asymptotic value of . Since this asymptotic value can only be reachedfor a long enough source, a realistic phase transition might not satisfy this. We will come to this pointlater. Similarly for MD, we can perform analogous manipulations and keep only the leading orderand also non-oscillatory term: G MD = 8 y cos (cid:104) ˜˜ ky − (cid:105) + · · · . (5.24)Upon integration, it gives the dominant contribution: Υ MD = 23 (cid:18) − y / (cid:19) . (5.25)– 35 – Figure 15 . The function Υ for radiation domination(blue solid) and matter domination(magenta dashed). For y (cid:29) , it approaches the previously observed asymptotic value of / . Thus barring otherdifferences for RD and MD, the different expansion behaviors lead to a suppression of gravitationalwave spectrum for MD, when compared with RD.With Υ( y ) obtained, the power spectrum as a function of y can be written in the following form P GW ( y, kR ∗ c ) = [16 πG (¯˜ (cid:15) + ¯˜ p ) ¯ U f ] π H H s y ( kR ∗ c ) × (cid:20)(cid:90) dy − cos (cid:16) ˜˜ ky − (cid:17) ˜Π ( kR ∗ c , β c | η − η | ) (cid:21) × Υ( y ) . (5.26)Here note that using Eq. 5.16, we have ˜˜ ky − = k ( η − η ) . The integral over y − is obtained byplugging the explicit expression of ˜Π , which results in a three-fold integral. The integration of y − over the three trigonometric functions result in a δ function, and makes the angle integration of ˜ q in Eq. 5.15 trivial. We are left eventually with a one fold integral over the magnitude of ˜ q , and thespectrum can be put in the following standard form: P GW ( y, kR ∗ c ) = 3Γ ¯ U f H R,s H H s ( a s R ∗ c ) ( kR ∗ c ) π ˜ P gw ( kR ∗ ) × y Υ( y ) , (5.27)where Γ = ¯˜ w/ ¯˜ e ≈ / , H R,s is defined to contain only the radiation energy density at t s : H R,s = H s √ − κ M , and the integral is hidden inside ˜ P gw ( kR ∗ ) : (cid:101) P GW ( kR ∗ ) = 14 πc s kR ∗ (cid:18) − c s c s (cid:19) (cid:90) z + z − dzz ( z − z + ) ( z − z − ) z + + z − − z ¯ P v ( z ) ¯ P v ( z + + z − − z ) . (5.28)– 36 – − kR ∗ c − − − − − ( k R ∗ c ) π ˜ P G W ( k R ∗ c ) k k − Figure 16 . The dimensionless gravitational wave power spectrum computed in the sound shell model. Thecalculation was performed for a weak phase transition with α = 0 . , v w = 0 . , and exponential bubblenucleation. The low and high frequency regimes follow the k and k − power law fits respectively (black solidlines) Here z = qR ∗ c , z ± = kR ∗ c c s (1 ± c s ) and ¯ P v ( z ) = π ¯ U f P v ( z ) z . Using Eq. 4.50, the explicit expressionfor ¯ P v is ¯ P v ( z ) = 164 π v w U f (cid:90) d ˜ T ˜ T ν ( ˜ T ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A (cid:32) z ˜ T (8 π ) / v w (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (5.29)Plugging in the explicit expressions of H and H R,s , we have P GW ( y, kR ∗ c ) = 3Γ ¯ U f ( H s a s R ∗ c ) ( kR ∗ c ) π ˜ P gw ( kR ∗ ) × (cid:40) (1 − κ M ) κ M y +1 − κ M (cid:41) × Υ( y ) . (5.30)For both RD and MD, the shape of the spectra are the same to a good approximation, and are thesame as that derived in the sound shell model and thus the properties of its shape [53] apply here forboth cases. In particular, the peak frequency of the spectrum is located at around kR ∗ c ≈ . Thismean a larger or smaller R ∗ can red or blue shift the spectrum respectively. For example, as shownin Fig. 7, increasing v w reduces R ∗ c and thus blue-shift the spectrum. For MD, it has a larger R ∗ andthus red-shift the spectrum.For RD, we recover the result found in Ref. [49], as long as Υ( y ) = 1 , which is only true for y (cid:29) . The reason only this asymptotic value is obtained in Ref. [49] is due to the over-simplifyingassumptions used (see Appendix B), in which case the second terms in both Eq. 5.23 and 5.25 aremissing. Whether or not the asymptotic values can be reached depends on how long the sourceremains active, and we continue in the next section on this question. As we saw earlier, the presence of an asymptotic value for Υ for large y in both cases is due tothe dilution of the source energy density. This asymptotic value was used in Ref. [49] to reach– 37 – Figure 17 . Time elapsed since t s in unit of Hubble time H − s at t ∗ . the conclusion that for RD the effective lifetime of the source is a Hubble time H − s for RD, i.e., τ sw = 1 /H s , which as we have seen is only true if Υ = 1 for y (cid:29) . The question is, however,whether this asymptotic value can be reached in a realistic time frame. In Fig. 17, we show the timeelapsed since the reference time t s , in unit of the Hubble time H − s . For RD, t − t s /H s = y − , (5.31)and for MD t − t s /H s = 23 ( y / − . (5.32)At about a Hubble time, Υ ≈ . for both RD and MD, which is less than a half of the asymptoticvalue for RD and for MD. We need many Hubble times for Υ to approach the asymptoticvalue. The problem is certain physical processes might prohibit the sound waves from being activefor such a long time, and thus the asymptotic value might never be reached. One such process isthe possible formation of shocks and turbulence. Another is the existence of possible dissipativeprocesses, whose presence damps the sound waves. If either of these processes quenches the soundwaves, the asymptotic value will not be achieved. In this case, the effective lifetime is shorter than theHubble time for RD, and the result obtained with an effective lifetime of a Hubble time overestimatesthe gravitational wave production. The time scale for turbulence is roughly [50, 107] τ sw ∼ L f ¯ U f ∼ R ∗ ¯ U f . (5.33)Therefore τ sw /H s ∼ H s R ∗ ¯ U f . (5.34)– 38 – . . . . . . v w . . . . . . . . α Deflagration Hybrid D e t o n a t i o n Forbidden − − − − U f . . . . . . v w . . . . . . . . α Deflagration Hybrid D e t o n a t i o n Forbidden − − − − U f Figure 18 . ¯ U f on the plane of ( v w , α ) . The left figure is ¯ U f of the fluid around a single bubble. The rightfigure is ¯ U f of the fluid calculated from the velocity power spectrum. As we have seen in Fig. 7, H s R s ∼ − and different expansion histories lead to larger or smallervalues. To delay the appearance of turbulence and thus approach the asymptotic value of Υ thusrequires smaller fluid velocity ¯ U f or larger bubble separation. While H s R ∗ depends on specificexpansion behavior adopted, the value of ¯ U f is more or less universal, and its value is shown in Fig. 18on the plane of ( v w , α ) . We show here two versions of it obtained using two different methods: oneby solving the velocity profile around a single bubble and the other by integrating over the velocitypower spectrum (see Ref. [53] for details). Thus whether or not above ratio becomes large enoughdepends on the details of the phase transition in a given cosmological context. Even in cases where theturbulence is delayed or not present, i.e., for sufficiently strong or weak phase transitions respectively,the damping of the sound waves caused by some weak processes could still shorten the lifetime inthe form of shear viscosity [49]. It seem unlikely for any scenario to be very close to the asymptoticvalue. We will mainly consider the case of RD as it is the most frequently encountered scenario. Denote thetemperature after the gravitational wave production as T e with the scale factor being a e . The amountof redshifting is described by the scale factor ratio a e /a . For radiation in thermal equilibrium and inadiabatic expansion, the relation between a e and a is governed by entropy conservation: g s ( T e ) a e T e = g s ( T γ ) a T γ , (5.35)where g s is the relativistic degrees of freedom for entropy; T γ is the temperature of the CMB photonwith T γ ≈ . K . At the present time, the relativistic species includes photons and decoupledneutrinos, thus g s = 2 + × N eff ( ) / ≈ . for N eff = 3 . . Using these, the ratio of thescale factor can be put into the following form: a e a = 1 . × − (cid:18) g s ( T e )100 (cid:19) / (cid:18) T e GeV (cid:19) (cid:18) Hz H e (cid:19) . (5.36)– 39 –or the peak frequency at kR ∗ = z p where z p ≈ [50], the frequency at t e is f p = z p πR ∗ ( t e ) , (5.37)where R ∗ ( t e ) is evaluated at the end of the gravitational wave production and note all previouslygenerated gravitational waves at higer frequencies at kR ∗ c = z p have all redshifted to the frequencyproduced at t e . Then the corresponding frequency today is f SW = 2 . × − Hz (cid:18) g s ( T e )100 (cid:19) / (cid:18) T e GeV (cid:19) (cid:16) z p (cid:17) (cid:18) H e R ∗ ( t e ) (cid:19) . (5.38)We can express R ∗ by β ( v w ) using Eq. 3.50, so that, H e R ∗ ( t e ) = (8 π ) − / a ( t f ) a ( t e ) 1 v w β ( v w ) H e = (8 π ) − / v w β ( v w ) H e × y . (5.39)Here we neglect the very small difference between t f , the time when all the bubbles have disappearedand t s , and we have shown explicitly the dependence of β on v w . Also note β is evaluated at t f when I ( t f ) = 1 . The factor y − is significant when the lifetime of the source is long. Then the presentpeak frequency becomes f SW = 8 . × − Hz v w (cid:18) g s ( T e )100 (cid:19) / (cid:18) T e GeV (cid:19) (cid:16) z p (cid:17) (cid:20) β ( v w ) /yH e (cid:21) . (5.40)For the energy fraction of gravitational waves, the dilution of gravitational waves leads to the fol-lowing connection: h Ω GW ( t , f ) = h (cid:18) a e a (cid:19) (cid:18) H e H (cid:19) Ω GW ( t e , a f /a e ) , = 1 . × − (cid:18) g s ( T e ) (cid:19) / Ω GW ( t e , a f /a e ) . (5.41)Here h ≈ . , the Hubble parameter today in unit of km /s/ Mpc. Then plugging the explicitexpression for P GW in Eq. 5.30, we have h Ω GW ( f ) = 4 . × − (cid:18) g s ( T e ) (cid:19) / Γ ¯ U f [ H s R ∗ ( t s )] A S SW ( f )Υ( y ) . (5.42)Here we have defined A S SW ( f ) to be ( kR ∗ ) ˜ P gw ( kR ∗ ) / (2 π ) with appropriate redshifting factorsincluded. One can either use the prediction from the sound shell model to determine A S SW ( f ) , oruse result from numerical simulations [50]. We choose the latter as it should give a more accurateresult, in which case A ≈ . and [2] S SW ( f ) = (cid:18) ff SW (cid:19) (cid:20)
74 + 3( f /f SW ) (cid:21) / . (5.43)For the term H s R ∗ ( t s ) , similar to Eq. 5.39, we can write H s R ∗ ( t s ) = (8 π ) / v w H s β ( v w ) . (5.44) We use a notation where k in kR ∗ is physical wavenumber, and k in kR ∗ c is a comoving wavenumber. – 40 – - - - - - - - - - - - T a iji T i a n Q i n B B O L I SA D E C I GO A sy m p t o t i c Figure 19 . The present day gravitational wave energy density spectra for H ∗ ∆ t = 0 . , and for H ∗ ∆ t (cid:29) when it takes the asymptotic form. Here ∆ t = t − t s and is the time elapsed since t s , the time when the sourcebecomes active. In all three cases, v w = 0 . , α = 0 . , T e = 100 GeV and β/ ( yH ∗ ) = 100 . The shaded regionsat the top are experimental sensitive regions for several proposed space-based detectors. Therefore the final spectrum is h Ω GW ( f ) = 8 . × − (cid:18) g s ( T e ) (cid:19) / Γ ¯ U f (cid:20) H s β ( v w ) (cid:21) v w S SW ( f ) × Υ( y ) . (5.45)For a long lifetime of the source, the main changes are the suppression factor Υ( y ) . In Fig. 19, weshow the spectra for several choices of H ∗ ∆ t , with z p = 10 (see caption for more details).For MD, apparently the extra dominant matter content will decay to radiation at some time later,which will inject entropy to the standard radiation sector. This can be studied using two methods. Inthe first method, one can assume a very quick and thus instantaneous decay of the matter, whichthen allows to use energy conservation to get the new heated radiation temperature. In the secondmethod, a more precise account of the matter decay is provided, with the conclusion that there is noheating up of the radiation but one gets a slower cooling of the radiation, as was firstly pointed out inRef. [94]. Therefore one needs to follow more closely the entropy evolution by taking into accountfinite matter decay width, following the procedure of Ref. [94] or a more closely related examplestudied in Ref. [81]. This however introduces extra model dependent varieties and is beyond thescope of this work. Note current simulations only probe relatively weak transitions and this spectrum might not be applicable for strongtransitions. As shown in a recent simulation [102], there is a strong deficit in the gravitational wave production for somestrong transitions. – 41 –
Summary
We studied in detail the cosmological first order phase transition and the calculation of resultingstochastic gravitational waves in an expanding universe, with radiation and matter dominated universeas two representative examples. Firstly we studied the changes to process of bubble formation andcollision, including important observables such as the mean bubble separation and its relation with β . We also derived the unbroken bubble wall area, the bubble conformal lifetime distribution whichare needed for the calculation of the gravitational wave spectrum. We then derived the full set ofdifferential equations as used in numerical simulations in an expanding universe. We found thatsimple rescalings works such that the equations governing the velocity profile around a single bubblemaintains the same form as in Minkowski spacetime and that the velocity profile remains the samewhen appropriate substitution of variables are used. We then generalized the sound shell model to theexpanding universe and derived the velocity power spectrum. This result is used to derive analyticallythe gravitational wave power spectrum from the sound waves, the dominant source. We found thatthe standard formula of the spectrum needs to include an additional suppression factor Υ , which isa function of the lifetime of the source. For radiation domination, the asymptotic value of Υ is when the lifetime of the source is very long, and corresponds to the usually adopted spectrum in theliterature. This asymptotic value however can not be reached as the onset of shocks and turbulencemay disrupt the sound waves and possible dissipative processes may further damp it. Thereforean additional suppression factor needs to be taken into account when using the gravitational wavespectrum from sound waves and we provided simple analytical expression for Υ . HG and KS are supported by DOE Grant desc0009956. TRIUMF receives federal funding via acontribution agreement with the National Research Council of Canada. We thank Mark Hindmarshfor helpful communications. We acknowledge Elizabeth Loggia for her early involvement in thisproject. May our field become more accommodating to mothers.
A The Example Effective Potential
Here we provide details of the example effective potential used in Sec. 3, so that those results can bereproduced more easily. The effective potential was originally used as a high temperature approxi-mation for the standard model, given by V ( φ, T ) = D ( T − T ) φ − ET φ + λ φ . (A.1)Here D > , E > , λ > and λ has a weak dependence on T . The first term has a positivecoefficient when T > T to restore the symmetry. The third, the cubic term, when is sufficientlysmaller, helps create a barrier together with the first term, and creates another minimum. Sincethis example is only used to provide a simple benchmark effective potential to show the effects ofthe expansion of the universe, we will take these parameters to be T independent. It should benoted that an effective potential of this form can characterize features of a wide class of beyond thestandard model scenarios in the high temperature approximation. We will use this effective potentialto calculate bounce solutions and corresponding parameters relevant for the phase transition.Though there are four free parameters for this simple effective potential, a rescaling of boththe coordinates and the scalar fields allows to reduce to only one dynamical parameter [108]. The– 42 –
10 20 30 400.00.51.01.52.02.53.0
0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.10 Figure 20 . Left panel: the bounce solutions for the example effective potential with rescaled fields andcoordinates used in this work for different choices of σ , with the color-map denoting values of σ . Right panel:comparison of the corresponding S ( T ) /T obtained with different packages and the analytical fit provided inRef. [108]. rescaled fields and coordinates are defined as Φ = 2
ET φ/M and X = M x . The Lagrangian thenbecomes L = M E T (cid:20)
12 ( ∂ X Φ) −
12 Φ + 12 Φ − σ Φ (cid:21) , (A.2)where σ ≡ λM / (2 E T ) . The behavior of the effective potential for the rescaled fields duringthe phase transition is solely controlled by σ . When σ < / , a second minimum develops. Thecorresponding temperature is: T = (cid:115) T − E λD . (A.3)When σ = 1 , this minimum is degenerate with the one at the origin, which corresponds to a criticaltemperature of T c = (cid:115) T − E λD . (A.4)Therefore for the rescaled field Φ and coordinate X , there is essentially one parameter σ that de-termines the shape of the potential. Calculating the bounce solution and S for all choices of σ issufficient to cover the full four parameter space. Define the S action for the rescaled fields andcoordinates as ˜ S ( σ ) , then the S ( T ) for the original four parameter theory can be obtained directlyas S ( T ) T = M E T ˜ S ( σ ) . (A.5)The bounce solutions for various choices of σ are shown in the left panel of Fig. 20 and the cor-responding S ( σ ) shown as red dotted and green dashed lines for solutions solved from Cosmo-Transitions [109] and BubbleProfiler [110] respectively. In this plot, there is also a purple curve,corresponding to the analytical fit in Ref. [108]: ˜ S ( σ ) = 4 × . × (cid:26) σ (cid:20) . − σ + 0 . − σ ) (cid:21)(cid:27) . (A.6) This is of course different from the σ defined in Eq. 4.25 – 43 –e can see in the whole region plotted, the three results agree very well with each other. So ourresults in previous sections can be followed by simply choosing above analytical fit. For the exampleused in Sec. 3, T = 75 GeV, E = D = 0 . and λ = 0 . , which gives T c = 106 . GeV.
B The Previously Derived Effective Lifetime of the Source
Here we revisit the deviation that led to the conclusion that the effective lifetime of the source isone Hubble time in a radiation dominated universe, as was originally obtained in Ref. [49]. We willfollow closely their notations, using the conformal time η as variable instead of y , and using a ∗ ratherthan a s . Also we study both RD and MD, though only RD is studied in Ref. [49].We start with Eq. 2.14 and do the integrals over ˜ η and ˜ η . We can keep only the leadingcontribution by neglecting the highly oscillatory part in the Green’s functions. This means for thetrigonometric function, we keep only the parts with argument ( ˜ η − ˜ η ) and find ∂G (˜ η, ˜ η ) ∂ ˜ η ∂G (˜ η, ˜ η ) ∂ ˜ η = ˜ η ˜ η × (cid:26) ˜ η − (1 + ˜ η − ) cos(˜ η − ˜ η ) , ˜ η − (1 + 3˜ η − + 9˜ η − )[(˜ η − ˜ η ) sin(˜ η − ˜ η ) + (1 + ˜ η ˜ η ) cos(˜ η − ˜ η )] , (B.1)where the upper and lower row applies to radiation and matter dominated universe respectively. Nowswitch integration variables from ˜ η and ˜ η to x ≡ (˜ η + ˜ η ) / and z = ˜ η − ˜ η . This results in therelation ˜ η ˜ η = x − z . Under these manipulations, the power spectral density of h (cid:48) becomes: P h (cid:48) = [16 πG (¯˜ (cid:15) + ¯˜ p ) ¯ U f ] L f (cid:26) ˜ η − (1 + ˜ η − )˜ η − (1 + 3˜ η − + 9˜ η − ) (cid:27) (cid:90) dx (cid:90) dz k ˜ η ˜ η a ∗ a ( η ) a ( η ) × (cid:26) cos zz sin z + (1 + x − z ) cos z (cid:27) ˜Π ( ˜ L f , ˜ η , ˜ η ) . (B.2)Here ˜ L f ≡ kL f . The expression can be reorganized to show the correct dependence on a ( η ) and wehave for the correlator of ˙ h : P ˙ h = a ∗ a ( η ) 1 k [16 πG (¯˜ (cid:15) + ¯˜ p ) ¯ U f ] L f (cid:26) η − η − + 9˜ η − (cid:27) (cid:90) ˜ η ˜ η ∗ dx (cid:90) dz × (cid:40) ˜ η ∗ x − z / η ∗ ( x − z / (cid:41) (cid:26) cos zz sin z + (1 + x − z ) cos z (cid:27) ˜Π ( ˜ L f , ˜ η , ˜ η ) . (B.3)As we have seen the source is largely stationary, that is, the correlator ˜Π ( ˜ L f , ˜ η , ˜ η ) depends onlyon z but not on x . Then it can be written as ˜Π ( ˜ L f , z ) . Also the autocorrelation time z is very smallcompared with the Hubble time, so we can neglect the z dependence on the denominators in the firstcurly bracket and keep only the x term for MD in the second curly bracket, which then allows theintegration over x , giving (cid:90) ˜ η ˜ η ∗ dx x = 1˜ η ∗ − η , (cid:90) ˜ η ˜ η ∗ dx x = 13 ( 1˜ η ∗ − η ) . (B.4)Here is where things become subtle. The second term for RD is neglected in Ref. [49]. This leads to aresult that corresponds to the asymptotic value Υ = 1 for RD, and as we have seen the short duration– 44 –f the source does not allow to neglect this term. Lets continue to reproduce the result of Ref. [49] bykeeping only the first term. This gives P ˙ h = a ∗ a ( η ) 1 k [16 πG (¯˜ (cid:15) + ¯˜ p ) ¯ U f ] L f (cid:26) η − (1 + 3˜ η − + 9˜ η − ) / (cid:27) ˜ η ∗ × (cid:90) dz cos( z )2 ˜Π ( ˜ L f , z )= a ∗ a ( η ) [16 πG (¯˜ (cid:15) + ¯˜ p ) ¯ U f ] L f (cid:26) η − (1 + 3˜ η − + 9˜ η − ) / (cid:27) ( a ∗ η ∗ )( a ∗ L f ) (cid:101) P GW ( kL f ) . (B.5)In the second line, the following definition is used: (cid:101) P GW ( kL f ) = 1 kL f (cid:90) dz cos z ( ˜ L f , z ) . (B.6)The variables appearing in above equations can further be reorganized so that we have a result similarto Eq.(A10) in Ref. [49]: P GW ( t, k ) = 3Γ ¯ U f (cid:18) a ∗ a H ∗ R H H ∗ (cid:19) (cid:26) η − (1 + 3˜ η − + 9˜ η − ) / (cid:27) × ( H ∗ a ∗ η ∗ )( H ∗ a ∗ L f ) ( kL f ) π (cid:101) P GW ( kL f ) . (B.7)For RD, H ∗ a ∗ η ∗ = 1 and a ∗ L f is the physical length scale ( L ∗ f in Ref. [49]). If we also neglect thevariation of the Hubble rate from H ∗ to H , and since in this case H ∗ R = H ∗ , and also neglect theterms suppressed by / ˜ η in the curly bracket due to the assumed relation ˜ η (cid:29) ˜ η ∗ , then the result forRD reduces to Eq.(A11) in Ref. [49]. Because H ∗ a ∗ η ∗ = 1 and also because the power spectrum inMinkowski spacetime is proportional to H ∗ τ sw , it is concluded in Ref. [49] that the effective lifetimeis a Hubble time. This is true if indeed ˜ η (cid:29) ˜ η ∗ , but as we have seen it requires many Hubble timesfor the asymptotic value to be reached. The sound wave, however, is likely to be disrupted by theonset of shocks or turbulence or damped by other dissipative processes, which certainly do not allowthe sound wave to remain active that long for the asymptotic value to be reached. So the main pointis we can not assume ˜ η (cid:29) ˜ η ∗ and neglect the second term in the first equation of Eq. B.4.While non-relevant here for MD, we can still compare its asymptotic value with what we alreadyfind in previous sections. From above equation we can see the quantity in the curly bracket is / forMD and for RD. But for MD, H ∗ a ∗ η ∗ = 2 , then the asymptotic value of Υ is / for MD, whichis consistent with our previous result. References [1] C. Caprini et al. , “Science with the space-based interferometer eLISA. II: Gravitational waves fromcosmological phase transitions,”
JCAP no. 04, (2016) 001, arXiv:1512.06239[astro-ph.CO] .[2] D. J. Weir, “Gravitational waves from a first order electroweak phase transition: a brief review,”
Phil.Trans. Roy. Soc. Lond.
A376 no. 2114, (2018) 20170126, arXiv:1705.01783 [hep-ph] .[3] A. Mazumdar and G. White, “Cosmic phase transitions: their applications and experimentalsignatures,” arXiv:1811.01948 [hep-ph] . – 45 –
4] G. Bertone et al. , “Gravitational wave probes of dark matter: challenges and opportunities,” arXiv:1907.10610 [astro-ph.CO] .[5] C. Caprini et al. , “Detecting gravitational waves from cosmological phase transitions with LISA: anupdate,”
JCAP (2020) 024, arXiv:1910.13125 [astro-ph.CO] .[6] E. Barausse et al. , “Prospects for Fundamental Physics with LISA,” arXiv:2001.09793[gr-qc] .[7] LIGO Scientific, Virgo
Collaboration, B. P. Abbott et al. , “Upper Limits on the StochasticGravitational-Wave Background from Advanced LIGO’s First Observing Run,”
Phys. Rev. Lett. no. 12, (2017) 121101, arXiv:1612.02029 [gr-qc] . [Erratum: Phys.Rev.Lett. 119, 029901(2017)].[8]
LIGO Scientific, Virgo
Collaboration, B. Abbott et al. , “Search for the isotropic stochasticbackground using data from Advanced LIGO’s second observing run,”
Phys. Rev. D no. 6, (2019)061101, arXiv:1903.02886 [gr-qc] .[9] P. Amaro-Seoan et al. , “Laser interferometer space antenna,” arxiv:1702.00786 , arXiv:1702.00786 [astro-ph.IM] .[10] K. Yagi and N. Seto, “Detector configuration of DECIGO/BBO and identification of cosmologicalneutron-star binaries,” Phys. Rev.
D83 (2011) 044011, arXiv:1101.3940 [astro-ph.CO] .[Erratum: Phys. Rev.D95,no.10,109901(2017)].[11] X. Gong et al. , “Descope of the ALIA mission,”
J. Phys. Conf. Ser. no. 1, (2015) 012011, arXiv:1410.7296 [gr-qc] .[12]
TianQin
Collaboration, J. Luo et al. , “TianQin: a space-borne gravitational wave detector,”
Class.Quant. Grav. no. 3, (2016) 035010, arXiv:1512.02076 [astro-ph.IM] .[13] C. Grojean and G. Servant, “Gravitational Waves from Phase Transitions at the Electroweak Scale andBeyond,” Phys. Rev.
D75 (2007) 043507, arXiv:hep-ph/0607107 [hep-ph] .[14] V. Vaskonen, “Electroweak baryogenesis and gravitational waves from a real scalar singlet,”
Phys. Rev.
D95 no. 12, (2017) 123515, arXiv:1611.02073 [hep-ph] .[15] G. Dorsch, S. Huber, T. Konstandin, and J. No, “A Second Higgs Doublet in the Early Universe:Baryogenesis and Gravitational Waves,”
JCAP (2017) 052, arXiv:1611.05874 [hep-ph] .[16] A. Beniwal, M. Lewicki, M. White, and A. G. Williams, “Gravitational waves and electroweakbaryogenesis in a global study of the extended scalar singlet model,” arXiv:1810.02380[hep-ph] .[17] Z. Kang, P. Ko, and T. Matsui, “Strong First Order EWPT and Strong Gravitational Waves in Z -symmetric Singlet Scalar Extension,” arXiv:1706.09721 [hep-ph] .[18] C. Delaunay, C. Grojean, and J. D. Wells, “Dynamics of Non-renormalizable Electroweak SymmetryBreaking,” JHEP (2008) 029, arXiv:0711.2511 [hep-ph] .[19] M. Chala, C. Krause, and G. Nardini, “Signals of the electroweak phase transition at colliders andgravitational wave observatories,” JHEP (2018) 062, arXiv:1802.02168 [hep-ph] .[20] S. A. Ellis, S. Ipek, and G. White, “Electroweak Baryogenesis from Temperature-Varying Couplings,” JHEP (2019) 002, arXiv:1905.11994 [hep-ph] .[21] A. Alves, D. Gonalves, T. Ghosh, H.-K. Guo, and K. Sinha, “Di-Higgs Production in the b Channeland Gravitational Wave Complementarity,”
JHEP (2020) 053, arXiv:1909.05268[hep-ph] .[22] A. Alves, T. Ghosh, H.-K. Guo, K. Sinha, and D. Vagie, “Collider and Gravitational WaveComplementarity in Exploring the Singlet Extension of the Standard Model,” JHEP (2019) 052, arXiv:1812.09333 [hep-ph] . – 46 –
23] J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen, “Gravitational wave energy budget in stronglysupercooled phase transitions,”
JCAP (2019) 024, arXiv:1903.09642 [hep-ph] .[24] A. P. Morais and R. Pasechnik, “Probing multi-step electroweak phase transition with multi-peakedprimordial gravitational waves spectra,” JCAP (2020) 036, arXiv:1910.00717 [hep-ph] .[25] A. Addazi, A. Marcian, and R. Pasechnik, “Probing Trans-electroweak First Order Phase Transitionsfrom Gravitational Waves,” MDPI Physics no. 1, (2019) 92–102, arXiv:1811.09074[hep-ph] .[26] C. Wainwright, S. Profumo, and M. J. Ramsey-Musolf, “Gravity Waves from a Cosmological PhaseTransition: Gauge Artifacts and Daisy Resummations,” Phys. Rev. D (2011) 023521, arXiv:1104.5487 [hep-ph] .[27] R. Zhou, L. Bian, and H.-K. Guo, “Probing the Electroweak Sphaleron with Gravitational Waves,” Phys. Rev. D no. 9, (2020) 091903, arXiv:1910.00234 [hep-ph] .[28] J. Bernon, L. Bian, and Y. Jiang, “A new insight into the phase transition in the early Universe withtwo Higgs doublets,”
JHEP (2018) 151, arXiv:1712.08430 [hep-ph] .[29] P. Schwaller, “Gravitational Waves from a Dark Phase Transition,” Phys. Rev. Lett. no. 18, (2015)181101, arXiv:1504.07263 [hep-ph] .[30] J. Jaeckel, V. V. Khoze, and M. Spannowsky, “Hearing the signal of dark sectors with gravitationalwave detectors,”
Phys. Rev. D no. 10, (2016) 103519, arXiv:1602.03901 [hep-ph] .[31] M. Chala, G. Nardini, and I. Sobolev, “Unified explanation for dark matter and electroweakbaryogenesis with direct detection and gravitational wave signatures,” Phys. Rev. D no. 5, (2016)055006, arXiv:1605.08663 [hep-ph] .[32] A. Addazi, “Limiting First Order Phase Transitions in Dark Gauge Sectors from Gravitational Wavesexperiments,” Mod. Phys. Lett. A no. 08, (2017) 1750049, arXiv:1607.08057 [hep-ph] .[33] W. Chao, H.-K. Guo, and J. Shu, “Gravitational Wave Signals of Electroweak Phase TransitionTriggered by Dark Matter,” JCAP no. 09, (2017) 009, arXiv:1702.02698 [hep-ph] .[34] I. Baldes, “Gravitational waves from the asymmetric-dark-matter generating phase transition,”
JCAP (2017) 028, arXiv:1702.02117 [hep-ph] .[35] A. Addazi and A. Marciano, “Gravitational waves from dark first order phase transitions and darkphotons,” Chin. Phys.
C42 no. 2, (2018) 023107, arXiv:1703.03248 [hep-ph] .[36] D. Croon, V. Sanz, and G. White, “Model Discrimination in Gravitational Wave spectra from DarkPhase Transitions,”
JHEP (2018) 203, arXiv:1806.02332 [hep-ph] .[37] I. Baldes and C. Garcia-Cely, “Strong gravitational radiation from a simple dark matter model,” JHEP (2019) 190, arXiv:1809.01198 [hep-ph] .[38] M. Fairbairn, E. Hardy, and A. Wickens, “Hearing without seeing: gravitational waves from hot andcold hidden sectors,” JHEP (2019) 044, arXiv:1901.11038 [hep-ph] .[39] D. Dunsky, L. J. Hall, and K. Harigaya, “Dark Matter, Dark Radiation and Gravitational Waves fromMirror Higgs Parity,” JHEP (2020) 078, arXiv:1908.02756 [hep-ph] .[40] P. Archer-Smith, D. Linthorne, and D. Stolarski, “Gravitational Wave Signals from Multiple HiddenSectors,” Phys. Rev. D no. 9, (2020) 095016, arXiv:1910.02083 [hep-ph] .[41] E. Hall, T. Konstandin, R. McGehee, and H. Murayama, “Asymmetric Matters from a DarkFirst-Order Phase Transition,” arXiv:1911.12342 [hep-ph] .[42] L. Bian, W. Cheng, H.-K. Guo, and Y. Zhang, “Gravitational waves triggered by B − L chargedhidden scalar and leptogenesis,” arXiv:1907.13589 [hep-ph] . – 47 –
43] X. Wang, F. P. Huang, and X. Zhang, “Phase transition dynamics and gravitational wave spectra ofstrong first-order phase transition in supercooled universe,”
JCAP (2020) 045, arXiv:2003.08892 [hep-ph] .[44] D. Croon, T. E. Gonzalo, and G. White, “Gravitational Waves from a Pati-Salam Phase Transition,” arXiv:1812.02747 [hep-ph] .[45] A. Greljo, T. Opferkuch, and B. A. Stefanek, “Gravitational Imprints of Flavor Hierarchies,” Phys.Rev. Lett. no. 17, (2020) 171802, arXiv:1910.02014 [hep-ph] .[46] W.-C. Huang, F. Sannino, and Z.-W. Wang, “Gravitational Waves from Pati-Salam Dynamics,” arXiv:2004.02332 [hep-ph] .[47] V. Brdar, L. Graf, A. J. Helmboldt, and X.-J. Xu, “Gravitational Waves as a Probe of Left-RightSymmetry Breaking,”
JCAP (2019) 027, arXiv:1909.02018 [hep-ph] .[48] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, “Gravitational waves from the sound ofa first order phase transition,” Phys. Rev. Lett. (2014) 041301, arXiv:1304.2433 [hep-ph] .[49] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, “Numerical simulations of acousticallygenerated gravitational waves at a first order phase transition,”
Phys. Rev.
D92 no. 12, (2015) 123009, arXiv:1504.03291 [astro-ph.CO] .[50] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, “Shape of the acoustic gravitationalwave power spectrum from a first order phase transition,”
Phys. Rev.
D96 no. 10, (2017) 103520, arXiv:1704.05871 [astro-ph.CO] .[51] J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, “Energy Budget of Cosmological First-orderPhase Transitions,”
JCAP (2010) 028, arXiv:1004.4187 [hep-ph] .[52] M. Hindmarsh, “Sound shell model for acoustic gravitational wave production at a first-order phasetransition in the early Universe,”
Phys. Rev. Lett. no. 7, (2018) 071301, arXiv:1608.04735[astro-ph.CO] .[53] M. Hindmarsh and M. Hijazi, “Gravitational waves from first order cosmological phase transitions inthe Sound Shell Model,”
JCAP no. 12, (2019) 062, arXiv:1909.10040 [astro-ph.CO] .[54] G. Kane, K. Sinha, and S. Watson, “Cosmological Moduli and the Post-Inflationary Universe: ACritical Review,”
Int. J. Mod. Phys. D no. 08, (2015) 1530022, arXiv:1502.07746[hep-th] .[55] R. Allahverdi et al. , “The First Three Seconds: a Review of Possible Expansion Histories of the EarlyUniverse,” arXiv:2006.16182 [astro-ph.CO] .[56] T. Banks, M. Berkooz, S. Shenker, G. W. Moore, and P. Steinhardt, “Modular cosmology,” Phys. Rev.D (1995) 3548–3562, arXiv:hep-th/9503114 .[57] T. Banks, M. Berkooz, and P. Steinhardt, “The Cosmological moduli problem, supersymmetrybreaking, and stability in postinflationary cosmology,” Phys. Rev. D (1995) 705–716, arXiv:hep-th/9501053 .[58] G. Coughlan, W. Fischler, E. W. Kolb, S. Raby, and G. G. Ross, “Cosmological Problems for thePolonyi Potential,” Phys. Lett. B (1983) 59–64.[59] T. Banks, D. B. Kaplan, and A. E. Nelson, “Cosmological implications of dynamical supersymmetrybreaking,”
Phys. Rev. D (1994) 779–787, arXiv:hep-ph/9308292 .[60] B. Dutta, L. Leblond, and K. Sinha, “Mirage in the Sky: Non-thermal Dark Matter, Gravitino Problem,and Cosmic Ray Anomalies,” Phys. Rev. D (2009) 035014, arXiv:0904.3773 [hep-ph] .[61] R. Allahverdi, B. Dutta, and K. Sinha, “Successful Supersymmetric Dark Matter with ThermalOver/Under-Abundance from Late Decay of a Visible Sector Scalar,” Phys. Rev. D (2013) 075024, arXiv:1212.6948 [hep-ph] . – 48 –
62] B. S. Acharya, P. Kumar, K. Bobkov, G. Kane, J. Shao, and S. Watson, “Non-thermal Dark Matter andthe Moduli Problem in String Frameworks,”
JHEP (2008) 064, arXiv:0804.0863[hep-ph] .[63] B. S. Acharya, G. Kane, S. Watson, and P. Kumar, “A Non-thermal WIMP Miracle,” Phys. Rev. D (2009) 083529, arXiv:0908.2430 [astro-ph.CO] .[64] A. L. Erickcek, K. Sinha, and S. Watson, “Bringing Isolated Dark Matter Out of Isolation: Late-timeReheating and Indirect Detection,” Phys. Rev. D no. 6, (2016) 063502, arXiv:1510.04291[hep-ph] .[65] M. Sten Delos, T. Linden, and A. L. Erickcek, “Breaking a dark degeneracy: The gamma-ray signatureof early matter domination,” Phys. Rev. D no. 12, (2019) 123546, arXiv:1910.08553[astro-ph.CO] .[66] A. L. Erickcek, “The Dark Matter Annihilation Boost from Low-Temperature Reheating,”
Phys. Rev.D no. 10, (2015) 103505, arXiv:1504.03335 [astro-ph.CO] .[67] M. Drees and F. Hajkarim, “Dark Matter Production in an Early Matter Dominated Era,” JCAP (2018) 057, arXiv:1711.05007 [hep-ph] .[68] C. Cosme, M. Dutra, T. Ma, Y. Wu, and L. Yang, “Neutrino Portal to FIMP Dark Matter with an EarlyMatter Era,” arXiv:2003.01723 [hep-ph] .[69] R. Allahverdi, B. Dutta, and K. Sinha, “Baryogenesis and Late-Decaying Moduli,” Phys. Rev. D (2010) 035004, arXiv:1005.2804 [hep-ph] .[70] C. Pallis, “Kination-dominated reheating and cold dark matter abundance,” Nucl. Phys. B (2006)129–159, arXiv:hep-ph/0510234 .[71] J. Lankinen, O. Kerppo, and I. Vilja, “Reheating via gravitational particle production in the kinationepoch,”
Phys. Rev. D no. 6, (2020) 063529, arXiv:1910.07520 [gr-qc] .[72] K. Nakayama and F. Takahashi, “Running Kinetic Inflation,”
JCAP (2010) 009, arXiv:1008.2956 [hep-ph] .[73] C. Pallis, “Quintessential kination and cold dark matter abundance,” JCAP (2005) 015, arXiv:hep-ph/0503080 .[74] D. Grin, T. L. Smith, and M. Kamionkowski, “Axion constraints in non-standard thermal histories,” Phys. Rev. D (2008) 085020, arXiv:0711.1352 [astro-ph] .[75] K. Dimopoulos and T. Markkanen, “Non-minimal gravitational reheating during kination,” JCAP (2018) 021, arXiv:1803.07399 [gr-qc] .[76] K. Redmond and A. L. Erickcek, “New Constraints on Dark Matter Production during Kination,” Phys. Rev. D no. 4, (2017) 043511, arXiv:1704.01056 [hep-ph] .[77] D. Bettoni, G. Domnech, and J. Rubio, “Gravitational waves from global cosmic strings inquintessential inflation,” JCAP (2019) 034, arXiv:1810.11117 [astro-ph.CO] .[78] D. Bettoni and J. Rubio, “Hubble-induced phase transitions: Walls are not forever,” JCAP (2020)002, arXiv:1911.03484 [astro-ph.CO] .[79] S. Bhattacharya, S. Mohanty, and P. Parashari, “Primordial black holes and gravitational waves innon-standard cosmologies,” arXiv:1912.01653 [astro-ph.CO] .[80] G. Barenboim and W.-I. Park, “Gravitational waves from first order phase transitions as a probe of anearly matter domination era and its inverse problem,” Phys. Lett.
B759 (2016) 430–438, arXiv:1605.03781 [astro-ph.CO] .[81] N. Bernal and F. Hajkarim, “Primordial Gravitational Waves in Nonstandard Cosmologies,”
Phys. Rev.
D100 no. 6, (2019) 063502, arXiv:1905.10410 [astro-ph.CO] . – 49 –
82] C. Caprini and D. G. Figueroa, “Cosmological Backgrounds of Gravitational Waves,”
Class. Quant.Grav. no. 16, (2018) 163001, arXiv:1801.04268 [astro-ph.CO] .[83] F. D’Eramo and K. Schmitz, “Imprint of a scalar era on the primordial spectrum of gravitationalwaves,” Phys. Rev. Research. (2019) 013010, arXiv:1904.07870 [hep-ph] .[84] M. Geller, A. Hook, R. Sundrum, and Y. Tsai, “Primordial Anisotropies in the Gravitational WaveBackground from Cosmological Phase Transitions,” Phys. Rev. Lett. no. 20, (2018) 201303, arXiv:1803.10780 [hep-ph] .[85] Y. Cui, M. Lewicki, D. E. Morrissey, and J. D. Wells, “Cosmic Archaeology with Gravitational Wavesfrom Cosmic Strings,”
Phys. Rev. D no. 12, (2018) 123505, arXiv:1711.03104 [hep-ph] .[86] J. Ellis, M. Lewicki, and J. M. No, “On the Maximal Strength of a First-Order Electroweak PhaseTransition and its Gravitational Wave Signal,” Submitted to: JCAP (2018) , arXiv:1809.08242[hep-ph] .[87] J. Ellis, M. Lewicki, and J. M. No, “Gravitational waves from first-order cosmological phasetransitions: lifetime of the sound wave source,” arXiv:2003.07360 [hep-ph] .[88] S. Weinberg,
Cosmology . 9, 2008.[89] A. D. Linde, “Fate of the False Vacuum at Finite Temperature: Theory and Applications,”
Phys. Lett. (1981) 37–40.[90] A. D. Linde, “Decay of the False Vacuum at Finite Temperature,”
Nucl. Phys.
B216 (1983) 421.[Erratum: Nucl. Phys.B223,544(1983)].[91] G. V. Dunne and H. Min, “Beyond the thin-wall approximation: Precise numerical computation ofprefactors in false vacuum decay,”
Phys. Rev.
D72 (2005) 125004, arXiv:hep-th/0511156[hep-th] .[92] A. Andreassen, D. Farhi, W. Frost, and M. D. Schwartz, “Precision decay rate calculations in quantumfield theory,”
Phys. Rev.
D95 no. 8, (2017) 085011, arXiv:1604.06090 [hep-th] .[93] S. Weinberg,
The quantum theory of fields. Vol. 2: Modern applications . Cambridge University Press,2013.[94] R. J. Scherrer and M. S. Turner, “Decaying Particles Do Not Heat Up the Universe,”
Phys. Rev.
D31 (1985) 681.[95] A. H. Guth and E. J. Weinberg, “Cosmological Consequences of a First Order Phase Transition in theSU(5) Grand Unified Model,”
Phys. Rev.
D23 (1981) 876.[96] G. D. Moore and T. Prokopec, “How fast can the wall move? A Study of the electroweak phasetransition dynamics,”
Phys. Rev.
D52 (1995) 7182–7204, arXiv:hep-ph/9506475 [hep-ph] .[97] R. Apreda, M. Maggiore, A. Nicolis, and A. Riotto, “Gravitational waves from electroweak phasetransitions,”
Nucl. Phys.
B631 (2002) 342–368, arXiv:gr-qc/0107033 [gr-qc] .[98] R.-G. Cai and S.-J. Wang, “Energy budget of cosmological first-order phase transition in FLRWbackground,”
Sci. China Phys. Mech. Astron. (2018) 080411, arXiv:1803.03002 [gr-qc] .[99] J. Ignatius, K. Kajantie, H. Kurki-Suonio, and M. Laine, “The growth of bubbles in cosmologicalphase transitions,” Phys. Rev. D (1994) 3854–3868, arXiv:astro-ph/9309059 .[100] H. Kurki-Suonio and M. Laine, “On bubble growth and droplet decay in cosmological phasetransitions,” Phys. Rev. D (1996) 7163–7171, arXiv:hep-ph/9512202 .[101] B.-H. Liu, L. D. McLerran, and N. Turok, “Bubble nucleation and growth at a baryon numberproducing electroweak phase transition,” Phys. Rev.
D46 (1992) 2668–2688.[102] D. Cutting, M. Hindmarsh, and D. J. Weir, “Vorticity, kinetic energy, and suppressed gravitationalwave production in strong first order phase transitions,” arXiv:1906.00480 [hep-ph] . – 50 – JCAP no. 09, (2014) 028, arXiv:1407.3132 [hep-ph] .[104] A. Brandenburg, K. Enqvist, and P. Olesen, “Large scale magnetic fields from hydromagneticturbulence in the very early universe,”
Phys. Rev.
D54 (1996) 1291–1300, arXiv:astro-ph/9602031 [astro-ph] .[105] R. M. Gailis, N. E. Frankel, and C. P. Dettmann, “Magnetohydrodynamics in the expanding Universe,”
Phys. Rev.
D52 no. 12, (1995) 6901.[106] H. Kurki-Suonio and M. Laine, “Supersonic deflagrations in cosmological phase transitions,”
Phys.Rev.
D51 (1995) 5431–5437, arXiv:hep-ph/9501216 [hep-ph] .[107] U.-L. Pen and N. Turok, “Shocks in the Early Universe,”
Phys. Rev. Lett. no. 13, (2016) 131301, arXiv:1510.02985 [astro-ph.CO] .[108] M. Dine, R. G. Leigh, P. Y. Huet, A. D. Linde, and D. A. Linde, “Towards the theory of theelectroweak phase transition,”
Phys. Rev. D (1992) 550–571, arXiv:hep-ph/9203203 .[109] C. L. Wainwright, “CosmoTransitions: Computing Cosmological Phase Transition Temperatures andBubble Profiles with Multiple Fields,” Comput. Phys. Commun. (2012) 2006–2013, arXiv:1109.4189 [hep-ph] .[110] P. Athron, C. Balzs, M. Bardsley, A. Fowlie, D. Harries, and G. White, “BubbleProfiler: finding thefield profile and action for cosmological phase transitions,” arXiv:1901.03714 [hep-ph] ..