GW_analytical

GW_analytical.py is a Python routine that contains analytical calculations and useful mathematical functions.

Adapted from the original GW_analytical in cosmoGW (https://github.com/AlbertoRoper/cosmoGW), created in Dec. 2021

Currently part of the cosmoGW code:

https://github.com/cosmoGW/cosmoGW/ https://github.com/cosmoGW/cosmoGW/blob/main/src/cosmoGW/GW_analytical.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/12/2021

• Updated: 31/08/2024

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

Contributors

• Antonino Midiri, Madeline Salomé

References

• [RoperPol:2022iel]: A. Roper Pol, C. Caprini, A. Neronov, D. Semikoz, “The gravitational wave signal from primordial magnetic fields in the Pulsar Timing Array frequency band,” Phys. Rev. D 105, 123502 (2022), arXiv:2201.05630

• [RoperPol:2023bqa]: A. Roper Pol, A. Neronov, C. Caprini, T. Boyer, D. Semikoz, “LISA and γ-ray telescopes as multi-messenger probes of a first-order cosmological phase transition,” arXiv:2307.10744 (2023)

• [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

• [Caprini:2024gyk]: A. Roper Pol, I. Stomberg, C. Caprini, R. Jinno, T. Konstandin, H. Rubira, “Gravitational waves from first-order phase transitions: from weak to strong,” JHEP 07 (2025) 217, arXiv:2409.03651.

• [Caprini:2024hue]: E. Madge, C. Caprini, R. Jinno, M. Lewicki, M. Merchand, G. Nardini, M. Pieroni, A. Roper Pol, V. Vaskonen, “Gravitational waves from first-order phase transitions in LISA: reconstruction pipeline and physics interpretation,” JCAP 10 (2024), 020, arXiv:2403.03723.

• [RoperPol:2025b]: A. Roper Pol, A. Midiri, M. Salomé, C. Caprini, “Modeling the gravitational wave spectrum from slowly decaying sources in the early Universe: constant-in-time and coherent-decay models,” in preparation

cosmoGW.GW_analytical.Iabn(a=5, b=0.6666666666666666, alp=2, n=0, norm=True, alpha2=False, piecewise=False)

Compute the n-th moment of the smoothed bPL spectra.

\[I_{ab; n} (\alpha) = \int K^n \zeta(K) dK = {\cal I}_{ab; n} (\alpha) \frac{\Gamma \bigl[\frac{a + n + 1}{\alpha (a + b)}\bigr] \Gamma\bigl[\frac{b - n - 1}{\alpha (a + b)}\bigr]} {\Gamma\bigl[\frac{1}{\alpha}\bigr]}\]

Parameters

a, bfloat

Slopes of the spectrum (default 5 and 2/3).

alpfloat

Smoothness of the transition \(\alpha\) (default 2).

nint

Moment of the integration (default 0).

normbool

Normalize the spectrum (default True).

alpha2bool

Use alternative convention (default False).

piecewisebool

Use piecewise broken power law (default False).

Returns

moment_valuefloat

Value of the n-th moment.

cosmoGW.GW_analytical.Kstarm(a=5, b=0.6666666666666666, alp=2, m=-1, norm=True, alpha2=False, piecewise=False)

Compute the location of the peak of the function

\[\zeta_m (K) = K^m \zeta(K),\]

where \(\zeta(K)\) is the original smoothed double broken power law (dbP).

Parameters

Same as zetam.

Returns

Kstarmfloat

Position of the peak of new function \(K^m \zeta(K)\).

cosmoGW.GW_analytical.calA(a=5, b=0.6666666666666666, alp=2, norm=True, alpha2=False, piecewise=False, dlogk=True)

Compute the parameter \({\cal A} = I_{ab;0}\) that relates the peak and the integrated values of the smoothed bPL spectrum.

References

RoperPol:2022iel, equation 8 (for dlogk False).

RoperPol:2025b, appendix A (for dlogk True).

Parameters

Same as Iabn function with n = -1 or 0 (if dlogk is True or False).

dlogkbool

Use logarithmic spacing (default True).

Returns

calAfloat

Value of the parameter \({\cal A}\).

cosmoGW.GW_analytical.calB(a=5, b=0.6666666666666666, alp=2, n=1, norm=True, alpha2=False, piecewise=False, dlogk=True)

Compute the parameter

\[{\cal B} = \left(\frac{I_{ab; -1 - n}}{I_{ab; -1}}\right)^{1/n}\]

that relates the peak and the integral scale (when \(n = 1\)), \({\cal B} = \xi\, k_{\text{peak}}\). When dlogk is False, the definition changes slightly, \(-1 \to 0\).

Parameters

Same as Iabn function.

nint

Moment of the integration (default 1).

dlogkbool

Use logarithmic spacing (default True).

Returns

calBfloat

Value of the parameter \({\cal B}\).

Reference

RoperPol:2025b, appendix A (for dlogk True)

cosmoGW.GW_analytical.calC(a=5, b=0.6666666666666666, alp=2, tp='vort', norm=True, alpha2=False, piecewise=False, dlogk=True)

Compute the parameter \(\cal C\) for the TT-projected stress spectrum.

It gives the value of the stresses spectrum at \(K \to 0\) limit as a prefactor pref that depends on the type of source and an integral over

\[ \begin{align}\begin{aligned}\int \frac{\zeta^2}{P^4} d\log k, \qquad {\rm when \ dlogk \ is \ True \ (see \ appendix \ B \ of \ RoperPol:2025b)}\\\int \frac{\zeta^2}{P^2} d K, \qquad {\rm when \ dlogk \ is \ False \ (RoperPol:2022iel \ for \ vortical \ and }\\\qquad \qquad {\rm eq. \ 46 \ of \ RoperPol:2023dzg \ for \ compressional \ fields)}\end{aligned}\end{align} \]

The spectrum of the stresses is then normalized such that (see appendix B of RoperPol:2025b):

\[E_\Pi (K) = K^3 E_\ast^2 \, {\cal C} \, \zeta_\Pi (K).\]

Parameters

a, bfloat

Slopes of the spectrum (default 5 and 2/3).

alpfloat

Smoothness of the transition \(\alpha\) (default 2).

tpstr

Type of sourcing field: ‘vort’, ‘comp’, ‘hel’, or ‘mix’.

normbool

Normalize the spectrum (default True).

alpha2bool

Use alternative convention (default False).

piecewisebool

Use piecewise broken power law (default False).

dlogkbool

Use per unit log(k) or per unit k (default True).

Returns

calCfloat

Value of the parameter \(\cal C\).

cosmoGW.GW_analytical.calIab_n_alpha(a=5, b=0.6666666666666666, alp=2, n=0, norm=True)

Compute the normalization factor used in the calculation of \(I_{ab; n}\) (see Iabn function),

\[{\cal I}_{ab; n} (\alpha) = \frac{1}{\alpha (a + b)} \left(\frac{a + b}{b} \right)^{1/\alpha} \left(\frac{a}{b}\right)^{-\frac{a + 1 + n}{\alpha (a + b)}}\]

Parameters

a, bfloat

Slopes of the smoothed_bPL function (default 5 and 2/3).

alpfloat

Smoothness parameter \(\alpha\) (default 2).

nint

n-moment of the integral (default 0).

normbool

Normalize the spectrum (default True).

Returns

calIfloat

Normalization parameter for the integral, \(\cal I_{ab; n} (\alpha)\).

Reference

Appendix A of RoperPol:2025b

cosmoGW.GW_analytical.smoothed_bPL(k, A=1.0, a=5, b=0.6666666666666666, kpeak=1.0, alp=2, norm=True, Omega=False, alpha2=False, piecewise=False, dlogk=True)

Return the value of the smoothed broken power law (bPL) model.

The spectrum is of the form:

\[\zeta(K) = A \cdot (b + \left| a \right|)^{1/\alpha} \frac{K^a}{\left[ b + c \cdot K^{\alpha(a + b)} \right]^{1/\alpha}}\]

where \(K = k/k_\ast\), \(c = 1\) if \(a = 0\) or \(c = \text{abs}(a)\) otherwise.

If norm is False:

\[\zeta(K) = A \cdot \frac{K^a}{ \left( 1 + K^{\alpha(a + b)} \right)^{1/\alpha}}\]

RoperPol:2022iel and RoperPol:2025b consider a spectral function \(\zeta(K)\) defined such that the average squared field corresponds to

\[ \begin{align}\begin{aligned}\langle v^2 \rangle \propto \Omega_v = A k_\ast \int \zeta(K) dK, \qquad {\rm in \ RoperPol:2022iel}\\\langle v^2 \rangle \propto \Omega_v = A \int \zeta(K) d \ln K, \qquad {\rm in \ RoperPol:2025b}\end{aligned}\end{align} \]

The first convention can be chosen in the following functions if dlogk is set to False, while the second one is assumed when dlogk is True.

An alternative parameterization is used when alpha2 is True:

\[\zeta(K) = A \cdot (b + \left| a \right|)^{(a + b)/\alpha} \frac{K^a}{\left[ b + c \cdot K^{\alpha} \right]^{(a + b)/\alpha}}\]

Parameters

karray_like

Array of wave numbers.

Afloat, optional

Amplitude of the spectrum (default 1).

afloat, optional

Slope at low wave numbers (default 5).

bfloat, optional

Slope at high wave numbers (default 2/3).

kpeakfloat, optional

Spectral peak position \(k_\ast\) (default 1).

alpfloat, optional

Smoothness of the transition \(\alpha\) (default 2).

normbool, optional

Normalize spectrum to peak at \(k_\ast\) with amplitude A (default True).

Omegabool, optional

Use integrated energy density \(\Omega_v\) as input A.

alpha2bool, optional

Use alternative convention for spectrum.

piecewisebool, optional

Return piecewise broken power law. Corresponds to \(\alpha \to \infty\).

dlogkbool, optional

Use spectrum defined per unit \(\log(k)\) or per unit k.

Returns

specndarray

Spectrum array.

cosmoGW.GW_analytical.smoothed_double_bPL(k, kpeak1, kpeak2, A=1.0, a=5, b=1, c=0.6666666666666666, alp1=2, alp2=2, kref=1.0, alpha2=False)

Return the value of the smoothed double broken power law (double bPL) model.

The spectrum is of the form:

\[\zeta(K) = A K^a \times \bigl(1 + (K/K1)^{(b - a) \alpha_1}\bigr)^{1/\alpha_1} \times \bigl(1 + (K/K2)^{(c + b) \alpha_2}\bigr)^{-1/\alpha_2},\]

where \(K = k/k_\ast\), K1 and K2 are the two position peaks, a is the low-k slope, b is the intermediate slope, and -c is the high-k slope. \(\alpha_1\) and \(\alpha_2\) are the smoothness parameters for each spectral transition.

The same broken power law with a slightly different form is used in Caprini:2024gyk, Caprini:2024hue and can be used setting alpha2 = True:

\[\zeta(K) = A K^a \times \bigl(1 + (K/K1)^{\alpha_1}\bigr)^{(b - a)/\alpha_1} \times \bigl(1 + (K/K2)^{\alpha_2}\bigr)^{-(c + b)/\alpha_2},\]

Parameters

karray_like

Array of wave numbers.

kpeak1, kpeak2float

Peak positions.

Afloat, optional

Amplitude of the spectrum.

afloat, optional

Slope at low wave numbers (default 5).

bfloat, optional

Slope at intermediate wave numbers (default 2/3).

cfloat, optional

Slope at high wave numbers (default 2/3).

alp1, alp2float, optional

Smoothness of the transitions.

kreffloat, optional

Reference wave number used to normalize the spectrum.

alpha2bool, optional

Use different normalization.

Returns

specndarray

Spectrum array.

References

RoperPol:2023dzg, equation 50 RoperPol:2023bqa, equation 7

cosmoGW.GW_analytical.zetam(a=5, b=0.6666666666666666, alp=2, m=-1, norm=True, alpha2=False, piecewise=False)

Compute the amplitude at the peak of the function

\[\zeta_m (K) = K^m \zeta(K),\]

where \(\zeta(K)\) is the original smoothed double broken power law (dbP).

Parameters

a, bfloat

Slopes of the spectrum (default 5 and 2/3).

alpfloat

Smoothness of the transition \(\alpha\) (default 2).

mint

Power of \(K^m\) (default -1).

normbool

Normalize the spectrum (default True).

alpha2bool

Use alternative convention (default False).

piecewisebool

Use piecewise broken power law (default False).

Returns

zetamfloat

Amplitude of new function \(K^m \zeta(K)\).