Gravitational waves from bubble collisions and fluid motion in strongly supercooled phase transitions

We estimate the gravitational wave spectra generated in strongly supercooled phase transitions by bubble collisions and fluid motion. We derive analytically in the thin-wall approximation the efficiency factor that determines the share of the energy released in the transition between the scalar field and the fluid. We perform numerical simulations including the efficiency factor as a function of bubble radius separately for all points on the bubble surfaces to take into account their different collision times. We find that the efficiency factor does not significantly change the gravitational wave spectra and show that the result can be approximated by multiplying the spectrum obtained without the efficiency factor by its value at the radius Reff≃5/β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{\textrm{eff}} \simeq 5/\beta $$\end{document}, where β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} is the approximate inverse duration of the transition. We also provide updated fits for the gravitational wave spectra produced in strongly supercooled transitions from both bubble collisions and fluid motion depending on the behaviour of the sources after the collision.

We estimate the gravitational wave spectra generated in strongly supercooled phase transitions by bubble collisions and fluid motion.We derive analytically in the thin-wall approximation the efficiency factor that determines the share of the energy released in the transition between the scalar field and the fluid.We perform numerical simulations including the efficiency factor as a function of bubble radius separately for all points on the bubble surfaces to take into account their different collision times.We find that the efficiency factor does not significantly change the gravitational wave spectra and show that the result can be approximated by multiplying the spectrum obtained without the efficiency factor by its value at the radius R eff 5/β, where β is the approximate inverse duration of the transition.We also provide updated fits for the gravitational wave spectra produced in strongly supercooled transitions from both bubble collisions and fluid motion depending on the behaviour of the sources after the collision.

