1r"""
  2$\newcommand{\bm}[1]{\boldsymbol{#1}}$
  3
  4Welcome to our visitors tour of the **Gradient Expansion Formalism Factory**!
  5
  6---
  7
  8# GEFF: Established in 2025
  9
 10This python package is designed to handle gauge-field production during cosmic inflation
 11using the *gradient expansion formalism* (GEF).
 12
 13If you are interested in axion inflation, the package comes with everything you need:
 14- **a variety of flavors**
 15    - pure axion inflation (PAI) : *good ol' axion inflation*
 16    - fermionic axion inflation (FAI) : *axion inflation with Standard Model fermions!*
 17- **useful tools**
 18    - Resolve the dynamics of axion inflation including homogeneous backreaction.
 19    - Analyze the gauge-field spectrum.
 20    - Determine the vacuum and induced tensor power spectrum.
 21    - Compute gravitational-wave spectra.
 22
 23
 24But we don't want to hold you back! The package provides a flexible framework to create your **own GEF flavor**, with all built-in tools at your disposable. 
 25It is indeed a true GEF *factory* !
 26
 27You can install this package using pip
 28
 29```bash
 30pip install cosmo-geff
 31```
 32
 33or using the `geff.yml` file found at the [GitHub repository](https://github.com/richard-von-eckardstein/GEFF) for this package,
 34
 35```bash
 36conda env create -f geff.yml
 37``` 
 38
 39---
 40
 41# The refreshing taste of GEF
 42
 43The GEF is a numerical technique to determine the dynamics and backreaction of gauge-fields during inflation
 44by directly evolving the time-dependent quantum expectation values of the gauge field, 
 45e.g., $\langle \bm{E}^2 \rangle$, $\langle \bm{B}^2 \rangle$, $\langle \bm{E} \cdot \bm{B} \rangle$, etc.
 46If this is the first time you encounter the GEF, here are some useful articles on the topic:
 47* [2109.01651](https://arxiv.org/abs/2109.01651)
 48* [2310.09186](https://arxiv.org/abs/2310.09186)
 49* [2408.16538](https://arxiv.org/abs/2408.16538)
 50
 51The strategy behind the GEF is to take Maxwell's equations in an expanding spacetime,
 52
 53<a name="max">$$\operatorname{div} \bm{E} = 0\, , \qquad \operatorname{div} \bm{B} = 0\, ,$$</a>
 54$$\dot{\bm{E}} + 2 H \bm{E} - \frac{1}{a}\operatorname{rot} \bm{B} + \bm{J} = 0 \, ,$$
 55$$\dot{\bm{B}}  + 2 H \bm{B} + \frac{1}{a}\operatorname{rot} \bm{E} = 0 \,$$
 56
 57and use them to formulate a tower of ODEs for the quantities
 58
 59$$ \mathcal{F}_\mathcal{E}^{(n)} = \frac{a^4}{k_{{\rm UV}}^{n+4}}\langle \bm{E} \cdot \operatorname{rot}^n \bm{E}\rangle = \int\limits_{0}^{k_{{\rm UV}}(t)}\frac{{\rm d} k}{k} \frac{a^2 k^{n+3}}{2 \pi^2 k_{{\rm UV}}^{n+4}}  \sum_{\lambda}\lambda^n |\dot{A}_\lambda(t,k)|^2\, ,$$
 60$$ \mathcal{F}_\mathcal{G}^{(n)} = -\frac{a^4}{2 k_{{\rm UV}}^{n+4}}\langle \bm{E} \cdot \operatorname{rot}^n \bm{B} + \bm{B} \cdot \operatorname{rot}^n \bm{E}\rangle = \int\limits_{0}^{k_{{\rm UV}}(t)} \frac{{\rm d} k}{k} \frac{a k^{n+4}}{2 \pi^2 k_{{\rm UV}}^{n+4}}\sum_{\lambda}\lambda^{n+1} \operatorname{Re}[\dot{A}_\lambda(t,k)A_\lambda^*(t,k)] \, ,$$
 61$$ \mathcal{F}_\mathcal{B}^{(n)} = \frac{a^4}{k_{{\rm UV}}^{n+4}}\langle \bm{B} \cdot \operatorname{rot}^n \bm{B}\rangle = \int\limits_{0}^{k_{{\rm UV}}(t)}\frac{{\rm d} k}{k} \frac{k^{n+5}}{2 \pi^{2}k_{{\rm UV}}^{n+4}} \sum_{\lambda}\lambda^n |A_\lambda(t,k)|^2 \, ,$$
 62
 63where, $k_{\rm UV}$ is a suitably chosen UV regulator which can vary with time. 
 64For completeness, we have also given the expression for $\mathcal{F}_\mathcal{X}^{(n)}$ in terms of the mode functions $A_\lambda(t,k)$.
 65
 66The ODE for the $\mathcal{F}_\mathcal{X}^{(n)}$'s are then given by
 67
 68$$\frac{\rm d}{{\rm d} t} \mathcal{F}_\mathcal{E}^{(n)} + (4+n)\frac{{\rm d} \ln k_{\rm UV}}{{\rm d} t} \mathcal{F}_\mathcal{E}^{(n)}  + 2\frac{k_{\rm UV}}{a}\mathcal{F}_\mathcal{G}^{(n+1)} + 2 \frac{a^4}{k_{\rm UV}^{n+4}} \langle \bm{J} \cdot \operatorname{rot}^n \bm{E} \rangle =  S_{\mathcal{E}}^{(n)}\, , $$
 69$$\frac{\rm d}{{\rm d} t} \mathcal{F}_\mathcal{G}^{(n)} + (4+n)\frac{{\rm d} \ln k_{\rm UV}}{{\rm d} t} \mathcal{F}_\mathcal{G}^{(n)} - \frac{k_{\rm UV}}{a}\left(\mathcal{F}_\mathcal{E}^{(n+1)} - \mathcal{F}_\mathcal{B}^{(n+1)}\right) - \frac{a^4}{k_{\rm UV}^{n+4}} \langle \bm{J} \cdot \operatorname{rot}^n \bm{B} \rangle= S_{\mathcal{G}}^{(n)}\, , $$
 70$$\frac{\rm d}{{\rm d} t} \mathcal{F}_\mathcal{B}^{(n)} + (4+n)\frac{{\rm d} \ln k_{\rm UV}}{{\rm d} t} \mathcal{F}_\mathcal{B}^{(n)} - 2\frac{k_{\rm UV}}{a}\mathcal{F}_\mathcal{G}^{(n+1)}  =  S_{\mathcal{B}}^{(n)}\, .$$
 71
 72Although these are infinitely many coupled ODE's, one can typically determine an analytical closing condition, such that they may be truncated at some order $n_{\rm tr}$.
 73
 74The ODEs for the $\mathcal{F}_\mathcal{X}^{(n)}$ can now be simply solved alongside those of the inflationary background.
 75This way, one can handle gauge-field backreaction onto the background dynamics of inflation.
 76
 77The GEFF package is designed to help the user in the process of solving these equations in the following way
 78 - pre-defined and ready-to-use [GEF flavors](#basics)
 79 - tailored [algorithm](#algorithm) to solve the inflationary background dynamics
 80 - options to [implement your own GEF flavor](#model_creation)
 81
 82 ---
 83
 84<a name="basics">
 85# Sampling the GEF flavors
 86
 87GEF models come in various flavors, some of the most intresting applications of the GEF are already implemented in this package.
 88As the first part of our tour, we explore the basic application of the GEFF code using these pre-defined models.
 89
 90## Choosing your flavor
 91
 92We start by sampling a first flavor of our choice; the model "pai":
 93
 94```python
 95from geff import compile_model
 96
 97# Create the GEF model of our choice
 98paiGEF = compile_model("pai")
 99``` 
100The object `paiGEF` is the compiled version of the model found under `.models.pai`. It defines an ODE solver 
101which needs to be initialized with information on the setup we want to study:
102The "pai" model expects information on the inflaton&ndash;vector coupling, $\beta / M_{\rm P}$, 
103initial conditions for the inflaton field, $\varphi(0)$, $\dot\varphi(0)$, and the shape of the inflaton potential, $V(\varphi)$.
104
105In this example, we configure the model to start on the slow-roll attractor of a chaotic inflation potential:
106
107$$ \varphi(0) = 15.55 M_{\rm P}, \qquad \dot{\varphi}(0) = -\sqrt{\frac{2}{3}} m M_{\rm P}, \qquad V(\varphi) = \frac{1}{2}m^2 \varphi^2$$
108
109where $m = 6.16 \times 10^{-6} M_{\rm P}$. We set the coupling to $\beta = 15$.
110
111The necessary information is passed to `paiGEF` as keyword arguments:
112```python
113import numpy as np
114
115m = 6.16e-6
116beta = 15
117
118phi = 15.55
119dphi = -np.sqrt(2/3)*m
120
121def V(x): return 0.5*m**2*x**2
122def dV(x): return m**2*x
123
124mod = paiGEF(beta=beta, phi=phi, dphi=dphi, V=V, dV=dV)
125```
126If you want to know what input is expected by the GEF model, use the `print_input` method of `paiGEF`.
127
128> **A note on units**:
129> The pre-defined GEF flavors in the GEFF package work in Planck units $M_{\rm P}=1$. 
130> From the input, the GEF determines the Hubble rate at initialization, $H_0$, (also in Planck units).
131> Internally, the numerical routines work with dimensionless quantities, e.g., $\bar{X} = X H_0^{-a} M_{\rm P}^{-b}$ with $a$ and $b$ indicating
132> how $X$ scales with an inverse timescale (e.g., $H_0$) and an energy scale (e.g., $M_{\rm P}$).
133> For example, the dimensionless inflaton velocity would be like $\dot{\bar{\varphi}} = \dot{\varphi}/(H_0 M_{\rm P})$,
134> i.e., $\dot{\varphi}$ scales like an amplitude, $\varphi \sim M_{\rm P}$, and a derivative, $\partial_t \sim H_0$.
135> Don't worry, this is happening under the hood, but if you want more details, see `.bgtypes`.
136
137## Getting a taste
138
139Now, that our model is initialized, we can solve the inflationary background evolution from these initial conditions:
140```python
141sol, spec, info = mod.run()
142```
143The `run` method returned three objects. The evolution of the background dynamics is contained in `sol`,
144the evolution on the gauge-field mode functions, $A_\lambda(t, k)$, are computed and returned as the
145object `spec`, while `info` is just a byproduct that contains full information on the ODE solution in `sol`.
146For basic applications, all information we actually want is in `sol` and `spec`.
147
148Let us focus on `sol`. It is a `BGSystem` object which we can use to access the time evolution of several important inflationary quantities defined by our "pai" model.
149For example, you can use it to make a basic plot showing the evolution of the energy densities during inflation as a function of $e$-folds
150
151```python
152import matplotlib.pyplot as plt
153
154# Plot the evolution of the inflaton amplitude as a function of e-folds:
155plt.plot(sol.N, sol.dphi**2/(6*sol.H**2)) # inflaton kinetic energy density
156plt.plot(sol.N, (sol.E + sol.B)/(6*sol.H**2)) # electromagnetic energy density
157plt.plot(sol.N, sol.V(sol.phi)/(3*sol.H**2)) # inflaton potential energy density
158
159plt.yscale("log")
160plt.ylim()
161
162plt.show()
163``` 
164
165How did we know that `sol` owns the attributes `N`,`phi`, etc.? You can find out using `sol.value_names()`.
166To print a full description of all available variables for "pai", you can use `paiGEF.print_ingredients()`. 
167Otherwise, this information can also be found at `.models.pai`. 
168If you need a brief description of any variable `X`, use `X.get_description()`.
169
170> **A note on variables:**
171> The variables encoded in a GEF model use a custom class called `Variable`. They work like a `numpy` array. 
172> Additionally, constants are realized using the class `Constant`, and work like a `float`.
173> Dimensionful functions like the inflaton potential use the `Func` class, and can be used like a regular function.
174> These classes are defined to take care of unit conversions and are collectively attatched to a `BGSystem`. If you are curious, have a look at `.bgtypes`.
175
176We did not only get back a `sol` object, but also `spec`. This object contains the gauge-field mode functions $A_\lambda(t,k)$.
177The code computes these mode functions after having solved the dynamics of the background system ($H(t)$, $\varphi(t)$ etc.),
178and uses them for estimating the convergence of the background solution. 
179We briefly discuss this algorithm in [the next section](#algorithm).
180
181The object `spec` is an instance of the `GaugeSpec` class, which defines useful methods for handling time-dependent gauge-field spectra.
182For details, see `.mbm.GaugeSpec`.
183
184> **A note on gauge-field spectra**
185> If you are only interested in the background evolution and not in the gauge-field spectrum, set `nmodes=None` in `run()`.
186> However, it is advised to compute gauge-field spectra to ensure consistency between the background evolution and the spectra.
187
188Next, let us store our GEF and mode solutions in a file:
189```python
190# some dummy paths for illustration
191gefpath = "some_gef_file.dat"
192mbmpath = "some_mbm_file.dat"
193
194sol.save_variables(gefpath)
195spec.save_spec(mbmpath)
196```
197
198The data can be restored from these files using
199```python
200from geff import GaugeSpec
201
202sol = mod.load_GEFdata(gefpath)
203spec = GaugeSpec.read_spec(mbmpath)
204```
205
206> **A note on storage:**
207> The `save_variables` method does not store information on constants or functions. So, to retrieve
208> the full information on our GEF run, we need to use an appropriately configured instance of
209> `paiGEF`. In the example above, this is achieved by reusing `mod`.
210
211## A rich palette
212
213Obtaining the inflationary dynamics is nice but it is not all the GEFF can do. For example, let's use our results to compute the corresponding gravitational-wave spectrum:
214```python
215from geff.tools import PowSpecT, omega_gw
216
217# Use sol to initialize the PowSpecT class
218pt_fai = PowSpecT(sol)
219
220# Compute the vacuum and induced power spectrum for 100 momentum modes k using spec:
221ks, pt_spec = pt_fai.compute_pt(100, spec)
222
223# from the power spectrum, we can then deduce the gravitational-wave spectrum:
224f, gwspec = omega_gw(ks, pt_spec["tot"], Nend=sol.N[-1], Hend=sol.H[-1])
225``` 
226We can just use  `sol` to extract the relevant information on the background solution.
227It's that easy!
228
229To finish this first part of our tour, let us sample a second GEF flavor, "fai_kh":
230```python
231# initialize the model
232faiGEF = compile_model("fai_kh", {"picture":"electric"})
233
234# chose initial conditions
235mod = faiGEF(beta=...)
236
237# solve the ODEs as before
238sol, spec, info = mod.run(...)
239
240...
241...
242```
243This model comes in three varieties ("pictures") corresponding to the effective treatment of fermions in the model.
244We specified the particular variety py passing a `settings` dictionary upon creating the model.
245
246For more details on the available models, see `geff.models`.
247
248<a name="algorithm">
249# On the factory floor
250
251Next on the tour of the GEF factory, we visit the production line.
252The following diagram sketches the basic algorithm behind the `run` method.
253
254```mermaid
255---
256config:
257    flowchart:
258        defaultRenderer: "elk"
259    theme: 'base'
260    themeVariables:
261        secondaryColor: '#CFCFC6'
262        tertiaryColor: '#FAFAFA'
263---
264graph TB
265    subgraph A["`**GEF Model**`"]
266        direction TB
267        St((initial data <br>from model))--> GS
268        subgraph GS["`**GEFSolver**`"]
269            direction LR
270            GS1[/initial data/] --> GS2[solve background ODEs] --> GS3[/background <br> dynamics/]
271        end
272        GS -.-> MbM
273        subgraph MbM["`**ModeSolver**`"]
274            direction LR
275            MbM1[/background <br> dynamics/] --> MbM2[compute <br>mode functions] --> MbM3[/spectrum/]
276        end
277        MbM -.-> A2[compare <br>background dynamics <br>& spectrum]
278        C[self-correction <br> using spectrum]
279        GS --> A2
280        A2 -.-> |disagreement| C -.->  GS
281    A2 --> |agreement|Fin[finalize] --> R1((background <br> dynamics))
282    Fin -.-> R2((spectrum))
283    end
284```
285
286As we have seen in the [first section](#basics), the `run` method is executed as part of a GEF Model.
287Each GEF Model consists of two components, the `GEFSolver` and the `ModeSolver`. The former determines the time-dependent inflationary background, 
288the latter computes the mode functions $A_\lambda(t,k)$.
289The two results can then be compared to eachother to assess the convergence of the background dynamics.
290If the two disagree, the GEF will attempt to self-correct using $A_\lambda(t,k)$.
291
292Note that, if you use `run(nmodes=None)`, the dotted lines in the diagram can be ignored; the background solution is immediately returned without
293computing the gauge-field spectrum.
294
295For more details on the `GEFSolver`, see `geff.solver`, while for the `ModeSolver` see `geff.mbm`.
296
297---
298
299<a name="model_creation">
300# Create your own flavor
301
302Having explored the potential of the GEFF code, you may be inclined to create your own GEF flavor.
303To help you in this process, we show how to implement an example toy model.
304
305> **Warning**: Before jumping into this section, we advise that you familiarize yourself with the GEF.
306> Also, this tutorial works best, if you have a basic understanding of the classes defined in `.bgtypes`.
307
308## The first step is the hardest
309
310First, we need to work out the mathematical formulation of our model.
311
312Let us consider the case of Abelian gauge-field production in de Sitter space ($H={\rm const.}$) by a current of the type $\bm{J} = 2 H \xi \bm{B}$, 
313where $\xi$ is a constant, which we refer to as instability parameter.The ODE tower for the gauge-field bilinears are then given by:
314
315$$\frac{\rm d}{{\rm d} t} \mathcal{F}_\mathcal{E}^{(n)} + (4+n)\frac{{\rm d} \ln k_\mathrm{h}}{{\rm d} t} \mathcal{F}_\mathcal{E}^{(n)}  + 2\frac{k_\mathrm{h}}{a}\mathcal{F}_\mathcal{G}^{(n+1)} - 4 H \xi \mathcal{F}_\mathcal{G}^{(n)}=  S_{\mathcal{E}}^{(n)}\, , $$
316$$\frac{{\rm d}}{{\rm d} t} \mathcal{F}_\mathcal{G}^{(n)} + (4+n)\frac{{\rm d} \ln k_\mathrm{h}}{{\rm d} t} \mathcal{F}_\mathcal{G}^{(n)} - \frac{k_\mathrm{h}}{a}\left(\mathcal{F}_\mathcal{E}^{(n+1)} - \mathcal{F}_\mathcal{B}^{(n+1)}\right) - 2 H \xi \mathcal{F}_\mathcal{B}^{(n)}= S_{\mathcal{G}}^{(n)}\, , $$
317$$\frac{{\rm d}}{{\rm d} t} \mathcal{F}_\mathcal{B}^{(n)} + (4+n)\frac{{\rm d} \ln k_\mathrm{h}}{{\rm d} t} \mathcal{F}_\mathcal{B}^{(n)} - 2\frac{k_\mathrm{h}}{a}\mathcal{F}_\mathcal{G}^{(n+1)}  =  S_{\mathcal{B}}^{(n)}\, .$$
318
319One can determine that a sensible regularization scale for this model is given by $k_{\rm h}(t) = 2aH\xi$. 
320
321The boundary terms $S_\mathcal{X}^{(n)}$ are a consequence of the time dependence of $k_{\rm h}$. They are expressed in terms of Whittaker functions,
322but the GEFF has a module that takes care of them. We will see this later.
323
324Once we have worked out the equations, we can turn to writing a new model file.
325
326## Only the best ingredients
327
328The first thing we should state in our model file is the name of the new GEF flavor:
329```python
330import numpy as np
331
332name = "tutorial"
333```
334> **A note on settings**: Following the model's name, you can also define a `settings` dictionary.
335> This dictionary should contain some key&ndash;value pair defining its name and the default setting.
336> The `settings` should be accompanied by a function called `interpret_settings`. This method will be called on model creation,
337> so you can use it to define how user input settings are handled. See, e.g., `.models.fai_kh` for how this is done in practice.
338
339Next, we need to define define and categorize the variables which appear in our model. This is taken care of by the functions `define_var`, `define_const`,  and `define_func` in `.bgtypes`.
340This module also contains some pre-defined variables, which are often encountered in GEF models.
341
342There are three variables which every GEF model needs to define:
343* $t$ - *cosmic time* : all ODE's are solved in terms of $t$ starting from $t=0$.
344* $N$ - *$e$-folds* : needed by most internal methods. The code expects that $N(t=0)=0$.
345* $H$ - *Hubble rate* : needed by most internal methods.
346
347Beyond these staples, some extra variables appear in our ODE's:
348* $\mathcal{F}_\mathcal{X}^{(n)}$ - *the gauge-field bilinears*
349* $k_{\rm h}$ - *the UV regulator*
350* $\xi$ -  *the instability parameter*
351
352To properly account for the gauge field, we also need to add in the following three variables:
353* $\mathcal{E}^{(0)} \equiv \langle \bm{E}^2 \rangle$ - *called `E` by the GEF*
354* $\mathcal{B}^{(0)} \equiv \langle \bm{B}^2 \rangle$ - *called `B` by the GEF*
355* $\mathcal{G}^{(0)} \equiv -\langle \bm{E} \cdot \bm{B} \rangle$ - *called `G` by the GEF*
356
357> **Note on gauge field bilinears** The time evolution of the variables $\mathcal{F}_\mathcal{X}^{(n)}$ will not be saved by the GEFF code,
358> since we are typically only interest in the quantities with $n=0$. The output `info` returned by the `run` method contains the full information on $\mathcal{F}_\mathcal{X}^{(n)}$,
359> but only the information on $n=0$ is passed to `sol` in the form of  `E`, `B` and `G`.
360
361All these variables need to be defined in our model file. This is done as follows:
362```python
363from geff.bgtypes import define_const
364
365# We make use of the fact that a lot of these variables are pre-defined:
366from geff.bgtypes import t, N, a, E, B, G, kh, GF
367
368# We also need to define some new objects:
369H = define_const("H", qu_omega=1, qu_mu=0) # Hubble rate (scales like inverse time)
370xi = define_const("xi", qu_omega=0, qu_mu=0) # instability parameter
371```
372We use `define_const` to define the constants for our model: $H$ and $\xi$.
373The Hubble rate has mass dimension one, and scales with an inverse time-scale $\omega$ as, $H = \bar{H} \omega$.
374It does not scale like an amplitude $\mu$. Hence, `qu_omega=1` and `qu_mu=0`.
375Similarly, $\xi$ is just a number, and we set `qu_omega=0` and `qu_mu=0`.
376This information needs to be passed to properly allow for unit conversions in the code.
377More information on units and scaling is given in `.bgtypes`.
378
379All the variables we have defined serve a specific purpose in our GEF model. To inform the code of this, we need to classify each of them in one of these categories:
380* **time**: a time variable (needs to be `t`)
381* **dynamical**: variables whose time-evolution is determined from a differential equation.
382* **gauge**: the gauge-field variable `GF` representing $\mathcal{F}_\mathcal{X}^{(n)}$.
383* **static**: variables whose time-evolution is computed from other variables.
384* **constant**: constants of time.
385* **function**: functions of the above.
386
387In our case, this would look as follows:
388```python
389quantities = {
390            "time":[t], # this is mandatory!
391            "dynamical":[kh], # kh is best evolved from an ODE
392            "gauge":[GF], # state the obvious
393            "static":[N, a, E, B, G], # directly computed from other variables
394            "constant":[H, xi], # we assume de-Sitter space, and xi is a constant
395            "function":[]  # our model does not need any functions
396            }
397```
398
399## Write a recipe
400
401With the variables defined, an important step towards our GEF model is already taken. 
402Next, let us set up the `GEFSolver`. (For more details, see `.solver.GEFSolver`.)
403
404Internally, the GEFF package uses `scipy.integrate.solve_ivp` to solve differential equations. 
405This requires that ODE's are formulated as $\dot{\vec{y}} = f(t, \vec{y})$ with $\vec{y}$ as a `numpy` array.
406However, the GEFF prefers the `BGSystem` class, which takes care of unit conversions.
407So, we need to define how to translate between the two.
408
409First up, we define how to interpret user input to initialize the array $\vec{y}$:
410```python
411def initial_conditions(sys, ntr):
412    yini = np.zeros(1 + 3*(ntr+1))
413    yini[0] = np.log(2*sys.xi*sys.H) #index 0 of yini is log(kh)
414
415    # initialize all F_X^n as zero at indices [1:]
416
417    return yini
418``` 
419Note how `sys` has attributes `xi` and `H`. These correspond to the variables we have defined in the previous step.
420
421Importantly, upon defining `initial_conditions`, we also make a choice; which entries in $\vec{y}$ correspond to which dynamical variable.
422In our simple toy model, there actually is no choice to be made: The GEFF package expects that the gauge-field variables are stored *after* all
423dynamical variables, as the number of the $\mathcal{F}_\mathcal{X}^{(n)}$'s will vary depending on $n_{\rm tr}$.
424Hence, $k_{\rm h}$ necessarily goes first. However, we do have the choice of evolving $\log k_{\rm h}$ instead of $k_{\rm h}$.
425
426
427Next, we write the recipe used to define our GEF ODE. The ODE evolution is computed in two steps:
4281. `update_values`: Update `sys` according to $\vec{y}(t)$.
4292. `timestep`: Compute $\dot{\vec{y}}(t)$ from `sys`.
430
431To define `update_values`, we can use all the variables which we have previously declared (they are assumed to be in numerical units):
432```python
433def update_values(t, y, sys):
434    # evolution of spacetime
435    sys.t.value = t
436    sys.N.value = sys.t*sys.H #perfect de Sitter
437
438    # watch out, you can use sys.X as an array,
439    # but to pass to other functions, its better to use sys.X.value
440    sys.a.value = np.exp(sys.N.value)
441
442    # define how kh is computed from xi:
443    sys.kh.value = np.exp(y[0]) # y[0] is log(kh)
444
445    # use that y[1] = F_E^0, y[2] = F_B^0, y[3] = F_G^0
446    rescale = (sys.kh/sys.a)**4
447    sys.E.value = rescale * y[1]
448    sys.B.value = rescale * y[2]
449    sys.G.value = rescale * y[3]
450
451    return
452```
453For `timestep`, we can make use of some pre-defined functions in the `.utility` module.
454
455```python
456from geff.utility.eom import gauge_field_ode
457from geff.utility.boundary import boundary_approx
458from geff.utility.general import heaviside
459
460def timestep(t, y, sys):
461    dydt = np.zeros_like(y)
462    
463    dlnkhdt = sys.H.value #dlnkhdt derivative
464    dydt[0] = dlnkhdt
465
466    xi = sys.xi.value
467
468    # compute boundary terms
469    W = boundary_approx(xi)
470
471    # reshape arrays to fit gauge_field_ode
472    Fcol = y[1:].shape[0]//3
473    F = y[1:].reshape(Fcol,3)
474
475    # compute the gauge-field ODEs
476    dFdt = gauge_field_ode(F, sys.a, sys.kh, 2*sys.H*xi, W, dlnkhdt)
477    # note that we can use 'a', 'kh' etc.;
478    # 'update_values' is called before 'timestep'
479
480    dydt[1:] = dFdt.reshape(Fcol*3) # reshape to fit dydt
481
482    return dydt
483``` 
484
485These are all the ingredients we need to formulate the GEF ODE's. We can combine them using the `.solver.GEFSolver` class factory:
486
487```python
488from geff.solver import GEFSolver
489
490solver = GEFSolver(initial_conditions, update_values, timestep, quantities)
491``` 
492
493> **Note**: This is just a basic `GEFSolver`. You can also define `Event` objects for a solver.
494> An `Event` will check for a certain condition while the ODEs are being solved, and can terminate the solver if the condition is met.
495> For example, you can define an `Event` to check if the end of inflation has been reached, or if some positive definite quantity has become negative.
496> A `GEFSolver` can be configured to check for any `Event` occurrences and react to them in user-specified ways.
497> For more details, see `geff.solver`.
498
499We also should define the `ModeSolver`. In this toy model, we can just use a pre-defined class. For more complex situations, use `.mbm.ModeSolver`.
500```python
501from geff.mbm import BaseModeSolver
502
503MbM = BaseModeSolver
504```
505
506> **Note on naming**: The names `solver` and `MbM` are not arbitrary. The `compile_model` function will look for objects with exactly these names.
507> Please stick to the naming convention we give here to ensure your model works as intended.
508
509## The finishing touch
510
511 The last thing we need to do is define how our new GEF model is initialized.
512
513First, we need to declare, what input our GEF model expects from the user.
514There are two constants, $H$ and $\xi$, and the user should tell us their value:
515```python
516model_input = [H, xi]
517``` 
518Our model does not require other input; $\mathcal{F}_{\mathcal{X}}^{(n)}$ is initially set to zero, and $k_h(0)$ is determined from $\xi$ and $H$. 
519
520The last step is to define the units of our GEF model based on the user input. This is achieved by the `define_units` function:
521```python
522def define_units(H):
523    # The characteristic inverse time scale is the constant Hubble rate in Planck units
524    freq = H
525    # The charateristic energy scale is the Planck mass (in Planck units)
526    amp = 1. 
527    return freq, amp
528```
529The arguments of `define_units` necessarily needs to be a subset of `model_input`. In this case, we only need the Hubble rate, `H` to define our unit system:
530The energy scale $\mu$ is the Planck mass in Planck units, while the inverse time scale is the constant Hubble rate $H$.
531
532We are finally done! We can put everything we defined above in a file, which we call "tutorial.py", and we are good to go!
533
534If all went well, you can now use your own GEF flavor just like the pre-defined ones:
535```python
536import numpy as np
537from geff import compile_model
538
539# Here, we assume you have saved your model as "tutorial.py"
540import tutorial
541
542# Pre-defined models can be initialized by strings, 
543# for your own model, intialize by passing the module itself to compile_model
544TutorialGEF = compile_model(tutorial)
545
546H = 5e-6
547xi = 5
548
549mod = TutorialGEF(H=H, xi=xi)
550
551sol, spec, info = mod.run()
552
553...
554```
555
556# Attribution
557
558If you use this software in your work, please cite:
559
560```
561von Eckardstein, R. (2025). GEFF: The Gradient Expansion Formalism Factory. Zenodo. https://doi.org/10.5281/zenodo.17356579
562
563@misc{geff,
564  author       = {von Eckardstein, Richard},
565  title        = {GEFF: The Gradient Expansion Formalism Factory},
566  month        = oct,
567  year         = 2025,
568  publisher    = {Zenodo},
569  doi          = {10.5281/zenodo.17356579},
570  url          = {https://doi.org/10.5281/zenodo.17356579},
571}
572
573von Eckardstein, R. (2025). GEFF: The Gradient Expansion Formalism Factory - A tool for inflationary gauge-field production. arXiv. https://arxiv.org/abs/2510.12644
574
575@misc{vonEckardstein:2025jug,
576    author        = {von Eckardstein, Richard},
577    title         = {GEFF: The Gradient Expansion Formalism Factory - A tool for inflationary gauge-field production},
578    eprint        = {2510.12644},
579    archivePrefix = {arXiv},
580    primaryClass  = {astro-ph.CO},
581    reportNumber  = {MS-TP-25-37},
582    month         = oct,
583    year          = 2025
584}
585```
586<script type="module">
587  import mermaid from 'https://cdn.jsdelivr.net/npm/mermaid@11/dist/mermaid.esm.min.mjs';
588  import elkLayouts from 'https://cdn.jsdelivr.net/npm/@mermaid-js/layout-elk@0/dist/mermaid-layout-elk.esm.min.mjs';
589  mermaid.registerLayoutLoaders(elkLayouts);
590</script>
591"""
592from .gef import compile_model
593from .mbm import GaugeSpec
594from .bgtypes import BGSystem
595
596__version__ = "0.1.1"