1 Introduction

A detection of signatures of primordial gravitational waves (GWs) in the cosmic microwave background was claimed by the BICEP2 experiment in 2014 [1]. But this was shown to be likely caused by interstellar dust soon thereafter [2], and the search for true signatures of primordial GWs is still ongoing. The detection of the GW signatures, if confirmed, would gain the greatest importance, among other things, from its link to cosmic “inflation”: primordial GWs are seen as the smoking gun for the “Big Bang” expansion. According to the inflation theory, the early universe experienced an extreme burst of expansion, which lasted a tiny fraction of a second, but smoothed out irregularities–inhomogeneities, anisotropies, and the curvature of space, and made the universe appear homogeneous and isotropic [35].

It has been suggested that inflationary models of the early universe most likely lead to a “multiverse” [6, 7]. One such model is “eternal inflation” [8]: it proposes that many bubbles of spacetime individually nucleate and grow inside an ever-expanding background multiverse. The nucleation and growth of such bubbles can be modeled by a Coleman–de Luccia (CDL) instanton, a type of quantum transition between two classically disconnected vacua at different energies; the higher energy (false vacuum), the lower energy (true vacuum) [911]. A scalar field initially in the false vacuum state may tunnel quantum mechanically to the true vacuum state. This nucleates bubbles of the true vacuum (new phase) inside of the false vacuum (old phase) background; through a first-order phase transition. These bubbles then expand and collide with each other. The mechanism of bubble collisions can be effectively modeled by the CDL instanton: as bubbles continue to collide repeatedly, the scalar field transitions back and forth repeatedly between the false vacuum and the true vacuum, eventually settling down in the true vacuum as the collision process is gradually terminated.

From the viewpoints of physical cosmology, bubble collisions and GWs emitted from the collisions are interesting in the following contexts: (1) Our primordial inflation would be completed by a second-order (not by a first-order) phase transition. However, there is a possibility that some weaker inflation could occur after the primordial inflation; for example, “thermal inflation” [12]. It is quite probable that the thermal inflation is completed by a first-order phase transition, and therefore bubble collisions could take place through a CDL instanton. Then there would be some signatures of bubble collisions, which would presumably be carried via GWs [13]. (2) Suppose that we live in a single large true vacuum bubble and that the boundary of our bubble would collide with another bubble that is located outside our observable universe [14]. Then there would exist some signatures of bubble collisions and these could be carried via GWs. For scenario (1), the mechanism of bubble collision–GW emission should be modeled stochastically. However, for scenario (2), the mechanism can be well approximated by a two-bubble collision model.

There were numerous studies about bubble collisions and GWs emitted from the collisions. Among others, Hawking et al. [15] and Wu [16] studied the mechanism of the collision of two bubbles using the thin-wall approximation. Johnson et al. [17] and Hwang et al. [18] investigated the collision of two bubbles in full General Relativity via numerical computations. Kosowsky et al. [19] computed the GW spectrum resulting from two-bubble collisions in first-order phase transitions in flat spacetime using numerical simulations. Caprini et al. [20] developed a model for the bubble velocity power spectrum to calculate analytically the GW spectrum generated by two-bubble collisions in first-order phase transitions in flat spacetime.

In this paper, we focus on collisions of two equal-sized bubbles and compute GWs emitted from the collisions in time domain. Largely, our analysis proceeds in two steps through Sects. 2 and 3. In Sect. 2, we study the mechanism of bubble collisions by means of a real scalar field and a quartic potential of this field, building the simplest possible model for a CDL instanton. The Einstein equations and a scalar field equation are derived for this system and are solved simultaneously for the full General Relativistic treatment of the collision dynamics. Hwang et al. [18] is closely reviewed for this purpose. In Sect. 3, using the scalar field model from Sect. 2, we compute GWs from the bubble collisions in a straightforward manner. In the quadrupole approximation, time-domain gravitational waveforms are directly obtained by integrating the energy-momentum tensors over the volume of the wave sources, where the energy-momentum tensors are expressed in terms of the scalar field, the local geometry and the potential; therefore, containing all necessary information about the bubble collisions. Part of the computational results from Ref. [18] are recycled here to build the energy-momentum tensors. In parallel with the scalar field solutions in Sect. 2, which have been obtained with various false vacuum field values [18], we present gravitational waveforms emitted during

(i) the initial-to-intermediate stage of strong collisions, and

(ii) the final stage of weak collisions.

The former is obtained numerically, in full General Relativity and the latter analytically, in the flat spacetime approximation. The thin-wall and quadrupole approximations are assumed to simplify our analysis and the next-to-leading order corrections beyond these approximations are disregarded. However, the approximations serve our purpose well: we aim to provide qualitative illustrations of the time-domain gravitational waveforms from the bubble collisions, which will be useful for constructing the templates for observation in the future. We adopt the unit convention, \(c=G=1\) for all our computations of GWs.

Fig. 1
figure 1

The quartic potential \(V( S) =\frac{3V_{ \mathrm {f}}}{\beta S_\mathrm{{f}}^{4}}S^{4}-\frac{2( \beta +3) V_\mathrm{{f}}}{\beta S_\mathrm{{f}}^{3}}S^{3}+\frac{ 3( \beta +1) V_\mathrm{{f}}}{\beta S_\mathrm{{f} }^{2}}S^{2}\); expressed in terms of the false vacuum field \(S_\mathrm{{f}}= \sqrt{4\pi }\phi _\mathrm{{f}}\), the vacuum energy of the false vacuum \(V_\mathrm{{f}}\) and a free parameter \(\beta \). The patterns of the potential are shown for \(\beta =0.1\), \(V_{ \mathrm {f}}=10^{-4}\) and \(S_\mathrm{{f}}=\sqrt{4\pi } \phi _\mathrm{{f}}=\) (1) \(0.1\), (2) \(0.2\), (3) \(0.3\), (4) \(0.4\) (credit: Hwang et al. [18])

2 Gravity–scalar field dynamics for colliding bubbles

The mechanism of two equal-sized colliding bubbles can be effectively modeled by means of a CDL instanton [810]. Basically, one can build a model for this, which consists of gravitation, a real scalar field and a potential of the field. In this section we introduce one such model from Hwang et al. [17], which is built with a quartic potential, the simplest possible one for the CDL instanton.

2.1 Dynamics of bubble collisions

A system of Einstein gravity coupled with a scalar field that is governed by a potential can be described by the following action:

$$\begin{aligned} \mathcal {S}=\int \mathrm{d}^{4}x\sqrt{-g}\left[ \frac{1}{16\pi }R-\frac{1}{2} \nabla _{\mu }\phi \nabla ^{\mu }\phi -V\left( \phi \right) \right] , \end{aligned}$$
(1)

where \(R\) denotes the Ricci scalar, \(\phi \) the scalar field, and \( V(\phi )\) the potential of the scalar field [18]. From this system the Einstein equations are derived:

$$\begin{aligned} R_{\mu \nu }-\frac{1}{2}Rg_{\mu \nu }=8\pi T_{\mu \nu }, \end{aligned}$$
(2)

where the energy-momentum tensors on the right-hand side are written as

$$\begin{aligned} T_{\mu \nu }=\phi _{;\mu }\phi _{;\nu }-\frac{1}{2}\phi _{;\rho }\phi _{;\sigma }g^{ \rho \sigma }g_{\mu \nu } -V(\phi )g_{\mu \nu }. \end{aligned}$$
(3)

