#  MiMeS: Misalignment Mechanism Solver

Karamitros Dimitrios

###### Abstract

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.

Program summary:

Program title: MiMeS.

Developer’s respository link: https://github.com/dkaramit/MiMeS.

Programming language: C++ and python.

Licensing provisions: MIT license.

Nature of problem: Solving numerically the axion (or axion-like-particle)
equation of motion, in order to determine the corresponding relic abundance.
The library is designed to be quite general, and can be used to obtain the
relic abundance in various cosmological scenarios, and axion-like-particle
models.

Solution method: Embedded Runge-Kutta for the numerical integration of the
equation of motion. The user may choose between explicit and Rosenborck
methods, or implement their own Butcher tableau. For the various
interpolations, the library uses cubic splines.

Restrictions: The the derivative of the axion-angle initially is assumed
vanish. This is hard-coded in the library, and there is no easy way for the
user to change it. Furthermore, any additional contribution from decays or
annihilation of plasma particles to the axion (or ALP) energy density is
assumed to be subdominant.

###### Contents

  1. 1 Introduction
  2. 2 Physics background
    1. 2.1 The WKB approximation
    2. 2.2 Notation
    3. 2.3 Adiabatic invariant and the anharmonic factor
  3. 3 MiMeS usage
    1. 3.1 First steps
      1. 3.1.1 Using MiMeS in C++
      2. 3.1.2 Using MiMeS in python
  4. 4 Assumptions and user input
    1. 4.1 Restrictions
    2. 4.2 Options at Compile-time
    3. 4.3 User input
      1. 4.3.1 Compile-time input
      2. 4.3.2 Run-time input
    4. 4.4 Complete Examples
      1. 4.4.1 complete example in C++
      2. 4.4.2 complete example in python
      3. 4.4.3 Results
  5. 5 Acknowledgements
  6. 6 Summary
  7. A Basics of embedded Runge-Kutta Mehtods
    1. A.1 Embedded RK methods
    2. A.2 Explicit embedded RK methods
    3. A.3 Rosenbrock methods
  8. B C++ classes
    1. B.1 Cosmo class
    2. B.2 AnharmonicFactor class
    3. B.3 AxionMass class
    4. B.4 AxionEOM class
    5. B.5 Axion class
  9. C MiMeS python interface
    1. C.1 Cosmo class
    2. C.2 AxionMass class
    3. C.3 Axion class
  10. D Other modules
    1. D.1 WKB module
    2. D.2 The ScanScript module
      1. D.2.1 The Scan class
      2. D.2.2 The ScanObs class
    3. D.3 FT module
  11. E Utilities
    1. E.1 FormatFile.sh
    2. E.2 timeit.sh
    3. E.3 Timer C++ class
    4. E.4 linspace function
    5. E.5 logspace function
    6. E.6 map function
  12. F Quick guide to the user input

##  1 Introduction

The axion is a hypothetical particle that was originally introduced in order
to solve the strong CP-problem 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
$100{\rm GeV}$), 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 non-standard contribution to the energy density is absent before
Big Bang Nucleosynthesis [12, 13] becomes active; for temperatures [14, 15,
16, 17] $T\gtrsim\mathcal{O}(10)~{}{\rm 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,
111MiMeS 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. 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, 222Especially 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.
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.

##  2 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, $\theta$, as $A=\theta\ f_{a}{}$, with $f_{a}{}$ the scale at which the
PQ symmetry breaks. 333If in a model under study, there is no $f_{a}{}$, this
parameter is still expected by MiMeS, but the user can set $f_{a}{}=1$. The
axion angle follows the EOM

$\left(\dfrac{d^{2}}{dt^{2}}+3H(t)\
\dfrac{d}{dt}\right)\theta(t)+\tilde{m}_{a}{}^{2}(t)\ \sin\theta(t)=0\;,$
(2.1)

with $H(t)$ the Hubble parameter (determined by the cosmology), and
$\tilde{m}_{a}{}(t)$ the time (temperature) dependent mass of the axion.
Usually the axion mass is written as

$\tilde{m}_{a}{}^{2}(T)=\dfrac{\chi(T)}{f_{a}{}^{2}}\;,$ (2.2)

with $\chi$ 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

Figure 1: The mass of the axion as a function of the temperature for
$f_{a}{}=10^{12}~{}{\rm GeV}$, using the data provided in ref. [25].

We assume that at very high temperatures $\tilde{m}_{a}{}\to 0$ (for the QCD
axion this is true; see Fig. (1)) – i.e. $\tilde{m}_{a}{}\ll H$ – at very
early times. Therefore, at the very early Universe, the EOM is444For the QCD
axion the PQ symmetry breaking scale determines whether there are domains with
different $\theta_{\rm 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].

$\left(\dfrac{d^{2}}{dt^{2}}+3H(t)\ \dfrac{d}{dt}\right)\theta(t)=0\;.$ (2.3)

which is solved by $\theta=\theta_{\rm
ini}{}+C\displaystyle\int_{0}^{t}dt^{\prime}\
\left(\dfrac{a(t^{\prime}=0)}{a(t^{\prime})}\right)^{3}$, with $\theta_{\rm
ini}{}$ a constant (some initial value) and $a$ the scale factor of the
Universe. That is, as the Universe expands, and as long as $\tilde{m}_{a}{}\ll
H$, $\theta\to\theta_{\rm ini}{}$, and $\dot{\theta}\to 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_{\rm ini}{}$) such that $\dot{\theta}|_{t=t_{\rm
ini}{}}=0$ and $\theta|_{t=t_{\rm ini}{}}=\theta_{\rm ini}{}$. This is the
most common case (there are exceptions to this; e.g. [27]), and it is what
MiMeS uses.

##### The shift symmetry

It is important to note that $\dot{\theta}\sim a^{-3}$ is expected from
symmetry grounds. The axion is assumed to be the phase of a complex scalar
field

$\Phi(x)=\dfrac{1}{\sqrt{2}}(f_{a}{}+\phi(x))\exp{(i\theta(x))}\;,$ (2.4)

charged under a global symmetry. If this symmetry is not broken, then the
Lagrangian is invariant under a shift $\theta\to\theta+z$. Then there exists a
Noether charge $Q\sim a^{3}\dot{\theta}=const.$ Thus, at high temperatures tha
shift symmetry is restored (i.e. the mass of the axion vanishes is much
smaller than its kinetic energy), and $\dot{\theta}\sim a^{-3}$.

On the other hand, if at high temperatures there exist shift symmetry breaking
interactions, then $a^{3}\dot{\theta}\neq const.$ This means that at low
temperatures, $\dot{\theta}$ can be sizeable. The explicit breaking of the
shift symmetry was introduced recently as a means to alter the axion
production mechanism (for details on how this affects the evolution of the
axion see e.g. [27, 28]). Practically, a non-vanishing $\dot{\theta}$ can
change the point where it starts to oscillate, which greatly affects its final
energy density.

From this point, we are not going to discuss a non-vanishing
$\dot{\theta}|_{t=t_{\rm ini}{}}$ further, since it is beyond the scope of
MiMeS.

###  2.1 The WKB approximation

In order to solve analytically eq. (2.1), we assume $\theta\ll 1$, which
results in the linearised EOM

$\left(\dfrac{d^{2}}{dt^{2}}+3H(t)\
\dfrac{d}{dt}+\tilde{m}_{a}{}^{2}(t)\right)\theta(t)=0\;.$ (2.5)

Using a trial solution $\theta_{\rm trial}=\exp\left[i\displaystyle\int dt\
\Big{(}\psi(t)+3/2\ i\ H(t)\Big{)}\right]$, and defining
$\Omega^{2}=\tilde{m}_{a}{}^{2}-\dfrac{9}{4}H^{2}-\dfrac{3}{2}\dot{H}$ we can
transform the eq. (2.5) to

$\psi^{2}=\Omega^{2}+i\ \dot{\psi}\;,$ (2.6)

which has a formal solution $\psi=\pm\sqrt{\Omega^{2}+i\dot{\psi}}$. In the
WKB approximation, we assume a slow time-dependence; that is
$\dot{\psi}\ll\Omega^{2}$ and $\dot{\Omega}\ll\Omega^{2}$. Then we can
approximate $\psi$ as

$\psi\approx\pm\Omega+\dfrac{i}{2}\dfrac{d\log\Omega}{dt}\;,$ (2.7)

which results in the general solution of eq. (2.5)

$\theta\approx\dfrac{1}{\sqrt{\Omega}}\exp\left(-\dfrac{3}{2}\int dt\
H\right)\left[A\cos\left(\int dt\ \Omega\right)+B\sin\left(\int dt\
\Omega\right)\right]\;.$ (2.8)

Applying, then, the initial conditions $\dot{\theta}|_{t=t_{\rm ini}{}}=0$ and
$\theta|_{t=t_{\rm ini}{}}\approx\theta_{\rm ini}{}$, we arrive at

$\theta(t)\approx\theta_{\rm ini}{}\sqrt{\dfrac{\Omega_{\rm
ini}}{\Omega(t)}}\left(\dfrac{a}{a_{\rm ini}{}}\right)^{-3/2}\
\cos\left(\int_{t_{\rm ini}{}}^{t}dt^{\prime}\ \Omega(t^{\prime})\right)\;.$
(2.9)

In order to further simplify this approximate result, we note that $\theta$
deviates from $\theta_{\rm ini}{}$ close to $t=t_{\rm osc}{}$ – corresponding
to $T=T_{\rm osc}{}$, the so-called “oscillation temperature" –
$\tilde{m}_{a}{}|_{t=t_{\rm osc}{}}=3H|_{t=t_{\rm osc}{}}$, which is defined
as the point at which the axion begins to oscillate. This observation allows
us to set $t_{\rm ini}{}=t_{\rm osc}{}$. Moreover, at $t>t_{\rm osc}{}$, we
approximate $\Omega\approx\tilde{m}_{a}{}$, as $H^{2}$ and $\dot{H}$ become
much smaller than $\tilde{m}_{a}{}^{2}$ quickly after $t=t_{\rm osc}{}$.
Finally, the axion angle takes the form

$\theta(t)\approx\theta_{\rm
osc}{}\left(\dfrac{3}{4}\right)^{1/4}\sqrt{\dfrac{\tilde{m}_{a}{}|_{t=t_{\rm
osc}{}}}{\tilde{m}_{a}{}(t)}}\left(\dfrac{a}{a_{\rm osc}{}}\right)^{-3/2}\
\cos\left(\int_{t_{\rm osc}{}}^{t}dt^{\prime}\
\tilde{m}_{a}{}(t^{\prime})\right)\;,$ (2.10)

