Gravitational waves from bubble collisions in FLRW spacetime

They adopt thin wall approximation and envelope approximation in their work. The so-called thin wall approximation means that we assume all energy of the bubble is located in the bubble wall with an infinitesimal width $\ell_{B}$. With the help of this approximation, they write the energy momentum tensor $T^{B}$ of the uncollided wall of a single bubble nucleated at $x_{N}\equiv(t_{N},\bm{x}_{N})$ as:

where

and

At here, $x\equiv(t,\bm{x})$ denotes the spacetime point, $\hat{\bullet}$ indicates the unit vector in the same direction of $\vec{\bullet}$, $v$ is the speed of the bubble wall, $\rho_{0}$ represents the energy density released by the phase transitions, and $\kappa$ denotes the fraction of the energy localized in the bubble wall to the released energy during the phase transition.

Envelope approximation means that as long as two or more bubbles collide, the colliding parts of bubble walls and corresponding energy-momentum tensor vanish at once. With the guarantee of envelope approximation, every spatial point($\bm{x}$) can be passed by bubble only once. Since once a spatial point enter the inner part of any single bubble, it will transform from the false vacuum state to the true vacuum state and the inverse process is forbidden. See FIG.1 for a straightforward understanding of this important approximation.

Besides these two approximations, they denote the bubble nucleation rate per unit time per volume as:

where $t_{\star}$ denotes a time point typically around the transition time, $\Gamma_{\star}$ is the transition rate at $t_{\star}$, and $\beta$ is a constant.

In order to obtain the GW spectrum, they define equal-time correlator of GW and unequal-time correlator of energy momentum tensor given by Eq.(6), which we’ll also take use of in our own derivation.

Where $\Pi_{ij}$ is the transverse-traceless part of the energy-momentum tensor $T_{ij}$.

According to the assumptions and approximations given by the previous subsection, they find the core quantity we need to obtain for describing the GW spectrum is the transverse and traceless part of $\langle T_{ij}(x)T_{ij}(y)\rangle$ i.e. $\langle\Pi_{ij}(x)\Pi_{ij}(y)\rangle$. After figuring it out, we can get $\Pi(t_{x},t_{y},k)$ with the help of Eq.(6) and Fourier transformation. As long as we have $\Pi(t_{x},t_{y},k)$, we can calculate the dimensionless GW energy density per logarithmic frequency $\Omega_{\text{GW}}$ by:

The superscript M indicates Minkowski spacetime. The definition of $\Delta(k/\beta,v)$ is given by

Here $\rho_{\text{tot}}=\rho_{0}+\rho_{\text{rad}}$, $\alpha:=\rho_{0}/\rho_{\text{rad}}$ represents the fraction of the released energy density from the phase transitions to the energy density of radiation, $H_{\star}$ denotes the typical value of Hubble parameter around the transition time, $t_{i}$ and $t_{f}$ denote the start time and end time of phase transitions. $\sigma$ is defined as the ratio of $\beta$ and $H_{\star}$, which can characterize the speed of phase transitions to the expansion speed of the universe. For the convenience of latter comparison, we define a new variable $\Delta^{\text{M}}$ at here, we’ll see that $\Delta^{\text{M}}$ is the specific variable that we’ll use to compare with our results in Sec.5 rather than $\Delta$ itself. To prevent possible confusion, we’ll explicitly write superscript M to distinguish $\Delta$ and $\Delta^{\text{M}}$.

