geff.models.fai_kh

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