Gravitational waves from first-order phase transitions: Ultra-supercooled transitions and the fate of relativistic shocks

TL;DR

This paper investigates gravitational wave generation during ultra-supercooled first-order phase transitions with dominant latent heat, focusing on relativistic shock waves and fluid shell dynamics to understand GW signals.

Contribution

It introduces an effective theory for relativistic fluid shell propagation in strong phase transitions, enhancing understanding of GW production mechanisms.

Findings

Relativistic fluid shells can produce significant gravitational waves.

The duration of fluid configurations impacts GW signal strength.

Shock wave behavior is crucial in ultra-supercooled transitions.

Abstract

We study the gravitational wave (GW) production in extremely strong first order phase transitions where the latent heat density dominates the plasma energy density, $\alpha \gtrsim 1$. In such transitions, bubbles develop extremely thin and relativistic fluid configurations, resulting in strong shock waves after collisions. We first propose a strategy to understand the GW production in such a system by separating the problem into the propagation part and the collision part. Focusing on the former, we next develop an effective theory for the propagation of the relativistic fluid shells. Using this effective theory, we finally calculate the expected duration of the relativistic fluid configurations and discuss its implications to the GW production.

Equations

$$\alpha$$

$$\boxed{\begin{array}{c}\mathrm{Extremely}\mathrm{thin}\mathrm{and}\mathrm{relativistic}\\ \mathrm{fluid}\mathrm{profiles}\end{array}}$$

$${\partial }_{\mu }{T}^{\mu \nu }$$

$$T_{\mu \nu }^{\mathrm{plasma}}$$

$$w_{+}{v}_{+}^2{\gamma }_{+}^2+p_{+}$$

$$p_{+}=\frac{1}{3}{a}_{+}{T}_{+}^4-ϵ,$$

$$p_{-}=\frac{1}{3}{a}_{-}{T}_{-}^4,$$

$$v_{+}{v}_{-}$$

$${\alpha }_{+}$$

$$(\xi -v)\frac{{\partial }_{\xi }\rho }{w}$$

$$(1-v\xi )\frac{{\partial }_{\xi }p}{w}$$

$${\partial }_{\xi }v$$

$$\mu (\xi ,v)$$

$$T_{\mu \nu }$$

$$p$$

$$c_s^2$$

$${\partial }_{t}u+{\partial }_{r}f+g$$

$$u$$

$${\partial }_t\left(\begin{array}{c}\rho \\ v\end{array}\right)+A{\partial }_r\left(\begin{array}{c}\rho \\ v\end{array}\right)+h$$

$$A$$

$$v_{C_{\pm }}$$

$$\pm \frac{c_s}{\rho +p}\left(\frac{\partial }{\partial t}+v_{C_{\pm }}\frac{\partial }{\partial r}\right)\rho +\frac{1}{1-v^2}\left(\frac{\partial }{\partial t}+v_{C_{\pm }}\frac{\partial }{\partial r}\right)v\pm \frac{v{c}_s}{1\pm v{c}_s}\frac{d-1}{r}$$

$$R_{\pm }$$

$${\frac{d}{dt}|}_{\pm }$$

$$R_{\pm }$$

$${\gamma }_s^2(t)$$

$${\rho }_{\mathrm{peak}}(t)$$

$$\ln {\rho }^{\prime }(t)$$

$$\ln {\gamma }^{2^{\prime }}(t)$$

$$p_{\mathrm{peak}}=\frac{p_0+{\rho }_0{v}_{\mathrm{peak}}{v}_s}{1-v_{\mathrm{peak}}{v}_s},v_s=\frac{(p_{\mathrm{peak}}+{\rho }_{\mathrm{peak}})v_{\mathrm{peak}}}{p_{\mathrm{peak}}{v}_{\mathrm{peak}}^2+{\rho }_{\mathrm{peak}}-{\rho }_0(1-v_{\mathrm{peak}}^2)}.$$

$${\gamma }_s^2$$

$$\frac{{\rho }_{\mathrm{peak}}}{{\rho }_0}$$

Full text

1. 1 Introduction

2. 2 Brief review on hydrodynamics around the bubble

3. 2.1 Basic equations

4. 2.1.1 Energy-momentum conservation

5. 2.1.2 Boundary condition at the scalar wall

6. 2.1.3 Fluid equations

7. 2.2 Detonation, deflagration and hybrid

8. 2.2.1 Detonation

9. 2.2.2 Deflagration

10. 2.2.3 Hybrid

11. 3 Problem and basic strategy

12. 4 Effective theory for relativistic shock propagation

13. 4.1 Assumptions and approximations

14. 4.2 Equations of motion and Riemann invariants

15. 4.2.1 Equations of motion

16. 4.2.2 Riemann invariants

17. 4.3 Effective theory for relativistic shock propagation

18. 4.3.1 Fundamental variables

19. 4.3.2 Equations

20. 4.3.3 Solution

21. 4.3.4 Initial conditions

22. 4.4 Comparison with numerical simulations

23. 5 Implications to gravitational wave production

24. 6 Discussion and conclusions

25. A Effective theory without the relativistic limit

26. B Details on numerical simulation

1 Introduction