I. INTRODUCTION
The first observations of gravitational waves (GWs) by LIGO/Virgo signified the beginning of a new era in astrophysics and cosmology.While up to now all observed events were produced by compact object binaries [1][2][3][4], this new messenger brings hope also for detection of primordial signals in the form of stochastic GW backgrounds.Given the tremendous advancements in sensitivity that are expected throughout a broad frequency spectrum with the upcoming experiments [5][6][7][8][9][10][11][12][13][14][15], the prospects for probing the early Universe processes are great even though the compact object binaries that will contribute to the stochastic GW background make the detection of its primordial components more difficult [16].Interestingly, the recent pulsar timing observations [17][18][19][20] feature a common spectrum-process which could be an early indication of the upcoming first detection of a stochastic GW background, potentially of primordial origin [21][22][23][24][25][26][27][28][29][30][31][32][33].
Many high-energy processes, including phase transitions [34,35], cosmic strings [36] and inflation [37], occurring in the early Universe may generate a detectable stochastic GW background.In this paper we focus on cosmological first-order phase transitions featured in various particle physics models.They are intensive processes where bubbles of the new phase nucleate, expand and eventually convert the whole Universe in the true vacuum phase [38].Interactions between the expanding bubble walls and the surrounding fluid cause motion and inhomogeneities in the fluid, and both the collisions of the bubble walls and the motion of fluid inhomogeneities source GWs [39,40].The resulting GW spectra from these components have been extensively studied with numerical and semi-analytical methods (see e.g.[41][42][43][44][45][46][47] for recent progress).These studies indicate that different sources active during the transition can produce different GW spectra.
In order to determine the GW spectrum generated in a phase transition in a given particle physics model, we need to estimate how much each of the GW sources contributes to the final GW spectrum.The vacuum energy released in the transition is split between the gradient energy of the scalar field bubble wall and motion in the fluid.How the total released energy is split depends on the strength of the interactions between the wall and the particles in the fluid, and on the strength of the transition.
In strongly supercooled phase transitions it is possible that the interactions of the bubble wall with fluid do not stop the wall from accelerating before it collides with other bubbles.In this case most of the released energy is in the bubble walls and the bubble collisions give the dominant contribution to the GW spectrum.This can happen in particular in quasi-conformal models [48][49][50][51][52][53][54][55][56][57][58][59][60][61].If the interactions instead are sufficiently strong, the bubble wall reaches a terminal velocity before the collisions and majority of the released energy goes into fluid motion.This is the typical case in extensions of the Standard Model featuring polynomial scalar potentials [62][63][64][65][66][67][68][69][70][71][72].
In this paper we derive analytically an efficiency factor that determines how large is the contribution from each of the GW sources.We perform numerical simulations of the phase transition, describing both of the GW sources, bubble walls and fluid motion, in the thin-wall limit, to show how the efficiency factor affects the final GW spectrum.Moreover, we derive analytically the probability density function for the radius at which a given point on the surface of a bubble collides with another bubble and verify the results against our numerical simulations.Finally we also provide updated fits to the spectral shapes of the GW signals that can be produced by all sources active in very strong phase transitions.We estimate how the energy released in the bubble expansion is shared between the scalar field gradients and the fluid motion in strongly supercooled phase transitions by studying the bubble expansion under the influence of pressure terms caused by the interactions of the wall with the ambient fluid.We perform the computation consistently in the thin-wall limit, which gives a good description of the system if the bubble reaches ultra-relativistic velocities.The following analysis improves earlier approximations used in the literature [58,71,73].
In the thin-wall limit the energy carried by the bubble walls can be modeled using a simple analytical prescription.This assumes that the bubble walls are spherical shells with a certain surface energy density and the interactions of the walls with the ambient fluid are local.In this limit, neglecting the expansion of the Universe,1 the evolution of the bubble radius R can be described by the equation of motion [74] that arises from energy conservation of the coupled system of fluid and the scalar field bubble.Here ∆P ( Ṙ) denotes the pressure difference across the bubble wall and σ is the surface tension of the wall.The latter is defined through the scalar potential V as [38] σ where ϕ c > 0 denotes the field value at which the potential energy is the the same as in the false vacuum that lies at the origin, V (ϕ c ) = V (0).In terms of the Lorentz factor of the bubble wall, γ = 1/ 1 − Ṙ2 , the equation of motion (1) is given by The bubble nucleates at rest, γ = 1, dγ/dR = 0, with an initial radius R 0 .By Eq. ( 3) we can relate the wall tension to the initial radius as R 0 = 2σ/∆P 0 , where ∆P 0 ≡ ∆P (γ = 1).For a constant pressure difference, ∆P = ∆P 0 , the solution of the equation of motion is The total pressure difference across the bubble wall, accounting for 1 → 1 scatterings and 1 → N splittings at the bubble wall, is given by where ∆V denotes the potential energy difference between the minima.The pressure arising from 1 → 1 scatterings quickly reaches a constant value in the relativistic limit [74,75].Subsequently, the γ dependence of the total pressure difference arises only from ∆P 1→N , for which we consider two forms.The first one, suggested in [76][77][78] is linear in the Lorentz factor ∆P 1→N = ∆P 1→N γ, and the second one, suggested in [79,80], is quadratic in the Lorentz factor, ∆P 1→N = ∆P 1→N γ 2 In both cases ∆P 1→N is a constant.By plugging ∆P (γ) into (3), we can easily solve γ as a function of R. When ∆P 1→N ∆V − ∆P 1→1 the solution can be approximated by Eq. ( 4).Assuming in addition that R R 0 , the Lorentz factor grows linearly with the radius, γ ≈ 2R/(3R 0 ).Eventually, as the bubble wall accelerates, γ becomes large enough for the 1 → N splittings to be important, ∆P 1→N ∼ ∆V − ∆P 1→1 , after which it asymptotically reaches the value where c = 1, 2 depending on the scaling of the 1 → N pressure, ∆P 1→N ∝ γ c .The change from the linear growth to asymptotically constant behaviour occurs when the radius reaches The solution γ(R), in the limit γ eq 1, can be approximated by a simple broken power-law In Fig. 1 we show the full solution γ(R) for different values of γ eq for both ∆P 1→N ∝ γ (solid) and ∆P 1→N ∝ γ 2 (dashed) together with the above approximation shown by the dotted black curves.In both cases the transition from linear growth, γ ∝ R to the constant value γ ≈ γ eq is quite fast, and the difference between the two cases is small.The main effect of the scaling of ∆P 1→N is that it changes γ eq and R eq .We define the efficiency factor κ as the fraction of the total released energy within a unit solid angle that goes into the bubble wall energy, The rest of the released energy, 1 − κ(R), goes into fluid motion.This is a good approximation for strongly supercooled transitions.In weaker transitions one also needs  3), while the black dashed curves correspond to the approximations ( 8) and (10).
to keep in mind that some of the energy going into the fluid will be lost on its heating which will reduce the overall GW signal from the fluid motion [71,81].
Using the approximation (8), we can express the efficiency factor as where is a constant, K < 1.The parameters α and α ∞ are defined by scaling with the radiation energy density ρ R as α = ∆V /ρ R and α ∞ = ∆P 1→1 /ρ R (see [71] for more details).Typically for strongly supercooled transitions K ≈ 1.As shown in the right panel of Fig. 1, the efficiency factor remains constant at R R eq and decreases as κ ∝ 1/R at R R eq .In the same way as for γ(R), the difference between the cases ∆P 1→N ∝ γ and ∆P 1→N ∝ γ2 is small.