Also, the scalar field equation for the system reads

$$\begin{aligned} \nabla ^{2}\phi =\frac{\mathrm{d}V}{\mathrm{d}\phi }. \end{aligned}$$
(4)

The Einstein equations (2) and the scalar field equation (4) constitute a scalar field model that effectively describes the mechanism of two colliding bubbles in curved spacetime [18]. Given a potential \( V( \phi ) \), the scalar field solution \(\phi \) and the geometry solution \(g_{\mu \nu }\) should be obtained by solving Eqs. (2) and (4) simultaneously.Footnote 1 To this end, we prescribe an ansatz for the geometry \(g_{\mu \nu }\) with the hyperbolic symmetry, using the double-null coordinates:

$$\begin{aligned} \mathrm{d}s^{2}=-\alpha _\mathrm{{h}}^{2}( u,v) \mathrm{d}u\mathrm{d}v+r_\mathrm{{h} }^{2}( u,v) \mathrm{d}H^{2}, \end{aligned}$$
(5)

where \(\mathrm{d}H^{2}=\mathrm{d}\chi ^{2}+\sinh ^{2}\chi \mathrm{d}\theta ^{2}\) with \(0\le \chi <\infty \), \(0\le \theta <2\pi \) [18], and \(\alpha _\mathrm{{h} }( u,v) \) and \(r_\mathrm{{h}}( u,v) \) are to be determined by solving Eqs. (2) and (4) simultaneously in the coordinates \(( u,v,\chi ,\theta ) \). In the flat spacetime limit, the double-null coordinates are defined as \(u\equiv \tau -x\) and \(v\equiv \tau +x\) with \(\tau ^{2}\equiv t^{2}-y^{2}-z^{2}\), \(t=\tau \cosh \chi \), \( y=\tau \sinh \chi \sin \theta \), \(z=\tau \sinh \chi \cos \theta \): in our analysis, the \(x\)-axis of Cartesian coordinates is chosen to coincide with a line adjoining the centers of the two bubbles, and the \(y\)-axis and the \(z\)-axis lie in a plane perpendicular to the \(x\)-axis.

2.2 Solving the scalar field equation

To build the simplest model of the CDL instanton for two identical colliding bubbles, we consider the potential in Sect. 2.1 to be

$$\begin{aligned} V( \phi ) =\frac{p_{4}}{4!}\phi ^{4}+\frac{p_{3}}{3!}\phi ^{3}+ \frac{p_{2}}{2!}\phi ^{2}, \end{aligned}$$
(6)

where \(p_{2}\), \(p_{3}\), and \(p_{4}\) are constants which can be appropriately chosen to tune the shape of the potential. Bubble collisions are represented by the scalar field moving along this potential: the field initially in the false vacuum state (at higher local minimum of potential) tunnels quantum mechanically to the true vacuum state (at lower local minimum of potential), repeating the transitions back and forth between the two states, eventually settling down in the true vacuum state.

Following Ref. [18], we may rescale the scalar field, \(S\equiv \sqrt{ 4\pi }\phi \) for computational convenience, and can specify \(p_{2}\), \(p_{3}\), and \(p_{4}\) in terms of the false vacuum field \(S_\mathrm{{f}}=\sqrt{4\pi } \phi _\mathrm{{f}}\), the vacuum energy of the false vacuum \(V_\mathrm{{f}}\) and a free parameter \(\beta \). The potential (6) can then be rewritten as

$$\begin{aligned} V(S)=\frac{3V_\mathrm{{f}}}{\beta S_\mathrm{{f}}^{4}}S^{4}-\frac{ 2(\beta +3)V_\mathrm{{f}}}{\beta S_\mathrm{{f}}^{3}}S^{3}+\frac{ 3(\beta +1)V_\mathrm{{f}}}{\beta S_\mathrm{{f}}^{2}}S^{2}. \end{aligned}$$
(7)

With this potential the scalar field equation (4), which is now rescaled, reads

$$\begin{aligned} \nabla ^{2}S=\frac{12V_\mathrm{{f}}}{\beta S_\mathrm{{f}}^{4}}S^{3}-\frac{ 6(\beta +3)V_\mathrm{{f}}}{\beta S_\mathrm{{f}}^{3}}S^{2}+\frac{ 6(\beta +1)V_\mathrm{{f}}}{\beta S_\mathrm{{f}}^{2}}S. \end{aligned}$$
(8)

This is a non-linear wave equation whose analytical solution is not generally known: we normally approach this type of problem with numerical methods.

Now, we solve the scalar field equation (8) simultaneously with the Einstein equations (2), using the ansatz given by (5), in the coordinates \(( u,v,\chi ,\theta )\). However, it turns out that our scalar field solution is independent of the coordinates \(\chi \) and \( \theta \) and is expressed in the coordinates \(( u,v)\) only; namely, \(S( u,v) =\sqrt{4\pi }\phi ( u,v)\) [15, 18, 19]. With the choice of the constants, \( \beta =0.1\), \(V_\mathrm{{f}}=10^{-4}\) and \(S_\mathrm{{f}}=\) (1) \(0.1\), (2) \( 0.2\), (3) \(0.3\), (4) \(0.4\) in Eq. (7), the potential takes the forms as given by Fig. 1 [18]. With this potential, our numerical solution \(S( u,v) \) is obtained as presented by Fig. 2 [18]. In each case of \(S_\mathrm{{f}}\), (1)–(4), the bubble wall has a different value of tension due to a different value of \(S_\mathrm{{f}}\) as shown by Fig. 1. In the top left of Fig. 2 the bubble has the lowest tension while in the bottom right it has the highest tension among the four cases of \(S_\mathrm{{f}}\). This results in the wall crossing regions in the top left being relatively wider than those in the bottom right.

Fig. 2
figure 2

The numerical solutions \(S(u,v)= \sqrt{4\pi }\phi (u,v)\) obtained with the potential from Fig. 1; with various false vacuum field values, \(S_{ \mathrm {f}}=\sqrt{4\pi }\phi _\mathrm{{f}}=\) (1) \(0.1\), (2) \(0.2\), (3) \(0.3\), (4) \(0.4\) (credit: Hwang et al. [18])

3 Gravitational waves from bubble collisions

In Sect. 2 we have built a system of two equal-sized colliding bubbles in curved spacetime by means of a CDL instanton model, considering the potential given by Eq. (7) [18]. Now, we consider that the a present observer lives in the true vacuum region of one of the two bubbles and that signatures of bubble collisions which took place in the distant past are being carried to the present observer via GWs. Here we assume that the distance between the center of collision region and the observer can be arbitrarily large (within the size of our universe), and thus that the time for a collision event, which is the retarded time to the present observer, can be quite far in the distant past; namely, \(t_\mathrm{{R }}=t_\mathrm{{P}}-r/c\ll t_\mathrm{{P}}\), where \(t_\mathrm{{R}}\) denotes the retarded time, \(t_\mathrm{{P}}\) the present time and \(r\) the distance. Therefore, our bubble collision may be regarded as a localized event, as long as the observer is reasonably far away from the collision region. In Fig. 3 a causal relationship between two colliding bubbles and an observer is depicted using a null-cone. Here collision events that took place in the distant past are placed within the intersection of the timelike zone (green-colored region) of a null-cone (green dashed lines) and the ‘diamond’ zone (region enclosed by black dashed lines): all the collision events as our GW sources, namely the collisions in Fig. 2 (strong collisions in the initial-to-intermediate stage) and the collisions to be discussed in Sect. 3.2 later (weak collisions in the final stage) should be considered to have taken place within this intersection.Footnote 2