where $\theta_{\rm osc}{}=\theta|_{t=t_{\rm osc}{}}$. This equation is further
simplified if we assume that $\theta_{\rm osc}{}\approx\theta_{\rm ini}{}$,
i.e.

$\theta(t)\approx\theta_{\rm
ini}{}\left(\dfrac{3}{4}\right)^{1/4}\sqrt{\dfrac{\tilde{m}_{a}{}|_{t=t_{\rm
osc}{}}}{\tilde{m}_{a}{}(t)}}\left(\dfrac{a}{a_{\rm osc}{}}\right)^{-3/2}\
\cos\left(\int_{t_{\rm osc}{}}^{t}dt^{\prime}\
\tilde{m}_{a}{}(t^{\prime})\right)\;.$ (2.11)

It is worth mentioning that the accuracy of this approximation depends, in
general, on $T_{\rm osc}{}$; it determines the difference between $\theta_{\rm
ini}{}$ and $\theta_{\rm osc}{}$, the deviation of $\dot{\theta}|_{t=t_{\rm
osc}{}}$ from $0$, and whether $\dot{\Omega}\ll\Omega^{2}$.

##### Axion energy density

In the small angle approximation, the energy density of the axion is

$\displaystyle\rho_{a}=\dfrac{1}{2}f_{a}{}^{2}\left[\dot{\theta}^{2}+\tilde{m}_{a}{}^{2}\theta^{2}\right]\;.$
(2.12)

For the relic abundance of axions, we need to calculate their energy density
at very late times. That is, $\dot{\tilde{m}}_{a}=0$, $\tilde{m}_{a}{}\gg H$
and $\dot{H}\ll H^{2}$. After some algebra, we obtain the approximate form of
the energy density as a function of the scale factor

$\displaystyle\rho_{a}\approx\dfrac{m_{a}{}}{2}\ f_{a}{}^{2}\ \theta_{\rm
ini}{}^{2}\ \tilde{m}_{a}{}(a_{\rm osc}{})\ \left(\dfrac{a_{\rm
osc}{}}{a}\right)^{3}\;,$ (2.13)

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_{\rm osc}{}$,
the axion energy density gets diluted, since

$a^{3}\ s=\gamma\ a_{\rm osc}{}^{3}\ s_{\rm osc}\Rightarrow\left(\dfrac{a_{\rm
osc}{}}{a}\right)^{3}=\gamma^{-1}\dfrac{s}{s_{\rm osc}}\;,$ (2.14)

with $\gamma$ the amount of entropy injection to the plasma between $t_{\rm
osc}{}$ and $t$. Therefore, the present (at $T=T_{0}$) energy density of the
axion, becomes

$\displaystyle\rho_{a,0}=\gamma^{-1}\dfrac{s_{0}}{s_{\rm osc}}\ \dfrac{1}{2}\
f_{a}{}^{2}\ m_{a}{}\ \tilde{m}_{a}{}_{,{\rm osc}}\ \theta_{\rm ini}{}^{2}\;,$
(2.15)

with $m_{a}{}$ the mass of the axion at $T=T_{0}$. Notice that the explicit
dependence on $f_{a}{}$ cancels if $\tilde{m}_{a}{}=\sqrt{\chi(T)}/f_{a}{}$.
That is, $f_{a}{}$ only affects the energy density of the axions through its
impact on $T_{\rm osc}{}$.

###  2.2 Notation

The EOM (2.1) depends on time, which is not useful variable in cosmology,
especially non-standard cosmologies. Therefore, we introduce

$\displaystyle u=\log\dfrac{a}{a_{\rm ini}{}}\;,$ (2.16)

which results in

$\displaystyle\dfrac{dF}{dt}$ $\displaystyle=H\dfrac{dF}{du}$ (2.17)
$\displaystyle\dfrac{d^{2}F}{dt^{2}}$ $\displaystyle=H^{2}\
\left(\dfrac{d^{2}F}{du^{2}}+\dfrac{1}{2}\dfrac{d\log
H^{2}}{du}\dfrac{dF}{du}\right)\;.$

The EOM in terms of $u$, then, becomes

$\dfrac{d^{2}\theta}{du^{2}}+\left[\dfrac{1}{2}\dfrac{d\log
H^{2}}{du}+3\right]\dfrac{d\theta}{du}+\
\left(\dfrac{\tilde{m}_{a}{}}{H}\right)^{2}\ \sin\theta=0\;.$ (2.18)

Notice that in a radiation dominated Universe

$\dfrac{d\log H^{2}}{du}=-\left(\dfrac{d\log g_{\rm eff}{}}{d\log
T}+4\right)\delta_{h}^{-1}\;,$

with $\delta_{h}=1+\dfrac{1}{3}\dfrac{d\log h_{\rm eff}{}}{d\log T}$. In a
general cosmological setting, if the expansion rate is dominated by an energy
density that scales as $\rho\sim a^{-c}$, $\dfrac{d\log H^{2}}{du}=-c$. We
notice that close to rapid particle annihilations and decays, the evolution of
the energy densities change, and $\dfrac{d\log H^{2}}{du}$ can only be
computed numerically.

Moreover, it is worth mentioning that the EOM 2.18 with the initial condition
$\theta(u=0)=\theta_{\rm ini}{}$ and $d\theta/du(u=0)=0$ can be written as a
system of first order ordinary differential equations

$\displaystyle\dfrac{d\zeta}{du}+\left[\dfrac{1}{2}\dfrac{d\log
H^{2}}{du}+3\right]\zeta+\ \left(\dfrac{\tilde{m}_{a}{}}{H}\right)^{2}\
\sin\theta=0\;.$ $\displaystyle\dfrac{d\theta}{du}-\zeta=0\;.$ (2.19)

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.

###  2.3 Adiabatic invariant and the anharmonic factor

The EOM 2.1 ca be solved analytically in the approximation $\theta\ll 1$.
Moreover, even for $\theta\ll 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 $\mathcal{H}(\theta,p;t)$, the equations of
motion are

$\dot{p}=-\dfrac{\partial\mathcal{H}}{\partial\theta}\;,\;\;\dot{\theta}=\dfrac{\partial\mathcal{H}}{\partial
p}\;.$ (2.20)

Moreover, we note that

$d\mathcal{H}=\dot{\theta}\ dp-\dot{p}\
d\theta+\dfrac{\partial\mathcal{H}}{\partial t}\ dt\;.$ (2.21)

If this system exhibits closed orbits (e.g. if it oscillates), we define

$J\equiv C\ \oint p\ d\theta\;,$ (2.22)

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,

$\dfrac{dJ}{dt}=C\ \oint\Big{(}\dot{p}\ d\theta+p\ d\dot{\theta}\Big{)}=C\
\displaystyle\int_{t}^{t+T}\dfrac{\partial\mathcal{H}}{\partial t^{\prime}}\
dt^{\prime}\approx T\ \dfrac{\partial\mathcal{H}(t^{\prime})}{\partial
t^{\prime}}\Big{|}_{t^{\prime}=t}\approx 0\;.$

##### Application to the axion

The Hamiltonian that results in the EOM of eq. (2.1) is

$\mathcal{H}=\dfrac{1}{2}\dfrac{p^{2}}{f_{a}{}^{2}\ a^{3}}+V(\theta)\
a^{3}\;,$ (2.23)

with

$\displaystyle p=f_{a}{}^{2}\ a^{3}\ \dot{\theta}$ (2.24) $\displaystyle
V(\theta)=\tilde{m}_{a}{}^{2}f_{a}{}^{2}(1-\cos\theta)\;.$ (2.25)

Notice that the Hamiltonian varies slowly if
$\dot{\tilde{m}}_{a}(T)/\tilde{m}_{a}{}\ll\tilde{m}_{a}{}$ and
$H\ll\tilde{m}_{a}{}$, which are the adiabatic conditions. When these
conditions are met, the adiabatic invariant for this system becomes

$J=\dfrac{\oint p\ d\theta}{\pi f_{a}{}^{2}}=\dfrac{1}{\pi
f_{a}{}^{2}}\oint\sqrt{2\left(\mathcal{H}(\theta)-V(\theta)\ a^{3}\right)\
f_{a}{}^{2}a^{3}\ }\ d\theta=\dfrac{2}{\pi f_{a}{}^{2}}\int_{-\theta_{\rm
peak}{}}^{\theta_{\rm peak}{}}\sqrt{2\left(\mathcal{H}(\theta_{\rm
peak}{})-V(\theta)\ a^{3}\right)\ f_{a}{}^{2}a^{3}\ }d\theta\;,$ (2.26)

where we note that $\theta_{\rm peak}{}$ denotes the maximum of $\theta$ – the
peak of the oscillation, which corresponds to $p=0$. That is,
$\mathcal{H}(\theta_{\rm peak}{})=V(\theta_{\rm peak}{})\ a^{3}$. Therefore,
the adiabatic invariant, takes the form

$\displaystyle J=$ $\displaystyle\dfrac{2\sqrt{2}}{\pi
f_{a}{}}\int_{-\theta_{\rm peak}{}}^{\theta_{\rm peak}{}}\sqrt{V(\theta_{\rm
peak}{})-V(\theta)}a^{3}d\theta=\dfrac{2\sqrt{2}}{\pi}\
\tilde{m}_{a}{}\,a^{3}\ \displaystyle\int_{-\theta_{\rm peak}{}}^{\theta_{\rm
peak}{}}\sqrt{\cos\theta-\cos\theta_{\rm peak}{}}\ d\theta\;,$ (2.27)

where, for the last equality. we have used the adiabatic conditions, i.e.
negligible change of $\tilde{m}_{a}{}$ and $a$ during one period. Usually, the
adiabatic invariant is written as [29, 30]

$J=a^{3}\ \tilde{m}_{a}{}\ \theta_{\rm peak}{}^{2}\,f(\theta_{\rm peak}{})\;,$
(2.28)

where

$f(\theta_{\rm peak}{})=\dfrac{2\sqrt{2}}{\pi\theta_{\rm
peak}{}^{2}}\displaystyle\int_{-\theta_{\rm peak}{}}^{\theta_{\rm
peak}{}}d\theta\sqrt{\cos\theta-\cos\theta_{\rm peak}{}}\;,$ (2.29)

is called the anharmonic factor, with $0.5\lesssim f(\theta_{\rm peak}{})\leq
1$ (see Fig. (2)).

Figure 2: The anharmonic factor for $0\leq\theta_{\rm peak}{}<\pi$.

##### 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
$\theta$ 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 $\theta$. Once the adiabatic conditions are fulfilled, we can stop
the integration at a peak, $\theta_{\rm peak}{}_{,*}$ – which corresponds to
$T=T_{*}$ and $a=a_{*}$. Then, the value of the maximum angle today
($\theta_{\rm peak}{}_{,0}\ll 1$) is related to $\theta_{\rm peak}{}_{,*}$ via