According to analyses, $\Delta$ can be divided into two parts as $\Delta=\Delta^{(s)}+\Delta^{(d)}$. $\Delta^{(s)}$ and $\Delta^{(d)}$ correspond to the contribution to $\Delta$ from the single bubble case and double bubble case respectively. The so-called single bubble case and double bubble case need to be understood by considering the relative positions of bubble walls and two spatial points $\bm{x}$ and $\bm{y}$. It is easy to realize that only if $T_{ij}(x)\neq 0$ and $T_{ij}(y)\neq 0$ can this situation contributes to $\langle T_{ij}(x)T_{ij}(y)\rangle$, which means that the spatial points $\bm{x}$ and $\bm{y}$ must be in the bubble wall(s) at the time points $t_{x}$ and $t_{y}$, respectively. So we can find that there are two possible situations. The first one is that $\bm{x}$ and $\bm{y}$ are passed by the same bubble at $t_{x}$ and $t_{y}$ respectively, and the second one is that $\bm{x}$ and $\bm{y}$ are passed by two different bubbles at $t_{x}$ and $t_{y}$ respectively. FIG.2 displays these two possible situations, the left panel and right panel correspond to the single bubble case and double bubble case, respectively. Because $t_{x}$ can be equal or not equal to $t_{y}$, for each case there are two more possibilities. The upper left panel shows the situation where $t_{x}=t_{y}$. The spatial points $\bm{x}$ and $\bm{y}$ are passed by the same bubble at $t=t_{x}=t_{y}$. The lower left panel shows the situation where the bubble firstly passes the spatial point $\bm{x}$ at $t=t_{x}$ and lately passes the spatial point $\bm{y}$ at $t=t_{y}$. The right panel can be understood similarly.

In order to calculate $\Delta^{(s)}$ and $\Delta^{(d)}$, we need to study the past cones of $x$ and $y$ to decide the permitted spacetime points for bubbles to nucleate. For each possible nucleation point, we need to figure out two quantities. The first one is the probability of such a situation, and the second one is the contribution to $\langle T_{ij}(x)T_{ij}(y)\rangle$ from this situation. By summing up all of the possible situations, we can get the final result of $\langle T_{ij}(x)T_{ij}(y)\rangle$. In a word, calculation process can be described by $\displaystyle{\langle\text{TT}\rangle(x,y)=\sum_{\text{situation}}\text{P}(\text{situation})\text{TT}(x,y)\Big{|}_{\text{situation}}}$. The detailed calculation idea can be found in their Appendix A.

As a result of their work, Jinno and Takimoto find that the GW spectrum can be described as two double integrations which can be calculated numerically. The shape of GW spectrum they’ve got is consistent with the result given by numerical simulations. They also give the fitting formula of $\Delta$ as:

with $\Delta_{\text{peak}}=0.043,f_{\text{peak}}/\beta=1.24/2\pi$ and $(c_{l},c_{h})=(0.064,0.48)$. We’ll use this formula to calculate $\Delta$ and compare it with our results in Sec.5.

In this work, we adopt thin wall approximation and envelope approximation as Jinno and Takimoto’s choice in their work. The concrete definitions of these two approximations have been discussed in Sec.2. According to Jinno and Takimoto’s setup, we can write the energy-momentum tensor of an uncollided bubble nucleated at $x_{N}\equiv(\eta_{N},\bm{x}_{N})$ as:

where

with $\Delta\bm{x}\equiv a(\eta)\absolutevalue{\bm{x}-\bm{x}_{N}}$. To remind, the expressions of $T_{ij}$ given by Eq.(13) is a modeling rather than a rigorous expression by derivation. The physical meanings of $\rho_{0}$ and $\kappa$ are both the same as which has been given in the previous section. In our work, we only consider the situation where the speed of bubble walls equals to the speed of the light, i.e. $v=c=1$. So we have:

Note that $\ell_{B}$ is the comoving width of the bubble wall. For the future convenience, we define $\Psi_{ij}(\eta,\bm{x})$ and $\Psi_{ij}(\eta,\bm{k})$ as:

where superscript TT is the abbreviation of Transverse-Traceless. $\Psi_{ij}(\eta,\bm{k})$ and $T_{ij}^{\text{TT}}(\eta,\bm{k})$ are the Fourier modes of $\Psi_{ij}(\eta,\bm{x})$ and $T_{ij}^{\text{TT}}(\eta,\bm{x})$ with wave vector $\bm{k}$, respectively. In this work, our conventions of Fourier transformation are given by $\displaystyle{\int\text{d}^{3}xe^{i\bm{k}\cdot\bm{x}}}$ and $\displaystyle{\int\text{d}^{3}k/(2\pi)^{3}e^{-i\bm{k}\cdot\bm{x}}}$. We can extract the TT parts of $T_{ij}(\eta,\bm{k})$ by contracting it with Lambda tensor whose definition is given by Eq.(22).