In the above scenario, our GWs from bubble collisions can be computed in a straightforward manner: by integrating the energy-momentum tensors combined with a Green function over the volume of the wave sources, where the energy-momentum tensors are expressed in terms of the scalar field, the local geometry and the potential by means of Eqs. (3), (5), and (8); therefore, they contain all necessary information about the bubble collisions. A mathematical description of this computation is given as follows. In the transverse trace-free (TT) gauge, GWs, as derived from the perturbed Einstein equations in linearized gravity, can be expressed as

$$\begin{aligned} h_{ij}^\mathrm{{TT}}( t,\mathbf {x})= & {} \overline{h}_{ij}^\mathrm{{TT }}( t,\mathbf {x}) =\frac{4G}{c^{4}}\Lambda _{ij,kl}( \mathbf { n})\nonumber \\&\times \int \mathrm{d}^{3}x^{\prime }\frac{1}{\vert \mathbf {x-x}^{\prime }\vert }T^{kl}\left( t-\frac{\vert \mathbf {x-x}^{\prime }\vert }{c},\mathbf {x}^{\prime }\right) , \end{aligned}$$
(9)

where \(\overline{h}_{ij}\equiv h_{ij}-\frac{1}{2}\delta _{ij}h_{k}^{k}\) and the unit vector \(\mathbf {n}\) denotes the propagation direction of the waves, and the projection tensor for gravitational radiation,

$$\begin{aligned} \Lambda _{ij,kl}( \mathbf {n}) \equiv P_{ik}P_{jl}-\frac{1}{2} P_{ij}P_{kl} \end{aligned}$$
(10)

with

$$\begin{aligned} P_{ij}=\delta _{ij}-n_{i}n_{j}. \end{aligned}$$
(11)
Fig. 3
figure 3

A causal relationship between two colliding bubbles and an observer: gravitational waves from bubble collisions can be observed by a distant observer (red spot) who is in the causal future of a collision event (black spot). The dotted vertical line inside the timelike zone (green-colored region) of a null-cone (green dashed lines) represents the time-axis, along which bubble collisions took place in the distant past; within the intersection of the timelike zone of the null-cone and the ‘diamond’ zone (region enclosed by black dashed lines). \(u\) and \(v\) represent the double-null coordinates as defined in Sect. 2. T and F denote a true vacuum and a false vacuum, respectively

We find that the computation would be technically quite difficult with the integral as it is in Eq. (9): the way the source point \(\mathbf {x} ^{\prime }\) (integration variable) is combined with the field point \(\mathbf { x}\) in the integrand would make our calculation quite intractable. However, expressing the integral in expansion, we obtain a more computationally favorable form:

$$\begin{aligned}&h_{ij}^\mathrm{{TT}}( t,\mathbf {x}) =\frac{4G}{c^{4}r}\Lambda _{ij,kl}( \mathbf {n}) \bigg [ \int \mathrm{d}^{3}x^{\prime }T^{kl}( t_{ \mathrm {R}},\mathbf {x}^{\prime }) \nonumber \\&\quad +\frac{1}{c}n^{m}\frac{\mathrm{d}}{\mathrm{d}t_\mathrm{{R}}}\int \mathrm{d}^{3}x^{\prime }T^{kl}( t_\mathrm{{R}},\mathbf {x}^{\prime }) x_{m}^{\prime } \nonumber \\&\quad +\frac{1}{2c^{2}}n^{m}n^{p}\frac{\mathrm{d}^{2}}{\mathrm{d}t_\mathrm{{R}}^{2}}\int \mathrm{d}^{3}x^{\prime }T^{kl}( t_\mathrm{{R}},\mathbf {x}^{\prime }) x_{m}^{\prime }x_{p}^{\prime }+\cdots \bigg ] _{t_\mathrm{{R}}=t-r/c},\nonumber \\ \end{aligned}$$
(12)

where \(r=\vert \mathbf {x}\vert \) and \(t_\mathrm{{R}}=t-r/c\) denotes the retarded time. In particular, the computation resulting from the first term alone in the square bracket in Eq. (12) is called the “quadrupole approximation”. The next terms will provide corrections to this computation.

3.1 Computation of gravitational waves in the quadrupole approximation

The complete information about the motion of the colliding two-bubble system is encoded in the scalar field solution \(S=\sqrt{4\pi }\phi \), as given by Fig. 2, and thus is carried by the energy-momentum tensors through Eq. (3): to be precise, the energy-momentum tensors are comprised of the scalar field \(S\) and the geometry \(g_{\mu \nu }\), which are obtained by solving Eqs. (2) and (4) simultaneously [18]. As described by Eqs. (9) and (12), GWs from the system are computed with the energy-momentum tensors being the sources. It is believed that the two bubbles will be in highly relativistic motion when they collide [15]. In view of this, corrections due to the next-to-leading order terms in Eq. (12) should not be disregarded if one aims to compute GWs from the system accurately. However, although not perfectly accurate, the leading order term alone in Eq. (12) provides the “quadrupole approximation” of GWs:

$$\begin{aligned} _\mathrm{{Q}}h_{ij}^\mathrm{{TT}}( t,\mathbf {x}) =\frac{4}{r} \Lambda _{ij,kl}( \mathbf {n}) I_{kl}( t_\mathrm{{R}}), \end{aligned}$$
(13)

where we have adopted the unit convention \(c=G=1\), and

$$\begin{aligned} I_{kl}( t_\mathrm{{R}}) \equiv \int \mathrm{d}^{3}x^{\prime }T^{kl}( t_\mathrm{{R}},\mathbf {x}^{\prime }). \end{aligned}$$
(14)

Throughout this paper our computation is carried out only based on this piece. Our main purpose is to provide qualitative insights into patterns of GWs from the colliding two-bubble system in time domain, and the next-to-leading order corrections in Eq. (12) are disregarded in our analysis.

Following Ref. [19], we can reduce the amount of computation in a great deal. As described in Sect. 2.1, we choose the \(x\)-axis to coincide with the line adjoining the centers of the two bubbles. With the axial symmetry about the \(x\)-axis, the off-diagonal components are zero and we can put \(I_{kl}\) in the form

$$\begin{aligned} I_{kl}=D\delta _{kl}+\triangle \delta _{kx}\delta _{lx}. \end{aligned}$$
(15)

Here, the first term turns out to be

$$\begin{aligned} D=\frac{1}{2}(I_{yy}+I_{zz}), \end{aligned}$$
(16)

which does not contribute to gravitational radiation due to Eqs. (10) and (11). The second term is given by

$$\begin{aligned} \triangle =I_{xx}-\frac{1}{2}(I_{yy}+I_{zz}). \end{aligned}$$
(17)

Therefore, \(I_{kl}\) is practically equivalent to \(\triangle \delta _{kx} \delta _{lx}\):

$$\begin{aligned} I_{kl}\sim \delta _{kx}\delta _{lx}\left[ I_{xx}-\frac{1}{2}(I_{yy}+I_{zz} )\right] . \end{aligned}$$
(18)

