geff.tools.pt

$\newcommand{\bm}[1]{\boldsymbol{#1}}$

This module is designed to compute the tensor power spectrum from a GEF result.

The polarized tensor power spectrum at time $t$ with momentum $k$ is given by $$\mathcal{P}_{T}(t, k) = \mathcal{P}_{T}^{\mathrm{vac}}(t, k) + \frac{1}{2}\sum_{\lambda=\pm 1} \mathcal{P}_{T,\lambda}^{\mathrm{ind}}(t, k)\, ,$$ where $\lambda$ labels the helicity of tensor modes.

The vacuum contribution is given by

$$\mathcal{P}_{T}^{\mathrm{vac}}(t, k) = \frac{4 k^3}{\pi^2 M_{\rm P}^2} |u_0(t, k)|^2\, ,$$

while the induced contribution is

$$\mathcal{P}_{T,\lambda}^{\mathrm{ind}}(t, k) = \frac{k^3}{2 \pi^2 M_{\rm P}^4} \int \frac{{\rm d}^3 \bm{p}}{(2 \pi)^3} \sum_{\alpha,\beta = \pm1} \left(1 + \lambda \alpha \frac{\bm{k} \cdot \bm{p}}{k p} \right)^2 \left(1 + \lambda \beta \frac{k^2 - \bm{k} \cdot \bm{p}}{kq} \right)^2 $$ $$ \qquad \qquad \qquad \times \left|\int_{-\infty}^\infty {\rm d} s \frac{G_k(t, s)}{a^3(s)} \left[A'_\alpha(s, p)A'_\beta(s, q) + \alpha \beta\, p q\, A_\alpha(s, p) A_\alpha(s, q) \right] \right|^2 \, , $$

with momentum $q = |\bm{p} + \bm{q}|$ and scale-factor $a$.

The vacuum modes $u_0(t,k)$ obey the mode equation $$\mathcal{D}_k {u_0} = \ddot{u}_0 + 3 H \dot{u}_0 + \frac{k^2}{a^2} {u_0} = 0$$ with the retarded Green function $G_k(t',t)$ defined for the operator $\mathcal{D}_k$.

The gauge-field mode functions $A_\lambda(t,k)$ are defined as in the geff.mode_by_mode module.

For details on the numerical computation, see the Appendix B of 2508.00798.

  1import numpy as np
  2from scipy.interpolate import CubicSpline, PchipInterpolator
  3from scipy.integrate import solve_ivp, trapezoid
  4
  5from geff._docs.docs_pt import DOCS
  6from geff.mbm import GaugeSpec
  7from geff.bgtypes import BGSystem, Constant
  8
  9from typing import Tuple
 10from types import NoneType
 11
 12__doc__ = DOCS["module"]
 13
 14class PowSpecT:
 15    r"""
 16    A class used to compute the tensor power spectrum including vacuum and gauge-field induced contributions.
 17
 18    These main method of this module is `compute_pt`, which computes both the vacuum and gauge-field induced contribution
 19      to the tensor power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}$ and $\mathcal{P}_{T,\lambda}^{\mathrm{ind}}$.
 20
 21    Results are internally computed using numerical units, but are returned in physical units.
 22    """
 23    def __init__(self, insys : BGSystem):
 24        """
 25        Initialise the class from a GEF solution.
 26
 27        Parameters
 28        ----------
 29        insys : BGSystem
 30            the GEF solution.
 31        """
 32        #Set GEF results to Hubble units.
 33        sys = BGSystem.from_system(insys, True)
 34        sys.units = False
 35        
 36        #import the background evolution
 37        N = sys.N.value
 38        a = np.exp(N)
 39        H = sys.H.value
 40
 41        self._omega = sys.omega
 42            
 43        #Set the range of modes
 44        self.maxk = max(a*H)
 45        self.mink = 1e4
 46
 47        #Define Useful quantities
 48
 49        self._t = sys.t.value
 50        self._N = N
 51        self._H = H
 52        self._xi= sys.xi.value
 53
 54        #define interpolated quantities as needed
 55        self._af = CubicSpline(self._t, a)
 56        if isinstance(self._H, Constant):
 57            self._Hf = lambda t: self._H
 58        else:
 59            self._Hf = CubicSpline(self._t, H)
 60        self._HN = CubicSpline(self._N, self._H)
 61        self._khN = CubicSpline(N, sys.kh)
 62
 63        #Obtain eta as a functio of time for phases
 64        def deta(t, y): return 1/self._af(t)
 65        
 66        soleta = solve_ivp(deta, [min(self._t), max(self._t)], np.array([0]), t_eval=self._t)
 67
 68        self._etaf = CubicSpline(self._t, soleta.y[0,:])
 69        return
 70    
 71    def compute_pt(self, nmodes : int, spec : GaugeSpec, FastGW : bool=True,
 72                    atols : list=[1e-3,1e-20], rtols : list=[1e-4,1e-4], momgrid : int=100
 73                    ) -> Tuple[np.ndarray,dict]:
 74        r"""
 75        Compute the full tensor power spectrum.
 76
 77        The method lodes data on gauge modes $A_\lambda(t,k)$ from a file indicated by `mbm_file`.
 78        The power spectrum is evaluated fpr `nmodes` log-spaced momenta $k \in [10^{4}a(0)H(0), 10 a(t_{\rm max}) H(t_{\rm max})]$
 79
 80        Parameters
 81        ----------
 82        nmodes : NDArray
 83            the number of momenta
 84        ModePath : str
 85            path to a file containing gauge-field mode functions
 86        FastGW : bool
 87            If `True`, only the expected dominant contributions to $\mathcal{P}_{T,\lambda}^{\mathrm{ind}}$ are computed
 88        atols : list
 89            absolute tolerance for `compute_homogeneous` (index 0) and `compute_green` (index 1)
 90        rtols : list
 91            relative tolerance for `compute_homogeneous` (index 0) and `compute_green` (index 1)
 92        momgrid : int
 93            passed to `compute_ptind`
 94        
 95        Returns
 96        -------
 97        ks : NDArray
 98            the momenta $k$
 99        PT : dict
100            a dictionary containing all contributions to the tensor power spectrum and its total.
101        """
102
103        k = np.logspace(np.log10(self.mink), np.log10(10*self.maxk), nmodes)
104        ks, tstarts = self._find_tinit_bd(k, mode="k")
105        
106        Ngrid = spec["N"]
107        tgrid = spec["t"]
108        pgrid = spec["k"]
109
110        GaugeModes = {"+":(spec["Ap"], spec["dAp"]), "-":(spec["Am"], spec["dAm"])}
111        
112        indend = np.argmin(abs(Ngrid-max(self._N)))
113
114        PT = {"tot":[], "vac":[], "ind+,++":[], "ind+,+-":[], "ind+,--":[], "ind-,++":[], "ind-,+-":[], "ind-,--":[]}
115
116        GWpols = [("+", 1), ("-",-1)]
117
118        gaugepols=[(("+",1),("+",1)),
119                    (("+",1),("-",-1)),
120                        (("-",-1),("-",-1))]
121
122        sign = np.sign(self._xi[0])
123
124        for i, k in enumerate(ks):
125            tstart = tstarts[i]
126
127            if k > 10**(5/2)*self.maxk:
128                for key in PT.keys():
129                    PT[key].append(0)
130            else:
131                uk, _ = self.compute_homogeneous(k, tstart, tgrid, atol=atols[0], rtol=rtols[0])
132                Green = self.compute_green(k, uk, tstart, tgrid, indend, atol=atols[1], rtol=rtols[1])
133                
134                PT["vac"].append(self.compute_ptvac(k, uk[indend]))
135
136                for lgrav in GWpols:
137                    for mu in gaugepols:
138                        GWpol = lgrav[1]
139
140                        Ax = GaugeModes[mu[0][0]][0]
141                        dAx = GaugeModes[mu[0][0]][1]
142                        lx = mu[0][1]
143                        Ay = GaugeModes[mu[1][0]][0]
144                        dAy = GaugeModes[mu[1][0]][1]
145                        ly = mu[1][1]
146
147                        if FastGW and (lx != sign or ly != sign):
148                            PT[f"ind{lgrav[0]},{mu[0][0]}{mu[1][0]}"].append(0.)
149                        else:
150                            PT[f"ind{lgrav[0]},{mu[0][0]}{mu[1][0]}"].append(
151                                self.compute_ptind(k, GWpol, Ngrid, indend, Green, pgrid, lx, Ax, dAx, ly, Ay, dAy, momgrid) )
152
153        PT["tot"] = np.zeros(ks.shape)
154        for key in PT.keys():
155            PT[key] = np.array(PT[key])
156            if ("+" in key) or ("-" in key):
157                PT["tot"] += 0.5*PT[key]
158            elif key=="vac":
159                PT["tot"] += PT[key]
160
161        #this returns k with units restored. Needs to be matched by omega_GW
162        return ks*self._omega, PT
163    
164    def compute_homogeneous(self, k : float, tvac : float, teval : np.ndarray|NoneType, atol : float=1e-3, rtol : float=1e-4) -> Tuple[np.ndarray, np.ndarray]:
165        r"""
166        Compute the evolution of a vacuum mode starting from Bunch–Davies vacuum.
167
168        For a mode $u_0(t,k)$ with momentum $k$, initialise in Bunch–Davies when $k = 10^{5/2} a(t_{\rm vac})H(t_{\rm vac})$.
169        Its evolution is obtained by solving `tensor_mode_eq` using `scipy.integrate.solve_ivp`.
170
171        Parameters
172        ----------
173        k : float
174           the momentum $k$.
175        tvac : float
176            the initial time, $t_{\rm init}$
177        teval : NDArray or None
178            time coordinates $t$ of $u_0(t,k)$ (`None`: same as class input `sys`) 
179        atol : float
180            the absolute tolerance for `solve_ivp`
181        rtol : float
182            the relative tolerance for `solve_ivp`
183
184        Returns
185        -------
186        u : NDArray
187            the vacuum tensor modes, $\sqrt{2k}u(t, k)$
188        duk : NDArray
189            the derivative of the vacuum tensor modes, $a\sqrt{2/k}\dot{u}(t, k)$
190        """
191
192        if len(teval)==0:
193            teval = self._t
194        tend = max(teval)
195
196        #conformal time needed for relative phases
197        eta = self._etaf(teval)
198
199        istart = 0
200        while teval[istart]<tvac:
201            istart+=1
202        
203        #define the ODE for the GW modes
204        def ode(t, y): return self.tensor_mode_eq( y, k, self._Hf(t), self._af(t) )
205        
206        #Initialise the modes in Bunch Davies
207        Zini = np.array([1, -10**(-5/2), 0, -1])
208
209        sol = solve_ivp(ode, [tvac, tend], Zini, t_eval=teval[istart:], method="RK45", atol=atol, rtol=rtol)
210        if not(sol.success):
211            print("Something went wrong")
212
213        #the mode was in vacuum before tvac
214        vac = list( np.exp(-1j*eta[:istart]*k) )
215        dvac = list( -1j*np.exp(-1j*eta[:istart]*k) )
216
217        #Create an array tracking a modes evolution from Bunch Davies to late times. Ensure equal length arrays for every mode k
218        uk = np.array( vac + list( (sol.y[0,:] + 1j*sol.y[2,:])*np.exp(-1j*k*eta[istart]) ) )/self._af(teval)
219        duk = np.array( dvac + list( (sol.y[1,:] + 1j*sol.y[3,:])*np.exp(-1j*k*eta[istart]) ) )/self._af(teval)
220
221        return uk, duk
222    
223    @staticmethod
224    def tensor_mode_eq(y : np.ndarray, k : float, H : float, a : float) -> np.ndarray:
225        r"""
226        Mode equation for vacuum tensor modes.
227
228        Parameters
229        ----------
230        y : NDArray
231            the vacuum tensor mode and its derivatives.
232            $a \sqrt{2k} u_0(t,k)$ and $ a^2\sqrt{2/k} \dot{u}_0(t, k)$ 
233        k : float
234            comoving momentum $k$
235        H : float
236            Hubble rate, $H(t)$
237        a : float
238            scale factor, $a(t)$
239
240        Returns
241        -------
242        dydt : NDArray
243            an array of time derivatives of y
244
245        """
246        dydt = np.zeros(y.shape)
247
248        #real
249        dydt[0] = (H*y[0] + k/a*y[1])
250        dydt[1] = ( -H*y[1] - k/a*y[0] )
251
252        #imaginary
253        dydt[2] = (H*y[2] + k/a*y[3])
254        dydt[3] = ( -H*y[3] - k/a*y[2] )
255
256        return dydt
257    
258    def compute_green(self, k : float, uk : np.ndarray, tvac : float, teval : np.ndarray|NoneType, ind : int, atol : float=1e-30, rtol : float=1e-4) -> np.ndarray:
259        r"""
260        Compute the evolution of the Green function $G_k(t',t)$ for fixed $t'$.
261
262        The evolution is obtained by solving `green_ode` from $t=t'$ backwards until $t_{\rm vac}$, defined through $k = 10^{5/2} a(t_{\rm end})H(t_{\rm vac})$.
263        From $t_{\rm end}$ onwards, the Green function is instead computed from the mode $u_0(t, k)$. The method uses `scipy.integrate.solve_ivp`.
264
265        Parameters
266        ----------
267        k : float
268           the momentum $k$.
269        uk : NDArray
270            the mode function $u_0(t,k)$
271        tvac : float
272            the final time, $t_{\rm vac}$
273        teval : NDArray or None
274            time coordinates $t$ of $G_k(t',t)$ (`None`: same as class input `sys`) 
275        ind : int
276            index of teval corresponding to $t'$
277        atol : float
278            the absolute tolerance for `solve_ivp`
279        rtol : float
280            the relative tolerance for `solve_ivp`
281
282        Return
283        ------
284        GreenN : NDArray
285            the Green function $-k G_k(t', t)$ 
286        """
287        if len(teval)==0:
288            teval = self._t
289            
290        istart = 0
291        while teval[istart]<tvac:
292            istart+=1
293
294        # G(k, t, t) = 0 by definition
295        Aini = np.array([0, 1])
296    
297
298        # Solve the EoM for G backwards in time starting from G(k, t, t)
299        def Aode(t, y): return -self.green_ode(y, k, self._Hf(-t), self._af(-t))
300        solA = solve_ivp(Aode, [-teval[ind], -tvac], Aini, t_eval=-teval[istart:ind+1][::-1],
301                         method="RK45", atol=atol, rtol=rtol)
302
303        #For numerical stability, only solve the EoM for the Green function until t' = tvac. Afterwards, compute it directly from the vacuum modes.
304        GreenN = np.zeros(teval.shape)
305        GreenN[istart:ind+1] = solA.y[0,:][::-1]
306        GreenN[:istart] = ( (uk[ind].conjugate()*uk).imag*self._af(teval)**2 )[:istart]
307
308        return GreenN
309    
310    @staticmethod
311    def green_ode(y : np.ndarray, k : float, H : float, a : float) -> np.ndarray:
312        r"""
313        Evolve the re-scaled Green function for tensor modes in time.
314
315        The re-scaled Green function is
316        $$B_{t'}^k(t) = 2 k a(t)^2 \operatorname{Im} \left[ u_0^*(t', k) \, u_0(t,k) \right] = k G_k(t', t) \,.$$
317        
318        whose ODE is coupled to 
319
320        $$C_{t'}^k(t) = 2 a(t)^3 \operatorname{Im} \left[ u_0^*(t', k) \, {\dot{u}_0}(t,k) \right] \,.$$
321
322        Parameters
323        ----------
324        y : NDArray
325            the functions B and C
326        k : float
327            comoving momentum $k$
328        H : float
329            Hubble rate, $H(t)$
330        a : float
331            scale factor, $a(t)$
332
333        Returns
334        -------
335        dydt : NDArray
336            an array of time derivatives of A
337
338        """
339        dydt = np.zeros(y.shape)
340        dydt[0] =(2*H*y[0] + k/a*y[1])
341        dydt[1] = -k/a*y[0]
342        return dydt
343
344    
345    def _find_tinit_bd(self, init : np.ndarray, mode : str="k") -> Tuple[np.ndarray, np.ndarray]:
346        """
347        Determines the pair of $k$ and $t$ satisfying $k = 10^(5/2)a(t)H(t)$.
348
349        Depending on `mode`, `init` may be a time coordinate (`mode='t'`), $e$-folds (`mode='N'`) or momentum (`mode='k'`).
350
351        Parameters
352        ----------
353        init : array
354            the input array (t, N, or k)
355        mode : str
356            indicate the type of `init`
357
358        Returns
359        -------
360        ks : NDarray
361            an array of momenta
362        tstart : NDarray
363            an array of times
364
365        Raises
366        ------
367        KeyError
368            if `mode` is not 't, 'k' or 'N'
369        """
370
371        pwr = 5/2
372
373        t = self._t
374        def logkH(t): return np.log(self._af(t)*self._Hf(t))
375        if mode=="t":
376            tstart = init
377            logks = logkH(tstart)
378            ks = 10**(pwr)*np.exp(logks)
379
380        elif mode=="k":
381            ks = init
382
383            tstart = []
384            for k in ks:
385                ttmp  = self._t[np.searchsorted(10**(pwr)*np.exp(logkH(self._t)), k, "right")]
386                tstart.append(ttmp)
387            tstart = np.array(tstart)
388
389        elif mode=="N":
390            tstart = CubicSpline(self._N, t)(init)
391            ks = 10**pwr*np.exp(logkH(tstart))
392
393        else:
394            raise KeyError("'mode' must be 't', 'k' or 'N'")
395
396        return ks, tstart
397        
398    def compute_ptvac(self, k : np.ndarray, uk : np.ndarray) -> np.ndarray:
399        r"""
400        Compute the vacuum power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}(t,k)$.
401
402        Parameters
403        ----------
404        k : NDArray
405            momentum, $k$
406        uk : NDArray
407            vacuum tensor mode $u_0(t,k)$
408
409        Return
410        ------
411        PTvac : NDArray
412            the vacuum tensor power spectrum
413        """
414        PTvac = 2*(k*self._omega)**2/( np.pi**2 ) * abs(uk)**2
415        return PTvac
416
417    def compute_ptind(self, k : float, lgrav : float, Ngrid : np.ndarray, ind: int, GreenN : np.ndarray, pgrid : np.ndarray,
418                                l1 : float, A1 : np.ndarray, dA1 : np.ndarray, l2 : float, A2 : np.ndarray, dA2 : np.ndarray,
419                                momgrid : int=100) -> float:
420        r"""
421        Compute the vacuum power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}(t,k)$.
422
423        The integral is computed by integrating over the inner integral over the gauge-mode functions $A_\mu(s,p)$ using `scipy.integrate.trapezoid`.
424        The external momentum integrals are also computed using `scipy.integrate.trapezoid` on a grid of momenta momgrid x momgrid. The mode functions
425        are interpolated over momentum to match this grid.
426
427        Parameters
428        ----------
429        k : float
430            momentum $k$
431        lgrav : NDArray
432            helicity $\lambda$
433        Ngrid : NDArray
434            $e$-folds $\log a$
435        ind : integer
436            the power spectrum is computed at $t$ corresponding to `Ngrid[ind]`
437        GreenN : NDArray
438            the Green function $kG_k(t, s)$ with $s$ corresponding to `Ngrid`
439        pgrid : NDArray
440            momenta $p$ for the mode functions $A_\mu(s,p)$
441        l1, l2 : float
442            the gauge-field helicities $\mu_1$, $\mu_2$
443        A1, A2 : float
444            the gauge-field mode functions $\sqrt{2k} A_{\mu_1, \mu_2}(s, p)$
445        dA1, dA2 : float
446            the gauge-field mode function derivatives $a(s)\sqrt{2/k} \dot{A}_{\mu_1, \mu_2}(s, p)$
447        momgrid : int
448            the grid size for the momentum integrals.
449
450        Returns
451        -------
452        PTind : float
453            the gauge-field induced tensor power spectrum
454        """
455
456        cutUV = self._khN(Ngrid[ind])/k
457        cutIR = min(pgrid)/k
458        HN = self._HN(Ngrid)
459
460        logAs = np.linspace(np.log(max(0.5, cutIR)), np.log(cutUV), momgrid)
461
462        Ax_rad = PchipInterpolator(np.log(pgrid), abs(A1))
463        Ax_phase = PchipInterpolator(np.log(pgrid), np.arccos(A1.real/abs(A1)))
464        
465        dAx_rad = PchipInterpolator(np.log(pgrid), abs(dA1))
466        dAx_phase = PchipInterpolator(np.log(pgrid), np.arccos(dA1.real/abs(dA1)))
467
468        Ay_rad = PchipInterpolator(np.log(pgrid), abs(A2))
469        Ay_phase = PchipInterpolator(np.log(pgrid), np.arccos(A2.real/abs(A2)))
470        
471        dAy_rad = PchipInterpolator(np.log(pgrid), abs(dA2))
472        dAy_phase = PchipInterpolator(np.log(pgrid), np.arccos(dA2.real/abs(dA2)))
473
474        def Afuncx(x): return Ax_rad(x)*np.exp(1j*Ax_phase(x))
475        def dAfuncx(x): return dAx_rad(x)*np.exp(1j*dAx_phase(x))
476
477        def Afuncy(x): return Ay_rad(x)*np.exp(1j*Ay_phase(x))
478        def dAfuncy(x): return dAy_rad(x)*np.exp(1j*dAy_phase(x))
479
480        IntOuter = []
481        for logA in logAs:
482            A = np.exp(logA)
483            Blow = (cutIR - A)
484            Bhigh = (A - cutIR)
485            if Bhigh>0.5:
486                Blow = (A - cutUV)
487                Bhigh = (cutUV - A)
488                if Bhigh>0.5:
489                    Blow = -0.5
490                    Bhigh = 0.5
491            Bs = np.linspace(Blow, Bhigh, momgrid)[1:-1]
492
493            IntInner = np.zeros(Bs.shape)
494
495            for j, B in enumerate(Bs):
496                Ax = Afuncx(np.log(k*(A+B)))
497                dAx = dAfuncx(np.log(k*(A+B)))
498                Ay = Afuncy(np.log(k*(A-B)))
499                dAy = dAfuncy(np.log(k*(A-B)))
500
501                mom =  abs( l1*l2 + 2*lgrav*( (l1+l2)*A + (l1-l2)*B ) + 4*(A**2 - B**2) + 8*lgrav*A*B*( (l1-l2)*A - (l1+l2)*B ) - 16*l1*l2*A**2*B**2 )
502                z = max(A+B,A-B)
503                mask = np.where(z<self._khN(Ngrid)/k, 1, 0)
504
505                val = (dAx*dAy + l1*l2*Ax*Ay)*mask*k/(np.exp(3*Ngrid)*HN)
506                timeintegrand = GreenN*val*mom
507            
508                timeintre = trapezoid(timeintegrand[:ind].real, Ngrid[:ind])
509                
510                timeintim = trapezoid(timeintegrand[:ind].imag, Ngrid[:ind])
511                
512                IntInner[j] = (timeintre**2 + timeintim**2)*A
513
514  
515            IntOuter.append(trapezoid(IntInner, Bs))
516        IntOuter = np.array(IntOuter)
517
518        PTind = trapezoid(IntOuter, logAs) / (16*np.pi**4)*(k*self._omega)**4
519
520        return PTind
521
class PowSpecT: 
 15class PowSpecT:
 16    r"""
 17    A class used to compute the tensor power spectrum including vacuum and gauge-field induced contributions.
 18
 19    These main method of this module is `compute_pt`, which computes both the vacuum and gauge-field induced contribution
 20      to the tensor power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}$ and $\mathcal{P}_{T,\lambda}^{\mathrm{ind}}$.
 21
 22    Results are internally computed using numerical units, but are returned in physical units.
 23    """
 24    def __init__(self, insys : BGSystem):
 25        """
 26        Initialise the class from a GEF solution.
 27
 28        Parameters
 29        ----------
 30        insys : BGSystem
 31            the GEF solution.
 32        """
 33        #Set GEF results to Hubble units.
 34        sys = BGSystem.from_system(insys, True)
 35        sys.units = False
 36        
 37        #import the background evolution
 38        N = sys.N.value
 39        a = np.exp(N)
 40        H = sys.H.value
 41
 42        self._omega = sys.omega
 43            
 44        #Set the range of modes
 45        self.maxk = max(a*H)
 46        self.mink = 1e4
 47
 48        #Define Useful quantities
 49
 50        self._t = sys.t.value
 51        self._N = N
 52        self._H = H
 53        self._xi= sys.xi.value
 54
 55        #define interpolated quantities as needed
 56        self._af = CubicSpline(self._t, a)
 57        if isinstance(self._H, Constant):
 58            self._Hf = lambda t: self._H
 59        else:
 60            self._Hf = CubicSpline(self._t, H)
 61        self._HN = CubicSpline(self._N, self._H)
 62        self._khN = CubicSpline(N, sys.kh)
 63
 64        #Obtain eta as a functio of time for phases
 65        def deta(t, y): return 1/self._af(t)
 66        
 67        soleta = solve_ivp(deta, [min(self._t), max(self._t)], np.array([0]), t_eval=self._t)
 68
 69        self._etaf = CubicSpline(self._t, soleta.y[0,:])
 70        return
 71    
 72    def compute_pt(self, nmodes : int, spec : GaugeSpec, FastGW : bool=True,
 73                    atols : list=[1e-3,1e-20], rtols : list=[1e-4,1e-4], momgrid : int=100
 74                    ) -> Tuple[np.ndarray,dict]:
 75        r"""
 76        Compute the full tensor power spectrum.
 77
 78        The method lodes data on gauge modes $A_\lambda(t,k)$ from a file indicated by `mbm_file`.
 79        The power spectrum is evaluated fpr `nmodes` log-spaced momenta $k \in [10^{4}a(0)H(0), 10 a(t_{\rm max}) H(t_{\rm max})]$
 80
 81        Parameters
 82        ----------
 83        nmodes : NDArray
 84            the number of momenta
 85        ModePath : str
 86            path to a file containing gauge-field mode functions
 87        FastGW : bool
 88            If `True`, only the expected dominant contributions to $\mathcal{P}_{T,\lambda}^{\mathrm{ind}}$ are computed
 89        atols : list
 90            absolute tolerance for `compute_homogeneous` (index 0) and `compute_green` (index 1)
 91        rtols : list
 92            relative tolerance for `compute_homogeneous` (index 0) and `compute_green` (index 1)
 93        momgrid : int
 94            passed to `compute_ptind`
 95        
 96        Returns
 97        -------
 98        ks : NDArray
 99            the momenta $k$
100        PT : dict
101            a dictionary containing all contributions to the tensor power spectrum and its total.
102        """
103
104        k = np.logspace(np.log10(self.mink), np.log10(10*self.maxk), nmodes)
105        ks, tstarts = self._find_tinit_bd(k, mode="k")
106        
107        Ngrid = spec["N"]
108        tgrid = spec["t"]
109        pgrid = spec["k"]
110
111        GaugeModes = {"+":(spec["Ap"], spec["dAp"]), "-":(spec["Am"], spec["dAm"])}
112        
113        indend = np.argmin(abs(Ngrid-max(self._N)))
114
115        PT = {"tot":[], "vac":[], "ind+,++":[], "ind+,+-":[], "ind+,--":[], "ind-,++":[], "ind-,+-":[], "ind-,--":[]}
116
117        GWpols = [("+", 1), ("-",-1)]
118
119        gaugepols=[(("+",1),("+",1)),
120                    (("+",1),("-",-1)),
121                        (("-",-1),("-",-1))]
122
123        sign = np.sign(self._xi[0])
124
125        for i, k in enumerate(ks):
126            tstart = tstarts[i]
127
128            if k > 10**(5/2)*self.maxk:
129                for key in PT.keys():
130                    PT[key].append(0)
131            else:
132                uk, _ = self.compute_homogeneous(k, tstart, tgrid, atol=atols[0], rtol=rtols[0])
133                Green = self.compute_green(k, uk, tstart, tgrid, indend, atol=atols[1], rtol=rtols[1])
134                
135                PT["vac"].append(self.compute_ptvac(k, uk[indend]))
136
137                for lgrav in GWpols:
138                    for mu in gaugepols:
139                        GWpol = lgrav[1]
140
141                        Ax = GaugeModes[mu[0][0]][0]
142                        dAx = GaugeModes[mu[0][0]][1]
143                        lx = mu[0][1]
144                        Ay = GaugeModes[mu[1][0]][0]
145                        dAy = GaugeModes[mu[1][0]][1]
146                        ly = mu[1][1]
147
148                        if FastGW and (lx != sign or ly != sign):
149                            PT[f"ind{lgrav[0]},{mu[0][0]}{mu[1][0]}"].append(0.)
150                        else:
151                            PT[f"ind{lgrav[0]},{mu[0][0]}{mu[1][0]}"].append(
152                                self.compute_ptind(k, GWpol, Ngrid, indend, Green, pgrid, lx, Ax, dAx, ly, Ay, dAy, momgrid) )
153
154        PT["tot"] = np.zeros(ks.shape)
155        for key in PT.keys():
156            PT[key] = np.array(PT[key])
157            if ("+" in key) or ("-" in key):
158                PT["tot"] += 0.5*PT[key]
159            elif key=="vac":
160                PT["tot"] += PT[key]
161
162        #this returns k with units restored. Needs to be matched by omega_GW
163        return ks*self._omega, PT
164    
165    def compute_homogeneous(self, k : float, tvac : float, teval : np.ndarray|NoneType, atol : float=1e-3, rtol : float=1e-4) -> Tuple[np.ndarray, np.ndarray]:
166        r"""
167        Compute the evolution of a vacuum mode starting from Bunch&ndash;Davies vacuum.
168
169        For a mode $u_0(t,k)$ with momentum $k$, initialise in Bunch&ndash;Davies when $k = 10^{5/2} a(t_{\rm vac})H(t_{\rm vac})$.
170        Its evolution is obtained by solving `tensor_mode_eq` using `scipy.integrate.solve_ivp`.
171
172        Parameters
173        ----------
174        k : float
175           the momentum $k$.
176        tvac : float
177            the initial time, $t_{\rm init}$
178        teval : NDArray or None
179            time coordinates $t$ of $u_0(t,k)$ (`None`: same as class input `sys`) 
180        atol : float
181            the absolute tolerance for `solve_ivp`
182        rtol : float
183            the relative tolerance for `solve_ivp`
184
185        Returns
186        -------
187        u : NDArray
188            the vacuum tensor modes, $\sqrt{2k}u(t, k)$
189        duk : NDArray
190            the derivative of the vacuum tensor modes, $a\sqrt{2/k}\dot{u}(t, k)$
191        """
192
193        if len(teval)==0:
194            teval = self._t
195        tend = max(teval)
196
197        #conformal time needed for relative phases
198        eta = self._etaf(teval)
199
200        istart = 0
201        while teval[istart]<tvac:
202            istart+=1
203        
204        #define the ODE for the GW modes
205        def ode(t, y): return self.tensor_mode_eq( y, k, self._Hf(t), self._af(t) )
206        
207        #Initialise the modes in Bunch Davies
208        Zini = np.array([1, -10**(-5/2), 0, -1])
209
210        sol = solve_ivp(ode, [tvac, tend], Zini, t_eval=teval[istart:], method="RK45", atol=atol, rtol=rtol)
211        if not(sol.success):
212            print("Something went wrong")
213
214        #the mode was in vacuum before tvac
215        vac = list( np.exp(-1j*eta[:istart]*k) )
216        dvac = list( -1j*np.exp(-1j*eta[:istart]*k) )
217
218        #Create an array tracking a modes evolution from Bunch Davies to late times. Ensure equal length arrays for every mode k
219        uk = np.array( vac + list( (sol.y[0,:] + 1j*sol.y[2,:])*np.exp(-1j*k*eta[istart]) ) )/self._af(teval)
220        duk = np.array( dvac + list( (sol.y[1,:] + 1j*sol.y[3,:])*np.exp(-1j*k*eta[istart]) ) )/self._af(teval)
221
222        return uk, duk
223    
224    @staticmethod
225    def tensor_mode_eq(y : np.ndarray, k : float, H : float, a : float) -> np.ndarray:
226        r"""
227        Mode equation for vacuum tensor modes.
228
229        Parameters
230        ----------
231        y : NDArray
232            the vacuum tensor mode and its derivatives.
233            $a \sqrt{2k} u_0(t,k)$ and $ a^2\sqrt{2/k} \dot{u}_0(t, k)$ 
234        k : float
235            comoving momentum $k$
236        H : float
237            Hubble rate, $H(t)$
238        a : float
239            scale factor, $a(t)$
240
241        Returns
242        -------
243        dydt : NDArray
244            an array of time derivatives of y
245
246        """
247        dydt = np.zeros(y.shape)
248
249        #real
250        dydt[0] = (H*y[0] + k/a*y[1])
251        dydt[1] = ( -H*y[1] - k/a*y[0] )
252
253        #imaginary
254        dydt[2] = (H*y[2] + k/a*y[3])
255        dydt[3] = ( -H*y[3] - k/a*y[2] )
256
257        return dydt
258    
259    def compute_green(self, k : float, uk : np.ndarray, tvac : float, teval : np.ndarray|NoneType, ind : int, atol : float=1e-30, rtol : float=1e-4) -> np.ndarray:
260        r"""
261        Compute the evolution of the Green function $G_k(t',t)$ for fixed $t'$.
262
263        The evolution is obtained by solving `green_ode` from $t=t'$ backwards until $t_{\rm vac}$, defined through $k = 10^{5/2} a(t_{\rm end})H(t_{\rm vac})$.
264        From $t_{\rm end}$ onwards, the Green function is instead computed from the mode $u_0(t, k)$. The method uses `scipy.integrate.solve_ivp`.
265
266        Parameters
267        ----------
268        k : float
269           the momentum $k$.
270        uk : NDArray
271            the mode function $u_0(t,k)$
272        tvac : float
273            the final time, $t_{\rm vac}$
274        teval : NDArray or None
275            time coordinates $t$ of $G_k(t',t)$ (`None`: same as class input `sys`) 
276        ind : int
277            index of teval corresponding to $t'$
278        atol : float
279            the absolute tolerance for `solve_ivp`
280        rtol : float
281            the relative tolerance for `solve_ivp`
282
283        Return
284        ------
285        GreenN : NDArray
286            the Green function $-k G_k(t', t)$ 
287        """
288        if len(teval)==0:
289            teval = self._t
290            
291        istart = 0
292        while teval[istart]<tvac:
293            istart+=1
294
295        # G(k, t, t) = 0 by definition
296        Aini = np.array([0, 1])
297    
298
299        # Solve the EoM for G backwards in time starting from G(k, t, t)
300        def Aode(t, y): return -self.green_ode(y, k, self._Hf(-t), self._af(-t))
301        solA = solve_ivp(Aode, [-teval[ind], -tvac], Aini, t_eval=-teval[istart:ind+1][::-1],
302                         method="RK45", atol=atol, rtol=rtol)
303
304        #For numerical stability, only solve the EoM for the Green function until t' = tvac. Afterwards, compute it directly from the vacuum modes.
305        GreenN = np.zeros(teval.shape)
306        GreenN[istart:ind+1] = solA.y[0,:][::-1]
307        GreenN[:istart] = ( (uk[ind].conjugate()*uk).imag*self._af(teval)**2 )[:istart]
308
309        return GreenN
310    
311    @staticmethod
312    def green_ode(y : np.ndarray, k : float, H : float, a : float) -> np.ndarray:
313        r"""
314        Evolve the re-scaled Green function for tensor modes in time.
315
316        The re-scaled Green function is
317        $$B_{t'}^k(t) = 2 k a(t)^2 \operatorname{Im} \left[ u_0^*(t', k) \, u_0(t,k) \right] = k G_k(t', t) \,.$$
318        
319        whose ODE is coupled to 
320
321        $$C_{t'}^k(t) = 2 a(t)^3 \operatorname{Im} \left[ u_0^*(t', k) \, {\dot{u}_0}(t,k) \right] \,.$$
322
323        Parameters
324        ----------
325        y : NDArray
326            the functions B and C
327        k : float
328            comoving momentum $k$
329        H : float
330            Hubble rate, $H(t)$
331        a : float
332            scale factor, $a(t)$
333
334        Returns
335        -------
336        dydt : NDArray
337            an array of time derivatives of A
338
339        """
340        dydt = np.zeros(y.shape)
341        dydt[0] =(2*H*y[0] + k/a*y[1])
342        dydt[1] = -k/a*y[0]
343        return dydt
344
345    
346    def _find_tinit_bd(self, init : np.ndarray, mode : str="k") -> Tuple[np.ndarray, np.ndarray]:
347        """
348        Determines the pair of $k$ and $t$ satisfying $k = 10^(5/2)a(t)H(t)$.
349
350        Depending on `mode`, `init` may be a time coordinate (`mode='t'`), $e$-folds (`mode='N'`) or momentum (`mode='k'`).
351
352        Parameters
353        ----------
354        init : array
355            the input array (t, N, or k)
356        mode : str
357            indicate the type of `init`
358
359        Returns
360        -------
361        ks : NDarray
362            an array of momenta
363        tstart : NDarray
364            an array of times
365
366        Raises
367        ------
368        KeyError
369            if `mode` is not 't, 'k' or 'N'
370        """
371
372        pwr = 5/2
373
374        t = self._t
375        def logkH(t): return np.log(self._af(t)*self._Hf(t))
376        if mode=="t":
377            tstart = init
378            logks = logkH(tstart)
379            ks = 10**(pwr)*np.exp(logks)
380
381        elif mode=="k":
382            ks = init
383
384            tstart = []
385            for k in ks:
386                ttmp  = self._t[np.searchsorted(10**(pwr)*np.exp(logkH(self._t)), k, "right")]
387                tstart.append(ttmp)
388            tstart = np.array(tstart)
389
390        elif mode=="N":
391            tstart = CubicSpline(self._N, t)(init)
392            ks = 10**pwr*np.exp(logkH(tstart))
393
394        else:
395            raise KeyError("'mode' must be 't', 'k' or 'N'")
396
397        return ks, tstart
398        
399    def compute_ptvac(self, k : np.ndarray, uk : np.ndarray) -> np.ndarray:
400        r"""
401        Compute the vacuum power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}(t,k)$.
402
403        Parameters
404        ----------
405        k : NDArray
406            momentum, $k$
407        uk : NDArray
408            vacuum tensor mode $u_0(t,k)$
409
410        Return
411        ------
412        PTvac : NDArray
413            the vacuum tensor power spectrum
414        """
415        PTvac = 2*(k*self._omega)**2/( np.pi**2 ) * abs(uk)**2
416        return PTvac
417
418    def compute_ptind(self, k : float, lgrav : float, Ngrid : np.ndarray, ind: int, GreenN : np.ndarray, pgrid : np.ndarray,
419                                l1 : float, A1 : np.ndarray, dA1 : np.ndarray, l2 : float, A2 : np.ndarray, dA2 : np.ndarray,
420                                momgrid : int=100) -> float:
421        r"""
422        Compute the vacuum power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}(t,k)$.
423
424        The integral is computed by integrating over the inner integral over the gauge-mode functions $A_\mu(s,p)$ using `scipy.integrate.trapezoid`.
425        The external momentum integrals are also computed using `scipy.integrate.trapezoid` on a grid of momenta momgrid x momgrid. The mode functions
426        are interpolated over momentum to match this grid.
427
428        Parameters
429        ----------
430        k : float
431            momentum $k$
432        lgrav : NDArray
433            helicity $\lambda$
434        Ngrid : NDArray
435            $e$-folds $\log a$
436        ind : integer
437            the power spectrum is computed at $t$ corresponding to `Ngrid[ind]`
438        GreenN : NDArray
439            the Green function $kG_k(t, s)$ with $s$ corresponding to `Ngrid`
440        pgrid : NDArray
441            momenta $p$ for the mode functions $A_\mu(s,p)$
442        l1, l2 : float
443            the gauge-field helicities $\mu_1$, $\mu_2$
444        A1, A2 : float
445            the gauge-field mode functions $\sqrt{2k} A_{\mu_1, \mu_2}(s, p)$
446        dA1, dA2 : float
447            the gauge-field mode function derivatives $a(s)\sqrt{2/k} \dot{A}_{\mu_1, \mu_2}(s, p)$
448        momgrid : int
449            the grid size for the momentum integrals.
450
451        Returns
452        -------
453        PTind : float
454            the gauge-field induced tensor power spectrum
455        """
456
457        cutUV = self._khN(Ngrid[ind])/k
458        cutIR = min(pgrid)/k
459        HN = self._HN(Ngrid)
460
461        logAs = np.linspace(np.log(max(0.5, cutIR)), np.log(cutUV), momgrid)
462
463        Ax_rad = PchipInterpolator(np.log(pgrid), abs(A1))
464        Ax_phase = PchipInterpolator(np.log(pgrid), np.arccos(A1.real/abs(A1)))
465        
466        dAx_rad = PchipInterpolator(np.log(pgrid), abs(dA1))
467        dAx_phase = PchipInterpolator(np.log(pgrid), np.arccos(dA1.real/abs(dA1)))
468
469        Ay_rad = PchipInterpolator(np.log(pgrid), abs(A2))
470        Ay_phase = PchipInterpolator(np.log(pgrid), np.arccos(A2.real/abs(A2)))
471        
472        dAy_rad = PchipInterpolator(np.log(pgrid), abs(dA2))
473        dAy_phase = PchipInterpolator(np.log(pgrid), np.arccos(dA2.real/abs(dA2)))
474
475        def Afuncx(x): return Ax_rad(x)*np.exp(1j*Ax_phase(x))
476        def dAfuncx(x): return dAx_rad(x)*np.exp(1j*dAx_phase(x))
477
478        def Afuncy(x): return Ay_rad(x)*np.exp(1j*Ay_phase(x))
479        def dAfuncy(x): return dAy_rad(x)*np.exp(1j*dAy_phase(x))
480
481        IntOuter = []
482        for logA in logAs:
483            A = np.exp(logA)
484            Blow = (cutIR - A)
485            Bhigh = (A - cutIR)
486            if Bhigh>0.5:
487                Blow = (A - cutUV)
488                Bhigh = (cutUV - A)
489                if Bhigh>0.5:
490                    Blow = -0.5
491                    Bhigh = 0.5
492            Bs = np.linspace(Blow, Bhigh, momgrid)[1:-1]
493
494            IntInner = np.zeros(Bs.shape)
495
496            for j, B in enumerate(Bs):
497                Ax = Afuncx(np.log(k*(A+B)))
498                dAx = dAfuncx(np.log(k*(A+B)))
499                Ay = Afuncy(np.log(k*(A-B)))
500                dAy = dAfuncy(np.log(k*(A-B)))
501
502                mom =  abs( l1*l2 + 2*lgrav*( (l1+l2)*A + (l1-l2)*B ) + 4*(A**2 - B**2) + 8*lgrav*A*B*( (l1-l2)*A - (l1+l2)*B ) - 16*l1*l2*A**2*B**2 )
503                z = max(A+B,A-B)
504                mask = np.where(z<self._khN(Ngrid)/k, 1, 0)
505
506                val = (dAx*dAy + l1*l2*Ax*Ay)*mask*k/(np.exp(3*Ngrid)*HN)
507                timeintegrand = GreenN*val*mom
508            
509                timeintre = trapezoid(timeintegrand[:ind].real, Ngrid[:ind])
510                
511                timeintim = trapezoid(timeintegrand[:ind].imag, Ngrid[:ind])
512                
513                IntInner[j] = (timeintre**2 + timeintim**2)*A
514
515  
516            IntOuter.append(trapezoid(IntInner, Bs))
517        IntOuter = np.array(IntOuter)
518
519        PTind = trapezoid(IntOuter, logAs) / (16*np.pi**4)*(k*self._omega)**4
520
521        return PTind

A class used to compute the tensor power spectrum including vacuum and gauge-field induced contributions.

These main method of this module is compute_pt, which computes both the vacuum and gauge-field induced contribution to the tensor power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}$ and $\mathcal{P}_{T,\lambda}^{\mathrm{ind}}$.

Results are internally computed using numerical units, but are returned in physical units.

PowSpecT(insys: geff.bgtypes.BGSystem) 
24    def __init__(self, insys : BGSystem):
25        """
26        Initialise the class from a GEF solution.
27
28        Parameters
29        ----------
30        insys : BGSystem
31            the GEF solution.
32        """
33        #Set GEF results to Hubble units.
34        sys = BGSystem.from_system(insys, True)
35        sys.units = False
36        
37        #import the background evolution
38        N = sys.N.value
39        a = np.exp(N)
40        H = sys.H.value
41
42        self._omega = sys.omega
43            
44        #Set the range of modes
45        self.maxk = max(a*H)
46        self.mink = 1e4
47
48        #Define Useful quantities
49
50        self._t = sys.t.value
51        self._N = N
52        self._H = H
53        self._xi= sys.xi.value
54
55        #define interpolated quantities as needed
56        self._af = CubicSpline(self._t, a)
57        if isinstance(self._H, Constant):
58            self._Hf = lambda t: self._H
59        else:
60            self._Hf = CubicSpline(self._t, H)
61        self._HN = CubicSpline(self._N, self._H)
62        self._khN = CubicSpline(N, sys.kh)
63
64        #Obtain eta as a functio of time for phases
65        def deta(t, y): return 1/self._af(t)
66        
67        soleta = solve_ivp(deta, [min(self._t), max(self._t)], np.array([0]), t_eval=self._t)
68
69        self._etaf = CubicSpline(self._t, soleta.y[0,:])
70        return

Initialise the class from a GEF solution.

Parameters
  • insys (BGSystem): the GEF solution.
def compute_pt( self, nmodes: int, spec: geff.mbm.GaugeSpec, FastGW: bool = True, atols: list = [0.001, 1e-20], rtols: list = [0.0001, 0.0001], momgrid: int = 100) -> Tuple[numpy.ndarray, dict]: 
 72    def compute_pt(self, nmodes : int, spec : GaugeSpec, FastGW : bool=True,
 73                    atols : list=[1e-3,1e-20], rtols : list=[1e-4,1e-4], momgrid : int=100
 74                    ) -> Tuple[np.ndarray,dict]:
 75        r"""
 76        Compute the full tensor power spectrum.
 77
 78        The method lodes data on gauge modes $A_\lambda(t,k)$ from a file indicated by `mbm_file`.
 79        The power spectrum is evaluated fpr `nmodes` log-spaced momenta $k \in [10^{4}a(0)H(0), 10 a(t_{\rm max}) H(t_{\rm max})]$
 80
 81        Parameters
 82        ----------
 83        nmodes : NDArray
 84            the number of momenta
 85        ModePath : str
 86            path to a file containing gauge-field mode functions
 87        FastGW : bool
 88            If `True`, only the expected dominant contributions to $\mathcal{P}_{T,\lambda}^{\mathrm{ind}}$ are computed
 89        atols : list
 90            absolute tolerance for `compute_homogeneous` (index 0) and `compute_green` (index 1)
 91        rtols : list
 92            relative tolerance for `compute_homogeneous` (index 0) and `compute_green` (index 1)
 93        momgrid : int
 94            passed to `compute_ptind`
 95        
 96        Returns
 97        -------
 98        ks : NDArray
 99            the momenta $k$
100        PT : dict
101            a dictionary containing all contributions to the tensor power spectrum and its total.
102        """
103
104        k = np.logspace(np.log10(self.mink), np.log10(10*self.maxk), nmodes)
105        ks, tstarts = self._find_tinit_bd(k, mode="k")
106        
107        Ngrid = spec["N"]
108        tgrid = spec["t"]
109        pgrid = spec["k"]
110
111        GaugeModes = {"+":(spec["Ap"], spec["dAp"]), "-":(spec["Am"], spec["dAm"])}
112        
113        indend = np.argmin(abs(Ngrid-max(self._N)))
114
115        PT = {"tot":[], "vac":[], "ind+,++":[], "ind+,+-":[], "ind+,--":[], "ind-,++":[], "ind-,+-":[], "ind-,--":[]}
116
117        GWpols = [("+", 1), ("-",-1)]
118
119        gaugepols=[(("+",1),("+",1)),
120                    (("+",1),("-",-1)),
121                        (("-",-1),("-",-1))]
122
123        sign = np.sign(self._xi[0])
124
125        for i, k in enumerate(ks):
126            tstart = tstarts[i]
127
128            if k > 10**(5/2)*self.maxk:
129                for key in PT.keys():
130                    PT[key].append(0)
131            else:
132                uk, _ = self.compute_homogeneous(k, tstart, tgrid, atol=atols[0], rtol=rtols[0])
133                Green = self.compute_green(k, uk, tstart, tgrid, indend, atol=atols[1], rtol=rtols[1])
134                
135                PT["vac"].append(self.compute_ptvac(k, uk[indend]))
136
137                for lgrav in GWpols:
138                    for mu in gaugepols:
139                        GWpol = lgrav[1]
140
141                        Ax = GaugeModes[mu[0][0]][0]
142                        dAx = GaugeModes[mu[0][0]][1]
143                        lx = mu[0][1]
144                        Ay = GaugeModes[mu[1][0]][0]
145                        dAy = GaugeModes[mu[1][0]][1]
146                        ly = mu[1][1]
147
148                        if FastGW and (lx != sign or ly != sign):
149                            PT[f"ind{lgrav[0]},{mu[0][0]}{mu[1][0]}"].append(0.)
150                        else:
151                            PT[f"ind{lgrav[0]},{mu[0][0]}{mu[1][0]}"].append(
152                                self.compute_ptind(k, GWpol, Ngrid, indend, Green, pgrid, lx, Ax, dAx, ly, Ay, dAy, momgrid) )
153
154        PT["tot"] = np.zeros(ks.shape)
155        for key in PT.keys():
156            PT[key] = np.array(PT[key])
157            if ("+" in key) or ("-" in key):
158                PT["tot"] += 0.5*PT[key]
159            elif key=="vac":
160                PT["tot"] += PT[key]
161
162        #this returns k with units restored. Needs to be matched by omega_GW
163        return ks*self._omega, PT

Compute the full tensor power spectrum.

The method lodes data on gauge modes $A_\lambda(t,k)$ from a file indicated by mbm_file. The power spectrum is evaluated fpr nmodes log-spaced momenta $k \in [10^{4}a(0)H(0), 10 a(t_{\rm max}) H(t_{\rm max})]$

Parameters
  • nmodes (NDArray): the number of momenta
  • ModePath (str): path to a file containing gauge-field mode functions
  • FastGW (bool): If True, only the expected dominant contributions to $\mathcal{P}_{T,\lambda}^{\mathrm{ind}}$ are computed
  • atols (list): absolute tolerance for compute_homogeneous (index 0) and compute_green (index 1)
  • rtols (list): relative tolerance for compute_homogeneous (index 0) and compute_green (index 1)
  • momgrid (int): passed to compute_ptind
Returns
  • ks (NDArray): the momenta $k$
  • PT (dict): a dictionary containing all contributions to the tensor power spectrum and its total.
def compute_homogeneous( self, k: float, tvac: float, teval: numpy.ndarray | None, atol: float = 0.001, rtol: float = 0.0001) -> Tuple[numpy.ndarray, numpy.ndarray]: 
165    def compute_homogeneous(self, k : float, tvac : float, teval : np.ndarray|NoneType, atol : float=1e-3, rtol : float=1e-4) -> Tuple[np.ndarray, np.ndarray]:
166        r"""
167        Compute the evolution of a vacuum mode starting from Bunch&ndash;Davies vacuum.
168
169        For a mode $u_0(t,k)$ with momentum $k$, initialise in Bunch&ndash;Davies when $k = 10^{5/2} a(t_{\rm vac})H(t_{\rm vac})$.
170        Its evolution is obtained by solving `tensor_mode_eq` using `scipy.integrate.solve_ivp`.
171
172        Parameters
173        ----------
174        k : float
175           the momentum $k$.
176        tvac : float
177            the initial time, $t_{\rm init}$
178        teval : NDArray or None
179            time coordinates $t$ of $u_0(t,k)$ (`None`: same as class input `sys`) 
180        atol : float
181            the absolute tolerance for `solve_ivp`
182        rtol : float
183            the relative tolerance for `solve_ivp`
184
185        Returns
186        -------
187        u : NDArray
188            the vacuum tensor modes, $\sqrt{2k}u(t, k)$
189        duk : NDArray
190            the derivative of the vacuum tensor modes, $a\sqrt{2/k}\dot{u}(t, k)$
191        """
192
193        if len(teval)==0:
194            teval = self._t
195        tend = max(teval)
196
197        #conformal time needed for relative phases
198        eta = self._etaf(teval)
199
200        istart = 0
201        while teval[istart]<tvac:
202            istart+=1
203        
204        #define the ODE for the GW modes
205        def ode(t, y): return self.tensor_mode_eq( y, k, self._Hf(t), self._af(t) )
206        
207        #Initialise the modes in Bunch Davies
208        Zini = np.array([1, -10**(-5/2), 0, -1])
209
210        sol = solve_ivp(ode, [tvac, tend], Zini, t_eval=teval[istart:], method="RK45", atol=atol, rtol=rtol)
211        if not(sol.success):
212            print("Something went wrong")
213
214        #the mode was in vacuum before tvac
215        vac = list( np.exp(-1j*eta[:istart]*k) )
216        dvac = list( -1j*np.exp(-1j*eta[:istart]*k) )
217
218        #Create an array tracking a modes evolution from Bunch Davies to late times. Ensure equal length arrays for every mode k
219        uk = np.array( vac + list( (sol.y[0,:] + 1j*sol.y[2,:])*np.exp(-1j*k*eta[istart]) ) )/self._af(teval)
220        duk = np.array( dvac + list( (sol.y[1,:] + 1j*sol.y[3,:])*np.exp(-1j*k*eta[istart]) ) )/self._af(teval)
221
222        return uk, duk

Compute the evolution of a vacuum mode starting from Bunch–Davies vacuum.

For a mode $u_0(t,k)$ with momentum $k$, initialise in Bunch–Davies when $k = 10^{5/2} a(t_{\rm vac})H(t_{\rm vac})$. Its evolution is obtained by solving tensor_mode_eq using scipy.integrate.solve_ivp.

Parameters
  • k (float): the momentum $k$.
  • tvac (float): the initial time, $t_{\rm init}$
  • teval (NDArray or None): time coordinates $t$ of $u_0(t,k)$ (None: same as class input sys)
  • atol (float): the absolute tolerance for solve_ivp
  • rtol (float): the relative tolerance for solve_ivp
Returns
  • u (NDArray): the vacuum tensor modes, $\sqrt{2k}u(t, k)$
  • duk (NDArray): the derivative of the vacuum tensor modes, $a\sqrt{2/k}\dot{u}(t, k)$
@staticmethod
def tensor_mode_eq(y: numpy.ndarray, k: float, H: float, a: float) -> numpy.ndarray: 
224    @staticmethod
225    def tensor_mode_eq(y : np.ndarray, k : float, H : float, a : float) -> np.ndarray:
226        r"""
227        Mode equation for vacuum tensor modes.
228
229        Parameters
230        ----------
231        y : NDArray
232            the vacuum tensor mode and its derivatives.
233            $a \sqrt{2k} u_0(t,k)$ and $ a^2\sqrt{2/k} \dot{u}_0(t, k)$ 
234        k : float
235            comoving momentum $k$
236        H : float
237            Hubble rate, $H(t)$
238        a : float
239            scale factor, $a(t)$
240
241        Returns
242        -------
243        dydt : NDArray
244            an array of time derivatives of y
245
246        """
247        dydt = np.zeros(y.shape)
248
249        #real
250        dydt[0] = (H*y[0] + k/a*y[1])
251        dydt[1] = ( -H*y[1] - k/a*y[0] )
252
253        #imaginary
254        dydt[2] = (H*y[2] + k/a*y[3])
255        dydt[3] = ( -H*y[3] - k/a*y[2] )
256
257        return dydt

Mode equation for vacuum tensor modes.

Parameters
  • y (NDArray): the vacuum tensor mode and its derivatives. $a \sqrt{2k} u_0(t,k)$ and $ a^2\sqrt{2/k} \dot{u}_0(t, k)$
  • k (float): comoving momentum $k$
  • H (float): Hubble rate, $H(t)$
  • a (float): scale factor, $a(t)$
Returns
  • dydt (NDArray): an array of time derivatives of y
def compute_green( self, k: float, uk: numpy.ndarray, tvac: float, teval: numpy.ndarray | None, ind: int, atol: float = 1e-30, rtol: float = 0.0001) -> numpy.ndarray: 
259    def compute_green(self, k : float, uk : np.ndarray, tvac : float, teval : np.ndarray|NoneType, ind : int, atol : float=1e-30, rtol : float=1e-4) -> np.ndarray:
260        r"""
261        Compute the evolution of the Green function $G_k(t',t)$ for fixed $t'$.
262
263        The evolution is obtained by solving `green_ode` from $t=t'$ backwards until $t_{\rm vac}$, defined through $k = 10^{5/2} a(t_{\rm end})H(t_{\rm vac})$.
264        From $t_{\rm end}$ onwards, the Green function is instead computed from the mode $u_0(t, k)$. The method uses `scipy.integrate.solve_ivp`.
265
266        Parameters
267        ----------
268        k : float
269           the momentum $k$.
270        uk : NDArray
271            the mode function $u_0(t,k)$
272        tvac : float
273            the final time, $t_{\rm vac}$
274        teval : NDArray or None
275            time coordinates $t$ of $G_k(t',t)$ (`None`: same as class input `sys`) 
276        ind : int
277            index of teval corresponding to $t'$
278        atol : float
279            the absolute tolerance for `solve_ivp`
280        rtol : float
281            the relative tolerance for `solve_ivp`
282
283        Return
284        ------
285        GreenN : NDArray
286            the Green function $-k G_k(t', t)$ 
287        """
288        if len(teval)==0:
289            teval = self._t
290            
291        istart = 0
292        while teval[istart]<tvac:
293            istart+=1
294
295        # G(k, t, t) = 0 by definition
296        Aini = np.array([0, 1])
297    
298
299        # Solve the EoM for G backwards in time starting from G(k, t, t)
300        def Aode(t, y): return -self.green_ode(y, k, self._Hf(-t), self._af(-t))
301        solA = solve_ivp(Aode, [-teval[ind], -tvac], Aini, t_eval=-teval[istart:ind+1][::-1],
302                         method="RK45", atol=atol, rtol=rtol)
303
304        #For numerical stability, only solve the EoM for the Green function until t' = tvac. Afterwards, compute it directly from the vacuum modes.
305        GreenN = np.zeros(teval.shape)
306        GreenN[istart:ind+1] = solA.y[0,:][::-1]
307        GreenN[:istart] = ( (uk[ind].conjugate()*uk).imag*self._af(teval)**2 )[:istart]
308
309        return GreenN

Compute the evolution of the Green function $G_k(t',t)$ for fixed $t'$.

The evolution is obtained by solving green_ode from $t=t'$ backwards until $t_{\rm vac}$, defined through $k = 10^{5/2} a(t_{\rm end})H(t_{\rm vac})$. From $t_{\rm end}$ onwards, the Green function is instead computed from the mode $u_0(t, k)$. The method uses scipy.integrate.solve_ivp.

Parameters
  • k (float): the momentum $k$.
  • uk (NDArray): the mode function $u_0(t,k)$
  • tvac (float): the final time, $t_{\rm vac}$
  • teval (NDArray or None): time coordinates $t$ of $G_k(t',t)$ (None: same as class input sys)
  • ind (int): index of teval corresponding to $t'$
  • atol (float): the absolute tolerance for solve_ivp
  • rtol (float): the relative tolerance for solve_ivp
Return

GreenN : NDArray the Green function $-k G_k(t', t)$

@staticmethod
def green_ode(y: numpy.ndarray, k: float, H: float, a: float) -> numpy.ndarray: 
311    @staticmethod
312    def green_ode(y : np.ndarray, k : float, H : float, a : float) -> np.ndarray:
313        r"""
314        Evolve the re-scaled Green function for tensor modes in time.
315
316        The re-scaled Green function is
317        $$B_{t'}^k(t) = 2 k a(t)^2 \operatorname{Im} \left[ u_0^*(t', k) \, u_0(t,k) \right] = k G_k(t', t) \,.$$
318        
319        whose ODE is coupled to 
320
321        $$C_{t'}^k(t) = 2 a(t)^3 \operatorname{Im} \left[ u_0^*(t', k) \, {\dot{u}_0}(t,k) \right] \,.$$
322
323        Parameters
324        ----------
325        y : NDArray
326            the functions B and C
327        k : float
328            comoving momentum $k$
329        H : float
330            Hubble rate, $H(t)$
331        a : float
332            scale factor, $a(t)$
333
334        Returns
335        -------
336        dydt : NDArray
337            an array of time derivatives of A
338
339        """
340        dydt = np.zeros(y.shape)
341        dydt[0] =(2*H*y[0] + k/a*y[1])
342        dydt[1] = -k/a*y[0]
343        return dydt

Evolve the re-scaled Green function for tensor modes in time.

The re-scaled Green function is $$B_{t'}^k(t) = 2 k a(t)^2 \operatorname{Im} \left[ u_0^*(t', k) \, u_0(t,k) \right] = k G_k(t', t) \,.$$

whose ODE is coupled to

$$C_{t'}^k(t) = 2 a(t)^3 \operatorname{Im} \left[ u_0^*(t', k) \, {\dot{u}_0}(t,k) \right] \,.$$

Parameters
  • y (NDArray): the functions B and C
  • k (float): comoving momentum $k$
  • H (float): Hubble rate, $H(t)$
  • a (float): scale factor, $a(t)$
Returns
  • dydt (NDArray): an array of time derivatives of A
def compute_ptvac(self, k: numpy.ndarray, uk: numpy.ndarray) -> numpy.ndarray: 
399    def compute_ptvac(self, k : np.ndarray, uk : np.ndarray) -> np.ndarray:
400        r"""
401        Compute the vacuum power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}(t,k)$.
402
403        Parameters
404        ----------
405        k : NDArray
406            momentum, $k$
407        uk : NDArray
408            vacuum tensor mode $u_0(t,k)$
409
410        Return
411        ------
412        PTvac : NDArray
413            the vacuum tensor power spectrum
414        """
415        PTvac = 2*(k*self._omega)**2/( np.pi**2 ) * abs(uk)**2
416        return PTvac

Compute the vacuum power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}(t,k)$.

Parameters
  • k (NDArray): momentum, $k$
  • uk (NDArray): vacuum tensor mode $u_0(t,k)$
Return

PTvac : NDArray the vacuum tensor power spectrum

def compute_ptind( self, k: float, lgrav: float, Ngrid: numpy.ndarray, ind: int, GreenN: numpy.ndarray, pgrid: numpy.ndarray, l1: float, A1: numpy.ndarray, dA1: numpy.ndarray, l2: float, A2: numpy.ndarray, dA2: numpy.ndarray, momgrid: int = 100) -> float: 
418    def compute_ptind(self, k : float, lgrav : float, Ngrid : np.ndarray, ind: int, GreenN : np.ndarray, pgrid : np.ndarray,
419                                l1 : float, A1 : np.ndarray, dA1 : np.ndarray, l2 : float, A2 : np.ndarray, dA2 : np.ndarray,
420                                momgrid : int=100) -> float:
421        r"""
422        Compute the vacuum power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}(t,k)$.
423
424        The integral is computed by integrating over the inner integral over the gauge-mode functions $A_\mu(s,p)$ using `scipy.integrate.trapezoid`.
425        The external momentum integrals are also computed using `scipy.integrate.trapezoid` on a grid of momenta momgrid x momgrid. The mode functions
426        are interpolated over momentum to match this grid.
427
428        Parameters
429        ----------
430        k : float
431            momentum $k$
432        lgrav : NDArray
433            helicity $\lambda$
434        Ngrid : NDArray
435            $e$-folds $\log a$
436        ind : integer
437            the power spectrum is computed at $t$ corresponding to `Ngrid[ind]`
438        GreenN : NDArray
439            the Green function $kG_k(t, s)$ with $s$ corresponding to `Ngrid`
440        pgrid : NDArray
441            momenta $p$ for the mode functions $A_\mu(s,p)$
442        l1, l2 : float
443            the gauge-field helicities $\mu_1$, $\mu_2$
444        A1, A2 : float
445            the gauge-field mode functions $\sqrt{2k} A_{\mu_1, \mu_2}(s, p)$
446        dA1, dA2 : float
447            the gauge-field mode function derivatives $a(s)\sqrt{2/k} \dot{A}_{\mu_1, \mu_2}(s, p)$
448        momgrid : int
449            the grid size for the momentum integrals.
450
451        Returns
452        -------
453        PTind : float
454            the gauge-field induced tensor power spectrum
455        """
456
457        cutUV = self._khN(Ngrid[ind])/k
458        cutIR = min(pgrid)/k
459        HN = self._HN(Ngrid)
460
461        logAs = np.linspace(np.log(max(0.5, cutIR)), np.log(cutUV), momgrid)
462
463        Ax_rad = PchipInterpolator(np.log(pgrid), abs(A1))
464        Ax_phase = PchipInterpolator(np.log(pgrid), np.arccos(A1.real/abs(A1)))
465        
466        dAx_rad = PchipInterpolator(np.log(pgrid), abs(dA1))
467        dAx_phase = PchipInterpolator(np.log(pgrid), np.arccos(dA1.real/abs(dA1)))
468
469        Ay_rad = PchipInterpolator(np.log(pgrid), abs(A2))
470        Ay_phase = PchipInterpolator(np.log(pgrid), np.arccos(A2.real/abs(A2)))
471        
472        dAy_rad = PchipInterpolator(np.log(pgrid), abs(dA2))
473        dAy_phase = PchipInterpolator(np.log(pgrid), np.arccos(dA2.real/abs(dA2)))
474
475        def Afuncx(x): return Ax_rad(x)*np.exp(1j*Ax_phase(x))
476        def dAfuncx(x): return dAx_rad(x)*np.exp(1j*dAx_phase(x))
477
478        def Afuncy(x): return Ay_rad(x)*np.exp(1j*Ay_phase(x))
479        def dAfuncy(x): return dAy_rad(x)*np.exp(1j*dAy_phase(x))
480
481        IntOuter = []
482        for logA in logAs:
483            A = np.exp(logA)
484            Blow = (cutIR - A)
485            Bhigh = (A - cutIR)
486            if Bhigh>0.5:
487                Blow = (A - cutUV)
488                Bhigh = (cutUV - A)
489                if Bhigh>0.5:
490                    Blow = -0.5
491                    Bhigh = 0.5
492            Bs = np.linspace(Blow, Bhigh, momgrid)[1:-1]
493
494            IntInner = np.zeros(Bs.shape)
495
496            for j, B in enumerate(Bs):
497                Ax = Afuncx(np.log(k*(A+B)))
498                dAx = dAfuncx(np.log(k*(A+B)))
499                Ay = Afuncy(np.log(k*(A-B)))
500                dAy = dAfuncy(np.log(k*(A-B)))
501
502                mom =  abs( l1*l2 + 2*lgrav*( (l1+l2)*A + (l1-l2)*B ) + 4*(A**2 - B**2) + 8*lgrav*A*B*( (l1-l2)*A - (l1+l2)*B ) - 16*l1*l2*A**2*B**2 )
503                z = max(A+B,A-B)
504                mask = np.where(z<self._khN(Ngrid)/k, 1, 0)
505
506                val = (dAx*dAy + l1*l2*Ax*Ay)*mask*k/(np.exp(3*Ngrid)*HN)
507                timeintegrand = GreenN*val*mom
508            
509                timeintre = trapezoid(timeintegrand[:ind].real, Ngrid[:ind])
510                
511                timeintim = trapezoid(timeintegrand[:ind].imag, Ngrid[:ind])
512                
513                IntInner[j] = (timeintre**2 + timeintim**2)*A
514
515  
516            IntOuter.append(trapezoid(IntInner, Bs))
517        IntOuter = np.array(IntOuter)
518
519        PTind = trapezoid(IntOuter, logAs) / (16*np.pi**4)*(k*self._omega)**4
520
521        return PTind

Compute the vacuum power spectrum, $\mathcal{P}_{T,\lambda}^{\mathrm{vac}}(t,k)$.

The integral is computed by integrating over the inner integral over the gauge-mode functions $A_\mu(s,p)$ using scipy.integrate.trapezoid. The external momentum integrals are also computed using scipy.integrate.trapezoid on a grid of momenta momgrid x momgrid. The mode functions are interpolated over momentum to match this grid.

Parameters
  • k (float): momentum $k$
  • lgrav (NDArray): helicity $\lambda$
  • Ngrid (NDArray): $e$-folds $\log a$
  • ind (integer): the power spectrum is computed at $t$ corresponding to Ngrid[ind]
  • GreenN (NDArray): the Green function $kG_k(t, s)$ with $s$ corresponding to Ngrid
  • pgrid (NDArray): momenta $p$ for the mode functions $A_\mu(s,p)$
  • l1, l2 (float): the gauge-field helicities $\mu_1$, $\mu_2$
  • A1, A2 (float): the gauge-field mode functions $\sqrt{2k} A_{\mu_1, \mu_2}(s, p)$
  • dA1, dA2 (float): the gauge-field mode function derivatives $a(s)\sqrt{2/k} \dot{A}_{\mu_1, \mu_2}(s, p)$
  • momgrid (int): the grid size for the momentum integrals.
Returns
  • PTind (float): the gauge-field induced tensor power spectrum