MiMeS: Misalignment Mechanism Solver

We introduce a C++ header-only library that is used to solve the axion equation of motion, MiMeS. MiMeS makes no assumptions regarding the cosmology and the mass of the axion, which allows the user to consider various cosmological scenarios and axion-like models. MiMeS also includes a convenient python interface that allows the library to be called without writing any code in C++, with minimal overhead.


Introduction
The axion is a hypothetical particle that was originally introduced in order to solve the strong CPproblem of the standard model (SM) [1,2,3]. Furthermore, the axion is assumed to acquire a non-zero vacuum expectation value (VEV), which results in the spontaneous breaking of a global symmetry, called Peccei-Quinn (PQ). That is the axion is a pseudo-Nambu-Goldstone boson. Moreover, it also appears to be a valid dark matter (DM) candidate [4,5,6], as it is electrically neutral and long-lived. The axion starts, at very early time, with a VEV close to the PQ breaking scale (usually much higher than 100GeV), but due to the expansion of the Universe it eventually ends oscillating around zero. This oscillation results in an apparent constant number of axion particles (or constant energy of the axion field) today, which can account for the observed [7] DM relic abundance. Apart form the axion, there are other hypothetical particles (for early examples of such particles see [8,9]; and [10] for a review), called axion-like particles (ALPs), which are not related to the strong CP-problem. That is, these ALPs, interact with the SM in a different way than the axion. However, they can still account for the DM of the Universe, and both axions and ALPs, follow the similar dynamics during the early universe; the follow the same equation of motion (EOM), with a different mass. The mass of the axion is dictated by QCD, while it is different (model specific) for ALPs. Therefore, the cosmological evolution of both axions and ALPs is often discussed together (see, for example ref. [11]). Moreover, deviation from the standard cosmological evolution is possible, as long as any nonstandard contribution to the energy density is absent before Big Bang Nucleosynthesis [12,13] becomes active; for temperatures [14,15,16,17] T O(10) MeV. For the QCD axion, such studies have been performed [18,19]. However, updated experimental data or new ALP models may require more similar systematic studies. To the best of our knowledge, a library for the calculation of the axion (or ALP) relic abundance does not exist. Therefore, in this article, we introduce MiMeS; a header-only C++ library, 1 that solves the axion (or ALP) EOM, where both the mass and the underlying cosmological evolution are treated as user inputs. That is, MiMeS can be used to compute the relic abundance of axions or ALPs, in a wide variety of scenarios.
There are several advantages of having a library that can calculate the relic abundance -in principle -fast. For example, one can perform a scan over different cosmological scenarios for various ALPs cases, automatically, writing only a few lines of code. Furthermore, the availability of a tool can help the community reproduce published results, without spending time duplicating the overall effort. MiMeS aims to be a useful addition in the list of available computational tools for physics, 2 since it is designed to be simple to use, simple to understand, and simple to modify.
This article is organized as follows: In section 2, we introduce the EOM and show how it can be solved approximately. Also, we introduce the notation that MiMeS follows, we derive the "adiabatic invariant" of the system, and discuss how MiMeS uses it. In the next section, we introduce MiMeS, by showing how it can be downloaded, compiled, and run for the first time. Also, we explain in detail all the parameters that MiMeS as a user input at run-time. In section 4, we discuss the few assumption that MiMeS makes, its default compile-time options and user input, and how the user can change them. Moreover, we provide a complete example in both C++ and python.

Physics background
Although there are several works in the literature (such as [24,11]) that can provide an insight on the cosmological evolution of axions, in this section we define, derive, and discuss various quantities we need, in order to understand how MiMeS works in detail.
The EOM The axion field, A, is usually expressed in terms of the so-called axion angle, θ, as A = θ f a , with f a the scale at which the PQ symmetry breaks. 3 The axion angle follows the EOM with χ a function of the temperature. For the QCD axion, this has been calculated using lattice simulations in [25]. MiMeS comes with the data provided by ref. [25]. However, the user is free to change them, or use another function for the mass.
Initial conditions We assume that at very high temperaturesm a → 0 (for the axion this is true; see Fig. (1))i.e.m a H -at very early times. Therefore, at the very early Universe, the EOM is 4 which is solved by θ = θ ini + C t 0 dt a(t = 0) a(t ) 3 , with θ ini a constant (some initial value) and a the scale factor of the Universe. That is, as the Universe expands, and as long asm a H, θ → θ ini , and 1 MiMeS is distributed under the MIT license. A copy of this license should be available in the MiMeS root directory. If you have not received one, you can find it at github.com/dkaramit/MiMeS/blob/master/LICENSE. 2 Especially for dark matter, where tools for both thermal and non-thermal production have been developed [20,21,22,23]; however, without the addition of the misalignment mechanism. 3 If in a model under study, there is no fa, this parameter is still expected by MiMeS, but the user can set fa = 1. 4 For the QCD axion the PQ symmetry breaking scale determines whether there are domains with different θini. However, one can still use a mean value of the angle as its initial condition. This is discussed in the literature; e.g. in [4,26,18]. The mass of the axion as a function of the temperature for f a = 10 12 GeV, using the data provided in ref. [25].
θ → 0. Since we would like to calculate the angle today, we can integrate eq. (2.1) from a time after inflation (call it t = t ini ) such thatθ| t=t ini = 0 and θ| t=t ini = θ ini . This is the most common case (there are exceptions to this; e.g. [27]), and it is what MiMeS uses.

The WKB approximation
In order to solve analytically eq. (2.1), we assume θ 1, which results in the linearised EOM Using a trial solution θ trial = exp i dt ψ(t) + 3/2 i H(t) , and defining Ω 2 =m a 2 − 9 4 H 2 − 3 2Ḣ we can transform the eq. (2.4) to ψ 2 = Ω 2 + iψ , (2.5) which has a formal solution ψ = ± Ω 2 + iψ. In the WKB approximation, we assume a slow timedependence; that isψ Ω 2 andΩ Ω 2 . Then we can approximate ψ as Applying, then, the initial conditionsθ| t=t ini = 0 and θ| t=t ini ≈ θ ini , we arrive at dt Ω(t ) . (2.8) In order to further simplify this approximate result, we note that θ deviates from θ ini close to t = t osc -corresponding to T = T osc , the so-called "oscillation temperature" -m a | t=tosc = 3H| t=tosc , which is defined as the point at which the axion begins to oscillate. This observation allows us to set t ini = t osc . Moreover, at t > t osc , we approximate Ω ≈m a , as H 2 andḢ become much smaller thañ m a 2 quickly after t = t osc . Finally, the axion angle takes the form where θ osc = θ| t=tosc . This equation is further simplified if we assume that θ osc ≈ θ ini , i.e. (2.10) It is worth mentioning that the accuracy of this approximation depends, in general, on T osc ; it determines the difference between θ ini and θ osc , the deviation ofθ| t=tosc from 0, and whetherΩ Ω 2 .
Axion energy density In the small angle approximation, the energy density of the axion is For the relic abundance of axions, we need to calculate their energy density at very late times. That is,ṁ a = 0,m a H andḢ H 2 . After some algebra, we obtain the approximate form of the energy density as a function of the scale factor 12) which shows that the energy density of axions at late times scales as the energy density of matter; i.e. the number of axion particles is conserved. If there is a period of entropy injection to the plasma for T < T osc , the axion energy density gets diluted, since with γ the amount of entropy injection to the plasma between t osc and t. Therefore, the present (at T = T 0 ) energy density of the axion, becomes with m a the mass of the axion at T = T 0 . Notice that the explicit dependence on f a cancels if m a = χ(T )/f a . That is, f a only affects the energy density of the axions through its impact on T osc .

Notation
The EOM (2.1) depends on time, which is not useful variable in cosmology, especially non-standard cosmologies. Therefore, we introduce The EOM in terms of u, then, becomes Notice that in a radiation dominated Universe In a general cosmological setting, if the expansion rate is dominated by an energy density that scales as ρ ∼ a −c , d log H 2 du = −c. We notice that close to rapid particle annihilations and decays, the evolution of the energy densities change, and d log H 2 du can only be computed numerically.
Moreover, it is worth mentioning that the EOM 2.17 with the initial condition θ(u = 0) = θ ini and dθ/du(u = 0) = 0 can be written as a system of first order ordinary differential equations This form of the EOM is what MiMeS uses, as it is suitable for integration via Runge-Kutta (RK) methods -which are briefly discussed in Appendix A.

