aainstitutetext: Departament de Física Quàntica i Astrofísica and Institut de Ciències del Cosmos (ICC), Universitat de Barcelona, Martí  i Franquès 1, ES-08028, Barcelona, Spain.bbinstitutetext: Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluís Companys 23, ES-08010, Barcelona, Spain.ccinstitutetext: CPHT, CNRS, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

Hydrodynamics of Relativistic Superheated Bubbles

Yago Bea a    Jorge Casalderrey-Solana a,b    David Mateos c    Mikel Sanchez-Garitaonandia yagobea@icc.ub.edu jorge.casalderrey@ub.edu dmateos@fqa.ub.edu mikel.sanchez@polytechnique.edu
Abstract

Relativistic, charged, superheated bubbles may play an important role in neutron star mergers if first-order phase transitions are present in the phase diagram of Quantum Chromodynamics. We describe the properties of these bubbles in the hydrodynamic regime. We find two qualitative differences with supercooled bubbles. First, the pressure inside an expanding superheated bubble can be higher or lower than the pressure outside the bubble. Second, some fluid flows develop metastable regions behind the bubble wall. We consider the possible role of a conserved charge akin to baryon number. The fluid flow profiles are unaffected by this charge if the speed of sound is constant in each phase, but they are modified for more general equations of state. We compute the efficiency factor relevant for gravitational wave production.

CPHT-RR050.062024

1 Introduction

A variety of arguments suggest that at least two first-order phase transitions (FOPT) may be present in the phase diagram of Quantum Chromodynamics (QCD) as a function of temperature and baryon chemical potential Stephanov:2004wx ; Alford:2007xm ; Fukushima:2010bq ; Guenther:2022wcr , as sketched in Fig. 1. One is the transition from hadronic matter to quark matter. The other is the transition from a non-superconducting to a color-superconducting phase. While these transitions are well motivated, rigorously establishing their existence is a fundamental open problem in nuclear and particle physics whose solution has resisted both theoretical and experimental attempts for decades.

Gravitational waves produced in neutron star (NS) mergers could provide direct experimental access to these phase transitions Casalderrey-Solana:2022rrn . Numerical simulations of NS mergers based on equations of state (EoS) with a hadronic-quark matter phase transition include Most:2018eaw ; Most:2019onn ; Ecker:2019xrw ; Prakash:2021wpz ; Weih:2019xvw ; Tootle:2022pvd . These studies have shown that the dynamics of the merger results in the formation of regions in which the matter is sufficiently heated and/or compressed that the thermodynamically preferred phase is the quark-matter phase. In other words, the matter in these regions is pushed along the black dotted curve in Fig. 1. For brevity, we will refer to these regions simply as “superheated regions”.

Refer to caption
Figure 1: Two possible phase transitions in QCD, indicated by the solid red curves. T𝑇Titalic_T and μ𝜇\muitalic_μ are the temperature and the baryon chemical potential, respectively. The dotted black curve on the left panel shows a possible evolution of a region of a NS merger as this region is heated and/or compressed. The points dubbed A𝐴Aitalic_A, Asuperscript𝐴A^{\prime}italic_A start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT and C𝐶Citalic_C correspond to the states shown in Fig. 2. See text and Ref. Casalderrey-Solana:2022rrn .

Although no simulation based on an EoS with a color-superconducting phase has been performed, the large baryon densities found in existing simulations make it conceivable that color-superconducting matter may also be formed in NS mergers.

The energy density along the black dotted curve in Fig. 1 displays the multivalued form characteristic of a FOPT, as shown in Fig. 2. The superheated region is pushed from left to right along the lower branch of the phase diagram. Once the superheating is large enough, namely once the region of interest is sufficiently deep into the metastable branch, bubbles of the stable phase begin to nucleate. The point where this happens is labeled “B𝐵Bitalic_B” in Fig. 2, the nucleated phase is labelled “C𝐶Citalic_C”, and the transition is indicated by a vertical, solid, black arrow. We will refer to these nucleated bubbles as “superheated bubbles”. The dynamics of these bubbles can produce gravitational waves that would provide direct experimental access to the QCD phase transition Casalderrey-Solana:2022rrn . The order-of-magnitude estimate in this reference shows that the frequency of these gravitational waves falls roughly in the MHz range, and that they may be potentially observable with future detectors.

Refer to caption
Figure 2: Energy density as a function of the position along the black dotted curve in Fig. 1. Both T𝑇Titalic_T and μ𝜇\muitalic_μ increase from left to right. The black, dotted, vertical line indicates the location of the phase transition, determined by the condition that the states A𝐴Aitalic_A and Asuperscript𝐴A^{\prime}italic_A start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT have the same free energy density. The green and blue curves indicate stable and metastable states, respectively. As some region of the NS merger is sufficiently heated and/or compressed, it enters the lower metastable branch. At the point B𝐵Bitalic_B bubbles of the preferred state C𝐶Citalic_C on the upper stable branch are nucleated. The direction of this phase transition is the opposite of that in a supercooled phase transition, which would take place as indicated by the vertical, dotted, red arrow.

The direction of the “superheated” transition is opposite to that of a cosmological FOPT. In this case the Universe is supercooled, namely pushed from right to left along the upper branch into the upper metastable region, until it transitions down to the lower stable branch, as indicated in Fig. 2 by the dashed, vertical, red arrow. We will refer to the nucleated bubbles in this case as “supercooled bubbles”. This type of transition will also occur in a NS merger if the superheated region cools down again.

The dynamics of relativistic supercooled bubbles has been extensively studied in the literature, motivated by their potential implications for cosmological FOPT — see e.g. Hindmarsh:2020hop for a review. Motivated by their possible role in NS mergers, in this paper we initiate the study of relativistic, charged, superheated bubbles, focusing on their hydrodynamic properties in the simple case of a bag-model, conformal EoS. We find two qualitative differences compared to supercooled bubbles. First, the pressure inside a superheated bubble may be higher or lower than the pressure outside the bubble. In contrast, in supercooled bubbles the inside pressure is always higher than the outside pressure. Second, some fluid flows develop regions behind the bubble wall that are metastable for any choice of parameters in the EoS, which would therefore decay at sufficiently long times. In contrast, in supercooled transitions there always exists a choice of parameters for which the flow behind the wall only explores stable states.

The physics of relativistic and non-relativistic bubbles is expected to be qualitatively similar if the two phases that determine the phase transition have comparable compressibilities in the range of pressures accessed by the bubble dynamics. For example, relativistic supercooled bubbles find a non-relativistic analog in non-relativistic chemical combustions, for this reason often referred to as “relativistic combustions”. Thus, one may expect that these similarities are also present for relativisitic superheated bubbles. However, a qualitative difference appears in a particular case of obvious phenomenological interest: the nucleation and expansion of vapor bubbles in a superheated liquid. In this case the liquid is effectively incompressible, and this leads to a bubble wall velocity that decreases as t1/2superscript𝑡12t^{-1/2}italic_t start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT at asymptotically late times. For completeness, we review this case in Appendix A.

Part of the contents of this paper overlap with those in Ref. Barni:2024lkj , which appeared while this paper was being typewritten.

2 Self-similar profiles

At late times, when the bubbles have become macroscopic in size, the lack of any macroscopic scale implies that these bubbles are expected to become self-similar. This means that the fluid flow profile is a function of only ξ=r/t𝜉𝑟𝑡\xi=r/titalic_ξ = italic_r / italic_t, with r𝑟ritalic_r the radius of the spherical bubble. Besides the bubble wall and possible shocks in the system, the very diluted gradients mean that a perfect-fluid description should be a valid approximation of the plasma. In this case, the constitutive equations for the stress tensor and the conserved charge take the form

Tμν=wuμuν+pgμν,jμ=nuμ,formulae-sequencesuperscript𝑇𝜇𝜈𝑤superscript𝑢𝜇superscript𝑢𝜈𝑝superscript𝑔𝜇𝜈superscript𝑗𝜇𝑛superscript𝑢𝜇T^{\mu\nu}=wu^{\mu}u^{\nu}+pg^{\mu\nu},\quad\quad j^{\mu}=nu^{\mu}\,,italic_T start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT = italic_w italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT italic_u start_POSTSUPERSCRIPT italic_ν end_POSTSUPERSCRIPT + italic_p italic_g start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT , italic_j start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT = italic_n italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT , (1)

with p𝑝pitalic_p the pressure, w=e+p𝑤𝑒𝑝w=e+pitalic_w = italic_e + italic_p the enthalpy, e𝑒eitalic_e the energy density, n𝑛nitalic_n the charge density of the conserved charge and uμ=γ(1,v)superscript𝑢𝜇𝛾1𝑣u^{\mu}=\gamma(1,\vec{v})italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT = italic_γ ( 1 , over→ start_ARG italic_v end_ARG ) the fluid four-velocity. The thermodynamic variables are not all independent but are related to the the temperature T𝑇Titalic_T and chemical potential of the plasma μ𝜇\muitalic_μ through the EoS. In this section we will not assume any specific EoS (but later in this section we will assume a constant speed of sound).

The equations of motion are simply the conservation equations of the stress-tensor and the current:

μTμν=0,μjμ=0.formulae-sequencesubscript𝜇superscript𝑇𝜇𝜈0subscript𝜇superscript𝑗𝜇0\partial_{\mu}T^{\mu\nu}=0\,,\quad\quad\partial_{\mu}j^{\mu}=0\,.∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_T start_POSTSUPERSCRIPT italic_μ italic_ν end_POSTSUPERSCRIPT = 0 , ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_j start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT = 0 . (2)

As in the neutral case (see e.g. Espinosa:2010hh ), we can project the motion of the flow along the direction parallel and orthogonal to the velocity of the fluid. For this purpose we multiply the conservation equations by uμ=γ(1,v)superscript𝑢𝜇𝛾1𝑣u^{\mu}=\gamma(1,\vec{v})italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT = italic_γ ( 1 , over→ start_ARG italic_v end_ARG ) and u¯μ=γ(v,v/v)superscript¯𝑢𝜇𝛾𝑣𝑣𝑣\bar{u}^{\mu}=\gamma(v,\vec{v}/v)over¯ start_ARG italic_u end_ARG start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT = italic_γ ( italic_v , over→ start_ARG italic_v end_ARG / italic_v ). Using that uμνuμ=0subscript𝑢𝜇subscript𝜈superscript𝑢𝜇0u_{\mu}\partial_{\nu}u^{\mu}=0italic_u start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT ∂ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT = 0, we obtain the following equations:

wμuμ+uμμe=0𝑤subscript𝜇superscript𝑢𝜇superscript𝑢𝜇subscript𝜇𝑒0\displaystyle w\,\partial_{\mu}u^{\mu}+u^{\mu}\,\partial_{\mu}e=0italic_w ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT + italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_e = 0 ,absent\displaystyle\,,, (3)
wu¯νuμμuνu¯μμp=0𝑤superscript¯𝑢𝜈superscript𝑢𝜇subscript𝜇subscript𝑢𝜈superscript¯𝑢𝜇subscript𝜇𝑝0\displaystyle w\,\bar{u}^{\nu}u^{\mu}\,\partial_{\mu}u_{\nu}-\bar{u}^{\mu}\,% \partial_{\mu}p=0italic_w over¯ start_ARG italic_u end_ARG start_POSTSUPERSCRIPT italic_ν end_POSTSUPERSCRIPT italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT - over¯ start_ARG italic_u end_ARG start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_p = 0 ,absent\displaystyle\,,,
nμuμ+uμμn=0𝑛subscript𝜇superscript𝑢𝜇superscript𝑢𝜇subscript𝜇𝑛0\displaystyle n\,\partial_{\mu}u^{\mu}+u^{\mu}\,\partial_{\mu}n=0italic_n ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT + italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_n = 0 .absent\displaystyle\,..

The self-similarity assumption implies that

uμμ=γt(ξv)ξ,superscript𝑢𝜇subscript𝜇𝛾𝑡𝜉𝑣subscript𝜉\displaystyle u^{\mu}\partial_{\mu}=-\frac{\gamma}{t}(\xi-v)\partial_{\xi}\,,italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT = - divide start_ARG italic_γ end_ARG start_ARG italic_t end_ARG ( italic_ξ - italic_v ) ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT , (4)
u¯μμ=γt(1ξv)ξ,superscript¯𝑢𝜇subscript𝜇𝛾𝑡1𝜉𝑣subscript𝜉\displaystyle\bar{u}^{\mu}\partial_{\mu}=\frac{\gamma}{t}(1-\xi v)\partial_{% \xi},over¯ start_ARG italic_u end_ARG start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT = divide start_ARG italic_γ end_ARG start_ARG italic_t end_ARG ( 1 - italic_ξ italic_v ) ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT ,
μuμ=(1ξv)γ3tξ+γvr(d1),subscript𝜇superscript𝑢𝜇1𝜉𝑣superscript𝛾3𝑡subscript𝜉𝛾𝑣𝑟𝑑1\displaystyle\partial_{\mu}u^{\mu}=(1-\xi v)\frac{\gamma^{3}}{t}\partial_{\xi}% +\frac{\gamma v}{r}(d-1),∂ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT = ( 1 - italic_ξ italic_v ) divide start_ARG italic_γ start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG start_ARG italic_t end_ARG ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT + divide start_ARG italic_γ italic_v end_ARG start_ARG italic_r end_ARG ( italic_d - 1 ) ,

where d𝑑ditalic_d is the number of spatial dimensions. Substituting these expressions into the equations of motion leads to the following set of differential equations:

(ξv)ξew𝜉𝑣subscript𝜉𝑒𝑤\displaystyle(\xi-v)\frac{\partial_{\xi}e}{w}( italic_ξ - italic_v ) divide start_ARG ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_e end_ARG start_ARG italic_w end_ARG =\displaystyle== γ2(1ξv)ξv+(d1)vξ,superscript𝛾21𝜉𝑣subscript𝜉𝑣𝑑1𝑣𝜉\displaystyle\gamma^{2}(1-\xi v)\partial_{\xi}v+(d-1)\frac{v}{\xi},italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( 1 - italic_ξ italic_v ) ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_v + ( italic_d - 1 ) divide start_ARG italic_v end_ARG start_ARG italic_ξ end_ARG , (5)
(1ξv)ξpw1𝜉𝑣subscript𝜉𝑝𝑤\displaystyle(1-\xi v)\frac{\partial_{\xi}p}{w}( 1 - italic_ξ italic_v ) divide start_ARG ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_p end_ARG start_ARG italic_w end_ARG =\displaystyle== γ2(ξv)ξv,superscript𝛾2𝜉𝑣subscript𝜉𝑣\displaystyle\gamma^{2}(\xi-v)\partial_{\xi}v,italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ξ - italic_v ) ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_v , (6)
(ξv)ξnn𝜉𝑣subscript𝜉𝑛𝑛\displaystyle(\xi-v)\frac{\partial_{\xi}n}{n}( italic_ξ - italic_v ) divide start_ARG ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_n end_ARG start_ARG italic_n end_ARG =\displaystyle== γ2(1ξv)ξv+(d1)vξ.superscript𝛾21𝜉𝑣subscript𝜉𝑣𝑑1𝑣𝜉\displaystyle\gamma^{2}(1-\xi v)\partial_{\xi}v+(d-1)\frac{v}{\xi}\,.italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( 1 - italic_ξ italic_v ) ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_v + ( italic_d - 1 ) divide start_ARG italic_v end_ARG start_ARG italic_ξ end_ARG . (7)

These can be recast into a more convenient form:

γ2(1ξv)(ν(ξ,v)2cs21)ξv=(d1)vξ,superscript𝛾21𝜉𝑣𝜈superscript𝜉𝑣2superscriptsubscript𝑐𝑠21subscript𝜉𝑣𝑑1𝑣𝜉\displaystyle\gamma^{2}(1-\xi v)\left(\frac{\nu(\xi,v)^{2}}{c_{s}^{2}}-1\right% )\partial_{\xi}v=(d-1)\frac{v}{\xi}\,,italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( 1 - italic_ξ italic_v ) ( divide start_ARG italic_ν ( italic_ξ , italic_v ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG - 1 ) ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_v = ( italic_d - 1 ) divide start_ARG italic_v end_ARG start_ARG italic_ξ end_ARG , (8a)
ξww=γ2(1+1cs2)ν(ξ,v)ξv,subscript𝜉𝑤𝑤superscript𝛾211superscriptsubscript𝑐𝑠2𝜈𝜉𝑣subscript𝜉𝑣\displaystyle\frac{\partial_{\xi}w}{w}=\gamma^{2}\left(1+\frac{1}{c_{s}^{2}}% \right)\nu(\xi,v)\partial_{\xi}v\,,divide start_ARG ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_w end_ARG start_ARG italic_w end_ARG = italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( 1 + divide start_ARG 1 end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) italic_ν ( italic_ξ , italic_v ) ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_v , (8b)
ξnn=11+cs2ξww,subscript𝜉𝑛𝑛11superscriptsubscript𝑐𝑠2subscript𝜉𝑤𝑤\displaystyle\frac{\partial_{\xi}n}{n}=\frac{1}{1+c_{s}^{2}}\frac{\partial_{% \xi}w}{w}\,,divide start_ARG ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_n end_ARG start_ARG italic_n end_ARG = divide start_ARG 1 end_ARG start_ARG 1 + italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG divide start_ARG ∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_w end_ARG start_ARG italic_w end_ARG , (8c)

where

cs2=ep+nwnpsuperscriptsubscript𝑐𝑠2subscript𝑒𝑝𝑛𝑤subscript𝑛𝑝c_{s}^{2}=\partial_{e}p+\frac{n}{w}\,\partial_{n}pitalic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = ∂ start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT italic_p + divide start_ARG italic_n end_ARG start_ARG italic_w end_ARG ∂ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT italic_p (9)

is the speed of sound and

ν(ξ,v)=ξv1ξv𝜈𝜉𝑣𝜉𝑣1𝜉𝑣\nu(\xi,v)=\frac{\xi-v}{1-\xi v}italic_ν ( italic_ξ , italic_v ) = divide start_ARG italic_ξ - italic_v end_ARG start_ARG 1 - italic_ξ italic_v end_ARG (10)

is the Lorentz-boosted velocity.

If each branch of the EoS is conformal, namely, if the pressure in each phase takes the form

p(T,μ)=Td+1f(μT)+const.,𝑝𝑇𝜇superscript𝑇𝑑1𝑓𝜇𝑇const.p(T,\mu)=T^{d+1}f\left(\frac{\mu}{T}\right)+\mbox{const.}\,,italic_p ( italic_T , italic_μ ) = italic_T start_POSTSUPERSCRIPT italic_d + 1 end_POSTSUPERSCRIPT italic_f ( divide start_ARG italic_μ end_ARG start_ARG italic_T end_ARG ) + const. , (11)

with f𝑓fitalic_f an arbitrary function, then the speed of sound is constant cs2=1/dsuperscriptsubscript𝑐𝑠21𝑑c_{s}^{2}=1/ditalic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = 1 / italic_d. In this case (8) becomes a nested system that can be sequentially solved for v𝑣vitalic_v, w𝑤witalic_w and n𝑛nitalic_n, similarly to what happens in the neutral case Espinosa:2010hh . The nested structure means that the addition of charge does not change the properties of self-similar flows for conformal EoS. In fact, the last equation simply implies that nw1/(1+cs2)proportional-to𝑛superscript𝑤11superscriptsubscript𝑐𝑠2n\propto w^{1/(1+c_{s}^{2})}italic_n ∝ italic_w start_POSTSUPERSCRIPT 1 / ( 1 + italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_POSTSUPERSCRIPT.

In contrast, if the EoS is non-conformal, the speed of sound will depend on the local properties of the plasma, cs=cs(e,n)subscript𝑐𝑠subscript𝑐𝑠𝑒𝑛c_{s}=c_{s}(e,n)italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ( italic_e , italic_n ), making equations (8) coupled to one another. In this case the effect of the charge is already present at the level of the self-similar flow. We leave studies of non-conformal theories for future work.

At the location of the bubble wall or of the possible shocks that will form during the expansion of the bubble, the perfect fluid description will break down. However, we can integrate the conservation laws, in the radial direction, in the rest frame of these discontinuities. Doing so leads to the following matching conditions:

w+v+2γ+2+p+subscript𝑤superscriptsubscript𝑣2superscriptsubscript𝛾2subscript𝑝\displaystyle w_{+}v_{+}^{2}\gamma_{+}^{2}+p_{+}italic_w start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_γ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_p start_POSTSUBSCRIPT + end_POSTSUBSCRIPT =wv2γ2+p,absentsubscript𝑤superscriptsubscript𝑣2superscriptsubscript𝛾2subscript𝑝\displaystyle=w_{-}v_{-}^{2}\gamma_{-}^{2}+p_{-}\,,= italic_w start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_γ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_p start_POSTSUBSCRIPT - end_POSTSUBSCRIPT , (12)
w+v+γ+2subscript𝑤subscript𝑣superscriptsubscript𝛾2\displaystyle w_{+}v_{+}\gamma_{+}^{2}italic_w start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT =wvγ2,absentsubscript𝑤subscript𝑣superscriptsubscript𝛾2\displaystyle=w_{-}v_{-}\gamma_{-}^{2}\,,= italic_w start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ,
n+v+γ+subscript𝑛subscript𝑣subscript𝛾\displaystyle n_{+}v_{+}\gamma_{+}italic_n start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT =nvγ,absentsubscript𝑛subscript𝑣subscript𝛾\displaystyle=n_{-}v_{-}\gamma_{-}\,,= italic_n start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ,

where +++(--) refers to the plasma in front (behind) the discontinuity, the velocities v±subscript𝑣plus-or-minusv_{\pm}italic_v start_POSTSUBSCRIPT ± end_POSTSUBSCRIPT are measured in the rest frame of the discontinuity, and γ𝛾\gammaitalic_γ is the usual Lorentz factor. Rearranging the first two matching conditions we arrive at

v+v=p+pe+e,v+v=e+p+e++p,n+v+γ+=nvγ.formulae-sequencesubscript𝑣subscript𝑣subscript𝑝subscript𝑝subscript𝑒subscript𝑒formulae-sequencesubscript𝑣subscript𝑣subscript𝑒subscript𝑝subscript𝑒subscript𝑝subscript𝑛subscript𝑣subscript𝛾subscript𝑛subscript𝑣subscript𝛾v_{+}v_{-}=\frac{p_{+}-p_{-}}{e_{+}-e_{-}}\,,\quad\quad\frac{v_{+}}{v_{-}}=% \frac{e_{-}+p_{+}}{e_{+}+p_{-}}\,,\quad\quad n_{+}v_{+}\gamma_{+}=n_{-}v_{-}% \gamma_{-}.italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT = divide start_ARG italic_p start_POSTSUBSCRIPT + end_POSTSUBSCRIPT - italic_p start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG italic_e start_POSTSUBSCRIPT + end_POSTSUBSCRIPT - italic_e start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG , divide start_ARG italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG = divide start_ARG italic_e start_POSTSUBSCRIPT - end_POSTSUBSCRIPT + italic_p start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG italic_e start_POSTSUBSCRIPT + end_POSTSUBSCRIPT + italic_p start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG , italic_n start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT = italic_n start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT . (13)

Eqs. (8) and (13) determine the entire fluid flow as a function of three input parameters: the wall velocity ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT, the nucleation energy density ensubscript𝑒𝑛e_{n}italic_e start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT and the nucleation charge density nnsubscript𝑛𝑛n_{n}italic_n start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT. In the following we will solve these equations in the case of spherical bubbles in three spatial dimensions (d=3𝑑3d=3italic_d = 3) and a constant speed of sound.

When cssubscript𝑐𝑠c_{s}italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT is constant we can solve for the velocity in (8) independently of the rest of the variables. By assuming different initial conditions v(ξ0)=v0𝑣subscript𝜉0subscript𝑣0v(\xi_{0})=v_{0}italic_v ( italic_ξ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) = italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT one can obtain all possible flows that are self-similar, shown in Fig. 3. Contrary to the supercooled case, we seek solutions where at ξ>0𝜉0\xi>0italic_ξ > 0 (expanding bubbles) the flow has negative speed v<0𝑣0v<0italic_v < 0 (energy inflow). In order to produce Fig. 3, we chose cs2=1/3superscriptsubscript𝑐𝑠213c_{s}^{2}=1/3italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = 1 / 3. We expect that different values of the speed of sound will change the picture quantitatively but not qualitatively.

Refer to caption
Figure 3: Self-similar flows for superheated bubbles in the case cs2=1/3superscriptsubscript𝑐𝑠213c_{s}^{2}=1/3italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = 1 / 3. The vertical dotted line lies at ξ=cs𝜉subscript𝑐𝑠\xi=c_{s}italic_ξ = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT, marking the beginning of detonations. The dashed blue curve corresponds to the locus ν(ξ,v)ξ=cs2𝜈𝜉𝑣𝜉superscriptsubscript𝑐𝑠2\nu(\xi,v)\xi=c_{s}^{2}italic_ν ( italic_ξ , italic_v ) italic_ξ = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, where the shock must be located . The dotted-dashed red curve indicates the points where the boosted velocity ν(ξ,v)=cs𝜈𝜉𝑣subscript𝑐𝑠\nu(\xi,v)=c_{s}italic_ν ( italic_ξ , italic_v ) = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT and ξvsubscript𝜉𝑣\partial_{\xi}v\rightarrow\infty∂ start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_v → ∞. Expanding bubbles in the shaded area are not allowed.

For v0𝑣0v\rightarrow 0italic_v → 0, all flows tend to ξ=cs𝜉subscript𝑐𝑠\xi=c_{s}italic_ξ = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT, doing so from low to high values of ξ𝜉\xiitalic_ξ. This implies that flows that start at ξ<cs𝜉subscript𝑐𝑠\xi<c_{s}italic_ξ < italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT can continuously connect to a v=0𝑣0v=0italic_v = 0 state, while flows starting at ξ>cs𝜉subscript𝑐𝑠\xi>c_{s}italic_ξ > italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT will inevitably have to develop a discontinuity to do so. Moreover, in the region ξ<cs𝜉subscript𝑐𝑠\xi<c_{s}italic_ξ < italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT another discontinuity may occur, for which the matching conditions reduce to ν(ξ,v)ξ=cs2𝜈𝜉𝑣𝜉superscriptsubscript𝑐𝑠2\nu(\xi,v)\xi=c_{s}^{2}italic_ν ( italic_ξ , italic_v ) italic_ξ = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. This discontinuity is a shock and its location is displayed in Fig. 3 as a dashed blue curve.

With this information we can already see the type of fluid flows that we can get. For low velocity we expect a so-called “deflagration”, namely a flow in front of the wall that smoothly connects with the asymptotic state at rest. For high velocity, we expect a “detonation”, namely a flow behind the wall that will develop a shock, connecting with the state at rest inside the bubble. Finally, we can also have hybrid solutions, which combine features from both deflagrations and detonations. Hybrid flows consist of a flow behind the wall (as in a detonation) ending in a shock, plus a flow ahead of the wall that smoothly connects with v=0𝑣0v=0italic_v = 0 (as in a deflagration).

There is yet another kind of flows, those that are contained in the shaded area of Fig. 3. Given that they would become multivalued, these flows cannot cross the dot-dashed line, and they must follow to higher ξ𝜉\xiitalic_ξ. As they cannot smoothly connect to a state with v=0𝑣0v=0italic_v = 0, they jump discontinously to it at the shock location. This shock is different to the one present in hybrids and detonations as the flow in motion is behind and not ahead. An entropic analysis shows that the shock front in Fig. 3 increases entropy only if the flow in motion is ahead. Hence, no flow starting in the mentioned area is realizable.

In the next subsection we explore in more detail the three type of solutions just mentioned.

3 Solutions for a conformal equation of state

In order to analyze the fluid flows in more detail, a choice of an EoS is needed. Assuming two branches such that cs2=1/3superscriptsubscript𝑐𝑠213c_{s}^{2}=1/3italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = 1 / 3 in each of them, we choose

pH(T,μ)subscript𝑝𝐻𝑇𝜇\displaystyle p_{H}(T,\mu)italic_p start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT ( italic_T , italic_μ ) =T4(aH+bHμ2T2+cHμ4T4)ϵ,absentsuperscript𝑇4subscript𝑎𝐻subscript𝑏𝐻superscript𝜇2superscript𝑇2subscript𝑐𝐻superscript𝜇4superscript𝑇4italic-ϵ\displaystyle=T^{4}\left(a_{H}+b_{H}\frac{\mu^{2}}{T^{2}}+c_{H}\frac{\mu^{4}}{% T^{4}}\right)-\epsilon\,,= italic_T start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ( italic_a start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT + italic_b start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT divide start_ARG italic_μ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_T start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_c start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT divide start_ARG italic_μ start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG start_ARG italic_T start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG ) - italic_ϵ , (14)
pL(T,μ)subscript𝑝𝐿𝑇𝜇\displaystyle p_{L}(T,\mu)italic_p start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( italic_T , italic_μ ) =T4(aL+bLμ2T2+cLμ4T4),absentsuperscript𝑇4subscript𝑎𝐿subscript𝑏𝐿superscript𝜇2superscript𝑇2subscript𝑐𝐿superscript𝜇4superscript𝑇4\displaystyle=T^{4}\left(a_{L}+b_{L}\frac{\mu^{2}}{T^{2}}+c_{L}\frac{\mu^{4}}{% T^{4}}\right)\,,= italic_T start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ( italic_a start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT + italic_b start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT divide start_ARG italic_μ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_T start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_c start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT divide start_ARG italic_μ start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG start_ARG italic_T start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG ) ,

where aH,Lsubscript𝑎𝐻𝐿a_{H,L}italic_a start_POSTSUBSCRIPT italic_H , italic_L end_POSTSUBSCRIPT, bH,Lsubscript𝑏𝐻𝐿b_{H,L}italic_b start_POSTSUBSCRIPT italic_H , italic_L end_POSTSUBSCRIPT and cH,Lsubscript𝑐𝐻𝐿c_{H,L}italic_c start_POSTSUBSCRIPT italic_H , italic_L end_POSTSUBSCRIPT reflect the number of degrees of freedom in each phase and “H𝐻Hitalic_H” and “L𝐿Litalic_L” refer to the high- and the low-energy phases, respectively. The constant ϵitalic-ϵ\epsilonitalic_ϵ is sometimes referred to as the “bag constant” and it measures the energy difference between the two phases in vacuum, that is, at T=μ=0𝑇𝜇0T=\mu=0italic_T = italic_μ = 0. For the choices aH=bH=cH=1subscript𝑎𝐻subscript𝑏𝐻subscript𝑐𝐻1a_{H}=b_{H}=c_{H}=1italic_a start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = italic_b start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = 1 and aL=bL=cL=1/2subscript𝑎𝐿subscript𝑏𝐿subscript𝑐𝐿12a_{L}=b_{L}=c_{L}=1/2italic_a start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT = italic_b start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT = 1 / 2, the phase diagram is shown in Fig. 4: There is a curve of FOPT transition at T2+μ2=2ϵsuperscript𝑇2superscript𝜇22italic-ϵT^{2}+\mu^{2}=\sqrt{2\epsilon}italic_T start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_μ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = square-root start_ARG 2 italic_ϵ end_ARG. We do not choose this EoS because we expect it to provide a realistic description of QCD but because it is the simplest example in which an investigation of superheated bubbles is possible. We will report on more realistic examples elsewhere.

Refer to caption
Figure 4: Phase diagram for the choice aH=bH=cH=1subscript𝑎𝐻subscript𝑏𝐻subscript𝑐𝐻1a_{H}=b_{H}=c_{H}=1italic_a start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = italic_b start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = 1 and aL=bL=cL=1/2subscript𝑎𝐿subscript𝑏𝐿subscript𝑐𝐿12a_{L}=b_{L}=c_{L}=1/2italic_a start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT = italic_b start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT = 1 / 2. “L𝐿Litalic_L” and “H𝐻Hitalic_H” refer to the low- and high-energy phases, respectively. All quantities are measured in units of the bag constant, ϵitalic-ϵ\epsilonitalic_ϵ, which is the energy difference between the two phases at T=μ=0𝑇𝜇0T=\mu=0italic_T = italic_μ = 0.

Having specified an EoS, we can now solve Eqs. (8) together with the matching conditions (13) in order to obtain the entire flow.

As mentioned earlier, for low, subsonic wall velocities ξw<cssubscript𝜉𝑤subscript𝑐𝑠\xi_{w}<c_{s}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT < italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT we have deflagration-type solutions. These are characterized by the fact that the fluid behind the wall is at rest. In contrast, the fluid in front of the wall gets into motion perturbed by the movement of the wall. The fluid in motion continuously connects at large distances with the superheated fluid at rest. An example of a deflagration is shown in the top panels of Figs. 5 and 6. Contrary to what happens in the supercooled case Espinosa:2010hh , the flow velocity is negative everywhere and, notably, there is no shock present in the system. The fluid flow in front of the wall continuously connects to the metastable state at ξ=cs𝜉subscript𝑐𝑠\xi=c_{s}italic_ξ = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT.

Refer to caption
Refer to caption
Refer to caption
Refer to caption
Refer to caption
Refer to caption
Figure 5: Fluid flows for superheated bubbles in the case of deflagrations (top), detonations (middle) and hybrids (bottom), with parameters ξw={0.35,0.8,0.5}subscript𝜉𝑤0.350.80.5\xi_{w}=\{0.35,0.8,0.5\}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT = { 0.35 , 0.8 , 0.5 } and αn={0.0614,0.143,0.0421}subscript𝛼𝑛0.06140.1430.0421\alpha_{n}=\{0.0614,0.143,0.0421\}italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT = { 0.0614 , 0.143 , 0.0421 }, respectively. The left panel shows the fluid velocity. The right panel shows the enthalpy (solid-black) and the charge density (dashed-black) normalized to their values at infinity. Grey regions correspond to the H𝐻Hitalic_H-phase and white regions to the L𝐿Litalic_L-phase, with the boundary being the bubble wall. Vertical dotted lines indicate the speed of sound.
Refer to caption
Refer to caption
Refer to caption
Refer to caption
Refer to caption
Refer to caption
Figure 6: Energy density (left) and pressure (right) profiles, normalized to their values at infinity, of the fluid flows shown in Fig. 5.

In the case of a high, supersonic wall velocity ξw>cssubscript𝜉𝑤subscript𝑐𝑠\xi_{w}>c_{s}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT > italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT, the solution will be of detonation type. In this case, the wall is moving too fast for the fluid ahead to respond and so it remains at rest. As the wall passes the fluid is set in motion and, as it cannot smoothly connect with a state of v=0𝑣0v=0italic_v = 0, it develops a shock down the stream. The shock then allows to connect with a fluid at rest inside the bubble. An example of a detonation is depicted in the middle panels of Figs. 5 and 6. As in the case of deflagrations, the detonation flow differs qualitatively from the supercooled case. In the latter, a rarefaction wave gets formed behind the wall, continuously connecting with the state at rest inside the bubble at ξw=cssubscript𝜉𝑤subscript𝑐𝑠\xi_{w}=c_{s}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT. Instead, in the current case the fluid dragged by the motion of the wall develops a shock that connects with the state inside the bubble.

Finally, for wall velocities close to but below the speed of sound, cs2ξwcssuperscriptsubscript𝑐𝑠2subscript𝜉𝑤subscript𝑐𝑠c_{s}^{2}\leq\xi_{w}\leq c_{s}italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ≤ italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ≤ italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT, a new type of flow can be constructed combining both the flows observed in detonations and deflagrations, giving rise to a hybrid solution. An example of this is displayed in the bottom panels of Figs. 5 and 6. The defining feature is that, in the rest frame of the bubble wall, the flow in front of it moves at the speed of sound v+=cssubscript𝑣subscript𝑐𝑠v_{+}=c_{s}italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT. As in any deflagration, the flow ahead of the wall continuously connects with the asymptotic state while the flow behind develops a shock that takes the fluid to rest. Hybrids appear in the current case for subsonic wall speeds, in contrast with the supercooled case Espinosa:2010hh .

The deflagration shown in Fig. 6(top) is an example of a superheated flow in which the pressure inside the bubble (at ξ=0𝜉0\xi=0italic_ξ = 0) is lower than the pressure outside the bubble (at ξ=1𝜉1\xi=1italic_ξ = 1), whereas the detonations (middle panel) and the hybrids (bottom panel) shown in Fig. 6 are examples in which the inside pressure is higher than the outside pressure. By constructing a large family of expanding superheated bubbles we have observed that all three types of flows can have either ordering of the pressures — see Fig. 7.

The specific detonation shown in Fig. 6(middle) has been chosen to illustrate that, for some flows, the pressure in a region just behind the wall, in this case in the region ξ0.8less-than-or-similar-to𝜉0.8\xi\lesssim 0.8italic_ξ ≲ 0.8, can be negative. This is necessarily lower than the pressure of the L𝐿Litalic_L-phase, which is always positive regardless of the temperature and of the values of the parameters aH,Lsubscript𝑎𝐻𝐿a_{H,L}italic_a start_POSTSUBSCRIPT italic_H , italic_L end_POSTSUBSCRIPT, bH,Lsubscript𝑏𝐻𝐿b_{H,L}italic_b start_POSTSUBSCRIPT italic_H , italic_L end_POSTSUBSCRIPT and cH,Lsubscript𝑐𝐻𝐿c_{H,L}italic_c start_POSTSUBSCRIPT italic_H , italic_L end_POSTSUBSCRIPT that specify the EoS. This is in stark contrast with supercooled flows, for which there is always a choice of these parameters for which the flow behind the wall only explores stable states. Instead, in the superheated case, at arbitrarily late times there is an arbitrarily large, approximately homogeneous region behind the wall that is thermodynamically disfavoured with respect to the L𝐿Litalic_L-phase at the same temperature. This region will eventually decay by nucleating bubbles of the L𝐿Litalic_L-phase if the bubble is allowed to expand for arbitrarily long times. Whether this happens in practice will depend on the nucleation rates in front and behind the bubble wall, on the bubble wall velocities, etc. Nevertheless, from the viewpoint of the dynamics of a single bubble, flows with negative pressure may regarded as thermodynamically unstable. In Fig. 7 these flows correspond to points above the blue curve labelled as p=0subscript𝑝0p_{-}=0italic_p start_POSTSUBSCRIPT - end_POSTSUBSCRIPT = 0. Note that deflagrations cannot suffer from this problem since the phase right behind the wall is directly the stable H𝐻Hitalic_H-phase.

We have also constructed detonations with ξw<cssubscript𝜉𝑤subscript𝑐𝑠\xi_{w}<c_{s}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT < italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT. In these, the wall connects the outside v=0𝑣0v=0italic_v = 0 state with some state below the blue line in Fig. 3, then the flow extends until the blue line and connects with the inside v=0𝑣0v=0italic_v = 0 state through a shock. We checked that the entropy production in these solutions can be positive. These solutions could in principle compete with the deflagrations and hybrids. However, they may be unstable and decay into the hybrids if perturbed, in analogy to what happens in the supercooled case with deflagrations and hybrids for ξw>cssubscript𝜉𝑤subscript𝑐𝑠\xi_{w}>c_{s}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT > italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT Espinosa:2010hh . We leave this stability analysis for future work.

All the flows above can be labeled by three independent parameters, that is, Tnsubscript𝑇𝑛T_{n}italic_T start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT, μnsubscript𝜇𝑛\mu_{n}italic_μ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT and ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT or, equivalently, ensubscript𝑒𝑛e_{n}italic_e start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT, nnsubscript𝑛𝑛n_{n}italic_n start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT and ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT. In the current case in which p=cs2e+constant𝑝superscriptsubscript𝑐𝑠2𝑒constantp=c_{s}^{2}e\,+\,\text{constant}italic_p = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_e + constant, one can solve for the enthalpy (or energy density) and the fluid velocity first, and then obtain the density for a given choice of nnsubscript𝑛𝑛n_{n}italic_n start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT. Therefore, for each choice of (en,ξw)subscript𝑒𝑛subscript𝜉𝑤(e_{n},\xi_{w})( italic_e start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT , italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ) we can extract a whole family of analogous solutions parameterized by nnsubscript𝑛𝑛n_{n}italic_n start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT. As in the case of supercooled bubbles, it is useful to characterize this choice in terms of the so-called “strength factor” αnsubscript𝛼𝑛\alpha_{n}italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT. Its general definition is111We have inverted the order in the numerator with respect to the supercooled case to ensure that it is positive. (see e.g. Hindmarsh:2020hop ; Ares:2020lbt )

αn13(eH3pH)(eL3pL)wL|(T,μ)=(Tn,μn).subscript𝛼𝑛evaluated-at13subscript𝑒𝐻3subscript𝑝𝐻subscript𝑒𝐿3subscript𝑝𝐿subscript𝑤𝐿𝑇𝜇subscript𝑇𝑛subscript𝜇𝑛\alpha_{n}\equiv\left.\frac{1}{3}\frac{(e_{H}-3p_{H})-(e_{L}-3p_{L})}{w_{L}}% \right|_{(T,\mu)=(T_{n},\mu_{n})}.italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ≡ divide start_ARG 1 end_ARG start_ARG 3 end_ARG divide start_ARG ( italic_e start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT - 3 italic_p start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT ) - ( italic_e start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT - 3 italic_p start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ) end_ARG start_ARG italic_w start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT end_ARG | start_POSTSUBSCRIPT ( italic_T , italic_μ ) = ( italic_T start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT , italic_μ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) end_POSTSUBSCRIPT . (15)

In the present case this reduces to

αn=ϵeL(Tn,μn)=ϵen,subscript𝛼𝑛italic-ϵsubscript𝑒𝐿subscript𝑇𝑛subscript𝜇𝑛italic-ϵsubscript𝑒𝑛\alpha_{n}=\frac{\epsilon}{e_{L}(T_{n},\mu_{n})}=\frac{\epsilon}{e_{n}}\,,italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT = divide start_ARG italic_ϵ end_ARG start_ARG italic_e start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( italic_T start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT , italic_μ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) end_ARG = divide start_ARG italic_ϵ end_ARG start_ARG italic_e start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT end_ARG , (16)

so choosing αnsubscript𝛼𝑛\alpha_{n}italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT is equivalent to choosing the nucleation energy density in units of ϵitalic-ϵ\epsilonitalic_ϵ. However, for general EoS both the strength factor and the energy density must be simultaneously specified in order to characterize a bubble solution.

The solutions presented in Figs. 5 and 6 depend on the values of ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT and αnsubscript𝛼𝑛\alpha_{n}italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT, but are independent on the values of the constants {aH,bH,cH,aL,bL,cL}subscript𝑎𝐻subscript𝑏𝐻subscript𝑐𝐻subscript𝑎𝐿subscript𝑏𝐿subscript𝑐𝐿\{a_{H},b_{H},c_{H},a_{L},b_{L},c_{L}\}{ italic_a start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT , italic_b start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT , italic_a start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT , italic_b start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT } in (14). However, the entropy production depends on these constants in (14), and the solutions will produce positive or negative entropy depending on the values of these constants. Below we comment further on this point.

4 Bounds

We will now see that the presence of a forbidden region in Fig. 3, together with the existence of a lower bound on the energy density in each branch, imposes bounds on the possible values of the parameters (ξw,αn)subscript𝜉𝑤subscript𝛼𝑛(\xi_{w},\alpha_{n})( italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT , italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) that characterize each possible fluid flow.

Let us start with detonations. Given a wall velocity ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT, consider increasing the strength αnsubscript𝛼𝑛\alpha_{n}italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT or, equivalently, decreasing the nucleation energy ensubscript𝑒𝑛e_{n}italic_e start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT relative to the bag constant ϵitalic-ϵ\epsilonitalic_ϵ. This decreases the energy behind the wall esubscript𝑒e_{-}italic_e start_POSTSUBSCRIPT - end_POSTSUBSCRIPT and increases the speed vsubscript𝑣v_{-}italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT. Clearly, the strongest transition we can consider is that for which the energy density in the fluid behind the wall reaches its vacuum energy, that is, e=ϵsubscript𝑒italic-ϵe_{-}=\epsilonitalic_e start_POSTSUBSCRIPT - end_POSTSUBSCRIPT = italic_ϵ. Precisely at this point we also reach the maximum possible velocity behind the wall, v=1subscript𝑣1v_{-}=1italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT = 1. Solving the matching conditions with these values results in the bound

αn,detξwcs21+ξw.subscript𝛼𝑛detsubscript𝜉𝑤superscriptsubscript𝑐𝑠21subscript𝜉𝑤\alpha_{n,\text{det}}\leq\frac{\xi_{w}-c_{s}^{2}}{1+\xi_{w}}.italic_α start_POSTSUBSCRIPT italic_n , det end_POSTSUBSCRIPT ≤ divide start_ARG italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT - italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 1 + italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT end_ARG . (17)

Detonations are therefore confined to the region simultaneously allowed by the conditions (17) and csξw1subscript𝑐𝑠subscript𝜉𝑤1c_{s}\leq\xi_{w}\leq 1italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ≤ italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ≤ 1. This region is displayed in pink in Fig. 7.

Refer to caption
Figure 7: Possible fluid flows for different strength factors αnsubscript𝛼𝑛\alpha_{n}italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT and wall velocities ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT. The bounds (17), (18) and (19) are represented as black solid curves. The flows allowed by the entropic bound (21) at zero chemical potential for different values of aL/aHsubscript𝑎𝐿subscript𝑎𝐻a_{L}/a_{H}italic_a start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT / italic_a start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT are those below the corresponding discontinuous black curves. The flows for which the pressure inside the bubble is lower (higher) than the pressure outside the bubble are those below (above) the red curve labelled as pin=pnsubscript𝑝insubscript𝑝𝑛p_{\textrm{in}}=p_{n}italic_p start_POSTSUBSCRIPT in end_POSTSUBSCRIPT = italic_p start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT. The flows for which the fluid pressure becomes negative at some point behind the wall are those above the blue curve labelled as p=0subscript𝑝0p_{-}=0italic_p start_POSTSUBSCRIPT - end_POSTSUBSCRIPT = 0.

In the case of deflagrations, a similar upper bound can be obtained for the strength, although it cannot be expressed in terms of a simple analytical formula. Fig. 3 shows that, given a ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT, the speed ahead of the wall can at most reach the boundary of the excluded region. This happens when the speed of the fluid in front of the wall moves at the speed of sound with respect to the wall, v+=cssubscript𝑣subscript𝑐𝑠v_{+}=c_{s}italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT. Notice that this is precisely the value at which hybrid solutions begin to exist. The matching conditions at the wall then lead to the following bound for the state inside the bubble:

α(csξw)2cs22csξw+1.subscript𝛼superscriptsubscript𝑐𝑠subscript𝜉𝑤2superscriptsubscript𝑐𝑠22subscript𝑐𝑠subscript𝜉𝑤1\alpha_{-}\leq\frac{(c_{s}-\xi_{w})^{2}}{c_{s}^{2}-2c_{s}\xi_{w}+1}.italic_α start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ≤ divide start_ARG ( italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT - italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 2 italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT + 1 end_ARG . (18)

Translating this bound to one on αnsubscript𝛼𝑛\alpha_{n}italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT can only be done numerically, and it requires solving for the entire fluid flow for deflagrations. The result can be seen in Fig. 7, where deflagrations are confined to the gray shaded region.

Finally, we can obtain a bound on hybrid solutions. Once again, we see in Fig. 3 that, given a ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT, the flow behind the wall has to fall below the excluded region. This puts the following constraint on the velocity behind the wall:

1v(ξw)cs2ξwξw(1cs2).1𝑣subscript𝜉𝑤superscriptsubscript𝑐𝑠2subscript𝜉𝑤subscript𝜉𝑤1superscriptsubscript𝑐𝑠2-1\leq v(\xi_{w})\leq-\frac{c_{s}^{2}-\xi_{w}}{\xi_{w}(1-c_{s}^{2})}.- 1 ≤ italic_v ( italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ) ≤ - divide start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT end_ARG start_ARG italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ( 1 - italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_ARG . (19)

The upper bound, which depends on the speed of sound, means that the wall must move faster than the shock, i.e. the flow must lie below the shaded region in Fig. 3. By solving for the whole flow, these two bounds translate into the two solid black lines that confine the green region in Fig. 7. The lower (upper) bound in (19) corresponds to the upper (lower) bound on ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT in Fig. 7. Notice that the bound on deflagrations happens to precisely match the lower bound on hybrids and therefore these two types of flows do not coexist. The matching of the two bounds could be expected as deflagrations are bounded by the flows that fulfill v+=cssubscript𝑣subscript𝑐𝑠v_{+}=c_{s}italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT, which is one of the defining conditions for hybrid flows. Thus deflagrations and hybrids do not compete in the same region of parameters.

We see that, as we increase the strength factor, the first solutions that cease to exist are the hybrid flows. The next ones to be forbidden are detonations. Therefore, for strong transitions only deflagrations survive, but they are confined to regions with very low wall velocities. This picture is opposite to that in the supercooled case (Espinosa:2010hh, ), where deflagrations cease to exist for strong transitions while hybrids and detonations become confined to higher and higher wall velocities.

In Fig. 7 we have added a long-dashed red curve separating the flows for which the pressure inside the bubble is higher (flows above the curve) or lower (flows below the curve) than the pressure outside the bubble. We see that, for a given wall velocity, these two cases occur for phase transitions with high or low strength factor, respectively. Motivated by these findings, we have explored the possible orderings of the pressures for a supercooled phase transition, where the standard picture is that the inside pressure is always higher than the outside pressure. We have constructed a large family of supercooled flows and we have found no counterexample to this picture.

The bounds above are model-independent, in the sense that they follow only from the allowed regions of Fig. 3 and from the existence of a minimum energy density in each phase. As we now discuss, there are two additional bounds that depend on some details of the model, meaning on the numerical coefficients aH,Lsubscript𝑎𝐻𝐿a_{H,L}italic_a start_POSTSUBSCRIPT italic_H , italic_L end_POSTSUBSCRIPT, bH,Lsubscript𝑏𝐻𝐿b_{H,L}italic_b start_POSTSUBSCRIPT italic_H , italic_L end_POSTSUBSCRIPT and cH,Lsubscript𝑐𝐻𝐿c_{H,L}italic_c start_POSTSUBSCRIPT italic_H , italic_L end_POSTSUBSCRIPT that specify the EoS.

The first bound follows from the existence of an upper limit on αnsubscript𝛼𝑛\alpha_{n}italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT given by the minimum energy density possible for nucleation, since nucleation can only happen in the stable region of the low-energy phase. This minimum will be attained at some point on the curve of first-order phase transitions, so we have

αnϵMin(eL(Tc,μc)).subscript𝛼𝑛italic-ϵMinsubscript𝑒𝐿subscript𝑇𝑐subscript𝜇𝑐\alpha_{n}\leq\frac{\epsilon}{\mathrm{Min}(e_{L}(T_{c},\mu_{c}))}.italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ≤ divide start_ARG italic_ϵ end_ARG start_ARG roman_Min ( italic_e start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( italic_T start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT , italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ) ) end_ARG . (20)

This bound would correspond to a horizontal line in Fig. 7. We do not show it in the figure since its position would depend on all the numerical coefficients in the EoS, but we note that it may compete with the other bounds presented in this section.

The second bound follows from entropic considerations. The flows presented above are described by ideal hydrodynamics, and therefore they do not produce entropy except at the discontinuities. At each discontinuity the entropy production must be non-negative. By Stoke’s theorem this amounts to the requirement that the entropy flow is non-negative across the discontinuity, which results in the condition

γ+v+s+γvs0.subscript𝛾subscript𝑣subscript𝑠subscript𝛾subscript𝑣subscript𝑠0\gamma_{+}v_{+}s_{+}-\gamma_{-}v_{-}s_{-}\geq 0\,.italic_γ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT + end_POSTSUBSCRIPT - italic_γ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ≥ 0 . (21)

The entropy densities on either side of the discontinuity, s±subscript𝑠plus-or-minuss_{\pm}italic_s start_POSTSUBSCRIPT ± end_POSTSUBSCRIPT, depend on the numerical coefficients that specify the EoS. The flows allowed by this bound for μ=0𝜇0\mu=0italic_μ = 0 and aL/aH={0.5,0.75,0.9}subscript𝑎𝐿subscript𝑎𝐻0.50.750.9a_{L}/a_{H}=\{0.5,0.75,0.9\}italic_a start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT / italic_a start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = { 0.5 , 0.75 , 0.9 } are those below the corresponding discontinuous curves in Fig. 7. We have verified that, for aL/aH=0.5subscript𝑎𝐿subscript𝑎𝐻0.5a_{L}/a_{H}=0.5italic_a start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT / italic_a start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = 0.5, the entropy production is positive for the flows in Figs. 5 and 6. In all the cases that we have examined we have found that the most restrictive bound is the entropic bound. It would be interesting to investigate if this is true in general.

Notice that flows that produce negative entropy are not physically meaningful as expanding, superheated bubbles. However, their time reversal gives rise to physically meaningful flows with positive entropy production that can be interpreted as contracting “drops” that were studied in Rezzolla:1995kv ; Rezzolla:2013dea .

5 Efficiency factor

Ref. Casalderrey-Solana:2022rrn showed that the dynamics of superheated bubbles in NS mergers can give rise to an interesting GW signal in the MHz range. A relevant quantity to estimate the amount of GW production is the so called efficiency factor Hindmarsh:2020hop

κ=1ϵ(4π/3)(ξwt)3wγ2v2d3x=3ϵξw3wγ2v2ξ2𝑑ξ.𝜅1italic-ϵ4𝜋3superscriptsubscript𝜉𝑤𝑡3𝑤superscript𝛾2superscript𝑣2superscript𝑑3𝑥3italic-ϵsuperscriptsubscript𝜉𝑤3𝑤superscript𝛾2superscript𝑣2superscript𝜉2differential-d𝜉\kappa=\frac{1}{\epsilon(4\pi/3)(\xi_{w}t)^{3}}\int w\gamma^{2}v^{2}d^{3}x=% \frac{3}{\epsilon\xi_{w}^{3}}\int w\gamma^{2}v^{2}\xi^{2}d\xi\,.italic_κ = divide start_ARG 1 end_ARG start_ARG italic_ϵ ( 4 italic_π / 3 ) ( italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT italic_t ) start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG ∫ italic_w italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_v start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_x = divide start_ARG 3 end_ARG start_ARG italic_ϵ italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG ∫ italic_w italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_v start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d italic_ξ . (22)

which is a measure of how much energy has been deposited in the motion of the fluid. We have computed the efficiency factor for several strength factors as a function of the wall velocity in the allowed regions where we can find bubble solutions. The result is summarized in Fig. 8. The empty regions in this figure arise form the disallowed regions, above the solid black curves, in Fig. 7.

Refer to caption
Figure 8: Efficiency factor (22) as a function of the bubble wall velocity for several choices of the strength factor.

The behavior of the efficiency factor for deflagrations and detonations is quite simple, since it is monotonous in both cases. For deflagrations the efficiency factor increases with the wall velocity at fixed strength, whereas for detonations it decreases. In both cases the efficiency factor increases with the strength factor at fixed wall velocity. Hybrids exhibit a richer, non-monotonic behaviour. For low αnsubscript𝛼𝑛\alpha_{n}italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT, hybrids connect the deflagrations with the detonations and the efficiency factor has a maximum at some intermediate value of the wall velocity. For moderate strengths, the maximum efficiency factor increases with the strength factor. As the strength is increased, the low-velocity hybrids cease to exist and the efficiency factor becomes a monotonically decreasing function of the wall velocity. At this stage, increasing the strength implies discarding hybrid solutions and the maximum efficiency factor decreases. In summary, the maximum efficiency factor corresponds to a hybrid at some intermediate values of both the wall velocity and the strength factor.

Fig. 8 includes all possible flows allowed by hydrodynamics. However, it must be kept in mind that there are more restrictive, model-dependent bounds that will exclude some of these flows, which in turn may reduce the possible maximum value of the efficiency factor.

6 Large jump in the number of degrees of freedom

We now turn our attention to the particular case in which the jump in the number of degrees of freedom between the two phases is large, that is, eLeHmuch-less-thansubscript𝑒𝐿subscript𝑒𝐻e_{L}\ll e_{H}italic_e start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ≪ italic_e start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT. In the supercooled case, this analysis led to an estimate for the wall velocity as a function of the nucleation temperature Sanchez-Garitaonandia:2023zqz . We will obtain an analogous result in the superheated case.

A convenient way to parameterize a large jump in the number of degrees of freedom is to rescale the L𝐿Litalic_L-phase pressure as pLxpLsubscript𝑝𝐿𝑥subscript𝑝𝐿p_{L}\to x\,p_{L}italic_p start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT → italic_x italic_p start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT, with x1much-less-than𝑥1x\ll 1italic_x ≪ 1. This immediately implies that all extensive quantities in the L𝐿Litalic_L-phase scale as x𝑥xitalic_x. In turn, this means that the strength factor diverges as

αn=ϵeLx1.subscript𝛼𝑛italic-ϵsubscript𝑒𝐿similar-tosuperscript𝑥1\alpha_{n}=\frac{\epsilon}{e_{L}}\sim x^{-1}\,.italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT = divide start_ARG italic_ϵ end_ARG start_ARG italic_e start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT end_ARG ∼ italic_x start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT . (23)

In view of Fig. 7, this immediately tells us that the flow must be a deflagration with very small wall velocity, since this is the only type of flow that exists for arbitrarily large strength factors. We have verified that the black solid curve on the left of Fig. 7 at small ξwsubscript𝜉𝑤\xi_{w}italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT is given by

ξwαn1.similar-tosubscript𝜉𝑤superscriptsubscript𝛼𝑛1\xi_{w}\sim\alpha_{n}^{-1}\,.italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ∼ italic_α start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT . (24)

Comparing the two equations above we conclude that the velocity must scale at small x𝑥xitalic_x as

ξwx.similar-tosubscript𝜉𝑤𝑥\xi_{w}\sim x\,.italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ∼ italic_x . (25)

This can also be seen from the following argument. In the limit x0𝑥0x\to 0italic_x → 0 the critical temperature and chemical potential approach finite, x𝑥xitalic_x-independent values Tc,μcsubscript𝑇𝑐subscript𝜇𝑐T_{c},\mu_{c}italic_T start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT , italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT. At temperatures and chemical potentials that are not parametrically separated from these the ratios of the following quantities in front and behind the wall obey

e+es+sw+wx.similar-tosubscript𝑒subscript𝑒subscript𝑠subscript𝑠similar-tosubscript𝑤subscript𝑤similar-to𝑥\frac{e_{+}}{e_{-}}\sim\frac{s_{+}}{s_{-}}\sim\frac{w_{+}}{w_{-}}\sim x\,.divide start_ARG italic_e start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG italic_e start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG ∼ divide start_ARG italic_s start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG italic_s start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG ∼ divide start_ARG italic_w start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG italic_w start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG ∼ italic_x . (26)

The pressure jump works slightly differently. By definition, pH(Tc,μc)=pL(Tc,μc)subscript𝑝𝐻subscript𝑇𝑐subscript𝜇𝑐subscript𝑝𝐿subscript𝑇𝑐subscript𝜇𝑐p_{H}(T_{c},\mu_{c})=p_{L}(T_{c},\mu_{c})italic_p start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT ( italic_T start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT , italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ) = italic_p start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( italic_T start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT , italic_μ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ). This means that near the critical curve not only pLsubscript𝑝𝐿p_{L}italic_p start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT but also pHsubscript𝑝𝐻p_{H}italic_p start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT are 𝒪(x)𝒪𝑥\mathcal{O}(x)caligraphic_O ( italic_x ), and therefore so is their difference. It then follows that

v+v=p+pex,v+v=ee++px1.formulae-sequencesubscript𝑣subscript𝑣subscript𝑝subscript𝑝subscript𝑒similar-to𝑥subscript𝑣subscript𝑣subscript𝑒subscript𝑒subscript𝑝similar-tosuperscript𝑥1v_{+}v_{-}=\frac{p_{+}-p_{-}}{e_{-}}\sim x\,,\quad\frac{v_{+}}{v_{-}}=\frac{e_% {-}}{e_{+}+p_{-}}\sim x^{-1}\,.italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT = divide start_ARG italic_p start_POSTSUBSCRIPT + end_POSTSUBSCRIPT - italic_p start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG italic_e start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG ∼ italic_x , divide start_ARG italic_v start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG = divide start_ARG italic_e start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG italic_e start_POSTSUBSCRIPT + end_POSTSUBSCRIPT + italic_p start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG ∼ italic_x start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT . (27)

From this we conclude that fluid velocity behind the wall in the rest frame of the wall is vxsimilar-tosubscript𝑣𝑥v_{-}\sim xitalic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ∼ italic_x. This is related to the velocity in the lab frame, v𝑣vitalic_v, through

v=vξw1vξw.subscript𝑣𝑣subscript𝜉𝑤1𝑣subscript𝜉𝑤v_{-}=\frac{v-\xi_{w}}{1-v\xi_{w}}\,.italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT = divide start_ARG italic_v - italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT end_ARG start_ARG 1 - italic_v italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT end_ARG . (28)

For deflagrations v=0𝑣0v=0italic_v = 0 behind the wall and hence

ξw|v|x,similar-tosubscript𝜉𝑤subscript𝑣similar-to𝑥\xi_{w}\sim|v_{-}|\sim x\,,italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ∼ | italic_v start_POSTSUBSCRIPT - end_POSTSUBSCRIPT | ∼ italic_x , (29)

in agreement with Eq. (25).

7 Discussion

Motivated by their possible role in NS mergers, in this paper we have studied the hydrodynamics of relativistic superheated bubbles in the presence of a conserved charge akin to baryon number. Equations (8) determine the possible flow profiles. In general, these equations are coupled to one another and one must solve for the velocity and the charge density simultaneously. Here we have observed that this problem simplifies if the speed of sound is constant in each branch of the EoS. In this case one can solve sequentially for v,w𝑣𝑤v,witalic_v , italic_w and n𝑛nitalic_n, meaning that the velocity profile is unaffected by the presence of the charge. For simplicity, in this paper we have restricted ourselves to this case by choosing a bag-model EoS for which cs2=1/3superscriptsubscript𝑐𝑠213c_{s}^{2}=1/3italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = 1 / 3 in both branches. One tantalizing result is that, in the case of a large jump in the number of degrees of freedom between the two phases, only deflagrations are possible and the bubble wall velocity is related to the strength of the transition through (24). Extrapolating to the hadronic-quark matter transition in QCD this would predict a value for the wall velocity around ξw0.1similar-tosubscript𝜉𝑤0.1\xi_{w}\sim 0.1italic_ξ start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ∼ 0.1, in agreement with the ballpark value suggested in Casalderrey-Solana:2022rrn . It would be interesting to verify these estimates with a more realistic approximation to the QCD EoS. In this case the velocity and the charge profiles will not decouple from one another and the fluid flows could be qualitatively different.

We have found two qualitative differences between superheated and supercooled flows. The first one is that some superheated detonations and hybrids can develop a necessarily metastable region behind the bubble wall, which at sufficiently long times will decay by nucleating bubbles of the L𝐿Litalic_L-phase. The second difference is that the pressure inside an expanding superheated bubble can be higher or lower than the pressure outside it. Even in the latter case, the bubble still expands by accreting energy from the outside. This raises an interesting dynamical question about how this steady state is reached from the slightly over-critical bubble that is initially nucleated. The reason is that, for the critical bubble, the inside pressure must be slightly higher than the outside pressure, since it must exactly compensate for the effect of the latter plus the surface tension. It would be interesting to see the time evolution of the inside pressure along the first instants after nucleation in a microscopic model, following the analysis in the supercooled case Bea:2021zsu . An additional insight provided by such a model would be a prediction of the bubble wall velocity from the microscopic dynamics. We will report on these aspects elsewhere.

Acknowledgements

We thank Alexandre Serantes for useful discussions. This work was supported by the “Center of Excellence Maria de Maeztu 2020-2023” award to the ICCUB (CEX2019-000918-M) funded by MCIN/AEI/10.13039/501100011033. YB, DM and JCS acknowledge support from grants PID2019- 105614GB-C21, PID2019-105614GB-C22, PID2022-136224NB-C21, PID2022-136224NB-C22, and 2021-SGR-872. The work of YB is also funded by a “Beatriu de Pinós” postdoctoral program under the Ministry of Research and Universities of the Government of Catalonia (2022 BP 00225) and by a “Maria Zambrano” postdoctoral fellowship from the University of Barcelona. The work of MSG is supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No758759).

Appendix A Non-relativistic, superheated bubbles with an incompressible phase

Superheated bubbles occur in many everyday situations, such as the formation of vapor bubbles in superheated water. In these cases, the matter surrounding the bubble walls is non-relativistic. Additionally, for the relevant range of pressures, the matter outside the superheated bubbles behaves as an incompressible fluid. These two features make the dynamics of those superheated bubbles very different from their relativistic analogues, as described in this work. Specifically, in the non-relativistic case, the bubble size does not grow at a constant rate; instead, the growth rate decreases with time t1/2superscript𝑡12t^{-1/2}italic_t start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT. For completeness, this appendix reviews the dynamics of non-relativistic superheated bubbles, and the discussion here follows closely the review vbblreview

We assume a non-relativistic system in its low-temperature phase that has been superheated past its critical temperature. In this system, characterized by a mass density ρLsubscript𝜌𝐿\rho_{L}italic_ρ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT, the low-temperature phase is incompressible, meaning ρLsubscript𝜌𝐿\rho_{L}italic_ρ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT remains constant. When a bubble of the thermodynamically stable high-temperature phase nucleates, it grows due to entropic preference. This phase has its own mass density, ρHsubscript𝜌𝐻\rho_{H}italic_ρ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT. However, similar to vapor bubbles in superheated water, we do not assume the fluid inside the bubble is incompressible. As with the relativistic case, once the bubble grows sufficiently large, the dynamics inside and outside the bubble can be described by non-relativistic hydrodynamics. The interface between the two phases can then be approximated as a discontinuity in the hydrodynamic fields.

Paralleling the discussion of the relativistic case, we begin by characterizing the fluid velocity in the vicinity of the bubble. This is analogous to the velocity profiles presented in fig. (3). In the low-temperature phase, the mass continuity equation combined with the incompressibility condition implies that the velocity field is divergenceless, i.e., 𝐯=0𝐯0{\bf\nabla v}=0∇ bold_v = 0. Under spherical symmetry, the velocity field outside the bubble is radial and can be expressed as:

𝐯(r,t)=F(t)r2r^rR(t)formulae-sequence𝐯𝑟𝑡𝐹𝑡superscript𝑟2^𝑟𝑟𝑅𝑡{\bf v}(r,t)=\frac{F(t)}{r^{2}}\hat{r}\,\quad\quad r\geq R(t)bold_v ( italic_r , italic_t ) = divide start_ARG italic_F ( italic_t ) end_ARG start_ARG italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG over^ start_ARG italic_r end_ARG italic_r ≥ italic_R ( italic_t ) (30)

where R(t)𝑅𝑡R(t)italic_R ( italic_t ) is the bubble radius and F(t)𝐹𝑡F(t)italic_F ( italic_t ) is a function determined by mass conservation, analogous to the charge conservation equation in non-relativistic systems. The velocity field vanishes far from the bubble, ensuring that the fluid at a distance remains at rest.

Since the bubble wall does not accumulate mass, the mass flux on both sides of the wall must be equal. This is equivalent to the last matching condition in Eq. (12).Imposing that the matter inside the bubble is at rest, and considering that as the bubble grows, the change in density between the high and low-temperature phases must be compensated by the influx of mass from the low-temperature phase outside the wall:

(ρLρH)R˙(t)=ρLv(R(t),t)F(t)=R2(t)R˙(t)(1ρHρL)subscript𝜌𝐿subscript𝜌𝐻˙𝑅𝑡subscript𝜌𝐿𝑣𝑅𝑡𝑡𝐹𝑡superscript𝑅2𝑡˙𝑅𝑡1subscript𝜌𝐻subscript𝜌𝐿\left(\rho_{L}-\rho_{H}\right)\dot{R}(t)=\rho_{L}v(R(t),t)\Rightarrow F(t)=R^{% 2}(t)\dot{R}(t)\left(1-\frac{\rho_{H}}{\rho_{L}}\right)( italic_ρ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT ) over˙ start_ARG italic_R end_ARG ( italic_t ) = italic_ρ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT italic_v ( italic_R ( italic_t ) , italic_t ) ⇒ italic_F ( italic_t ) = italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_t ) over˙ start_ARG italic_R end_ARG ( italic_t ) ( 1 - divide start_ARG italic_ρ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT end_ARG ) (31)

This expression for the velocity field, combined with the Navier-Stokes (NS) equation for incompressible fluids,

(pρL)=𝐯t+(𝐯)𝐯ν2𝐯𝑝subscript𝜌𝐿𝐯𝑡𝐯𝐯𝜈superscript2𝐯-{\bf\nabla}\left(\frac{p}{\rho_{L}}\right)=\frac{\partial\bf{v}}{\partial t}+% \left({\bf v}\cdot{\bf\nabla}\right){\bf v}-\nu\nabla^{2}{\bf v}- ∇ ( divide start_ARG italic_p end_ARG start_ARG italic_ρ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT end_ARG ) = divide start_ARG ∂ bold_v end_ARG start_ARG ∂ italic_t end_ARG + ( bold_v ⋅ ∇ ) bold_v - italic_ν ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_v (32)

where ν=η/ρL𝜈𝜂subscript𝜌𝐿\nu=\eta/\rho_{L}italic_ν = italic_η / italic_ρ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT is the kinematic viscosity, determines the form of the pressure field in front of the bubble. By substituting Eqs. (30) and (31) into the NS equation and integrating in the radial direction, we obtain the pressure difference in the low-temperature phase both far away from the bubble, psubscript𝑝p_{\infty}italic_p start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT and next to the bubble wall pL(R)subscript𝑝𝐿𝑅p_{L}(R)italic_p start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( italic_R ):

pL(R)pρF=(1Q)(12(3+Q)R˙(t)2+R(t)R¨(t))subscript𝑝𝐿𝑅subscript𝑝subscript𝜌𝐹1𝑄123𝑄˙𝑅superscript𝑡2𝑅𝑡¨𝑅𝑡\frac{p_{L}(R)-p_{\infty}}{\rho_{F}}=(1-Q)\left(\frac{1}{2}(3+Q)\dot{R}(t)^{2}% +R(t)\ddot{R}(t)\right)divide start_ARG italic_p start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( italic_R ) - italic_p start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT end_ARG = ( 1 - italic_Q ) ( divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( 3 + italic_Q ) over˙ start_ARG italic_R end_ARG ( italic_t ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_R ( italic_t ) over¨ start_ARG italic_R end_ARG ( italic_t ) ) (33)

with Q=ρH/ρLQ=\equiv\rho_{H}/\rho_{L}italic_Q = ≡ italic_ρ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT / italic_ρ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT. The pressure outside the bubble can be related to the pressure inside the bubble through the non-relativistic analog of the matching conditions for the momentum flux across the wall, given by the first equation in Eq. (12). Taking into account the effect of the bubble’s surface tension, γ𝛾\gammaitalic_γ, the fluid stress on both sides of the wall is given by:

pL(R)=pH2γR4(1Q)νRR˙subscript𝑝𝐿𝑅subscript𝑝𝐻2𝛾𝑅41𝑄𝜈𝑅˙𝑅p_{L}(R)=p_{H}-2\frac{\gamma}{R}-4(1-Q)\frac{\nu}{R}\dot{R}italic_p start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT ( italic_R ) = italic_p start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT - 2 divide start_ARG italic_γ end_ARG start_ARG italic_R end_ARG - 4 ( 1 - italic_Q ) divide start_ARG italic_ν end_ARG start_ARG italic_R end_ARG over˙ start_ARG italic_R end_ARG (34)

where the last term comes from the sheer stress contribution of the velocity field, Eq. (31). For very large bubbles, the effect of the surface tension and the viscosity becomes negligible; this is the reason why we have not considered those contribution in the relativistic analysis. Nevertheless, keeping those terms, and together with Eq. (33), we obtain the so called Rayleigh-Plesset (RP) equation

pHpρL=(1Q)(12(3+Q)R˙(t)2+R(t)R¨(t)+4νRR˙)+2γRsubscript𝑝𝐻subscript𝑝subscript𝜌𝐿1𝑄123𝑄˙𝑅superscript𝑡2𝑅𝑡¨𝑅𝑡4𝜈𝑅˙𝑅2𝛾𝑅\frac{p_{H}-p_{\infty}}{\rho_{L}}=(1-Q)\left(\frac{1}{2}(3+Q)\dot{R}(t)^{2}+R(% t)\ddot{R}(t)+4\frac{\nu}{R}\dot{R}\right)+2\frac{\gamma}{R}divide start_ARG italic_p start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT - italic_p start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT end_ARG = ( 1 - italic_Q ) ( divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( 3 + italic_Q ) over˙ start_ARG italic_R end_ARG ( italic_t ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_R ( italic_t ) over¨ start_ARG italic_R end_ARG ( italic_t ) + 4 divide start_ARG italic_ν end_ARG start_ARG italic_R end_ARG over˙ start_ARG italic_R end_ARG ) + 2 divide start_ARG italic_γ end_ARG start_ARG italic_R end_ARG (35)

which governs the dynamics of non-relativistic bubbles in incompressible fluids.

To complete the analysis, we also must consider the effect of the energy flux across the bubble wall, which in the relativistic case is given by the middle equation in Eq. (12). In general, in the non-relativistic limit, the energy flux is given by

=ρ(12v2+w¯)𝐯η𝐯σκT𝜌12superscript𝑣2¯𝑤𝐯𝜂𝐯𝜎𝜅𝑇{\bf\mathcal{F}}=\rho\left(\frac{1}{2}v^{2}+{\bar{w}}\right){\bf v}-\eta{\bf v% }\cdot\sigma-\kappa{\bf\nabla}Tcaligraphic_F = italic_ρ ( divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_v start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + over¯ start_ARG italic_w end_ARG ) bold_v - italic_η bold_v ⋅ italic_σ - italic_κ ∇ italic_T (36)

with w¯=e¯+pρ¯𝑤¯𝑒𝑝𝜌\bar{w}=\bar{e}+\frac{p}{\rho}over¯ start_ARG italic_w end_ARG = over¯ start_ARG italic_e end_ARG + divide start_ARG italic_p end_ARG start_ARG italic_ρ end_ARG and e¯¯𝑒{\bar{e}}over¯ start_ARG italic_e end_ARG are the enthalpy and energy per unit mass, σ𝜎\sigmaitalic_σ is the sheer tensor and κ𝜅\kappaitalic_κ the heat conductivity. For spherical bubbles and using Eqs. (30)and (31), the radial component of the flux right outside of the bubble wall is given by

outsider=12(1Q)3ρR˙3+ρw¯(1Q)R˙4η(1Q)2R˙2RκTrsuperscriptsubscriptoutside𝑟12superscript1𝑄3𝜌superscript˙𝑅3𝜌¯𝑤1𝑄˙𝑅4𝜂superscript1𝑄2superscript˙𝑅2𝑅𝜅𝑇𝑟\mathcal{F}_{\rm outside}^{r}=\frac{1}{2}(1-Q)^{3}\rho\dot{R}^{3}+\rho\bar{w}(% 1-Q)\dot{R}-4\eta(1-Q)^{2}\frac{\dot{R}^{2}}{R}-\kappa\frac{\partial T}{% \partial r}caligraphic_F start_POSTSUBSCRIPT roman_outside end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( 1 - italic_Q ) start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_ρ over˙ start_ARG italic_R end_ARG start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT + italic_ρ over¯ start_ARG italic_w end_ARG ( 1 - italic_Q ) over˙ start_ARG italic_R end_ARG - 4 italic_η ( 1 - italic_Q ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT divide start_ARG over˙ start_ARG italic_R end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_R end_ARG - italic_κ divide start_ARG ∂ italic_T end_ARG start_ARG ∂ italic_r end_ARG (37)

This expression shows that non-relativistic bubbles cannot grow at a constant rate, as their relativistic analogues do, at least for normal fluids, such that the density of the fluid phase is larger than the density of the vapour phase, Q<1𝑄1Q<1italic_Q < 1. For a constant and positive R˙˙𝑅\dot{R}over˙ start_ARG italic_R end_ARG, the first to terms of the expression above dominate and the energy flux is positive, which implies that the fluid motion outside the wall takes away energy from the interior of the bubble. This can also inferred from the combination of Eqs. (30)and (31), which show that the velocity of the fluid outside of the bubble is positive, unlike the relativistic case. Therefore, the only way that the latent heat can be compensated is by reducing the internal energy of the fluid outside; when this energy arrives to the critical energy, the pressure difference in the inside and outside of the bubble become equal and the fluid velocity must vanish, as inferred from Eq. (35). Therefore, superheated bubbles cannot expand at constant velocities for arbitrarily long times.

For bubbles to continue growing at asymptotic late times, energy must be transferred toward the interior of the bubble. This is achieved by heat transfer, represented by the last term in Eq. (37), which becomes the dominant mechanism at these late times. In this regime, integrating over the bubble’s surface, the energy balance between the heat flux and the latent heat can be expressed as follows:

Lddt(43πR2ρH)=4πR3κTr|r=R(t)𝐿𝑑𝑑𝑡43𝜋superscript𝑅2subscript𝜌𝐻evaluated-at4𝜋superscript𝑅3𝜅𝑇𝑟𝑟𝑅𝑡L\frac{d}{dt}\left(\frac{4}{3}\pi R^{2}\rho_{H}\right)=4\pi R^{3}\kappa\left.% \frac{\partial T}{\partial r}\right|_{r=R(t)}italic_L divide start_ARG italic_d end_ARG start_ARG italic_d italic_t end_ARG ( divide start_ARG 4 end_ARG start_ARG 3 end_ARG italic_π italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ρ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT ) = 4 italic_π italic_R start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_κ divide start_ARG ∂ italic_T end_ARG start_ARG ∂ italic_r end_ARG | start_POSTSUBSCRIPT italic_r = italic_R ( italic_t ) end_POSTSUBSCRIPT (38)

with L𝐿Litalic_L the latent heat of the transition. In this limit, the temperature gradient outside the bubble is primarily governed by heat diffusion, yielding TrΔT/(Dt)1/2similar-to𝑇𝑟Δ𝑇superscript𝐷𝑡12\frac{\partial T}{\partial r}\sim\Delta T/(Dt)^{1/2}divide start_ARG ∂ italic_T end_ARG start_ARG ∂ italic_r end_ARG ∼ roman_Δ italic_T / ( italic_D italic_t ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT with ΔTΔ𝑇\Delta Troman_Δ italic_T, where ΔTΔ𝑇\Delta Troman_Δ italic_T is the temperature difference between the wall and the ambient, and D𝐷Ditalic_D is the thermal diffusivity of the fluid outside the wall. Consequently, in this regime, the bubble radius R𝑅Ritalic_R scales as Rt1/2similar-to𝑅superscript𝑡12R\sim t^{1/2}italic_R ∼ italic_t start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT. For a comprehensive derivation of this expression, refer to vbblreview and the references therein.

References