HydroGrav: Precise hydrodynamics and gravitational waves for cosmological phase transitions

In this work, we demonstrate our improved hydrodynamic calculation using the simplest extension of the Standard Model Higgs sector. Namely, we consider a first-order phase transition driven by a real scalar singlet coupled to the Higgs [19, 36], whose tree-level potential is

where the Higgs doublet and scalar singlet are

Here, $\phi_{h}$ and $\phi_{s}$ are background fields, $h$ and $s$ are the dynamic degrees of freedom of the Higgs and scalar fields, respectively, and $G^{\pm,0}$ are the Goldstone bosons. After symmetry-breaking, we are left with two singly charged $G^{\pm}$ and a neutral CP-odd $G^{0}$ Goldstone bosons, corresponding to the electroweak $W^{\pm}$ and $Z$ bosons. In addition, we have the CP-even scalars $h$ and $s$, with masses

where $h$ is the SM Higgs, and we have assumed the EWSB vacuum $(\phi_{h},\phi_{s})=(v_{h},0)$. In regions of the parameter space, the final vacuum is reached by a two-step transition

and we consider GWs generated by the second step of the transition.

The effective potential is constructed in the context of three-dimensional effective field theory (3dEFT) using the code DRalgo [26]²²2DRalgo utilises the package GroupMath [29] during model creation.. Specifically, we consider the effective field theory obtained by integrating out the heavy thermal scale, yielding a dimensionally reduced three-dimensional theory. Within this framework, the one-loop effective potential is given by

where the subscript ‘$3$’ indicates quantities in the dimensionally reduced theory. The tree-level potential $V_{3,0}$ is obtained by replacing all quantities in $V_{0}(\phi)$ with their 3D counterparts. Through matching to the parent UV theory, all 3D parameters become functions of both temperature and the parameters of the underlying 4D Lagrangian. The scalar mass eigenstates contributing to the sum in 3.9 are obtained using the same approach.

For a generic theory that undergoes a single first-order phase transition, the effective potential possesses two distinct minima: $\phi_{m}^{+}$, associated with the high-temperature symmetric phase, and $\phi_{m}^{-}$, corresponding to the low-temperature broken phase. Utilising equations (3.3) and (3.4), the equations of state in the symmetric ($s$) and broken ($b$) phases are expressed as

while the associated sound speed in each vacuum phase is given by

We denote the pressure and energy density on either side of the wall as $p_{\pm}$ and $e_{\pm}$, defined by

For hydrodynamic configurations that develop a shock front, the local thermodynamic quantities at the interface are given by

However, a standard approach in the literature approximates the exact equation of state (4.15), derived from the effective potential, via the so-called bag model to simplify the hydrodynamics of cosmological phase transitions [68, 47, 45, 63, 28, 62, 52, 54, 55, 46]. Within this framework, the thermal contribution to the effective potential is expressed as $(1/3)a(T)T^{4}$, where $a(T)=g_{\text{eff}}(T)\pi^{2}/30$ and $g_{\text{eff}}(T)$ denote the effective number of relativistic degrees of freedom at temperature $T$. This assumes a purely radiation-dominated equation of state $p=(1/3)e$. Consequently, the bag equation of state is characterised by a constant sound speed of $c_{s}^{2}=1/3$ in both vacuum phases, yielding [18]

where $\epsilon$ is the false-vacuum energy of the Higgs potential (often referred to as the bag constant), defined to vanish in the broken phase at $T=0$, and $a_{\pm}=g_{\pm}\pi^{2}/30$. The quantities $g_{\pm}$ are the effective degrees of freedom in the symmetric and broken phases, which satisfy $g_{+}>g_{-}$.

The phase transition strength is conventionally parametrised by

which defines the ratio of the latent vacuum energy released during the transition to the background radiation energy density. Here, the trace anomaly of the fluid’s stress-energy tensor is defined as

with $\Delta\theta=\theta_{+}-\theta_{-}$.