Adiabatic invariant and the anharmonic factor
The EOM 2.1 ca be solved analytically in the approximation θ 1. Moreover, even for θ 1, the WKB approximation fails to capture the dynamics before the adiabatic conditions are met, and result in an inaccurate axion relic abundance. Therefore, a numerical integration should be preferred. Furthermore, in order to reduce the computation time, the numerical integration needs to stop as soon as the axion begins to evolve adiabatically. After this point, we can correlate its energy density at later times, using an "adiabatic invariant", which can be defined for oscillatory systems with varying period.
Definition of the adiabatic invariant Given a system with Hamiltonian H(θ, p; t), the equations of motion areṗ Moreover, we note that dH =θ dp −ṗ dθ + ∂H ∂t dt . (2.20) If this system exhibits closed orbits (e.g. if it oscillates), we define where the integral is over a closed path (e.g. a period, T ), and C indicates that J can always be rescaled with a constant. This quantity is the adiabatic invariant of the system, if the Hamiltonian varies slowly during a cycle. That is, Application to the axion The Hamiltonian that results in the EOM of eq. (2.1) is with Notice that the Hamiltonian varies slowly ifṁ a (T )/m a m a and H m a , which are the adiabatic conditions. When these conditions are met, the adiabatic invariant for this system becomes (2.25) where we note that θ peak denotes the maximum of θ -the peak of the oscillation, which corresponds to p = 0. That is, H(θ peak ) = V (θ peak ) a 3 . Therefore, the adiabatic invariant, takes the form where, for the last equality. we have used the adiabatic conditions, i.e. negligible change ofm a and a during one period. Usually, the adiabatic invariant is written as [28,29] J = a 3m a θ peak 2 f (θ peak ) , is called the anharmonic factor, with 0.5 f (θ peak ) ≤ 1 (see Fig. (2)).
The role of the adiabatic invariant in the axion relic energy density The adiabatic invariant allows us to calculate the maximum value of the angle θ at late times from its corresponding value at some point just after the adiabatic conditions were fulfilled. In order to do this, we can numerically integrate eq. (2.1), and identify the maxima of θ. Once the adiabatic conditions are fulfilled, we can stop the integration at a peak, θ peak, * -which corresponds to T = T * and a = a * . Then, the value of the maximum angle today (θ peak,0 1) is related to θ peak, * via Plugging this into eq. (2.11) (withθ = 0, i.e. at today's peak), we arrive at the energy density today where γ is the entropy injection coefficient between T * and T 0 , defined from Notice that eq. (2.30) is similar to the corresponding WKB result (2.14) at t osc → t * , multiplied by the anharmonic factor f (θ peak, * ). So, the numerical integration is needed in order to correctly identify T * and θ peak, * , which are greatly affected by the underlying cosmology (especially in cases where entropy injection is active close to T osc ) and whether θ ini 1 [29,19].
It is worth mentioning that MiMeS identifies the maxima in real time. Then, integration stops as soon as J becomes almost constant. That is, adiabaticity is assumed to be reached once J does not change significantly between consecutive peaks. MiMeS expects from the user to decide for how many peaks and how little J needs to change, in order for the system to be considered adiabatic. More details are given in the next section, where we discuss hw MiMeS is used.

MiMeS usage
The latest stable version of MiMeS is available at github.com/dkaramit/MiMeS/tree/stable, which can also be obtained by running: It is important to note that MiMeS relies on NaBBODES [30] and SimpleSplines [31]. These are two libraries developed independently, and in order to get MiMeS with the latest version of these libraries, one needs to run following commands 1 git clone https://github.com/dkaramit/MiMeS.git This downloads the master branch of MiMeS; and NaBBODES [30] and SimpleSplines [31] as submodules. This guaranties that MiMeS uses the most updated version of these libraries, although it may not be stable.
Once everything is downloaded successfully, we can go inside the MiMeS directory, and run "bash configure.sh" and "make". The bash script configure.sh, just writes some paths to some files, formats the data files provided in an acceptable format (in section 4.2 the format is explained), and makes some directories. The makefile is responsible for compiling some examples and checks, as well as the shared libraries that needed for the python interface. If everything runs successfully, there should be two new directories exec and lib. Inside exec, there are several executables that ready to run, in order to ensure that the code runs (e.g. no segmentation fault occurs). For example, MiMeS/exec/AxionSolve_check.run, should print the values of the parameters θ ini and f a , the oscillation temperature and the corresponding value of θ, the evolution of the axion (e.g. temperature, θ, ρ a , etc.), and the values of various quantities on the peaks of the oscillation. In the directory lib, there are several shared libraries for the python interface.
Although there are various options available at compile-time, we first discuss how MiMeS can be used, in order for the role of these options to be clear.

First steps
There are several examples in C++ (MiMeS/UserSpace/Cpp) and python (MiMeS/UserSpace/Python), as well as jupyter notebooks (MiMeS/UserSpace/JupyterNotebooks), that show in detail how MiMeS can be used. Here, we discuss the various functions one may use as it will provide some insight for the following discussions.

Using MiMeS in C++
The class that is responsible for the solution of the EOM is mimes::Axion<LD,Solver,Method>, located in MiMeS/src/Axion/AxionSolve.hpp. However, in order to use it, we first have to define the mass of the axion as a function of the temperature and f a . The axion mass is defined as an instance of the mimes::AxionMass<LD>, which is defined in the header file MiMeS/src/AxionMass/AxionMass.hpp. We should note that LD is the numerical type to be used (e.g. long double). The other template arguments are related to the differential equation solver, and their role will be explained in later sections.
In order to start, at the top of the main .cpp file, we need to include 1 #include "src/AxionMass/AxionMass.hpp" 2 #include "src/Axion/AxionSolve.hpp" Notice that if the this .cpp file is not in the root directory of MiMeS, we need to compile it using the flag -Ipath-to-root, "path-to-root" the relative path to the root directory of MiMeS; e.g. if the .cpp is in the MiMeS/UserSpace/Cpp/Axion directory, this flag should be -I../../../.
The mass of the axion can be defined in two ways; either via a data file or as a user defined function.
Axion mass form data file In many cases, axion mass cannot be written in a closed form. In these cases, the mass is assumed to be of the form of eq. (2.2), and the user has to provide a data file that tabulates the function χ(T ). Then, the axion mass can be defined as 1 mimes::AxionMass<LD> axionMass(chi_PATH, minT, maxT); The template parameter LD is a numeric type (e.g. double or long double). The argument chi_PATH is a (relative or absolute) path to a file with two columns; T (in GeV) and χ (in GeV 4 ), with increasing T . In this case, the axion mass is interpolated between the temperatures minT and maxT. These two parameters are just suggestions, and the actual interpolation limits, TMin and TMax, are chosen as the closest temperatures in the data file. That is, in general, TMin ≥ minT and TMax ≤ maxT. Beyond these limits, the mass is assumed to constant by default. However, one can change this by using with ma2\_MAX and ma2\_MIN the axion mass squared as functions (or any other callable objects) with signatures LD ma2_MAX(LD T, LD fa) and LD ma2_MIN(LD T, LD fa). These functions are called for T ≥TMax and T ≤TMin, respectively. Usually, in order to ensure continuity of the axion mass, one needs to know TMin, TMax, χ(TMin), and χ(TMax); which can be found by calling axionMass.getTMin(), axionMass.getTMax(), axionMass.getChiMin(), and axionMass.getChiMax(), respectively.
Axion mass function In some cases, the dependence of the axion mass on the temperature can be expressed analytically. If this is the case, the user can define the axion mass via with ma2 the axion mass squared as a callable object with signature LD ma2(LD T, LD fa).
The EOM solver Once the axion mass is defined, we can declare a variable that will be used to solve the EOM, as Here, LD should be the numeric type to be used; it is recommended to use long double, but other choices are also available as we discuss later. Moreover Solver and Method depend on the type of Runge-Kutta (RK) the user chooses. The available choices are shown in table5. A discussion on RK methods is given in Appendix A. The various parameters are as follows: 1. theta_i: Initial angle.
3. umax : If u >umax the integration stops (remember that u = log(a/a i )). Typically, this should be a large number (∼ 1000), in order to avoid stopping the integration before the axion begins to evolve adiabatically.

TSTOP:
If the temperature drops below this, integration stops. In most cases this should be around 10 −4 GeV, in order to be sure that any entropy injection has stopped before integration stops (since BBN bounds [12,13] should not be violated).

ratio_ini:
Integration starts when 3H/m a ≈ratio_ini (the exact point depends on the file "inputFile", which we will see later).
6. N_convergence_max and convergence_lim: Integration stops when the relative difference between two consecutive peaks is less than convergence_lim for N_convergence_max consecutive peaks. This is the point beyond which adiabatic evolution is assumed.
7. inputFile: Relative (or absolute) path to a file that describes the cosmology. the columns should be: u T [GeV] log H, sorted so that u increases. 6 It is important to remember that MiMeS assumes that the entropy injection has stopped before the lowest temperature given in inputFile. Since MiMeS is unable to guess the cosmology beyond what is given in this file, the user has to make sure that there are data between the initial temperature (which corresponds to ratio_ini), and TSTOP.
12. absolute_tolerance (optional): Absolute tolerance of the RK solver 13. relative_tolerance (optional): Relative tolerance of the RK solver. Generally, both absolute and relative tolerances should be 10 −8 . In some cases, however, one may need more accurate result (e.g. if f_a is extremely high, the oscillations happen violently, and the system destabilizes). In any case, if the tolerances are below 10 −8 , LD should be long double. MiMeS by default uses long double variables, in order to change it see the options available in section 4.2.
14. beta (optional): Controls how agreesive the adaptation is. Generally, it should be around but less than 1.
15. fac_max, fac_min (optional): The stepsize does not increase more than fac_max, and less than fac_min. This ensures a better stability. Ideally, fac_max= ∞ and fac_min= 0, but in reality one must tweak them in order to avoid instabilities. 16. maximum_No_steps (optional): Maximum steps the solver can take. Quits if this number is reached even if integration is not finished.
In order to understand the role of the optional parameters, some basic techniques of Runge-Kutta methods are discussed in Appendix A.
The EOM (2.17), then can be solved using Once the EOM is solved, we can access T osc , θ osc , and Ωh 2 via ax.T_osc, ax.theta_osc, and ax.relic. The entire evolution (the points the integrator took) of the axion angle is stored in ax.points, which is a two-dimensional std::vector<LD>, with the columns corresponding to a, T [GeV], θ, dθ/du, ρ a . Moreover, the peaks of he oscillation are stored in another two-dimensional std::vector<LD>, with the columns corresponding to a, T [GeV], θ peak , dθ/du = 0, ρ a , J. We should note that the peaks are identified using linear interpolation between integration points, in order to ensure that dθ/du = 0. That is, the values stored in ax.peaks do not exist in ax.points. Finally, the local errors (see Appendix A) of the integration points of θ and ζ are stored in ax.dtheta and ax.dzeta.
Changing axion mass definition The axion mass definiton can change at any time by changing the data file (or ma2_MIN and ma2_MAX) or using using This function remakes the interpolations, clears all vector, and sets all variables to 0.
Changing initial condition The final member function is mimes::Axion::setTheta_i, which allows the user to set a different θ ini without generating another instance. 7 This function is used as where new_theta_ini is the new value of θ ini . Running this function resets all variables to 0 (except T_osc and a_osc, since they should not change), and clears all std::vector<LD> variables, which allows the user to simply run "ax.solveAxion();" as if ax was a freshly defined instance.

Using MiMeS in python
The modules for the python interface are located in MiMeS/src/interfacePy. Although the usage of the AxionMass and Axion classes is similar to the C++ case, it is worth showing explicitly how the python interface works. One should keep in mind that the various template arguments discussed in the C++ case have to be chosen at compile-time. That is, for the python interface, one needs to choose the numeric type, and RK method to be used when the shared libraries are compiled. This is done by assigning the relevant variable in MiMeS/Definitions.mk before running "make". The various options are discussed in section 4.1, and outlined in table 6.
The two relevant classes are defined in the modules interfacePy.AxionMass and interfacePy.Axion, and can be loaded in a python script as 1 from sys import path as sysPath 2 sysPath.append('path_to_src') 3 from interfacePy.AxionMass import AxionMass 4 from interfacePy.Axion import Axion It is important that 'path_to_src' provides the relative path to the MiMeS/src directory. For example, if the script is located in MiMeS/UserSpace/Python, 'path_to_src' should be '../../src'.
Axion mass definition via a data file As before, we first need to define the axion mass. In order to define the axion mass via a file, we use 1 AxionMass axionMass(chi_PATH, minT, maxT) Here, the constructor requires the same parameters as in C++. Moreover, the axion mass beyond the interpolation limits can be changed via Although the naming is the same as in the C++ case, there is an important difference. Namely, ma2_MAX and ma2_MIN have to be functions (that take T and f a as arguments and returnm a 2 ), and cannot be any other callable object. The reason is that MiMeS uses the ctypes module, which only works with objects compatible with C. Moreover, the values of TMin, TMax, χ(TMin), and χ(TMax) can be obtained by axionMass.getTMin(), axionMass.getTMax(), axionMass.getChiMin(), and axionMass.getChiMax(), respectively.
Axion mass definition via a function Again this can be done as 1 AxionMass axionMass(ma2) with ma2 the axion mass squared, which should be a function (not any callable object) of T and f a .
Importand note Once an AxionMass is no longer required, the destructor must be called. In this case, we can run 1 del axionMass The reason is that MiMeS constructs a pointer for every instance of the class, which needs to be deleted manually.
Moreover, in order to get the various quantities on the peaks of the oscillation, we can run 1

ax.getPeaks()
This makes numpy arrays that contain the scale factor (ax.a_peak), temperature (ax.T_peak), θ (ax.theta_peak), its derivative with respect to u (ax.zeta_peak, which should be equal to 0), the energy density of the axion (ax.rho_axion_peak), and the values of the adiabatic invariant on the peaks (ax.adiabatic_invariant). Moreover, we can run the following 1 ax.getErrors() in order to store the local errors (see Appendix A) of θ and ζ in ax.dtheta and ax.dzeta, respectively.
Changing axion mass definition We can change the axion mass by changing the file, ma2_MIN and ma2_MAX, or using with ma2 a function that takes T and f a , and returns them a 2 . As in the C++ case, if the definition of the mass of the axion changes (including the definitions beyond the interpolation limits), we have to call 1 ax.restart() in order to remake the interpolation of cosmological quantities and reset the various variables.
Changing the initial condition The initial condition θ ini can be changed using 1 ax.setTeta_i (new_theta_ini) which is faster than running the constructor again, since all the interpolations are reused. However, running this function, erases all the arrays, and resets all variables to 0 (except T_osc and a_osc, as they should not change).
Importand note The Axion class constructs a pointer to an instance of the underlying mimes::Axion class, which has to be manually deleted. Therefore, once ax is used it must be deleted, i.e. we need to run 1 del ax We should note that this must run even if we assign another instance to the same variable ax, otherwise we risk running out of memory.

Assumptions and user input
MiMeS only makes a few, fairly general, assumptions. First of all, it is assumed that the axion energy density is always subdominant compared to radiation or any other component of the Universe, and that decays and annihilations of particles have a negligible effect on the axion energy density. Moreover, the initial condition is always assumed to be θ t=t ini = θ ini andθ| t=t ini = 0. Furthermore, it is also assumed that 3H/m a increases monotonically at high temperatures. Also, it is assumed that the entropy of the plasma resumes its conserved status at a temperature higher than the minimum temperature of inputFile (which is required by the constructor of the mimes::Axion<LD,Solver,Method> class). Finally, MiMeS does not try to predict anything regarding the cosmology. Therefore, the temperatures in inputFile must cover the entire region of integration; i.e. the maximum temperature has to be larger than the one required to reach ratio_ini, while the minimum one should be lower than TSTOP.

Options at Compile-time
The user has a number of options regarding different aspects for the code. If MiMeS is used without using the available makefiles, then they must use the correct values for the various template arguments, If omitted, the type of the numeric values is double (double precision). On the other hand, if LONG=long, the type is long double. Generally, using double should be enough. For the sake of numerical stability, however, it is advised to always use LONG=long, as it a safer option. The reason is that the axion angle redshifts, and can become very small, which introduces "rounding errors". Moreover, if the parameters absolute_tolerance or absolute_tolerance are chosen to be below ∼ 10 −8 , then double precision numbers may not be enough, and LONG=long is preferable. This choice comes at the cost of speed; double precision operations are usually preformed much faster. It is important to note that LONG defines a macro with the same name (in the C++ examples), which then defines the macro (again in the C++ examples) as #define LD LONG double. The macro LD, then is used as the corresponding template argument in the various classes. We point out again that if one chooses not to use the makefile files, the template arguments need to be known at compile-time. So the user has to define them in the code.
3. LONGpy: the same as LONG, but for the python interface. One should keep in mind thtat this cannot be changed inside python scripts. It just instructs ctypes what numeric type to use. Since the preferred way to compile the shared libraries is via running "make" in the root directory of MiMeS, this variable needs to be defined inside MiMeS/Definitions.mk. By default, this variable is set to long, since this is the most stable choice in general.

SOLVER:
MiMeS uses the ordinary differential equation (ODE) integrators of ref. [30]. Currently, there are two available choices; Rosenbrock and RKF. The former is a general embedded Rosenbrock implementation and it is used if SOLVER=1, while the latter is a general explicit embedded Runge-Kutta implementation and can be chosen by using SOLVER=2 (a brief description of how these algorithms are implemented can be found in Appendix A). By default inside the Definitions.mk files SOLVER=1, because the axion EOM tends to oscillate rapidly. However, in some cases, a high order explicit method may also work. Note that this variable defines a macro that is then used as the second template argument of the mimes::Axion<LD,Solver,Method> class. The preferred way to do it in the shared libraries is via the MiMeS/Definitions.mk file, however, the user if free to compile everything in a different way. In this case, the various Definitions.mk files, are not being used, and the user must define the relevant arguments in the code where MiMeS is used.

METHOD:
Depending on the type of solver, there are some available methods. 8 • For SOLVER=1, the available methods are METHOD=RODASPR2 and METHOD=ROS34PW2. The RODASPR2 choice is a fourth order Rosenbrock-Wanner method (more information can be found in ref. [33]). The ROS34PW2 choice corresponds to a third order Rosenbrock-Wanner method [34].
• For SOLVER=2, the only reliable method available in NaBBODES is the Dormand-Prince [35] chosen if METHOD=DormandPrince, which is an explicit Runge-Kutta method of seventh order.
This variable defines a macro (with the same name) that is passed as the third template parameter of mimes::Axion<LD,Solver,Method> (i.e. METHOD<LD> in the place of Method). If the compilation is not done via the makefile files, the user must define the relevant template arguments in the code.
6. CC: The C++ compiler that one chooses to use. The default option is CC=g++, which is the GNU C++ compiler, and is available for all systems. Another option is to use the clang compiler, which is chosen by CC=clang -lstdc++. MiMeS is mostly tested using g++, but clang also seems to work (and the resulting executables are sometimes faster), but the user has to make sure that their version of the compiler of choice supports the C++ 17 standard, otherwise MiMeS probably will not work.

7.
OPT: Optimization level of the compiler. By default, this is OPT=O3, which produces executables that are marginally faster than OPT=O1 and OPT=O2, but significantly faster than OPT=O0. There is another choice, OPT=Ofast, but it can cause various numerical instabilities, and is generally considered dangerous -although we have not observed any problems when running MiMeS.
It is important to note, once again, that the variables that correspond to template arguments must be known at compile time. Thus, if the compilation is done without the help of the various makefile files, the template arguments must be given, otherwise compilation will fail. 9 For example, the choice LONG=long, SOLVER=1, and METHOD=RODASPR2 will be used to compile the shared libraries (and C++ example in MiMeS/UserSpace/Cpp/Axion) with mimes::Axion<long double,1,RODASPR2<long double>>.
In order to fully understand the template arguments, the signatures of all classes and functions are given in Appendix B.

Compile-time input
Files MiMeS requires files that provide data for the relativistic degrees of freedom (RDOF) of the plasma, and the anharmonic factor. Although MiMeS is shipped with the standard model RDOF found in [36], and a few points for f (θ peak ) introduced in eq. (2.28), the user is free change them via the corresponding variables in MiMeS/Paths.mk. Moreover, there is a set of data for the QCD axion mass as calculated in ref. [25]. The variables pointing to these data files are cosmoDat, axMDat, and anFDat, for the RDOF, axion mass, and the anharmonic factor; respectively. The format of the files has to be the following: • The RDOF data must be given in three columns; T (in GeV), h eff , and g eff .
• The axion mass data must be given in two columns; T (in GeV), χ (in GeV 4 ). Here, χ is defined as in eq. (2.2). The user can provide a function instead of data for the axion mass, by leaving the axMDat variable empty.
• The data for the anharmonic factor must be give in two columns θ peak f (θ peak ); with increasing θ peak .
The paths to these files should be given at compile time. That is, once Paths.mk changes, "bash configure.sh" and then "make" must be ran. The user can change the content of the data files (without changing their paths), in order to use them without compiling MiMeS again. However, the user has to make sure that all the files are sorted so that the values of first column increase (with no duplicates or empty lines). In order to ensure this, it is advised to run "bash FormatFile.sh path-to-file" (in Appendix E there are some details on MiMeS/src/FormatFile.sh), in order to format the file (that should exist in "path-to-file") so that it complies with the requirements of MiMeS.
These paths are stored as strings in MiMeS/src/misc_dir/path.hpp at compile-time (they are defined as constexpr), and can be accessed once this header file is included. The corresponding variables are cosmo_PATH, chi_PATH, and anharmonic_PATH, for the path to data file of the plasma quantities, χ(T ), and f (θ peak ); respectively. Although, the axion mass data file may be omitted -since the axion mass is defined by the user, the variable chi_PATH is still useful if the axion mass is defined via a data file, as it is automatically converted to an absolute path. 10

Run-time input
The run-time user input is described in sec. 3. The user has to provide the parameters that describe the axion evolution, θ ini and f a .
Moreover, the maximum allowed value of u and the minimum value of T , allow the user to decide when integration stops even if the axion has not reached its adiabatic evolution. Ideally, umax= ∞ and TSTOP= 0, but MiMeS is designed to be as general as possible and there may be cases where one needs to stop the integration abruptly. 11 Furthermore, ratio_ini allows the user to choose a desired point at which the interpolation of the data in inputFile begins. This can save valuable computation time as well as memory, as only the necessary data are stored and searched. Generally, for ratio_ini> 10 3 , the relic abundance becomes independent from ratio_ini, but one has to choose it carefully, in order to find a balance between accuracy and computation time.
Finally, the convergence conditionsi.e. N_convergence_max and convergence_lim -allow the user to decide when the adiabatic evolution of the axion begins. Generally, the relic abundance does not have a strong dependence on these parameters as long as N_convergence_max> 3 and convergence_lim< 10 −1 (i.e. the adiabatic invariant does not vary more that 10% for three consecutive peaks of the oscillation). However, we should note that greedy choices (e.g. N_convergence_max= 100 and convergence_lim= 10 −5 ) are dangerous, as θ tends to oscillate rapidly and destabilize the differential equation solver. Therefore, these parameters should be chosen carefully, in order to ensure that integration stops when axion has reached its adiabatic evolution, without destabilising the EOM.

Complete Examples
Although one can modify the examples provided in MiMeS/UserSpace, in this section we show a complete example in both C++ and python. The underlying cosmology is assumed to be an EMD scenario, 12 . with the evolution of the energy densities of the plasma and the matter field (Φ) shown in Fig. (3). The various regime changes are indicated as T E 1,2 -corresponding to the first and second time where ρ Φ = ρ R , and T D 1,2 -which correspond to the start and end of entropy injection (for more precise definition see [38]). Regarding the axion, this cosmological scenario alters its evolution, since both the Hubble parameter and the entropy of the plasma change significantly. This results in a shifted oscillation temperature, and the dilution of the axion energy density.
Moreover, we will use the QCD axion mass of ref. [25]. For the axion mass beyond the minimum and maximum temperatures, T min, max , in the corresponding data file, we will takẽ (4.1)

complete example in C++
In order to write a C++ program that uses MiMeS in order to solve the EOM 2.18, we must include the header files src/AxionMass/AxionMss.hpp and src/Axion/AxionSolve.hpp. In the example at hand, the main function should contain the definition of the axion mass 1 // use chi_PATH to interpolate the axion mass. / * This is the axion mass squared beyond the interpolation limits 5 for the current data. If you don't specify them, the axion mass 6 is taken to be constant beyond these limits * / It should be noted that, since we have used chi_PATH, we need to include the path.hpp header file from MiMeS/src/misc_dir. Moreover, we note that since the second and third arguments are 0 and mimes::Cosmo<LD>::mP (the Planck mass; see B), interpolation of the axion mass will be between the minimum and maximum values that appear in the data file.
Following, the axion EOM can be defined and solved as 1e-1, 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7)); 7 8 ax.solveAxion(); In order to print anything, however, we need to write a few more lines of code. For example, the relic abundance is printed by adding std::cout<<ax.relic<<"\n"; after ax.solveAxion();. It should be noted that, in order to make sure that the code is compiled successfully, one should add at the top of the file any other header they need. In particular, in this case, we need to add #include<iostream>, in order to be able to use std::cout.
Parameter choice It should be noted that the string that is assigned to the variable inputFile, is a path of a file that exists in MiMeS/UserSpace/InputExamples. Also, rootDir is a constant static (char[]) variable that is defined in MiMeS/src/path.hpp, and it is automatically generated when bash configure.sh is executed. Moreover, the other parameters used in the constructor are: • relative_tolerance= 10 −11 • beta= 0.9 • fac_max= 1.2 • fac_min= 0.8 • maximum_No_steps= 10 7 Template arguments Furthermore, The first template parameter is long double, which means that all the numeric types (apart from integers like maximum_No_steps) have "long double" precision. which is useful when considering low tolerances as in this example. The second and third template parameters are responsible for the Runge-Kutta method we use. That is, in this case, the second template parameter -1 -means that we use a Rosenbrock method, with the method being RODASPR2.
Notice that the method also needs a template parameter that is used to declare all member variables of the RODASPR2 class as long double.
Compilation Before we compile, we have to make sure that bash configure.sh has been executed. Assuming that we name the file that contains the code is axionExample.cpp, and it is located in MiMeS/UserSpace, an executable can be produced as Both of these commands should create an executable that solves the axion EOM. That is, assuming we have a terminal open in MiMeS/UserSpace, we can run ./axion, and get the results we chose. It should be noted that all paths that MiMeS uses by default are written as absolute paths in MiMeS/src/misc_dir/paths.hpp when we run bash configure.sh. Therefore, if the inputFile variable is also an absolute path, the executable can be copied and used in any other directory of the same system. However, it should be preferred that executables are kept under the MiMeS directory, in order to be able to compile them with different data file paths if needed.
Alternative axion mass definition For this particular example, we could have used the approximate definition of the axion mass