In the upcoming decades, gravitational waves will prove to be an invaluable tool in probing the early universe. They will provide us opportunities to test high-energy physics theories such as inflation [1], cosmological first-order phase transitions [2, 3], preheating [4], topological defects [5], and so on. To observe gravitational waves, the LIGO Scientific Collaboration [6] and Virgo Collaboration [7] are now in operation, and KAGRA [8] is also expected to join soon. In addition, space-based detectors such as LISA [9] (and its possible successor BBO [10]) and DECIGO [11] have been proposed.

Among various sources of gravitational waves, cosmological first-order phase transitions are of our interest in this paper. Although both the electroweak and QCD phase transitions are known to be crossovers in the Standard Model (SM) [12, 13, 14], first-order phase transitions often occur in a wide range of extensions of the SM, and may play a crucial role in explaining the baryon asymmetry of the universe [15, 16, 17, 18]. During a first-order phase transition, bubbles of the broken phase nucleate, expand, and eventually collide with each other (Fig. 2).♢♢\diamondsuit1♢♢\diamondsuit1$\diamondsuit$1 See e.g. Refs. [19, 20, 21, 22] for pioneering works on the GW production in such systems.

The properties of the transitions are characterised by several quantities, one of which is the ratio between the latent heat and the plasma energy density just before the transition:

[TABLE]

For small $\alpha$ ($\ll 1$), the dynamics of the bubbles and the resulting GW production are relatively well-known (though still far from complete): the released latent heat turns into heat and the kinetic energy of the surrounding plasma, which in turn propagates in the form of sound waves characterised by the linearity of the fluid equation $(\partial_{t}^{2}-c_{s}^{2}\nabla^{2})\vec{v}=0$. These sound waves propagate even after bubble collisions and damping of the scalar configurations, and can be a long-lasting source of gravitational waves [23, 24, 25, 26]. At late times, these sound waves turn into turbulence which produce another characteristic form of the GW spectrum [22, 27, 28, 29, 30, 31] (see also reviews e.g. Refs. [32, 33, 34]).

In contrast, the dynamics of the transitions with large $\alpha$ ($\gtrsim 1$) – what we call the ultra-supercooled transitions in this paper – is rather unknown, despite that they are quite interesting from an observational viewpoint.♢♢\diamondsuit2♢♢\diamondsuit2$\diamondsuit$2 See Refs. [35, 36] for a recent discussion on the maximal possible strength of the transitions in specific models.

Indeed, first-order phase transitions in a certain class of models are known to have this property [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55] and we, at least naively, expect a large amount of gravitational wave to be produced because of the huge latent heat released in the universe. Previously, it had been thought that the scalar walls start to run away with sufficiently large $\alpha$ [56], but recently it has been pointed out that particle-splitting processes prohibit runaway solutions at least when the scalar field is gauged [57].♢♢\diamondsuit3♢♢\diamondsuit3$\diamondsuit$3 Nevertheless, the scalar field can play an important role in the GW production in extremely strong transitions. In such cases, the resulting GW spectrum has been estimated by the so-called envelope approximation [19, 58, 59, 60], but the effects beyond the envelope approximation can also be important [61, 62, 63, 64].

This means that an extremely strong detonation is realized in such a large $\alpha$ transition, and therefore the initial fluid profile after bubble collision is quite different from that of sound waves (characterised by linearity $(\partial_{t}^{2}-c_{s}^{2}\nabla^{2})\vec{v}=0$) or turbulence. Even if the system were to develop into sound waves or turbulence at late times, we must first understand its evolution just after bubble collision to know the resulting GW spectrum from this type of transition:

[TABLE]

Indeed, the spectrum of the gravitational waves generated from thin and long-lived bubbles can be quite different from that from the overlapping sound shells [26, 60]. Therefore, in order to correctly predict the GW spectrum from ultra-supercooled transitions, we need to understand the evolution of thin and relativistic fluid profiles. One may consider numerical simulations, but difficulties arise from (see also the profiles in Fig. 2 or 6):

• •

Strong energy and momentum concentration in the fluid and the resulting huge hierarchy between the thickness of the profile and the size of the whole system (typically $\gtrsim{\mathcal{O}}(10)$ $\times$ (bubble size at collisions))

• •

Strong shock waves forming around the propagating fluid

Given this, in this paper we tackle the ultra-supercooling regime and the resulting GW production with an analytic approach, in line with the present authors’ investigation so far [59, 60, 62].

The paper is organized as follows. In Sec. 2, we review the hydrodynamics around the expanding bubble walls before collisions. In Sec. 3, we make clear the problems we tackle and give an overview of the strategy to deal with them. In Sec. 4, we develop an effective theory for the relativistic shock propagation and check its validity by comparing it with the result of numerical simulations. In Sec. 5, we discuss the implications of our results to the GW production in large $\alpha$ transitions. In Sec. 6, we summarize. In Appendix A, we derive the equations without the relativistic limit, while in Appendix B we explain the details of numerical simulations.

2 Brief review on hydrodynamics around the bubble

In this section, we briefly review the fluid dynamics around an expanding bubble. This section closely follows Ref. [65]. See this reference for full derivations and more detailed explanations.

2.1 Basic equations

2.1.1 Energy-momentum conservation

All the necessary equations for the fluid profile before collision are obtained from the energy-momentum conservation

[TABLE]