The bag equation of state is generally a robust approximation in the symmetric phase, where Standard Model particles remain effectively massless during an electroweak-scale transition and behave dynamically akin to radiation. However, this approximation typically breaks down in the broken phase, where thermal corrections induce deviations in the sound speed from $c_{s}^{2}=1/3$ [53, 69]. A phenomenological remedy involves modifying the equation of state (4.19) to accommodate different sound speeds in the respective vacuum phases [53, 32, 31, 77, 1]. To this end, we adopt the $\mu\nu$ (or improved bag) model introduced in [77]³³3The formulation of the $\mu\nu$ model described in [31, 1] is physically equivalent to (4.22) when defining the temperature-dependent coefficients in equation (2.15) of [31] as $a_{+}(T)=a_{+}T^{4-\mu}$ and $a_{-}(T)=a_{-}T^{4-\nu}$, with $\mu=1+1/c_{s,+}^{2}$ and $\nu=1+1/c_{s,-}^{2}$., which reads

where $c_{s,+}=c_{s,s}(T_{+})$ and $c_{s,-}=c_{s,b}(T_{-})$ denote the sound speeds in the symmetric and broken vacuum phases, respectively, and are treated as constants.

In this generalised framework, the trace anomaly takes the form $\bar{\theta}=e-p/c_{s,-}$, yielding the modified strength parameter

In the limit $c_{s,+}^{2}=c_{s,-}^{2}=1/3$, the bag model definition (4.20) is recovered.

By combining this strength parameter with the matching conditions (4.12) and the $\mu\nu$ equation of state (4.22), the fluid velocity immediately in front of the wall is

where $\alpha_{+}$ is the strength parameter evaluated at $T_{+}$, and the positive (negative) root is taken for $|v_{+}|>|v_{-}|$ ($|v_{+}|<|v_{-}|$). For hydrodynamic profiles exhibiting a shock front, the equation of state across the shock (4.18) takes the form

so the matching conditions (4.14) are

Using either of these equation of state approximations significantly simplifies the relativistic hydrodynamics required to compute the gravitational wave spectrum. This permits a straightforward mapping of a specific effective potential onto the fluid dynamics by extracting the relevant phase transition parameters and evaluating either (4.20) or (4.23). Nevertheless, this methodology remains intrinsically ‘model-independent’, as the parametrised equation of state relies purely on macroscopic properties rather than the microscopic structure of a particular particle physics model.

In the remainder of this work, we quantify the impact of these simplifying assumptions on both the fluid dynamics and the resulting gravitational wave spectrum. We achieve this by performing the hydrodynamic calculations utilising the exact, fully model-dependent equation of state derived directly from the effective potential. While partially model-dependent hydrodynamic treatments have been explored previously [72, 78, 75]⁴⁴4In [72], the exact equation of state was used within the broken phase, while the $\mu\nu$ equation of state was used in the symmetric phase. In [75], the authors extend the bag model with additional temperature-dependent terms in the broken phase that reflect the structure of the effective potential., our fully model-dependent approach operates independently of both the bag and $\mu\nu$ approximations.

Hydrodynamic solutions for a perfect fluid can be categorised by the direction of flow of the fluid across the bubble wall [51]. For deflagration solutions, the fluid flows outwards from the bubble wall, so the pressure decreases ($p_{+}>p_{-}$) and velocity increases ($|v_{+}|<|v_{-}|$). For detonation solutions, the fluid flows inwards from the bubble wall, so the pressure increases ($p_{+}<p_{-}$) and velocity decreases ($|v_{+}|>|v_{-}|$). Here, we use the ‘$+$’ subscript to denote quantities in front of the bubble wall and the ‘$-$’ subscript for quantities behind the wall. The fluid velocities $v_{\pm}$ are given in the bubble wall frame, and we use a tilde to denote the fluid velocities $\tilde{v}_{\pm}$ in the centre-of-bubble frame (also called the fluid frame).