complete example in python
In order to be able use the AxionMass and Axion classes in python, we need to import the corresponding modules from MiMeS/src/interfacePy. That is, assuming that the script from which we intend to import MiMeS/src/interfacePy/Axion/Axion.py and MiMeS/src/interfacePy/Axion/AxionMass.py is in MiMeS/UserSpace, on top of the script, we need to write # This is the axion mass squared beyond the interpolation limits for the current data.

5
# If you don't specify them, the axion mass is taken to be constant beyond these limits Again, the parameters passed to the constructor are the same as in the C++ example. The axion EOM, then, is solved using 1 ax.solveAxion () In contrast to the C++ usage of this function, this only stores the θ ini , f a , θ osc , T osc , and Ωh 2 in the variables ax.theta_i, ax.fa, ax.theta_osc, ax.T_osc, and ax.relic; respectively. Therefore, we can print the relic abundance by calling print(ax.relic). In order to get the integration points (i.e. the evolution of the angle and the other quantities), the quantities at the peaks, and the local integration errors, we need to call The documentation of any python function can be read directly inside the script using ? as a prefix. For example, in order to see what the functionality and usage of the getPeaks function, we can call "?ax.getPeaks", and its documentation will be printed. As already mentioned, it is important to always delete any instance of the AxionMass and Axion classes once they are not needed. In this case this is done by calling The entire code As in the C++ example, this example consists of only a few lines of code. The script we described here is  One can find a complete example, including the option to create several plots, in the script MiMeS/UserSpace/Python/Axion.py. Also, the same example can be used interactively, in jupyter notebook environment [39], that can be found in MiMeS/UserSpace/JupyterNotebooks/Axion.ipynb. One can read the comments, and change all different parameters, in order to examine how the results are affected.
Alternative axion mass definition In order to use the approximate axion mass as defined in eq.  Figure 4: The evolution of the axion angle, θ, with the temperature in early matter dominated (a) and radiation dominated (b) cases. For the EMD case, the θ ini is chosen so that Ωh 2 = 0.12.
In Fig. (4a) we show the evolution of θ for temperatures T ∈ [0.025, 0.15] GeV, where the vertical line indicates T osc , while the horizontal one the corresponding value of θ, θ osc . Also, the blue curve connects the peaks of the oscillation. For comparison, in Fig. (4b) we also show the evolution of the angle in a radiation dominated Universe with constant entropy (i.e. standard cosmological scenario). From these figures, we can see that the effect of an early matter domination reduces the amplitude of oscillation -due to the injection of entropy -as well the oscillation temperature -due to the increase of the Hubble parameter. Furthermore, the angle at T = T osc in the EMD scenario is much smaller that the corresponding value in the standard cosmological case, since the entropy injection causes the scale factor to be larger compared to the scale factor at the same temperature. Moreover, in Fig. (5), we show the relative local error of integration as well as a histogram of the number of integration steps for 0.025 GeV < T < 0.15 GeV. The local relative errors, defined in Appendix A, are shown in Fig. (5a). The black and red lines correspond to δθ/θ and δζ/ζ, respectively. This figure indicates that the local errors are relatively well behaved for T > T osc , and they only start to oscillate violently once the oscillations start. However, the adaptation of the integration step seems to work, as the errors are kept below ∼ 10 −7 (the initial relative error of ζ is large because ζ ≈ 0 before the T osc ). In order to examine how the step-size is adapted to the difficulty of the problem, we also show a histogram which shows how many integration steps are taken for fixed temperature intervals. 13 In this figure, we see that the number of integration steps increases rapidly for temperatures below the oscillation temperature. This is expected, since integration becomes more difficult as the frequency of the oscillation increases -m a /H increases rapidly for T < T osc . That is, the local integration error tends to increase. Thus, in order to reduce the local error, the embedded RK method we employ, reduces the step-size. This means that for T T osc the number of integration steps increase drastically. This is the general picture of how adaptation happens. However, one has to experiment with all the available parameters, in order to solve the axion EOM as accurately and fast as possible.

Acknowledgements
The author acknowledges support by the Lancaster-Manchester-Sheffield Consortium for Fundamental Physics, under STFC research grant ST/T001038/1.

Summary
We have introduced MiMeS; a header-only library written C++ that is used to compute the axion (or ALP) relic abundance, by solving the corresponding EOM, in a user defined underlying cosmology. MiMeS makes only a few assumptions, which allows the user to explore a wide range ALP and cosmological scenarios.
In this manuscript, we have provided a detailed explanation on how to use MiMeS, by showing examples in both C++ and python (paragraphs (3.1) and (4.3)). We have described the user input that is expected, and the various options available (sections (3) and (4)). Moreover, in the Appendix, we provide a detailed review of all the internal components that comprise MiMeS. We briefly discuss the Runge-Kutta methods that MiMeS uses, and show how the user can implement their own. We explain the functionality of all the classes, modules, and utilities. Also, we provide a detailed input and option guide.
In the future, MiMeS will be extended in several ways. First, we should implement new functionality that will allow the user to automatically compare against various experimental data (although, there is already an available module [40] that can be used for this goal). This is going to be helpful, as the user will only need to use one program (or script) to compute what is needed. We also aim to supplement MiMeS with the option to produce ALPs via interactions of the plasma (e.g. freeze-out/in), which may be useful in certain cases. Moreover, a later version of MiMeS may allow the user to define different initial condition for ζ as well, since there are cases where this is needed (e.g. [27]). Finally, MiMeS will continue to improve by correcting mistakes, or implementing suggestions by the community.