III. BUBBLE NUCLEATION AND COLLISIONS
Soon after the bubble has nucleated, we can neglect its initial radius, and, if the friction terms are sufficiently small (γ eq 1), we can approximate that the bubble radius grows as R = t − t n , where t n denotes the nucleation time of the bubble.Moreover, assuming that the bubbles are much smaller than the Hubble horizon, we can neglect the expansion of the Universe.The expected number of bubbles reaching a given point is then given by where Γ(t) denotes the bubble nucleation rate per unit time and volume, and the probability that the given point still is in the false vacuum at time t is Let us consider a bubble nucleation rate Γ(t) = Ce A(t) .Around the time t * when the transition proceeds, we can expand A(t) to get Γ(t) = Ce A(t * )+β(t−t * ) = Γ 0 e βt , where β ≡ d ln Γ/dt| t=t * and Γ 0 ≡ Ce A(t * )−βt * .As the transition is not an instantaneous process, the choice of t * includes some ambiguity.It is convenient to choose t * by requiring that P (t * ) = 1/e, which gives Γ 0 = β 4 /(8π), and N (t) = e βt .
Next, let us consider a point on the surface of a bubble that nucleated at time t n .If the point is still in the false vacuum when the radius of the bubble is R, then it has stayed in the false vacuum for the whole time 0 ≤ t − t n < R. The probability for this is P (t n + R).So, the probability that a bubble nucleated within time t n < t < t n + dt n in a volume V , and that a point on its surface is still in the false vacuum at radius R, is given by dt n V Γ(t n )P (t n + R).By integrating this over the nucleation time t n we get the probability density function for the radius at which a bubble surface element collides with the surface of another bubble, which we normalize to unity, dR c p(R c ) = 1.For the exponential bubble nucleation rate, Γ(t) ∝ e βt , this gives (independently of the prefactor Γ 0 The above result provides a good cross-check for the numerical simulations that we will use for the GW computation.In Fig. 2 the solid curve and the gray band indicate the mean and variance of the R c distribution obtained from 90 simulations with simulation volume (16/β) 3 and each of the simulation including at least 70 bubbles.In these simulations we nucleate thin-wall bubbles according to the exponential rate inside a cubic box with periodic boundary conditions, evolve them according to R = t − t n , discretise the bubble surfaces, and find the radius at which each of the points on the bubble surface collide with a wall of another bubble using the cosine rule.We label the bubbles with index j and denote the position vectors of the bubble centers by x j .Consider a point defined by the angles θ and φ on the surface of the bubble j = j .The radius at which that point collides with a surface of another bubble is given by where the minimum is taken over all bubbles, d 2 j ≡ | x j − x j | 2 is the distance between the bubble nucleation centers, ∆t j ≡ t n,j − t n,j is the time between their nucleation, and θ j is the angle between the vector x j − x j and the vector corresponding to the angles θ and φ.As shown in Fig. 2, the simulation result agrees well with the analytical result (15).
A widely used approximation for the bubble radius upon collision comes from the bubble number density n bubbles = dt n Γ(t n )P (t n ) which leads to R * = n −1/3 bubbles = (8π) 1/3 /β.From p(R c ) we can calculate moments of the bubble radius when a bubble surface element collides with the surface of another bubble, R n c = dR c R n c p(R c ).For the exponential bubble nucleation rate this gives Given that the released energy scales with the radius to the third power, this leads to a different estimate of the average bubble radius R 3 c 1/3 = 6 1/3 /β more appropriate for computation of the GW spectrum.

IV. GRAVITATIONAL WAVES
The energy released in the bubble expansion is divided between the bubble wall and the fluid shell that follows right behind the wall.Both the bubble walls and the fluid shells source GWs.We model these sources in the thinwall limit and calculate the GW spectrum accounting for the efficiency factor κ(R) for the bubble collisions and 1 − κ(R) for the fluid motion.The modeling of the fluid motion in the thin-wall limit is based on the assumption that the released energy going to fluid motion is strongly localized in a thin shell.Before collision this fluid shell is right behind the bubble wall and after the collision it propagates to the same direction as before the collision,however , depending on how strong the interaction are between the fluid and the scalar field, it's velocity can slow down to the speed of sound.
We calculate the GW spectrum as e.g. in Refs.[43,44] assuming that, as in the previous section, the bubble nucleation follows exponential rate per unit volume, Γ ∝ e βt .Each of the contributions (l =bubbles, fluid) to the resulting energy spectrum of GWs can be expressed as where ) encodes the spectral shape of the signal.The integral is over the wavevector k directions, and the integrand is ∝ V s /β 5 if the volume V s over which we average the GW energy spectrum is sufficiently big.
Using the thin-wall limit, the functions C l,+ and C l,× in the direction k = (0, 0, 1), can be expressed as The sum runs over all the bubbles nucleated in the volume V s , t n,j is the time of nucleation of the bubble j, z j is the z coordinate of its center, and R j = v l (t − t n,j ), where v l is the bubble wall/fluid shell velocity, denotes the radius of the bubble/fluid shell j at time t.For the bubble walls we use v bubbles = 1 both before and after the collision, whereas for the fluid shells we use v fluid = 1 before the collision and after the collision we consider two cases: v fluid = 1 and v fluid = c s = 1/ √ 3. The former is appropriate for very strong transitions [82], whereas the latter is realized for weaker transitions [45].The functions g +,× read g + (φ) = cos(2φ) and g × (φ) = sin(2φ).
The function f l (R) encodes the scaling of the GW source [44].For the bubble collisions contribution, we follow the results of lattice simulations [43,44], which showed that the maximum of the stress-energy tensor scales as T rr ∝ R −ξ after the collision.The power ξ in general depends on the underlying particle physics model.In particular, it was shown in [43] that breaking of a global symmetry corresponds to ξ = 2 while in [44] it was shown that in models where the phase transition breaks a gauge symmetry correspond to ξ = 3. Accounting also for the efficiency factor κ, the f l function for bubble collisions is given by where R c denotes the bubble radius at the moment of collision, t = t c .In contrast with Refs.[43,44], where R c was determined numerically by the bisection method, we find R c using Eq. ( 16).Also for the fluid motion we assume that the maximum of the stress-energy tensor scales as R −ξ after the collision.The function f l for fluid motion then reads In the perfect fluid description, that assumes the fluid to remain in local equilibrium at all times, the transversetraceless part of the stress energy tensor of the fluid reads T ij = γ 2 v i v j w, where v is the fluid velocity and w is its enthalpy density.By the interactions of the fluid with the wall, an overdense fluid shell with radial velocity v r > 0 builds up around the bubble wall.If the wall reaches a terminal velocity, the fluid shell settles into a self-similar profile [81].We expect that this shell continues to propagate after the bubble wall collides with the wall of another bubble without changing its shape significantly in the collision.This behaviour was first observed for weaker transitions, α < O(0.1), in lattice simulations [83,84].For our case of very strong transitions, α ≥ O(10), we used a simplified simulation involving only the fluid and assuming extra symmetry of the system to retain only one spatial dimension as in [45]. 3We begun with simulating collisions of two planar shells and verified that they are not significantly modified and instead simply propagate onward.We next simulated the evolution of spherically symmetric fluid shells after the collision.Fig. 3 shows an illustrative example of our results.We find that the maximum of the rr component of the stress energy tensor scales as T rr ∝ R −3 .This matches to the same scaling found in [45] for weak transitions and motivates us to consider ξ = 3 for the scaling of the fluid related GW source.For comparison, we consider also ξ = 2 which corresponds to the bulk flow model [86].
In principle, after the GW generation through relativistic fluid shells finishes, one would expect to enter the period of sound waves [41,83,84,87] and perhaps also turbulence [88][89][90][91].However, in the very strong transitions, that are our primary interest, we expect that the fluid will remain in the relativistic shells for a long time after the transition (at least until the shell radius has grown by O(1) factor).The energy carried by the inhomogeneities has then significantly diluted once the sound wave and turbulence periods begin and, therefore, we expect the main contribution on the GW spectrum in very strong transitions to arise from the relativistic fluid shells or the scalar field bubbles themselves.Thus, in the present analysis we neglect the periods of sound waves and turbulence.
For certain simple forms of the f l function the time integral in Eq. ( 20) can be done analytically, which makes the simulation significantly faster.In particular, it can be done analytically if f l is a broken power-law with integer powers.We consider the form (10) for the efficiency factor κ with c = 1.Strictly speaking our results then hold for the case that ∆P 1→N ∝ γ.However, since the difference in κ(R) for c = 1 and c = 2 is very small, our results give a good approximation also of the case ∆P 1→N ∝ γ 2 .The pressure ∆P 1→N mainly just determines the asymptotic radius R eq through Eqs. ( 7) and (6).In our simulations R eq is an input parameter, and we perform the numerical simulations for several values  I. Fitted values for the parametrization of the spectral shape (24) and fitted value of βR eff in Eq. ( 23).The corresponding spectra are shown in Fig 5.
of R eq .We also assume that K ≈ 1, which typically holds for strongly supercooled transitions, so that V. RESULTS We perform 90 simulations with simulation volume (16/β) 3 , each including at least 70 bubbles, for a range of R eq values including both signals due to bubble walls and the surrounding fluid in each of the cases described in the previous section.From the simulations we compute the spectral shape function (19).In each case we fit the data combined from the 90 simulations with a broken power-law spectrum of the form where a, b > 0 determine the low and high frequency power-law tails of the spectrum, c > 0 the width of the transition between these power-laws, while f p and A the peak frequency and amplitude of the spectrum respectively.The resulting GW spectra are shown in Fig. 4 with the solid curves.The color coding indicates different values of R eq .For the solid curves in Fig. 4 the efficiency factor is directly included in the simulation as in Eqs. ( 21) and (22).A commonly used approximation for the effect of the efficiency factor on the GW spectrum is to multiply the spectra obtained for bubble collisions and fluid motion without any efficiency factor by κ(R eff ) 2 and by (1 − κ(R eff )) 2 , respectively.To check this, we have computed the amplitude of the GW spectrum in each case relative to the R eq case that gives the largest amplitude and fitted R eff .The data points and resulting fits for all cases are shown in Fig. 5 and the corresponding fitted values of R eff in the last line of Table I.We find that the effect of the efficiency factor is almost independent of the behaviour of the GW source after the collisions.In all cases our results give R eff 5/β, showing that the often used approximation with R eff = (8π) 1/3 /β ≈ 2.9/β slightly underestimates R eff .Moreover, the results of applying the fitted efficiency factor are shown in Fig. 4 by the dashed curves.For these curves we have used the mean values given in Table I that are obtained by averaging over the fits with different equilibrium radius, except for the amplitude for which we use the strongest signal for each source.We see that the dashed curves agree very well with the fully numerical results shown by the solid curves.This shows that the efficiency factor does not change the shape of the GW spectrum but gives only an overall suppression factor.
To summarize, we have shown that for very strong transitions, α α ∞ , the GW spectrum from bubble collisions and from fluid motion, accounting for the distribution of energy between these sources, is given by where the efficiency factor is given by Eq. (23).The fitted values of the parameters A, a, b, c, f p and R eff are given in Table I.For weaker transitions, α α ∞ , also the prefactor K, given in Eq. ( 11), needs to be accounted for, as well as the suppression arising from heating of the fluid around the bubble wall [81].In the limit of large wall velocity appropriate for strong transitions this reduction takes a simple form [71] where The GW spectrum today can be obtained from ( 25) by accounting for the scaling of the amplitude and frequency  with the scale factor [43]:4 where h denotes the dimensionless Hubble constant, h = 0.674 [93], and h * the inverse Hubble time at the transition redshifted to its value today Here T * denotes the temperature after the transition (including reheating) and g * the effective number of relativistic degrees of freedom at that temperature.At scales larger than the horizon scale at the time of the transition the source is not coherent and, consequently, in standard radiation domination the slope of the spectrum changes to Ω GW ∝ f 3 for f < h * [94,95]. 5n this paper we have revisited the energy budget of strong first-order phase transitions to verify its impact on the produced gravitational wave spectra.We have gone beyond the current state-of-art by including the efficiency as a function of radius of the bubble accounting for the collision time of each point on the bubble surface.We have utilised numerical simulations randomly nucleating bubbles in a three dimensional box with periodic boundaries and used these to compute the GW spectra.This has allowed us to confirm that a simplified treatment of simply scaling entire spectra with an efficiency factor computed at some characteristic radius is accurate as the spectral shapes do not change due to the efficiency factor significantly.We did, however, find that in order to accurately describe the results the characteristic radius used in the simplified calculation should be around R eff ≈ 5/β rather than the usually employed average bubble separation R * = (8π) 1 3 /β ≈ 2.9/β.In each simulation we have also took into account the scaling of the GW sources after the collision in order to provide new fits for the resulting spectra from strongly supercooled transitions.The results are shown in Table I and, starting from strongest transitions, include bubble collision spectra for both T rr ∝ R −3 , appropriate for gauge symmetry breaking, and T rr ∝ R −2 , appropriate for global symmetry breaking.Going towards slightly weaker transitions, we have provided the spectrum generated by fluid motion with the scaling T rr ∝ R −3 and assuming the fluid remains in the form of relativistic shocks v fluid = 1 after the transition.For transitions which are not extremely strong, we have show results closer to the sound wave picture in which the velocity of the fluid quickly relaxes to the speed of sound v fluid = c s , again assuming the scaling T rr ∝ R −3 .Finally, for illustration, we also provide fluid spectra assuming the scaling T rr ∝ R −2 .
Taking into account that for very relativistic walls the fluid profiles are extremely peaked, we have thus show that the final GW spectrum will be indistinguishable from an even stronger transition where bubble collisions would be the main source.Only for weaker transitions where the hydrodynamical effects change the propagation speed of the fluid shells, the spectrum diverges from the spectrum arising from bubble collisions.