$\displaystyle\theta_{\rm peak}{}_{,0}^{2}$
$\displaystyle=\left(\dfrac{a_{*}}{a_{0}}\right)^{3}\
\dfrac{\tilde{m}_{a}{}_{,*}}{m_{a}{}}\ f(\theta_{\rm peak}{}_{,*})\
\theta_{\rm peak}{}_{,*}^{2}=\gamma^{-1}\ \dfrac{s_{0}}{s_{*}}\
\dfrac{\tilde{m}_{a}{}_{,*}}{m_{a}{}}\ f(\theta_{\rm peak}{}_{,*})\
\theta_{\rm peak}{}_{,*}^{2}\;.$ (2.30)

Plugging this into eq. (2.12) (with $\dot{\theta}=0$, i.e. at today’s peak),
we arrive at the energy density today

$\rho_{a,0}=\gamma^{-1}\ \dfrac{s_{0}}{s_{*}}\ m_{a}{}\ \tilde{m}_{a}{}_{,*}\
\dfrac{1}{2}\ f_{a}{}^{2}\ \theta_{\rm peak}{}_{,*}^{2}\;\ f(\theta_{\rm
peak}{}_{,*})\;,$ (2.31)

where $\gamma$ is the entropy injection coefficient between $T_{*}$ and
$T_{0}$, defined from

$a^{3}(T_{0})\ s(T_{0})=\gamma\ a^{3}(T_{*})\ s(T_{*})\;$ (2.32)