Then by Eqs. (14) and (18) we may express

$$\begin{aligned}&I_{kl}(t_\mathrm{{R}})=\delta _{kx}\delta _{lx}\int \mathrm{d}^{3}x^{\prime }\nonumber \\&\times \left[ T^{xx}(t_\mathrm{{R}},\mathbf {x}^{\prime })-\frac{1}{2} [T^{yy}(t_\mathrm{{R}},\mathbf {x}^{\prime })+T^{zz} (t_{ \mathrm {R}},\mathbf {x}^{\prime })]\right] . \end{aligned}$$
(19)

Now, recall from Sect. 2.1 that in the flat spacetime we define the hyperbolic coordinates \(\tau \), \(\chi \), \(\theta \) by

$$\begin{aligned}&t_\mathrm{{R}} = \tau \cosh \chi , \end{aligned}$$
(20)
$$\begin{aligned}&y =\tau \sinh \chi \sin \theta , \end{aligned}$$
(21)
$$\begin{aligned}&z =\tau \sinh \chi \cos \theta , \end{aligned}$$
(22)

so that

$$\begin{aligned} \tau ^{2}=t_\mathrm{{R}}^{2}-\rho ^{2}, \end{aligned}$$
(23)

where \(0\le \chi <\infty \), \(0\le \theta <2\pi \) and \(\rho \equiv \sqrt{ y^{2}+z^{2}}\). In these coordinates the flat spacetime metric takes the form

$$\begin{aligned} \mathrm{d}s^{2}=-\mathrm{d}\tau ^{2}+\mathrm{d}x^{2}+\tau ^{2}( \mathrm{d}\chi ^{2}+\sinh ^{2}\chi \mathrm{d}\theta ^{2}). \end{aligned}$$
(24)

In this geometry, however, the scalar field solution is independent of the coordinates \(\chi \) and \(\theta \) and is expressed in the coordinates \( ( \tau ,x) \) only; namely, \(\phi ( \tau ,x) \) [15, 19]. Taking this into account, we should rewrite the volume element for the integral (19) as

$$\begin{aligned} \mathrm{d}^{3}x=\mathrm{d}x\rho \mathrm{d}\rho \mathrm{d}\theta =-\mathrm{d}x\tau \mathrm{d}\tau \mathrm{d}\theta , \end{aligned}$$
(25)

which is defined at the instant \(t_\mathrm{{R}}\) by means of Eq. (23).

From Eq. (23) we see that \(\tau \) has an upper bound \(t_\mathrm{{R}}\) with \(\rho =0\). This represents the exterior surface of the bubble walls, i.e.

$$\begin{aligned} \tau =t_\mathrm{{R}}. \end{aligned}$$
(26)

However, the interior surface is found from Eq. (23) to be

$$\begin{aligned} \tau =t_\mathrm{{R}}-\frac{\eta ^{2}}{8t_\mathrm{{R}}}+\mathcal {O}\left( \frac{\eta ^{4}}{t_\mathrm{{R}}^{3}}\right) , \end{aligned}$$
(27)

given a wall thickness \(\eta \sim 2\rho \ll t_\mathrm{{R}}\). Then from Eqs. (19), (25), (26), and (27) we can compute \( I_{kl}( t_\mathrm{{R}}) \) effectively out of a volume piece \( \mathcal {V}\):

$$\begin{aligned} I_{kl}( t_\mathrm{{R}})= & {} \delta _{kx}\delta _{lx}\int _{\mathcal { V}}\mathrm{d}^{3}x^{\prime }\left[ T_{xx}-\frac{1}{2}( T_{yy}+T_{zz}) \right] \nonumber \\= & {} 2\pi \delta _{kx}\delta _{lx}\int _{t_\mathrm{{R}}-\eta ^{2}/( 8t_{ \mathrm {R}}) }^{t_\mathrm{{R}}}\mathrm{d}\tau ^{\prime }\tau ^{\prime }\int _{-x_\mathrm{{o}}}^{x_\mathrm{{o}}}\mathrm{d}x^{\prime }\nonumber \\&\times \left[ T_{xx}-\frac{1}{2} ( T_{yy}+T_{zz}) \right] +\mathcal {O}\left( \frac{\eta ^{4}}{t_{ \mathrm {R}}^{3}}\right) \nonumber \\= & {} \frac{\pi }{4}\delta _{kx}\delta _{lx}\eta ^{2}\int _{-x_\mathrm{{o}}}^{x_{ \mathrm {o}}}\mathrm{d}x^{\prime }\left[ T_{xx}-\frac{1}{2}( T_{yy}+T_{zz}) \right] _{\tau =t_\mathrm{{R}}}\nonumber \\&+\,\mathcal {O}\left( \frac{\eta ^{4}}{t_{ \mathrm {R}}^{3}}\right) , \end{aligned}$$
(28)

where the volume piece \(\mathcal {V}\) is defined from a thin cylindrical shell in motion, having the thickness \(\eta ^{2}/( 8t_\mathrm{{R} })\), extending along the \(x\)-axis, by means of Eqs. (26) and (27): from \(t_\mathrm{{R}}-\eta ^{2}/( 8t_\mathrm{{R}}) + \mathcal {O}( \eta ^{4}/t_\mathrm{{R}}^{3}) \le \tau \le \) \(t_{ \mathrm {R}}\) we find \(\mathcal {V}=\Delta x\vert \tau \Delta \tau \vert \Delta \theta =\frac{\pi }{4}\eta ^{2}\Delta x[ 1+\mathcal {O }( \eta ^{2}/t_\mathrm{{R}}^{2}) ] \), and the limit of the integral \(x_\mathrm{{o}}=\Delta x/2\) should be chosen to be sufficiently large such that collision effects be fully covered in numerical integration.Footnote 3 Following Ref. [15], we estimate a bubble wall thickness \(\eta \), assuming that the walls will be highly relativistic when they collide, having the Lorentz factor \(\gamma \):

$$\begin{aligned} \eta \sim \frac{3}{2}\phi _\mathrm{{f}}/( \xi ^{2}\gamma ) \sim \frac{3}{2}\phi _\mathrm{{f}}^{2}/( b\epsilon ^{4}), \end{aligned}$$
(29)

where \(\phi _\mathrm{{f}}\) denotes the scalar field value at the false vacuum and \(\xi ^{4}\) the effective height of the potential barrier between the two minima, and the Lorentz factor \(\gamma =b\epsilon ^{4}/( \xi ^{2}\phi _\mathrm{{f}}) \) with \(2b\) representing the separation of the bubbles and \(\epsilon ^{4}\) the potential difference between the two minima (which is equivalent to \(V_\mathrm{{f}}\) in our analysis in Sect. 2.2) [15].

However, as described in Sect. 2.2, our scalar field \(S=\sqrt{ 4\pi }\phi \) is obtained by solving Eqs. (2) and (8) simultaneously, using the ansatz (5), in the coordinates \(( u,v,\chi ,\theta ) \). Then by Eq. (3) the energy-momentum tensors should be expressed in the same coordinates. Now, due to the definitions of \(u\) and \(v\) in the flat spacetime limit, and by Eqs. (20), (21), and (22) we have