The energy-momentum tensor $T_{\mu\nu}$ has two contributions, one from the scalar field $\phi$ driving the transition $T^{\phi}_{\mu\nu}$ and the other from the plasma particles $T^{\rm plasma}_{\mu\nu}$. Since we are interested in the dynamics of the cosmological-size bubbles, the plasma particles can be well-approximated by the fluid description. Assuming the perfect fluid form, we may write the corresponding part as

[TABLE]

where $w=\rho+p$ is the enthalpy density with $\rho$ and $p$ the energy density and pressure, respectively, $u^{\mu}=\gamma(1,\vec{v})$ is the fluid four-velocity with $\gamma=1/\sqrt{1-|\vec{v}|^{2}}$ the gamma factor of the fluid, and $g_{\mu\nu}={\rm diag}(+,-,-,-)$ is the metric tensor.

2.1.2 Boundary condition at the scalar wall

The bubbles that produce gravitational waves are of cosmological scale. In contrast, the scalar field dynamics is restricted within the “particle-physics scales” (say, $\sim$ TeV*-1* size for a TeV scale transition), a configuration that can be regarded as an infinitesimal boundary separating the symmetric phase $+$ and the broken phase $-$. The energy-momentum conservation across the boundary gives the following two conditions:

[TABLE]

where $v_{\pm}$ are the incoming and outgoing fluid velocities measured from the wall, respectively.

Using the relativistic gas approximation for simplicity, the pressure and energy density can be written in terms of temperature $T$ and vacuum energy $\epsilon$ of the symmetric phase relative to that of the broken phase:

[TABLE]

where $a_{\pm}$ are the number of light degrees of freedom in the symmetric and broken phase, respectively. Eqs. (2.3)–(2.5) give us two relations between $v_{\pm}$ and two dimensionless quantities from three dimensionful quantities, $T_{\pm}$ and $\epsilon$:

[TABLE]

where

[TABLE]

Note that the first equation gives $v_{+}v_{-}=1/3$ for $\alpha_{+}=0$. This will be important when we discuss the position of the shock front later.

2.1.3 Fluid equations

Except at the wall position, the fluid dynamics is governed by the energy-momentum conservation involving the fluid only. Since there is no characteristic length scale in the present system, the conservation equation $\partial^{\mu}T^{\rm fluid}_{\mu\nu}=0$ can be written in terms of the velocity of a coordinate point in the wave profile, $\xi\equiv r/t$, where $r$ is the radial coordinate and $t$ is the time elapsed since the nucleation of the bubble. The point at the wall position reaches the terminal velocity $\xi=\xi_{\rm wall}$. $v(\xi)$ is then the fluid velocity at $\xi$ as measured from the bubble center. After projecting onto the directions along and perpendicular to the flow of the fluid, we obtain the following equations

[TABLE]

The two derivatives $\partial_{\xi}\rho$ and $\partial_{\xi}p$ are related through the speed of sound in the plasma $c_{s}^{2}=(dp/dT)/(d\rho/dT)$♢♢\diamondsuit4♢♢\diamondsuit4$\diamondsuit$4In this paper, we use $c_{s}^{2}=1/3$ in both the symmetric and broken phases assuming that all the qualitative features of the system do not depend significantly on the possible deviations from this value., and hence we get the following equation describing the velocity profile:

[TABLE]

where $\mu$ is the Lorentz-transformed fluid velocity measured from $\xi$:

[TABLE]

Solving Eq. (2.10) for $\xi$ as a function of $v$ gives the trajectories (gray lines) shown in the left panel of Fig. 1. The boundary conditions at the wall positions are denoted by stars. The $v=\xi$ line separates the blue and white regions, while the red region is enclosed by $\mu\xi=c_{s}^{2}$ (the upper boundary) and $\mu=c_{s}$ (the lower boundary). Note that the derivative $\partial_{\xi}v$ diverges at the lower boundary along the gray lines. As seen from this figure, this system allows several types of solutions which we turn to next.

2.2 Detonation, deflagration and hybrid

There are three types of solutions for the motion of plasma in the present system: detonation, deflagration, and hybrid.

2.2.1 Detonation

First, let us take the boundary condition (i.e. the wall position) to be somewhere below the red region in Fig. 1. An example profile is outlined in blue in Fig. 1. In this case, the differential equation (2.10) can be easily solved to give solutions of detonation type. In detonations, the supersonic wall precedes the fluid in motion, and a nonzero fluid velocity arises from $c_{s}<\xi<\xi_{\rm wall}$ where $\xi_{\rm wall}$ is the wall velocity as mentioned previously.

2.2.2 Deflagration

Next, let us take the boundary condition in the white region satisfying $\xi<c_{s}$. An example is shown in yellow in Fig. 1, which has the fluid moving in front of the wall. These solutions are of deflagration type.

The deflagration profiles have a shock front. To see this, first note that all the gray lines corresponding to deflagrations cross the $\mu=c_{s}$ line (the lower boundary between the red and white regions) where the derivative $\partial_{\xi}v$ diverges. As we follow the lines further down to the attractor point $(\xi,v)=(c_{s},0)$, the profiles become double-valued. Since this is unphysical, we need to devise a method to have the fluid velocities jump to zero. The solution is to consider a shock front developing in the front end of such profile. The location of the shock is given by a consideration similar to Eq. (2.6): we can redefine the quantities with the subscripts $\pm$ to describe the fluid in front of and behind the shock front, respectively, with $\alpha_{+}$ taken to be [math] since there is no energy injection at the shock position. This gives us $v_{+}v_{-}=1/3$, which translates into $\mu\xi=1/3$ (the upper boundary between the red and white regions) in the plasma frame. Therefore, the fluid velocity is nonzero behind the shock front, i.e. the crossing point of the gray line and the boundary $\mu\xi=1/3$, while it drops to zero in front of the shock front.