Notice that eq. (2.31) is similar to the corresponding WKB result (2.15) at
$t_{\rm osc}{}\to t_{*}$, multiplied by the anharmonic factor $f(\theta_{\rm
peak}{}_{,*})$. So, the numerical integration is needed in order to correctly
identify $T_{*}$ and $\theta_{\rm peak}{}_{,*}$, which are greatly affected by
the underlying cosmology (especially in cases where entropy injection is
active close to $T_{\rm osc}{}$) and whether $\theta_{\rm ini}{}\gtrsim 1$
[30, 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.

##  3 MiMeS usage

The latest stable version of MiMeS is available at
github.com/dkaramit/MiMeS/tree/stable, which can also be obtained by running:
555Instructions on how git can be installed can be found in
https://github.com/git-guides/install-git.

⬇

1 git clone -b stable https://github.com/dkaramit/MiMeS.git

Moreover, one can download MiMeS from mimes.hepforge.org, or by running

⬇

1 wget -c https://mimes.hepforge.org/downloads/MiMeS-x.y.z.tar.gz -O - | tar -xz

where x.y.z define the version of the library (currently this is 1.0.0).

It is important to note that MiMeS relies on NaBBODES [31] and SimpleSplines
[32]. 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

2 cd MiMeS

3 git submodule init

4 git submodule update \--remote

This downloads the master branch of MiMeS; and NaBBODES [31] and SimpleSplines
[32] 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.3 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 $\theta_{\rm ini}{}$ and $f_{a}{}$, the oscillation
temperature and the corresponding value of $\theta$, the evolution of the
axion (e.g. temperature, $\theta$, $\rho_{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.

###  3.1 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.

####  3.1.1 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 $\chi(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 ${\rm GeV}$) and $\chi$ (in ${\rm 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 $\geq$ minT and TMax $\leq$ maxT. Beyond these
limits, the mass is assumed to constant by default. However, one can change
this by using

⬇

1 axionMass.set_ma2_MAX(ma2_MAX);

2 axionMass.set_ma2_MIN(ma2_MIN);

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\geq$TMax and $T\leq$TMin,
respectively. Usually, in order to ensure continuity of the axion mass, one
needs to know TMin, TMax, $\chi($TMin$)$, and $\chi($TMax$)$; which can be
found by calling axionMass.getTMin(), axionMass.getTMax(),
axionMass.getChiMin(), and axionMass.getChiMax(), respectively.

##### Axion mass from a 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

⬇

1 mimes::AxionMass<LD> axionMass(ma2);

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

⬇

1 mimes::Axion<LD,Solver,Method> ax(theta_i, fa, umax, TSTOP, ratio_ini,

2 N_convergence_max, convergence_lim, inputFile, &axionMass,
initial_step_size,

3 minimum_step_size, maximum_step_size, absolute_tolerance,
relative_tolerance,

4 beta, fac_max, fac_min, maximum_No_steps);

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.

The various parameters are as follows:

  1. 1.

theta_i: Initial angle.

  2. 2.

fa: The PQ scale.

  3. 3.

umax : If $u>$umax the integration stops (remember that $u=\log(a/a_{i})$).
Typically, this should be a large number ($\sim 1000$), in order to avoid
stopping the integration before the axion begins to evolve adiabatically.

  4. 4.

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

  5. 5.

ratio_ini: Integration starts when $3H/\tilde{m}_{a}{}\approx$ratio_ini (the
exact point depends on the file ‘‘inputFile", which we will see later).

  6. 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. 7.

inputFile: Relative (or absolute) path to a file that describes the cosmology.
the columns should be: $u$ $T~{}[{\rm GeV}]$ $\log H$, sorted so that $u$
increases. 666One can run “bash MiMeS/src/FormatFile.sh inputFile” in order to
sort it and remove any unwanted duplicates. See Appendix E for details of
MiMeS/src/FormatFile.sh. 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.

  8. 8.

axionMass: An instance of the mimes::AxionMass<LD> class, passed by pointer.

  9. 9.

initial_stepsize (optional): Initial step the solver takes.

  10. 10.

maximum_stepsize (optional): This limits the step-size to an upper limit.

  11. 11.

minimum_stepsize (optional): This limits the step-size to a lower limit.

  12. 12.

absolute_tolerance (optional): Absolute tolerance of the RK solver

  13. 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.3.

  14. 14.

beta (optional): Controls how agreesive the adaptation is. Generally, it
should be around but less than 1.

  15. 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$=\infty$ and fac_min$=0$, but in reality one must tweak them in order
to avoid instabilities.

  16. 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 RK methods are discussed in Appendix A.

The EOM (2.18), then can be solved using

⬇

1 ax.solveAxion();

Once the EOM is solved, we can access $T_{\rm osc}{}$, $\theta_{\rm osc}{}$,
and $\Omega 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~{}[{\rm GeV}]$, $\theta$, $d\theta/du$, $\rho_{a}$.
Moreover, the peaks of he oscillation are stored in another two-dimensional
std::vector<LD>, with the columns corresponding to $a$, $T~{}[{\rm GeV}]$,
$\theta_{\rm peak}{}$, $d\theta/du=0$, $\rho_{a}$, $J$. We should note that
the peaks are identified using linear interpolation between integration
points, in order to ensure that $d\theta/du=0$. That is, the values stored in
ax.peaks do not exist in ax.points.

As already mentioned, MiMeS uses embedded RK methods in order to solve the
axion EOM. That is, each integration point is calculated twice using two
estimates of different order. The difference between the two estimates is then
interpreted as the local integration error (see Appendix A for details). These
local errors at the integration points of $\theta$ and $\zeta$ 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

⬇

1 ax.set_ma2(ma2);

with ma2 a callable object with signature LD ma2(LD T, LD fa).

However, since integration starts at a temperature determined by ratio_ini, if
the mass changes (including the definitions beyond the interpolation limits),
we have to remake the interpolation of the underlying cosmology described by
the parameter inputFile. Thus, if the definition of the mass changes, we need
to call

⬇

1 ax.restart();

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 $\theta_{\rm ini}{}$ without generating another instance.
777Since the interpolations of the data of inputFile are made inside the
constructor of the mimes::Axion<LD,Solver,Method> class,
mimes::Axion<LD,Solver,Method>::setTheta_i is a faster choice if ones needs to
solve the EOM for a different initial condition. This function is used as

⬇

1 ax.setTeta_i(new_theta_ini);

where new_theta_ini is the new value of $\theta_{\rm 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.

####  3.1.2 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.2, 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

⬇

1 axionMass.set_ma2_MAX(ma2_MAX)

2 axionMass.set_ma2_MIN(ma2_MIN)

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 return $\tilde{m}_{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, $\chi($TMin$)$, and $\chi($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.

##### The EOM solver

We can define an Axion instance as follows

⬇

1 ax=Axion(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max,
convergence_lim,

2 inputFile, axionMass, initial_step_size, minimum_step_size,
maximum_step_size,

3 absolute_tolerance, relative_tolerance, beta, fac_max, fac_min,
maximum_No_steps)

Here the input parameters are the same as in the C++ case, and outlined in
table 2. Moreover, the usage of the class can be found by running Axion after
loading the module. The only slight difference compared to the C++ case is
that the axionMass instance is not passed as a pointer; it is done internally
using ctypes.

Using the defined variable (ax in this example), we can simply run

⬇

1 ax.solveAxion()

in order to solve the EOM of the axion. In contrast to the C++ implementation,
this only gives us access to $T_{\rm osc}{}$, $\theta_{\rm osc}{}$, and
$\Omega h^{2}$; the corresponding variables are ax.T_osc, ax.theta_osc, and
ax.relic. In order to get the evolution of the axion field, we need to run

⬇

1 ax.getPoints()

This will make numpy [33] arrays that contain the scale factor (ax.a),
temperature (ax.T), $\theta$ (ax.theta), its derivative with respect to $u$
(ax.zeta), and the energy density of the axion (ax.rho_axion).

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), $\theta$ (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 $\theta$ and $\zeta$ 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

⬇

1 ax.set_ma2(ma2)

with ma2 a function that takes $T$ and $f_{a}{}$, and returns the
$\tilde{m}_{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 $\theta_{\rm 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.

##  4 Assumptions and user input

###  4.1 Restrictions

MiMeS only makes a few, fairly general, assumptions. If these assumptions are
violated, then MiMeS might not work as expected.

First of all, it is assumed that the axion energy density is always
subdominant compared to radiation or any other dominant component of the
Universe, and that decays and annihilations of particles have a negligible
effect on the axion energy density.

Moreover, the initial condition -- discussed in detail in section 2 \-- is
always assumed to be $\theta_{t=t_{\rm ini}{}}=\theta_{\rm ini}{}$ and
$\dot{\theta}|_{t=t_{\rm ini}{}}=0$. This cannot change, since it is an
integral part of how MiMeS finds a suitable starting point. Furthermore, it is
also assumed that $3H/\tilde{m}_{a}{}$ increases monotonically at high
temperatures. 888This assumption is connected to the initial condition, as a
sizeable mass at high temperatures can break the shift symmetry mentioned in
section 2.

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.

###  4.2 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, explained in Appendix B.
The various choices we for the shared libraries used by the python interface
are given in MiMeS/Definitions.mk while the corresponding options for the C++
examples are in the Definitions.mk files inside the subdirectories of
MiMeS/UserSpace/Cpp. The options correspond to different variables, which are

  1. 1.

rootDir: The relative path of root directory of MiMeS.

  2. 2.

LONG: This sets the numeric types for the C++ examples. It should be either
long or omitted. 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 $\sim 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. 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.

  4. 4.

SOLVER: MiMeS uses the ordinary differential equation (ODE) integrators of
ref. [31]. 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.

  5. 5.

METHOD: Depending on the type of solver, there are some available methods.
999It is worth mentioning that NaBBODES is built in order to be a template for
all possible Rosenbrock and explicit Runge-Kutta embedded methods, and one can
provide their own Butcher tableau if they want to use another method, as shown
in Appendix A.

     * •

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. [34]). The ROS34PW2 choice corresponds to a
third order Rosenbrock-Wanner method [35].

     * •

For SOLVER=2, the only reliable method available in NaBBODES is the Dormand-
Prince [36] 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. 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. 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. 101010In C++ the template
arguments are part of the definition of a class; if the template arguments are
not known, the class is not even constructed. 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.

###  4.3 User input

####  4.3.1 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 [37], and a few points for
$f(\theta_{\rm peak}{})$ introduced in eq. (2.29), 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 ${\rm GeV}$), $h_{\rm
eff}{}$, and $g_{\rm eff}{}$.

  * •

The axion mass data must be given in two columns; $T$ (in ${\rm GeV}$), $\chi$
(in ${\rm GeV}^{4}$). Here, $\chi$ 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 $\theta_{\rm
peak}{}$ $f(\theta_{\rm peak}{})$; with increasing $\theta_{\rm peak}{}$.

The paths to these files should be given at compile time. That is, once
Paths.mk changes, we must run ‘‘bash configure.sh" and then ‘‘make" in order
to make sure that they will be used. 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,
$\chi(T)$, and $f(\theta_{\rm 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. 111111Absolute
paths have the advantage to be accessible from everywhere else in the system.
Thus, executables that seek the corresponding files can be called and copied
easily.

####  4.3.2 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, $\theta_{\rm 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$=\infty$ 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. 121212These two variables are not optional,
because the user must be aware of them, in order to choose them according to
their needs.

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 conditions -- i.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 $\theta$ tends to oscillate
rapidly which destabilizes the differential equation solver. Therefore, these
parameters should be chosen carefully, in order to ensure that integration
stops when the axion has reached its adiabatic evolution, without
destabilising the EOM.

###  4.4 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, 131313In terms of the
parametrisation introduced in ref. [38, 39], for this scenario we choose
$T_{\rm end}=10^{-2}~{}{\rm GeV},\;c=3,\;T_{\rm ini}=10^{12}~{}{\rm
GeV},\;\text{and}\;r=0.1$.. with the evolution of the energy densities of the
plasma and the matter field ($\Phi$) shown in Fig. (3).

Figure 3: The energy densities of the plasma and a hypothetical decaying
fluid. The parameters of the cosmological scenario $T_{\rm END}=10^{-2}~{}{\rm
GeV}$, $c=3$, $T_{\rm ini}=10^{12}~{}{\rm GeV}$, and $r=0.1$; the parameters
follow the definitions in ref. [39].

The various regime changes are indicated as $T_{{\rm E}_{1,2}}$ \--
corresponding to the first and second time where $\rho_{\Phi}=\rho_{\rm R}$,
and $T_{{\rm D}_{1,2}}$ \-- which correspond to the start and end of entropy
injection (for more precise definition see [39]). 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_{\rm min,\,max}$, in the
corresponding data file, we will take

$\tilde{m}_{a}{}=\Bigg{\\{}\begin{matrix}\dfrac{\chi(T_{\rm
min})}{f_{a}{}^{2}}&\text{for }T<T_{\rm min}\\\ \\\ \dfrac{\chi(T_{\rm
max})}{f_{a}{}^{2}}\left(\dfrac{T}{T_{\rm max}}\right)^{-8.16}&\text{for
}T>T_{\rm max}\end{matrix}\;.$ (4.1)

####  4.4.1 complete example in C++

In order to write a C++ program that uses MiMeS in order to solve the EOM
2.19, 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.

2 mimes::AxionMass<long double> axionMass(chi_PATH,0,mimes::Cosmo<long
double>::mP);

3

4 /* 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 */

7

8 /*set ma2 for T>TMax*/

9 long double TMax=axionMass.getTMax();

10 long double chiMax=axionMass.getChiMax();

11 axionMass.set_ma2_MAX(

12 [&chiMax,&TMax](long double T, long double fa){

13 return chiMax/fa/fa*std::pow(T/TMax,-8.16);

14 }

15 );

16

17 /*set ma2 for T<TMin*/

18 long double TMin=axionMass.getTMin();

19 long double chiMin=axionMass.getChiMin();

20 axionMass.set_ma2_MIN(

21 [&chiMin,&TMin](long double T, long double fa){

22 return chiMin/fa/fa;

23 }

24 );

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

⬇

1 std::string inputFile = std::string(rootDir)+

2 std::string("/UserSpace/InputExamples/MatterInput.dat");

3

4 mimes::Axion<long double, 1, RODASPR2<long double> >

5 ax(0.1, 1e16, 500, 1e-4, 1e3, 15, 1e-2, inputFile, &axionMAss,

6 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:

  * •

theta_i$=0.1$ • fa$=10^{16}~{}{\rm GeV}$ • umax $=500$ • TSTOP$=10^{-4}~{}{\rm
GeV}$ • ratio_ini$=10^{3}$ • N_convergence_max$=15$

  * •

convergence_lim$=10^{-2}$ • initial_stepsize$=10^{-1}$ •
maximum_stepsize$=10^{-1}$ • minimum_stepsize$=10^{-8}$ •
absolute_tolerance$=10^{-11}$ • 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

⬇

1 g++ -O3 -std=c++17 -lm -I../ -o axion axionExample.cpp

or

⬇

1 clang -lstdc++ -O3 -std=c++17 -lm -I../ -o axion axionExample.cpp

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.

##### The entire code

This example consists of only a few lines of code, which, including the change
in the axion mass beyond the interpolation limit, is

⬇

1 #include<iostream>

2 //include everything you need from MiMeS

3 #include"src/Axion/AxionSolve.hpp"

4 #include"src/AxionMass/AxionMass.hpp"

5 #include"src/Cosmo/Cosmo.hpp"

6 #include"src/misc_dir/path.hpp"

7

8 int main(){

9 // use chi_PATH to interpolate the axion mass.

10 mimes::AxionMass<long double> axionMass(chi_PATH,0,mimes::Cosmo<long
double>::mP);

11

12 /*set $\tilde{m}_{a}{}^{2}$ for $T\geq T_{\rm max}$*/

13 long double TMax=axionMass.getTMax();

14 long double chiMax=axionMass.getChiMax();

15 axionMass.set_ma2_MAX(

16 [&chiMax,&TMax](long double T, long double fa){

17 return chiMax/fa/fa*std::pow(T/TMax,-8.16);

18 }

19 );

20

21 /*set $\tilde{m}_{a}{}^{2}$ for $T\leq T_{\rm min}$*/

22 long double TMin=axionMass.getTMin();

23 long double chiMin=axionMass.getChiMin();

24 axionMass.set_ma2_MIN(

25 [&chiMin,&TMin](long double T, long double fa){

26 return chiMin/fa/fa;

27 }

28 );

29

30 std::string inputFile = std::string(rootDir)+

31 std::string("/UserSpace/InputExamples/MatterInput.dat");

32

33 mimes::Axion<long double, 1, RODASPR2<long double> >

34 ax(0.1, 1e16, 500, 1e-4, 1e3, 15, 1e-2, inputFile, &axionMass,

35 1e-1, 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7));

36

37 ax.solveAxion();

38

39 std::cout<<ax.relic<<"\n";

40

41 return 0;

42 }

We should point out that another example is given in
MiMeS/UserSpace/Cpp/Axion, where all parameters are taken as command-line
inputs, the various compilation-time options can be given
MiMeS/UserSpace/Cpp/Axion/Definitions.mk, and then compiled using make.
Therefore, the user can just modify this code in order to meet their needs.
That is, using this example, one only needs to add their preferred definition
of ma2_MAX and ma2_MIN -- or change the axionMass variable (in
MiMeS/src/static.hpp) using a function as mentioned in section 4.3.1 without
writing the entire program themselves.

##### Alternative axion mass definition

For this particular example, we could have used the approximate definition of
the axion mass

$\tilde{m}_{a}{}^{2}\approx
m_{a}{}^{2}\times\Bigg{\\{}\begin{matrix}\left(\dfrac{T_{\rm
QCD}}{T}\right)^{8.16}&\quad\text{for }T\geq T_{\rm QCD}\\\ \\\
1&\quad\text{for }T\leq T_{\rm QCD}\end{matrix}\;,$ (4.2)

where $T_{\rm QCD}\approx 150~{}{\rm MeV}$ and $m_{a}{}=\dfrac{3.2\times
10^{-5}~{}{\rm GeV}^{4}}{f_{a}{}^{2}}$. This can be done by substituting the
axionMass definition (lines $9$-$28$) with

⬇

1 auto ma2 = [](long double T,long double fa){

2 long double TQCD=150*1e-3;

3 long double ma20=3.1575e-05/fa/fa;

4 if(T<=TQCD){return ma20;}

5 else{return ma20*std::pow((TQCD/T),8.16);}

6 };

7 mimes::AxionMass<long double> axionMass(ma2);

####  4.4.2 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

⬇

1 #add the relative path for MiMeS/src

2 from sys import path as sysPath

3 sysPath.append(’../src’)

4

5 from interfacePy.AxionMass import AxionMass #import the AxionMass class

6 from interfacePy.Axion import Axion #import the Axion class

7 from interfacePy.Cosmo import mP #import the Planck mass

Once everything we need is imported, we can simply follow the steps outlined
in section 3.1.2. For the example at hand, we can create an instance of the
AxionMass class as

⬇

1 # AxionMass instance

2 axionMass = AxionMass(r’../src/data/chi.dat’,0,mP)

3

4 # 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

6 TMax=axionMass.getTMax()

7 chiMax=axionMass.getChiMax()

8 TMin=axionMass.getTMin()

9 chiMin=axionMass.getChiMin()

10

11 axionMass.set_ma2_MAX( lambda T,fa: chiMax/fa/fa*pow(TMax/T,8.16))

12 axionMass.set_ma2_MIN( lambda T,fa: chiMin/fa/fa )

Similar to the C++ case, the second and third arguments are $0$ and mP (the
Planck mass). However, in the python interface, mP needs to be imported from
the Cosmo module.

Then we can simply create an instance of the Axion class as

⬇

1 #in python it is more convenient to use relative paths

2 inputFile = "./InputExamples/MatterInput.dat"

3

4 ax = Axion(0.1, 1e16, 500, 1e-4, 1e3, 15, 1e-2, inputFile, axionMass,

5 1e-1, 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7))

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
$\theta_{\rm ini}{}$, $f_{a}{}$, $\theta_{\rm osc}{}$, $T_{\rm osc}{}$, and
$\Omega 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

⬇

1 #this gives you all the points of integration

2 ax.getPoints()

3

4 #this gives you the peaks of the oscillation

5 ax.getPeaks()

6

7 #this gives you local integration errors

8 ax.getErrors()

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

⬇

1 del ax

2 del axionMass

##### Compilation of the shared library

As described in section 4.3.1, we may need to change the default data file
paths, or the various compilation options. This is done through the variables
in MiMeS/Definitions.mk and MiMeS/Paths.mk described in section 4.2.

Once we have chosen everything according to our needs, the library can be
created by opening a terminal inside the root directory of MiMeS and running

⬇

1 bash configure.sh

2 make lib/Axion_py.so

##### The entire code

As in the C++ example, this example consists of only a few lines of code. The
script we described here is

⬇

1 #add the relative path for MiMeS/src

2 from sys import path as sysPath

3 sysPath.append(’../src’)

4

5 from interfacePy.AxionMass import AxionMass #import the AxionMass class

6 from interfacePy.Axion import Axion #import the Axion class

7 from interfacePy.Cosmo import mP #import the Planck mass

8

9 # AxionMass instance

10 axionMass = AxionMass(r’../src/data/chi.dat’,0,mP)

11

12 # define $\tilde{m}_{a}{}^{2}$ for $T\leq T_{\rm min}$

13 TMin=axionMass.getTMin()

14 chiMin=axionMass.getChiMin()

15 axionMass.set_ma2_MIN( lambda T,fa: chiMin/fa/fa )

16

17 # define $\tilde{m}_{a}{}^{2}$ for $T\geq T_{\rm max}$

18 TMax=axionMass.getTMax()

19 chiMax=axionMass.getChiMax()

20 axionMass.set_ma2_MAX( lambda T,fa: chiMax/fa/fa*pow(TMax/T,8.16))

21

22 #in python it is more convenient to use relative paths

23 inputFile = inputFile="./InputExamples/MatterInput.dat"

24

25 ax = Axion(0.1, 1e16, 500, 1e-4, 1e3, 15, 1e-2, inputFile, axionMass,

26 1e-1, 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7))

27

28 ax.solveAxion()

29

30 print(ax.relic)

31

32 #once we are done we should run the destructor

33 del ax

34 del axionMass

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 [40], 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. (4.2), we could
replace the definition of the axionMass variable (lines $10$-$21$) with

⬇

1 def ma2(T,fa):

2 TQCD=150*1e-3

3 ma20=3.1575e-05/fa/fa

4 if T<=TQCD:

5 return ma20;

6 return ma20*(TQCD/T)**8.16

7 axionMass = AxionMass(ma2)

####  4.4.3 Results

(a)

(b)

Figure 4: The evolution of the axion angle, $\theta$, with the temperature in
early matter dominated (a) and radiation dominated (b) cases. For the EMD
case, the $\theta_{\rm ini}{}$ is chosen so that $\Omega h^{2}=0.12$.

In Fig. (4(a)) we show the evolution of $\theta$ for temperatures
$T\in[0.025,0.15]~{}{\rm GeV}$, where the vertical line indicates $T_{\rm
osc}{}$, while the horizontal one the corresponding value of $\theta$,
$\theta_{\rm osc}{}$. Also, the blue curve connects the peaks of the
oscillation. For comparison, in Fig. (4(b)) 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_{\rm 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~{}{\rm
GeV}<T<0.15~{}{\rm GeV}$.

(a)

(b)

Figure 5: (a) Local relative errors of $\theta$ and $\zeta$ for the
temperature range $0.025{\rm GeV}\leq T\leq 0.15{\rm GeV}$. (b) Histogram of
number of points in the same temperature range.

The local relative errors, defined in Appendix A, are shown in Fig. (5(a)).
The black and red lines correspond to $\delta\theta/\theta$ and
$\delta\zeta/\zeta$, respectively. This figure indicates that the local errors
are relatively well behaved for $T>T_{\rm 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 $\sim
10^{-7}$ (the initial relative error of $\zeta$ is large because $\zeta\approx
0$ before the $T_{\rm 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.
141414More precisely, the histogram is made by dividing the temperature range
$0.025~{}{\rm GeV}\lesssim T\lesssim 0.15~{}{\rm GeV}$ to $30$ bins of equal
size. 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 -- $\tilde{m}_{a}{}/H$ increases rapidly for $T<T_{\rm 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\lesssim T_{\rm 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.

##  5 Acknowledgements

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

##  6 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.4)). 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 [41] 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 $\zeta$ 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

##  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.
151515Boundary value problems, and higher order differential equations are
expressed as first order ODEs, and then solved. Similarly to eq. (2.19).
Although there are some very insightful sources in the literature (e.g. [42,
43, 44]) 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

$\dfrac{d\vec{y}}{dt}=\vec{f}(\vec{y},t)\;,$ (A.1)

with given initial condition $\vec{y}(0)$. Also, the components of $\vec{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

$\vec{y}_{n+1}=\vec{y}_{n}+h\sum_{i=1}^{s}b_{i}\ \vec{k}_{i}\;,$ (A.2)

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 $\vec{k}_{i}$
define the corresponding RK method. For example, the classic Euler method is
an RK method with $s=1$, $b=1$, and $\vec{k}_{1}=\vec{f}(\vec{y_{n},t_{n}})$.
Methods with $\vec{k}$ that depends on previous step (i.e.
$\vec{k}=\vec{k}\left[y_{n},t_{n}\right]$), are called explicit, while the
ones that try to also predict next step (i.e.
$\vec{k}=\vec{k}\left[y_{n},t_{n};y_{n+1},t_{n+1}\right]$) are called
implicit. 161616Generally, by substituting implicit methods $y_{n+1}$ as in
eq. (A.2), we end up with a system of equations that need to be solved in
order to compute $\vec{k}$.

###  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 $\vec{k}$ many times within the same iteration. Therefore, together
with the iteration of eq. (A.2), a second estimate is given by

$\vec{y}_{n+1}^{*}=\vec{y}_{n}+h\sum_{i=1}^{s}b_{i}^{*}\vec{k}_{i}\;,$ (A.3)

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

$\Delta\equiv\sqrt{\frac{1}{N}\sum_{d=1}^{N}\left(\frac{y_{n+1,d}-y^{*}_{n+1,d}}{\Lambda_{d}}\right)^{2}}\;,$
(A.4)

where $n$ the iteration number, $N$ the number of ODEs, $d$ the component of
$\vec{y}$, and $\Lambda_{d}$ defined as

$\Lambda_{d}={\tt Atol}+\max\left(|y_{n+1,d}|,|y^{*}_{n+1,d}|\right)\;{\tt
Rtol}\;,$ (A.5)

with Atol and Rtol the absolute and relative tolerances that characterise the
desirable accuracy we want to achieve; user defined values, typically
Atol$=$Rtol$\ll 1$. With these definitions, the desirable error is reached
when $\Delta\lesssim 1$.

##### Step-control

The definition A.4, allows us to adjust the step-size $h$ in such way that
$\Delta\approx 1$. That is, we take trial steps, and $h$ is adapted until
$\Delta\approx 1$. A simple adaptive strategy adjusts step-size, using

$\displaystyle h\to\beta\;h\;{\rm max}\left[f_{\rm min},{\rm
min}\left(\Delta^{-\frac{1}{p+1}},f_{\rm max}\right)\right]\;,$ (A.6)

with $p$ the order of the RK method ($p+1$ is the order of the embedded one),
$\beta$ a bias factor of the adaptive strategy (typically $\beta$ 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 $\beta\ 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.3. The
correspondence between them and the RK parameters is given in table 1.

MiMeS | Runge-Kutta
---|---
absolute_tolerance | Atol
relative_tolerance | Rtol
b | $\beta$
fac_min | $f_{\rm min}$
fac_max | $f_{\rm max}$
Table 1: Correspondence between the user MiMeS user input and the RK
parameters described in the text.

###  A.2 Explicit embedded RK methods

Explicit methods use only the information of the previous step in order to
compute $\vec{k}_{i}$ from

$\vec{k}_{i}=\vec{f}\Bigg{(}\vec{y}_{n}+h\Big{(}\sum_{j=1}^{i-1}a_{ij}\vec{k}_{j}\Big{)},t_{n}+c_{i}\;h\Bigg{)}\;,\quad\forall
i=1,2,\dots,s\;,$ (A.7)

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

$\begin{array}[]{c|cccccc}0&0&0&0&\dots&0&0\\\ c_{2}&a_{21}&0&0&\dots&0&0\\\
c_{3}&a_{31}&a_{32}&0&\dots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\ddots&\vdots\\\
c_{s}&a_{s1}&a_{s2}&a_{s3}&\dots&a_{ss-1}&0\\\ \hline\cr
p&b_{1}&b_{2}&b_{3}&\dots&b_{s-1}&b_{s}\\\
p+1&b_{1}^{\star}&b_{2}^{\star}&b_{3}^{\star}&\dots&b_{s-1}^{\star}&b_{s}^{\star}\end{array}$
(A.8)

It should be noted that $c_{i}=\displaystyle\sum_{j=1}^{s}a_{ij}$.

###  A.3 Rosenbrock methods

Explicit RK methods, encounter instabilities when a system is ‘‘stiff" 171717A
definition of stiffness can be found in [43, 42].; 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 $\vec{k}$; i.e. in
implicit methods. However, then one has to solve a non-linear set of equations
in order to compute $\vec{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 [43]). In these methods, $\vec{k}$ is determined
by

$\left(\hat{I}-\gamma
h\hat{J}\right)\cdot\vec{k}_{i}=h\vec{f}\Big{(}\vec{y}_{n}+\sum_{j=1}^{i-1}a_{ij}\vec{k}_{j},t_{n}+c_{i}h\Big{)}+h^{2}\left(\gamma+\sum_{j=1}^{i-1}\gamma_{ij}\right)\dfrac{\partial\vec{f}}{\partial
t}+h\hat{J}\cdot\sum_{j=1}^{i-1}\gamma_{ij}\vec{k}_{j}\;,$ (A.9)

which is written in such a way that everything is evaluated at $t=t_{n}$. In
eq. (A.9), $\hat{J}=\dfrac{\partial\vec{f}}{\partial\vec{y}}$ the Jacobian of
the system of ODEs, $\hat{I}$ the unit matrix with dimension equal to the
number of ODEs. Moreover, $\gamma$ and $\gamma_{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

⬇

1 /*-----Implementation of the Heun-Euler method-----*/

2 /*-----It shouldn’t be used for MiMeS, since it is likely to fail-----*/

3

4 //LD is the numeric type (\ie double, or long double).

5 template<class LD>

6 //use struct because its variables are public by default.

7 struct HeunEuler{

8 static constexpr unsigned int s=2; // HeunEuler is a 2-stage RK

9 static constexpr unsigned int p=1; // first order, with embedded $p+1=2$

10

11 //these are aliases for the arrays

12 using arr=std::array<LD,s>;

13 using arr2=std::array<std::array<LD,s>,s>;

14

15 static constexpr arr b={0.5,0.5}; // this is $b_{i}$

16 static constexpr arr bstar={0.5,0.5}; // this is $b_{i}^{*}$

17 static constexpr arr c={0,1}; // remember that
$c_{i}=\displaystyle\sum_{j=1}^{s}a_{ij}$

18

19 // this is $a_{ij}$

20 static constexpr arr2 a={

21 arr{0,0},

22 arr{1,0}

23 };

24 };

In order to implement the ROS3w [35] method, one can add the following code in
the header file MiMeS/src/NaBBODES/Rosenbrock/METHOD.hpp

⬇

1 /*-----Implementation of the ROS3w method-----*/

2 /*-----It shouldn’t be used, since it is likely to fail to produce a
result-----*/

3

4

5 //LD is the numeric type (i.e. double, or long double).

6 template<class LD>

7 struct ROS3w{

8 static constexpr unsigned int s=3; // 3-stage

9 static constexpr unsigned int p=2; // second order (embedded order is $p+1$)

10

11 //aliases for the arrays

12 using arr=std::array<LD,s>;

13 using arr2=std::array<std::array<LD,s>,s>;

14

15 static constexpr arr b={0.25,0.25,0.5 }; // this is $b_{i}$

16 // this is $b_{i}^{*}$

17 static constexpr arr bstar={
0.746704703274,0.1144064078371,0.1388888888888};

18

19 static constexpr arr c={0,2/3.,4/3.}; // remember that
$c_{i}=\displaystyle\sum_{j=1}^{s}a_{ij}$

20 static constexpr LD gamma=0.4358665215084; // this is $\gamma$

21

22 // this is $a_{ij}$

23 static constexpr arr2 a={

24 arr{0,0,0},

25 arr{2/3.,0,0},

26 arr{2/3.,2/3.,0}

27 };

28

29 // this is $\gamma_{ij}$

30 static constexpr arr2 g={

31 arr{0,0,0},

32 arr{0.3635068368900,0,0},

33 arr{-0.8996866791992,-0.1537997822626,0}

34 };

35 };

##### 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. 181818One needs to make sure
to use the correct Solver template argument (or SOLVER variable in the
Definitions.mk files), otherwise compilation will fail.

##  Appendix 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
${\rm GeV}$), $h_{\rm eff}{}$, $g_{\rm eff}{}$, with increasing $T$. The
parameters minT and maxT are minimum and maximum interpolation temperatures.
These temperatures are just limits, and the actual interpolation is done
between the closest temperatures in the data file. Moreover, beyond the
interpolation temperatures, both $h_{\rm eff}{}$ and $g_{\rm 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=\dfrac{2\pi^{2}}{45}h_{\rm eff}{}T^{3}$. These quantities are given as the
member functions:

  * •

template<class LD> LD mimes::Cosmo<LD>::heff(LD T): $h_{\rm eff}{}$ as a
function of $T$.

  * •

template<class LD> LD mimes::Cosmo<LD>::geff(LD T): $g_{\rm eff}{}$ as a
function of $T$.

  * •

template<class LD> LD mimes::Cosmo<LD>::dheffdT(LD T): $dh_{\rm eff}{}/dT$ as
a function of $T$.

  * •

template<class LD> LD mimes::Cosmo<LD>::dgeffdT(LD T): $dg_{\rm eff}{}/dT$ as
a function of $T$.

  * •

template<class LD> LD mimes::Cosmo<LD>::dh(LD T):
$\delta_{h}=1+\frac{1}{3}\frac{d\log h_{\rm eff}{}}{d\log T}$ as a function of
$T$.

  * •

template<class LD> LD mimes::Cosmo<LD>::s(LD T): The entropy density of the
plasma as a function of $T$.

  * •

template<class LD> LD mimes::Cosmo<LD>::rhoR(LD T): The energy density of the
plasma as a function of $T$.

  * •

template<class LD> LD mimes::Cosmo<LD>::Hubble(LD T): The Hubble parameter
assuming radiation dominated expansion as a function of $T$.

Moreover, there are several cosmological quantities are given as members
variables:

  * •

template<class LD> constexpr static LD mimes::Cosmo<LD>::T0: CMB temperature
today [45] in ${\rm GeV}$.

  * •

template<class LD> constexpr static LD mimes::Cosmo<LD>::h_hub: Dimensionless
Hubble constant [45].

  * •

template<class LD> constexpr static LD mimes::Cosmo<LD>::rho_crit: Critical
density [45] in ${\rm GeV}^{3}$.

  * •

template<class LD> constexpr static LD mimes::Cosmo<LD>::relicDM_obs: Central
value of the measured DM relic abundance [7].

  * •

template<class LD> constexpr static LD mimes::Cosmo<LD>::mP: Planck mass [45]
in ${\rm GeV}$.

###  B.2 AnharmonicFactor class

The class mimes::AnharmonicFactor<LD> is responsible interpolating the
anharmonic factor as defined in eq. (2.29). 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 $\theta_{\rm peak}{}$ (which
should be in increasing order) and $f(\theta_{\rm 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 $\theta_{\rm
peak}{}=$theta_peak. Although, there is no need to call this function beyond
the interpolation limits (as long as the data file contains $0\leq\theta_{\rm
peak}{}\leq\pi$), 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.4). 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 ${\rm GeV}$) and $\chi$ (in ${\rm 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}\geq$minT and $T_{\max}\leq$maxT. Beyond these limits, by default,
the axion mass is assumed to be constant. However, this can be changed by
using the member functions

⬇

1 void set_ma2_MIN(std::function<LD(LD,LD)> ma2_MIN)

2 void set_ma2_MAX(std::function<LD(LD,LD)> ma2_MAX)

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}$, $\chi(T_{\rm min})$, and
$\chi(T_{\rm 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_{\rm min}$.

  * •

template<class LD> LD mimes::AxionMass<LD>::getTMax(): This function returns
the maximum interpolation temperature, $T_{\rm max}$.

  * •

template<class LD> LD mimes::AxionMass<LD>::getChiMin(): This function returns
$\chi(T_{\rm min})$.

  * •

template<class LD> LD mimes::AxionMass<LD>::getChiMax(): This function returns
$\chi(T_{\rm max})$.

An alternative way to define the axion mass is via the constructor

⬇

1 template<class LD>

2 mimes::AxionMass<LD>(std::function<LD(LD,LD)> ma2)

Here, the only argument is the axion mass squared, $\tilde{m}_{a}{}$, defined
as a callable object.

Once an instance of the class is defined, we can get $\tilde{m}_{a}{}^{2}$
using the member function

⬇

1 template<class LD> LD mimes::AxionMass<LD>::ma2(LD T, LD fa)

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:

⬇

1 template<class LD> void
mimes::AxionMass<LD>::set_ma2(std::function<LD(LD,LD)> ma2)

###  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.19, which is passed to the ODE solver of
NaBBODES.

The constructor of the class is

⬇

1 template<class LD>

2 mimes::AxionEOM<LD>(LD fa, LD ratio_ini, std::string inputFile,
AxionMass<LD> *axionMass)

The role of the arguments are discussed in section (3.1.1) and (4.3.2) as well
as in table 2. One the instance is created, the interpolations are constructed
by calling the member function

⬇

1 template<class LD> void mimes::AxionEOM<LD>::makeInt()

Then, the temperature as a function of $u=\log a/a_{\rm ini}{}$, is given via
the member function

⬇

1 template<class LD> LD mimes::AxionEOM<LD>::Temperature(LD u)

Another useful member function is

⬇

1 template<class LD> LD mimes::AxionEOM<LD>::logH2(LD u)

This function returns $\log H^{2}$ as a function of $u$. Moreover, its
derivative, $\frac{d\log H^{2}}{du}$, is computed using

⬇

1 template<class LD> LD mimes::AxionEOM<LD>::dlogH2du(LD u)

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

⬇

1 template<class LD>

2 void mimes::AxionMass<LD>::operator()(std::array<LD,2> &lhs,
std::array<LD,2> &y, LD u)

Here, the inputs are $u=\log a/a_{\rm ini}{}$, y[0]$=\theta$, and
y[1]$=\zeta$; which are used to calculate the components of the EOM, with
lhs[0]$=\frac{d\theta}{du}$ and lhs[1]$=\frac{d\zeta}{du}$.

###  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.19. Its header file is
MiMeS/src/Axion/AxionSolve.hpp, and its constructor is

⬇

1 template<class LD, const int Solver, class Method>

2 mimes::Axion< LD, Solver, Method >(LD theta_i, LD fa, LD umax, LD TSTOP,

3 LD ratio_ini, unsigned int N_convergence_max, LD convergence_lim,

4 std::string inputFile, AxionMass<LD> *axionMass, LD initial_step_size=1e-2,

5 LD minimum_step_size=1e-8, LD maximum_step_size=1e-2, LD
absolute_tolerance=1e-8,

6 LD relative_tolerance=1e-8, LD beta=0.9, LD fac_max=1.2, LD fac_min=0.8,

7 unsigned int maximum_No_steps=10000000)

The various arguments are discussed in section (3.1.1) and (4.3.2); and
outlined in table 2.

The member function responsible for solving the EOM is

⬇

1 template<class LD, const int Solver, class Method>

2 void mimes::Axion< LD, Solver, Method>::solveAxion()

Once this function finishes, the results are stored in several member
variables.

The quantities $a/a_{\rm ini}{},\ T,\ \theta,\ \zeta,\rho_{a}$, at the
integration steps are stored in

⬇

1 template<class LD, const int Solver, class Method>

2 std::vector< std::vector<LD> > mimes::Axion< LD, Solver, Method>::points

The quantities $a/a_{\rm ini}{},\ T,\ \theta,\ \zeta,\rho_{a},\ J$, at the
peaks of the oscillation are stored in

⬇

1 template<class LD, const int Solver, class Method>

2 std::vector< std::vector<LD> > mimes::Axion< LD, Solver, Method>::peaks

Note that these points are computed using linear interpolation between two
integration points with a change in the sign of $\zeta$.

The local integration errors for $\theta$ and $\zeta$ are stored in

⬇

1 template<class LD, const int Solver, class Method>

2 std::vector<LD> mimes::Axion< LD, Solver, Method>::dtheta

3

4 template<class LD, const int Solver, class Method>

5 std::vector<LD> mimes::Axion< LD, Solver, Method>::dzeta

Moreover, the oscillation temperature, $T_{\rm osc}{}$, and the corresponding
values of $a/a_{\rm ini}{}$ and $\theta$ are given in

⬇

1 template<class LD, const int Solver, class Method>

2 LD mimes::Axion< LD, Solver, Method>::T_osc

3

4 template<class LD, const int Solver, class Method>

5 LD mimes::Axion< LD, Solver, Method>::a_osc

6

7 template<class LD, const int Solver, class Method>

8 LD mimes::Axion< LD, Solver, Method>::theta_osc

Also, the entropy injection between the last peak ($T=T_{\rm peak}$) and today
($T=T_{0}$), $\gamma$ (defined as in eq. (2.32)), is given in

⬇

1 template<class LD, const int Solver, class Method>

2 LD mimes::Axion< LD, Solver, Method>::gamma

The relic abundance is stored in the following member variable

⬇

1 template<class LD, const int Solver, class Method>

2 LD mimes::Axion< LD, Solver, Method>::relic

We can set another initial condition, $\theta_{\rm ini}{}$, using

⬇

1 template<class LD, const int Solver, class Method>

2 void mimes::Axion< LD, Solver, Method>::setTheta_i(LD theta_i)

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

⬇

1 template<class LD, const int Solver, class Method>

2 void mimes::Axion< LD, Solver, Method>::restart()

Again, this function clears all member variables. So it should be used with
caution.

Finally, there is static mimes::Cosmo<LD> member variable

⬇

1 template<class LD, const int Solver, class Method>

2 static mimes::Cosmo<LD> mimes::Axion< LD, Solver, Method>::plasma

This variable can be used without an instance of the
mimes::Axion<LD,Solver,Method> class.

##  Appendix 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

2 sysPath.append(’path_to_src’)

The available models are Cosmo, AxionMass, and Axion, each defines a class
with the same name.

###  C.1 Cosmo class

The Cosmo module defines the Cosmo class, which contains information about the
plasma. The relevant shared library (lib/libCosmo.so) is obtained by compiling
MiMeS/src/Cosmo/Cosmo.cpp using make lib/libCosmo.so.

The class can be imported by running

⬇

1 from interfacePy.Cosmo import Cosmo

Its constructor is

⬇

1 Cosmo(cosmo_PATH, minT=0, maxT=mP)

The argument cosmo_PATH is the path (a string) of a data file that contains
$T$ (in ${\rm GeV}$), $h_{\rm eff}{}$, $g_{\rm 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_{\rm
eff}{}$ and $g_{\rm 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

⬇

1 del cosmo

The member functions of this class are:

  * •

Cosmo.heff(T): $h_{\rm eff}{}$ as a function of $T$.

  * •

Cosmo.geff(T): $g_{\rm eff}{}$ as a function of $T$.

  * •

Cosmo.dheffdT(T): $dh_{\rm eff}{}/dT$ as a function of $T$.

  * •

Cosmo.dgeffdT(T): $dg_{\rm eff}{}/dT$ as a function of $T$.

  * •

Cosmo.dh(T): $\delta_{h}=1+\frac{1}{3}\frac{d\log h_{\rm eff}{}}{d\log T}$ as
a function of $T$.

  * •

Cosmo.s(T): The entropy density of the plasma as a function of $T$.

  * •

Cosmo.rhoR(T): The energy density of the plasma as a function of $T$.

  * •

Cosmo.Hubble(T): The Hubble parameter assuming radiation dominated expansion
as a function of $T$.

The several cosmological quantities are given as members variables:

  * •

Cosmo.T0: CMB temperature today [45] in ${\rm GeV}$.

  * •

Cosmo.h_hub: Dimensionless Hubble constant [45].

  * •

Cosmo.rho_crit: Critical density [45] in ${\rm GeV}^{3}$.

  * •

Cosmo.relicDM_obs: Central value of the measured DM relic abundance [7].

  * •

Cosmo.mP: Planck mass [45] in ${\rm GeV}$.

Note that these values can be directly imported from the module, without
declaring an instance of the class, as

⬇

1 from interfacePy.Cosmo import T0, h_hub, rho_crit, relicDM_obs, mP

###  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.4).
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 ${\rm GeV}$) and $chi$ (in ${\rm 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}\geq$minT and
$T_{\max}\leq$maxT. Beyond these limits, by default, the axion mass is assumed
to be constant. However, this can be changed by using the member functions

⬇

1 AxionMass.set_ma2_MIN(ma2_MIN)

2 AxionMass.set_ma2_MAX(ma2_MAX)

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}$, $\chi(T_{\rm min})$, and
$\chi(T_{\rm max})$. These values can be obtained using the member functions

  * •

AxionMass.getTMin(): This function returns the minimum interpolation
temperature, $T_{\rm min}$.

  * •

AxionMass.getTMax(): This function returns the maximum interpolation
temperature, $T_{\rm max}$.

  * •

AxionMass.getChiMin(): This function returns $\chi(T_{\rm min})$.

  * •

AxionMass.getChiMax(): This function returns $\chi(T_{\rm max})$.

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
${\rm GeV}$) and $f_{a}$, and returns $m_{a}^{2}(T)$ (in ${\rm 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 returns $\tilde{m}_{a}{}^{2}$ is

⬇

1 AxionMassma2(T,fa)

We should note that another ma2 can be changed using the following member
function:

⬇

1 AxionMass.set_ma2(ma2)

Again, ma2(T,fa) is a function (not any callable object) that takes $T$ (in
${\rm GeV}$) and $f_{a}$, and returns $m_{a}^{2}(T)$ (in ${\rm GeV}$).

###  C.3 Axion class

The class Axion, solves the axion EOM 2.19. 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.

Its constructor is

⬇

1 Axion(theta_i, fa, umax, TSTOP,

2 ratio_ini, unsigned int N_convergence_max, convergence_lim,

3 inputFile, axionMass, initial_step_size=1e-2,

4 minimum_step_size=1e-8, maximum_step_size=1e-2, absolute_tolerance=1e-8,

5 relative_tolerance=1e-8, beta=0.9, fac_max=1.2, fac_min=0.8,

6 unsigned int maximum_No_steps=10000000)

Again, the various arguments are discussed in section (3.1.1) and (4.3.2); and
can also be found in table 2. Notice one important difference between this and
the C++ verison of this class; the instance of AxionMass, axionMass, is passed
by value (there are no pointers in python). However, internally, the
constructor converts this instance to a pointer, which is then passed to the
underlying C function responsible for creating the relevant instance.

The member function responsible for solving the EOM is

⬇

1 Axion.solveAxion()

Once this function is finished, the following member funcions are available

  * •

Axion.T_osc: the oscillation temperature, $T_{\rm osc}{}$, in ${\rm GeV}$.

  * •

Axion.a_osc: $\dfrac{a}{a_{\rm ini}{}}$ at the oscillation temperature.

  * •

Axion.theta_osc: $\theta_{\rm osc}{}$, i.e. $\theta$ at $T_{\rm osc}{}$.

  * •

Axion.gamma: the entropy injection between the last peak ($T=T_{\rm peak}$)
and today ($T=T_{0}$), $\gamma$ (defined as in eq. (2.32)).

  * •

Axion.relic The relic abundance of the axion.

The evolution of $a/a_{\rm ini}{},\ T,\ \theta,\ \zeta,\rho_{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, $\dfrac{a}{a_{\rm ini}{}}$.

  * •

Axion.T: The temperature in ${\rm GeV}$.

  * •

Axion.theta: The axion angle, $\theta$.

  * •

Axion.zeta: The derivative of $\theta$,
$\zeta\equiv\dfrac{d\theta}{d\log(a/a_{\rm ini}{})}$.

  * •

Axion.rho_axion: The axion energy density in ${\rm GeV}^{4}$.

Moreover, the function

⬇

1 Axion.getPeaks()

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_{\rm ini}{},\ T,\ \theta,\ \zeta,\rho_{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 $\zeta$.

The local integration errors for $\theta$ and $\zeta$ are stored in
Axion.dtheta and Axion.dzeta, after the following function is run

⬇

1 Axion.getErrors()

Another initial condition, $\theta_{\rm 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.

##  Appendix 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.14)).

Moreover, there is a function that helps to determine $T_{\rm osc}{}$ and
$\gamma$, 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 $\gamma$ (the
entropy injection between Tosc and today) and $T_{\rm osc}{}$.

###  D.2 The ScanScript module

The python interface of MiMeS has two simple classes that help to scan over
$\theta_{\rm 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 $\theta_{\rm ini}{}$, $f_{a}{}$ (${\rm GeV}$), $\theta_{\rm
osc}{}$, $T_{\rm osc}{}$ (in ${\rm GeV}$), and $\Omega h^{2}$, for every
combination of $\theta_{\rm ini}{}$ and $f_{a}{}$ that are passed as input.

This class is imported using

⬇

1 from interfacePy.ScanScript import Scan

Its constructor is

⬇

1 ScanScript.Scan(cpus,table_fa,table_theta_i,umax,TSTOP,ratio_ini,

2 N_convergence_max,convergence_lim,inputFile,

3 PathToCppExecutable, break_after=5*60,break_time=5,break_command=’’,

4 initial_step_size=1e-2, minimum_step_size=1e-8, maximum_step_size=1e-2,

5 absolute_tolerance=1e-8, relative_tolerance=1e-8,

6 beta=0.9, fac_max=1.2, fac_min=0.8, maximum_No_steps=int(1e7))

Some arguments are the same as in the Axion class, and their definition can be
found in table 2. The other arguments are

  * •

cpus: number of points to run simultaneously (number of cpus available).

  * •

table_fa: table of $f_{a}{}$ to scan.

  * •

table_theta_i: table of $\theta_{\rm ini}{}$ to scan

  * •

PathToCppExecutable: path to an executable that takes theta_i, fa, umax,
TSTOP, ratio_ini, N_convergence_max, convergence_lim, and inputFile; and
prints theta_i, fa, theta_osc, T_osc, and relic.

  * •

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 $\theta_{\rm ini}{}$ closer to the observed DM relic
abundance. The result file is a time-coded file with columns correspond that
to $\theta_{\rm ini}{}$, $f_{a}{}$ (in ${\rm GeV}$), $\theta_{\rm osc}{}$,
$T_{\rm osc}{}$ (in ${\rm GeV}$), and $\Omega h^{2}$. It is imported using

⬇

1 from interfacePy.ScanScript import ScanObs

Its constructor is

⬇

1 ScanScript.ScanObs(cpus,table_fa,len_theta,umax,TSTOP,ratio_ini,

2 N_convergence_max,convergence_lim,inputFile,axionMass,

3 PathToCppExecutable, relic_obs,relic_err_up,relic_err_low,

4 cosmo=Cosmo(_PATH_+r’/src/data/eos2020.dat’,0,Cosmo.mP),

5 break_after=5*60,break_time=5,break_command=’’,

6 initial_step_size=1e-2, minimum_step_size=1e-8, maximum_step_size=1e-2,

7 absolute_tolerance=1e-8, relative_tolerance=1e-8,

8 beta=0.9, fac_max=1.2, fac_min=0.8, maximum_No_steps=int(1e7)

Here, len_theta is the number of different values of $\theta$ 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 $\Omega 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 $\theta_{\rm ini}{}\ll 1$, and finds the value of $\theta_{\rm ini}{}$
that would result in $\Omega h^{2}=$relic_obs; $\theta_{\rm
ini}{}^{\rm(approx)}$. This is easy to do, since for $\theta_{\rm ini}{}\ll
1$, the EOM becomes independent of $\theta_{\rm ini}{}$.

Then, it creates an array, table_theta, of size len_theta that contains values
of $\theta_{\rm ini}{}$ between $Min(0.85\ \theta_{\rm
ini}{}^{\rm(approx)},1)$ and $Max(1.2\ \theta_{\rm ini}{}^{\rm(approx)},\pi)$.
Then it calls ScanObs.run_batch(), which works as Scan.run_batch().

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 [46]. The class can be
imported using

⬇

1 from interfacePy.FT import FT

The constructor of the class is

⬇

1 FT(_M_xticks,_M_yticks, _M_xticks_exception,_M_yticks_exception,
_m_xticks,_m_yticks,

2 xmin,xmax,ymin,ymax,xscale,yscale)

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.

  * •

xmin, xmax (ymin,ymax): Minimum and maximum of the x (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

⬇

1 FT.format_ticks(plt,sub)

It is important to not that plt is the matplotlib.pyplot module that is
usually imported as

⬇

1 import matplotlib.pyplot as plt

#### Example of FT

Consider as an example the plot of $f(x)=e^{-\frac{1}{2}x^{2}}$ for $-5\leq
x\leq 5$. In order to do this, we need to run

⬇

1 import numpy as np#you usually need numpy

2

3 #---these are for plots---#

4 import matplotlib.pyplot as plt

5

6 #load the FT module

7 from sys import path as sysPath

8 sysPath.append(’../../src’)

9 from interfacePy.FT import FT

10

11 fig=plt.figure(figsize=(9,4))

12 fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.95,
right=0.9,wspace=0.0,hspace=0.0)

13 sub = fig.add_subplot(1,1,1)

14

15 x=np.linspace(-5,5,500)

16 sub.plot(x,np.exp(-1/2. * x**2))

Now we can format the ticks using FT. For the shake of this example, for the
x-axis, we choose as major ticks all integrs except $3$ and $-3$, and minor
ticks the halves of each major tick. Also, for the y-axis, we choose to only
show major ticks for every $0.1$ in the interval of interest. Finally, we set
the plot limits as $x\in[-5,5]$ ad $y\in[0,1]$, usng linear scale for both.
This can be done as

⬇

1 #set major ticks

2 _M_xticks=[i for i in range(-5,6)]

3 _M_yticks=[i/10. for i in range(0,11)]

4

5 #set major ticks that will not have a label

6 _M_xticks_exception=[-3,3]

7 _M_yticks_exception=[]

8

9 _m_xticks=[i+j/2. for i in range(-5,6) for j in [0,1]]

10 _m_yticks=[]

11 ft=FT(_M_xticks, _M_yticks, _M_xticks_exception, _M_yticks_exception,
_m_xticks, _m_yticks,

12 xmin=-5, xmax=5, ymin=0, ymax=1, xscale=’linear’, yscale=’linear’)

13

14 ft.format_ticks(plt,sub)

Figure 6: Example plot obtained by using the FT class.

The resulting plot is shown in Fig. (6).

##  Appendix 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

⬇

1 bash MiMeS/src/util/FormatFile.sh path_to_dat_file

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>
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, $\chi(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. 191919For 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:
3
foo

The script is called as

⬇

1 bash MiMeS/src/util/timeit.sh path_to_executable path_to_dat_file

This script is useflul for running the example code in
MiMeS/UserSpace/Cpp/Axion/Axion.cpp, which is compiled by calling ‘‘make
examples" in the root directory of MiMeS. The executable that results from the
compilation of this source, MiMeS/UserSpace/Cpp/Axion/Axion.run, expects the
same arguments as the constructor of the mimes::Axion<LD,Solver,Method> class
(see section B). This means that we can write a file (call it INPUT) that
contains all these arguments in each line, and simply run

⬇

1 bash MiMeS/src/util/timeit.sh MiMeS/UserSpace/Cpp/Axion/Axion.run INPUT

The files KinationInputs, MatterInputs, and RDInputs; in
MiMeS/UserSpace/Cpp/Axion, show these files should look like in different
cases.

###  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-time -- i.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 $\sum_{i=0}^{10^{4}}\ i$ can be determined as

⬇

1 #include "src/util/timeit.hpp"

2

3 int main(){

4 auto sum = 0;

5 {

6 Timer _timer_;

7 for(auto i=0; i<=10000; ++i){ ++sum;}

8 }

9

10 return 0;

11 }

###  E.4 linspace function

The function linspace, defined in MiMeS/src/util/linspace.hpp, is used to
generate linearly spaced numbers between some boundary. This function is
overloaded; there are two functions with the same name and different
signature. The first one is

⬇

1 /*linspace where the list is pushed in a vector that is pass by reference*/

2 template<class LD>

3 void mimes::util::linspace(LD min, LD max, unsigned int length,
std::vector<LD> &X)

Here, X is cleared, and gets filled with length linearly spaced numbers
between min and max.

The overloaded version of this function is

⬇

1 /*linspace where the list returned as a vector*/

2 template<class LD>

3 std::vector<LD> mimes::util::linspace(LD min, LD max, unsigned int length)

This function does not get a reference to a std::vector<LD>, instead it return
one.

###  E.5 logspace function

The function logspace, defined in MiMeS/src/util/logspace.hpp, is used to
generate $\log_{10}$-spaced points between a given boundary. The signatures of
this function are

⬇

1 /*logspace where the list is pushed in a vector that is pass by reference*/

2 template<class LD>

3 void mimes::util::logspace(LD min, LD max, unsigned int length,
std::vector<LD> &X)

4

5 /*logspace where the list returned as a vector*/

6 template<class LD>

7 std::vector<LD> mimes::util::logspace(LD min, LD max, unsigned int length)

Both of these functions work as linspace.

###  E.6 map function

The function map, defined in MiMeS/src/util/map.hpp, applys a function of an
std::vector, and fills another std::vector with the result. The signatures of
this function are

⬇

1 /*The list is pushed in a vector that is pass by reference*/

2 template<class LD>

3 void map(const std::vector<LD> &X, std::function<LD(LD)> func,
std::vector<LD> &Y )

4

5 /*The list returned as a vector*/

6 template<class LD>

7 std::vector<LD> map(const std::vector<LD> &X, std::function<LD(LD)> func)

Both of these functions evaluate func (a callable object) for every element of
X. If another std::vector is given, then it is cleared and filled with the
results, otherwise the function returns a std::vector with the values of
func$($X$)$.

##  Appendix 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.19.
---
theta_i | Initial angle.
fa | The PQ scale.
umax |  Once $u=\log a/a_{i}>$umax, the integration stops. Typical value: $\sim 500$.
TSTOP |  Once $T<$TSTOP, integration stops. Typical value: $10^{-4}~{}{\rm GeV}$.
ratio_ini |  Integration starts at $u$ with $3H/\tilde{m}_{a}{}\approx$ratio_ini. Typical value: $\sim 10^{3}$.
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~{}[{\rm 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.
initial_stepsize |  Initial step-size of the solver. Default value: $10^{-2}$.
minimum_stepsize |  Lower limit of the step-size. Default value: $10^{-8}$.
maximum_stepsize |  Upper limit of the step-size. Default value: $10^{-2}$.
absolute_tolerance | Absolute tolerance of the RK solver (see also table 1). Default value: $10^{-8}$.
relative_tolerance | Relative tolerance of the RK solver (see also table 1). Default value: $10^{-8}$.
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 table 1). Default values: $1.2$ and $0.8$, respectively.
maximum_No_steps | If integration needs more than maximum_No_steps integration stops. Default value: $10^{7}$.
Table 2: Table of the constructor arguments of the
mimes::AxionSolve<LD,Solver,Method> class. 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 ${\rm GeV}$), $\chi(T)$ (in ${\rm 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}\geq$minT and $T_{\max}\leq$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 returns $\tilde{m}_{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 ${\rm GeV}$), $h_{\rm eff}{}$, $g_{\rm eff}{}$. If the path changes one must run bash configure.sh and make.
axMDat | Relative path to data file with $T$ (in ${\rm GeV}$), $\chi$ (defined from eq. (2.2)). 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 $\theta_{\rm peak}{}$, $f(\theta_{\rm peak}{})$. If the path changes one must run bash configure.sh and make.
Table 4: Paths to the required data files. Variables defined in the
MiMeS/Paths.mk files, and used when running bash configure.sh in the root
directory of MiMeS. Template arguments.
---
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 $\theta$.
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.
Table 5: Template arguments of the various MiMeS classes. User compile-time
options. Variables in the various Definitions.mk files.
---
rootDir | The relative path of root directory of MiMeS. Relevant only when compiling using make. Available in all Definitions.mk.
LONG | long for long double or empty for double. This is defines a macro in the source files of the various C++ examples. Available in Definitions.mk inside the various subdirectories of MiMeS/UserSpace/Cpp.
LONGpy | long or empty. Same as LONG, applies in the python modules. Available in MiMeS/Definitions.mk.
SOLVER | In order to use a Rosenbrock method SOLVER=$1$. For explicit RK method, SOLVER=$2$. This defines a macro that is passes as the second template argument of mimes::Axion<LD,Solver,Method>. The corresponding variable in MiMeS/Definitions.mk applies to the python modules. The variable in MiMeS/UserSpace/Cpp/Axion/Definitions.mk applies to the example in the same directory.