$$\begin{aligned}&u =\tau -x=\sqrt{t_\mathrm{{R}}^{2}-( y^{2}+z^{2}) }-x, \end{aligned}$$
(30)
$$\begin{aligned}&v =\tau +x=\sqrt{t_\mathrm{{R}}^{2}-( y^{2}+z^{2}) }+x, \end{aligned}$$
(31)
$$\begin{aligned}&\chi =\tanh ^{-1}\left( \sqrt{\frac{y^{2}+z^{2}}{t_\mathrm{{R}}^{2}}} \right) . \end{aligned}$$
(32)

Using these relations, we find

$$\begin{aligned} T_{xx}= & {} T_{uu}\left( \frac{\partial u}{\partial x}\right) ^{2}+2T_{uv} \frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\,T_{vv}\left( \frac{\partial v}{\partial x}\right) ^{2}\nonumber \\&+\,T_{\chi \chi }\left( \frac{ \partial \chi }{\partial x}\right) ^{2} =T_{uu}-2T_{uv}+T_{vv},\end{aligned}$$
(33)
$$\begin{aligned} T_{yy}= & {} T_{uu}\left( \frac{\partial u}{\partial y}\right) ^{2}+2T_{uv} \frac{\partial u}{\partial y}\frac{\partial v}{\partial y}\nonumber \\&+\,T_{vv}\left( \frac{\partial v}{\partial y}\right) ^{2}+T_{\chi \chi }\left( \frac{ \partial \chi }{\partial y}\right) ^{2} \nonumber \\= & {} \frac{y^{2}}{t_\mathrm{{R}}^{2}-( y^{2}+z^{2}) }(T_{uu}+2T_{uv}+T_{vv}) \nonumber \\&+\,\frac{t_\mathrm{{R}}^{2}y^{2}}{[ t_{ \mathrm {R}}^{2}-( y^{2}+z^{2}) ] ^{2}( y^{2}+z^{2}) }T_{\chi \chi },\end{aligned}$$
(34)
$$\begin{aligned} T_{zz}= & {} T_{uu}\left( \frac{\partial u}{\partial z}\right) ^{2}+2T_{uv} \frac{\partial u}{\partial z}\frac{\partial v}{\partial z}+T_{vv}\left( \frac{\partial v}{\partial z}\right) ^{2}\nonumber \\&+\,T_{\chi \chi }\left( \frac{ \partial \chi }{\partial z}\right) ^{2} \nonumber \\= & {} \frac{z^{2}}{t_\mathrm{{R}}^{2}-( y^{2}+z^{2}) }( T_{uu}+2T_{uv}+T_{vv})\nonumber \\&+\,\frac{t_\mathrm{{R}}^{2}z^{2}}{[ t_{ \mathrm {R}}^{2}-( y^{2}+z^{2}) ] ^{2}( y^{2}+z^{2}) }T_{\chi \chi }. \end{aligned}$$
(35)

Substituting Eqs. (33), (34), and (35) into Eq. (28), we obtain

$$\begin{aligned} I_{kl}( t_\mathrm{{R}})= & {} \frac{\pi }{4}\delta _{kx}\delta _{lx}\eta ^{2}\int _{-x_\mathrm{{o}}}^{x_\mathrm{{o}}}\mathrm{d}x^{\prime }\nonumber \\&\times \left[ T_{uu}-2T_{uv}+T_{vv}-\frac{1}{2t_\mathrm{{R}}^{2}}T_{\chi \chi }\right] _{t_\mathrm{{R}}}+\,\mathcal {O}\left( \frac{\eta ^{4}}{t_\mathrm{{R}}^{3}} \right) ,\nonumber \\ \end{aligned}$$
(36)

where the subscript \(t_\mathrm{{R}}\) outside the square bracket means that the double-null coordinates \(( u,v) \) are defined at \(\tau =t_{ \mathrm {R}}\); namely, \(u=t_\mathrm{{R}}-x\) and \(v=t_\mathrm{{R}}+x\). In the actual computation of Eq. (36), we integrate \(T_{uu}( u,v) \), \(T_{uv}( u,v) \), \(T_{vv}( u,v) \) and \(T_{\chi \chi }( u,v) \), which are constructed out of the scalar field solution \(S( u,v) =\sqrt{4\pi }\phi ( u,v) \), the geometry solution \(g_{uv}( u,v) \), \(g_{\chi \chi }( u,v) \), \( g_{\theta \theta }( u,v) \) and the potential \(V( S) \) via Eq. (3). Then we need to change the variable of integration, from \( x\) to \(u\) or \(v\). Using the relations \(u=t_\mathrm{{R}}-x\) and \(v=t_\mathrm{{ R}}+x\), we can convert

$$\begin{aligned} \mathrm{d}x=-\mathrm{d}u\,\,\text { or }\,\,\mathrm{d}x=\mathrm{d}v. \end{aligned}$$
(37)

Then we may rewrite

$$\begin{aligned} \int _{-x_\mathrm{{o}}}^{x_\mathrm{{o}}}\mathrm{d}x^{\prime }~T_{ab}( u,v)= & {} \int _{t_\mathrm{{R}}-x_\mathrm{{o}}}^{t_\mathrm{{R}}+x_\mathrm{{o} }}\mathrm{d}u~T_{ab}( u,2t_\mathrm{{R}}-u)\nonumber \\= & {} \int _{t_\mathrm{{R}}-x_{ \mathrm {o}}}^{t_\mathrm{{R}}+x_\mathrm{{o}}}\mathrm{d}v~T_{ab}( 2t_\mathrm{{R} }-v,v) \nonumber \\= & {} 2\int _{0}^{x_\mathrm{{o}}}\mathrm{d}u~T_{ab}( t_\mathrm{{R}}+u,t_\mathrm{{R} }-u) \nonumber \\= & {} 2\int _{0}^{x_\mathrm{{o}}}\mathrm{d}v~T_{ab}( t_\mathrm{{R}}-v,t_{ \mathrm {R}}+v), \end{aligned}$$
(38)

where \(T_{ab}\) represents any of \(T_{uu}\), \(T_{uv}\), \(T_{vv}\) and \(T_{\chi \chi }\), and the expressions in the second line have been obtained via translations, \(u\rightarrow u-t_\mathrm{{R}}\) and \(v\rightarrow v-t_\mathrm{{ R}}\). Then by Eqs. (36) and (38) \(I_{kl}( t_\mathrm{{R} }) \) can be expressed as

$$\begin{aligned} I_{kl}( t_\mathrm{{R}})= & {} \frac{\pi }{2}\delta _{kx}\delta _{lx}\eta ^{2}\bigg [ \int _{0}^{x_\mathrm{{o}}}\mathrm{d}u~T_{uu}( t_\mathrm{{R} }\pm u,t_\mathrm{{R}}\mp u) \nonumber \\&-\,2\int _{0}^{x_\mathrm{{o}}}\mathrm{d}u~T_{uv}( t_\mathrm{{R}}\pm u,t_\mathrm{{R}}\mp u) \nonumber \\&+\int _{0}^{x_\mathrm{{o}}}\mathrm{d}u~T_{vv}( t_\mathrm{{R}}\pm u,t_\mathrm{{R}}\mp u) \nonumber \\&-\frac{1}{2t_\mathrm{{R}}^{2}}\int _{0}^{x_{ \mathrm {o}}}\mathrm{d}u~T_{\chi \chi }( t_\mathrm{{R}}\pm u,t_\mathrm{{R}}\mp u) \bigg ] +\mathcal {O}\left( \frac{\eta ^{4}}{t_\mathrm{{R}}^{3}} \right) . \end{aligned}$$
(39)
Fig. 4
figure 4