These solutions can be further categorised as weak (fluid velocities on either side of the wall are both subsonic or both supersonic), strong (one of $v_{\pm}$ is subsonic and the other is supersonic), or Chapman-Jouguet ($v_{-}=c_{s,-}$ or $v_{+}=c_{s,+}$). The sixteen types of hydrodynamic solutions are shown in Table 1. Note that weak supersonic deflagrations and weak subsonic detonations belong to the class of ‘inverse’ phase transitions [8, 9], which we do not consider in this work. These involve the vacuum undergoing a phase transition from a low-temperature phase to a high-temperature phase, so the false vacuum energy is negative (i.e. $\alpha_{\bar{\theta}}<0$) and fluid flows outwards from the bubble in the centre-of-bubble frame.

For deflagration solutions, the fluid inside the bubble is at rest in the centre-of-bubble frame ($\tilde{v}_{-}=0$), so $v_{-}=\xi_{w}$. The bubble wall is therefore subsonic ($\xi_{w}<c_{s,-}$) for weak subsonic deflagrations, or equal to the speed of sound ($\xi_{w}=c_{s,-}$) for Chapman-Jouguet deflagrations. Strong deflagrations are known to be mechanically unstable [51] and would quickly decay into a rarefaction wave and a weak subsonic deflagration [49]. The subsonic expansion of the bubble compresses the fluid in front of the wall, forming a discontinuous shock front $\xi_{\text{sh}}>c_{s,+}$ that reheats the fluid behind it (i.e. $T_{+}>T_{N}$).

The condition $|v_{+}|<|v_{-}|$ for deflagrations gives both a lower and upper bound for the strength parameter $\alpha_{+}$ such that a shock can form. Inverting equation (4.24), we find that $\alpha_{+}^{\text{min}}=0$ for $v_{+}=v_{-}$, corresponding to the case in which no energy is transferred from the phase transition to the fluid, and $\alpha_{+}^{\text{max}}=1/3$ for $v_{+}=0$, which represents the maximum phase transition strength compatible with a subsonic wall. For $\alpha_{+}<\alpha_{+}^{\text{min}}$ or $\alpha_{+}>\alpha_{+}^{\text{max}}$, weak subsonic deflagration solutions are not possible.

For strongly supercooled transitions, detonation solutions are possible [35, 27]. Here, the fluid ahead of the bubble is undisturbed ($\tilde{v}_{+}=0$), so $v_{+}=\xi_{w}$. Since the pressure in front of the wall is greater than behind it, fluid flows into the bubble as a rarefaction wave until it smoothly comes to a halt at $\xi=c_{s,-}$. Given the convention $v(\xi)\geq 0$ for direct transitions, this implies $\partial_{\xi}v\geq 0$ across the interval $c_{s,-}\leq\xi\leq\xi_{w}$. From (4.7a), we must therefore have

At the bubble wall, $\mu(\xi=\xi_{w},\tilde{v}_{-})=v_{-}$, so for (4.27) to be satisfied, we require $v_{-}\geq c_{s,-}$. This is only possible for weak supersonic and Chapman-Jouguet detonations. Furthermore, strong detonations are prohibited as they cannot satisfy the boundary conditions at the wall and the centre of the bubble [68, 49, 8].

We note that for inverse transitions, $v(\xi)\leq 0$, so this condition becomes $v_{-}\leq c_{s,-}$, permitting only weak subsonic and Chapman-Jouguet inverse detonations. The minimum speed of the fluid in front of the wall is determined by the Chapman-Jouguet condition, $v_{-}=c_{s,-}$ ($v_{+}=c_{s,+}$ for inverse transitions). In general, this can be calculated by solving the matching conditions at the wall for $v_{+}$. For the $\mu\nu$ equation of state, we define

The positive root is called the Jouguet detonation velocity, $v_{J}^{\text{det}}$, and implies the wall is supersonic ($\xi_{w}>v_{J}^{\text{det}}>c_{s,+}$). The thermodynamic properties of the fluid beyond the wall are just those of the symmetric phase, so $T_{+}=T_{N}$ and $w_{+}=w_{N}$. From here on, we refer to weak subsonic deflagrations and weak supersonic detonations as just deflagrations and detonations, unless otherwise stated.