The transition ratio per unit conformal time per comoving volume we use in our paper is given byCaprini_2008 :

where $\eta_{\star}$ denotes the conformal time corresponding to the start of phase transitions, and $\Gamma_{\star}$ is the value of $\Gamma(\eta)$ at this specific time. Note that the definition of $\widetilde{\beta}$ is given by $\widetilde{\beta}:=a(\eta)\beta\Big{|}_{\eta=\eta_{\star}}=$const. where $\beta$ is more common to be used in the literature e.g. Ref.Kosowsky_1993 . At here, it is necessary to explain that Eq.(17) only holds when $\eta\sim\eta_{\star}$, so in our work, we consider a finite length of phase transitions to make sure this relationship can be satisfied. Actually, later we’ll see the choice of FLRW spacetime itself gives us a relatively natural way to determine the duration of phase transitions which will directly influence the final GW spectrum. In most literature, people estimate the duration of the phase transitions by $\beta^{-1}$, and in Jinno and Takimoto’s work they set the phase transitions start at $t_{i}=-\infty$ and end at $t_{f}=+\infty$ which itself means that they regard the duration of phase transitions as infinite long in their actual derivation and numerical calculation. Although they give an argument to this point that this treatment doesn’t change their final results, in our view, a finite duration might be a better choice and could make the results more reliable. The exact method to determine the length of phase transitions will be discussed in the latter section.

Without loss of generality, we will set $\widetilde{\beta}\equiv 1$ in the rest of this paper which is just a choice of the unit like setting the speed of light $c\equiv 1$. Besides, it’s also worth mentioning that although $\widetilde{\beta}$ is a constant given a certain value of $\beta$ and $a$ i.e. given a specific physical scenario, we are not allowed to compare any dimensional physical quantity by comparing the number of them in different scenarios. Because the physical meaning of “1” in each particular scenario is different. However, by the virtue of dimensionless property of the quantity $h_{0}^{2}\Omega_{\text{GW}}$, we can compare it safely in different physical scenarios.

In this part we begin to derive the analytical expressions of GW spectrum using the model proposed by Ref.Jinno_2017 .

Let’s consider the equation of motion of the metric perturbation satisfy at first. As we have stated already, our spacetime background is FLRW spacetime. Including tensor perturbations, we have the spacetime metric as:

Here $h_{ij}$ are all traceless and transverse. By calculating the Christoffel symbol, Riemman tensor and Ricci tensor, we have the linearized Einstein’s equation to the first order in $h_{ij}$, over the FLRW spacetime as Caprini_2018 :

Where $H$ is Hubble parameter at $t$, $\nabla^{2}\equiv\partial^{i}\partial_{i}$ indicates the Laplacian associated with comoving coordinates $x^{i}$ and $\dot{\bullet}$ denotes the derivative of $\bullet$ respective to cosmic time $t$. Define $H_{ij}(\eta,\bm{x})=a(\eta)h_{ij}(\eta,\bm{x})$ and Eq.(19) in Fourier space becomes

where $\bullet^{\prime}$ denotes the derivative of $\bullet$ with respect to conformal time $\eta$. For the convenience of subsequent derivation, we define $\zeta\equiv k\eta$ at here. Since we assume that the phase transitions happen in RD era, we can use the relationship of scale factor and conformal time $a(\eta)=\chi\eta$ to simplify the equation of motion and obtain Eq.(21) where $\chi$ is a constant.

As we have stated before, $\Psi_{ij}(\eta,\bm{x})$ is related to the TT part of energy-momentum tensor $T_{ij}$, and we can use Lambda tensor to extract the TT part of $T_{ij}$. The definition of Lambda tensor is given by:

It’s easy to check that the Lambda tensor has a basic property which is of great use:

So we can rewrite $\Psi_{ij}$ in terms of $T_{ij}$ with the help of Lambda tensor:

Assuming that the phase transitions begin at $\zeta_{i}=k\eta_{i}$ and end at $\zeta_{f}=k\eta_{f}$. $\Psi_{ij}=0$ before the phase transitions. Let’s consider a time point $\zeta_{i}<\zeta<\zeta_{f}$ during the phase transitions, we can get the solution of $H_{ij}$ by Green’s function method:

When phase transitions finish, $\Psi_{ij}=0$ (no bubbles anymore, all spacetime points stay in the true vacuum state), so we have

We require $H_{ij}$ to satisfy the connection requirements at $\zeta=\zeta_{f}$, i.e. $H_{ij}|_{\zeta_{f}}$ and $H^{\prime}_{ij}|_{\zeta_{f}}$ given by Eq.(25) and Eq.(26) should be equal to each other. Using these requirements, we can arrive at

Insert Eq.(27) into Eq.(26), we can get the final expression of $H_{ij}$. Note that this solution only holds during the RD era, so our results cannot be extended to the phase transitions that happen in MD (Matter Dominated) or $\Lambda$D ($\Lambda$ Dominated) eras.

Now let’s discuss the power spectrum of GW. First of all, we define the equal-time two-point correlator of GW by

Where $\langle...\rangle$ denotes ensemble average and we will use spatial and time average to substitute ensemble average in the actual calculation. $(2\pi)^{3}$ is a coefficient putting at here out of convenience for future derivation. $\delta^{(3)}(\bm{k}-\bm{q})$ shows up out of the homogeneity of the system. Interested readers could refer to Chapter 7 of Ref.Maggiore_2008 to know more details and more properties of stochastic gravitational wave background.

We define the unequal-time correlator of energy-momentum tensor by²²2For simplicity, we call $T_{ij}$ and $\Psi_{ij}$ both as energy-momentum tensor:

we can expand the left hand side of Eq.(29) as:

Where we have defined $\bm{X}:=(\bm{x}+\bm{y})/2$ and $\bm{r}:=\bm{x}-\bm{y}$ to simplify our derivation. Note that the correlator itself doesn’t rely on $\bm{x}$ or $\bm{y}$, on the contrary, it only depends on $(\bm{x}-\bm{y})\equiv\bm{r}$. So we have

which tells us the factor $\Psi(\eta_{x},\eta_{y},k)$ appeared in Eq.(29) can be written as:

We know that the energy density of GW is given byMTW :

Note that the $t/\eta$ above means $t$ or $\eta$ rather than $t$ by $\eta$, and here $h_{ij}$ are all transverse and traceless. Now insert Eq.(26) and Eq.(27) into Eq.(28) and make an inverse Fourier transformation. Substitute the numerator of Eq.(33) by the result we’ve got, finally we can arrive at:

Now we define the dimensionless energy density per logarithmic comoving wave number of GW by:

The definition of $\rho_{\text{tot}}$ is unchanged comparing with the definition given in Sec.2. Since we consider the spatial curvature equals to zero, the total energy of the universe specifically equals to the energy needs to close the universe i.e. $\rho_{\text{cr}}=\rho_{\text{tot}}$. We assume that the phase transitions happen during RD era, so we can neglect the contribution of energy density by matter $\rho_{\text{mat}}$ and cosmological constant $\rho_{\Lambda}$ at here.

Substitute $\rho_{\text{GW}}$ by Eq.(34), we can rewrite $\Omega_{\text{GW}}$ as:

Here, we use superscript F to indicate FLRW spacetime. The dimensionless quantity $\sigma:=\frac{\widetilde{\beta}}{\mathcal{H}}=\frac{\beta}{H}$ denotes the fraction of speed of the phase transitions to the expansion speed of the universe as the same with the definition given in the previous section. Define $\tau$ as the effective duration of phase transitions, since $\tau$ itself is dependent on $\sigma$, different values of $\sigma$ will influence the upper limit of integral of $\eta_{x}$ and $\eta_{y}$ in Eq.(38), so $\Delta$ is also a function of $\sigma$.

At this stage, we must remind readers again that our definition of $\Delta^{\text{F}}$ is a little bit different with the definition of $\Delta$ given by the first line of Eq.(7). Our $\Delta^{\text{F}}$ is a function of $k/\widetilde{\beta},\sigma$ and $\eta$, while the $\Delta$ defined by Jinno and Takimoto is a function of $k/\beta$ and the speed of bubble wall $v$. The difference appears because of the different treatments of the duration of the phase transitions. In their previous work, they regard the phase transitions start at $t_{i}=-\infty$ and end at $t_{f}=+\infty$ in their calculation, which indeed simplify their analyses and decouple $\sigma$ from $\Delta$, but also could bring about a practical problem in latter numerical integration. In the actual integration, we cannot use $\infty$ as the upper limit, on the contrary, we have to artificially specify a finite value which is large enough to confirm the convergence of the integral. However, if we don’t know the specific relationship between $\tau$ and $\sigma$, then we may specify a too large number as upper limit of integration which can definitely lead to a longer time for computation. On the other hand, if we pick a too small value as the upper limit of integration, the result we get will be inaccurate and lack of value for reference. Besides the reasons above, there is another important motivation to consider a finite duration of phase transitions. Guo et al. consider a finite life time of the sources in their work to study the GW spectrum from sound waves during the first order cosmological phase transitions in an expanding universe.Guo_2020 As a result, they find an additional suppression factor $\Upsilon(y)$ should be included in the new expression of GW spectrum which could decrease GW spectrum significantly when the life time of sources is quite short. This point can be easily checked according to the FIG.15 of Ref.Guo_2020 . Be inspired by their work, we think it is necessary to take the finite duration of phase transitions into account. So in our work, $\tau$ itself is a function of $\sigma$ which induces that $\sigma$ cannot be decoupled from $\Delta^{\text{F}}$ anymore. As a result, the value of $\Delta^{\text{F}}$ in our work cannot compare with $\Delta$ in Ref.Jinno_2017 directly. In fact $\Delta^{\text{F}}$ and $\Delta^{\text{M}}$ are two quantities we will make a comparison later, where $\Delta^{\text{F}}$ represents our result and $\Delta^{\text{M}}$ represents Jinno’s result.

Using Friedmann EquationWeinberg :

$\Delta^{\text{F}}(k/\widetilde{\beta};\sigma;\eta)$ can be rewritten as:

Now, by comparing our Eq.(38) with Eq.(9), we can easily find the difference between these two equations out of two different choices of spacetime background. We’ll see the difference more clearly later from Eq.(57) and Eq.(64). At here, we adopt the Eq.(22) of Ref.Kosowsky_1994 to describe the functional relationship between $\kappa$ and $\alpha$ for latter numerical estimation as a benchmark:

Note that this equation only holds for the case of Jouguet detonation, readers can adopt the expressions given by Ref.Espinosa_2010 to do more estimations when $\alpha$ and $v$ take different values. In Ref.Espinosa_2010 , Espinosa et al. studied all of the bubble expansion regimes without specifying any particle physics model and gave the fitting formulas of $\kappa(v,\alpha)$ in their Appendix A. Note that since our derivation is only valid for phase transitions happen during the RD era, so we artificially restrict $\alpha\in[0,1]$. The work considering other eras and other ranges of $\alpha$ needs to be done in the future.

Now we can clearly discover that the most important quantity we need to focus on is $\Delta^{\text{F}}(k/\widetilde{\beta};\sigma;\eta)$. The coefficients $\kappa$ and $\alpha$ should be uniquely specified by choosing a specific model. In our work, we don’t choose any specific model and regard $\alpha$ as a free variable, $\kappa$ is a function of $\alpha$ described by Eq.(39). Through this approach, we can obtain a general understanding of the behavior of GW spectrum without being restricted by any specific model.