The numerical plots of \(_\mathrm{{Q}}h^\mathrm{{TT} }(t)r/(2\pi \eta ^{2})\) with various false vacuum field values, \(S_\mathrm{{f}}=\sqrt{4\pi } \phi _\mathrm{{f}}=\) (1) \(0.1\), (2) \(0.2\), (3) \(0.3\), (4) \(0.4\). As \( \eta ^{2}\sim \phi _\mathrm{{f}}^{4}\), the amplitude of \(_{ \mathrm {Q}}h^\mathrm{{TT}}(t)\) should scale as (1) \(1\), (2) \( 2^{4} \), (3) \(3^{4}\), (4) \(4^{4}\)

If the wall thickness \(\eta \) can be taken sufficiently small in Eq. (39), then by Eq. (13) we can compute the bubble-collision-induced GWs in the quadrupole approximation as

$$\begin{aligned} _\mathrm{{Q}}h_{ij}^\mathrm{{TT}}(t,\mathbf {x})\approx & {} \frac{ 2\pi }{r}\eta ^{2}\Lambda _{ij,kl}(\mathbf {n})\delta _{kx}\delta _{lx} \nonumber \\&\times \bigg [\int _{0}^{x_\mathrm{{o}}}\mathrm{d}u~T_{uu}(t_\mathrm{{R} }\pm u,t_\mathrm{{R}}\mp u)\nonumber \\&-\,2\int _{0}^{x_\mathrm{{o} }}\mathrm{d}u~T_{uv}(t_\mathrm{{R}}\pm u,t_\mathrm{{R}}\mp u) \nonumber \\&+\int _{0}^{x_\mathrm{{o}}}\mathrm{d}u~T_{vv}(t_\mathrm{{R}}\pm u,t_{ \mathrm {R}}\mp u)\nonumber \\&-\frac{1}{2t_\mathrm{{R}}^{2}}\int _{0}^{x_\mathrm{{o} }}\mathrm{d}u~T_{\chi \chi }(t_\mathrm{{R}}\pm u,t_\mathrm{{R}}\mp u)\bigg ]. \end{aligned}$$
(40)

Now, without loss of generality we may choose

$$\begin{aligned} \mathbf {n=}{(n_{x},n_{y},n_{z})=(\cos \vartheta ,\sin \vartheta ,0)}, \end{aligned}$$
(41)

where \(\vartheta \) denotes the angle of propagation taken from the \(x\)-axis. From this it follows that

$$\begin{aligned} \Lambda _{ij,kl}(\mathbf {n})\delta _{kx}\delta _{lx}= & {} \delta _{ix} \delta _{jx}-2\delta _{ix}n_{j} \cos \vartheta \nonumber \\&+\,\frac{1}{2}n_{i}n_{j}\left( 1+ \cos ^{2}\vartheta \right) -\frac{1}{2}\delta _{ij}\sin ^{2}\vartheta , \end{aligned}$$
(42)

due to Eqs. (10) and (11). Substituting this into Eq. (40), we finally find the expression

$$\begin{aligned} _\mathrm{{Q}}h_{ij}^\mathrm{{TT}}(t,\mathbf {x})\approx & {} \frac{ 2\pi }{r}\eta ^{2}\bigg [\delta _{ix}\delta _{jx}-2\delta _{ix}n_{j}\cos \vartheta \nonumber \\&+\,\frac{1}{2}n_{i}n_{j} \left( 1+\cos ^{2}\vartheta \right) -\frac{1}{2} \delta _{ij}\sin ^{2}\vartheta \bigg ] \nonumber \\&\times \bigg [\int _{0}^{x_\mathrm{{o}}}\mathrm{d}u~T_{uu}(t_\mathrm{{R} }\pm u,t_\mathrm{{R}}\mp u)\nonumber \\&-2\int _{0}^{x_\mathrm{{o} }}\mathrm{d}u~T_{uv}(t_\mathrm{{R}}\pm u,t_\mathrm{{R}}\mp u) \nonumber \\&+\int _{0}^{x_\mathrm{{o}}}\mathrm{d}u~T_{vv}(t_\mathrm{{R}}\pm u,t_{ \mathrm {R}}\mp u)\nonumber \\&-\frac{1}{2t_\mathrm{{R}}^{2}}\int _{0}^{x_\mathrm{{o} }}\mathrm{d}u~T_{\chi \chi }(t_\mathrm{{R}}\pm u,t_\mathrm{{R}}\mp u)\bigg ], \end{aligned}$$
(43)

where \(t_\mathrm{{R}}=t-r=t-\vert \mathbf {x}\vert \), and the wall thickness \(\eta \) can be specified by means of (29); namely, in terms of the quantities for the bubble collision profiles, such as the false vacuum field \(\phi _\mathrm{{f}}\) (equivalent to \(S_\mathrm{{f}}/\sqrt{4\pi }\)), the potential difference between the two minima \(\epsilon ^{4}\) (equivalent to \(V_\mathrm{{f}}\)) and half the separation of the bubbles \(b\) [15].

Result 1: The numerical computations of (43) are presented in Fig. 4; with various false vacuum field values, \(S_\mathrm{{f} }=\sqrt{4\pi }\phi _\mathrm{{f}}=\) (1) \(0.1\), (2) \(0.2\), (3) \(0.3\), (4) \(0.4\), in accordance with the scalar field solutions as presented by Fig. 2. Due to Eqs. (29) and (43), the amplitude of our GWs \(_\mathrm{{Q}}h^\mathrm{{TT}}( t) \) scales as \(\phi _\mathrm{{ f}}^{4}\) if the other conditions, \(\epsilon ^{4}\) and \(b\) are kept the same. Thus, with \(S_\mathrm{{f}}=\) (1) \(0.1\), (2) \(0.2\), (3) \(0.3\), (4) \(0.4\), the amplitude scales as (1) \(1\), (2) \(2^{4}\), (3) \(3^{4}\), (4) \(4^{4}\). The frequency of the waves is modulating due to the non-linearity of the collision dynamics in the all four cases of \(S_\mathrm{{f}}\), (1)–(4). However, the modulating frequency increases overall as \(S_\mathrm{{f}}\) increases, which is analogous to the tendency exhibited by \(S( u,v) =\sqrt{4\pi }\phi ( u,v) \) as shown in Fig. 2. In Fig. 4, we present \(_\mathrm{{Q}}h^\mathrm{{TT} }( t) r/( 2\pi \eta ^{2}) \) instead of \(_\mathrm{{Q} }h^\mathrm{{TT}}( t) \), and thus all the waveforms are plotted in the same scale. One should note here that our actual numerical data of the energy-momentum tensors \(T_{ab}\) for Eq. (43) have been obtained via Eq. (3) after solving Eqs. (2) and (4) simultaneously [18]. Therefore, our \(T_{ab}\) contain the full physical information about the bubble collisions in terms of the scalar field \(S=\sqrt{4\pi }\phi \), the geometry \(g_{ab}\) and the potential \(V( S)\); with the radiation reaction effects included in \(S\) and \(g_{ab}\).Footnote 4

3.2 A simplified method to compute gravitational waves in the quadrupole approximation