So far, we have constructed consistent solutions for bubble wall velocities $\xi_{w}<c_{s,+}$ (weak subsonic deflagrations) and $\xi_{w}>v_{J}^{\text{det}}>c_{s,+}$ (weak supersonic detonations), which are shown in Figure 2. In fact, another stable solution can be constructed for bubble walls that are supersonic and less than the Jouguet detonation velocity. Hydrodynamic simulations [48] suggest that strong deflagrations, which we previously mentioned were unstable, tend to decay to form a rarefaction wave similar to a detonation. This gives rise to a type of weak supersonic deflagration called a hybrid. The rarefaction wave causes the fluid behind the bubble wall to be disturbed, so these differ from usual deflagration solutions and are more appropriately described as a superposition of a Chapman-Jouguet deflagration and weak supersonic detonation [49] with $v_{-}=c_{s,-}$ and $v_{+}<c_{s,+}$, whose bubble wall velocity satisfies $c_{s,+}\leq\xi_{w}\leq v_{J}^{\text{det}}$.

The bag model is just a special case of the $\mu\nu$ model with $c_{s,+}^{2}=c_{s,-}^{2}=1/3$, so the methods used to construct the fluid profiles in either case are identical. We first initialise a
PTParams_Bag object that creates a class of phase transition parameters: the wall velocity $\xi_{w}$, nucleation temperature $T_{N}$, strength parameter at nucleation temperature $\alpha_{N}=\alpha_{\bar{\theta}}(T_{N})$, inverse phase transition duration $\beta$, mean bubble separation $R_{*}$, nucleation type (exponential or simultaneous), the universe constants class, and the speed of sound in the broken and symmetric phases $c_{s,\pm}$ (for the bag model, $c_{s,\pm}^{2}=1/3$).

The hydrodynamic mode is simply determined by computing the Jouguet velocity (4.28). The solver algorithm for the simplified equations of state is implemented in the
solve_profile function. For deflagration solutions, the fluid is at rest inside the bubble and beyond the shock, so we construct a solution across $\xi_{w}<\xi<\xi_{\text{sh}}$. Given some initial condition $v(\xi_{w})=\tilde{v}_{+}$, we solve the equations of motion, initially using $w_{+}=w_{N}$ and $T_{+}=T_{N}$. We integrate from $v(\xi_{w})=\tilde{v}_{+}$ to $v(\xi_{\text{end}})=0$, where $\xi_{\text{end}}$ is some arbitrary end-point of the integration. The shock front exists at some $\xi_{\text{sh}}\in(c_{s,+},\xi_{\text{end}})$, which is found by iteratively stepping through the solution until the first shock matching condition in (4.26), with $v_{2}=\xi_{\text{sh}}$, is satisfied. The solution is then truncated to the position of the shock to obtain the profiles across the interval $\xi_{w}<\xi<\xi_{\text{sh}}$. To obtain correctly scaled enthalpy and temperature profiles, we compute $w_{1}/w_{N}$ and $T_{1}/T_{N}$ from $\xi_{\text{sh}}$ and normalise the solution such that $w(\xi_{\text{sh}})=w_{1}/w_{N}$ and $T(\xi_{\text{sh}})=T_{1}/T_{N}$. The enthalpy of the fluid behind the shock is

from the first matching condition in (4.13). From (4.25), the temperature behind the shock is simply

To determine the initial condition $v(\xi_{w})=\tilde{v}_{+}$ for deflagrations, we construct a residual function alN_residual that computes the strength parameter $\alpha_{N}$ from the matching conditions at the wall and compares it to the input value of $\alpha_{N}$. To compute $\alpha_{N}$ using our ‘guess’ $\tilde{v}_{+}$, we first solve the fluid equations of motion with the guess value to obtain $w_{+}/w_{N}$, the enthalpy in front of the wall. We then compute $\alpha_{+}$ from the matching condition at the wall by inverting (4.24),

where $v_{-}=\xi_{w}$ for deflagrations. The definition of the strength parameter (4.23) gives us the relation