Note that the $\Omega_{\text{GW}}$ given by Eq.(35) denotes the dimensionless energy density of GW per logarithmic comoving wave number when GW has just been generated after the phase transitions rather than nowadays. Considering the effect of cosmological redshift yields:

where $k$ is the magnitude of comoving wave vector $\bm{k}$, $T_{\star}$ is the temperature after the phase transitions and $g_{\star}$ is the relativistic degrees of freedom in the universe corresponding to the temperature $T_{\star}$.

For the simplicity of comparison, we adopt the notation system defined in Ref.Jinno_2017 . The only difference are listed below. We define two new time variables $\Xi$ and $\xi$ rather than the choice of $T$ and $t_{\text{d}}$:

Besides these two variables, we only need to replace $t$ with $\eta$ in our derivation. The specific physical meaning of these quantities can refer to the FIG.3 and FIG.5 of Ref.Jinno_2017 . Out of the special property of FLRW metric, we find that light still move along the 45 degree line in the $\eta-x$ plane of the spacetime diagram, which indicates that we can still use the physical picture given by Jinno and Takimoto directly and make our analyses easier to follow.

According to Eq.(42) and Eq.(43) of Ref.Jinno_2017 , we have the probability of two spacetime points staying in the false vacuum state as:

Where the integral of $I(x,y)$ can be expanded as two terms:

After integration, we have $I(x,y)=I_{x}^{(y)}+I_{y}^{(x)}$:

Before we go ahead, let’s explain the reason why we need to discuss the probability of spacetime points staying in the false vacuum. Remember that under the so-called envelope approximation, any spatial point can be passed by a bubble for only once. Assuming that a bubble wall passes spatial point $\bm{x}$ at $\eta=\eta_{x}$, we must confirm that there are no bubbles nucleated inside the past cones of $x$, otherwise the spatial point $\bm{x}$ must have already transformed from the false vacuum state to the true vacuum state before $\eta=\eta_{x}$ and cannot be passed by another bubble by the second time. For the sake of the reason above, we have to figure out the probability of spacetime points staying in the false vacuum at the very first step of our calculation.

So far, we have every ingredient prepared to analyze the energy spectrum of GW. Just as the definition given by Jinno and Takimoto, $\Delta^{\text{F}}$ can be decomposed as $\Delta^{\text{F}(s)}+\Delta^{\text{F}(d)}$, where $\Delta^{\text{F}(s)}$ indicates the contribution to $\Delta^{\text{F}}$ from single bubble case and $\Delta^{\text{F}(d)}$ denotes the contribution from double bubble case. The specific meaning of these two cases have been discussed in Sec.2. Here we only simply introduce the calculation procedure taking the single bubble case for an example and directly give our results, the detailed calculation ideas can refer to the Appendix A of Ref.Jinno_2017 . Remember that our ultimate target is to get the GW spectrum which can be calculated by Eq.(35) and Eq.(38). In order to use these two equations we must figure out $\Psi(\eta_{x},\eta_{y},k)$ defined by Eq.(32) at first, which requires us to find all of the possible spacetime points for bubbles to nucleate which can lead to non-vanishing energy-momentum tensor at two spacetime points $x$ and $y$ to make $\langle\Psi_{ij}(x)\Psi_{ij}(y)\rangle$ nonzero. When we consider the single bubble case, we’ll easily find that only if the bubble nucleates at some special points which are belonging to the intersection part of the past cones of $x$ and $y$ (in Ref.Jinno_2017 , they call it $\delta V_{xy}$), may this situation contribute to $\Delta^{\text{F}(s)}$. By taking the possibility for bubbles to nucleate in these spacetime points into consideration and calculating the contributions given by every single configuration, we can write the two-point correlator of energy-momentum tensor as below:

where $F_{0},F_{1}$ and $F_{2}$ are given by: