geff.models.pai

  1r"""
  2$\newcommand{\bm}[1]{\boldsymbol{#1}}$
  3Defines the GEF model "pai" corresponding to pure axion inflation.
  4
  5For more details on this model, see e.g., [arXiv:2109.01651](https://arxiv.org/abs/2109.01651).
  6
  7---
  8
  9The model knows the following variables:
 10* time variable: `t` - *cosmic time*, $t$ 
 11* dynamical variables:
 12    * `N` - *$e$-folds*,  $N$
 13    * `phi`, `dphi` - *inflaton amplitude, $\varphi$, and velocity, $\dot{\varphi}$* 
 14    * `kh` -  *the instability scale, $k_{\rm h}$*
 15* static variables:
 16    * `a` - *scale factor, $a$* 
 17    * `H` - *Hubble rate, $H$* 
 18    * `ddphi` - *inflaton acceleration, $\ddot{\varphi}$*
 19    * `E`, `B`, `G` - *gauge-field expectation values, $\langle \bm{E}^2 \rangle$, $\langle \bm{B}^2 \rangle$, -$\langle \bm{E} \cdot \bm{B} \rangle$*
 20    * `xi` - *instability parameter, $\xi$* 
 21* constants: 
 22    * `beta` - *coupling strength, $\beta$*
 23* functions: 
 24    * `V`,`dV` - *inflaton potential, $V(\varphi)$, and its derivative, $V_{,\varphi}(\varphi)$*
 25* gauge field: 
 26    * `GF` - *tower of gauge bilinears, $\mathcal{F}_{\mathcal X}^{(n)}$*
 27
 28The model expects the following input:
 29* `phi`, `dphi` - *initial data on the inflaton, $\varphi$, $\dot\varphi$*
 30* `rhoChi` - *initial data on the fermion energy density, $\rho_{\chi}$*
 31* `beta` - *coupling strength, $\beta$*
 32* `V`, `dV` - *potential shape, $V(\varphi)$, $V_{,\varphi}(\varphi)$*
 33
 34The model tracks the following events:
 35* end of inflation - terminate solver when $\ddot{a} < 0$
 36* negative norms - return an error when $\langle \bm{E}^2 \rangle$ or  $\langle \bm{B}^2 \rangle$ are negative 
 37"""
 38import numpy as np
 39
 40from geff.bgtypes import t, N, a, H, phi, dphi, ddphi, V, dV, E, B, G, xi, kh, beta, GF
 41from geff.solver import TerminalEvent, ErrorEvent, GEFSolver
 42from geff.mbm import ModeSolver
 43
 44from geff.utility.eom import friedmann, gauge_field_ode, dlnkh, klein_gordon, check_accelerated_expansion
 45from geff.utility.mode import bd_classic, mode_equation_classic
 46from geff.utility.boundary import boundary_pai
 47from geff.utility.general import heaviside
 48from geff._docs import generate_docs, docs_models
 49
 50
 51name : str = "PAI"
 52"""The models name."""
 53
 54quantities : dict={
 55            "time":[t], #time coordinate according to which EoMs are expressed
 56            "dynamical":[N, phi, dphi, kh], #variables which evolve in time according to an EoM
 57            "static":[a, H, ddphi, xi, E, B, G], #variables which are derived from dynamical variables
 58            "constant":[beta], #constant quantities in the model
 59            "function":[V, dV], #functions of variables such as scalar potentials
 60            "gauge":[GF] #Gauge fields whose dynamics is given in terms of bilinear towers of expectation values
 61            }
 62
 63#State which variables require input for initialisation
 64model_input = [ beta, phi, dphi, V, dV]
 65
 66#this functions is called upon initialisation of the GEF class
 67def define_units(phi, dphi, V):
 68    #compute Hubble rate at t0
 69    rhoK = 0.5*dphi**2
 70    rhoV = V(phi)
 71    H0 = friedmann( rhoK, rhoV )
 72    
 73    omega = H0 #Characteristic frequency is the initial Hubble rate in Planck units
 74    mu = 1. #Charatcterisic amplitude is the Planck mass (in Planck units)
 75
 76    return omega, mu
 77
 78#define sys_to_yini for GEFSolver
 79def initial_conditions(sys, ntr):
 80    yini = np.zeros((ntr+1)*3+4)
 81
 82    #from the 'input' dictionary
 83    yini[0] = sys.N.value
 84    yini[1] = sys.phi.value
 85    yini[2] = sys.dphi.value
 86
 87    #needs to be computed
 88    sys.initialise("kh")( abs(sys.dphi)*sys.beta )
 89
 90    #all gauge-field expectation values are assumed to be 0 at initialisation
 91    return yini
 92
 93#define update_sys for GEFSolver
 94def update_values(t, y, sys):
 95    #spacetime variables
 96    sys.t.value = t
 97    sys.N.value = y[0]
 98    sys.a.value = np.exp(y[0])
 99
100    #parse for convenience
101    sys.phi.value = y[1]
102    sys.dphi.value =  y[2]
103    sys.kh.value = np.exp(y[3])
104
105    #the gauge-field terms in y are not stored, save these values here
106    sys.E.value = y[4]*np.exp(4*(y[3]-y[0]))
107    sys.B.value = y[5]*np.exp(4*(y[3]-y[0]))
108    sys.G.value = y[6]*np.exp(4*(y[3]-y[0]))
109
110    #Hubble rate
111    sys.H.value = friedmann(0.5*sys.dphi**2, sys.V(sys.phi), 0.5*(sys.E+sys.B)*sys.omega**2)
112
113    #boundary term parameter
114    sys.xi.value = sys.beta*(sys.dphi/(2*sys.H))
115
116    #acceleration for convenience
117    sys.ddphi.value = klein_gordon(sys.dphi, sys.dV(sys.phi), sys.H, -sys.G*sys.beta*sys.omega**2)
118    return
119
120#define timestep for GEFSolver
121def compute_timestep(t, y, sys):
122
123    dydt = np.zeros(y.shape)
124
125    #odes for N and phi
126    dydt[0] = sys.H.value
127    dydt[1] = sys.dphi.value
128    dydt[2] = sys.ddphi.value
129
130    #achieving dlnkhdt is monotonous requires some care
131    dlnkhdt = dlnkh( sys.kh, sys.dphi, sys.ddphi, sys.beta,
132                       0., sys.xi, sys.a, sys.H )
133    
134    logfc = y[0] + np.log( 2*abs(sys.xi)*dydt[0])
135    dlnkhdt *= heaviside(dlnkhdt,0)*heaviside(logfc,y[3]*(1-1e-5))
136    dydt[3] = dlnkhdt
137
138    #compute boundary terms and then the gauge-field bilinear ODEs
139    Fcol = y[4:].shape[0]//3
140    F = y[4:].reshape(Fcol,3)
141    W = boundary_pai(sys.xi.value)
142    dFdt = gauge_field_ode(F, sys.a, sys.kh, 2*sys.H*sys.xi, W, dlnkhdt, L=20)
143
144    #reshape to fit dydt
145    dydt[4:] = dFdt.reshape(Fcol*3)
146
147    return dydt
148
149#Event 1: Track the end of inflation:
150def condition_EndOfInflation(t, y, sys):
151    #compute energy densities
152    rhoV = sys.V(y[1])
153    rhoK = y[2]**2/2
154    rhoEM = 0.5*(y[4]+y[5])*(sys.omega/sys.mu)**2*np.exp(4*(y[3]-y[0]))
155
156    #compute pressure
157    pV = -rhoV
158    pK = rhoK
159    pEM = 1/3*rhoEM
160
161    return check_accelerated_expansion([rhoV, rhoK, rhoEM], [pV, pK, pEM])/(sys.H)**2
162
163def consequence_EndOfInflation(sys, occurrence):
164    if occurrence:
165        #stop solving once the end of inflation is reached
166        return "finish", {}
167    else:
168        #increase tend given the current Hubble rate
169        tdiff = np.round(5/sys.H, 0)
170        #round again, sometimes floats cause problems in t_span and t_eval.
171        tend  = np.round(sys.t + tdiff, 0)
172        return "proceed", {"tend":tend}
173
174EndOfInflation : TerminalEvent = TerminalEvent("End of inflation", condition_EndOfInflation, -1, consequence_EndOfInflation)
175"""Defines the 'End of inflation' event."""
176
177#Event 2: ensure that E^2 and B^2 are positive
178def condition_NegativeNorms(t, y, sys):
179    return min(y[4], y[5])
180    
181NegativeNorms : ErrorEvent = ErrorEvent("Negative norms", condition_NegativeNorms, -1, "Negative value for E^2 or B^2.")
182"""Defines the 'Negative norms' event."""
183
184events = [EndOfInflation, NegativeNorms]
185
186#gather all information in the solver
187solver = GEFSolver(initial_conditions, update_values, compute_timestep, quantities, events)
188"""The solver used by the GEF model."""
189
190#define mode-by-mode solver
191MbM = ModeSolver(mode_equation_classic, {"a":a,"xi":xi, "H":H}, bd_classic, {})
192"""The mode solver used by the GEF model."""
193
194
195#define default docs for the above functions
196generate_docs(docs_models.DOCS)