2.2.3 Hybrid

Finally, let us take the boundary condition in the white region above the red region satisfying $\xi>c_{s}$. An example is shown in red in Fig. 1. These solutions are called the hybrid type because a nonzero fluid motion arises both in front of and behind the wall.

The fluid profile in front of the wall is similar to that of deflagrations: the fluid velocity takes its maximum value at the wall position and vanishes at the shock front position, $\mu\xi=1/3$. On the other hand, the profile behind the wall is clearly different from that of deflagrations: a nonzero fluid motion also arises behind the wall, a rarefaction wave of Jouguet type, starting from the line $\mu=c_{s}$. In fact, the single-valuedness of the fluid profile requires that the profile to be on or below the $\mu=c_{s}$ line, but it must also have its endpoint on or above the $\mu=c_{s}$ line lest the fluid velocity measured from the wall exceed the speed of sound, presenting another discontinuity.

3 Problem and basic strategy

Now let us focus on the fluid dynamics and the GW production in large $\alpha$ transitions. A sketch of the system is shown in Fig. 2. Just after bubble collisions, the relativistic fluid is pushed inside the regions of broken phase. To understand the evolution of such fluid, we first propose reducing the problem to two parts:

• •

Fluid propagation

• •

Fluid collision

each of which deforms the initial profile. We focus on the former in this paper and will report on the latter elsewhere. Note that even the first part is highly nontrivial because of the nonlinearity of the fluid equations of motion and the strong shock waves developing in the front end of the propagating fluid.

Let us take a closer look at the initial fluid profiles. As explained in Sec. 2, the left panel of Fig. 1 shows the fluid velocity around an expanding bubble. The horizontal axis is the distance measured from the bubble center normalized by the time elapsed from nucleation. The wall positions, denoted by stars, determine the location of the energy release and the strength of the fluid dragging. As explained previously, the profiles can be classified into three types depending on the wall position. From a particle physics point of view, the wall position is determined by the balance between the (1) pressure on the wall exerted by the energy release, and (2) friction on the wall arising from the coupling to the plasma particles.

The left panel shows that the walls with larger energy release $\alpha$ are located closer to the top-right corner. Since we are concentrating on the large $\alpha$ cases, the relevant combustion modes are detonations and hybrids. In the following we mainly work with detonations to make the analysis simple. We discuss hybrids in Secs. 5 and 6.

Now let us state the problems we tackle in this paper more clearly. Our initial relativistic fluid profile is a strong detonation. Assuming that the first fluid collision – which is simultaneous with the wall collision – does not change the profile significantly, we expect that this false-vacuum fluid would be pushed inside another expanding bubble, or the true vacuum. Fig. 3 is an illustration of the setup. Starting from this initial condition, we ask the following questions:

• •

How does the relativistic fluid evolve in time?

• •

What is the implication to the GW production?

We tackle the former in Sec. 4, while we discuss the latter in Sec. 5 based on our findings in Sec. 4.

As briefly mentioned above, shock waves are one of the main sources of complication in the present system. Generally a well-known phenomenon, they form when fluid is accelerated faster than the local speed of sound. However, calculating their time evolution is difficult because the shock front is a discontinuity at which derivatives diverge. While numerical methods are available in the literature, strong shock waves created in the phase transitions in the early universe are far beyond numerical tolerance. Therefore, we develop an effective theory for the shock front in Sec. 4 that allows us to describe its time evolution in the relativistic limit. We then compare its result with that of numerical simulations in an intermediate regime (see Fig. 4). Lastly, using the effective theory, we estimate how long such shock waves can last to source gravitational waves as thin, spherical objects before turning into sound waves or turbulence.

4 Effective theory for relativistic shock propagation

4.1 Assumptions and approximations

In this section, we model the system with a few quantities and derive the governing equations for propagation. As mentioned in Sec. 2, we consider the perfect-fluid form of the energy-momentum tensor:

[TABLE]

Note that we have omitted and will omit the label “plasma” in the following because the system consists of fluid only after bubble collision. We consider the equation of state

[TABLE]

to simplify the analysis. This gives the speed of sound

[TABLE]

The property of the fluid can change depending on the microphysics, but we expect that it does not change our qualitative results significantly.

Before moving on, we summarize our notations for subscripts here:

• •

“wall”: quantities pertaining to the bubble wall. For example, $\gamma_{\rm wall}$ is the $\gamma$ factor of the bubble wall.

• •

“max”: fluid quantities at the wall position during bubble expansion. For example, $\rho_{\rm max}$ and $\gamma_{\rm max}$ are the fluid energy density and $\gamma$ factor at the wall position during bubble expansion.

• •

“peak”: fluid quantities at the shock front after bubble collision. For example, $\rho_{\rm peak}$ and $\gamma_{\rm peak}$ are the fluid energy density and $\gamma$ factor at the shock front after bubbles have collided and shock profiles are pushed into the true vacuum fluid.

4.2 Equations of motion and Riemann invariants

The evolution equations for the system simply come from the conservation of the energy-momentum tensor. However, due to the discontinuities in the fluid profile, an ambiguity arises when we write down the evolution equations in differential forms. As we will see later, we have to impose the Rankine-Hugoniot relations to keep the conserved quantities intact across the discontinuities.

4.2.1 Equations of motion

Let us start from $\partial_{\mu}T^{\mu\nu}=0$. In the following, we assume $d$-dimensional spherical symmetry (e.g. $d=1$ is planar, $d=2$ is cylindrical, and $d=3$ is spherical). Our final interest is the spherical case, but we keep $d$ as a free parameter since we will check the validity of our theory using numerical simulation for the $d=1$ case later. The evolution equations reduce to (see e.g. Ref. [66])

[TABLE]

where

[TABLE]

In terms of $\rho$ and $v$, we have

[TABLE]

where $A$ and $h$ are

[TABLE]

4.2.2 Riemann invariants

It is known that two conserved quantities exist along the eigenvalues

[TABLE]

of the matrix $A$. Indeed, Eq. (4.6) can be rewritten as

[TABLE]

This implies that, for $d=1$, the quantities called the Riemann invariants

[TABLE]

are conserved along the directions defined by

[TABLE]

Note that these directions are paths of sound wave propagations since $v_{C_{\pm}}=(v\pm c_{s})/(1\pm vc_{s})$ are just the relativistic sums of $v$ and $c_{s}$. They are called the forward and backward characteristics, respectively, for which we have used the label $C_{\pm}$. For the equation of state (4.2), the conserved quantities become

[TABLE]

where we took the relativistic limit in the last equality.

4.3 Effective theory for relativistic shock propagation

4.3.1 Fundamental variables

In this section we derive the effective theory for the shock wave propagation. However, before doing so, let us talk about result of a numerical simulation for $d=1$ with $(\gamma_{\rm wall},\alpha)=(10,10)$ (Fig. 6). Note that the setup in this figure is not of our final interest because the initial profile is obtained from an expanding spherical bubble while the evolution is calculated with the assumption of planar symmetry. The aim of this figure is to get a rough idea of the shock propagation after bubble collision and to check the validity of the effective theory. The length scale of the horizontal axis is chosen so that the collision time $t_{c}$ is unity (with the nucleation time set to be zero as mentioned above). We observe that

• •

The initial position of the fluid front is $r/t_{c}\simeq 1$ because the wall is highly relativistic.

• •

Soon after the fluid is pushed into the low-energy plasma, the peak values of $\rho$ and $\gamma^{2}$ get rearranged to certain values. After that, the peak values evolve gradually.

• •

At the same time, a discontinuity detaches from the peak and evolves toward the rear side of the fluid profile.

As far as the GW production is concerned, we are not much interested in the tail part of the fluid profile. Rather, we focus on the time evolution of the high-energy peak where the energy-momentum tensor takes large values. Therefore, we model the strong shock front using the following five variables:

• •

$\gamma$-factor of the shock wave

[TABLE]

• •

Peak values of the fluid energy density and their $\gamma$-factors

[TABLE]

• •

Derivatives of the fluid energy density and the $\gamma$-factor at the shock front

[TABLE]

where $r_{s}(t)$ is the position of the shock front. In the following, we derive the governing equations for these variables.

4.3.2 Equations

Rankine-Hugoniot relations.

As discussed in Sec. 4.2, equations of motion in differential forms pose an ambiguity due to the discontinuity in the fluid profile. Therefore, we resort to the Rankine-Hugoniot conditions which guarantee the conservation law on both sides of the shock. In the present system, we have two Rankine-Hugoniot conditions arising from two-component equations of motion. They can be read off as (see e.g. Ref. [66]),

[TABLE]

Assuming the equation of state (4.2), these reduce to

[TABLE]

They further get simplified in the relativistic limit:

[TABLE]

Note that these are constraint equations.

Time evolution equations.

In order to derive the time evolution equations, we use the time evolution of the Riemann invariants (Eq. (4.9)) along the two characteristics $C_{\pm}$. Fig. 5 is a zoom-in around the shock front with the red and blue lines describing the direction of the characteristics. The rough idea is to relate the time evolution of the peak values (A $\to$ B) with the spatial derivative of the profile (C $\to$ A or D $\to$ B) by connecting B and C or A and D with the characteristics. We take the relativistic limit in the following.

Let us first use $C_{+}$. By considering A, B, and C, we have

[TABLE]

where we used the time evolution of $R_{+}$ (Eq. (4.9)) from the second to third line. Since $C_{+}$, $C_{-}$, fluid motion and shock front are all extremely inclined to the forward direction in the relativistic limit (see Fig. 5), we also used $v_{\rm peak}\simeq 1$. Furthermore, we replaced $r_{s}$ with $t$ in the last line by taking the relativistic limit. The relation between $\Delta t$ and $\Delta r_{\rm AC}$ is obtained from Eq. (4.20) and

[TABLE]

as

[TABLE]

Therefore, we get

[TABLE]

The other evolution is obtained from the consideration of A, B, and D with A and D connected by $C_{-}$:

[TABLE]

where we used the time evolution of $R_{-}$ (Eq. (4.9)) from the second to third line and $v_{\rm peak}\simeq 1$ from the third to fourth line. The relation between $\Delta t$ and $\Delta r_{\rm BD}$ is calculated from Eq. (4.20) and

[TABLE]

as

[TABLE]

The result is

[TABLE]

Here, we notice one interesting thing. Since $\rho_{\rm peak}$ is related to $\gamma_{\rm peak}^{2}$ via Eq. (4.20), $\ln\rho^{\prime}$ must be related to $\ln\gamma^{2^{\prime}}$ to make Eqs. (4.24) and (4.28) consistent. After some calculations, we find that

[TABLE]

Substituting this into either Eq. (4.24) or (4.28) and using Eq. (4.20), we obtain

[TABLE]

Energy and momentum domination by shock front.

So far, we have derived only four constraint or evolution equations for five variables. The last condition we add is the energy and momentum dominance by the shock front. In order to do so, we need a specific ansatz for the fluid profile. In our analysis, we adopt

[TABLE]

Note that these profiles are only estimates and the choice is rather arbitrary. However, as long as the energy and momentum of the system are dominated by the shocks, we expect that the relation “(energy or momentum at the peak) $\times$ (thickness of the fluid) $\simeq$ const.” holds well. Indeed, in the relativistic limit, the energy or momentum calculated from the above ansatz (for $d=2$ and $3$, energy or momentum per unit angle and solid angle, respectively) becomes

[TABLE]

and satisfies (energy or momentum at the peak) $=\rho_{\rm peak}\gamma_{\rm peak}^{2}$ and (thickness of the fluid) $=1/(\ln\rho^{\prime}+\ln\gamma^{2^{\prime}})$. Note that $\sigma$ has different mass dimensions depending on the value of $d$. We use Eq. (4.32) with $\sigma$ calculated from the initial profile as the last ingredient of our effective theory. Though this constraint is rather qualitative in contrast to the other four equations, our effective theory (4.20), (4.24), (4.28), and (4.32) describes the duration of the shock rather correctly as we will see later.

4.3.3 Solution

The five equations (4.20), (4.24), (4.28), and (4.32) can be solved analytically. For $d=1$, the solution becomes

[TABLE]

[TABLE]

[TABLE]

with the following constraints among the initial values

[TABLE]

On the other hand, for $d=3$, the solution becomes

[TABLE]

[TABLE]

with the terms with the exponent $\delta\equiv 10/13$ becoming numerically negligible at late times compared to other terms. The constraints among the initial values are given by

[TABLE]

This solution implies that the lifetime of the relativistic fluid profile is qualitatively given by

[TABLE]

as long as we take only propagation effects into account. Note that the thickness of the fluid ($\sim 1/\ln\rho^{\prime}$ or $1/\ln{\gamma^{2}}^{\prime}$) is also determined by this timescale as

[TABLE]

In Sec. 4.3.4, we discuss how to determine the initial conditions (4.36) and (4.41) just after bubble collision from the profile of expanding bubbles.

4.3.4 Initial conditions

So far, we have discussed the effective theory without mentioning the initial conditions. For detonations, which we have used in this paper, the initial values $\rho_{\rm peak}$ or $\gamma_{\rm peak}^{2}$ can be estimated by the Riemann invariant $R_{+}$ along the characteristic $C_{+}$ as shown in the left panel of Fig. 5.♢♢\diamondsuit5♢♢\diamondsuit5$\diamondsuit$5 For hybrids, we cannot use the estimate (4.44) because the shock does not start from the wall position $W$ in the left panel of Fig. 5. However, in Sec. 5, we slightly extrapolate Eq. (4.44) to estimate the GW production for the hybrid regime as well.

In the relativistic limit, the relation is given by

[TABLE]

With the Rankine-Hugoniot relation (4.20), we have

[TABLE]

In Sec. 4.4, we indeed see that this relation gives a good estimate on the initial values.

4.4 Comparison with numerical simulations

Let us check the validity of our effective theory using numerical simulations with moderate relativisticity. Adopting the numerical scheme proposed in Ref. [67], we perform a numerical simulation for $(\gamma_{\rm wall},\alpha)=(10,10)$, values with which taking the relativistic limit is valid (see Fig. 4). The result is shown in Fig. 6. The initial profiles, denoted by thick blue lines, are calculated per Ref. [65], and we evolved the system from $t=t_{c}$ to $t=2t_{c}$. We observe several features explained in Sec. 4.3.1. From this fluid evolution, we extracted the five quantities introduced in Sec. 4.3.1 and plotted them in Fig. 7. Though the numerical scheme in Ref. [67] tracks the shock evolution quite well, numerical viscosity still plays a role. Therefore, we obtain the peak quantities by extrapolation as explained in Appendix. B. In Fig. 7, the blue points are the quantities obtained by this extrapolation, while the red lines are the prediction from our effective theory. For $t_{c}<t<1.1\times t_{c}$, the blue points are not shown because the peak structure is not clear during that time since the substructures have not completely detached from the peak (e.g. the top-left panel of Fig. 6). We see that our theory predicts the shock evolution qualitatively well, while it does not coincide with the data perfectly. There are two possible reasons for this deviation:

• •

The Rankine-Hugoniot relations (4.20) and the time evolution equations (4.24) and (4.28) are derived in the relativistic limit, while our input parameters $(\gamma_{\rm wall},\alpha)=(10,10)$ are not as relativistic.

• •

The last ingredient of our effective theory (the energy and momentum domination by the shock front) is not a strict condition but rather a qualitative one.