Reference [19] presents a simplified method to compute the GWs \(_{ \mathrm {Q}}h_{ij}^\mathrm{{TT}}(t,\mathbf {x})\) of Eq. (13) by neglecting the gravitational effects on the bubbles: namely, \(g_{\mu \nu }\) in Eq. (3) is replaced by \(\eta _{\mu \nu }\), assuming that the bubbles are in flat spacetime. Then Eq. (14) can be simplified as

$$\begin{aligned} I_{kl}(t_\mathrm{{R}})= & {} \int \mathrm{d}^{3}x^{\prime }T^{kl}(t_\mathrm{{ R}},\mathbf {x}^{\prime })\nonumber \\= & {} \int \mathrm{d}^{3}x^{\prime }\partial _{k}\phi (t_{ \mathrm {R}},\mathbf {x}^{\prime })\partial _{l}\phi (t_\mathrm{{R}}, \mathbf {x}^{\prime }), \end{aligned}$$
(44)

where the energy-momentum tensors from Eq. (3) have been reduced; \( T_{ij}\rightarrow \partial _{i}\phi \partial _{j}\phi \) because the terms proportional to \(\delta _{ij}\) in \(T_{ij}\) makes no contribution to gravitational radiation (13) due to the property of Eq. (10); namely, \(\Lambda _{ij,kl}\delta _{ij}=0\) [19].

By Eqs. (18) and (44) we have

$$\begin{aligned} I_{kl}= & {} \delta _{kx}\delta _{lx}\int \mathrm{d}^{3}x^{\prime }\left[ \left( \frac{ \partial \phi }{\partial x^{\prime }}\right) ^{2}-\frac{1}{2}\left( \frac{ \partial \phi }{\partial y^{\prime }}\right) ^{2}-\frac{1}{2}\left( \frac{ \partial \phi }{\partial z^{\prime }}\right) ^{2}\right] \nonumber \\= & {} \delta _{kx}\delta _{lx}\int \mathrm{d}^{3}x^{\prime }\left[ \left( \frac{\partial \phi }{\partial x^{\prime }}\right) ^{2}-\frac{1}{2}\left( \frac{\partial \phi }{\partial \rho ^{\prime }}\right) ^{2}\right] . \end{aligned}$$
(45)

Using Eq. (23), we can modify

$$\begin{aligned} \left( \frac{\partial \phi }{\partial \rho }\right) ^{2}=\left( \frac{ \partial \phi }{\partial \tau }\right) ^{2}\left( \frac{\partial \tau }{ \partial \rho }\right) ^{2}=\frac{t_\mathrm{{R}}^{2}-\tau ^{2}}{\tau ^{2}} \left( \frac{\partial \phi }{\partial \tau }\right) ^{2}. \end{aligned}$$
(46)

Then in the same manner as described above by Eq. (28), the integral \( I_{kl}( t_\mathrm{{R}}) \) is computed out of the volume piece \( \mathcal {V}=\Delta x\vert \tau \Delta \tau \vert \Delta \theta = \frac{\pi }{4}\eta ^{2}\Delta x[ 1+\mathcal {O}( \eta ^{2}/t_{ \mathrm {R}}^{2}) ]\):

$$\begin{aligned} I_{kl}( t_\mathrm{{R}})= & {} \delta _{kx}\delta _{lx}\int _{\mathcal { V}}\mathrm{d}^{3}x^{\prime }\left[ \left( \frac{\partial \phi }{\partial x^{\prime }} \right) ^{2}-\frac{1}{2}\left( \frac{\partial \phi }{\partial \rho ^{\prime } }\right) ^{2}\right] \nonumber \\= & {} 2\pi \delta _{kx}\delta _{lx}\int _{t_\mathrm{{R}}-\eta ^{2}/( 8t_{ \mathrm {R}}) }^{t_\mathrm{{R}}}\mathrm{d}\tau ^{\prime }\tau ^{\prime }\int _{-x_\mathrm{{o}}}^{x_\mathrm{{o}}}\mathrm{d}x^{\prime }\nonumber \\&\times \left[ \left( \frac{ \partial \phi }{\partial x^{\prime }}\right) ^{2}-\frac{t_\mathrm{{R} }^{2}-\tau ^{\prime 2}}{2\tau ^{\prime 2}}\left( \frac{\partial \phi }{ \partial \tau ^{\prime }}\right) ^{2}\right] +\mathcal {O}\left( \frac{\eta ^{4}}{t_\mathrm{{R}}^{3}}\right) \nonumber \\= & {} \frac{\pi }{2}\delta _{kx}\delta _{lx}\eta ^{2}\int _{0}^{x_\mathrm{{o} }}\mathrm{d}x^{\prime }\left( \frac{\partial \phi }{\partial x^{\prime }}\right) _{\tau =t_\mathrm{{R}}}^{2}+\mathcal {O}\left( \frac{\eta ^{4}}{t_\mathrm{{R} }^{3}}\right) .\nonumber \\ \end{aligned}$$
(47)

If the wall thickness \(\eta \) can be taken sufficiently small in Eq. (47), by Eqs. (13) and (47) we can compute the bubble-collision-induced GWs in the quadrupole approximation as

$$\begin{aligned} {}_\mathrm{{Q}}h_{ij}^\mathrm{{TT}}( t,\mathbf {x}) \approx \frac{ 2\pi }{r}\Lambda _{ij,kl}( \mathbf {n}) \delta _{kx}\delta _{lx}\eta ^{2}\int _{0}^{x_\mathrm{{o}}}\mathrm{d}x^{\prime }\left( \frac{\partial \phi }{\partial x^{\prime }}\right) _{\tau =t_\mathrm{{R}}}^{2}. \end{aligned}$$
(48)

Substituting Eq. (42) into Eq. (48), we finally express

$$\begin{aligned} _\mathrm{{Q}}h_{ij}^\mathrm{{TT}}( t,\mathbf {x})\approx & {} \frac{ 2\pi }{r}\eta ^{2}\bigg [ \delta _{ix}\delta _{jx}-2\delta _{ix}n_{j}\cos \vartheta \nonumber \\&+\,\frac{1}{2}n_{i}n_{j}\left( 1+\cos ^{2}\vartheta \right) -\frac{1 }{2}\delta _{ij}\sin ^{2}\vartheta \bigg ] \nonumber \\&\times \int _{0}^{x_\mathrm{{o}}}\mathrm{d}x^{\prime }\left( \frac{\partial \phi }{ \partial x^{\prime }}\right) _{\tau =t_\mathrm{{R}}}^{2}, \end{aligned}$$
(49)

where \(t_\mathrm{{R}}=t-r=t-\vert \mathbf {x}\vert \), and \(\eta \) is specified by (29).

Result 2: Toward the end of the bubble collisions, \(\tau \gg 1\), the scalar field oscillates around the true vacuum state, i.e. \(\vert \phi \vert \ll 1\), being nearly monochromatic. Then we can approximate Eq. (8) as

$$\begin{aligned} \square \phi ( \tau ,x) \approx \frac{6( \beta +1) V_{ \mathrm {f}}}{\beta S_\mathrm{{f}}^{2}}\phi ( \tau ,x), \end{aligned}$$
(50)

