hydro_bubbles

hydro_bubbles.py is a Python routine that contains functions to study the 1D hydrodynamic solutions of expanding bubbles produced from first-order phase transitions.

Currently part of the cosmoGW code:

https://github.com/cosmoGW/cosmoGW/ https://github.com/cosmoGW/cosmoGW/blob/main/src/cosmoGW/hydro_bubbles.py

Note

For full documentation, visit Read the Docs

To use it, first install cosmoGW:

pip install cosmoGW

Author

• Alberto Roper Pol (alberto.roperpol@unige.ch)

Dates

• Created: 01/02/2023

• Updated: 21/08/2025 (release cosmoGW 1.0: https://pypi.org/project/cosmoGW)

Contributors

• Antonino Midiri, Simona Procacci

References

• [RoperPol:2025a]: A. Roper Pol, S. Procacci, A. S. Midiri, C. Caprini, “Irrotational fluid perturbations from first-order phase transitions,” in preparation

• [Espinosa:2010hh]: J. R. Espinosa, T. Konstandin, J. M. No, G. Servant, “Energy Budget of Cosmological First-order Phase Transitions,” JCAP 06 (2010) 028, arXiv:1004.4187.

• [Hindmarsh:2016lnk]: M. Hindmarsh, “Sound shell model for acoustic gravitational wave production at a first-order phase transition in the early Universe,” Phys. Rev. Lett. 120 (2018) 7, 071301, arXiv:1608.04735.

• [Hindmarsh:2019phv]: M. Hindmarsh, M. Hijazi, “Gravitational waves from first order cosmological phase transitions in the Sound Shell Model,” JCAP 12 (2019) 062, arXiv:1909.10040.

• [RoperPol:2023dzg]: A. Roper Pol, S. Procacci, C. Caprini, “Characterization of the gravitational wave spectrum from sound waves within the sound shell model,” Phys. Rev. D 109, 063531 (2024), arXiv:2308.12943.

Comments

The solutions are based in the 1D hydrodynamic descriptions of Espinosa:2010hh and Hindmarsh:2019phv.

The details accompanying the code are provided in RoperPol:2023dzg and RoperPol:2025a (appendix A).

cosmoGW.hydro_bubbles.Chapman_Jouget(alp)

Compute the Chapman-Jouget wall velocity.

The Chapman-Jouget velocity is the wall velocity at which the relative speed behind the wall becomes that of the speed of sound. It separates detonations and supersonic deflagrations.

\[v_{\rm CJ} = \frac{1}{\sqrt{3}(1 + \alpha)} \left[1 + \sqrt{\alpha(2 + 3\alpha)}\right].\]

Reference: Eq. 39 of Espinosa:2010hh

Parameters

alpfloat

Alpha at the + side of the wall (\(\alpha_{+}\)).

Returns

vcJfloat

Chapman-Jouget speed \(v_{\rm CJ}\).

References

Espinosa:2010hh (eq. 39) and Hindmarsh:2019phv (eq. B19). Note that Eq. B19 of Hindmarsh:2019phv has a typo.

cosmoGW.hydro_bubbles.Rstar_beta(vws=1.0, cs2=0.3333333333333333, corr=True)

Compute the ratio of the mean bubble separation \(R_\ast\) to the inverse nucleation rate parameter \(\beta\).

\[R_\ast \beta = (8\pi)^{1/3} v_w\]

If corr is True, applies a correction for the speed of sound:

\[R_\ast \beta \to R_\ast \beta \times \max(1, c_s/v_w)\]

Parameters

vwsfloat or array_like

Bubble wall velocity \(v_w\).

cs2float, optional

Square of the speed of sound \(c_s^2\) (default is 1/3).

corrbool, optional

If True, apply correction for the speed of sound (default is True).

Returns

Rbetafloat or ndarray

Ratio \(R_* \beta\).

cosmoGW.hydro_bubbles.compute_alphan(vw=0.5, alpha_obj=0.263, tol=1e-05, cs2=0.3333333333333333, quiet=False, max_it=30, Nxi=10000, ty='def')

Iteratively compute the value of \(\alpha_+\) corresponding to a target \(\alpha\).

Uses the 1D profile of w and Newton-Raphson update to find \(\alpha_+\) that gives the correct \(\alpha\).

Parameters

vwfloat

Wall velocity.

alpha_objfloat

Target value of alpha at nucleation temperature.

tolfloat, optional

Relative tolerance for convergence (default is 1e-5).

cs2float, optional

Speed of sound squared (default is 1/3).

quietbool, optional

If True, suppress debug output.

max_itint, optional

Maximum number of iterations (default is 30).

Nxiint, optional

Number of discretization points in xi.

tystr, optional

Type of solution (‘def’ or ‘hyb’).

Returns

xis0ndarray

Array of xi.

vvs0ndarray

Array of velocities.

wws0ndarray

Array of enthalpies.

xi_shfloat

Position of the shock.

shbool

Boolean, True if a shock forms.

w_plfloat

Plus value of the enthalpy.

w_mfloat

Minus value of the enthalpy.

alpha_nfloat

Converged alpha at nucleation.

alp_plusfloat

Value of alpha_+ leading to alpha_obj.

convbool

True if the algorithm has converged.

cosmoGW.hydro_bubbles.compute_def(vw=0.5, alpha=0.263, cs2=0.3333333333333333, Nxi=10000, shock=True)

Compute the solutions for a subsonic deflagration 1D profile given vw and alpha.

Uses _def_sol() to compute the velocity and enthalpy profiles.

Parameters

vwfloat, optional

Wall velocity (default is 0.5).

alphafloat, optional

Strength of the phase transition (default is 0.263).

cs2float, optional

Speed of sound squared (default is 1/3).

Nxiint, optional

Number of discretization points in xi.

shockbool, optional

Stop calculation once a shock is formed (default is True).

Returns

xisndarray

Array of xi.

vsndarray

Array of 1D velocities.

wsndarray

Array of 1D enthalpies.

xi_shfloat

Position of the shock.

shbool

Boolean, True if a shock forms.

w_plfloat

Plus value of the enthalpy across the bubble wall.

w_mfloat

Minus value of the enthalpy across the bubble wall.

cosmoGW.hydro_bubbles.compute_det(vw=0.77, alpha=0.091, cs2=0.3333333333333333, Nxi=10000)

Compute the solutions for a detonation 1D profile given vw and alpha.

Uses _det_sol() to compute velocity and enthalpy profiles.

Parameters

vwfloat, optional

Wall velocity (default is 0.77).

alphafloat, optional

Strength of the phase transition (default is 0.091).

cs2float, optional

Speed of sound squared (default is 1/3).

Nxiint, optional

Number of discretization points in xi.

Returns

xisndarray

Array of xi.

vsndarray

Array of 1D velocities.

wsndarray

Array of 1D enthalpies.

xi_shfloat

Position of the shock (set to vw for detonations).

shbool

Boolean, always False for detonations.

w_plfloat

Plus value of the enthalpy across the bubble wall.

w_mfloat

Minus value of the enthalpy across the bubble wall.

cosmoGW.hydro_bubbles.compute_hyb(vw=0.7, alpha=0.052, cs2=0.3333333333333333, Nxi=10000, shock=True)

Compute the solutions for a supersonic deflagration 1D profile given vw and alpha.

Uses _det_sol() and _def_sol() to compute velocity and enthalpy profiles.

Parameters

vwfloat, optional

Wall velocity (default is 0.7).

alphafloat, optional

Strength of the phase transition (default is 0.052).

cs2float, optional

Speed of sound squared (default is 1/3).

Nxiint, optional

Number of discretization points in xi.

shockbool, optional

Stop calculation once a shock is formed (default is True).

Returns

xisndarray

Array of xi.

vsndarray

Array of 1D velocities.

wsndarray

Array of 1D enthalpies.

xi_shfloat

Position of the shock.

shbool

Boolean, True if a shock forms.

w_plfloat

Plus value of the enthalpy across the bubble wall.

w_mfloat

Minus value of the enthalpy across the bubble wall.

cosmoGW.hydro_bubbles.compute_profiles_vws(alpha, vws=None, cs2=0.3333333333333333, Nvws=20, Nxi=10000, Nxi2=10, plot=False, plot_v='v', cols=None, alphan=True, quiet=True, tol=1e-05, max_it=30, ls='solid', alp=1.0, lam=False, legs=False, fs_lg=14, st_lg=2, eff=False, save=False, dec_vw=1, ress='results/1d_profiles', strs_vws=None, str_alp=None)

Compute velocity and enthalpy profiles for a given alpha and a range of wall velocities.

This function solves the hydrodynamic equations for bubble expansion and returns the velocity and enthalpy (or energy density perturbation) profiles for each wall velocity.

Parameters

alphafloat

Nucleation temperature alpha.

vwsarray_like, optional

Array of wall velocities. If None, uses a default linspace.

cs2float, optional

Square of the speed of sound (default is 1/3).

Nvwsint, optional

Number of wall velocities if vws is not provided.

Nxiint, optional

Number of discretization points in xi.

Nxi2int, optional

Number of discretization points in xi out of the profiles.

plotbool, optional

If True, plot the resulting profiles.

plot_vstr, optional

What to plot: ‘v’, ‘w’, or ‘both’ (default ‘v’).

colslist, optional

List of colors to use for plotting.

alphanbool, optional

If True, input alpha is at nucleation temperature; if False, input is alpha+.

quietbool, optional

If True, suppress debug output.

tolfloat, optional

Tolerance for convergence of alpha+ (default is 1e-5).

max_itint, optional

Maximum number of iterations for alpha+ convergence.

lsstr, optional

Line style for plots.

alpfloat, optional

Opacity for plots.

lambool, optional

If True, compute energy perturbations lambda instead of enthalpy.

legsbool, optional

Whether to show legend in plots.

fs_lgint, optional

Font size for legend.

st_lgint, optional

Legend style.

effbool, optional

If True, also compute efficiency factors kappa and omega.

savebool, optional

If True, save results to CSV files.

dec_vwint, optional

Decimal precision for wall velocity in filenames.

ressstr, optional

Directory to save results.

strs_vwslist, optional

Custom strings for wall velocities in filenames.

str_alplist, optional

Custom strings for alpha in filenames.

Returns

xisndarray

Array of xi positions.

vvsndarray

Array of velocities.

wwsndarray

Array of enthalpies (or energy density perturbations if lam is True).

alphas_nndarray

Array of nucleation alphas (if input is alpha+, returns alpha+).

convndarray

Array of booleans for alpha+ convergence.

shocksndarray

Array of booleans for shock formation.

xi_shocksndarray

Positions of shocks (or xi_front if no shock).

wmsndarray

Enthalpy values at - side of the wall.

kappasndarray, optional

Ratio of kinetic to vacuum energy density (if eff is True).

omegasndarray, optional

Ratio of energy density perturbations to vacuum energy density (if eff is True).

cosmoGW.hydro_bubbles.compute_profiles_vws_multalp(alphas, vws=None, cs2=0.3333333333333333, Nvws=20, Nxi=10000, Nxi2=10, alphan=True, quiet=True, tol=1e-05, max_it=30, lam=False, eff=False)

Compute velocity and enthalpy profiles for an array of alpha and wall velocities.

This function wraps compute_profiles_vws() for multiple alpha values and wall velocities, optionally computing efficiency factors.

Parameters

alphasarray_like

Array of nucleation temperature alpha values.

vwsarray_like, optional

Array of wall velocities. If None, uses a default linspace.

cs2float, optional

Square of the speed of sound (default is 1/3).

Nvwsint, optional

Number of wall velocities if vws is not provided.

Nxiint, optional

Number of discretization points in xi.

Nxi2int, optional

Number of discretization points in xi out of the profiles.

alphanbool, optional

If True, input alpha is at nucleation temperature; if False, input is alpha+.

quietbool, optional

If True, suppress debug output.

tolfloat, optional

Tolerance for convergence of alpha+ (default is 1e-5).

max_itint, optional

Maximum number of iterations for alpha+ convergence.

lambool, optional

If True, compute energy perturbations lambda instead of enthalpy.

effbool, optional

If True, also compute efficiency factors kappa and omega.

Returns

xisndarray

Array of xi positions.

vvsndarray

Array of velocities.

wwsndarray

Array of enthalpies (or energy density perturbations if lam is True).

alphas_nndarray

Array of nucleation alphas (if input is alpha+, returns alpha+).

convndarray

Array of booleans for alpha+ convergence.

shocksndarray

Array of booleans for shock formation.

xi_shocksndarray

Positions of shocks (or xi_front if no shock).

wmsndarray

Enthalpy values at - side of the wall.

kappasndarray, optional

Ratio of kinetic to vacuum energy density (if eff is True).

omegasndarray, optional

Ratio of energy density perturbations to vacuum energy density (if eff is True).

cosmoGW.hydro_bubbles.compute_w(v, xi, cs2=0.3333333333333333)

Compute the enthalpy from the solution of the velocity profile.

Parameters

xiarray_like

Self-similar r/t.

varray_like

1D velocity profile.

cs2float, optional

Square of the speed of sound (default 1/3).

Returns

wndarray

Enthalpy profile.

cosmoGW.hydro_bubbles.compute_xi_from_v(v, xi0, cs2=0.3333333333333333, shock=False)

Compute the solution \(\xi(v)\) using a 4th-order Runge-Kutta scheme.

Since \(dv/d\xi\) has a singularity, it is necessary to compute \(\xi(v)\) and then invert each of the solutions to have the full dynamical solution. For the physical solution, computing \(v(\xi)\) is more practical.

Parameters

varray_like

Velocity array.

xi0float

Position where boundary condition is known (first value in velocity array).

cs2float, optional

Square of the speed of sound (default 1/3).

shockbool, optional

Option to stop the integration when a shock is found (as in deflagrations).

Returns

xindarray

Self-similar array.

shbool

Boolean determining if a shock is formed.

indshint

Index determining the position of the shock (-1 if no shock).

cosmoGW.hydro_bubbles.fp_z(xi, vs, z, lz=False, ls=None, multi=True, quiet=False)

Compute the functions \(f'(z)\) and \(l(z)\) for the Fourier transform of the velocity field.

These functions provide the initial conditions used to compute the kinetic spectrum in the sound-wave regime according to the Sound-Shell Model.

Parameters

xindarray

Array of xi positions.

vsndarray

Array of velocity profiles.

zndarray

Array of z = k (t - t_n) where f’ and l functions are to be computed.

lzbool, optional

If True, compute l(z) (default False).

lsndarray, optional

Array of energy density perturbations (used if lz is True).

multibool, optional

If True, use an array of wall velocities.

quietbool, optional

If True, suppress debug output.

Returns

fpzsndarray

Array of f’(z) values.

lzsndarray, optional

Array of l(z) values (if lz is True).

cosmoGW.hydro_bubbles.kappas_Esp(vw, alp, cs2=0.3333333333333333)

Compute the efficiency in converting vacuum to kinetic energy density.

Uses semiempirical fits from Espinosa:2010hh, appendix A, following the bag equation of state.

Numerical values can be computed from the 1d profiles using compute_profiles_vws() setting eff to True.

Parameters

vwfloat or array_like

Wall velocity.

alpfloat or array_like

Strength of the phase transition at nucleation temperature.

cs2float, optional

Square of the speed of sound (default is 1/3).

Returns

kappafloat

Ratio of kinetic to vacuum energy density, computed in the bag equation of state.

cosmoGW.hydro_bubbles.kappas_from_prof(vw, alpha, xis, ws, vs)

Compute the kinetic energy density efficiency kappa and thermal factor omega from 1D profiles.

\[\kappa = \frac{4}{v_w^3 \alpha} \int \xi^2 \frac{w}{1-v^2} v^2 d\xi \omega = \frac{3}{v_w^3 \alpha} \int \xi^2 (w-1) d\xi\]

Parameters

vwfloat

Wall velocity.

alphafloat

Strength of the phase transition.

xisndarray

Array of xi positions.

wsndarray

Array of enthalpy profiles.

vsndarray

Array of velocity profiles.

Returns

kappafloat

Ratio of kinetic to vacuum energy density.

omegafloat

Ratio of energy density perturbations to vacuum energy density.

cosmoGW.hydro_bubbles.type_nucleation(vw, alp, cs2=0.3333333333333333)

Determine the type of bubble solution.

Returns ‘def’ for subsonic deflagrations, ‘hyb’ for supersonic deflagrations, or ‘det’ for detonations.

Parameters

vwfloat or array_like

Bubble wall speed.

alpfloat or array_like

Strength of the phase transition at the + side of the wall (\(\alpha_{+}\)).

cs2float, optional

Square of the speed of sound (default 1/3).

Returns

tystr or ndarray

Type of solution (‘def’, ‘hyb’, ‘det’).

cosmoGW.hydro_bubbles.v_shock(xi)

Compute the shock velocity at the - side of the shock.

Computed from the \(v_+ v_- = 1/3\) condition.

Parameters

xifloat or array_like

Self-similar variable.

Returns

vshfloat or ndarray

Shock velocity.

cosmoGW.hydro_bubbles.vplus_vminus(alpha, vw=1.0, ty='det', cs2=0.3333333333333333)

Compute \(v_+\) and \(v_-\) (in the wall frame) for different types of solutions.

This function returns the boundary velocities for deflagrations, detonations, and hybrids, as required for hydrodynamic matching.

Parameters

alphafloat

Value of alpha at the + side of the wall.

vwfloat, optional

Wall velocity.

tystr, optional

Type of solution (‘def’, ‘det’, ‘hyb’).

cs2float, optional

Speed of sound squared (default is 1/3).

Returns

vplusfloat or ndarray

Velocity on the + side of the wall (in wall frame), \(v_+\).

vminusfloat or ndarray

Velocity on the - side of the wall (in wall frame), \(v_-\).

cosmoGW.hydro_bubbles.w_shock(xi)

Compute the ratio of enthalpies w-/w+ across the shock.

Computed from the \(v_+ v_- = 1/3\) condition.

Parameters

xifloat or array_like

Self-similar variable.

Returns

wshfloat or ndarray

Ratio of enthalpies across the shock.

cosmoGW.hydro_bubbles.w_to_lam(xis, ws, vw, alphan)

Compute the energy density perturbations from the enthalpy profile using the bag equation of state.

\[\lambda(\xi) = \frac{3}{4} (w(\xi) - 1) - \frac{3}{4} \alpha_n \text{ for } \xi < v_w\]

Parameters

xisndarray

Array of xi positions.

wsndarray

Array of enthalpy profiles.

vwfloat

Wall velocity.

alphanfloat

Value of nucleation alpha.

Returns

lamndarray

Array of energy density perturbations.

Reference

Appendix A of RoperPol:2025a.