Regarding the first point, we plot Fig. 9 to check the Rankine-Hugoniot relations in our simulation. In the left panel, the blue and red points show $\gamma_{s}^{2}$ and $2\gamma_{\rm peak}^{2}$, respectively, which should coincide with each other in the relativistic limit as in the first relation in Eq. (4.20). The yellow points are the Rankine-Hugoniot relations without the relativistic limit (Eq. (A.1) in Appendix A). We see that all the blue, red, and yellow points coincide well with each other. On the other hand, the right panel of Fig. 9 is to check the second relation in Eq. (4.20). The blue and red points are $\rho_{\rm peak}/\rho_{0}$ and $(16/3)\gamma_{\rm peak}^{2}$, respectively, while the yellow points are the corresponding quantities without the relativistic limit (Eq. (A.2) in Appendix A). We see that the blue and yellow points coincide well with each other. Also, Fig. 9 checks the relation (4.29). The blue and red points are $\ln\rho^{\prime}$ and $(7/6)\ln{\gamma^{2}}^{\prime}$, respectively. As it is clear from the derivation, this relation is also subject to the first source of error listed above. However, we see that the difference between the two quantities $\ln\rho^{\prime}$ and $(7/6)\ln{\gamma^{2}}^{\prime}$ is only about $10\%$.

To summarize this section, we found that our effective theory of shock propagation describes the evolution of the system rather correctly. In the next section, we use it to discuss the implications to the GW production.

5 Implications to gravitational wave production

Now we discuss the implications of our effective theory to the GW production in ultra-supercooled transitions. In this section, we present our results in the $\gamma_{\rm wall}$-$\alpha$ plane. In order to translate the $\gamma_{\rm wall}$-$\gamma$ plane (the right panel of Fig. 1) into this plane, we follow the procedure described in Ref. [65]. The result is the parameter regions shown in Figs. 11–13. In these figures, the region below the solid blue line is for detonations, while the one between the solid red and blue lines is for hybrids.♢♢\diamondsuit6♢♢\diamondsuit6$\diamondsuit$6 Note that the transition from a detonation to a hybrid is continuous in $\alpha$ in contrast the right panel of Fig. 1 where there is no consistent solution inside the red region.

Though the estimate on the initial condition (4.44) applies only to the detonation profiles (see footnote $\diamondsuit$5), we slightly extrapolate our results to the hybrid region, interpreting $\rho_{\rm max}$ and $\gamma_{\rm max}$ in Eq. (4.44) as the maximal fluid energy density and $\gamma$ factor in the hybrid case.

Now let us explain each figure. Fig. 11 is the contour plots for $\rho_{\rm max}/\rho_{\infty}$ (left) and $\gamma_{\rm max}^{2}$ (right). The discontinuities across the blue line (corresponding to the boundary between hybrids and detonations) are already known [65] and we successfully reproduce them. One important observation from this figure is that the contours are parallel to the $\gamma_{\rm wall}$-axis in the detonation region. This means that

• $\rho_{\rm max}$ and $\gamma_{\rm max}$, which are the input parameters of our effective theory through Eq. (4.44), are $\gamma_{\rm wall}$ independent.

This fact may be helpful when estimating the GW production from particle physics models since it is often difficult to estimate the wall velocity from microphysics. Actually, this was already expected from Fig. 1: the limit $\xi_{\rm wall}\to 1$ does not cause any divergence in the differential equation (2.10), and therefore the dynamics does not depend on $\xi_{\rm wall}$ as long as it is close to unity.♢♢\diamondsuit7♢♢\diamondsuit7$\diamondsuit$7 The authors thank T. Konstandin for pointing this out.

Fig. 11 is the contour plots for $\rho_{0}/\rho_{\infty}$ (left) and $4\pi\sigma/\frac{4\pi}{3}\rho_{\infty}t_{c}^{3}$ (right). From the left panel, we see that, for a fixed $\gamma_{\rm wall}$, the plasma energy density in the broken phase gets smaller as $\alpha$ gets larger. This can be understood from the behavior of friction: the larger the $\alpha$, the larger the pressure on the wall, and therefore the friction on the wall increases accordingly in order to keep $\gamma_{\rm wall}$ constant. The increased friction results in a stronger drag of the fluid behind the wall, thus the energy density inside the bubble gets smaller. On the other hand, the quantity shown in the right panel is almost the same as $\alpha$. This is a cross-check of our calculation: $4\pi\sigma$ is the total energy in the whole bubble at $t=t_{c}$ (note that we are interested in the relativistic limit and the difference between energy and momentum does not matter), while $\frac{4\pi}{3}\rho_{\infty}t_{c}^{3}$ is the plasma energy density in the symmetric phase multiplied by the volume of the bubble at the time of collision. Therefore, this quantity should coincide with $\alpha$, and it does indeed.

Fig. 12 is the contour plots for $\rho_{\rm peak}/\rho_{\infty}$ (left) and $\gamma_{\rm peak}^{2}$ (right) calculated from $\rho_{\rm max}$ and $\gamma_{\rm max}^{2}$ using Eq. (4.44) and Eq. (4.20). In these plots, we extrapolated our results to the hybrid region by interpreting $\rho_{\rm max}$ and $\gamma_{\rm max}^{2}$ in Eq. (4.44) as the maximal energy density and $\gamma$ factor squared of the fluid at the wall position. Note that, as explained in the beginning of this section, this procedure is not justified in the strict sense because the point $W$ in the left panel of Fig. 5 does not coincide with the starting point of the shock wave. Therefore, the contours in the hybrid region should be taken as reference values.

