geff.models.fai_basic

  1r"""
  2$\newcommand{\bm}[1]{\boldsymbol{#1}}$
  3Module defining the GEF model "fai_basic" corresponding to fermionic axion inflation without a heuristic scale dependence.
  4
  5For more details on this model, see e.g., [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    * `delta` - *cumulative electric damping, $\Delta$*
 16    * `rhoChi` - *fermion energy density, $\rho_{\chi}$*
 17* static variables:
 18    * `a` - *scale factor, $a$* 
 19    * `H` - *Hubble rate, $H$* 
 20    * `ddphi` - *inflaton acceleration, $\ddot{\varphi}$*
 21    * `E`, `B`, `G` - *gauge-field expectation values, $\langle \bm{E}^2 \rangle$, $\langle \bm{B}^2 \rangle$, -$\langle \bm{E} \cdot \bm{B} \rangle$*
 22    * `xi` - *instability parameter, $\xi$* 
 23    * `sigmaE`, `sigmaB` - *electric and magnetic conductivities, $\sigma_{\rm E}$, $\sigma_{\rm B}$*
 24    * `xieff` - *effective instability parameter, $\xi_{\rm eff}$*
 25    * `s` - *electric damping parameter, $s = \sigma_{\rm E}/(2H)$*
 26* constants: 
 27    * `beta` - *coupling strength, $\beta$*
 28* functions: 
 29    * `V`,`dV` - *inflaton potential, $V(\varphi)$, and its derivative, $V_{,\varphi}(\varphi)$*
 30* gauge field: 
 31    * `GF` - *tower of gauge bilinears, $\mathcal{F}_{\mathcal X}^{(n)}$*
 32
 33The model expects the following input:
 34* `phi`, `dphi` - *initial data on the inflaton, $\varphi$, $\dot\varphi$*
 35* `rhoChi` - *initial data on the fermion energy density, $\rho_{\chi}$*
 36* `beta` - *coupling strength, $\beta$*
 37* `V`, `dV` - *potential shape, $V(\varphi)$, $V_{,\varphi}(\varphi)$*
 38
 39The model tracks the following events:
 40* end of inflation - terminate solver when $\ddot{a} < 0$
 41* negative norms - return an error when $\langle \bm{E}^2 \rangle$ or  $\langle \bm{B}^2 \rangle$ are negative 
 42"""
 43import numpy as np
 44
 45from geff.bgtypes import t, N, a, H, phi, dphi, ddphi, V, dV, E, B, G, xi, kh, beta, GF, define_var
 46from geff.solver import TerminalEvent, ErrorEvent, GEFSolver
 47from geff.mbm import ModeSolver
 48
 49from geff.utility.eom import (klein_gordon, friedmann, dlnkh_schwinger, ddelta, drhoChi, gauge_field_ode_schwinger,
 50                                        conductivities_collinear, conductivities_mixed, check_accelerated_expansion)
 51from geff.utility.boundary import boundary_fai
 52from geff.utility.general import heaviside
 53from geff.utility.mode import mode_equation_SE_no_scale, damped_bd
 54from geff._docs import generate_docs, docs_models
 55
 56name = "FAI basic"
 57"""The models name."""
 58
 59settings = {"picture":"mixed"}
 60"""The model settings.
 61
 62Possible settings are "mixed", "electric", "magnetic".
 63
 64Determines if conductivities are computed assuming collinear E&M fields 
 65("electric", "magnetic") or not ("mixed").
 66"""
 67
 68# parse settings
 69def define_conductivity():
 70    if settings["picture"]=="mixed":
 71        conductivity = conductivities_mixed
 72    elif settings["picture"]=="electric":
 73        def conductivity(a, H, E, B, G, omega):
 74            return conductivities_collinear(a, H, E, B, G, -1, omega)
 75    elif settings["picture"]=="magnetic":
 76        def conductivity(a, H, E, B, G, omega):
 77            return conductivities_collinear(a, H, E, B, G, 1, omega)
 78    else:
 79        raise KeyError(f"{settings['picture']} is an unknown choice for the setting 'picture'")
 80    return np.vectorize(conductivity)
 81
 82def interpret_settings():
 83    global conductivity
 84    conductivity = define_conductivity()
 85    return
 86
 87#Define all additional variables
 88sigmaE=define_var("sigmaE", 1, 0, "electric damping")
 89sigmaB=define_var("sigmaB", 1, 0, "magnetic damping")
 90delta=define_var("delta", 0, 0, "cumulative electric damping") 
 91xieff=define_var("xieff", 0, 0, "effective instabilty parameter") 
 92s=define_var("s", 0, 0, "electric damping parameter")
 93rhoChi=define_var("rhoChi", 4, 0, "fermion energy density")
 94
 95#Assign quantities to a dictionary, classifying them by their role:
 96quantities={
 97            "time":[t],
 98            "dynamical":[N, phi, dphi, kh, delta, rhoChi],
 99            "static":[a, H, xi, E, B, G, ddphi, sigmaE, sigmaB, xieff, s],
100            "constant":[beta],
101            "function":[V, dV],
102            "gauge":[GF]
103            }
104
105#State which variables require input for initialisation
106model_input = [phi, dphi, rhoChi, beta, V, dV]
107
108#this functions is called upon initialisation of the GEF class
109def define_units(phi, dphi, V, rhoChi):
110    #compute Hubble rate at t0
111    rhoK = 0.5*dphi**2
112    rhoV = V(phi)
113    H0 = friedmann( rhoK, rhoV, rhoChi )
114    
115    omega = H0 #Characteristic frequency is the initial Hubble rate
116    mu = 1. #Charatcterisic amplitude is the Planck mass
117    
118    return omega, mu
119
120#the new function for sys_to_yini in GEFSolver
121def initial_conditions(sys, ntr):
122    yini = np.zeros((ntr+1)*3+6)
123
124    #from the 'input' dictionary
125    yini[0] = sys.N.value
126    yini[1] = sys.phi.value
127    yini[2] = sys.dphi.value
128
129    #needs to be computed
130    sys.initialise("kh")( abs(sys.dphi)*sys.beta )
131    yini[3] = np.log(sys.kh.value)
132
133    #initialise delta and rhoChi
134    yini[4] = sys.delta.value
135    yini[5] = sys.rhoChi.value
136
137    #all gauge-field expectation values are assumed to be 0 at initialisation
138    return yini
139
140#define update_sys for GEFSolver
141def update_values(t, y, sys, atol=1e-20, rtol=1e-6):
142    #spacetime variables
143    sys.t.value = t
144    sys.N.value = y[0]
145    sys.a.value = np.exp(y[0])
146
147    #parse for convenience
148    sys.phi.value = y[1]
149    sys.dphi.value = y[2]
150    sys.kh.value = np.exp(y[3])
151    sys.delta.value = y[4]
152    sys.rhoChi.value = y[5]
153
154    #the gauge-field terms in y are not stored, save these values here
155    sys.E.value = y[6]*np.exp(4*(y[3]-y[0]))
156    sys.B.value = y[7]*np.exp(4*(y[3]-y[0]))
157    sys.G.value = y[8]*np.exp(4*(y[3]-y[0]))
158
159    #Hubble rate
160    sys.H.value = ( friedmann(0.5*sys.dphi**2, sys.V(sys.phi), 
161                                0.5*(sys.E+sys.B)*sys.omega**2, sys.rhoChi*sys.omega**2) )
162
163    #conductivities
164    sigmaE, sigmaB, ks = conductivity(sys.a.value, sys.H.value, sys.E.value,
165                                       sys.B.value, sys.G.value, sys.omega)
166
167    GlobalFerm = heaviside(np.log(ks),np.log(sys.a*sys.H))
168    sys.sigmaE.value = (GlobalFerm*sigmaE)
169    sys.sigmaB.value = (GlobalFerm*sigmaB)
170
171    #boundary term parameters
172    sys.s.value = sys.sigmaE/(2*sys.H)
173    sys.xi.value = sys.beta*(sys.dphi/(2*sys.H))
174    sys.xieff.value = sys.xi + sys.sigmaB/(2*sys.H)
175
176    #acceleration for convenience
177    sys.ddphi.value = klein_gordon(sys.dphi, sys.dV(sys.phi), sys.H, -sys.G*sys.beta*sys.omega**2)
178    return
179
180def compute_timestep(t, y, sys, atol=1e-20, rtol=1e-6):
181    dydt = np.zeros(y.shape)
182
183    #odes for N and phi
184    dydt[0] = sys.H.value
185    dydt[1] = sys.dphi.value
186    dydt[2] = sys.ddphi.value
187
188    #achieving dlnkhdt is monotonous requires some care
189    dlnkhdt = dlnkh_schwinger( sys.kh, sys.dphi, sys.ddphi, sys.beta,
190                                        0., sys.xieff, sys.s, sys.a, sys.H )
191    xieff = sys.xieff.value
192    s = sys.s.value
193    sqrtterm = np.sqrt(xieff**2 + s**2 + s)
194    r = (abs(xieff) + sqrtterm)
195    logfc = y[0] + np.log(r*dydt[0]) 
196    dlnkhdt *= heaviside(dlnkhdt,0)*heaviside(logfc,y[3]*(1-1e-5))
197    dydt[3] = dlnkhdt
198
199    #the other derivatives are straight forwards
200    dydt[4] = ddelta( sys.delta, sys.sigmaE )
201    dydt[5] = drhoChi( sys.rhoChi, sys.E, sys.G, sys.sigmaE, sys.sigmaB, sys.H )
202
203    #compute boundary terms and then the gauge-field bilinear ODEs
204    Fcol = y[6:].shape[0]//3
205    F = y[6:].reshape(Fcol,3)
206    W = boundary_fai(sys.xieff.value, sys.s.value)
207    dFdt = gauge_field_ode_schwinger( F, sys.a, sys.kh, 2*sys.H*sys.xieff,
208                                            sys.sigmaE, sys.delta,
209                                                W, dlnkhdt )
210    #reshape to fit dydt
211    dydt[6:] = dFdt.reshape(Fcol*3)
212
213    return dydt
214
215#Event 1: Track the end of inflation:
216def condition_EndOfInflation(t, y, sys):
217   #compute energy densities
218    rhoV = sys.V(y[1])
219    rhoK = y[2]**2/2
220    rhoEM = 0.5*(y[6]+y[7])*(sys.omega/sys.mu)**2*np.exp(4*(y[3]-y[0]))
221    rhoF = y[5]*(sys.omega/sys.mu)**2
222
223    #compute pressure
224    pV = -rhoV
225    pK = rhoK
226    pEM = 1/3*rhoEM
227    pF = 1/3*rhoF
228
229    return check_accelerated_expansion([rhoV, rhoK, rhoEM, rhoF], [pV, pK, pEM, pF])/(sys.H)**2
230
231def consequence_EndOfInflation(sys, occurance):
232    if occurance:
233        #stop solving once the end of inflation is reached
234        return "finish", {}
235    else:
236        #increase tend given the current Hubble rate
237        tdiff = np.round(5/sys.H, 0)
238        #round again, sometimes floats cause problems in t_span and t_eval.
239        tend  = np.round(sys.t + tdiff, 0)
240        return "proceed", {"tend":tend}
241    
242EndOfInflation = TerminalEvent("End of inflation", condition_EndOfInflation, -1, consequence_EndOfInflation)
243"""Defines the 'End of inflation' event."""
244
245#Event 2: ensure that E^2 and B^2 are positive
246def condition_NegativeNorms(t, y, sys):
247    return min(y[6], y[7])
248    
249NegativeNorms : ErrorEvent = ErrorEvent("Negative norms", condition_NegativeNorms, -1, "Negative value for E^2 or B^2.")
250"""Defines the 'Negative norms' event."""
251
252events = [EndOfInflation, NegativeNorms]
253
254
255
256#gather all information in the solver
257solver = GEFSolver(initial_conditions, update_values, compute_timestep, quantities, events)
258"""The solver used by the GEF model."""
259
260#define mode-by-mode solver
261MbM = ModeSolver(mode_equation_SE_no_scale, {"a":a, "xieff":xieff, "H":H, "sigmaE":sigmaE},
262                         damped_bd, {"a":a, "sigmaE":sigmaE, "delta":delta}, default_atol=1e-5)
263"""The mode solver used by the GEF model."""
264
265
266#define default docs for the above functions
267generate_docs(docs_models.DOCS)