Appendix A Basics of embedded Runge-Kutta Mehtods
Runge-Kutta (RK) methods are employed in order to solve an ordinary differential equation (ODE), or a system of ODEs of first order. 14 Although there are some very insightful sources in the literature (e.g. [41,42,43]) we give a brief overview of them in order to help the user to make appropriate decisions when using MiMeS.
The general form of a system of first order of ODEs is with given initial condition y(0). Also, the components of y denote the unknown functions. Note that we can always shift t to start at 0, which simplifies the notation. In order to solve the system of eq. (A.1), an RK method uses an iteration of the form with n denoting the iteration number, h the "step-size" that is used to progress t; t n+1 = t n + h. Moreover, s, b i and k i define the corresponding RK method. For example, the classic Euler method is an RK method with s = 1, b = 1, and k 1 = f ( y n , t n ). Methods with k that depends on previous step (i.e. k = k [y n , t n ]), are called explicit, while the ones that try to also predict next step (i.e. k = k [y n , t n ; y n+1 , t n+1 ]) are called implicit. 15

A.1 Embedded RK methods
A large category of RK methods are the so-called embedded RK methods. These methods make two estimates for the same step simultaneously -without evaluating k many times within the same iteration. Therefore, together with the iteration of eq. (A.2), a second estimate is given by with b * i is an extra parameter that characterise the "embedded" method of different order (typically, one order higher that the estimate A.2). The local (for the step t n+1 ) error divided by the scale of the solution, then, estimated as 14 Boundary value problems, and higher order differential equations are expressed as first order ODEs, and then solved. Similarly to eq. (2.18). 15 Generally, by substituting implicit methods yn+1 as in eq. (A.2), we end up with a system of equations that need to be solved in order to compute k.