Finally, Fig. 13 is a contour plot for the lifetime of shock waves (4.42). This is the main result of this paper: the propagating fluid in the broken phase remains relativistic for quite a long time ($\sim 10^{1,2,3}\times t_{c}$). This means that the thickness of the shock wave remains much smaller than the typical bubble size until late times, and that the onset of the sound wave regime (characterized by the linearity of the fluid equation) is delayed. However, please note that this conclusion is subject to change once we include the effects of shock wave collisions.

6 Discussion and conclusions

In this paper, we studied the fluid dynamics and the resulting GW production in ultra-supercooled first-order phase transitions characterized by $\alpha\gtrsim 1$. Though such transitions are observationally and theoretically motivated, their evolution and GW production mechanisms are quite unclear. Such system is likely to proceed with a strong detonation profile [57] in which energy gets highly concentrated behind the wall (see Fig. 2). It is extremely hard to tackle this system with numerical simulations because of the huge hierarchy between the energy concentration scale and the simulation volume, and also of the strong shock waves forming around the propagating fluid. Especially, it is unclear when the system develops into the sound wave or turbulence regimes, the former of which is characterized by the linearity of the fluid equation and is widely believed to serve as a long-lasting source of gravitational waves [23, 24, 25, 26]. It is of great importance to clarify this because the GW spectrum of overlapping sound shells can be quite different from that of thin and long-lived bubbles [26, 60].

In order to study the system, we divided the problem into two parts: the fluid propagation and collision during bubble expansion (Sec. 3). The aim of this paper was to tackle the former. For this purpose, we developed an effective theory of shock wave propagation (Sec. 4) valid in the relativistic limit described by a few degrees of freedom (Sec. 4.3.1). The theory has two constraint equations (the Rankine-Hugoniot relations) and two time evolution equations (resulting from the two characteristics), and becomes closed by assuming the energy or momentum dominance by the shock front (Sec. 4.3.2). The effective theory rendered the system analytically solvable (Sec. 4.3.3), and we checked the validity of the theory by comparing it with numerical simulations in the regime of moderate relativisticity (Sec. 4.4).

Equipped with the theory, we finally discussed the implications to the GW production in ultra-supercooled transitions in Sec. 5. We first found that the initial conditions of the effective theory depend on $\alpha$, but not on $\gamma_{\rm wall}$. This fact can be helpful since it is often difficult to calculate the wall velocity starting from particle physics models. We also found that the lifetime of the fluid relativisticity can be orders of magnitude larger than the typical bubble size (Fig. 13). However, in estimating this lifetime, we took only the propagation effects into account. To further understand the GW production in large $\alpha$ transitions, we need to incorporate the effects of fluid collisions. This will be addressed in future work.

Acknowledgments

The authors are grateful to Thomas Konstandin and Géraldine Servant for helpful comments on the manuscript, and particularly grateful to Sangjun Lee for helpful discussions at the initial stage of this project. The work of RJ was supported by Grants-in-Aid for JSPS Overseas Research Fellow (No. 201960698). The work of RJ was supported by IBS under the project code, IBS-R018-D1. The work of RJ is supported by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306. The work of MT was supported by JSPS Research Fellowships for Young Scientists.

Appendix A Effective theory without the relativistic limit

In this appendix, we summarize the Rankine-Hugoniot conditions and time evolution equations without taking the relativistic limit. For the constraint (4.32), we do not show the corresponding equation since Eq. (4.32) is already an approximate relation.

Rankine-Hugoniot conditions

The Rankine-Hugoniot conditions are $d$-independent and given by

[TABLE]

Time evolution equations

The time evolution equations are derived in almost the same way as Sec. 4.3.2 except that we do not take the relativistic limit. For $C_{+}$ and $C_{-}$, we have

[TABLE]

[TABLE]

These give the time evolution equations without the relativistic limit as

[TABLE]

where

[TABLE]

Note that $r_{s}$ in Eqs. (A.5) and (A.6) is related to $v_{s}$ through $v_{s}=dr_{s}/dt$ and the initial condition at the bubble nucleation time given by $r_{s}(t=0)=0$.

Appendix B Details on numerical simulation

In this appendix, we explain how we read off the shock position and peak values from numerical data. Fig. 14 is a snapshot of the numerical simulation explained in Sec. 4.4 at $t=1.5t_{c}$. As seen from the figure, both quantities are smeared around the peak due to the numerical scheme we adopted and the finiteness of the number of grids. We first identified the position where the absolute value of the spatial derivative of $\rho$ or $\gamma^{2}$ becomes the largest as the true shock position (yellow points). The shock position thus estimated is nearly identical for $\rho$ or $\gamma^{2}$. We next identify the two points where the value of $\rho/\rho_{\infty}$ or $\gamma^{2}$ is smaller than the peak value by $1$ or $2$ (red points). Then we finally extrapolate these points to the estimated shock position to estimate the true peak values of $\rho$ and $\gamma^{2}$ as well as the derivatives $\ln\rho^{\prime}$ and $\ln{\gamma^{2}}^{\prime}$.