FIG. 1 .
FIG. 1.The Lorentz factor of the bubble wall (left panel) and the efficiency factor κ (right panel) as a function of the bubble radius R for γeq = 10, 100, 10 3 , 10 4 from light to dark.The solid black curves show the case ∆P1→N ∝ γ and the orange dashed curves the case ∆P1→N ∝ γ 2 .In the left panel the vertical dashed lines indicate the values of Req/R0.The colored curves show the full solution of the equation of motion (3), while the black dashed curves correspond to the approximations (8) and(10).

FIG. 3 .
FIG.3.Time evolution of the perfect fluid stress energy tensor (solid lines) together with the r −3 scaling for comparison (dashed line).The red profile is the initial condition just after the bubble collision tc and darker colours show the profile at later times.The parameters for this example profile are α = 20 and γw = 50 corresponding to a very strong transition such that the velocity of the profile remains nearly constant and only very slowly changes to the speed of sound.
FIG.4.Fitted GW spectral shape sourced by bubble walls and fluid motion assuming different scalings of the source after the collisions, Trr ∝ R −ξ , and different velocities of the fluid shell after the collision.The solid curves show the results obtained by directly including the factor κ(R) to the simulation and the dashed curves the result obtained by scaling the result obtained without it factor by κ(R eff ).

FIG. 5 .
FIG.5.The blue and red points with error bars show the amplitude of the GW spectrum from bubble collisions and from fluid motion, respectively, relative to the Req value that gives the largest amplitude.The blue and red curves show κ bubbles and κ fluid , respectively.The parameters read ξ = 3 and v fluid = 1 for both solid lines and the points of corresponding colour.Dashed and dotted lines and their corresponding points show the remaining cases (ξ = 2 and v fluid = cs) which as we see largely overlap with the previous two.
(15)2.Probability distribution function of bubble radius at collision.The red dashed curve shows the analytical result(15).The black curve and the gray band show the mean and the variance of the result from our numerical simulations of bubble nucleation.