MiMeS
Runge-Kutta absolute_tolerance Atol relative_tolerance Rtol b β fac_min f min fac_max f max Table 1: Correspondence between the user MiMeS user input and the RK parameters described in the text.
where n the iteration number, N the number of ODEs, d the component of y, and Λ d defined as with Atol and Rtol the absolute and relative tolerances that characterise the desirable accuracy we want to achieve; user defined values, typically Atol=Rtol 1. With these definitions, the desirable error is reached when ∆ 1.
Step-control The definition A.4, allows us to adjust the step-size h in such way that ∆ ≈ 1. That is, we take trial steps, and h is adapted until ∆ ≈ 1. A simple adaptive strategy adjusts step-size, using with p the order of the RK method (p + 1 is the order of the embedded one), β a bias factor of the adaptive strategy (typically β is close but below 1), used to adjust the tendency of h to be somewhat smaller than what the step-control predicts. Also, f min and f max are the minimum and maximum allowed factors, respectively, that can multiply β h, used in order to avoid large fluctuations that can destabilise the process. All these parameters are chosen by the user, in order to make the step-control process as aggressive or safe as needed.
Correspondence between MiMeS parameters and RK ones The various parameters that MiMeS are described in section 4.2. The correspondence between them and the RK parameters is given in table 1.

A.2 Explicit embedded RK methods
Explicit methods use only the information of the previous step in order to compute k i from with a ij , c i , together with b i and b * i , consist the so-called Butcher tableau of the corresponding method. For explicit methods this is usually presented as c s a s1 a s2 a s3 . . . a ss−1 0 It should be noted that c i = s j=1 a ij .

A.3 Rosenbrock methods
Explicit RK methods, encounter instabilities when a system is "stiff" 16 ; e.g. when it oscillates rapidly, or has different elements at different scales. These problems are somewhat resolved by trying to predict the next step inside k; i.e. in implicit methods. However, then one has to solve a non-linear set of equations in order to compute y n+1 , which is generally a slow process, as the Newton method (or some variation) needs to be applied. However, there exist another way, a compromise between explicit and implicit methods. Linearly implicit RK methods, usually called Rosenbrock methods -with popular improvements as Rosenbrock-Wanner methods, introduce parameters in the diagonal of the Butcher tableau A.8 and linearise the system of non-linear equations (for details, see [42]). In these methods, k is determined by which is written in such a way that everything is evaluated at t = t n . In eq. (A.9),Ĵ = ∂ f ∂ y the Jacobian of the system of ODEs,Î the unit matrix with dimension equal to the number of ODEs. Moreover, γ and γ ij are parameters that characterise the method (along with a ij , c i , b i , and b * i ).
Implementing a new Butcher tableau in NaBBODES As already mentioned, MiMeS uses NaBBODES, which supports the implementation of new Butcher tableaux. This is done by adding a new class (or struct) inside the header file METHOD.hpp that can be found in MiMeS/src/NaBBODES/RKF for the explicit RK and MiMeS/src/NaBBODES/Rosenbrock for the Rosenbrock embedded methods. All the new parameters must be public, constexpr static variables, of a type that is the template parameter of the method. For example, the Heun-Euler method can be implemented by adding the following code in MiMeS/src/NaBBODES/RKF/METHOD.hpp In order to implement the ROS3w [34] method, one can add the following code in the header file MiMeS/src/NaBBODES/Rosenbrock/METHOD.hpp  How to compile MiMeS in order to use the newly implemented method Once a new Butcher tableau is implemented, the mimes::Axion class can use it. This class, just needs the name of method assigned to the corresponding template argument; Method. A convenient way to do this, is to define a macro using the -D flag of the compiler, and use macro as the corresponding template parameter of the mimes::Axion class. Alternatively, if one uses the makefile files, a method can be chosen by adding it in the corresponding Definitions.mk file; e.g. as METHOD=ROS3w for the ROW3 method. 17

B C++ classes
MiMeS is designed as an object-oriented header-only library. That is, all the basic components of the library are defined as classes inside header files. All the classes relevant to the use of MiMeS are under the namespace mimes.

B.1 Cosmo class
The mimes::Cosmo<LD> class is responsible for interpolation of the various quantities of the plasma. Its header file is MiMeS.src/Cosmo/Cosmo.hpp, and needs to be included in order to use this class. The template parameter LD is the numeric type that will be used, e.g. double. The constructor of this class is 1 template<class LD> 2 mimes::Cosmo<LD>(std::string cosmo_PATH, LD minT=0, LD maxT=mimes::Cosmo<LD>::mP) The argument cosmo_PATH is the path of the data file that contains T (in GeV), h eff , g eff , with increasing T . The parameters minT and maxT are minimum and maximum interpolation temperatures. These temperatures are just limits, and the action interpolation is done between the closest temperatures in the data file. Moreover, beyond the interpolation temperatures, both h eff and g eff are assumed to be constants.
Interpolation of the RDOF, allows us to define various quantities related to the plasma; e.g. the entropy density is defined as s = 2π 2 45 h eff T 3 . These quantities are given as the member functions: • template<class LD> LD mimes::Cosmo<LD>::heff(LD T): h eff as a function of T .

B.2 AnharmonicFactor class
The class mimes::AnharmonicFactor<LD> is responsible interpolating the anharmonic factor as defined in eq. (2.28). The corresponding header file is MiMeS/src/AnharmonicFactor/AnharmonicFactor.hpp. The constructor of this class is 1 template<class LD> 2 mimes::AnharmonicFactor<LD>(std::string anharmonic_PATH) Again, the template argument LD is a numeric type, and the anharmonic_PATH string is the path of the data file with data for θ peak (which should be in increasing order) and f (θ peak ). The member function that MiMeS uses is the overloaded call operator 1 template<class LD> LD mimes::AnharmonicFactor<LD>::operator()(LD theta_peak) This function returns the value of the anharmonic factor at θ peak =theta_peak. Although, there is no need to call this function beyond the interpolation limits (as long as the data file contains 0 ≤ θ peak ≤ π), it is important to note that the anharmonic factor is taken to be constant beyond these limits.

B.3 AxionMass class
The mimes::AxionMass<LD> class is responsible for the definition of the axion mass. The header file of this class is MiMeS/src/AxionMass/AxionMass.hpp. Its usage and member functions are described in the examples given in sections (3.1) and (4.3). However, it would be helpful to outline them here. The class has two constructors. The first one is 1 template<class LD> 2 mimes::AxionMass<LD>(std::string chi_PATH, LD minT=0, LD maxT=mimes::Cosmo::mP) The first argument, chi_PATH, is the path to a data file that contains two columns; T (in GeV) and χ (in GeV 4 ), with increasing T . The arguments minT and maxT are the interpolation limits. These limits are used in order to stop the interpolation in the closest temperatures that exist in the data file. That is the actual interpolation limits are T min ≥minT and T max ≤maxT. Beyond these limits, by default, the axion mass is assumed to be constant. However, this can be changed by using the member functions Here, ma2_MIN and ma2_MAX are functors that define the axion mass squared beyond the interpolation limits. In order to ensure that the axion mass is continuous, usually we need T min , T max , χ(T min ), and χ(T max ). These values can be obtained using the member functions • template<class LD> LD mimes::AxionMass<LD>::getTMin(): This function returns the minimum interpolation temperature, T min .
An alternative way to define the axion mass is via the constructor Here, the only argument is the axion mass squared,m a , defined as a callable object.
Once an instance of the class is defined, we can getm a 2 using the member function We should note that ma2 is a public std::function<LD(LD,LD)> member variable. Therefore, it can be assigned using the assignment operator. However, in order to change its definition, we can also use the following member function:

B.4 AxionEOM class
The mimes::AxionEOM<LD> class is not useful for the user. However, it is responsible for the interpolation of the underlying cosmology, and the definition of the axion EOM 2.18, which is passed to the ODE solver of NaBBODES.
The constructor of the class is It should be noted that the highest interpolation temperature is determined by ratio_ini while the lower interpolation temperature is the one given in the data file inputFile. Beyond these limits, all functions are assumed to be constant. Therefore, one should be careful, and choose an appropriate ratio_ini, and provide a lower temperature at which any entropy injection has stopped and the axion has reached its adiabatic evolution. Finally, the actual EOM is given an overloaded call operator

B.5 Axion class
The mimes::Axion<LD,Solver,Method> class is the class that combines all the others, and actually solves the axion EOM 2.18. Its header file is MiMeS/src/Axion/AxionSolve.hpp, and its constructor is We should note that running this function all variables are cleared. So we lose all information about the last time axionSolve() ran.
In case the mass of the axion is changed, we also need to remake the interpolation (i.e. run mimes::AxionEOM::makeInt()). This is done using This variable can be used without an instance of the mimes::Axion<LD,Solver,Method> class.

C MiMeS python interface
The various python modules, classes, and functions are designed to work exactly in the same way as the ones in C++. All the modules are located in src/interfacePY, so it is helpful to add the MiMeS/src path to the system path at the top of every script that uses MiMeS. This is done by adding 1 from sys import path as sysPath The available models are Cosmo, AxionMass, and Axion, each defines a class with the same name.

C.1 Cosmo class
The The argument cosmo_PATH is the path (a string) of a data file that contains T (in GeV), h eff , g eff , with accenting T . The second and third arguments, minT and maxT, are minimum and maximum interpolation temperatures, with the interpolation being between the closest temperatures in the data file. Moreover, beyond these limits, both h eff and g eff are assumed to be constants. It is important to note that the class creates a void pointer that gets recasted to mimes::Cosmo<LD> in order to call the various member functions. This means that once an instance of Cosmo is no longer needed, it must be deleted, in order to free the memory that it occupies. An instance, say cosmo, is deleted using The several cosmological quantities are given as members variables: • Cosmo.T0: CMB temperature today [44] in GeV.
Note that these values can be directly imported from the module, without declaring an instance of the class, as

C.2 AxionMass class
The AxionMass class is defined in the module with the same name that can be found in the directory MiMeS/src/interfacePy/AxionMass. This class is responsible for the definition of the axion mass. This module loads the corresponding shared library from MiMeS/lib/libma.so, which is created by compiling MiMeS/src/AxionMass/AxionMass.cpp using "make lib/libma.so". Its usage is described in the examples given in sections (3.1) and (4.3). Moreover, this class is used in the same way as mimes::AxionMass<LD>. However, we should append in this section the definition of its member functions in python.
The class is imported using 1 from interfacePy.AxionMass import AxionMass The constructor is 1

AxionMass(*args)
and can be used in two different ways. First, one can pass three arguments, i.e. 1

AxionMass(chi_PAT, minT, maxT)
The first argument is the path to a data file that contains two columns; T (in GeV) and chi (in GeV 4 ), with increasing T . The arguments minT and maxT are the interpolation limits. These limits are used in order to stop the interpolation in the closest temperatures in chi_PATH. That is the actual interpolation limits are T min ≥minT and T max ≤maxT. Beyond these limits, by default, the axion mass is assumed to be constant. However, this can be changed by using the member functions Here, ma2_MIN(T,fa) and ma2_MAX(T,fa), are functions (not any callable object), should take as arguments T and fa, and return the axion mass squared beyond the interpolation limits. In order to ensure that the axion mass is continuous, usually we need T min , T max , χ(T min ), and χ(T max ). These values can be obtained using the member functions An alternative way to define the axion mass is via the constructor 1

AxionMass(ma2)
Here, ma2(T,fa) is a function (not any callable object) that takes T (in GeV) and f a , and returns m 2 a (T ) (in GeV). As in the other python classes, once the instances of this class are no longer needed, they must be deleted using the destructor, del.
The member function that returnsm a 2 is 1

AxionMassma2(T,fa)
We should note that another ma2 can be changed using the following member function: Again, ma2(T,fa) is a function (not any callable object) that takes T (in GeV) and f a , and returns m 2 a (T ) (in GeV).

C.3 Axion class
The class Axion, solves the axion EOM 2.18. This class is defined in MiMeS/interfacePy/Axion/Axion.py, and imports the corresponding shared library. This library is compiled by running make lib/Axion_py.so, and its source file is MiMeS/src/Axion/Axion-py.cpp. As in the previous classes, its usage is similar to the C++ version.

Axion.solveAxion()
Once this function is finished, the following member funcions are available • Axion.T_osc: the oscillation temperature, T osc , in GeV.
• Axion.a_osc: a a ini at the oscillation temperature.
• Axion.relic The relic abundance of the axion.
The evolution of a/a ini , T, θ, ζ, ρ a at the integration steps, is not automatically accessible to user, but they can be made so using 1 Axion.getPoints() Then, the following member variables are filled • Axion.a: The scale factor over its initial value, a a ini .
• Axion.T: The temperature in GeV.
Moreover, the function fills the (numpy) arrays Axion.a_peak, Axion.T_peak, Axion.theta_peak, Axion.zeta_peak, Axion.rho_axion_peak, and Axion.adiabatic_invariant with the quantities a/a ini , T, θ, ζ, ρ a , J, at the peaks of the oscillation. These points are computed using linear interpolation between two integration points with a change in the sign of ζ. The local integration errors for θ and ζ are stored in Axion.dtheta and Axion.dzeta, after the following function is run 1 Axion.getErrors() Another initial condition, θ ini , can be used without declaring a new instance using 1

Axion.setTheta_i(theta_i)
We should note that running this function all variables are cleared.
As in the previous python classes, once an instance of the Axion class is no longer needed, it needs to be deleted, by calling the destructor, del.
Important difference between the C++ version Since the axion mass is passed by value in the constructor, a change of the AxionMass instance has no effect on the Axion instance that uses it. Therefore, if the definition of the axion mass changes, one has to declare a new instance of the Axion class. The new instance can be named using the name of the previous one, if the latter is deleted by running its destructor.

D Other modules
There are several modules available, that may help the user scan to obtain approximate results using the WKB approximation, sca the parameter space, and plot some results. These modules are not integral to MiMeS, and the user is free to ignore them.

D.1 WKB module
The WKB module can be used to calculate the axion relic abundance using the WKB approximation discussed in section 2. The module can be imported using 1 from interfacePy import WKB It contains the definition of a function that returns the relic abundance using the WKB approximation, which is 1 WKB.relic(Tosc, theta_osc ,ma2, gamma=1, 2 cosmo=Cosmo(_PATH_+r'/src/data/eos2020.dat',0,Cosmo.mP)) Here, Tosc is the oscillation temperate, ma2(T,fa) is m a 2 (T ) as a function that takes T and f a as arguments, cosmo an instance of the Cosmo class, and gamma the entropy injection (as defined in eq. (2.13)).
Moreover, there is a function that helps to determine T osc and γ, which is 1 WKB.getPoints(fa, ma2, inputFile, cosmo=Cosmo(_PATH_+r'/src/data/eos2020.dat',0,Cosmo.mP)) The arguments fa, inputFile,f a , and cosmo the path to a file that describes the cosmology (described  in table 2), f a (in GeV ), an instance of the Cosmo class, respectively. This function returns γ (the entropy injection between Tosc and today) and T osc .

D.2 The ScanScript module
The python interface of MiMeS has two simple classes that help to scan over θ ini and f a , in parallel. These two classes can be imported from the module ScanScript. Tis module is based on the bash script MiMeS/src/util/parallel_scan.sh, which automatically performs scans in parallel. This script is used as 1 bash src/util/parallel_scan.sh executable cpus inputFile Here, executable is the path to an executable, cpus the number of instances of the executable to launch simultaneously, and inputFile a file that contains the arguments the executable expects. This file should contain arguments for the executable in each line. This script, then, separates all in arguments in the inputFile in batches of size cpus, and runs one batch at a time.

