Welcome to CosmoGW documentation!
CosmoGW is a Python library for the study of cosmological gravitational wave (GW) backgrounds from different sources in the early Universe, focusing on the production of GWs from sound waves and MHD turbulence from cosmological phase transitions.
The repository is stored on GitHub: https://github.com/CosmoGW/cosmoGW.
It is also publicly available on PyPI: https://pypi.org/project/cosmoGW/
To install the code, you can use pip:
pip install cosmoGW
Note
This project is under active development.
For any comments, suggestions, contributions, or collaborations feel free to contact me at alberto.roperpol@unige.ch
For a full project overview, installation instructions, and usage examples, see the README file:
Below is the table of contents for the documentation
Contents
- CosmoGW v1.0
- GW_analytical
- GW_back
- GW_models
- cosmology
- GW_templates
- hydro_bubbles
- Author
- Dates
- Contributors
- References
- Comments
Chapman_Jouget()Rstar_beta()compute_alphan()compute_def()compute_det()compute_hyb()compute_profiles_vws()compute_profiles_vws_multalp()compute_w()compute_xi_from_v()fp_z()kappas_Esp()kappas_from_prof()type_nucleation()v_shock()vplus_vminus()w_shock()w_to_lam()
- interferometry
- utils
- analysis
- Tutorial: LISA and gamma-ray telescopes as multi-messenger probes of a first-order cosmological phase transition
- Tutorial. LISA and :math:`gamma`-ray telescopes as multi-messenger probes of a first-order cosmological phase transition
- Tutorial that shows how to use CosmoGW to compute the GW spectrum from sound waves and turbulence in a first-order phase transition.
- It follows the description presented in [RoperPol:2023bqa]: “LISA and γ-ray telescopes as multi-messenger probes of a first-order cosmological phase transition,” A. Roper Pol, A. Neronov, C. Caprini, T. Boyer and D. Semikoz, arXiv:2307.10744 (2023)
- Created by Alberto Roper Pol (alberto.roperpol@unige.ch), 23/07/2025
- Comprehensive documentation of CosmoGW is available at Read the Docs.
- We compute the GW background using the models presented in RoperPol:2023bqa for sound waves and turbulence, available in CosmoGW by the functions OmGW_spec_sw and OmGW_spec_turb
- For illustration, let’s generate the plot of the GW spectra for a particular phase transition, described by the parameters \(\alpha = 0.5\), \(\beta/H_\ast = 10\), \(v_w =0.95\), \(\varepsilon_{\rm turb} = 1\), \(T_\ast = 100\) GeV
- Tutorial. Gravitational wave production from phase transitions: sound waves and decaying compressional motion
- Tutorial that shows how to use CosmoGW to compute the GW production according to the model of locally stationary UETC validated by Higgsless simulations, which extends the sound-shell model used for weak phase transitions to strong phase transitions where non-linearities develop and the kinetic energy density decays in time
- It follows the description presented in [Caprini:2024gyk]: “Gravitational waves from first-order phase transitions: from weak to strong transitions,” C. Caprini, R. Jinno, T. Konstandin, H. Rubira, A. Roper Pol, I. Stomberg, JHEP 07 (2025) 217, arXiv:2409.03651
- The results are summarized and presented in [Stomberg:2025kxf]: “Gravitational wave spectra for cosmological phase transitions with non-linear decay of the fluid motion,” I. Stomberg and A. Roper Pol, Gravitation session of the 59th Rencontres de Moriond proceedings (2025), arXiv:2508.04263
- In general, this tutorial also shows how to produce the GW background from phase transitions considering previous templates in the literature.
- Created by Alberto Roper Pol (alberto.roperpol@unige.ch), 16/04/2025
- Comprehensive documentation of CosmoGW is available at Read the Docs.
- The GW spectrum computed from sound waves (compressional motion) can be in general decomposed in the following way:
- where the prefactor pref depends on the PT parameters: length scale \(R_\ast\) and ratio of kinetic to total energy density \(K = \rho_{\rm kin}/\rho_{\rm total}\).
- In the first place, we directly generate the GW spectrum obtained from the Higgsless simulations of Caprini:2024gyk as presented in Stomberg:2025, where we use the numerical values as interpolated from the results of the numerical simulations. We compare the resulting GW
- In the following of this notebook, we will describe each of the terms that come into play to determine the GW spectrum and will review the different templates developed in the recent years in the literature, showing how CosmoGW can be used to reproduce each case and generate GW spectrum from sound waves and, in general, from compressional motion accounting for the decay of the source.
- The prefactor in the GW production is generalized in Caprini:2024gyk based on the locally stationary UETC, such that it can be expressed as
- where \(\kappa\) is the efficiency factor in the bag equation of state, which can be computed using kappas_Esp function in hydro_bubbles.py.
- for an equation of state \(p = c_{\rm s}^2 \rho\) where \(c_{\rm s}^2 = 1/3\) is assumed to be a constant.
- Alternatively, we can consider the value of \(K(\tau)\) at the end of the phase transition \(K_0 = K(\tau_0)\) obtained from the numerical simulations of Caprini:2024gyk.
- Similarly, Caprini:2024hue uses \(K = 0.6 K_\xi\) based on the simulations of Jinno:2022mie. To incorporate a constant factor (like 0.6), the GW spectrum can just be multiplied by the squared of this factor, so this case does not require a particular option in the different functions of CosmoGW.
- We can compare the resulting \(K_0\) from both methods. Both are available to compute the prefactor of the GW production using the function OmGW_spec_sw in GW_templates.py with the options model_K0 = ‘Espinosa’ or ‘higgsless’
- We then compute the factor \(K_{\rm exp, int}^2\) assuming a power-law in time \(K(\tau) = K_0 (\tau/\tau_0)^{-b}\), where \(b\) is estimated based on the numerical results of Caprini:2024gyk.
- These results can be used calling the function pref_GWB_sw with model = ‘sound_waves’, or directly calling the function OmGW_spec_sws with the option model_decay = ‘sound_waves’. To take into account expansion of the Universe, these functions should be called with the argument expansion = True
- Let us check the evolution with the source duration of the GW amplitude for different values of \(\beta/H_\ast = (10, 100, 1000)\) and compare to the shock formation time for the chosen values of \(\alpha = 0.0046, 0.05\) and \(0.5\) when there is no decay of the source
- Alternatively, we can take into account the decay of the source calling the function pref_GWB_sw with the option model = ‘decaying’ and providing the value of the decay exponent \(b\), or directly calling the function OmGW_spec_sws with the option model_decay = ‘decaying’.
- We can also use the values of \(b\) interpolated from the numerical results of Caprini:2024gyk, where we note that for some of the phase transition parameters, especially for weak \(\alpha\), the values of \(b\) might deviate from the physical ones due to underresolution.
- Taking into account these values of \(b\) for the decay, we can again reevaluate the evolution with the source duration of the GW amplitude and compare to the shock formation time for the chosen values of \(\alpha = 0.0046, 0.05\) and \(0.5\) including the effect from the decay.
- Hence, the decay of the source will in general lead to smaller amplitudes of the GW spectrum for the same duration of the source. However, since this model applies after non-linearities develop, it allows to extend the duration to times longer than that of the shock formation time.
- We can use the function K2int provided in the GW_models.py routine of CosmoGW. We show in dots the time of shock formation for each phase transition parameter point \((\alpha, v_w)\). Note that the simulations of Caprini:2024gyk only last until \(\delta \tau_{\rm sim} \beta = 22\) so any results extrapolated beyond that might not be correct.
- We can now include the effect of Universe expansion for different values of \(\beta/H_\ast\)
- With these results, we can estimate what is the effect of the decay in the amplitude of the GW spectrum in an expanding background. On one hand, the simulations of Caprini:2024gyk evolve the system also beyond the time of shock formation (at least for large \(\alpha\)). Therefore, the amplitude will continue growing beyond the shock formation time following the non-linear dynamics. For how long we can extrapolate these results beyond the ending time of the simulations is an open problem as we expect non-linearities and turbulence to eventually overcome this production.
- In the first place, let’s compare the amplitude at the time of shock formation. We can see that the inclusion of a decay always leads to a smaller amplitude of the GW spectrum. Below, we show the ratio of the \(\Upsilon_b\) factor, taking into account the decay, to the \(\Upsilon\) factor (with no decay). This ratio increases as \(\beta/H_\ast\) increases, as the shock formation time is \(\tau_{\rm sh} {\cal H}_\ast = (R_\ast \beta) (\beta/H_\ast)^{-1}/v_f\) and becomes smaller for larger \(b\).
- However, as time evolves, the amplitude keeps growing, and the simulations of Caprini:2024gyk show that the development of non-linearities does not lead instantly to a saturation of the amplitude. Hence, if the GW amplitude continues to grow beyond this shock formation time, eventually its amplitude can become larger to the one estimated when decay of the source is ignored and the GW amplitude saturates at one shock formation time.
- Below we show the results for a duration of the sourcing 50 times that of the shock formation time, taking into account the decay of the kinetic energy. When the decay is strong, the growth is not as large, as the amplitude of the GW spectrum reaches saturation faster.
- Finally, we focus now on the spectral shape, for which we consider a few different options.
- where \(s = 2 \pi f {\cal R}_\ast/10\) and \(S_0\) is a normalization constant such that \(\int S(s) {\rm d}\ln s = 1\).
- This spectral shape can be considered calling the function Sf_shape_sw in the GW_templates module or the function OmGW_spec_sw using the argument model_shape = sw_LISAold
- where \(s = f {\cal R}_\ast \Delta_w\) with \(\Delta_w = \frac{|v_w - c_{\rm s}|}{v_w}\) being the sound-shell thickness, which depends on the wall velocity, and \(m = (9 \Delta_w^4)/(\Delta_w^4 + 1)\), such that \(S/S_0 = 1\) at the peak \(s = 1\).
- This model can be chosen taking the argument model_shape = sw_SSM.
- defined such that \(S/S_0 = 1\) at \(s = 1\) and, again, \(s = {\cal R}_\ast f \Delta_w\).
- This model can be chosen taking the argument model_shape = sw_HL.
- Before going over the more recent templates, let us first compare the spectral shape based on these different models
- Let us now consider the recent template used by the LISA Cosmology working group (Caprini:2024hue) and the one developed in Caprini:2024gyk.
- where \(S_0\) is a normalization constant such that \(\int S(f) \, {\rm d} \ln f = 1\). This type of spectrum are available using the function smoothed_double_bPL in GW_analytical using the argument alpha2 = True (default is a slightly different convention used in other work).
- In Caprini:2024hue, the following values are considered: slopes \(a_1 = 3\), \(a_2 = 1\), \(a_3 = 3\) as in previous work, and smoothness parameters \(\alpha_1 = 2\) and \(\alpha_2 = 4\). The positions of the frequencies are such that \(f/f_1 = 5f {\cal R}_\ast\) and \(f/f_2 = 2 f {\cal R}_\ast \Delta_w\). However, the sound-shell thickness is in this case computed using \(\Delta_w = \xi_{\rm shell}/\max(v_w, c_{\rm s})\).
- This template is available in CosmoGW using the option model_shape = sw_LISA
- Note that previous templates (sw_SSM and sw_HL) can also be expressed in this generic form (available in smoothed_double_bPL) using \(a = 9\), \(b = 1\), \(c = 4\), \(\alpha_1 = 4\), \(\alpha_2 = 2\), \(f/f_1 = f {\cal R}_\ast\), and \(f/f_2 = \sqrt{m/(5 -m)} f {\cal R}_\ast \Delta_w\) (sw_SSM) and \(a = 3\), \(b = 1\), \(c = 3\), \(\alpha_1 = 1.5\), \(\alpha_2 = 0.5\), \(f/f_1 = f {\cal R}_\ast\), and \(f/f_2 = c_{\rm s} f {\cal R}_\ast \Delta_w\).
- 4. Caprini:2024gyk presents the same fit for a doubly broken power law and uses the results of the Higgsless numerical simulations to compute the spectral parameters. The simulations are ran for \(\alpha = 0.0046\), \(\alpha = 0.05\), and \(\alpha = 0.5\), and finds, after averaging over statistical realizations and wall velocities, the parameters \(\alpha_1 = 3.6\), \(\alpha_2 = 2.4\), \(f/f_1 = 2.5 f {\cal R}_\ast\), and \(f/f_2 = 2 f {\cal R}_\ast \Delta_w\) for weak phase transitions (0.0046), while no dependence on \(\Delta_w\) is found for intermediate (0.05) and strong (0.5) phase transitions, \(f/f_2 = f {\cal R}_\ast\) and \(f/f_2 = 2 f {\cal R}_\ast\), respectively.
- This template is available in CosmoGW using the option model_shape = sw_HLnew with the argument strength = ‘weak’, ‘interm’, or ‘strong’
- Both templates consider the dependence on the actual sound-shell thickness \(\Delta_w = \xi_{\rm shell}/\max(v_w, c_s)\), which depends on both the wall velocity and the strength of the phase transition. It can be computed using the hydro_bubbles module in CosmoGW. Below we show the spectral shape for a few cases.
- We can observe how the resulting spectral shape are rather compatible and the recent results from the simulations of Caprini:2024gyk justify the validity of the analysis in Caprini:2024hue. However, for intermediate and strong phase transitions, a second peak does not seem to depend any longer on \(\Delta_w\) in the Higgsless simulations and the ratio between \(f_1\) and \(f_2\) becomes constant, \(f_2/f_1 = 2.5\) for strong phase transitions and \(f_2/f_1 = 1.25\) for intermediate ones. These can be used calling the spectral shape function in CosmoGW using the arguments model_shape=’sw_HLnew’ and strength=’interm’ or ‘strong’
- One last option is available in CosmoGW, which corresponds to using the values of the spectral peak positions \(k_1\) and \(k_2\) interpolating the results from the Higgsless simulations of Caprini:2024gyk. Note that the position of the knee \(k_1\) is close to the IR wave number of the simulations and hence, its extraction is difficult. Similarly, for weak phase transitions, the value of \(k_2\) for confined hybrids (small sound-shell thickness) becomes unresolved in the UV for some cases. The interpolated values can be taken using the argument interpolate_HL = True when calling the spectral shape with model_shape = sw_HL_new.
- Let us have a look at the numerical results and their interpolated values, starting with \(k_1\), and compared to the averaged value found \(f_1 {\cal R}_\ast = 0.39\)
- The position of the second spectral peak \(k_2\) is now shown, comparing the numerical values found in Caprini:2024gyk, their averaged values used in the template shown above sw_HLnew, and the interpolated results for different values of \(\alpha\).
- We first compute \(\Delta_w\) for \(\alpha = 0.0046\) to compare with the numerical results
- One additional result found in Caprini:2024gyk is that for some intermediate and strong phase transitions vortical motion starts to develop. Therefore, the spectral slope at large frequencies can deviate with respect to the one from purely compressional turbulence (-3) and become shallower (as the one found in Kolmogorov turbulence, -8/3). Indeed, in some cases a shallower value of \(n_3\) is directly found in the GW spectrum. The current data from the simulations can also be interpolated. However, take these results with a grain of salt as they require further studies of the kinetic spectrum.
- To use these values of \(n_3\) one can call the spectral shape function with interpolate_HL = True and interpolate_HL_n3 = True
- We can now compare the spectral shape found from interpolating the spectral parameters with those using the averaged ones (see plots below)