geff.tools.gw
This module is intended to compute a gravitational-wave spectrum from a tensor power spectrum.
The gravitational-wave spectrum, $h^2\Omega_{\rm GW}$ is computed using omega_gw, which evaluates the formula
$$\Omega_{\rm GW}(f) \equiv \frac{1}{3 H_0^2 M_{\rm P}^2}\frac{{\rm d} \rho_{\rm GW} (f)}{{\rm d} \ln{f}} = \frac{\pi^2}{3 H_0^2} f^2 |\mathcal{T}_{\rm GW}(f)|^2 \mathcal{P}_T(k_f), k_f) \, , \quad k_f = 2 \pi a_0 f \, ,$$
where $\mathcal{P}_T(k)$ is the tensor power spectrum with momentum $k$ at the end of inflation, and $H_0$ is the Hubble rate today. The transfer function $|\mathcal{T}_{\rm GW}(f)|^2$ is given by
$$|\mathcal{T}_{\rm GW}(f)|^2 \simeq \frac{H_0^2\Omega_r}{8 \pi^2 f^2} \frac{g_{*}(T_f)}{g_{*}(T_0)} \left(\frac{g_{*,S}(T_0)}{g_{*,S}(T_f)}\right)^{4/3} \left(1 + \frac{9}{16}\left(\frac{f_{\rm eq}}{\sqrt{2} f} \right)^2\right) |\mathcal{T}_{\rm reh}(f)|^2 \,.$$
It accounts for the evolution of $\Omega_{\rm GW}(f)$ from the end of inflation until today.
Here, $T_f$ is the temperature corresponding to the frequency $f$.
For $g_{*}$, $g_{*,S}$, $f_{\rm eq}$, $\Omega_r$, and $T_0$, we use the corresponding functions in geff.utility.cosmo.
The term $|\mathcal{T}_{\rm reh}(f)|^2$ accounts for the transition through reheating. For instantaneous reheating, $|\mathcal{T}_{\rm reh}(f)|^2 = 1$.
Otherwise, we assume that reheating proceeds via coherent oscillations of the inflaton field, such that
$$ |\mathcal{T}_{\rm reh}(f)|^2 = \frac{ \theta(f_{\rm end} - f)}{1 - 0.22 \left(\frac{f}{f_{\rm reh}}\right)^{1.5} + 0.65 \left(\frac{f}{f_{\rm reh}}\right)^{2}} \, , $$
as given in arXiv:1407.4785 and arXiv:2011.03323. Here, $f_{\rm reh}$ and $f_{\rm end}$ are, respectively, the frequencies at the end of reheating, and inflation.
The frequency $f$ can be computed from a comoving momentum $k$ using k_to_f, which evaluates
$$ f = \frac{k_f}{2 \pi a_0} =\frac{k_f}{2 \pi a_{\rm end}} e^{-N_{reh}} \left( \frac{g_{*,S}(T_0)}{g_{*,S}(T_{{\rm reh}})}\right)^{1/3} \frac{T_0}{T_{{\rm reh}}} \, . $$
where
$$ N_{\rm reh} = \frac{1}{3(1 + w_{\rm reh})} \ln \left(\frac{90 M_{\rm P} H_{\rm end}^2}{\pi^2 g_*(T_{\rm reh}) T_{\rm reh}^4} \right) \, . $$ and $T_{\rm reh}$ is the temperature of the SM plasma at the end of reheating.
1import numpy as np 2 3from geff._docs.docs_gw import DOCS 4from geff.utility.cosmo import g_rho, g_rho_freq, g_rho_0, g_s, g_s_freq, g_s_0, T_0, M_pl, gev_to_hz, omega_r, h, feq 5 6from numpy.typing import ArrayLike 7from typing import Tuple 8 9__doc__ = DOCS["module"] 10 11def omega_gw(k:np.ndarray, PT:np.ndarray, Nend:float, Hend:float, Trh:None|float=None) -> Tuple[np.ndarray, np.ndarray]: 12 r""" 13 Compute $h^2 \Omega_{\rm GW}(f)$ from a tensor power spectrum. 14 15 Parameters 16 ---------- 17 k : NDarray 18 momenta in Planck units 19 PT : NDarray 20 the tensor power spectrum at the end of inflation as a function of momentum 21 Nend : float 22 the number of e-folds at the end of inflation 23 Hend : float 24 the Hubble rate at the end of inflation (in Planck units) 25 Trh : None or float 26 the reheating temperature in GeV. If None, instantaneous reheating is assumed. 27 28 Returns 29 ------- 30 f : NDarray 31 frequencies today (in Hz) 32 h2OmegaGw : NDarray 33 the gravitational-wave spectrum as a function of frequency today 34 """ 35 36 f = k_to_f(k, Nend, Hend, Trh) 37 if Trh is None: 38 TransferRH=1 39 else: 40 frh = 1/(2*np.pi) * (g_s_0/g_s(Trh))**(1/3) * (np.pi**2*g_rho(Trh)/90)**(1/2) * (Trh/M_pl) * T_0*gev_to_hz 41 fend = 1/(2*np.pi) * (g_s_0/g_s(Trh))**(1/3) * (np.pi**2*g_rho(Trh)/90)**(1/3) * (Trh/M_pl)**(1/3) * (Hend)**(1/3) * T_0*gev_to_hz 42 43 TransferRH = 1/(1. - 0.22*(f/frh)**1.5 + 0.65*(f/frh)**2) 44 TransferRH = np.where(np.log(f/fend) < 0, TransferRH, np.zeros(f.shape)) 45 46 TransferMD = 1 + 9/32*(feq/f)**2 47 48 h2OmegaGW = h**2*omega_r/24 * PT * (g_rho_freq(f)/g_rho_0) * (g_s_0/g_s_freq(f))**(4/3) * TransferMD * TransferRH 49 50 return f, h2OmegaGW 51 52def k_to_f(k:np.ndarray, Nend:float, Hend:float, Trh:None|float=None) -> ArrayLike: 53 r""" 54 Compute frequency today from momentum at the end of inflation 55 56 Parameters 57 ---------- 58 k : NDarray 59 momenta in Planck units 60 Nend : float 61 the number of e-folds at the end of inflation 62 Hend : float 63 the Hubble rate at the end of inflation 64 Trh : None or float 65 the reheating temperature in GeV. If None, instantaneous reheating is assumed. 66 67 Return 68 ------ 69 f : NDarray 70 frequencies today (in Hz) 71 """ 72 73 if Trh is None: 74 Trh = np.sqrt(3*Hend/np.pi)*(10/106.75)**(1/4)*M_pl 75 Trh = Trh*(106.75/g_rho(Trh))**(1/4) 76 Nrh = 0 77 else: 78 wrh = 0 79 Nrh = np.log( 90*(Hend*M_pl**2)**2 / (np.pi**2*g_rho(Trh)*Trh**4 ) ) / ( 3 * (1 + wrh) ) 80 81 82 f = k*M_pl*gev_to_hz/(2*np.pi*np.exp(Nend)) * T_0/Trh * (g_s_0/g_s(Trh))**(1/3) * np.exp(-Nrh) 83 84 return f
12def omega_gw(k:np.ndarray, PT:np.ndarray, Nend:float, Hend:float, Trh:None|float=None) -> Tuple[np.ndarray, np.ndarray]: 13 r""" 14 Compute $h^2 \Omega_{\rm GW}(f)$ from a tensor power spectrum. 15 16 Parameters 17 ---------- 18 k : NDarray 19 momenta in Planck units 20 PT : NDarray 21 the tensor power spectrum at the end of inflation as a function of momentum 22 Nend : float 23 the number of e-folds at the end of inflation 24 Hend : float 25 the Hubble rate at the end of inflation (in Planck units) 26 Trh : None or float 27 the reheating temperature in GeV. If None, instantaneous reheating is assumed. 28 29 Returns 30 ------- 31 f : NDarray 32 frequencies today (in Hz) 33 h2OmegaGw : NDarray 34 the gravitational-wave spectrum as a function of frequency today 35 """ 36 37 f = k_to_f(k, Nend, Hend, Trh) 38 if Trh is None: 39 TransferRH=1 40 else: 41 frh = 1/(2*np.pi) * (g_s_0/g_s(Trh))**(1/3) * (np.pi**2*g_rho(Trh)/90)**(1/2) * (Trh/M_pl) * T_0*gev_to_hz 42 fend = 1/(2*np.pi) * (g_s_0/g_s(Trh))**(1/3) * (np.pi**2*g_rho(Trh)/90)**(1/3) * (Trh/M_pl)**(1/3) * (Hend)**(1/3) * T_0*gev_to_hz 43 44 TransferRH = 1/(1. - 0.22*(f/frh)**1.5 + 0.65*(f/frh)**2) 45 TransferRH = np.where(np.log(f/fend) < 0, TransferRH, np.zeros(f.shape)) 46 47 TransferMD = 1 + 9/32*(feq/f)**2 48 49 h2OmegaGW = h**2*omega_r/24 * PT * (g_rho_freq(f)/g_rho_0) * (g_s_0/g_s_freq(f))**(4/3) * TransferMD * TransferRH 50 51 return f, h2OmegaGW
Compute $h^2 \Omega_{\rm GW}(f)$ from a tensor power spectrum.
Parameters
- k (NDarray): momenta in Planck units
- PT (NDarray): the tensor power spectrum at the end of inflation as a function of momentum
- Nend (float): the number of e-folds at the end of inflation
- Hend (float): the Hubble rate at the end of inflation (in Planck units)
- Trh (None or float): the reheating temperature in GeV. If None, instantaneous reheating is assumed.
Returns
- f (NDarray): frequencies today (in Hz)
- h2OmegaGw (NDarray): the gravitational-wave spectrum as a function of frequency today
53def k_to_f(k:np.ndarray, Nend:float, Hend:float, Trh:None|float=None) -> ArrayLike: 54 r""" 55 Compute frequency today from momentum at the end of inflation 56 57 Parameters 58 ---------- 59 k : NDarray 60 momenta in Planck units 61 Nend : float 62 the number of e-folds at the end of inflation 63 Hend : float 64 the Hubble rate at the end of inflation 65 Trh : None or float 66 the reheating temperature in GeV. If None, instantaneous reheating is assumed. 67 68 Return 69 ------ 70 f : NDarray 71 frequencies today (in Hz) 72 """ 73 74 if Trh is None: 75 Trh = np.sqrt(3*Hend/np.pi)*(10/106.75)**(1/4)*M_pl 76 Trh = Trh*(106.75/g_rho(Trh))**(1/4) 77 Nrh = 0 78 else: 79 wrh = 0 80 Nrh = np.log( 90*(Hend*M_pl**2)**2 / (np.pi**2*g_rho(Trh)*Trh**4 ) ) / ( 3 * (1 + wrh) ) 81 82 83 f = k*M_pl*gev_to_hz/(2*np.pi*np.exp(Nend)) * T_0/Trh * (g_s_0/g_s(Trh))**(1/3) * np.exp(-Nrh) 84 85 return f
Compute frequency today from momentum at the end of inflation
Parameters
- k (NDarray): momenta in Planck units
- Nend (float): the number of e-folds at the end of inflation
- Hend (float): the Hubble rate at the end of inflation
- Trh (None or float): the reheating temperature in GeV. If None, instantaneous reheating is assumed.
Return
f : NDarray frequencies today (in Hz)