which we use to determine $\alpha_{N}^{\text{guess}}$ and then calculate the corresponding residual $|\alpha_{N}^{\text{guess}}-\alpha_{N}|$. The residual is computed in a bracket $(\tilde{v}_{+}^{\text{min}},\tilde{v}_{+}^{\text{max}})$ (most straightforwardly taken to be $(0,1)$) and minimised using a golden section minimisation method to find the correct $\tilde{v}_{+}$. In practice, the residual only exists in a neighbourhood of the true value for $\tilde{v}_{+}$. Outside of this neighbourhood, there are no solutions that satisfy both of the matching conditions at the bubble wall and shock front.

The energy density fluctuation profile $\lambda(\xi)$, defined as

is also computed. This quantity will be used in (7.8b) in the sound shell model. Here, $e(\xi)=e_{+}(T(\xi))$ for deflagrations. Using the $\mu\nu$ equation of state, we note that the relation (5.3) holds across the entire deflagration profile. That is, $T(\xi)/T_{N}=(w(\xi)/w_{N})^{1/4}$. The energy density fluctuation profile is therefore

Detonation solutions are considerably simpler. The fluid is at rest beyond the wall and comes to a halt inside the bubble at $\xi=c_{s,-}$, so we construct a solution across $c_{s,-}<\xi<\xi_{w}$. We use the initial conditions $(\xi_{0},w_{0},T_{0})=(\xi_{w},w_{-}/w_{N},T_{-}/T_{N})$ and integrate from $v(\xi_{w})=\tilde{v}_{-}$ to $v(c_{s,-})=0$, where $\tilde{v}_{-}=\mu(\xi_{w},|v_{-}|)$ and $v_{-}$ is determined from inverting the matching condition (4.24),

where $v_{+}=\xi_{w}$ and $\alpha_{+}=\alpha_{N}$ for detonations. The enthalpy just behind the wall is determined from the matching condition (4.11)

with $w_{+}=w_{N}$ for detonations. From the bag equation of state (4.19), $w_{\pm}=(1+c_{s,\pm}^{2})a_{\pm}T_{\pm}^{4}$. The temperature behind the wall is then

and the ratio $a_{+}/a_{-}$ is approximated to $1$ in this work. These initial conditions are then used to integrate the equations of motion (5.1) from $v(\xi_{w})=\tilde{v}_{-}$ to $v(c_{s,_{-}})=0$. The energy density fluctuation profile for detonations is slightly more complex with $e(\xi)=e_{-}(T(\xi))$, since the bag constant $\epsilon$ remains in the numerator of (5.6). To express $\lambda(\xi)$ in terms of the enthalpy, sound speed, and strength parameter, we solve (4.23) for $\epsilon$, and we find

Hybrid solutions require us to stitch together both a deflagration solution ($\xi_{w}<\xi<\xi_{\text{sh}}$) and a detonation solution ($c_{s,-}<\xi<\xi_{w}$). We do this using the usual deflagration/detonation solvers described above, but with the boundary condition $v_{-}=c_{s,-}$ in the deflagration part of the profiles and using the matching conditions at the wall (4.11) to obtain the initial condition for the detonation part of the profiles.

To construct the fluid profiles directly from the effective potential, we first initialise a PTParams_Veff object which takes in the same phase transition parameters as the simplified equations of state, except for the speed of sound in each phase. Instead, an EquationOfState object is passed in, which contains the pressure $p$ and energy density $e$ in each phase as a function of temperature from some effective potential. The equation of state data is stored in the class, and the speed of sound $c_{s}^{2}(T)=\mathrm{d}p/\mathrm{d}e$ in each phase is evaluated numerically.

The PTParams_Veff object is then used to initialise the FluidProfile class. Constructing fluid profiles for the exact equation of state is more involved than for simplified equations of state, since the algebraic forms of the thermodynamic quantities $p(T)$, $e(T)$, $w(T)$, and $s(T)$ are unknown. Consequently, the matching conditions in equation (4.12) must be imposed numerically using a root-finding algorithm. Throughout most of the code, we use a bounded 2D Newton’s method solver to find the root of

for matching at the wall and

for matching at the shock. Here, $p_{\pm}$ and $e_{\pm}$ are defined in (4.17), and $p_{1,2}$ and $e_{1,2}$ are defined similarly in (4.18).