where we have replaced the curved spacetime Laplacian \(\nabla ^{2}\) by the flat spacetime d’Alembertian \(\square \equiv -\partial ^{2}/\partial \tau ^{2}-( 2/\tau ) \partial /\partial \tau +\partial ^{2}/\partial x^{2}\), neglecting the gravitational effects on the bubbles to simplify the problem.Footnote 5 With the help of Refs. [22, 23], we obtain a solution for Eq. (50):

$$\begin{aligned} \phi ( \tau ,x) =\phi _\mathrm{{o}}\frac{J_{1}\left( \omega _{ \mathrm {t}}\sqrt{\tau ^{2}-x^{2}}\right) }{\omega _\mathrm{{t}}\sqrt{\tau ^{2}-x^{2}}}, \end{aligned}$$
(51)

where \(J_{n}\) denotes the Bessel function of the first kind, \(\omega _{ \mathrm {t}}\equiv \sqrt{6( \beta +1) V_\mathrm{{f}}/( \beta S_\mathrm{{f}}^{2}) }\) represents the ‘terminal’ frequency of the bubble collisions, and \(\phi _\mathrm{{o}}\) is the amplitude which is determined by the initial conditions of the field. Substituting Eq. (51) into Eq. (49), the GWs emitted from the bubble collisions in the final stage can be computed. Figure 5 shows the GWs \(_{ \mathrm {Q}}h^\mathrm{{TT}}( t) r/( 2\pi \eta ^{2}) \) computed with various false vacuum field values, \(S_\mathrm{{f}}=\sqrt{4\pi } \phi _\mathrm{{f}}=\) (1) \(0.1\), (2) \(0.2\), (3) \(0.3\), (4) \(0.4\). Corresponding to the field values are the frequencies, \(\omega _\mathrm{{t} }\simeq \) (1) \(0.8124\), (2) \(0.4062\), (3) \(0.2708\), (4) \(0.2031\) \(\sim \) (1) \(1\), (2) \(1/2\), (3) \(1/3\), (4) \(1/4\), as can be seen from Fig. 5. Also, the amplitude of \(_\mathrm{{Q}}h^\mathrm{{TT}}( t) r/( 2\pi \eta ^{2}) \) scales as (1) \(1\), (2) \(2\), (3) \(3\), (4) \(4\) due to the factor \(\omega _\mathrm{{t}}^{-1}\sim S_\mathrm{{f}}\), as can be seen from Fig. 5. Then the amplitude of \(_\mathrm{{Q}}h^\mathrm{{TT} }( t) \) should scale as \(\eta ^{2}\omega _\mathrm{{t}}^{-1}\sim S_\mathrm{{f}}^{5}\sim \) (1) \(1\), (2) \(2^{5}\), (3) \(3^{5}\), (4) \(4^{5}\).

Fig. 5
figure 5

Plots of \(_\mathrm{{Q}}h^\mathrm{{TT}}(t)r/(2 \pi \eta ^{2})\) computed via the expression (49), with \(\phi \) given by (51). Corresponding to the false vacuum field values, \(S_\mathrm{{f}}=\sqrt{4\pi } \phi _\mathrm{{f}}=\) (1) \(0.1\), (2) \(0.2\), (3) \(0.3\), (4) \(0.4\), are the frequencies, \(\omega _\mathrm{{t}}\simeq \) (1) \(0.8124\), (2) \( 0.4062\), (3) \(0.2708\), (4) \(0.2031\) \(\sim \) (1) \(1\), (2) \(1/2\), (3) \(1/3\), (4) \(1/4\). The amplitude of \(_\mathrm{{Q}}h^\mathrm{{TT} }(t)r/(2\pi \eta ^{2})\) scales as (1) \(1\), (2) \(2\), (3) \(3\), (4) \(4\) due to the factor \(\omega _\mathrm{{t} }^{-1}\sim S_\mathrm{{f}}\). Then the amplitude of \(_\mathrm{{Q}}h^\mathrm{{TT }}(t)\) should scale as \(\eta ^{2}\omega _\mathrm{{t} }^{-1}\sim S_\mathrm{{f}}^{5}\sim \) (1) \(1\), (2) \(2^{5}\), (3) \(3^{5}\), (4) \( 4^{5}\)

4 Conclusions

We have computed GWs emitted from collisions of two equal-sized bubbles in time domain. The waveforms have been obtained for:

  1. (i)

    the initial-to-intermediate stage of strong collisions, and

  2. (ii)

    the final stage of weak collisions, in full General Relativity and in the flat spacetime approximation, using numerical and analytical methods, respectively.

During (i), the waveforms show the non-linearity of the collisions, characterized by a modulating frequency and cusp-like bumps, whereas during (ii), the waveforms exhibit the linearity of the collisions, featured by a constant frequency and smooth oscillations, as can be checked from Figs. 4 and 5, respectively. Also, depending on the false vacuum field value \(\phi _\mathrm{{f}}\), the waveforms have different scales of frequency. During (i), the modulating frequency increases overall as the false vacuum field value \(\phi _\mathrm{{f}}\) increases, whereas during (ii), the frequency \(\omega _\mathrm{{t}}\) decreases as the false vacuum field value \(\phi _\mathrm{{f}}\) increases in an inversely proportional relationship, i.e. \(\omega _\mathrm{{t}}\sim \phi _\mathrm{{f} }^{-1}\). It is interesting to note that the relationship between the false vacuum field value and the frequency during (i) changes almost inversely during (ii). In addition, the false vacuum field value \(\phi _\mathrm{{f}}\) affects the amplitude of the waveforms. During (i), the amplitude scales as \( \eta ^{2}\sim \phi _\mathrm{{f}}^{4}\), whereas during (ii), the amplitude scales as \(\eta ^{2}\omega _\mathrm{{t}}^{-1}\sim \phi _\mathrm{{f}}^{5}\), where \(\eta \) is a bubble wall thickness.

One of the notable differences between the waveforms emitted during (i) and during (ii) is the sign, as can be seen from Figs. 4 and 5. This is due to the difference between Eqs. (43) and (49): the integral in (49) is always positive while its counterpart in (43) is not necessarily. This has to do with the composition of the integrands in the two expressions. The integrand in (43) consists of the energy-momentum tensors \(T_{ab}\) which have been obtained via Eq. (3) after solving Eqs. (2) and (4) simultaneously [18]: thus \(T_{ab}\) contain the full physical information of bubble collisions in terms of the scalar field \(\phi \), the geometry \(g_{ab}\) and the potential \(V(\phi )\); with the radiation reaction effects included in \(\phi \) and \(g_{ab}\). However, as explained in the beginning of Sect. 3.2, the integrand in (49) comes only from the first term, with the second and third terms being disregarded in (3) as the gravitational effects on the bubbles are assumed to be neglected, following Ref. [19]. This, combined with the thin-wall approximation, results in the integrand in (49) being positive, which leads to the integral being also positive. But this is not the case for the integral in (43) due to the minus signs appearing in (3) and in the integrand in (43).

Throughout the paper, we used the thin-wall and quadrupole approximations to simplify our computations. These approximations served our purpose well in that we were able to gain some qualitative insights into the time-domain gravitational waveforms emitted from bubble collisions. However, to obtain more physically reasonable waveforms, taking into account a generic thickness and relativistic motion of bubble wall, it will be inevitable to include in our computations the next-to-leading order corrections beyond each approximation. A huge amount of computation will be involved in this task, and we leave it for follow-up studies.