D.2.1 The Scan class
The Scan class writes a time-coded file (so it would be unique) with columns that correspond to θ ini , f a (GeV), θ osc , T osc (in GeV), and Ωh 2 , for every combination of θ ini and f a that are passed as input.
This class is imported using 1 from interfacePy.ScanScript import Scan • break_after, break_time: take a break after break_after seconds for break_time seconds.
• break_command (optional): before it takes a break, run this system command (this may be a script to send the results via e-mail, or back them up).
The scan runs using the method 1

Scan.run()
For every value of f a , this method calls another one Scan.run_batch(), which in turn calls the bash script MiMeS/src/util/parallel_scan.sh with arguments PathToCppExecutable, cpus, and a file that contains all the inputs for PathToCppExecutable for all values in table_theta.
As the scan runs, it prints the number of batches that have been evaluated, the mean time it takes to evaluate one batch, and an estimate of the remaining time. It should be noted that if the scan exits before it finishes, the next run will continue from the point at which it was stopped even if the inputs have changed. In order to start from the beginning, the user must delete the file count._mimes_.

D.2.2 The ScanObs class
The ScanObs class scans for different values of f a (given as a table) and finds the value of θ ini closer to the observed DM relic abundance. The result file is a time-coded file with columns correspond that to θ ini , f a (in GeV), θ osc , T osc (in GeV), and Ωh 2 . It is imported using Here, len_theta is the number of different values of θ to be used in the search of a value closer to the central value of the observed DM relic abunance. The arguments relic_obs, relic_err_up, relic_err_low are the central value of Ωh 2 , and its upper and lower error; respectively. All the other arguments have already been defined previously.
The scan runs using the method 1

ScanObs.run()
For every value of f a , this method first calculates the relic abundance for θ ini 1, and finds the value of θ ini that would result in Ωh 2 =relic_obs; θ ini (approx) . This is easy to do, since for θ ini 1, the EOM becomes independent of θ ini .
Similarly to the Scan class, if the scan exits before it finishes, the next run will continue from the point at which it was stopped. In order to start from the beginning, delete the file count._mimes_.

D.3 FT module
The FT class, defined in FT module that can be found in src/interfacePy/FT, is used to format the ticks when plotting using matplotlib [45]. The class can be imported using The various arguments are • _M_xticks,_M_yticks: A list for the major ticks in the x and y axes, respectively.
• _M_xticks,_M_yticks: A list for the major ticks in the x and y axes, respectively, for which no number is printed.
• _m_xticks,_m_yticks: A list for the minor ticks in the x and y axes, respectively.
• xscale, yscale: The scale of the x and y axes, respectively. The available values are "linear", "log", and "symlog".
Instances of FT should be defined after subplots have been defined, because the code seems clearer. Once this is done, we can format the ticks of a subplot, sub, using

Example of FT
Consider as an example the plot of f (x) = e − 1 2 x 2 for −5 ≤ x ≤ 5. In order to do this, we need to run

E Utilities
There are various utilities, that can be found MiMeS/src/util, which can be used to make the use of MiMeS easier. In this section we discuss how they work.

E.1 FormatFile.sh
This is a bash script that formats a data file so that it is compatible with the interpolation assumptions. It takes a path to a file as an argument, and it returns the same file, sorted (in ascenting order) with respect to the first column, with all duplicate or empty lines as well as the last "new line" (i.e. \n) removed. The script is called as Here path_to_dat_file is a path to the data file we would like to replace with the formatted one. Notice that FormatFile.sh must be run in order to ensure that the cosmological input needed by the mimes::Axion<LD,SOLVER,METHOD<LD>> class, is acceptable. Moreover, note that FormatFile.sh runs every time configure.sh is called, in order to format the data files for the RDOFs, χ(T ), and the anharmonic factor, in order to ensure that they comply with what MiMeS expects.

E.2 timeit.sh
This script takes two arguments. The first should be a path to an executable, the second a file that contains arguments to passed to the executable. The script, then, runs the executable, and prints in stderr the time it took (is seconds) to run it. It is important to note that each argument of the executable should be in a different line of the file. 18 The script is called as

E.3 Timer C++ class
The class mimes::util::Timer, defined in MiMeS/src/util/timeit.hpp, can be used in order to time processes in C++. An instance of this class prints its life-timei.e. time from the moment it was created until it went out of scope -in stderr. Instances of this class are intended to be used inside a scope in order to count the time it took to run a block of code. For example the time it takes to compute 10 4 i=0 i can be determined as 1 #include "src/util/timeit.hpp" For example , if we have an executable that takes two arguments, let's assume the first is 3 and the second is "foo", the file should read:

foo
Timer _timer_; F Quick guide to the user input We present tables (2), (3), (4), and (5), with the various available run-time inputs, required files, and template arguments. In table 6 we show the available compile-time options, that can be used when compiling using the various makefile files.
User run-time input for solving the EOM 2.18.

theta_i
Initial angle. fa The PQ scale.
umax Once u = log a/a i >umax, the integration stops. Typical value: ∼ 500.

N_convergence_max convergence_lim
Integration stops when the relative difference between two consecutive peaks is less than convergence_lim for N_convergence_max consecutive peaks.
inputFile Relative (or absolute) path to a file that describes the cosmology. The columns should be: u T [GeV] log H, with acceding u. Entropy injection should have stopped before the lowest temperature of given in inputFile. axionMass Instance of mimes::AxionMass<LD> class. In C++ this instance is passed as a pointer to the constructor of the mimes::Axion<LD,Solver,Method> class, while in python it is simply passed as a variable.
absolute_tolerance Absolute tolerance of the RK solver (see also beta Aggressiveness of the adaptation strategy (see also table 1). Default value: 0.9.

fac_max, fac_min
The step-size does not change more than fac_max and less than fac_min within a trial step (see also  Input required for the definition of the axion mass via a data file. chi_Path Relative or absolute path to data file with T (in GeV), χ(T ) (in GeV 4 ).
minT and maxT Interpolation limits. These are used in order to stop the interpolation at the closest temperatures that exist in the data file. This means that the actual interpolation limits are T min ≥minT and T max ≤maxT. Beyond these limits hat axion mass is assumed to be constant. However, this can be changed using mimes::AxionMass<LD>::set_ma2_MIN(std::function<LD(LD,LD)> ma2_MIN) and mimes::AxionMass<LD>::set_ma2_MAX(std::function<LD(LD,LD)> ma2_MAX).
Input required for the definition of the axion mass via a function.

ma2
A function (or a callable object) with signature LD ma2(LD T, LD fa) that returnsm a 2 (T ). Table 3: Arguments of the contructors of the mimes::AxionMass<LD> class.
Required data files, with corresponding variables in MiMeS/Paths.mk.
cosmoDat Relative path to data file with T (in GeV), h eff , g eff . If the path changes one must run bash configure.sh and make.
If the path changes one must run bash configure.sh and make. This variable can be omitted, since the user can define an AxionMass instance using any path.
anFDat Relative path to data file with θ peak , f (θ peak ). If the path changes one must run bash configure.sh and make.

LD
This template argument appears in all classes of MiMeS. The preferred choice is long double. However, in many cases double can be used. The user should be careful, as the later can lead to an inaccurate result; especially for low tolerances, and small values of θ.
Solver This is the second template argument of the mimes::Axion<LD,Solver,Method> class. The available choices are Solver=1 for Rosenbrock method, and Solver=2 for explicit RK method.
Method The third template argument of the mimes::Axion<LD,Solver,Method> class. Its value depends on the choice of Solver; For Solver=1, Method can be either RODASPR2<LD> (fourth order) or ROS34PW2<LD> (third order). For Solver=2, Method can only be DormandPrince<LD> (seventh order). Notice that the definitions of the various method classes, also need a template argument,LD, that must be the same as the first template argument of the mimes::Axion<LD,Solver,Method> class. If one defines their own Butcher table (see Appendix A), then they would have to follow their definitions and assumptions. OPT Available options are OPT=O1, O2, O3 (be default). This variable defines the optimization level of the compiler. The variable can be changed in all Definitions.mk files. In the root directory of MiMeS, the optimization level applies to the python modules (i.e. the shared libraries), while in the subdirectories of MiMeS/UserSpace/Cpp it only applies to example inside them. Table 6: User compile-time input and options. These are available in the various Definitions.mk files, which are used when compiling using make.