diff --git a/Documentation/Tex/MiMeS.tex b/Documentation/Tex/MiMeS.tex
index 1a0e080..f77bf1d 100644
--- a/Documentation/Tex/MiMeS.tex
+++ b/Documentation/Tex/MiMeS.tex
@@ -1,1443 +1,1579 @@
\documentclass[11pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage[left=2cm,right=2cm,top=2cm,bottom=2cm]{geometry}
\usepackage{slashed}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{float}
\usepackage{subcaption}
\usepackage{multirow}
\usepackage{minted}
\usepackage[dvipsnames]{xcolor}
\usepackage{multicol}
\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}
\include{macros}
%%%%%%%-----------------------Rules-----------------------%%%%%%%
% 1. Put any new macros in macros.tex.
% 2. Section should start with:
%\section{This is a section}\label{sec:Intro}
%\setcounter{equation}{0}
% 3. labels for Figs should start with fig:, for equations should start with eq:, for sections with sec:, etc.
%%%%%%%--------------------------------------------------%%%%%%%
\author{Karamitros Dimitrios}
\title{{\tt MiMeS}: Misalignment Mechanism Solver}
\begin{document}
\maketitle
\begin{abstract}
We introduce a \CPP header-only library that solves the Axion equation of motion, \mimes.
\mimes makes no assumptions regarding the cosmology and the thermal mass of the axion, which allows the user
to consider various cosmological scenarios and axion-like models.
\mimes also includes a convenient \PY interface, which allows the user to use it without writing any code \CPP.
\end{abstract}
\section{Introduction}\label{sec:intro}
\setcounter{equation}{0}
About the axion...
\section{Physics background}\label{sec:abundance}
\setcounter{equation}{0}
%
Although there are several works in the literature (such as~\cite{Chang:1998ys}) that can provide a insight on the solution of the axion equation of motion (EOM), in this section we define, derive, and discuss the various quantities we need, in order to understand how \mimes works in detail.
\paragraph{The EOM} %The axion field is written as
%
%\begin{equation}
% A \equiv \fa \ \theta,
% \label{eq:fa_def}
%\end{equation}
%%
%with $\fa$ the scale of the axion that determines the the scale at which the PQ symmetry breaks.
%
The EOM in terms of the axion angle, $\theta$, is
%
\begin{equation}
\lrb{\dfrac{d^2}{d t^2} + 3 H(t) \ \dfrac{d}{d t} } \theta(t) + \maT^2(t) \ \sin \theta(t) = 0 \; ,
\label{eq:eom}
\end{equation}
%
with $H(t)$ the Hubble parameter (determined by the cosmology), and $\maT(t)$ the time (temperature) dependent mass of the axion, with
\begin{equation}
\maT^2(T) = \dfrac{\chi(T)}{\fa^2}\;,
\label{eq:axion_mass_def}
\end{equation}
%
and $\chi$ a function of the temperature, and $\fa$ the scale at which the Peccei Quinn symmetry breaks \DK{ref??}. For the QCD axion, this has been calculated using lattice simulations in~\cite{Borsanyi:2016ksw}.~\footnote{The data provided by ref.~\cite{Borsanyi:2016ksw} are used by \mimes out-of-the-box. However, the user is free to change them.}
\paragraph{Initial conditions}
%
\begin{figure}[h!]
\includegraphics[width=1\textwidth]{figs/axion_mass.pdf}
\caption{The mass of the axion as a function of the temperature for $\fa=10^{12}~\GeV$, using the data provided in ref.~\cite{Borsanyi:2016ksw}.}
\label{fig:axion_mass}
\end{figure}
%
Assuming that the PQ symmetry breaks before inflation, the initial conditions (\ie at some $t=0$, after inflation) for the EOM is random. However, we note that $\maT \to 0$ (see \Figs{fig:axion_mass}) -- \ie $\maT \ll H$ -- at very early times. Therefore, after inflation, the EOM is simply
%
\begin{equation}
\lrb{\dfrac{d^2}{d t^2} + 3 H(t) \ \dfrac{d}{d t} } \theta(t) = 0 \; ,
\label{eq:massless_eom}
\end{equation}
%
which is solved by $\theta = \thetai + C \dint_{0}^t d t' \ \lrb{ \dfrac{a(t'=0)}{a(t')} }^3$, with $\thetai$ a constant and $a$ the scale factor of the Universe. That is, as the Universe expands, $\theta \approx \thetai$. Since we would like to calculate the relic abundance of axions, we can integrate \eqs{eq:eom} from a time after inflation (call it $t = \ti$) such that $ \dot \theta |_{t=\ti} = 0$ and $\theta|_{t=\ti} \approx \thetai$.
\subsection{The WKB approximation}
%
In order to solve analytically \eqs{eq:eom}, we assume $\theta \ll 1$, which results in the linearised EOM
%
\begin{equation}
\lrb{\dfrac{d^2}{d t^2} + 3 H(t) \ \dfrac{d}{d t} + \maT^2(t) } \theta(t) = 0 \; .
\label{eq:linear_eom}
\end{equation}
Using a trial solution $\theta_{\rm trial} = \exp\lrsb{ i \dint d t \ \lrBigb{u(t) +3/2 \ i \ H(t)} }$, and defining $\Omega^2 = \maT^2 - \dfrac{9}{4} H^2 - \dfrac{3}{2} \dot H $ we can transform the \eqs{eq:linear_eom} to
%
\begin{equation}
u^2 = \Omega^2 + i \ \dot u \; ,
\label{eq:eom_of_u}
\end{equation}
%
which has a formal solution $u = \pm \sqrt{\Omega^2 + i \dot u}$. Assuming that $\dot u \ll \Omega^2$ and $\dot \Omega \ll \Omega^2$, we can approximate $u$ as
%
\begin{equation}
u \approx \pm \Omega + \dfrac{i}{2} \dfrac{d \log \Omega}{d t} \;,
\label{eq:u_approx}
\end{equation}
%
which results in the general solution of \eqs{eq:linear_eom}
%
\begin{equation}
\theta \approx \dfrac{1}{\sqrt{\Omega}} \exp\lrb{-\dfrac{3}{2} \int d t \ H} \lrsb{ A \cos\lrb{ \int d t \ \Omega} + B \sin\lrb{ \int d t \ \Omega} } \;.
\label{eq:general_solution_eom_approx}
\end{equation}
Applying, then, the initial conditions $ \dot \theta |_{t=\ti} = 0$ and $\theta|_{t=\ti} \approx \thetai$, we arrive at
%
\begin{equation}
\theta(t) \approx \thetai \sqrt{ \dfrac{ \Omegai }{\Omega (t)} } \lrb{\dfrac{a}{\ai}}^{-3/2} \ \cos\lrb{ \int_{\ti}^t d t^\prime \ \Omega(t^\prime)} \;.
\label{eq:solution_eom_approx}
\end{equation}
In order to further simplify this approximate result, we note that $\theta$ deviates from $\thetai$ close to $t=\tosc$ -- corresponding to $T = \Tosc$, the so-called ``oscillation temperature" -- $\maT|_{t = \tosc} = 3 H|_{t = \tosc}$, which is defined as the point at which the axion begins to oscillate.
%
This observation allows us to set $\ti = \tosc$. Moreover, at $t > \tosc$, we approximate $\Omega \approx \maT$, as $H^2$ and $\dot H$ become much smaller than $\maT^2$ quickly after $t=\tosc$. Finally, the axion angle takes the form
%
\begin{equation}
\theta(t) \approx \thetaosc \lrb{\dfrac{3}{4}}^{1/4} \sqrt{ \dfrac{ \maT|_{t=\tosc} }{\maT (t)} } \lrb{\dfrac{a}{\aosc}}^{-3/2} \ \cos\lrb{ \int_{\tosc}^t d t^\prime \ \maT(t^\prime)} \;,
\label{eq:solution_eom_approx_theta_osc}
\end{equation}
%
where $\thetaosc = \theta|_{t=\tosc}$. This equation is further simplified if we assume that $\thetaosc \approx \thetai$, \ie
%
\begin{equation}
\theta(t) \approx \thetai \lrb{\dfrac{3}{4}}^{1/4} \sqrt{ \dfrac{ \maT|_{t=\tosc} }{\maT (t)} } \lrb{\dfrac{a}{\aosc}}^{-3/2} \ \cos\lrb{ \int_{\tosc}^t d t^\prime \ \maT(t^\prime)} \;.
\label{eq:solution_eom_approx_final}
\end{equation}
%
It is worth mentioning that the accuracy of this approximation depends, in general, on $\Tosc$; it determines the difference between $\thetai$ and $\thetaosc$, the deviation of $\dot \theta|_{t=\tosc}$ from $0$, and whether $\dot \Omega \ll \Omega^2$.
\paragraph{Axion energy density}
%
The energy density of the axion is
%
\begin{eqnarray}
\rho_{a} = \dfrac{1}{2} \fa^2 \lrsb{ \dot{\theta}^2 + \maT^2 \theta^2 } \;.
\label{eq:rho_a_def}
\end{eqnarray}
%
For the relic abundance of axions, we need to calculate their energy density at very late times. That is, $\dot{\tilde{m}}_a = 0$, $\maT \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 $a$)
%
\begin{eqnarray}
\rho_{a} \approx \dfrac{\ma }{2} \ \fa^2 \ \thetai^2 \ \maT(\aosc) \ \lrb{\dfrac{\aosc}{a}}^3 \;,
\label{eq:rho_a0}
\end{eqnarray}
%
which shows that the energy density of axions at late times scales as the energy density of matter. If there is a period of entropy injection to the plasma for $T<\Tosc$, the axion energy density gets diluted, since
%
\begin{equation}
a^3 \ s = \gamma \ \aosc^3 \ s_{\rm osc} \Rightarrow \lrb{\dfrac{\aosc}{a}}^3 = \gamma^{-1} \dfrac{s}{s_{\rm osc}} \;,
\end{equation}
%
with $\gamma$ the amount of entropy injection to the plasma between $\tosc$ and $t$. Therefore, the present (at $T=T_0$) energy density of the axion, becomes
%
\begin{eqnarray}
\rho_{a,0} = \gamma^{-1} \dfrac{s_0}{s_{\rm osc}} \ \dfrac{1 }{2} \ \fa^2 \ \ma \ \maT_{,{\rm osc}} \ \thetai^2 \;,
\label{eq:rho_a_approx}
\end{eqnarray}
with $\ma$ the mass of the axion at $T=T_0$. Notice that the explicit dependence on $\fa$ cancels, since $\maT \sim 1/\fa$. That is, $\fa$ only affects the energy density of the axions through its impact on $\Tosc$.
%
\subsection{Notation}\label{sec:notation}
%
%The WKB approximation is very useful in order to understand the evolution of the axion field. However, it fails to explain how the oscillation begins before $\dot \Omega \ll \Omega^2$ is reached. In this section we will try to understand the evolution of the axion as generally as possible.
%
The EOM~(\ref{eq:eom}) depends on time, which is not very useful variable in a non-standard comsological setting. Therefore, we introduce
%
\begin{eqnarray}
u = \log \dfrac{a}{\ai} \;,
\label{eq:natation}
\end{eqnarray}
%
which results in
%
\begin{eqnarray}
&\dfrac{d F}{dt} &= H \dfrac{d F}{du}
\nonumber \\
&\dfrac{d^2 F}{dt^2} &= H^2 \ \lrb{ \dfrac{d^2 F}{du^2} + \dfrac{1}{2} \dfrac{d \log H^2}{du} \dfrac{d F}{du} }\;.
\label{eq:deriv_u}
\end{eqnarray}
The EOM in terms of $u$, then, becomes
%
\begin{equation}
\dfrac{d^2 \theta}{du^2} + \lrsb{\dfrac{1}{2} \dfrac{d \log H^2}{du} + 3 } \dfrac{d \theta}{d u} + \ \lrb{\dfrac{\maT}{H}}^2 \ \sin \theta
=0 \;.
\label{eq:eom_u}
\end{equation}
Notice that in a radiation dominated Universe
%
$$
\dfrac{d \log H^2}{du} = -\lrb{ \dfrac{d \log \geff}{d \log T} +4 } \delta_h^{-1}\;,
$$
with $ \delta_h = 1+ \dfrac{1}{3} \dfrac{d \log \heff}{d \log T} $.
%
In a general cosmological setting, the expansion rate is dominated by an energy density that scales as $\rho \sim a^{-c}$, which results in $\dfrac{d \log H^2}{du} = -c$. However, 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.
It is also worth mentioning that the EOM~\ref{eq:eom_u} with the initial condition $\theta(u=0)=\thetai$ and $d\theta/du (u=0)=0$ can be written as a system of first order
%
\begin{eqnarray}
& \dfrac{d \zeta}{du} + \lrsb{\dfrac{1}{2} \dfrac{d \log H^2}{du} + 3 } \zeta + \ \lrb{\dfrac{\maT}{H}}^2 \ \sin \theta
=0 \;. \\ \nonumber
&\dfrac{d \theta}{d u} - \zeta=0 \;.
\label{eq:eom_sys}
\end{eqnarray}
%
This form of the EOM is suitable for integration via Runge-Kutta methods, which are briefly discussed in Appendix~\ref{app:RK}.
\subsection{Beyond the WKB approximation}\label{sec:beyond_WKB}
%
The WKB approximation is very useful, as it helps us understand the evolution of the axion field after it begins to oscillate adiabatically. However, it fails to capture the dynamics the adiabatic conditions are met, result in inaccurate axion relic abundance result. In this section, we examine the deviation of $\thetaosc$ from the initial value of $\theta$, as well how to treat cases where the EOM cannot be linearized.
\subsubsection{Behaviour close to the initial condition}\label{sec:close_to_ini}
%
The accuracy of the approximate evolution of the axion angle given in \eqs{eq:solution_eom_approx_final} depends on the deviation between $\thetai$ and $\thetaosc$. In order to examine their deviation, we expand \eqs{eq:eom} at a time $t =\ti + \delta t$ with $\delta t \to 0$. That is, we use the following approximations
%
\begin{equation*}
\ddot{\theta} \approx \dfrac{\dot \theta (\ti+ \delta t) - \dot \theta (\ti)}{ \delta t } = \dfrac{\dot \theta (\ti + \delta t) }{\delta t } \;,
\qquad
\dot \theta(\ti + \delta t) = \dfrac{\theta(\ti + \delta t) - \thetai}{\delta t} \;.
\end{equation*}
%
and, solve the EOM~\ref{eq:eom} for $\theta(\ti + \delta t)$. This results in
%
\begin{equation}
\theta(\ti + \delta t) \approx \thetai -\dfrac{\delta t^2}{1+3\ \delta t \ H(\ti + \delta t) } \ \maT^2(\ti + \delta t) \sin \theta(\ti + \delta t)
\approx \thetai - \delta t^2 \ \maT^2(\ti) \ \sin \thetai \ + \mathcal{O}(\delta t^3)\;,
\label{eq:theta_dt}
\end{equation}
%
which indicates that the angle decreases (increases) as the temperature drops if $\thetai>0$ ($\thetai<0$). Using the notation introduced in section~\ref{sec:notation}, \eqs{eq:theta_dt} takes the form
%
\begin{eqnarray}
\theta \approx \thetai - \delta u^2 \ \lrb{\dfrac{\maT}{H}}_{t=\ti}^2 \ \sin \thetai \;,
\label{eq:theta_du}
\end{eqnarray}
%
which can be used to estimate the angle at the oscillation temperature
%
\begin{eqnarray}
\thetaosc \approx \thetai - \lrb{\dfrac{\maT}{H}}_{t=\ti}^2
\lrsb{1- \dfrac{\Ti}{\Tosc}\lrb{\dfrac{\heff_{,\rm ini}}{\heff_{,\rm osc}}\gamma_{\rm osc}}^{1/3}}^{2} \ \sin \thetai \;,
\label{eq:theta_osc}
\end{eqnarray}
%
where $\gamma_{\rm osc}$ is the entropy injection between $\Ti$ and $\Tosc$. Notice that in the derivation of \eqs{eq:solution_eom_approx_final} we used $\theta_{\rm osc} = \thetai$ as our first approximation. Thus, \eqs{eq:theta_osc} provides a correction that takes into account the deviation between $\theta_{\rm osc} $ and $ \thetai$, and \eqs{eq:rho_a_approx} becomes (for $\thetai \ll 1$)
%
\begin{eqnarray}
\rho_{a,0} = \gamma^{-1} \dfrac{s_0}{s_{\rm osc}} \ \dfrac{1 }{2} \ \fa^2 \ \ma \ \maT_{,{\rm osc}} \ \thetai^2 \lrBiggcb{
1 - \ \lrb{\dfrac{\maT}{H}}_{t=\ti}^2 \ \lrb{1- \dfrac{\Ti}{\Tosc}\lrb{\dfrac{\heff_{,\rm ini}}{\heff_{,\rm osc}}\gamma_{\rm osc}}^{1/3}}^{2} }^2 \;,
\label{eq:rho_a_NLO}
\end{eqnarray}
%
which implies that the WKB approximation overestimates the energy density of the axion, since $\gamma_{\rm osc}$, $H_{\rm ini}$, and $\Tosc$ can vary.
\subsubsection{Example}\label{sec:comparison_between_WKB_approaches}
In order to demonstrate the effect of the deviation between $\thetai$ and $\thetaosc$, we show \Figs{fig:WKB_diff_i,fig:WKB_diff_osc}, for different cosmological scenarios. In both figures, the red and black (dashed) lines correspond to standard cosmological scenario and early matter dominated (EMD) Universe, respectively.~\footnote{For the matter dominated case (following the notation of refs.~\cite{Arias:2019uol,Arias:2020qty}), we have used
$T_{\rm end}=10^{-2} ~\GeV,\; c=3, \; T_{\rm ini}=10^{12} ~\GeV, \; \text{and}\; r=0.1$.} For the matter dominated case, entropy is injected to the plasma both before and after $\Tosc$. The entropy injection factor, $\gamma_{\rm osc}$, as a function of $\fa$ is shown in \Figs{fig:gamma_osc}.
In \Figs{fig:WKB_diff_i} we show the relative difference between $\Omega h^2$ using \eqs{eq:solution_eom_approx_final} and the result obtained from \mimes. It should be apparent that the accuracy of the estimate based on \eqs{eq:solution_eom_approx_final} depends of $\fa$, since for high values of $\fa$, $\Tosc$ is low. That is, the deviation between $\thetai$ and $\thetaosc$ is greater than the corresponding difference for lower $\fa$, especially for the matter dominated case, where the entropy injection is great (see \Figs{fig:gamma_osc}).
In \Figs{fig:WKB_diff_osc} we show the relative difference between $\Omega h^2$ obtained by integrating numerically \eqs{eq:eom_u} and using \eqs{eq:solution_eom_approx_theta_osc} with \eqs{eq:theta_osc}. The relative difference between the two estimates, once again, depends $\fa$. It should be noted that the estimate for $\thetaosc$ depends on the choice of $\Ti$. For \Figs{fig:WKB_diff_osc}, we define $\Ti$ as the temperature at which $3H/\maT = 2$. We note that this estimate is more accurate, its accuracy still depends on whether the axion evolves adiabatically at $T<\Tosc$. Moreover, an increase of $\gamma_{\rm osc}$ invalidates the estimate of $\thetaosc$; higher powers of $\delta u$ may be needed.
%
\begin{figure}[t]
\begin{subfigure}[]{0.5\textwidth}
\includegraphics[width=1\textwidth]{figs/WKB_diff_i.pdf}
\caption{}
\label{fig:WKB_diff_i}
\end{subfigure}
\begin{subfigure}[]{0.5\textwidth}
\includegraphics[width=1\textwidth]{figs/WKB_diff_osc.pdf}
\caption{}
\label{fig:WKB_diff_osc}
\end{subfigure}
\begin{center}
\begin{subfigure}[]{0.5\textwidth}
\includegraphics[width=1\textwidth]{figs/gamma_osc.pdf}
\caption{}
\label{fig:gamma_osc}
\end{subfigure}
\end{center}
\caption{...}
\label{fig:WKB_diff}
\end{figure}
%
Nevertheless, neither estimates based on \eqs{eq:solution_eom_approx_theta_osc,eq:solution_eom_approx_final} seem to provide a result comparable to numerically integrating \eqs{eq:eom_u}. Therefore, numerical integration should be preferred.
\subsubsection{Adiabatic invariant and the anharmonic factor}\label{sec:an_fac}
%
In oscillatory systems with varying period, the energy is not conserved, and it is usually useful to define an ``adiabatic invariant", which is an approximate constant of motion.
\paragraph{Definition of the adiabatic invariant}
%
Given a system with Hamiltonian $\mathcal{H}(\theta,p;t)$, the equations of motion are
%
\begin{equation}
\dot p = - \dfrac{\partial \mathcal{H}}{\partial \theta} \;, \;\;
\dot \theta = \dfrac{\partial \mathcal{H}}{\partial p} \;.
\label{eq:hamiltonian_eoms}
\end{equation}
Moreover, we note that
%
\begin{equation}
d \Ham = \dot \theta \ d p - \dot p \ d \theta + \dfrac{\partial \Ham}{\partial t} \ d t \;.
\label{eq:total_dH}
\end{equation}
If this system exhibits closed orbits (\eg if it oscillates), we define
%
\begin{equation}
J \equiv C \ \oint p \ d \theta \;,
\label{eq:adiabatic_inv_def}
\end{equation}
%
where the integral is over a closed path (\eg 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{d J}{d t} = C \ \oint \lrBigb{\dot p \ d \theta + p \ d \dot \theta} = C \ \dint_{t}^{t+T} \dfrac{\partial \Ham}{\partial t^\prime} \ d t^\prime \approx T \ \dfrac{\partial \Ham(t^{\prime})}{\partial t^{\prime}}\Big|_{t^{\prime}=t} \approx 0
\;.
\]
%
\paragraph{Application to the axion}
%
The Hamiltonian that results in the EOM of \eqs{eq:eom} is
%
\begin{equation}
\Ham = \dfrac{1}{2} \dfrac{p^2}{\fa^2 \ a^3} + V(\theta) \ a^3\;,
\label{eq:axion_H}
\end{equation}
%
with
%
\begin{eqnarray}
& p = \fa^2 \ a^3 \ \dot \theta \\
\label{eq:momentum}
& V(\theta) = \maT^2 \fa^2 (1-\cos \theta) \;.
\label{eq:potential}
\end{eqnarray}
Notice that the Hamiltonian varies slowly if $\dot \maT/\maT \ll \maT$ and $H \ll \maT$, which are the adiabatic conditions. When these conditions are met, the adiabatic invariant for this system becomes
%
\begin{equation}
J = \dfrac{\oint p \ d \theta}{\pi \fa^2} = \dfrac{1}{\pi \fa^2} \oint \sqrt{ 2\lrb{ \Ham(\theta) - V(\theta) \ a^3} \ \fa^2 a^3 \ } \ d \theta =
\dfrac{2}{\pi \fa^2} \int_{-\thetamax}^{\thetamax} \sqrt{ 2\lrb{ \Ham(\thetamax) - V(\theta) \ a^3} \ \fa^2 a^3 \ } d \theta \;,
\label{eq:J_axion_definition}
\end{equation}
%
where we note that $\thetamax$ denotes the maximum of $\theta$ -- the peak of the oscillation, which corresponds to $p=0$. That is, $\Ham(\thetamax) = V(\thetamax) \ a^3$. Therefore, the adiabatic invariant, takes the form
\begin{eqnarray}
J=& \dfrac{2 \sqrt{2} }{\pi \fa} \int_{- \thetamax} ^{\thetamax} \sqrt{ V(\thetamax) - V(\theta) } a^{3} d \theta =
\dfrac{2 \sqrt{2} }{\pi} \ \maT \, a^3 \ \dint_{- \thetamax}^{\thetamax} \sqrt{\cos \theta - \cos \thetamax} \ d \theta
\;,
\label{eq:J_axion_derivation}
\end{eqnarray}
%
where, for the last equality. we have used the adiabatic conditions, \ie negligible change of $\maT$ and $a$ during one period. Usually, the adiabatic invariant is written as~\cite{Lyth:1991ub,Bae:2008ue}
%
\begin{equation}
J = a^3 \ \maT \ \thetamax^2 \, f(\thetamax) \;,
\label{eq:J_axion_final_form}
\end{equation}
%
where
\begin{equation}
f(\thetamax) =\dfrac{ 2 \sqrt{2}}{\pi \thetamax^2 } \dint_{- \thetamax}^{\thetamax} d \theta \sqrt{ \cos \theta - \cos \thetamax } \;,
\label{eq:anharmonic_f}
\end{equation}
%
is called the anharmonic factor, with $ 0.5 \lesssim f(\thetamax) \leq 1$ (see \Figs{fig:anharmonic_factor}).
\begin{figure}[t]
\includegraphics[width=1\textwidth]{figs/anharmonic_factor.pdf}
\caption{The anharmonic factor for $0 \leq \thetamax < \pi $.}
\label{fig:anharmonic_factor}
\end{figure}
\paragraph{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 after the adiabatic conditions where met. It allows us to calculate the energy density of the axions field for $\theta \gtrsim 1$ since we have taken into account the exact potential.
In order to do this, we can numerically integrate \eqs{eq:eom}, and identify the maxima of $\theta$. Once the adiabatic conditions are fulfilled, we can stop the integration at a peak, $\thetamax_{,*}$ -- which corresponds to $T=T_{*}$ and $a=a_{*}$. Then, the value of the maximum angle today ($\thetamax_{,0} \ll 1$) is related to $\thetamax_{,*}$ via
%
\begin{eqnarray}
\thetamax_{,0}^2 &= \lrb{\dfrac{a_*}{a_0}}^3 \ \dfrac{\maT_{,*}}{\ma} \ f(\thetamax_{,*}) \ \thetamax_{,*}^2 =
\gamma^{-1} \ \dfrac{s_0}{s_*} \ \dfrac{\maT_{,*}}{\ma} \ f(\thetamax_{,*}) \ \thetamax_{,*}^2
\; .
\label{eq:theta_relation}
\end{eqnarray}
%
Using this, and since we evaluate the energy density at the maximum of $\theta$ (\ie $\dot \theta = 0$), we can determine the energy density of axions today from \eqs{eq:rho_a_def}. That is,
%
\begin{equation}
\rho_{a,0} = \gamma^{-1} \ \dfrac{s_0}{s_*} \ \ma \ \maT_{,*} \ \dfrac{1}{2} \ \fa^2 \ \thetamax_{,*}^2 \; \ f(\thetamax_{,*}) \;.
\label{eq:rho_axion_exact}
\end{equation}
%
Notice that this is form of $\rho_{a,0}$ is similar to the the WKB result~(\ref{eq:rho_a_approx}) at $\tosc \to t_*$, multiplied by the anharmonic factor $f(\thetamax_{,*})$. That is, if the axion starts to evolve adiabatically close to $t=\tosc$, and $\thetaosc \ll 1$, the WKB approximation is valid.
-
-
-
-\section{How to start using \mimes}\label{sec:start}
+\section{First steps with \mimes}\label{sec:first_steps}
\setcounter{equation}{0}
%
\begin{listing}
\begin{bash}
git clone git@github.com:dkaramit/MiMeS.git
cd MiMeS
git submodule init
git submodule update --remote
\end{bash}
%
\caption{Commands in order to download \mimes, the differential equation solvers, and the interpolation libraries.}
\label{list:git_download}
\end{listing}
%
The library can downloaded from \href{https:/github.com/dkaramit/MiMeS}{github.com/dkaramit/MiMeS}. Since the library depends on {\tt NaBBODES}~\cite{NaBBODES} {\tt SpleSplines}~\cite{SimpleSplines}, one has to download the, too. However, in order to get \mimes ready, one can run the commands shown in listing~\ref{list:git_download}. Although this method should be preferred in order to get the latest version, \mimes can also be downloaded from \href{https:/mimes.hepforge.org/downloads/}{mimes.hepforge.org/downloads}.
-Once everything is downloaded successfully, we can go inside the \mimes directory, and run \run{bash configure.sh}~\footnote{It is recommended to use {\tt bash}.} and \run{make}. The {\tt bash} script {\tt configure.sh}, just writes some paths to some files, formats the data files provided in a acceptable format (in section~\ref{sec:input} the format is explained), and makes some directories.
+Once everything is downloaded successfully, we can go inside the \mimes directory, and run \run{bash configure.sh}~\footnote{It is recommended to use {\tt bash}.} and \run{make}. The {\tt bash} script {\tt configure.sh}, just writes some paths to some files, formats the data files provided in an acceptable format (in section~\ref{sec:input} the format is explained), and makes some directories.
%
The {\tt makefile} is responsible for compiling some examples and checks, as well as the shared libraries that needed for the \PY interface. If everything runs successfully, there should be two new directories {\tt exec} and {\tt lib}. Inside {\tt exec}, there are several executables that ready to run, in order to ensure that the code runs (\eg no segmentation fault occurs). For example, {\tt exec/AxionSolve\_check.run}, should print the values of the parameters $\thetai$ and $\fa$, the oscillation temperature and the corresponding value of $\theta$, the evolution of the axion (\eg temperature, $\theta$, $\rho_{a}$, etc.), and the values of various quantities on the peaks of the oscillation. In the directory {\tt lib}, there are several shared libraries for the \PY interface.
-\subsection{First steps}\label{sec:first_steps}
+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.
+
+
+\subsection{First examples}\label{sec:First_examples}
%
There are several examples in \CPP ({\tt MiMeS/UserSpace/Cpp}) and \PY ({\tt MiMeS/UserSpace/Python}), as well as \JUPY notebooks ({\tt MiMeS/UserSpace/JupyterNotebooks}), that show in detail how \mimes can be used.
\subsubsection{Using \mimes in \CPP}
%
-In order to use the {\tt mimes::Axion<LD,Solver,Method>} class in a \CPP program, we first need to have access to the {\tt mimes::AxionMass<LD>} class, in order to define the mass of the axion as a function of the temperature and $\fa$. Therefore, we need to include the header file of this class; \ie
+The class that is responsible for the solution of the EOM is {\tt mimes::Axion<LD,Solver,Method>}, located in {\tt MiMeS/src/Axion}. However, in order to use it, we have to define the mass of the axion as a function of the temperature and $\fa$. The axion mass is defined as an instance of the {\tt mimes::AxionMass<LD>}, which is defined in the header file {\tt \mimes/src/AxionMass/AxionMass.hpp}. We should note that {\tt LD} is the numerical type we use (\eg { \tt long double}). The other template arguments are related to the differential equation solver, and their role will be explained in later sections.
+
+Therefore, at the top of the main {\tt .cpp} file, we need to include
%
\begin{cpp}
#include "MiMeS/src/AxionMass/AxionMass.hpp"
+ #include "MiMeS/src/Axion/AxionSolve.hpp"
\end{cpp}
%
-Furthermore, wen eed to include the header file {\tt AxionSolve.hpp} located in {\tt MiMeS/src/Axion}, in order to be able to use the {\tt mimes::Axion<LD,Solver,Method>} class that solves the axion EOM~\ref{eq:eom_sys}. That is, we need to include the following line in our code
+Notice that if the this {\tt .cpp} file is not in the root directory of \mimes, we need to compile it using the flag {\tt -Ipath-to-root}, "path-to-root" the relative path to the root directory of \mimes; \eg if the {\tt .cpp} is in the {\tt MiMeS/UserSpace/Cpp/Axion} directory, this flag should be {\tt -I../../../}.
+
+The mass of the axion can be defined in two ways; either via a data file or by a user defined function.
+
+\paragraph{Axion mass form data file}
+In many cases, axion mass cannot be written in a closed form. In these cases, the user may provide a data file. Then, the axion mass can be defined via this file as
%
\begin{cpp}
- #include "MiMeS/src/Axion/AxionSolve.hpp"
+ mimes::AxionMass<LD> axionMass(chi_PATH, minT, maxT);
\end{cpp}
+%
+Here {\tt chi\_PATH} is a (relative or absolute) path to a file with two columns; $T$ (in $\GeV$) and $\chi$ (in $\GeV^4$), with increasing $T$. In this case, the axion mass is interpolated between the temperatures {\tt minT} and {\tt maxT}. These two parameters are just suggestions, and the actual interpolation limits, {\tt TMin} and {\tt TMax}, are chosen as the closest temperatures in the data file. That is, in general, {\tt TMin} $\geq$ {\tt minT} and {\tt TMax} $\leq$ {\tt maxT}. Beyond these limits, the mass is assumed to constant by default. However, one can change this by using
+%
+\begin{cpp}
+ axionMass.set_ma2_MAX(ma2_MAX);
+ axionMass.set_ma2_MIN(ma2_MIN);
+\end{cpp}
+%
+with {\tt ma2\_MAX} and {\tt ma2\_MIN} the axion mass squared as functions (or any other callable objects) with signatures {\tt LD ma2\_MAX(LD T, LD fa)} and {\tt LD ma2\_MIN(LD T, LD fa)}. These functions are called for $T\geq${\tt TMax} and $T\leq${\tt TMin}, respectively. Usually, in order to ensure continuity of the axion mass, one needs to know {\tt TMin}, {\tt TMax}, $\chi(${\tt TMin}$)$, and $\chi(${\tt TMax}$)$; which can be found by calling {\tt axionMass.getTMin()}, {\tt axionMass.getTMax()}, {\tt axionMass.getChiMin()}, and {\tt axionMass.getChiMax()}, respectively.
-Notice that if the {\tt .cpp} is not in the root directory of \mimes, we need to compile it using the flag {\tt -Ipath-to-root}, "path-to-root" the relative path to the root directory of \mimes; \eg if the {\tt .cpp} is in the {\tt MiMeS/UserSpace/Cpp/Axion} directory, this flag should be {\tt -I../../../}.
-
-Since \mimes is designed in order to be as general as possible, one has to provide the axion mass. Therefore, before we define solve the EOM, we need to define the axion mass, as an instance of the {\tt mimes::AxionMass<LD>} class. For this example -- and the examples included in the current version -- we are going to use the mass of the QCD axion given in ref.~\cite{Borsanyi:2016ksw}. The {\tt mimes::AxionMass<LD>} class offers two ways of doing this. The first one interpolates $\chi(T)$ (defined as in \eqs{eq:axion_mass_def}) using a data file that the user provides, and the second one uses a user-defined function.
-
-For this example, we consider the QCD axion, with $\chi$ computed in ref.~\cite{Borsanyi:2016ksw}.
+\paragraph{Axion mass function}
+In some cases, the dependence of the axion mass on the temperature is known. If this is the case, the user can define the axion mass via
+%
+\begin{cpp}
+ mimes::AxionMass<LD> axionMass(ma2);
+\end{cpp}
+%
+with {\tt ma2} the axion mass squared as a callable object with signature {\tt LD ma2(LD T, LD fa)}.
-Then, we can assign an instance of {\tt mimes::Axion<LD,Solver,Method>} in a variable using
+\paragraph{The EOM solver}
+Once the axion mass is defined, we can declare a variable that will be used to solve the EOM, as
%
\begin{cpp}
mimes::Axion<LD,Solver,Method> ax(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max,
convergence_lim, inputFile, &axionMass, initial_step_size,minimum_step_size,
maximum_step_size, absolute_tolerance, relative_tolerance, beta, fac_max, fac_min,
maximum_No_steps);
\end{cpp}
%
-Here, {\tt LD} should be the desirable numeric type to be used; it is recommended to use {\tt long double}, but other choices are also available as we discuss on. Moreover {\tt Solver} and {\tt Method} depend on the type of Runge-Kutta the user chooses (see section~\ref{sec:options} asa well as the Appendices~\refs{app:RK,app:usr_input} for more details). and
+Here, {\tt LD} should be the desirable numeric type to be used; it is recommended to use {\tt long double}, but other choices are also available as we discuss later. Moreover {\tt Solver} and {\tt Method} depend on the type of Runge-Kutta the user chooses (see section~\ref{sec:options} as well as the Appendices~\refs{app:RK,app:usr_input} for more details).
- The various parameters are as follows:
+The various parameters are as follows:
%
\begin{enumerate}
\item {\tt theta\_i}: initial angle.
\item {\tt fa}: the PQ scale.
\item {\tt umax }: if $u>${\tt 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.
\item {\tt TSTOP}: if the temperature drops below this, integration stops. In most cases this should be around
$10^{-4}~\GeV$, in order to be sure that any entropy injection has stopped before integration stops (since BBN bounds~\cite{Kolb:206230,Peebles:1993} should not be violated).
\item {\tt ratio\_ini}: integration starts when $3H/\maT \approx${\tt ratio\_ini} (the exact point depends on the file ``{\tt inputFile}", which we will see later).
\item {\tt N\_convergence\_max} and {\tt convergence\_lim}: integration stops when the relative difference
between two consecutive peaks is less than {\tt convergence\_lim} for {\tt N\_convergence\_max}
consecutive peaks.
\item {\tt inputFile}: relative (or absolute) path to a file that describes the cosmology. the columns should be: $u$ $T ~[\GeV]$ $\log H$, sorted so that $u$ increases.~\footnote{One can run \run{bash MiMeS/src/FormatFile.sh inputFile} in order to sort it and remove any unwanted duplicates. See Appendix~\ref{app:util} for details of {\tt MiMeS/src/FormatFile.sh}.}
%
It is important to remember that \mimes assumes that the entropy injection has stopped before the lowest temperature of given in {\tt 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 {\tt ratio\_ini}, and {\tt TSTOP}).
\item {\tt axionMass}: an instance of the {\tt mimes::AxionMass<LD>} class, passed by pointer.
\item {\tt initial\_stepsize} (optional): initial step the solver takes.
\item {\tt maximum\_stepsize} (optional): This limits the step-size to an upper limit.
\item {\tt minimum\_stepsize} (optional): This limits the step-size to a lower limit.
\item {\tt absolute\_tolerance} (optional): absolute tolerance of the RK solver
\item {\tt 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 (\eg if {\tt f\_a} is extremely high,
the oscillations happen violently, and the system destabilizes). Whatever the case, if the
tolerances are below $10^{-8}$, long doubles have be used. \mimes by default uses {long double} variables,
in order to change it see the options available in section~\ref{sec:input}.
\item {\tt beta} (optional): controls how agreesive the adaptation is. Generally, it should be around but less than 1.
\item {\tt fac\_max}, {\tt fac\_min} (optional): the stepsize does not increase more than fac\_max, and less than fac\_min.
This ensures a better stability. Ideally, {\tt fac\_max}$=\infty$ and {\tt fac\_min}$=0$, but in reality one must
tweak them in order to avoid instabilities.
\item {\tt maximum\_No\_steps} (optional): maximum steps the solver can take Quits if this number is reached even if integration
is not finished.
\end{enumerate}
%
in order to understand the usage of the optional parameters, some basic techniques of Runge-Kutta methods are discussed in Appendix~\ref{app:RK}.
The EOM~(\ref{eq:eom_u}), then can be solved using
%
\begin{cpp}
ax.solveAxion();
\end{cpp}
%
-This, if this runs successfully (in simple cases, such as a radiation dominated Universe, this should take a fraction of a second), then $\Tosc$, $\thetaosc$, and $\Omega h^2$ are accessed via {\tt ax.T\_osc}, {\tt ax.theta\_osc}, and {\tt ax.relic}. The entire evolution (the points the integrator took) of the axion angle is stored in {\tt ax.points}, which is a two-dimensional {\tt std::vector}, with the columns corresponding to $a$, $T~[\GeV]$,
-$\theta$, $d\theta/du$, $\rho_a$. Moreover, the peaks of he oscillation are stored in another two-dimensional {\tt std::vector}, with the columns corresponding to $a$, $T~[\GeV]$, $\thetamax$, $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 {\tt mimes::Axion<LD,Solver,Method>::peaks} do not exist in {\tt mimes::Axion<LD,Solver,Method>::points}. Finally, the local errors (see Appendix~
-\ref{app:RK}) of the integration points of $\theta$ and $\zeta$ are stored in {\tt mimes::Axion<LD,Solver,Method>::dtheta} and {\tt mimes::Axion<LD,Solver,Method>::dzeta}.
+Once the EOM is solved, we can access $\Tosc$, $\thetaosc$, and $\Omega h^2$ via {\tt ax.T\_osc}, {\tt ax.theta\_osc}, and {\tt ax.relic}. The entire evolution (the points the integrator took) of the axion angle is stored in {\tt ax.points}, which is a two-dimensional {\tt std::vector}, with the columns corresponding to $a$, $T~[\GeV]$,
+$\theta$, $d\theta/du$, $\rho_a$. Moreover, the peaks of he oscillation are stored in another two-dimensional {\tt std::vector}, with the columns corresponding to $a$, $T~[\GeV]$, $\thetamax$, $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 {\tt ax.peaks} do not exist in {\tt ax.points}. Finally, the local errors (see Appendix~\ref{app:RK}) of the integration points of $\theta$ and $\zeta$ are stored in {\tt ax.dtheta} and {\tt ax.dzeta}.
-The final member function is {\tt mimes::Axion<LD,Solver,Method>::setTheta\_i}, which allows the user to set a different $\thetai$ without generating another instance.~\footnote{Since the interpolations of the data of {\tt inputFile} are made inside the constructor of the {\tt mimes::Axion<LD,Solver,Method>} class, {\tt 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
+
+\paragraph{Changing axion mass definition}
+%
+The axion mass definiton can change at any time by changing the data file (or {\tt ma2\_MIN} and {\tt ma2\_MAX}) or using using
+%
+\begin{cpp}
+ ax.set_ma2(ma2);
+\end{cpp}
+%
+with {\tt ma2} a callable object with signature {\tt LD ma2(LD T, LD fa)}.
+
+However, since integration starts at a temperature determined by {\tt 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 {\tt inputFile}. Thus, if we change the definition of the mass, we need to call
+%
+\begin{cpp}
+ ax.restart();
+\end{cpp}
+%
+This function remakes the interpolations, clears all vector, and sets all variable to $0$.
+
+
+\paragraph{Changing initial condition}
+%
+The final member function is {\tt mimes::Axion::setTheta\_i}, which allows the user to set a different $\thetai$ without generating another instance.~\footnote{Since the interpolations of the data of {\tt inputFile} are made inside the constructor of the {\tt mimes::Axion<LD,Solver,Method>} class, {\tt 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
%
\begin{cpp}
ax.setTeta_i(new_theta_ini);
\end{cpp}
%
where {\tt new\_theta\_ini} is the new value of $\thetai$. Running this function resets all variables to $0$ (except {\tt T\_osc} and {\tt a\_osc}, since they should not change), and clears all {\tt std::vector} variables, which allows the user to simply run \run{ax.solveAxion();} as if {\tt ax} was a freshly defined instance.
+
+
\subsubsection{Using \mimes in \PY}\label{sec:begin_py}
-The modules for the \PY interface are located in {\tt MiMeS/src/interfacePy}. Since, the most useful module is the one that actually solves the EOM~\ref{eq:eom_u}, we will explain here it can be used, and describe the others in Appendix~\ref{app:modules}.
+%
+The modules for the \PY interface are located in {\tt MiMeS/src/interfacePy}. Although the usage of the {\tt AxionMass} and {\tt Axion} classes is similar to the \CPP case, it is worth showing explicitly how the \PY interface works.
-The {\tt Axion} class in the module {\tt interfacePy.Axion} can be loaded in a \PY script as
+
+
+
+The two classes the modules {\tt interfacePy.AxionMass} and {\tt interfacePy.Axion}, and can be loaded in a \PY script as
%
\begin{py}
from sys import path as sysPath
from os import path as osPath
sysPath.append(osPath.join(osPath.dirname(__file__), 'path_to_src'))
+ from interfacePy.AxionMass import AxionMass
from interfacePy.Axion import Axion
\end{py}
%
It is important that \mintinline{python}{'path_to_src'} provides the relative path to the {\tt MiMeS/src} directory. For example, if the script is located in {\tt MiMeS/UserSpace/Python}, \mintinline{python}{'path_to_src'} should be \mintinline{python}{'../../src'}.
-Once the module is imported, we can define an {\tt Axion} instance as follows
+\paragraph{Axion mass definition via a data file}
+As before, we first need to define the axion mass. Similarly to the previous case, in order to define the axion mass via a file
+In many cases, axion mass cannot be written in a closed form. In these cases, the user may provide a data file. Then, the axion mass can be defined via this file as
+%
+\begin{py}
+ AxionMass axionMass(chi_PATH, minT, maxT)
+\end{py}
+%
+Here, the constructor requires the same parameters as in \CPP. Moreover, the axion mass beyond the interpolation limits can be changed via
+%
+\begin{py}
+ axionMass.set_ma2_MAX(ma2_MAX)
+ axionMass.set_ma2_MIN(ma2_MIN)
+\end{py}
+%
+Although the naming is the same as in the \CPP case, there on important difference. Namely, {\tt ma2\_MAX} and {\tt ma2\_MIN} have to be functions (that take $T$ and $\fa$ as arguments and return $\maT^2$), and cannot be any other callable object. The reason is that \mimes uses the {\tt ctypes} module, which only works with objects compatible with {\tt C}.
+%
+Moreover, the values of {\tt TMin}, {\tt TMax}, $\chi(${\tt TMin}$)$, and $\chi(${\tt TMax}$)$ can be obtained by {\tt axionMass.getTMin()}, {\tt axionMass.getTMax()}, {\tt axionMass.getChiMin()}, and {\tt axionMass.getChiMax()}, respectively.
+
+\paragraph{Axion mass definition via a function}
+%
+Again this can be done as
+%
+\begin{py}
+ AxionMass axionMass(ma2)
+\end{py}
+%
+with {\tt ma2} the axion mass squared, which should be a function (not any callable object) of $T$ and $\fa$.
+
+
+\paragraph{The EOM solver}
+%
+We can define an {\tt Axion} instance as follows
%
\begin{py}
- ax=Axion(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max, convergence_lim, inputFile,
- initial_step_size,minimum_step_size, maximum_step_size, absolute_tolerance,
- relative_tolerance, beta, fac_max, fac_min, maximum_No_steps)
+ ax=Axion(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max, convergence_lim,
+ inputFile, axionMass, initial_step_size, minimum_step_size, maximum_step_size,
+ absolute_tolerance, relative_tolerance, beta, fac_max, fac_min, maximum_No_steps)
\end{py}
%
-Here, {\tt Axion} is the constructor of the {\tt Axion} class, and {\tt ax} the variable -- which can have any name, and the input parameters are exactly the same as in the \CPP case (the usage of the class can be found by running {\tt ?Axion} after loading the module).
+Here the input parameters are the same as in the \CPP case (the usage of the class can be found by running {\tt ?Axion} after loading the module). The only slight difference is that the {\tt axionMass} is not passed as a pointer; which is done internally using {\tt ctypes}.
+
Using the defined variable ({\tt ax} in this example), we can simply run
%
\begin{py}
ax.solveAxion()
\end{py}
%
in order to solve the EOM of the axion. In contrast to the \CPP implementation, this only gives us access to $\Tosc$, $\thetaosc$, and $\Omega h^2$; the corresponding variables are {\tt ax.T\_osc}, {\tt ax.theta\_osc}, and {\tt ax.relic}. In order to get the evolution of the axion field, we need to run
%
\begin{py}
ax.getPoints()
\end{py}
%
This will make {\tt numpy} arrays that contain the scale factor ({\tt ax.a}), temperature ({\tt ax.T}), $\theta$ ({\tt ax.theta}), its derivative with respect to $u$ ({\tt ax.zeta}), and the energy density of the axion ({\tt ax.rho\_axion}).
Moreover, in order to get the various quantities on the peaks of the oscillation, we can run
%
\begin{py}
ax.getPeaks()
\end{py}
%
This makes {\tt numpy} arrays that contain the scale factor ({\tt ax.a\_peak}), temperature ({\tt ax.T\_peak}), $\theta$ ({\tt ax.theta\_peak}), its derivative with respect to $u$ ({\tt ax.zeta\_peak}, which should be equal to $0$), the energy density of the axion ({\tt ax.rho\_axion\_peak}), and the values of the adiabatic invariant on the peaks ({\tt ax.adiabatic\_invariant}).
Moreover, we can run the following
%
\begin{py}
ax.getErrors()
\end{py}
%
in order to store the local errors (see Appendix~\ref{app:RK}) of $\theta$ and $\zeta$ in {\tt ax.dtheta} and {\tt ax.dzeta}, respectively.
+\paragraph{Changing axion mass definition}
+%
+We can change the axion mass by changing the file, {\tt ma2\_MIN} and {\tt ma2\_MAX}, or using
+%
+\begin{cpp}
+ ax.set_ma2(ma2)
+\end{cpp}
+%
+with {\tt ma2} a function that takes $T$ and $\fa$, and returns the $\maT^2$. As in the \CPP case, if the definition of the mass of the axion changes (including the definitions beyond the interpolation limits), we have to call
+%
+\begin{cpp}
+ ax.restart()
+\end{cpp}
+%
+in order to remake the interpolation of cosmological quantities and reset the various variables.x
+\paragraph{Changing the initial condition}
+%
The initial condition $\thetai$ can be changed using
%
\begin{py}
ax.setTeta_i(new_theta_ini)
\end{py}
%
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 {\tt T\_osc} and {\tt a\_osc}, as they should not change).
-\paragraph{Importand} Th {\tt Axion} class constructs a {\tt mimes::Axion<LD,Solver,Method>} pointer, which has to be manually deleted. Therefore, once {\tt ax} is used it must be deleted, \ie we need to run
+\paragraph{Importand note} The {\tt Axion} class constructs a pointer to an instance of the underlying {\tt mimes::Axion} class, which has to be manually deleted. Therefore, once {\tt ax} is used it must be deleted, \ie we need to run
%
\begin{py}
del ax
\end{py}
%
We should note that this must run even if we assign another instance to the same variable {\tt ax}, otherwise we risk running out of memory.
+
\section{Assumptions and user input}\label{sec:assumptions}
\setcounter{equation}{0}
%
\mimes only makes a few, fairly general, assumptions. First of all, it is assumed that the axion energy density is always subdominant compared to radiation or any other component of the Universe, and that decays and annihilations of particles have a negligible effect on the axion energy density. Moreover, the initial condition is always assumed to be $\theta_{t=\ti} = \thetai$ and $\dot \theta|_{t=\ti}=0$.
%
Furthermore, it is also assumed that $3H/\maT$ increases monotonically at high temperatures.
%
Also, it is assumed that the entropy is resumes its conserved status at a temperature higher than the minimum temperature of {\tt inputFile} (which is required by the constructor of the {\tt mimes::Axion<LD,Solver,Method>} class).
%
Finally, \mimes does not try to predict anything regarding the cosmology. Therefore, the temperatures in {\tt inputFile} {\em must} cover the entire region of integration; \ie the maximum temperature has to be larger than the one required to reach {\tt ratio\_ini}, while the minimum one should be lower than {\tt TSTOP}.
\subsection{Options at Compile-time}\label{sec:options}
%
The user has a number of options regarding different aspects for the code. The various choices for the shared libraries used by the \PY interface are given in {\tt \mimes/Definitions.mk}. while the corresponding options for the \CPP examples are in the {\tt Definitions.mk} files in the directories inside {\tt \mimes/UserSpace/Cpp}. The options correspond to different variables, which are
%
\begin{enumerate}
\item {\tt rootDir}: the relative path of root directory of \mimes.
%
\item {\tt LONG}: this should be either {\tt long} or omitted. If omitted, the type of the numeric values is {\tt double} (double precision). On the other hand, if {\tt LONG=long}, the type is {\tt long double}. Generally, using {\tt double} should enough. For the sake of numerical stability, then, it is advised to always use {\tt 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 {\tt absolute\_tolerance} or {\tt absolute\_tolerance} are chosen to be below $\sim 10^{-8}$, then double precision numbers may not be enough, and {\tt 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 {\tt LONG} defines a macro with the same name, which then defines the macro {\tt LD} (in all the {\tt .cpp} files) as {\tt \#define LD LONG double}.
%
\item {\tt LONGpy}: the same as {\tt LONG}, but for the \PY interface.
%
\item {\tt SOLVER}: \mimes uses the ordinary differential equation (ODE) integrators of ref.~\cite{NaBBODES}. Currently, there are two available choices; {\tt Rosenbrock} and {\tt RKF}. The former is a general embedded Rosenbrock implementation and it is used if {\tt Solver=1}, while the latter is a general explicit embedded Runge-Kutta implementation and can be chose by using {\tt SOLVER=2} (a brief description of how these algorithms are implemented can be found in Appendix~\ref{app:RK}). The default option for that \mimes uses is {\tt SOLVER=1}, because the axion EOM tends to oscillate rapidly. However, in some cases, a high order explicit method may also work. This variable, is used to define a macro with the same name and value in {\tt \mimes/src/Axion/AxionSolve.hpp}.
%
\item {\tt METHOD}: Depending on the type of solver is chosen, there are some available methods.~\footnote{It is worth mentioning that {\tt 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~\ref{app:RK}. This variable, is used to define a macro with the same name and value in {\tt \mimes/src/Axion/AxionSolve.hpp}.}
- \begin{itemize}
- \item For {\tt SOLVER=1}, the available methods are
- {\tt METHOD=RODASPR2} and {\tt METHOD=ROS34PW2}. The {\tt RODASPR2} choice is a fourth order Rosenbrock-Wanner method (more information can be found in ref.~\cite{RANG2015128}). The {ROS34PW2} choice corresponds to a third order Rosenbrock-Wanner method~\cite{RangAngermann2005}.
- %
- \item For {\tt SOLVER=2}, the only reliable method available in {\tt NaBBODES} is the Dormand-Prince~\cite{DORMAND198019} chosen if {\tt METHOD=DormandPrince}, which is an explicit Runge-Kutta method of seventh order.
- \end{itemize}
+ \begin{itemize}
+ \item For {\tt SOLVER=1}, the available methods are
+ {\tt METHOD=RODASPR2} and {\tt METHOD=ROS34PW2}. The {\tt RODASPR2} choice is a fourth order Rosenbrock-Wanner method (more information can be found in ref.~\cite{RANG2015128}). The {ROS34PW2} choice corresponds to a third order Rosenbrock-Wanner method~\cite{RangAngermann2005}.
+ %
+ \item For {\tt SOLVER=2}, the only reliable method available in {\tt NaBBODES} is the Dormand-Prince~\cite{DORMAND198019} chosen if {\tt METHOD=DormandPrince}, which is an explicit Runge-Kutta method of seventh order.
+ \end{itemize}
%
\item {\tt CC}: the \CPP compiler that one chooses to use. The default option is {\tt CC=g++}, which is the {\tt GNU} \CPP compiler, and is available for all systems. Another option is to use the {\tt clang} compiler, which is chosen by {\tt CC=clang -lstdc++}. \mimes is mostly tested using {\tt g++}, but {\tt 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 \CPP17 standard, otherwise \mimes probably will not work.
%
\item {\tt OPT}: Optimization level of the compiler. By default, this is {\tt OPT=O3}, which produces executables that are marginally faster than {\tt OPT=O1} and {\tt OPT=O2}, but significantly faster than {\tt OPT=O0}. There is another choice, {\tt OPT=Ofast}, but it can cause various numerical instabilities, and is generally considered dangerous -- although we have not observed any problems when running \mimes.
\end{enumerate}
%
These options can only be used when the compilation is done using the various {\tt makefile} files. For compilation of a user provided {\tt .cpp} file, these have to be define by the user inside this file. If one compiles the examples without {\tt makefile}, the variables {\tt LONG}, {\tt LONGpy}, {\tt SOLVER}, and {\tt METHOD} correspond to macros with the same name.~\footnote{Macros can be defined using the {\tt -D} flag
-in the compiler. For example {\tt METHOD=DormandPrince} is the same as passing {\tt -DMETHOD=DormandPrince} to the compiler.} It should be pointed-out that, in the examples, {\tt LONG} is used to define a macro {\tt LD} which becomes {\tt double} or {\tt long double} accordingly. the actual choice of numeric type is passed as a template parameter in all classes.
+ in the compiler. For example {\tt METHOD=DormandPrince} is the same as passing {\tt -DMETHOD=DormandPrince} to the compiler.} It should be pointed-out that, in the examples, {\tt LONG} is used to define a macro {\tt LD} which becomes {\tt double} or {\tt long double} accordingly. the actual choice of numeric type is passed as a template parameter in all classes.
\subsection{User input}\label{sec:input}
%
\subsubsection{Compile-time input}\label{sec:compile_time_input}
%
-\paragraph{Files} \mimes requires files that provide data for the relativistic degrees of freedom (RDOF) of the plasma, as well as its mass as a function of the temperature. Although \mimes is shipped with the standard model RDOF found in~\cite{Saikawa:2020swg}, and the axion mass from~\cite{Borsanyi:2016ksw}, the user can change the corresponding variables in {\tt Paths.mk} (found in the root directory of \mimes) to another file path. The variables pointing to these data files are {\tt cosmoDat} and {\tt axMDat}, for the RDOF and axion mass, respectively.~\footnote{One can change the data file for anharmonic factor introduced in \eqs{eq:anharmonic_f} using the variable {\tt anFDat}.}
+\paragraph{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~\cite{Saikawa:2020swg}, and a few points for $f(\thetamax)$ introduced in \eqs{eq:anharmonic_f}, the user is free change them via the corresponding variables in {\tt Paths.mk} (in the root directory of \mimes). Moreover, there is a set of data for the QCD axion mass as calculated in ref.\cite{Borsanyi:2016ksw}.
+%
+The variables pointing to these data files are {\tt cosmoDat}, {\tt axMDat}, and {\tt anFDat}, for the RDOF, axion mass, and the anharmonic factor; respectively.
-The format of the files is the following:~\footnote{The anharmonic factor data are given as columns: $\thetamax$ $f(\thetamax)$; with increasing $\thetamax$.}
+The format of the files has to be the following:
%
\begin{itemize}
\item The RDOF data must be given in three columns; $T$ (in $\GeV$), $\heff$, and $\geff$.
%
\item The axion mass data must be given in two columns; $T$ (in $\GeV$), $\chi$ (in $\GeV^4$). Here, $\chi$ is defined as in \eqs{eq:axion_mass_def}.
The user can provide a function instead of data for the axion mass, by leaving the {\tt axMDat} variable empty.
+ %
+ \item The data for the anharmonic factor must be give in two columns $\thetamax$ $f(\thetamax)$; with increasing $\thetamax$.
\end{itemize}
%
The paths to these files should be given at compile time. That is, once {\tt Paths.mk} changes, \run{bash configure.sh} and then \run{\tt make} must be ran. The user can change the content of the data files (without changing their paths), in order to use them without compiling \mimes again. However, the user has to make sure that all the files are sorted so that the values of first column increase (with no duplicates or empty lines). In order to ensure this, it is advised to run \run{bash FormatFile.sh path-to-file} (in Appendix~\ref{app:util} there are some details on {\tt MiMeS/src/FormatFile.sh}), in order to format the file (that should exist in \run{path-to-file}) complies with the requirements of \mimes.
-\paragraph{Static variables} The header file {\tt MiMeS/src/static.hpp} contains definitions of various static variables that are used in the \PY modules and the {\tt mimes::Axion<LD,Solver,Method>} class. These variables are instances of {\tt mimes::Cosmo<LD>}, {\tt mimes::AxionMass<LD>}, and {\tt mimes::AnharmonicFactor<LD>}, and they are responsible for interpolation of the various plasma parameters (\ie $\heff(T)$, $\geff(T)$, $s(T)$, etc.), the axion mass, and the anharmonic factor, respectively.
-
-\subparagraph{{\tt cosmo<LD>}} The static variable {\tt cosmo<LD>} is an instance of {\tt mimes::Cosmo<LD>}, with {\tt LD} being a numeric type (\eg {\tt double}, {\tt long double}). It is expected to be defined in {\tt MiMeS/src/static.hpp}, as
-%
-\begin{cpp}
- template<class LD> static mimes::Cosmo<LD> cosmo(cosmo_PATH,minT,maxT);
-\end{cpp}
-%
-The arguments of the constructor are
-%
-\begin{enumerate}
- \item {\tt cosmo\_PATH}: a string of the absolute path of the relativistic degrees of freedom (RDOF). This macro is automatically
- written in {\tt MiMeS/src/misc\_dir/path.hpp} when \run{bash configure.sh} is ran. The RDOF data must be given in three columns; $T$ (in $\GeV$), $\heff$, and $\geff$.
- %
- \item {\tt minT} (optional): The minimum interpolation temperature. Bellow this limit, $\heff$ and $\geff$ are takes to be constant. Its default value is $0$.
- %
- \item {\tt maxT} (optional): The maximum interpolation temperature. Above this limit, $\heff$ and $\geff$ are takes to be constant. Its default value is $m_P$; the Planck mass $m_P = 1.22 \times 10^{19}~\GeV$.
-\end{enumerate}
-%
-It should be noted that the interpolation limits {\tt minT} and {\tt maxT} are suggestions; the actual limits will be set to their closest value that exist in the data file. Apart form the various member functions, there are several static variables; \eg the CMB temperature today ({\tt mimes::Cosmo<LD>::T0}), the Planck mass ({\tt mimes::Cosmo<LD>::mP}), etc. The values of all these parameters are taken as compile-time constants. However, the user can change their values directly from the class, which is located in {\tt \mimes/src/Cosmo/Cosmo.hpp}.
-
-
-\subparagraph{{\tt axionMass<LD>}} The static variable {\tt axionMass<LD>} is an instance of {\tt mimes::AxionMass<LD>}. It is expected to be defined in {\tt MiMeS/src/static.hpp}. The first way one can define the variable {\tt axionMass<LD>} is
-%
-\begin{cpp}
- template<class LD> static mimes::AxionMass<LD> axionMass(chi_PATH,minT,maxT);
-\end{cpp}
-%
-The arguments of the constructor are
-%
-\begin{enumerate}
- \item {\tt chi\_PATH}: a string of the absolute path of the data for $\chi$ (defined as in \eqs{eq:axion_mass_def}). This macro is automatically
- written in {\tt MiMeS/src/misc\_dir/path.hpp} when \run{bash configure.sh} is ran. The data must be given in two columns; $T$ (in $\GeV$), $\chi$ (in $\GeV^4$).
- %
- \item {\tt minT} (optional): The minimum interpolation temperature. Its default value is $0$.
- %
- \item {\tt maxT} (optional): The maximum interpolation temperature.Its default value is $m_P$.
-\end{enumerate}
-%
-As previously, the interpolation limits {\tt minT} and {\tt maxT} are suggestions, as the actual limits will be set to the closest value that exist in the data file. Beyond these limits, the axion mass by default is taken to be
-%
-\begin{equation}
- \maT = \Bigg\{
- \begin{matrix}
- \dfrac{\chi(T_{\rm min})}{\fa^2} & \text{for } T<T_{\rm min}
- \\ \\
- \dfrac{\chi(T_{\rm max})}{\fa^2} \lrb{\dfrac{T}{T_{\rm max}}}^{-8.16} & \text{for } T>T_{\rm max}
- \end{matrix} \;.
- \label{eq:axM-limits}
-\end{equation}
-%
-This is a good approximation for the QCD axion (the default data \mimes uses), but the user can change them inside the class, which is located in{\tt \mimes/src/AxionMass/AxionMass.hpp} (see section~\ref{sec:cpp_example} or Appendix~\ref{app:classes} for details).
+These paths are stored as strings in {\tt \mimes/src/misc\_dir/path.hpp} at compile-time (they are defined as {\tt constexpr}), and can be accessed once this header file is included. The corresponding variables are {\tt cosmo\_PATH}, {\tt chi\_PATH}, and {\tt anharmonic\_PATH}, for the path to data file of the plasma quantities, $\chi(T)$, and $f(\thetamax)$; respectively. Although, the axion mass data file may be omitted -- since the axion mass is defined by the user, the variable {\tt chi\_PATH} is still useful if the axion mass is defined via a data file, as it is automatically converted to an absolute path.~\footnote{Absolute 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.}
-The second way we can declare the {\tt axionMass<LD>} variable is
-%
-\begin{cpp}
- template<class LD> static mimes::AxionMass<LD> axionMass(ma2);
-\end{cpp}
-%
-where {\tt ma2} -- the axion mass squared in $\GeV^2$ -- is a function object with signature {\tt LD ma2 (LD T, LD fa)}. In order to do this,
-the user may omit a data file path for the {\tt axMDat} variable in {\tt \mimes/src/Definition.mk}.
-\subparagraph{{\tt anharmonicFactor<LD>}} The static variable {\tt anharmonicFactor<LD>} is an instance of the class {\tt mimes::AnharmonicFactor<LD>}. It is expected to be defined in {\tt MiMeS/src/static.hpp}, as
-%
-\begin{cpp}
- template<class LD> static mimes::AnharmonicFactor<LD> anharmonicFactor(anF_PATH);
-\end{cpp}
-%
-Here,~{\tt anF\_PATH} is a macro of a string of the absolute path of the data for the anharmonic factor. This macro is automatically
-written in {\tt MiMeS/src/misc\_dir/path.hpp} when \run{bash configure.sh} is ran. The data should be given in two columns; $\thetamax$, $f(\thetamax)$.
\subsubsection{Run-time input}\label{sec:run_time_input}
%
The run-time user input is described in sec.~\ref{sec:first_steps}. The user has to provide at least the parameters that describe the axion evolution, $\thetai$ and $\fa$.
Moreover, the maximum allowed value of $u$ and the minimum value of $T$, allow the user to decide when the integration has to stop even if the axion has not reached its adiabatic evolution. Ideally, {\tt umax}$=\infty$ and {\tt 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.~\footnote{These two variables are not optional, because the user must be aware of them, in order to choose them according to their needs.}
Furthermore, {\tt ratio\_ini} allows the user to choose a desired point at which the interpolation of the data in {\tt inputFile} begins. This can save valuable computation time as well as memory, as only the necessary data are stored and searched. Generally, for {\tt ratio\_ini}$>10^{3}$, the relic abundance becomes
independent from {\tt ratio\_ini}, but one has to choose it carefully, in order to find a balance between accuracy and computation time.
Finally, the convergence conditions -- \ie~{\tt N\_convergence\_max} and {\tt 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 log as they are {\tt N\_convergence\_max}$>3$ and {\tt convergence\_lim}$<10^{-1}$ (\ie the adiabatic invariant does not vary more that $10\%$ for three consecutive peaks of the oscillation).
%
However, we should note that greedy choices (\eg~{\tt N\_convergence\_max}$=100$ and {\tt convergence\_lim}$=10^{-5}$) are dangerous, as $\theta$ tends to oscillate rapidly and destabilize the differential equation solver. Therefore, these parameters should be chosen carefully, in order to ensure that integration stops when axion has reached its adiabatic evolution, without destabilising the EOM.
-
-
\subsection{Complete Examples}\label{sec:complete_examples}
%
Although one can modify the examples provided in {\tt \mimes/UserSpace}, in this section we show a complete example in both \CPP and \PY.
%
The underlying cosmology is assumed to be an EMD scenario. In terms of the parametrisation of ref.~\cite{Arias:2020qty}, it is defined through the same parameters as in the example shown in section~\ref{sec:comparison_between_WKB_approaches}.
%
and the evolution of the energy densities of the plasma and a fluid $\Phi$ is shown in \Figs{fig:energy_densities}.
%
\begin{figure}[t]
\includegraphics[width=1\textwidth]{figs/energy_densities.pdf}
\caption{The energy densities of the plasma and a hypothetical decaying fluid. The parameters of the cosmological scenario $T_{\rm END}=10^{-2}~\GeV$, $c=3$, $T_{\rm ini}=10^{12}~\GeV$, and $r=0.1$; the parameters follow the definitions in ref.~\cite{Arias:2020qty}. }
\label{fig:energy_densities}
\end{figure}
%
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~\cite{Arias:2020qty}). Regarding the axion, this cosmological scenario alters it evolution, since both the Hubble parameter and the entropy of the plasma change significantly. This results in a shift in the oscillation temperature, and the dilution of the axion energy density.
+Moreover, we will use the QCD axion mass of ref.~\cite{Borsanyi:2016ksw}. For the axion mass beyond the minimum and maximum temperatures, $T_{\rm min, \, max}$, in the corresponding data file, we will take
+%
+\begin{equation}
+ \maT = \Bigg\{
+ \begin{matrix}
+ \dfrac{\chi(T_{\rm min})}{\fa^2} & \text{for } T<T_{\rm min}
+ \\ \\
+ \dfrac{\chi(T_{\rm max})}{\fa^2} \lrb{\dfrac{T}{T_{\rm max}}}^{-8.16} & \text{for } T>T_{\rm max}
+ \end{matrix} \;.
+ \label{eq:axM-limits}
+\end{equation}
+
+
\subsubsection{complete example in \CPP}\label{sec:cpp_example}
%
-In order to write a \CPP program that uses \mimes in order to solve the EOM~\ref{eq:eom_sys}, we must include the header file {\tt src/Axion/AxionSolve.hpp}. Once this is done, we can declare various instances, and solve the axion EOM. For example in order to to do so in the cosmological scenario at hand, the {\tt main} function should contain
+In order to write a \CPP program that uses \mimes in order to solve the EOM~\ref{eq:eom_sys}, we must include the header files {\tt src/AxionMass/AxionMss.hpp} and {\tt src/Axion/AxionSolve.hpp}. In the example at hand, the {\tt main} function should contain the definition of the axion mass
+%
+\begin{cpp}
+ // use chi_PATH to interpolate the axion mass.
+ mimes::AxionMass<LD> axionMass(chi_PATH,0,mimes::Cosmo<LD>::mP);
+
+ // This is the axion mass squared beyond the interpolation limits for the current data.
+ // If you don't specify them, the axion mass is taken to be constant beyond these limits
+
+ /*set ma2 for T>TMax*/
+ long double TMax=axionMass.getTMax();
+ long double chiMax=axionMass.getChiMax();
+ axionMass.set_ma2_MAX(
+ [&chiMax,&TMax](long double T, long double fa){
+ return chiMax/fa/fa*std::pow(T/TMax,-8.16);
+ }
+ );
+
+ /*set ma2 for T<TMin*/
+ long double TMin=axionMass.getTMin();
+ long double chiMin=axionMass.getChiMin();
+ axionMass.set_ma2_MIN(
+ [&chiMin,&TMin](long double T, long double fa){
+ return chiMin/fa/fa;
+ }
+ );
+\end{cpp}
+%
+It should be noted that, since we have used {\tt chi\_PATH}, we need to include the {\tt path.hpp} header file from {\tt \mimes/src/misc\_dir}. Moreover, we note that we have used as parameter {\tt mimes::Cosmo<LD>::mP}, which is the Planck mass, as a static member of the {\tt Cosmo} class (described in Appendix~\ref{app:classes}), in order to interpolate the mass up to $T=m_P$, \ie the entire possible mass range.
+
+Following, the axion EOM can be solve as
%
\begin{cpp}
std::string inputFile = std::string(rootDir)+
std::string("/UserSpace/InputExamples/MatterInput.dat");
mimes::Axion<long double, 1, RODASPR2<long double> > ax(0.1, 1e16, 500, 1e-4, 1e3,
- 15, 1e-2, inputFile, 1e-1, 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7));
+ 15, 1e-2, inputFile, &axionMAss, 1e-1, 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7));
ax.solveAxion();
\end{cpp}
%
-This program should solve the axion EOM. However, one needs to print anything that they need. For example, the relic abundance is printed by adding {\tt std::cout<<ax.relic<<"\\n";} after {\tt 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 {\tt \#include<iostream>}, in order to be able to use {\tt std::cout}.
+This program should solve the axion EOM. However, one needs to print anything that they need. For example, the relic abundance is printed by adding {\tt std::cout<<ax.relic<<"\\n";} after {\tt 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 {\tt \#include<iostream>}, in order to be able to use {\tt std::cout}.
\paragraph{Parameter choice}
%
It should be noted that the string that is assigned to the variable {\tt inputFile}, is the path of file that exists in {\tt \mimes/UserSpace/InputExamples/MatterInput.dat}. Also, {\tt rootDir} is a constant static ({\tt char[]}) variable that is defined in {\tt \mimes/src/path.hpp} that , and it is automatically generated when {\tt bash configure.sh} is executed. Moreover, the other parameters used in the constructor are:
%
\begin{itemize}
\begin{minipage}{0.3\linewidth}
\item {\tt theta\_i}$=0.1$
\item {\tt fa}$=10^{16}~\GeV$
\item {\tt umax }$=500$
\item {\tt TSTOP}$=10^{-4}~\GeV$
\item {\tt ratio\_ini}$=10^3$
\item {\tt N\_convergence\_max}$=15$
\end{minipage}
\begin{minipage}{0.35\linewidth}
\item {\tt convergence\_lim}$=10^{-2}$
\item {\tt initial\_stepsize}$=10^{-1}$
\item {\tt maximum\_stepsize}$=10^{-1}$
\item {\tt minimum\_stepsize}$=10^{-8}$
\item {\tt absolute\_tolerance}$=10^{-11}$
\item {\tt relative\_tolerance}$=10^{-11}$
\end{minipage}
\begin{minipage}{0.3\linewidth}
\item {\tt beta}$=0.9$
\item {\tt fac\_max}$=1.2$
\item {\tt fac\_min}$=0.8$
\item {\tt maximum\_No\_steps}$=10^7$
\item[]{\vfill} \item[]{\vfill} \item[]{\vfill}
\end{minipage}
\end{itemize}
\paragraph{Template arguments}
%
Furthermore, The first template parameter is {\tt long double}, which means that all the numeric types (apart from integers like {\tt 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, n this case, the second template parameter -- $1$ -- means that we use a Rosenbrock method, with the method being {\tt RODASPR2}. Notice that the method also needs a template parameter that declares all member variables of the {\tt RODASPR2} class.
\paragraph{Compilation}
%
Assuming that we have already executed {\tt bash configure.sh}, compilation of the above example is trivial. Assuming that we name the file that contains the code {\tt axionExample.cpp}, and it is located in {\tt \mimes/UserSpace}, the only thing we need to do is run
%
\begin{pseudo}
g++ -O3 -std=c++17 -lm -I../ -o axion axionExample.cpp
\end{pseudo}
%
or
%
\begin{pseudo}
clang -lstdc++ -O3 -std=c++17 -lm -I../ -o axion axionExample.cpp
\end{pseudo}
Both of these commands should create an executable that solves the axion EOM. That is, assuming we have a terminal open in {\tt \mimes/UserSpace}, we can run {\tt ./axion}, and get the results we chose. It should be noted that all paths that \mimes uses by default are written as absolute paths when we run {bash ./configure.sh}. Therefore, if the {\tt 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.
-\paragraph{Axion mass limits}
-%
-In this example, we have used the default data files static variables for the plasma quantities and the axion mass. As already mentioned, one can change the variables {\tt cosmoDat} or {\tt axMDat} in {\tt \mimes/Paths.mk}, in order to change the path of the corresponding data file.
-
-However, using a data file for the axion mass, requires some functions that will return the axion mass for temperatures beyond the interpolation limits. \mimes, by default, uses \eqs{eq:axM-limits} that correspond to the QCD axion, and are valid approximate limits for the default data file. In the constructor, these two functions are defined as
-%
-\begin{cpp}
- /*$\maT^2$ for $T>T_{\rm max}$*/
- ma2_MAX=[this](LD T, LD fa){return this->chiMax/fa/fa*std::pow(T/this->TMax,-8.16);};
- /*$\maT^2$ for $T<T_{\rm max}$*/
- ma2_MIN=[this](LD T, LD fa){return this->chiMin/fa/fa;};
-\end{cpp}
-%
-where {\tt ma2\_MAX} and {\tt ma2\_MIN} are the axion mass squared for temperature above the maximum and below the minimum interpolation limit; $T_{\rm max}$ and $T_{\rm min}$, respectively. Frthermore, {\tt chiMax} and {\tt chiMin} correspond to $\chi(T_{\rm max})$ and $\chi(T_{\rm min})$, respectively.
-
-Therefore, if we use a different set of data for $\chi$ it is possible that these definitions would also need to be altered. This can be done directly in the constructor, by changing this aforementioned lines. Alternatively, since {\tt ma2\_MAX} and {\tt ma2\_MIN} are public variables, they can be changed inside the {\tt main} function -- before the declaration of the {\tt mimes::Axion<LD, Solver, Method>} instance we intend to use -- as any lambda function is \CPP. For example, in order to implement $\maT^2 = 0$ for $T>T_{\rm max}$, we simply need to change the definition of {\tt ma2\_MAX} in the constructor to
-%
-\begin{cpp}
- /*$\maT^2$ for $T>T_{\rm max}$*/
- ma2_MAX=[this](LD T, LD fa){return 0;};
-\end{cpp}
-%
-or add the following line before the declaration of {\tt ax}
-%
-\begin{cpp}
- /*$\maT^2$ for $T>T_{\rm max}$*/
- axionMass<long double>.ma2_MAX=[](LD T, LD fa){return 0;};
-\end{cpp}
-
-Then everything else -- \eg compilation, calling the member functions of the {\tt mimes::Axion<LD, Solver, Method>} class -- is exactly the same as before.
\paragraph{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
%
\begin{cpp}
#include<iostream>
+ //include everything you need from \mimes
#include"src/Axion/AxionSolve.hpp"
+ #include"src/AxionMass/AxionMass.hpp"
+ #include"src/Cosmo/Cosmo.hpp"
+ #include"src/misc_dir/path.hpp"
int main(){
- /*$\maT^2$ for $T>T_{\rm max}$*/
- axionMass<long double>.ma2_MAX=[](LD T, LD fa){return 0;};
-
+
+ // use chi_PATH to interpolate the axion mass.
+ mimes::AxionMass<LD> axionMass(chi_PATH,0,mimes::Cosmo<LD>::mP);
+
+ /*set $\maT^2$ for $T\geq T_{\rm max}$*/
+ long double TMax=axionMass.getTMax();
+ long double chiMax=axionMass.getChiMax();
+ axionMass.set_ma2_MAX(
+ [&chiMax,&TMax](long double T, long double fa){
+ return chiMax/fa/fa*std::pow(T/TMax,-8.16);
+ }
+ );
+
+ /*set $\maT^2$ for $T\leq T_{\rm min}$*/
+ long double TMin=axionMass.getTMin();
+ long double chiMin=axionMass.getChiMin();
+ axionMass.set_ma2_MIN(
+ [&chiMin,&TMin](long double T, long double fa){
+ return chiMin/fa/fa;
+ }
+ );
+
std::string inputFile = std::string(rootDir)+
std::string("/UserSpace/InputExamples/MatterInput.dat");
mimes::Axion<long double, 1, RODASPR2<long double> > ax(0.1, 1e16, 500, 1e-4, 1e3,
15, 1e-2, inputFile, 1e-1, 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7));
ax.solveAxion();
std::cout<<ax.relic<<"\n";
return 0;
}
\end{cpp}
%
-\paragraph{Final note}
%
-We should point out that another example is given in {\tt \mimes/UserSpace/Cpp/Axion}, where all parameters are taken as command-line inputs, the various compilation-time options can be given {\tt \mimes/UserSpace/Cpp/Axion/Definitions.mk}, and then compiled using {\tt 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 {\tt ma2\_MAX} and {\tt ma2\_MIN} -- or change the {\tt axionMass<LD>} static variable (in {\tt \mimes/src/static.hpp}) using a function as mentioned in section~\ref{sec:compile_time_input} without writing the entire program themselves.
+We should point out that another example is given in {\tt \mimes/UserSpace/Cpp/Axion}, where all parameters are taken as command-line inputs, the various compilation-time options can be given {\tt \mimes/UserSpace/Cpp/Axion/Definitions.mk}, and then compiled using {\tt 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 {\tt ma2\_MAX} and {\tt ma2\_MIN} -- or change the {\tt axionMass} variable (in {\tt \mimes/src/static.hpp}) using a function as mentioned in section~\ref{sec:compile_time_input} without writing the entire program themselves.
+
+\paragraph{Alternative axion mass definition}
+%
+For this particular example, we could have used the approximate definition of the axion mass
+%
+\begin{equation}
+ \maT^2 \approx \Bigg\{
+ \begin{matrix}
+ \ma^2 \times \lrb{\dfrac{T_{\rm QCD}}{T}}^{8.16} & \quad \text{for } T\geq T_{\rm QCD}
+ \\ \\
+ \ma^2 & \quad \text{for } T\leq T_{\rm QCD}
+\end{matrix} \;,
+ \label{eq:axM_approx}
+\end{equation}
+%
+where $T_{\rm QCD} \approx 150~\MeV$ and $\ma=\dfrac{3.2 \times 10^{-5} ~\GeV^4}{\fa^2}$. This can be done by substituting the {\tt axionMass} definition (lines $10$-$29$) with
+%
+\begin{cpp}
+ auto ma2 = [](long double T,long double fa){
+ long double TQCD=150*1e-3;
+ long double ma20=3.1575e-05/fa/fa;
+ if(T<=TQCD){return ma20;}
+ else{return ma20*std::pow((TQCD/T),8.16);}
+ };
+ mimes::AxionMass<LD> axionMass(ma2);
+\end{cpp}
\subsubsection{complete example in \PY}
%
-In order to be able use the {\tt Axion} class in \PY, we need to import the corresponding module from {\tt \mimes/src/interfacePy}. That is, assuming that the script from which we intend to import {\tt \mimes/src/interfacePy/Axion/Axion.py} is in {\tt \mimes/UserSpace}, we can import the {\tt Axion} class by writing the following
+In order to be able use the {\tt AxionMass} and {\tt Axion} classes in \PY, we need to import the corresponding module from {\tt \mimes/src/interfacePy}. That is, assuming that the script from which we intend to import {\tt \mimes/src/interfacePy/Axion/Axion.py} is in {\tt \mimes/UserSpace}, we can import the {\tt Axion} class by writing the following
%
\begin{py}
#add the relative path for MiMeS/src
from sys import path as sysPath
from os import path as osPath
sysPath.append(osPath.join(osPath.dirname(__file__), '../src'))
- #import the Axion class from the Axion module
- from interfacePy.Axion import Axion
+
+ from interfacePy.AxionMass import AxionMass #import the AxionMass class
+ from interfacePy.Axion import Axion #import the Axion class
+ from interfacePy.Cosmo import mP #import the Planck mass
\end{py}
%
-Once the class is imported, we can simply follow the steps outlined in section~\ref{sec:begin_py}. For the example at hand, we can create an instance of the axion class as
+Once everything we need is imported, we can simply follow the steps outlined in section~\ref{sec:begin_py}. For the example at hand, we can create an instance of the {\tt AxionMAss} class as
+%
+\begin{py}
+ # AxionMass instance
+ axionMass = AxionMass(r'../src/data/chi.dat',0,mP)
+
+ # This is the axion mass squared beyond the interpolation limits for the current data.
+ # If yoy don't specify them, the axion mass is taken to be constant beyond these limits
+ TMax=axionMass.getTMax()
+ chiMax=axionMass.getChiMax()
+ TMin=axionMass.getTMin()
+ chiMin=axionMass.getChiMin()
+
+ axionMass.set_ma2_MAX( lambda T,fa: chiMax/fa/fa*pow(T/TMax,-8.16) )
+ axionMass.set_ma2_MIN( lambda T,fa: chiMin/fa/fa )
+\end{py}
+%
+
+Then we can simply create an instance of the {\tt Axion} class as
+%
\begin{py}
#in python it is more convenient to use relative paths
inputFile = inputFile="./InputExamples/MatterInput.dat"
- ax = Axion(0.1, 1e16, 500, 1e-4, 1e3, 15, 1e-2, inputFile, 1e-1,
- 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7))
+ ax = Axion(0.1, 1e16, 500, 1e-4, 1e3, 15, 1e-2, inputFile, axionMass,
+ 1e-1, 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7))
\end{py}
%
Again, the parameters for the constructor are the same as in the \CPP example. The axion EOM, then, is solved using
%
\begin{py}
ax.solveAxion()
\end{py}
%
In contrast to the \CPP usage of this function, this only stores the $\thetai$, $\fa$, $\thetaosc$, $\Tosc$, and $\Omega h^2$ in the variables {\tt ax.theta\_i}, {\tt ax.fa}, {\tt ax.theta\_osc}, {\tt ax.T\_osc}, and {\tt ax.relic}; respectively. Therefore, we can print the relic abundance by calling \run{print(ax.relic)}. In order to get the integration points (\ie the evolution of the angle and the other quantities), the quantities at the peaks, and the local integration errors, we need to call
%
\begin{py}
#this gives you all the points of integration
ax.getPoints()
#this gives you the peaks of the oscillation
ax.getPeaks()
#this gives you local integration errors
ax.getErrors()
\end{py}
%
The documentation of any \PY function can be read directly inside the script using {\tt ?} as a prefix. For example, in order to see what the functionality and usage of the {\tt getPeaks} function, we can call \run{?ax.getPeaks}, and its documentation will be printed.
As already mentioned, it is important to always delete any instance of the {\tt Axion} class that is not going to be used. In this case this is done by calling
%
\begin{py}
del ax
\end{py}
-\paragraph{Axion mass limits}
-%
-Before we use the \PY interface of \mimes, we may need to change the default definitions of the axion mass beyond the interpolation points which can only be done inside the definition of the {\tt mimes::AxionMass} constructor (see the \CPP example in section~\ref{sec:cpp_example}), since we cannot reassign class variables inside header files.
-
-
\paragraph{Compilation of the shared library}
%
As described in section~\ref{sec:compile_time_input}, we may need to change the default data file paths, or the various compilation options. This is done through the variables in {\tt \mimes/Definitions.mk} and {\tt \mimes/Paths.mk} described in section~\ref{sec:options}.
-Once have chosen everything according to our needs, the library can be created by opening a terminal inside the root directory of \mimes and running
+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
%
\begin{py}
bash configure.sh
make lib/Axion_py.so
\end{py}
%
\paragraph{The entire code}
%
As in the \CPP example, this example consists of only a few lines of code. The script we described here is
%
\begin{py}
#add the relative path for MiMeS/src
from sys import path as sysPath
from os import path as osPath
sysPath.append(osPath.join(osPath.dirname(__file__), '../src'))
- #import the Axion class from the Axion module
- from interfacePy.Axion import Axion
+ from interfacePy.AxionMass import AxionMass #import the AxionMass class
+ from interfacePy.Axion import Axion #import the Axion class
+ from interfacePy.Cosmo import mP #import the Planck mass
+
+ # AxionMass instance
+ axionMass = AxionMass(r'../src/data/chi.dat',0,mP)
+
+ # define $\maT^2$ for $T\leq T_{\rm min}$
+ TMin=axionMass.getTMin()
+ chiMin=axionMass.getChiMin()
+ axionMass.set_ma2_MIN( lambda T,fa: chiMin/fa/fa )
+
+ # define $\maT^2$ for $T\geq T_{\rm max}$
+ TMax=axionMass.getTMax()
+ chiMax=axionMass.getChiMax()
+ axionMass.set_ma2_MAX( lambda T,fa: chiMax/fa/fa*pow(T/TMax,-8.16) )
#in python it is more convenient to use relative paths
inputFile = inputFile="./InputExamples/MatterInput.dat"
- ax = Axion(0.1, 1e16, 500, 1e-4, 1e3, 15, 1e-2, inputFile, 1e-1,
- 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7))
+ ax = Axion(0.1, 1e16, 500, 1e-4, 1e3, 15, 1e-2, inputFile, axionMass,
+ 1e-1, 1e-8, 1e-1, 1e-11, 1e-11, 0.9, 1.2, 0.8, int(1e7))
ax.solveAxion()
print(ax.relic)
#once we are done we should run the destructor
del ax
\end{py}
-\paragraph{Final note}
-%
One can find a complete example, including the option to create several plots, in {\tt \mimes/UserSpace/Python/Axion.py}. Also, the same example can be used interactively, in the jupyter notebook, that can be found in {\tt \mimes/UserSpace/JupyterNotebooks/Axion.ipynb}. One can read the comments, and change all different parameters, in order to examine how the results are affected.
-
+\paragraph{Alternative axion Mass definition}
+%
+In order to use the approximate axion mass as defined in \eqs{eq:axM_approx}, we could replace the definition of the {\tt axionMass} variable (lines $10$-$21$) with
+%
+\begin{py}
+ def ma2(T,fa):
+ TQCD=150*1e-3
+ ma20=3.1575e-05/fa/fa
+ if T<=TQCD:
+ return ma20;
+ return ma20*pow((TQCD/T),8.16)
+ axionMass = AxionMass(ma2)
+\end{py}
\subsubsection{Results}
%
In both \CPP and \PY one should get the same results (if the same numeric type, and RK parameters and method are used), which should be similar to \Figs{fig:theta_evolution-EMD}.
%
It should be noted that this particular example is chosen such that the axion relic abundance is close to its observed value~\cite{Planck:2018vyg} ($\Omega h^2 \approx 0.12$), and is indicative of the features one can expect if entropy is injected to the plasma for $T \ll \Tosc$.
%
\begin{figure}[h]
\begin{subfigure}[]{0.5\textwidth}
\includegraphics[width=1\textwidth]{figs/theta_evolution-EMD.pdf}
\caption{}
\label{fig:theta_evolution-EMD}
\end{subfigure}
\begin{subfigure}[]{0.5\textwidth}
\includegraphics[width=1\textwidth]{figs/theta_evolution-RD.pdf}
\caption{}
\label{fig:theta_evolution-RD}
\end{subfigure}
\caption{...}
\label{fig:results}
\end{figure}
%
In \Figs{fig:theta_evolution-EMD} we show the evolution of $\theta$ for temperatures $T \in [0.025,0.15]~\GeV$, where the vertical line indicates $\Tosc$, while the horizontal one the corresponding value of $\theta$, $\thetaosc$. Also, the blue curve connects the peaks of the oscillation. For comparison, in \Figs{fig:theta_evolution-RD} we also show the evolution of the angle in a radiation dominated Universe with constant entropy (\ie 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, as expected from \eqs{eq:theta_osc}, the angle at $T=\Tosc$ in the EMD scenario is much smaller that the corresponding value in the standard cosmological case. This shows again that approximating $\thetaosc \approx \thetai$ results in an inaccurate result, especially in cases where entropy is injected to the plasma (as expected by \eqs{eq:theta_osc}).
%
Moreover, in \Figs{fig:RK_response}, we show the relative local error of integration as well as a histogram of the number of integration steps for $0.025~\GeV < T < 0.15 ~\GeV$.
%
\begin{figure}[h]
\begin{subfigure}[]{0.5\textwidth}
\includegraphics[width=1\textwidth]{figs/local_errors-EMD.pdf}
\caption{}
\label{fig:local_errors-EMD}
\end{subfigure}
\begin{subfigure}[]{0.5\textwidth}
\includegraphics[width=1\textwidth]{figs/histogram-EMD.pdf}
\caption{}
\label{fig:histogram-EMD}
\end{subfigure}
\caption{...}
\label{fig:RK_response}
\end{figure}
%
The local relative errors, defined in Appendix~\ref{app:RK}, are shown in \Figs{fig:local_errors-EMD}. 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>\Tosc$, 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 relative error of $zeta$ is large because $\zeta \approx 0$ before the $\Tosc$). In order to examine how the step-size is adapted to the difficulty of the problem, we also show a histogram which show how many integration steps are taken for fixed temperature intervals.~\footnote{More precisely, the histogram is made by dividing the the temperature range $0.025~\GeV\lesssim T\lesssim 0.15~\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 -- $\maT/H$ increases rapidly for $T<\Tosc$. 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 \Tosc$ 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPENDICES
\pagebreak
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setcounter{section}{0}
\section*{Appendix}
\appendix
\renewcommand{\theequation}{\Alph{section}.\arabic{equation}}
\setcounter{equation}{0} % reset counter
%%%%%%%%%%%%%%%%%%%%%%%
\section{Basics of embedded Runge-Kutta Mehtods}\label{app:RK}
\setcounter{equation}{0}
%
Runge-Kutta (RK) methods are employed in order to solve an ordinary differential equation (ODE), or a system of ODEs of first order.~\footnote{Boundary value problems, and higher order differential equations are expressed as first order ODEs, and then solved. Similarly to \eqs{eq:eom_sys}.} Although there are some very insightful sources in the literature (\eg~\cite{Hairer,hairer2010solving,10.5555/1403886}) 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
%
\begin{equation}
\dfrac{d\vec{y}}{dt}=\vec{f}(\vec{y},t) \;,
\label{eq:ODE_definition}
\end{equation}
%
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 the system of \eqs{eq:ODE_definition}, an RK method uses an iteration of the form
%
\begin{equation}
\vec{y}_{n+1}=\vec{y}_{n}+ h\sum_{i=1}^{s} b_i \ \vec{k}_i \;,
\label{eq:RK_iter}
\end{equation}
%
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 (\ie $\vec{k} =\vec{k}\lrsb{y_{n},t_{n}}$), are called explicit, while the ones that try to also predict next step (\ie $\vec{k} =\vec{k}\lrsb{y_{n},t_{n}; y_{n+1},t_{n+1}}$) are called implicit.~\footnote{Generally, by substituting implicit methods $y_{n+1}$ as in \eqs{eq:RK_iter}, we end up with a system of equations that need to be solved in order to compute $\vec{k}$.}
\subsection{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 \eqs{eq:RK_iter}, a second estimate is given by
%
\begin{equation}
\vec{y}_{n+1}^{(*)}=\vec{y}_{n}+ h\sum_{i=1}^{s} b_i^{(*)} \vec{k}_i \;,
\label{eq:RK_Embedded_iter}
\end{equation}
%
with $b_i^{(*)}$ is an extra parameter that characterise the ``embedded" method of different order (typically, one order higher that the estimate~\ref{eq:RK_iter}). The local (for the step $t_{n+1}$) error divided by the scale of the solution, then, estimated as
%
\begin{equation}
\Delta \equiv \sqrt{\frac{1}{N}\sum_{d=1}^{N}\lrb{\frac{y_{n+1 ,d}-y^{*}_{n+1,d}}{ \Lambda_d }}^{2}} \;,
\label{eq:error_estimate}
\end{equation}
%
where $n$ the iteration number, $N$ the number of ODEs, $d$ the component of $\vec{y}$, and $\Lambda_d$ defined as
%
\begin{equation}
\Lambda_d = {\tt Atol} + \max\lrb{|y_{n+1, d}|,|y^{*}_{n+1 ,d}|} \; {\tt Rtol} \;,
\label{eq:RK_scale}
\end{equation}
%
with {\tt Atol} and {\tt Rtol} the absolute and relative tolerances that characterise the desirable accuracy we want to achieve; user defined values, typically {\tt Atol}$=${\tt Rtol}$\ll 1$. With these definitions, the desirable error is reached when $\Delta \lesssim 1$.
\paragraph{Step-control} The definition~\ref{eq:error_estimate}, 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
%
\begin{eqnarray}
h \to \beta \; h \; {\rm max} \lrsb{f_{\rm min}, {\rm min} \lrb{\Delta^{-\frac{1}{p+1}}, f_{\rm max}}} \;,
\label{eq:step-control}
\end{eqnarray}
%
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.
\paragraph{Correspondence between \mimes parameters and RK ones}
%
The various parameters that \mimes are described in section~\ref{sec:input}. The correspondence between them and the RK parameters is given in table~\ref{tab:RK_mimes_params}.
%
\begin{table}[t!]
\centering
\begin{tabular}{|c|c|}
\mimes & Runge-Kutta \\
\hline
{\tt absolute\_tolerance } & {\tt Atol} \\
\hline
{\tt relative\_tolerance } & {\tt Rtol} \\
\hline
{\tt b} & $\beta$ \\
\hline
{\tt fac\_min} & $f_{\rm min}$ \\
\hline
{\tt fac\_max} & $f_{\rm max}$ \\
\hline
\end{tabular}
\caption{Correspondence between the user \mimes user input and the RK parameters described in the text.}
\label{tab:RK_mimes_params}
\end{table}
\subsection{Explicit embedded RK methods}
%
Explicit methods use only the information of the previous step in order to compute $\vec{k}_i$ from
%
\begin{equation}
\vec{k}_{i}=\vec f\lrBiggb{ \vec{y}_{n}+h \lrBigb{ \sum_{j=1}^{i-1}a_{ij}\vec{k}_{j} }, t_{n}+ c_{i} \; h}\;, \quad \forall i=1,2,\dots, s\;,
\label{eq:explicit_RK_k}
\end{equation}
%
with $a_{ij}$, $c_i$, together with $b_i$ and $b_i^{(*)}$, consist the so-called Butcher tableau of teh corresponding method. For explicit methods this is usually presented as
%
\begin{equation}
\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_{s s-1}& 0 \\
\hline
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}
\label{eq:Butcher}
\end{equation}
%
It should be noted that $c_i = \displaystyle\sum_{j=1}^{s} a_{ij} $.
\subsection{Rosenbrock methods}
%
Explicit RK methods, encounter instabilities when a system is ``stiff"~\footnote{A definition of stiffness can be found in~\cite{hairer2010solving,Hairer}.}; when is oscillates rapidly, or has different elements at different scales. This problems are somewhat resolved by trying to predict the next step inside $\vec{k}$; \ie 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~\ref{eq:Butcher} and linearise the system of non-linear equations (for details, see~\cite{hairer2010solving}). In these methods, $\vec{k}$ is determined by
%
\begin{equation}
\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}\; ,
\label{eq:Ros_k}
\end{equation}
%
which is written in such a way that everything is evaluated at $t = t_{n}$. In \eqs{eq:Ros_k}, $\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^{(*)}$).
\paragraph{Implementing a new Butcher tableau in {\tt NaBBODES}} As already mentioned, \mimes uses {\tt NaBBODES}, which supports the implementation of new Butcher tableaux. This is done by adding a new class (or struct) inside the header file {\tt METHOD.hpp} that can be found in {\tt \mimes/src/NaBBODES/RKF} for the explicit RK and {\tt \mimes/src/NaBBODES/Rosenbrock} for the Rosenbrock embedded methods. All the new parameters must be {\tt public}, {\tt constexpr static} variables, of a type that is the template parameter of the method. For example, the Heun-Euler method can be implemented, adding the following code in {\tt \mimes/src/NaBBODES/RKF/METHOD.hpp}
%
\begin{cpp}
/*-----Implementation of the Heun-Euler method-----*/
/*-----It shouldn't be used for MiMeS, since it is likely to fail-----*/
//LD is the numeric type (\ie double, or long double).
template<class LD>
//use struct because its variables are public by default.
struct HeunEuler{
static constexpr unsigned int s=2; // HeunEuler is a 2-stage RK
static constexpr unsigned int p=1; // order is 1 (its embedded order is p+1=2)
//these are aliases for the arrays
using arr=std::array<LD,s>;
using arr2=std::array<std::array<LD,s>,s>;
static constexpr arr b={0.5,0.5}; // this is $b_i$
static constexpr arr bstar={0.5,0.5}; // this is $b_i^{(*)}$
static constexpr arr c={0,1}; // remember that $c_i = \displaystyle\sum_{j=1}^{s} a_{ij}$
// this is $a_{ij}$
static constexpr arr2 a={
arr{0,0},
arr{1,0}
};
};
\end{cpp}
In order to implement the ROS3w~\cite{RangAngermann2005} method, one can add the following code in the header file {\tt \mimes/src/NaBBODES/Rosenbrock/METHOD.hpp}
%
\begin{cpp}
/*-----Implementation of the ROS3w method-----*/
/*-----It shouldn't be used for MiMeS, since it is likely to fail-----*/
//LD is the numeric type (\ie double, or long double).
template<class LD>
struct ROS3w{
static constexpr unsigned int s=3; // 3-stage
static constexpr unsigned int p=2; // second order (embedded order is $p+1$)
//aliases for the arrays
using arr=std::array<LD,s>;
using arr2=std::array<std::array<LD,s>,s>;
static constexpr arr b={0.25,0.25,0.5 }; // this is $b_i$
// this is $b_i^{(*)}$
static constexpr arr bstar={ 0.746704703274,0.1144064078371,0.1388888888888};
static constexpr arr c={0,2/3.,4/3.}; // remember that $c_i = \displaystyle\sum_{j=1}^{s} a_{ij}$
static constexpr LD gamma=0.4358665215084; // this is $\gamma$
// this is $a_{ij}$
static constexpr arr2 a={
arr{0,0,0},
arr{2/3.,0,0},
arr{2/3.,2/3.,0}
};
// this is $\gamma_{ij}$
static constexpr arr2 g={
arr{0,0,0},
arr{0.3635068368900,0,0},
arr{-0.8996866791992,-0.1537997822626,0}
};
};
\end{cpp}
\paragraph{How to compile \mimes in order to use the newly implemented method} One a new Butcher tableau is implemented, the {\tt mimes::Axion} class can use it. This class, needs the name of method assigned to the macro {\tt METHOD}. That is, before you include the header file {\tt src/Axion/AxionSolve.hpp}, this macro needs to be defined. A convenient way to do this, is to define the macro using the {\tt -D} flag of the compiler; \eg in order to use the ROS3w method, this definition ca be done by adding {\tt -DMETHOD=ROS3w} when compiling the code. Alternatively, if one uses the {\tt makefile} files, a method can be chosen by adding it in the corresponding {\tt Definitions.mk} file as {\tt METHOD=ROS3w}.~\footnote{One make sure to use the correct {\tt SOLVER} macro -- {\tt SOLVER=1} for Rosenbrock method and {\tt SOLVER=2} for explicit RK method-- otherwise the compilation will fail.}
\section{\CPP classes}\label{app:classes}
\setcounter{equation}{0}
\DK{Describe all the classes.}
\DK{Mention that one has access to both the local error, {\tt mimes::Axion<LD,Solver,Method>::System::error[eq][step]}, and $\Delta$ defined in \eqs{eq:RK_scale} {\tt mimes::Axion<LD,Solver,Method>::System::Delta[step]}. }
\section{\PY modules}\label{app:modules}
\setcounter{equation}{0}
\DK{Describe the modules.}
\section{Utilities}\label{app:util}
\setcounter{equation}{0}
\DK{Describe {\tt FormatFile.sh} and {\tt timeit.sh}.}
\section{Quick guide to user input}\label{app:usr_input}
\setcounter{equation}{0}
In table~\ref{tab:input}, we provide a reference for all the user input, and compile-time options.
%
\begin{table}[h!]
\centering
\begin{tabular}{l l}
\hline\\[-0.4cm]
\multicolumn{2}{c}{\bf User run-time input} \\
\hline\\[-0.4cm]
{\tt theta\_i} & initial angle. \\
\hline\\[-0.4cm]
{\tt fa} & the PQ scale.\\
\hline\\[-0.4cm]
{\tt umax } & Once $u=\log a/a_i>${\tt umax}, the integration stops. Typical value: $\sim 500$.\\
\hline\\[-0.4cm]
{\tt TSTOP} & Once $T<${\tt TSTOP}, integration stops. Typical value: $10^{-4}~\GeV$.\\
\hline\\[-0.4cm]
{\tt ratio\_ini}& Integration starts at $u$ with $3H/\maT \approx${\tt ratio\_ini}. Typical value: $\sim 10^{3}$.\\
\hline\\[-0.4cm]
\multirow{1}{4cm}{{\tt N\_convergence\_max} {\tt convergence\_lim}} & \multirow{1}{12cm}{Integration stops when the relative difference
between two consecutive peaks is less than {\tt convergence\_lim} for {\tt N\_convergence\_max}
consecutive peaks. } \\ \\ \\
\hline\\[-0.4cm]
{\tt inputFile} & \multirow{1}{12cm}{Relative (or absolute) path to a file that describes the cosmology. The columns should be: $u$ $T ~[\GeV]$ $\log H$, with acceding $u$. Entropy injection should have stopped before the lowest temperature of given in {\tt inputFile}.} \\ \\ \\ \\
\hline\\[-0.4cm]
{\tt initial\_stepsize} & Initial step-size of the solver. Default value: $10^{-2}$.\\
\hline\\[-0.4cm]
{\tt minimum\_stepsize} & Lower limit of the step-size. Default value: $10^{-8}$.\\
\hline\\[-0.4cm]
{\tt maximum\_stepsize} & Upper limit of the step-size. Default value: $10^{-2}$.\\
\hline\\[-0.4cm]
{\tt absolute\_tolerance} & \multirow{1}{12cm}{Absolute tolerance of the RK solver (see also table~\ref{tab:RK_mimes_params}). Default value: $10^{-8}$.}\\\\
\hline\\[-0.4cm]
{\tt relative\_tolerance} & \multirow{1}{12cm}{Relative tolerance of the RK solver (see also table~\ref{tab:RK_mimes_params}). Default value: $10^{-8}$.}\\\\
\hline\\[-0.4cm]
{\tt beta} & \multirow{1}{12cm}{Aggressiveness of the adaptation strategy (see also table~\ref{tab:RK_mimes_params}). Default value: $0.9$.}\\\\
\hline\\[-0.4cm]
{\tt fac\_max}, {\tt fac\_min} &\multirow{1}{12cm}{The step-size does not change more than {\tt fac\_max} and less than {\tt fac\_min} within a trial step (see also table~\ref{tab:RK_mimes_params}). Default values: $1.2$ and $0.8$, respectively.} \\ \\ \\
\hline\\[-0.4cm]
{\tt maximum\_No\_steps} & \multirow{1}{12cm}{If integration needs more than {\tt maximum\_No\_steps} integration stops. Default value: $10^7$.}\\\\
\hline\\[-0.4cm]
\end{tabular}
\caption{Table of run-time user input.}
\label{tab:run_time-input}
\end{table}
\begin{table}[h!]
\centering
\begin{tabular}{l l}
\multicolumn{2}{c}{\bf Required data files, with corresponding variables in {\tt \mimes/Paths.mk}} \\
\hline\\[-0.4cm]
{\tt cosmoDat}& \multirow{1}{12cm}{Relative path to data file with $T$ (in $\GeV$), $\heff$, $\geff$. If the path changes one must run
{\tt bash configure.sh} and {\tt make}.}\\\\
\hline\\[-0.4cm]
{\tt axMDat}& \multirow{1}{12cm}{Relative path to data file with $T$ (in $\GeV$), $\chi$ (defined from \eqs{eq:axion_mass_def}). If the path changes one must run {\tt bash configure.sh} and {\tt make}. If a function of the axion mass is used, this variable can be omitted.}\\\\\\
\hline\\[-0.4cm]
{\tt anFDat}& \multirow{1}{12cm}{Relative path to data file with $\thetamax$, $f(\thetamax)$. If the path changes one must run
{\tt bash configure.sh} and {\tt make}.}\\\\
\hline\\[-0.4cm]
\end{tabular}
\caption{Paths to the required data files.}
\label{tab:input}
\end{table}
\begin{table}[h!]
\centering
\begin{tabular}{l l}
\multicolumn{2}{c}{\bf User compile-time input} \\
\hline\\[-0.4cm]
{\tt rootDir}& \multirow{1}{12cm}{The relative path of root directory of \mimes. Relevant only when compiling using {\tt make}. Available in all {\tt Definitions.mk}.}\\\\
\hline\\[-0.4cm]
{\tt LONG}& \multirow{1}{12cm}{{\tt long} for {\tt long double }or empty for {\tt double}. This should be a macro at the top of the {\tt .cpp} file. Available in all {\tt Definitions.mk}.}\\\\
\hline\\[-0.4cm]
{\tt LONGpy}& \multirow{1}{12cm}{{\tt long} or empty. Same as {\tt LONG}, applies in the python modules. Available in {\tt \mimes/Definitions.mk}.}\\\\
\hline\\[-0.4cm]
{\tt SOLVER}& \multirow{1}{12cm}{In order to use a Rosenbrock method {\tt SOLVER}=$1$. For explicit RK method, {\tt SOLVER}=$1$.
This should be a macro defined before we include {\tt \mimes/src/Axion/AxionSolve.hpp}. Therefore it is preferable to
define it with the compiler flag {\tt -DSOLVER=}$1$ or $2$. The corresponding variable in {\tt \mimes/Definitions.mk} applies to the python modules. The variable in {\tt \mimes/UserSpace/Cpp/Axion/Definitions.mk} applies to the example in the same directory.}\\\\\\\\\\\\\\
\hline\\[-0.4cm]
{\tt METHOD}& \multirow{1}{12cm}{Depending of the solver, this macro should name one of its available methods.
For {\tt SOLVER}=$1$, {\tt METHOD}=RODASPR2 (fourth order) or ROS34PW2 (third order).
For {\tt SOLVER}=$1$, {\tt METHOD}=DormandPrince (seventh order). This should be a macro defined before we include {\tt \mimes/src/Axion/AxionSolve.hpp}. Therefore it is preferable to
define it with the compiler flag {\tt -DMETHOD=}RODASPR2, ROS34PW2, or DormandPrince.
The corresponding variable in {\tt \mimes/Definitions.mk} applies to the python modules. The variable in {\tt \mimes/UserSpace/Cpp/Axion/Definitions.mk} applies to the example in the same directory.}\\\\\\\\\\\\\\\\\\\\
\hline\\[-0.4cm]
\multicolumn{2}{c}{\bf Compiler options} \\
\hline\\[-0.4cm]
{\tt CC} & \multirow{1}{12cm}{The preferred \CPP compiler ({\tt g++} by default). Corresponding variable in all {\tt Definitions.mk} files.} \\\\
\hline\\[-0.4cm]
{\tt OPT} & \multirow{1}{12cm}{Optimization level of the compiler. Available options {\tt OPT}={\tt O1}, {\tt O2}, {\tt O3} (be default). By Corresponding variable in all {\tt Definitions.mk} files.} \\\\
\hline\\[-0.4cm]
\end{tabular}
\caption{User compile-time input.}
\label{tab:compile_time-input}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\bibliography{refs}{}
\bibliographystyle{JHEP}
\end{document}
diff --git a/UserSpace/Python/Axion.py b/UserSpace/Python/Axion.py
index 2377bf2..45ff781 100644
--- a/UserSpace/Python/Axion.py
+++ b/UserSpace/Python/Axion.py
@@ -1,193 +1,193 @@
from time import time
from sys import stderr
from sys import path as sysPath
from os import path as osPath
sysPath.append(osPath.join(osPath.dirname(__file__), '../../src'))
#load the module
from interfacePy.Axion import Axion
from interfacePy.AxionMass import AxionMass
from interfacePy.Cosmo import Hubble,mP
theta_i, fa=0.94435, 1e12
umax=500
TSTOP=1e-4
ratio_ini=1e3
N_convergence_max, convergence_lim=5, 1e-2 #this is fine, but you can experiment a bit.
#radiation dominated example
inputFile="../InputExamples/RDinput.dat"
# Matter Domination example.
# the NSC parameters (using the notation of 2012.07202) are:
# T_end=1e-2 (GeV), c=3, T_ini=1e12 (GeV), and r=1e-1
# inputFile="../InputExamples/MatterInput.dat"
# Kination Domination example.
# the NSC parameters (using the notation of 2012.07202) are:
# T_end=0, c=6, T_ini=1e3 (GeV), and r=1e10
# inputFile="../InputExamples/KinationInput.dat"
# options for the solver
# These variables are optional. Yoou can use the Axion class without them.
initial_step_size=1e-2; #initial step the solver takes.
minimum_step_size=1e-8; #This limits the sepsize to an upper limit.
maximum_step_size=1e-2; #This limits the sepsize to a lower limit.
absolute_tolerance=1e-8; #absolute tolerance of the RK solver
relative_tolerance=1e-8; #relative tolerance of the RK solver
beta=0.9; #controls how agreesive the adaptation is. Generally, it should be around but less than 1.
#The stepsize does not increase more than fac_max, and less than fac_min.
#This ensures a better stability. Ideally, fac_max=inf and fac_min=0, but in reality one must
#tweak them in order to avoid instabilities.
fac_max=1.2;
fac_min=0.8;
maximum_No_steps=int(1e7); #maximum steps the solver can take Quits if this number is reached even if integration is not finished.
_=time()
# AxionMass instance
-axionMass = AxionMass(r'../../src/data/chi.dat',0,mP)
-
-
-#------------------------------------------------------------------------------#
-# this is the axion mass squared beyond the interpolation limits for the current data
-# if yo don't specify it, the axion mass is taken to be constant beyond these limits
-TMax=axionMass.getTMax()
-chiMax=axionMass.getChiMax()
-TMin=axionMass.getTMin()
-chiMin=axionMass.getChiMin()
-
-axionMass.set_ma2_MAX( lambda T,fa: chiMax/fa/fa*pow(T/TMax,-8.16) )
-axionMass.set_ma2_MIN( lambda T,fa: chiMin/fa/fa )
-#------------------------------------------------------------------------------#
-# def ma2(T,fa):
- # TQCD=150*1e-3;
- # ma20=3.1575e-05/fa/fa;
- # if T<=TQCD:
- # return ma20;
- # return ma20*pow((TQCD/T),8.16)
-# axionMass = AxionMass(ma2)
+# axionMass = AxionMass(r'../../src/data/chi.dat',0,0.6)
+# #------------------------------------------------------------------------------#
+# # this is the axion mass squared beyond the interpolation limits for the current data
+# # if yo don't specify it, the axion mass is taken to be constant beyond these limits
+# TMax=axionMass.getTMax()
+# chiMax=axionMass.getChiMax()
+# TMin=axionMass.getTMin()
+# chiMin=axionMass.getChiMin()
+
+# axionMass.set_ma2_MAX( lambda T,fa: chiMax/fa/fa*pow(T/TMax,-8.16) )
+# axionMass.set_ma2_MIN( lambda T,fa: chiMin/fa/fa )
+# #------------------------------------------------------------------------------#
+def ma2(T,fa):
+ TQCD=150*1e-3
+ ma20=3.1575e-05/fa/fa
+ if T<=TQCD:
+ return ma20
+ return ma20*pow((TQCD/T),8.16)
+axionMass = AxionMass(ma2)
# Axion instance
ax=Axion(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max, convergence_lim, inputFile, axionMass,
initial_step_size,minimum_step_size, maximum_step_size, absolute_tolerance,
relative_tolerance, beta, fac_max, fac_min, maximum_No_steps)
# solve the EOM (this only gives you the relic, T_osc, theta_osc, and a_osc)
ax.solveAxion()
print(ax.theta_i, ax.fa, ax.theta_osc, ax.T_osc ,ax.relic)
print(round(time()-_,3),file=stderr)
+
+
# change to True in order to mkae the plots
if False:
import matplotlib.pyplot as plt
from numpy import array as np_array
from numpy import abs as np_abs
ax.getPeaks()#this gives you the peaks of the oscillation
ax.getPoints()#this gives you all the points of integration
ax.getErrors()#this gives you local errors of integration
#########-----\theta-----#########
fig=plt.figure(figsize=(9,4))
fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.9, right=0.9,wspace=0.0,hspace=0.25)
sub = fig.add_subplot(1,1,1)
#this plot shows the peaks of the oscillation
sub.plot(ax.T_peak,ax.theta_peak,linestyle=':',marker='+',color='xkcd:blue',linewidth=2)
#this plot shows all the points
sub.plot(ax.T,ax.theta,linestyle='-',linewidth=2,alpha=1,c='xkcd:black')
sub.set_yscale('linear')
sub.set_xscale('linear')
sub.set_xlabel(r'$T ~[{\rm GeV}]$')
sub.xaxis.set_label_coords(0.5, -0.1)
sub.set_ylabel(r'$\theta$')
sub.yaxis.set_label_coords(-0.1,0.5)
sub.axhline(ax.theta_osc,linestyle=':',color='xkcd:red',linewidth=1.5)
sub.axvline(ax.T_osc,linestyle='--',color='xkcd:gray',linewidth=1.5)
fig.savefig('theta-T_examplePlot.pdf',bbox_inches='tight')
#########-----\zeta-----#########
fig=plt.figure(figsize=(9,4))
fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.9, right=0.9,wspace=0.0,hspace=0.25)
sub = fig.add_subplot(1,1,1)
sub.plot(ax.T,ax.zeta,linestyle='-',linewidth=2,alpha=1,c='xkcd:black')
sub.plot(ax.T_peak,ax.zeta_peak,linestyle=':',marker='+',color='xkcd:blue',linewidth=2)
sub.set_yscale('linear')
sub.set_xscale('linear')
sub.set_xlabel(r'$T ~[{\rm GeV}]$')
sub.xaxis.set_label_coords(0.5, -0.1)
sub.set_ylabel(r'$\zeta $')
sub.yaxis.set_label_coords(-0.1,0.5)
sub.axvline(ax.T_osc,linestyle='--',color='xkcd:gray',linewidth=1.5)
fig.savefig('zeta-T_examplePlot.pdf',bbox_inches='tight')
#########-----errors-----#########
fig=plt.figure(figsize=(9,4))
fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.9, right=0.9,wspace=0.0,hspace=0.25)
sub = fig.add_subplot(1,1,1)
sub.plot(ax.T,np_abs(ax.dtheta/ax.theta),linestyle='-',linewidth=2,alpha=1,c='xkcd:black',label=r'$\dfrac{\delta \theta}{\theta}$')
sub.plot(ax.T,np_abs(ax.dzeta/ax.zeta),linestyle='-',linewidth=2,alpha=1,c='xkcd:red',label=r'$\dfrac{\delta \zeta}{\zeta}$')
sub.set_yscale('log')
sub.set_xscale('linear')
sub.set_xlabel(r'$T ~[{\rm GeV}]$')
sub.xaxis.set_label_coords(0.5, -0.1)
sub.set_ylabel(r'local errors')
sub.yaxis.set_label_coords(-0.1,0.5)
sub.legend(bbox_to_anchor=(1, 0.0),borderaxespad=0.,
borderpad=0.05,ncol=1,loc='lower right',fontsize=14,framealpha=0)
sub.axvline(ax.T_osc,linestyle='--',color='xkcd:gray',linewidth=1.5)
fig.savefig('errors_examplePlot.pdf',bbox_inches='tight')
#########-----rho_a-----#########
fig=plt.figure(figsize=(9,4))
fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.9, right=0.9,wspace=0.0,hspace=0.25)
sub = fig.add_subplot(1,1,1)
sub.plot(ax.T,ax.rho_axion,linestyle='-',linewidth=2,alpha=1,c='xkcd:black')
sub.plot(ax.T_peak,ax.rho_axion_peak,linestyle=':',linewidth=2,alpha=1,c='xkcd:blue')
sub.set_yscale('log')
sub.set_xscale('linear')
sub.set_xlabel(r'$T ~[{\rm GeV}]$')
sub.xaxis.set_label_coords(0.5, -0.1)
sub.set_ylabel(r'$\rho_{a}(T) ~[{\rm GeV}^{4}]$')
sub.yaxis.set_label_coords(-0.1,0.5)
sub.axvline(ax.T_osc,linestyle='--',color='xkcd:gray',linewidth=1.5)
fig.savefig('rho_a-T_examplePlot.pdf',bbox_inches='tight')
#don't forget the destructor
del ax
\ No newline at end of file
diff --git a/src/Axion/Axion-py.cpp b/src/Axion/Axion-py.cpp
index 11b25e5..ebcac93 100644
--- a/src/Axion/Axion-py.cpp
+++ b/src/Axion/Axion-py.cpp
@@ -1,88 +1,89 @@
#include<string>
#include "src/Axion/AxionSolve.hpp"
#include"src/AxionMass/AxionMass.hpp"
// macros for the solver
#ifndef SOLVER
#define SOLVER 1
#define METHOD RODASPR2
#endif
// macros for the numeric type
#ifndef LONGpy
#define LONGpy
#endif
#ifndef LD
#define LD LONGpy double
#endif
// use this to cast void* to Axion*
#define Cast(ax) static_cast<mimes::Axion<LD,SOLVER,METHOD<LD>>*>(ax)
extern "C"{
//constructor
void* INIT(LD theta_i, LD fa, LD umax, LD TSTOP, LD ratio_ini, unsigned int N_convergence_max, LD convergence_lim, char* inputFile,
void *axionMass,
LD initial_step_size, LD minimum_step_size, LD maximum_step_size,
LD absolute_tolerance, LD relative_tolerance, LD beta, LD fac_max,
LD fac_min, unsigned int maximum_No_steps){
return new mimes::Axion<LD,SOLVER,METHOD<LD>>(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max,convergence_lim, static_cast<std::string>(inputFile),
static_cast<mimes::AxionMass<LD>*>(axionMass),
initial_step_size,minimum_step_size, maximum_step_size, absolute_tolerance, relative_tolerance, beta,
fac_max, fac_min, maximum_No_steps);
}
// destructor to delete the void*
void DEL(void *Ax){ delete static_cast<mimes::Axion<LD,SOLVER,METHOD<LD>>*>(Ax); Ax=nullptr; }
// set theta_i
void setTheta_i(LD theta_i, void * Ax){ Cast(Ax) -> setTheta_i(theta_i); }
void MAKE(void * Ax){ Cast(Ax) -> solveAxion(); }
+ void restart(void * Ax){ Cast(Ax) -> restart(); }
unsigned int getPointSize(void * Ax){ return Cast(Ax) -> pointSize; }
unsigned int getPeakSize(void * Ax){ return Cast(Ax) -> peakSize; }
void getResults(LD *results, void *Ax){
results[0]=Cast(Ax) ->T_osc;
results[1]=Cast(Ax)->a_osc;
results[2]=Cast(Ax) ->theta_osc;
results[3]=Cast(Ax) ->relic;
}
void getPoints(LD *a, LD *T, LD *theta, LD *zeta, LD *rho_a, void *Ax){
unsigned int N = Cast(Ax) ->pointSize;
for (unsigned int i = 0; i < N; i++){
a[i]=Cast(Ax) ->points[i][0];
T[i]=Cast(Ax) ->points[i][1];
theta[i]=Cast(Ax) ->points[i][2];
zeta[i]=Cast(Ax) ->points[i][3];
rho_a[i]=Cast(Ax) ->points[i][4];
}
}
void getErrors(LD *dtheta, LD *dzeta, void *Ax){
unsigned int N = Cast(Ax) ->pointSize;
for (unsigned int i = 0; i < N; i++){
dtheta[i]=Cast(Ax) ->dtheta[i];
dzeta[i]=Cast(Ax) ->dzeta[i];
}
}
void getPeaks(LD *a_peak, LD *T_peak, LD *theta_peak, LD*zeta_peak, LD *rho_a_peak, LD *ad_inv_peak, void *Ax){
unsigned int Npeaks = Cast(Ax) ->peakSize;
for (unsigned int i = 0; i < Npeaks; i++){
a_peak[i]=Cast(Ax) ->peaks[i][0];
T_peak[i]=Cast(Ax) ->peaks[i][1];
theta_peak[i]=Cast(Ax) ->peaks[i][2];
zeta_peak[i]=Cast(Ax) ->peaks[i][3];
rho_a_peak[i]=Cast(Ax) ->peaks[i][4];
ad_inv_peak[i]=Cast(Ax) ->peaks[i][5];
}
}
}
diff --git a/src/Axion/AxionEOM.hpp b/src/Axion/AxionEOM.hpp
index 2e45008..397b27d 100644
--- a/src/Axion/AxionEOM.hpp
+++ b/src/Axion/AxionEOM.hpp
@@ -1,167 +1,171 @@
#ifndef SYSTEM_AxionEOM
#define SYSTEM_AxionEOM
#include <iostream>
#include <fstream>
#include <cmath>
#include <array>
#include <string>
#include "src/SimpleSplines/Interpolation.hpp"
#include"src/AxionMass/AxionMass.hpp"
namespace mimes{
constexpr unsigned int Neqs=2;
template<class LD> using Array = std::array<LD,Neqs>;
template<class LD>
class AxionEOM{
public:
LD fa,ratio_ini;
LD T_stop, logH2_stop, u_stop;//temperature, logH^2, and u where we stop interpolation (the end of the file). They should be AFTER entropy conservation resumes (I usually stop any intergation at T<1 MeV or so, where the Universe expands in a standard way)!
LD T_ini, logH2_ini;//temperature and logH^2 (t=0 by definition) where we start interpolation (the end of the file). They should be AFTER entropy conservation resumes (I usually stop any intergation at T<1 MeV or so, where the Universe expands in a standard way)!
LD u_osc, T_osc;//temperature and logH^2 (t=0 by definition) where we start interpolation (the end of the file). They should be AFTER entropy conservation resumes (I usually stop any intergation at T<1 MeV or so, where the Universe expands in a standard way)!
std::vector<LD> u_tab,T_tab,logH2_tab;//these will store the data from inputFile
CubicSpline<LD> T_int; //interpolation of the temperature
CubicSpline<LD> logH2_int; //interpolation of logH^2
AxionMass<LD> *axionMass; //pointer to an instance of AxionMass
AxionEOM()=default;
~AxionEOM()=default;
// constructor of AxionEOM.
/*
fa: PQ scale in GeV (the temperature dependent mass is defined as m_a^2(T) = \chi(T)/f^2)
ratio_ini: interpolations start when 3H/m_a<~ratio_ini
inputFile: file that describes the cosmology. the columns should be: u T[GeV] logH
axionMass: pointer to an instance of AxionMass
*/
AxionEOM(LD fa, LD ratio_ini, std::string inputFile, AxionMass<LD> *axionMass){
this->fa=fa;
this->ratio_ini=ratio_ini;
this->axionMass=axionMass;
LD u=0,T=0,logH=0;//current line in file
LD u_prev=0,T_prev=0,logH_prev=0;//previous line in file
+ /*clear the data just in case*/
+ u_tab.clear();
+ T_tab.clear();
+ logH2_tab.clear();
std::ifstream data_file(inputFile,std::ios::in);
//quit if the file does not exist.
if(not data_file.good()){
std::cerr << inputFile<< " does not exist.";
std::cerr<<" Please make sure to provide a valid path.\n";
exit(1);
}
bool ini_check=true; //check when ratio_ini is reached
bool osc_check=true; //check when T_osc is reached
LD u_ini=0; //check when ratio_ini is reached in order to rescale u to start at 0 in the interpolations
LD ratio=0;// I will use ratio to store 3H/m_a as I read the file
//read the file and fill u_tab, T_tab, and logH2_tab; and find T_osc.
unsigned int N=0;
while (not data_file.eof()){
data_file>>u;
data_file>>T;
data_file>>logH;
//if there is an empty line, u does not change, so do skip this line :).
if(N>1 and u==u_prev){continue;}
ratio = 3*std::exp(logH) / std::sqrt(axionMass->ma2(T,fa));
if(ratio <= ratio_ini ){
//we need the following check in order to find u_ini (it is better to start at the point before ratio_ini is reached)
if(ini_check){
T_ini=T_prev;
logH2_ini=2*logH_prev;
u_ini=u_prev;
u_tab.push_back(0);
T_tab.push_back(T_prev);
logH2_tab.push_back(2*logH_prev);
ini_check=false;
}
u_tab.push_back(u - u_ini);
T_tab.push_back(T);
logH2_tab.push_back(2*logH);
}
//find T and u when oscillation begins.
if(ratio<=1){if(osc_check){osc_check=false; u_osc=u_prev-u_ini;T_osc=T_prev;}}
u_prev=u;
T_prev=T;
logH_prev=logH;
N++;
}
u_stop=u - u_ini;
T_stop=T;
logH2_stop=2*logH;
data_file.close();
};
// make the interpolations
void makeInt(){
T_int=CubicSpline<LD>{&u_tab,&T_tab};
logH2_int=CubicSpline<LD>{&u_tab,&logH2_tab};
}
//the temperature at t. in order to avoid getting out of bounds, we take T to be constant for
//t<0 and t>t_stop (not the best idea, but it works with the solver)
LD Temperature(LD u){
if(u<=0){return T_ini;}
if(u>=u_stop){return T_stop;}
return T_int(u);
}
//logH^2 at u. in order to avoid getting out of bounds, we take logH2 to be constant for
//u<0 and u>u_stop (not the best idea, but it works with the solver)
LD logH2(LD u){
if(u<=0){return logH2_ini;}
if(u>=u_stop){return logH2_stop;}
return logH2_int(u);
}
//the temperature at u. in order to avoid getting out of bounds, we take it to be 0 for
//t<0 and t>t_stop (not the best idea, but it works with the solver)
LD dlogH2du(LD u){
if(u<=0){return 0;}
if(u>=u_stop){return 0;}
return logH2_int.derivative_1(u);
}
//overload operator() in order to call it from the solver
void operator()(Array<LD> &lhs, Array<LD> &y, LD u){
//remember that t=log{a/a_i}
LD theta=y[0];
LD zeta=y[1];
// you need to load a file that contains u, T, logH, and interpolate T and logH^2
LD _T= Temperature(u);//temperature
LD _H2= std::exp(logH2(u)) ;// e^{log(H^2)}=H^2
LD _dlogH2du=dlogH2du(u);// dlogH^2/du
// in radiation dominated Universe, you can also use:
// LD _dlogH2du=-(cosmo.dgeffdT(_T)*_T/cosmo.geff(_T) + 4 )/cosmo.dh(_T);
lhs[0]=zeta;//dtheta/du
lhs[1]=-(0.5*_dlogH2du +3)*zeta - axionMass->ma2(_T,fa)/_H2*std::sin(theta);//dzeta/du
};
};
};
#endif
\ No newline at end of file
diff --git a/src/Axion/AxionSolve.hpp b/src/Axion/AxionSolve.hpp
index 8d728b3..9518add 100644
--- a/src/Axion/AxionSolve.hpp
+++ b/src/Axion/AxionSolve.hpp
@@ -1,386 +1,398 @@
#ifndef SolveAxion_included
#define SolveAxion_included
#include <cmath>
#include <string>
#include <vector>
//----Solver-----//
#include "src/NaBBODES/RKF/RKF_class.hpp"
#include "src/NaBBODES/RKF/RKF_costructor.hpp"
#include "src/NaBBODES/RKF/RKF_calc_k.hpp"
#include "src/NaBBODES/RKF/RKF_sums.hpp"
#include "src/NaBBODES/RKF/RKF_step_control-simple.hpp"
// #include "src/NaBBODES/RKF/RKF_step_control-PI.hpp"
#include "src/NaBBODES/RKF/RKF_steps.hpp"
#include "src/NaBBODES/RKF/METHOD.hpp"
//----Solver-----//
#include "src/NaBBODES/Rosenbrock/LU/LU.hpp"
#include "src/NaBBODES/Rosenbrock/Ros_class.hpp"
#include "src/NaBBODES/Rosenbrock/Ros_costructor.hpp"
#include "src/NaBBODES/Rosenbrock/Ros_LU.hpp"
#include "src/NaBBODES/Rosenbrock/Ros_calc_k.hpp"
#include "src/NaBBODES/Rosenbrock/Ros_sums.hpp"
#include "src/NaBBODES/Rosenbrock/Ros_step_control-simple.hpp"
// #include "src/NaBBODES/Rosenbrock/Ros_step_control-PI.hpp"
#include "src/NaBBODES/Rosenbrock/Ros_steps.hpp"
#include "src/NaBBODES/Rosenbrock/Jacobian/Jacobian.hpp"
#include "src/NaBBODES/Rosenbrock/METHOD.hpp"
//---Get what you need for the axion--//
#include"src/misc_dir/path.hpp"
#include"src/Axion/AxionEOM.hpp"
#include"src/AxionMass/AxionMass.hpp"
#include"src/Cosmo/Cosmo.hpp"
#include"src/AnharmonicFactor/AnharmonicFactor.hpp"
/*================================*/
// If statement, in order to use the desired solver
namespace mimes{
template<bool C, typename T1, typename T2>
struct IF{using type=T1;};
template<typename T1, typename T2>
struct IF<false,T1,T2>{using type=T2;};
};
namespace mimes{
template<class LD, const int Solver, class Method>
class Axion{
using RKSolver = typename IF<Solver==1,
Ros<Neqs, Method, Jacobian<Neqs, LD>, LD>,
typename IF<Solver==2,RKF<Neqs, Method, LD>,void>::type>::type;
LD umax,TSTOP,ratio_ini;
unsigned int N_convergence_max;
LD convergence_lim;
std::string inputFile;
- AxionEOM<LD> axionEOM;
LD initial_step_size, minimum_step_size, maximum_step_size, absolute_tolerance, relative_tolerance;
LD beta,fac_max,fac_min;
unsigned int maximum_No_steps;
constexpr static LD theta_small=1e-3;
LD theta_ini; //use this to rescale theta_i is it is less than theta_small
+ /*internal axionEOM variable*/
+ AxionEOM<LD> axionEOM;
+ /*pointer to AxionMass instance*/
+ AxionMass<LD> *axionMass;
public:
LD theta_i,fa;
LD theta_osc, T_osc, a_osc;
LD gamma;//gamma is the entropy injection from the point of the last peak until T_stop (the point where interpolation stops)
LD relic;
std::vector<std::vector<LD>> points;
std::vector<std::vector<LD>> peaks;
std::vector<LD> dtheta;
std::vector<LD> dzeta;
unsigned int pointSize,peakSize;
static AnharmonicFactor<LD> anharmonicFactor;
static Cosmo<LD> plasma;
// constructor of Axion.
/*
theta_i: initial angle
fa: PQ scale in GeV (the temperature dependent mass is defined as m_a^2(T) = \chi(T)/f^2)
umax: if u>umax the integration stops (rempember that u=log(a/a_i))
TSTOP: if the temperature drops below this, integration stops
ratio_ini: integration starts when 3H/m_a<~ratio_ini (this is passed to AxionEOM,
in order to make the interpolations start at this point)
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
inputFile: file that describes the plasmalogy. the columns should be: u T[GeV] logH
axionMass: pointer to an instance of AxionMass
-----------Optional arguments------------------------
initial_stepsize: initial step the solver takes.
maximum_stepsize: This limits the sepsize to an upper limit.
minimum_stepsize: This limits the sepsize to a lower limit.
absolute_tolerance: absolute tolerance of the RK solver
relative_tolerance: relative tolerance of the RK solver
Note:
Generally, both absolute and relative tolerances should be 1e-8.
In some cases, however, one may need more accurate result (eg if f_a is extremely high,
the oscillations happen violently, and the ODE destabilizes). Whatever the case, if the
tolerances are below 1e-8, long doubles *must* be used.
beta: controls how agreesive the adaptation is. Generally, it should be around but less than 1.
fac_max, fac_min: the stepsize does not increase more than fac_max, and less than fac_min.
This ensures a better stability. Ideally, fac_max=inf and fac_min=0, but in reality one must
tweak them in order to avoid instabilities.
maximum_No_steps: maximum steps the solver can take Quits if this number is reached even if integration
is not finished.
*/
Axion(LD theta_i, LD fa, LD umax, LD TSTOP, LD ratio_ini,
unsigned int N_convergence_max, LD convergence_lim, std::string inputFile,
AxionMass<LD> *axionMass,
LD initial_step_size=1e-2, LD minimum_step_size=1e-8, LD maximum_step_size=1e-2,
LD absolute_tolerance=1e-8, LD relative_tolerance=1e-8,
LD beta=0.9, LD fac_max=1.2, LD fac_min=0.8, unsigned int maximum_No_steps=10000000){
this->theta_i=theta_i;
this->fa=fa;
this->umax=umax;
this->TSTOP=TSTOP;
this->ratio_ini=ratio_ini;
this->N_convergence_max=N_convergence_max;
this->convergence_lim=convergence_lim;
this->inputFile=inputFile;
this->initial_step_size=initial_step_size;
this->minimum_step_size=minimum_step_size;
this->maximum_step_size=maximum_step_size;
this->absolute_tolerance=absolute_tolerance;
this->relative_tolerance=relative_tolerance;
this->beta=beta;
this->fac_max=fac_max;
this->fac_min=fac_min;
this->maximum_No_steps=maximum_No_steps;
- axionEOM = AxionEOM<LD>(fa, ratio_ini, inputFile, axionMass);
+ this->axionMass=axionMass;
+ axionEOM = AxionEOM<LD>(fa, ratio_ini, inputFile, this->axionMass);
axionEOM.makeInt();//make the interpolations of u,T,logH from inputFile
}
- Axion(){}
+ Axion()=default;
+ ~Axion()=default;
void solveAxion();
+
//use setTheta to set another initial condition
// Warning: this resets public variables (except fa)!
void setTheta_i(LD theta_i){
this->theta_i=theta_i;
- this->theta_osc=0;
points.clear();
peaks.clear();
dtheta.clear();
dzeta.clear();
pointSize=points.size();
peakSize=peaks.size();
theta_osc=theta_i;
gamma=1;
relic=0;
+ }
+ /*reset all variables and remake the interpolations of axionEOM*/
+ void restart(){
+ setTheta_i(theta_i);
+ axionEOM = AxionEOM<LD>(fa, ratio_ini, inputFile, axionMass);
+ axionEOM.makeInt();//make the interpolations of u,T,logH from inputFile
}
+
};
template<class LD, const int Solver, class Method>
AnharmonicFactor<LD> Axion<LD,Solver,Method>::anharmonicFactor(anharmonic_PATH);
template<class LD, const int Solver, class Method>
Cosmo<LD> Axion<LD,Solver,Method>::plasma(cosmo_PATH,0,mimes::Cosmo<LD>::mP);
template<class LD, const int Solver, class Method>
void Axion<LD,Solver,Method>::solveAxion(){
//when theta_i<<1, we can refactor the eom so that it becomes independent from theta_i.
// So, we can solve for theta_i=1e-3, and rescale the result to our desired initial theta.
// This helps avoid any roundoff errors when the amplitude of the oscillation becomes very small.
theta_ini=theta_i;
if(theta_ini<1e-3){theta_ini=1e-3;}
/*================================*/
Array<LD> y0={theta_ini, 0.};//initial conditions
/*================================*/
// instance of the solver
RKSolver System(axionEOM, y0, umax,
initial_step_size, minimum_step_size, maximum_step_size, maximum_No_steps,
absolute_tolerance, relative_tolerance, beta, fac_max,fac_min);
//You find these as you load the data from inputFile
//(it is done automatically in the constructor of axionEOM)
T_osc=axionEOM.T_osc;
a_osc=std::exp(axionEOM.u_osc);
// these parameters are helpful..
unsigned int current_step=0;//count the steps the solver takes
//check is used to identify the peaks, osc_check is used to find the oscillation temperature
bool check=true, osc_check=true;
//use these to count the peaks in order to apply the stopping conditions
unsigned int Npeaks=0, N_convergence=0;
//these variables will be used to fill points (the points the solver takes)
LD theta=0,zeta=0,u=0,a=0,T=0,H2=0,ma2=0;
LD rho_axion=0;
// we will need these for the peaks
LD zeta_prev=0,u_prev=0,theta_prev=0;
LD theta_peak=0,zeta_peak=0,u_peak=0,a_peak=0,T_peak=0,ma2_peak=0;
LD adInv_peak=0,rho_axion_peak=0;
LD an_diff=0;
// this is the initial step (ie points[0])
theta=y0[0]/theta_ini*theta_i;
zeta=y0[1]/theta_ini*theta_i;
T=axionEOM.Temperature(0);
H2=std::exp(axionEOM.logH2(0));
- ma2=axionEOM.axionMass->ma2(T,fa);
+ ma2=axionMass->ma2(T,fa);
if(std::abs(theta)<1e-3){rho_axion=fa*fa*(ma2*0.5*theta*theta);}
else{rho_axion=fa*fa*(ma2*(1 - std::cos(theta)));}
points.push_back(std::vector<LD>{1,T,theta,0,rho_axion});
dtheta.push_back(0);
dzeta.push_back(0);
// the solver identifies the peaks in theta by finding points where zeta goes from positive to negative.
// Once the peaks are identified, we can calculate the adiabatic invariant.
// The solver stops when the idiabatic invariant is almost constant.
while (true){
current_step++;// incease number of steps that the solver takes
//stop if you exceed umax or if the solver takes more than maximum_No_steps, stop
if(System.tn>=System.tmax or current_step==maximum_No_steps){break;}
//take the next step
System.next_step();
// store the local error
dtheta.push_back(System.ynext[0] - System.ynext_star[0]);
dzeta.push_back(System.ynext[1] - System.ynext_star[1]);
//update y (y[0]=theta, y[1]=zeta)
for (unsigned int eq = 0; eq < Neqs; eq++){
System.yprev[eq]=System.ynext[eq];
}
// increase tn
System.tn+=System.h;
//rescale theta and zeta which will be stored in points
//(we rescale them if theta_i<1e-3 in order to avoid roundoff errors)
theta=System.ynext[0]/theta_ini*theta_i;
zeta=System.ynext[1]/theta_ini*theta_i;
//remember that u=log(a/a_i)
u=System.tn;
a=std::exp(u);
//get the temperature which corresponds to u
T=axionEOM.Temperature(u);
//get the Hubble parameter squared which corresponds to u
H2=std::exp(axionEOM.logH2(u));
// If you use as T_osc the one provided by axionEOM, you will not have theta_osc.
// So use this (it doesn't really matter, as T_osc is not a very well defined parameter)
// We could use interpolation to get more accurate T_osc and theta_osc, but it's not worth it,
// because the oscillation temperature does not affect the results.
if(T<=T_osc and osc_check){
T_osc=T;
a_osc=a;
theta_osc=theta;
osc_check=false;
}
//if the temperature is below TSTOP, stop the integration
if(T<TSTOP){break;}
//axion mass squared
- ma2=axionEOM.axionMass->ma2(T,fa);
+ ma2=axionMass->ma2(T,fa);
// If theta<~1e-8, we have roundoff errors due to cos(theta)
// The solution is this (use theta<1e-3; it doesn't matter):
if(std::abs(theta)<1e-3){rho_axion=fa*fa*(0.5* H2*zeta*zeta+ ma2*0.5*theta*theta);}
else{rho_axion=fa*fa*(0.5* H2*zeta*zeta+ ma2*(1 - std::cos(theta)));}
//store current step in points
points.push_back(std::vector<LD>{a,T,theta,zeta,rho_axion});
//check if current point is a peak
if(zeta>0){check=false;}
if(zeta<=0 and check==false){
//increase the total number of peaks
Npeaks++;
//set check=true (this resets the check until the next peak)
check=true;
// use linear interpolation to find u at zeta=0
u_prev=std::log(points[current_step-1][0]);
theta_prev=points[current_step-1][2];
zeta_prev=points[current_step-1][3];
// these are all the quantities at the peak
u_peak=(-zeta_prev*u+zeta*u_prev)/(zeta-zeta_prev);
a_peak=std::exp(u_peak);
theta_peak=((theta-theta_prev)*u_peak+(theta_prev*u-theta*u_prev))/(u-u_prev);
zeta_peak=((zeta-zeta_prev)*u_peak+(zeta_prev*u-zeta*u_prev))/(u-u_prev);
T_peak=axionEOM.Temperature(u_prev);
- ma2_peak=axionEOM.axionMass->ma2(T_peak,fa);
+ ma2_peak=axionMass->ma2(T_peak,fa);
//compute the adiabatic invariant
adInv_peak=anharmonicFactor(theta_peak)*theta_peak*theta_peak *std::sqrt(ma2_peak) * a_peak*a_peak*a_peak ;
if(std::abs(theta)<1e-3){rho_axion_peak=fa*fa*(ma2_peak*0.5*theta_peak*theta_peak);}
else{rho_axion_peak=fa*fa*(ma2_peak*(1 - std::cos(theta_peak)));}
peaks.push_back(std::vector<LD>{a_peak,T_peak,theta_peak,zeta_peak,rho_axion_peak,adInv_peak});
// if the total number of peaks is greater than 2, then you can check for convergence
if(Npeaks>=2){
//difference of adiabatic invariant between current and previous peak
an_diff=std::abs(adInv_peak-peaks[Npeaks-2][5]);
if(adInv_peak>peaks[Npeaks-2][5]){
an_diff=an_diff/adInv_peak;
}else{
an_diff=an_diff/peaks[Npeaks-2][5];
}
//if the difference is smaller than convergence_lim increase N_convergence
//else, reset N_convergence (we stop if the difference is small between consecutive steps )
if(an_diff<convergence_lim){N_convergence++;}else{N_convergence=0;}
}
}
//if N_convergence>=N_convergence_max, compute the relic and stop the integration
if(N_convergence>=N_convergence_max){
//entropy injection from the point of the last peak until T_stop (the point where interpolation stops)
//the assumption is that at T_stop the universe is radiation dominated with constant entropy.
gamma=plasma.s(axionEOM.T_stop)/plasma.s(T_peak)*std::exp(3*(axionEOM.u_stop-u_peak));
//the relic of the axion
relic=Cosmo<LD>::h_hub*Cosmo<LD>::h_hub/Cosmo<LD>::rho_crit*plasma.s(Cosmo<LD>::T0)/plasma.s(T_peak)/gamma*0.5*
- std::sqrt(axionEOM.axionMass->ma2(Cosmo<LD>::T0,1)*axionEOM.axionMass->ma2(T_peak,1))*
+ std::sqrt(axionMass->ma2(Cosmo<LD>::T0,1)*axionMass->ma2(T_peak,1))*
theta_peak*theta_peak*anharmonicFactor(theta_peak);
break;
}
}
pointSize=points.size();
peakSize=peaks.size();
};
};
#endif
\ No newline at end of file
diff --git a/src/Axion/checks/AxionSolve_check.cpp b/src/Axion/checks/AxionSolve_check.cpp
index ed7d14c..09653b2 100644
--- a/src/Axion/checks/AxionSolve_check.cpp
+++ b/src/Axion/checks/AxionSolve_check.cpp
@@ -1,137 +1,149 @@
#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>
#include"src/Axion/AxionSolve.hpp"
#include"src/AxionMass/AxionMass.hpp"
#include"src/misc_dir/path.hpp"
// define here what you want to print.
#define preintPoints
#define printPeaks
#define printError
// macros for the solver
#ifndef SOLVER
#define SOLVER 1
#define METHOD RODASPR2
#endif
// macros for the numeric type
#ifndef LONG
#define LONG
#endif
#ifndef LD
#define LD LONG double
#endif
int main(int argc, char **argv){
// you can change the axion mass function for T above the interpolation range like this (these are not correct; they only show how it can be done):
// axionMass<LD>.ma2_MAX = [](LD T, LD fa)->LD{ return 0; };
// axionMass<LD>.ma2_MIN = [](LD T, LD fa)->LD{return 0;};//this gives Omegah^2 = 0 (but you are sure that it works)
//model parameters
LD theta_i = 0.94435;
LD fa = 1e12;
// solver parameters
LD umax = 500; //t at which the integration stops
LD TSTOP = 1e-3; // temperature at which integration stops
LD ratio_ini=1e3; // 3H/m_a at which integration begins (should be larger than 500 or so)
// stopping conditions.
// integration stops after the adiabatic invariant hasn't changed
// more than convergence_lim% for N_convergence_max consecutive peaks
unsigned int N_convergence_max=5;
LD convergence_lim=1e-2;
std::string inputFile=std::string(rootDir)+std::string("/UserSpace/InputExamples/RDinput.dat");
/*options for the solver*/
LD initial_step_size=1e-2; //initial step the solver takes.
LD minimum_step_size=1e-8; //This limits the sepsize to an upper limit.
LD maximum_step_size=1e-2; //This limits the sepsize to a lower limit.
LD absolute_tolerance=1e-8; //absolute tolerance of the RK solver
LD relative_tolerance=1e-8; //relative tolerance of the RK solver
LD beta=0.9; //controls how agreesive the adaptation is. Generally, it should be around but less than 1.
/*
the stepsize does not increase more than fac_max, and less than fac_min.
This ensures a better stability. Ideally, fac_max=inf and fac_min=0, but in reality one must
tweak them in order to avoid instabilities.
*/
LD fac_max=1.2;
LD fac_min=0.8;
int maximum_No_steps=int(1e7); //maximum steps the solver can take Quits if this number is reached even if integration is not finished.
+
mimes::AxionMass<LD> axionMass(chi_PATH,0,mimes::Cosmo<LD>::mP);
/*set ma2 for T>TMax*/
LD TMax=axionMass.getTMax();
LD chiMax=axionMass.getChiMax();
axionMass.set_ma2_MAX( [&chiMax,&TMax](LD T, LD fa){return chiMax/fa/fa*std::pow(T/TMax,-8.16);});
/*set ma2 for T<TMin*/
LD TMin=axionMass.getTMin();
LD chiMin=axionMass.getChiMin();
axionMass.set_ma2_MIN( [&chiMin,&TMin](LD T, LD fa){return chiMin/fa/fa;});
// std::function<LD(LD,LD)> ma2 = [](LD T,LD fa){
// LD TQCD=150*1e-3;
// LD ma20=3.1575e-05/fa/fa;
// if(T<=TQCD){return ma20;}
// else{return ma20*std::pow((TQCD/T),8.16);}
// };
// mimes::AxionMass<LD> axionMass(ma2);
mimes::Axion<LD,SOLVER,METHOD<LD>> Ax(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max,
convergence_lim,inputFile, &axionMass,
initial_step_size,minimum_step_size, maximum_step_size, absolute_tolerance, relative_tolerance, beta,
fac_max, fac_min, maximum_No_steps);
Ax.solveAxion();
std::cout<<std::setprecision(5)
<<"theta_i="<<theta_i<<" "<<"f_a="<< fa<<" GeV\n"<<"theta_osc~="<<Ax.theta_osc<<" "
<<"T_osc~="<<Ax.T_osc<<"GeV \n"
<<"Omega h^2="<<Ax.relic<<"\n";
+ /*if you want to change the mass, you need to run restart*/
+ // axionMass.ma2= [](LD T, LD fa){return 1/fa/fa;};
+ // Ax.restart();
+ // Ax.solveAxion();
+
+ // std::cout<<std::setprecision(5)
+ // <<"theta_i="<<theta_i<<" "<<"f_a="<< fa<<" GeV\n"<<"theta_osc~="<<Ax.theta_osc<<" "
+ // <<"T_osc~="<<Ax.T_osc<<"GeV \n"
+ // <<"Omega h^2="<<Ax.relic<<"\n";
+
+
// print all the points
#ifdef preintPoints
std::cout<<"---------------------points:---------------------\n";
std::cout<<"a/a_i\tT [GeV]\ttheta\tzeta\trho_a [GeV^4]"<<std::endl;
for(size_t i=0; i<Ax.pointSize; ++i ){
for(int j=0; j<5; ++j){
std::cout<<Ax.points[i][j];
if(j==4){std::cout<<"\n";}else{std::cout<<"\t";}
}
}
#endif
// print all the peaks
#ifdef printPeaks
std::cout<<"---------------------peaks:---------------------\n";
std::cout<<"a/a_i\tT [GeV]\ttheta\tzeta\trho_a [GeV^4]\tadiabatic_inv [GeV]"<<std::endl;
for(size_t i=0; i<Ax.peakSize; ++i ){
for(int j=0; j<6; ++j){
std::cout<<Ax.peaks[i][j];
if(j==5){std::cout<<"\n";}else{std::cout<<"\t";}
}
}
#endif
// print the local errors
#ifdef printError
std::cout<<"---------------------local errors:---------------------\n";
std::cout<<"dtheta\tdzeta"<<std::endl;
for(size_t i=0; i<Ax.pointSize; ++i ){
std::cout<<Ax.dtheta[i]<<"\t"<<Ax.dzeta[i]<<"\n";
}
#endif
return 0;
}
diff --git a/src/AxionMass/AxionMass.cpp b/src/AxionMass/AxionMass.cpp
index c66d954..43fe63c 100644
--- a/src/AxionMass/AxionMass.cpp
+++ b/src/AxionMass/AxionMass.cpp
@@ -1,49 +1,51 @@
#include<iostream>
#include<functional>
#include"src/AxionMass/AxionMass.hpp"
// macros for the numeric type
#ifndef LONGpy
#define LONGpy
#endif
#ifndef LD
#define LD LONGpy double
#endif
// use this to cast void* to axionMass*
#define Cast(axionMass) static_cast<mimes::AxionMass<LD>*>(axionMass)
typedef LD(* funcType )(LD,LD);
extern "C"{
//constructor
void* INIT_interpolation(char* path, LD minT, LD maxT){
return new mimes::AxionMass<LD>(path,minT,maxT);
}
void* INIT_function(funcType ma2){
return new mimes::AxionMass<LD>( ma2 );
}
// destructor to delete the void*
void DEL(void *axionMass){ delete Cast(axionMass); axionMass=nullptr; }
/*get the values of Tmin, Tmax, chi(Tmin), and chi(Tmax)*/
LD getTMin(void *axionMass){return Cast(axionMass)->getTMin();}
LD getTMax(void *axionMass){return Cast(axionMass)->getTMax();}
LD getChiMin(void *axionMass){return Cast(axionMass)->getChiMin();}
LD getChiMax(void *axionMass){return Cast(axionMass)->getChiMax();}
+ void set_ma2(funcType ma2,void *axionMass){ Cast(axionMass)->set_ma2(ma2); }
+
/*set the functions for ma2 beyond the interpolation limits*/
void set_ma2_MAX(funcType ma2_MAX,void *axionMass){ Cast(axionMass)->set_ma2_MAX(ma2_MAX); }
void set_ma2_MIN(funcType ma2_MIN,void *axionMass){ Cast(axionMass)->set_ma2_MIN(ma2_MIN); }
// axion mass squared
LD ma2(LD T, LD fa, void* axionMass){return Cast(axionMass)->ma2(T,fa);}
}
diff --git a/src/AxionMass/AxionMass.hpp b/src/AxionMass/AxionMass.hpp
index 8d74b00..8eef144 100644
--- a/src/AxionMass/AxionMass.hpp
+++ b/src/AxionMass/AxionMass.hpp
@@ -1,131 +1,133 @@
#ifndef CHI_head
#define CHI_head
#include<iostream>
#include<fstream>
#include<vector>
#include<cmath>
#include<string>
#include<functional>
#include"src/Cosmo/Cosmo.hpp"
#include "src/SimpleSplines/Interpolation.hpp"
namespace mimes{
template<class LD>
class AxionMass {
private:
using VecLD=std::vector<LD>;
using func=std::function<LD(LD,LD)>;
CubicSpline<LD> chi;
VecLD Ttab,chitab;
// these are the functions you can use in order to get the mass and its derivative beyond the
// interpolation limits.
// ma2_MAX and dma2dT_MAX denote the approximation for T>TMAX,
// and ma2_MIN and dma2dT_MIN the approximation for T<TMIN.
LD TMin, TMax, chiMax, chiMin;
func ma2_MAX,ma2_MIN;
public:
func ma2;
AxionMass()=default;
AxionMass(func ma2);
AxionMass(std::string path, LD minT, LD maxT);
AxionMass& operator=(const AxionMass& axionMass)=delete;//delete it explicitely, in order to prevent copies
/*get the values of Tmin, Tmax, chi(Tmin), and chi(Tmax)*/
LD getTMin(){return TMin;}
LD getTMax(){return TMax;}
LD getChiMin(){return chiMin;}
LD getChiMax(){return chiMax;}
+
+ void set_ma2(func ma2){ this->ma2=ma2;}
/*set the functions for ma2 beyond the interpolation limits*/
void set_ma2_MAX(func ma2_MAX){ this->ma2_MAX=ma2_MAX;}
void set_ma2_MIN(func ma2_MIN){ this->ma2_MIN=ma2_MIN;}
LD operator()(LD T, LD fa){return ma2(T,fa);}//overload this just in case you need it.
};
/*
The constructor of the class.
path: a string with a valid path of the data file. The columns of the file must be:
T [GeV] chi [GeV^4]
Here chi is defined through m_a^2(T) = chi(T)/f_a^2, ie chi=m_a^2_{f_a=1 GeV}.
minT: minimum temperature of interpolation. Below this temperature m_a assumed to be constant.
maxT: maximum temperature of interpolation. Above this temperature m_a assumed to follow
m_a^2(T)=chiMax/f_a^2*(T/TMax)^{-8.16},
with chiMax=chi(TMax), and TMax the temberature closest to maxT.
Note: If the file has only one point, the mass is assumed to be constant!
*/
template<class LD>
AxionMass<LD>::AxionMass(std::string path, LD minT, LD maxT){
unsigned int N=0;
LD T,chi;
std::ifstream data_file(path,std::ios::in);
//quit if the file does not exist.
if(not data_file.good()) {
std::cerr << path << " does not exist.";
std::cerr <<" Please make sure to provide a valid path for the AxionMass datafile";
std::cerr <<" in MiMeS/Definitions.mk\n";
exit(1);
}
while (not data_file.eof()){
data_file>>T;
data_file>>chi;
if(T>=minT and T<=maxT){
Ttab.push_back(T); //temperature in GeV
chitab.push_back(chi); //chi in GeV**4
N++;
}
}
data_file.close();
TMin=Ttab[0];
TMax=Ttab[N-1];
chiMax=chitab[N-1];
chiMin=chitab[0];
this->ma2_MAX=[this](LD T, LD fa){return this->chiMax/fa/fa;};
this->ma2_MIN=[this](LD T, LD fa){return this->chiMin/fa/fa;};
this->chi=CubicSpline<LD>(&Ttab,&chitab);
ma2=[this](LD T, LD fa){
// axion mass squared at temperature T and f_\alpha=fa
if(T>=TMax){return ma2_MAX(T,fa);}//use this beyond the upper limit
if(T<=TMin){return ma2_MIN(T,fa);}//use this beyond the lower limit
return this->chi(T)/fa/fa;
};
}
/*
Constructor.
LD ma2(LD T,LD fa): the axion mass squared, which is then accessed by AxionMass::ma2.
The derivative wrt the temperature (AxionMass::dma2dT) is numerically approximated.
Note: LD here is a macro for "double" of "long double"
*/
template<class LD>
AxionMass<LD>::AxionMass(func ma2){
// the mass of the axion squared
this->ma2=ma2;
}
};
#endif
\ No newline at end of file
diff --git a/src/interfacePy/Axion/Axion.py b/src/interfacePy/Axion/Axion.py
index 2a52596..bd0349b 100644
--- a/src/interfacePy/Axion/Axion.py
+++ b/src/interfacePy/Axion/Axion.py
@@ -1,321 +1,330 @@
from ctypes import CDLL, c_longdouble, c_double, c_char_p, c_void_p, c_int, POINTER
from numpy import array as np_array
from time import time
from sys import path as sysPath
from os import path as osPath
sysPath.append(osPath.join(osPath.dirname(__file__), '../../'))
from misc_dir.path import _PATH_
from misc_dir.type import cdouble
from interfacePy.AxionMass import AxionMass
cint=c_int
void_p=c_void_p
char_p=c_char_p
#load the library
axionLib = CDLL(_PATH_+'/lib/Axion_py.so')
# set the argumet and return types of the different functions
axionLib.INIT.argtypes= cdouble, cdouble, cdouble, cdouble, cdouble, cint, cdouble, char_p, void_p, cdouble, cdouble, cdouble, cdouble, cdouble, cdouble, cdouble, cdouble, cint
axionLib.INIT.restype = void_p
axionLib.DEL.argtypes= void_p,
axionLib.DEL.restype = None
axionLib.setTheta_i.argtypes= cdouble, void_p
axionLib.setTheta_i.restype = None
axionLib.MAKE.argtypes= void_p,
axionLib.MAKE.restype = None
+axionLib.restart.argtypes= void_p,
+axionLib.restart.restype = None
+
axionLib.getPointSize.argtypes= void_p,
axionLib.getPointSize.restype = cint
axionLib.getPeakSize.argtypes= void_p,
axionLib.getPeakSize.restype = cint
axionLib.getResults.argtypes=POINTER(cdouble), void_p
axionLib.getResults.restype=None
axionLib.getErrors.argtypes=POINTER(cdouble),POINTER(cdouble), void_p
axionLib.getErrors.restype=None
axionLib.getPoints.argtypes=POINTER(cdouble),POINTER(cdouble),POINTER(cdouble),POINTER(cdouble),POINTER(cdouble),void_p
axionLib.getPoints.restype=None
axionLib.getPeaks.argtypes=POINTER(cdouble),POINTER(cdouble),POINTER(cdouble),POINTER(cdouble),POINTER(cdouble),POINTER(cdouble),void_p
axionLib.getPeaks.restype=None
class Axion:
'''
class that solves the axion eom and stores the resulting evolution of different quantities (eg the angle),
as well as the relic.
Methods:
solveAxion(): solves the EOM and stores the T_osc, \\theta_osc, and \\Omega h^2;
the corresponding variables are self.T_osc, self.theta_osc, and self.relic
getPoints(): This will make numpy arrays that contain the scale factor (self.a),
temperature (self.T), \\theta (self.theta), its derivative with respect to u (self.zeta),
and the energy density of the axion (self.rho_axion).
getPeaks(): This will make numpy arrays that contain the scale factor (self.a),
temperature (self.T), \\theta (self.theta), its derivative with respect to u (self.zeta),
the energy density of the axion (self.rho_axion), and the adiabatic invariant (self.adiabatic_invariant)
on the peaks of the oscillation.
setTheta_i(theta_i): This sets another initial value for theta, without rebuilding
the interpolations (should be faster than declaring another instance).
Under the hood, the constructor gets a new (pointer to an) instance of
mimes::Axion and casts it to void*. Then every member function of this class casts this void* to
mimes::Axion* and calls the corresponding mimes::Axion member from the shared library.
However, the only thing one needs to remember is to delete any instance of this class
once the calculations are done.
'''
def __init__(self, theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max,convergence_lim,inputFile,
axionMass,
initial_step_size=1e-2, minimum_step_size=1e-8, maximum_step_size=1e-2,
absolute_tolerance=1e-8, relative_tolerance=1e-8,
beta=0.9, fac_max=1.2, fac_min=0.8, maximum_No_steps=int(1e7)):
'''
The constructor of the Axion class.
theta_i: initial angle
fa: PQ scale in GeV (the temperature dependent mass is defined as m_a^2(T) = \chi(T)/f^2)
umax: if u>umax the integration stops (rempember that u=log(a/a_i))
TSTOP: if the temperature drops below this, integration stops
ratio_ini: integration starts when 3H/m_a<~ratio_ini (this is passed to AxionEOM,
in order to make the interpolations start at this point)
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
inputFile: file that describes the cosmology. the columns should be: u T[GeV] logH
axionMass: instance of AxionMass
-----------Optional arguments------------------------
initial_stepsize: initial step the solver takes.
maximum_stepsize: This limits the sepsize to an upper limit.
minimum_stepsize: This limits the sepsize to a lower limit.
absolute_tolerance: absolute tolerance of the RK solver
relative_tolerance: relative tolerance of the RK solver
Note:
Generally, both absolute and relative tolerances should be 1e-8.
In some cases, however, one may need more accurate result (eg if f_a is extremely high,
the oscillations happen violently, and the ODE destabilizes). Whatever the case, if the
tolerances are below 1e-8, long doubles *must* be used.
beta: controls how agreesive the adaptation is. Generally, it should be around but less than 1.
fac_max, fac_min: the stepsize does not increase more than fac_max, and less than fac_min.
This ensures a better stability. Ideally, fac_max=inf and fac_min=0, but in reality one must
tweak them in order to avoid instabilities.
maximum_No_steps: maximum steps the solver can take Quits if this number is reached even if integration
is not finished.
'''
self.theta_i, self.fa=theta_i, fa
self.voidAx=void_p()
_file_=char_p(bytes(inputFile, encoding='utf-8'))
self.axionMass=axionMass
self.voidAx=axionLib.INIT(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max, convergence_lim, _file_,self.axionMass.pointer(),
initial_step_size,minimum_step_size, maximum_step_size, absolute_tolerance, relative_tolerance, beta,
fac_max, fac_min, maximum_No_steps)
self.a_peak=[]
self.T_peak=[]
self.theta_peak=[]
self.zeta_peak=[]
self.adiabatic_invariant=[]
self.rho_axion_peak=[]
self.a=[]
self.T=[]
self.theta=[]
self.zeta=[]
self.rho_axion=[]
self.dtheta=[]
self.dzeta=[]
def __del__(self):
'''
Overloaded del.
Also deallocates self.voidAx.
'''
axionLib.DEL(self.voidAx)
del self.voidAx
del self.axionMass
del self.T_osc
del self.a_osc
del self.theta_osc
del self.relic
del self.a_peak
del self.T_peak
del self.theta_peak
del self.zeta_peak
del self.adiabatic_invariant
del self.rho_axion_peak
del self.a
del self.T
del self.theta
del self.zeta
del self.rho_axion
del self.dtheta
del self.dzeta
+
def setTheta_i(self,theta_i):
'''set another initial value of \\theta without rebuilding the interpolations'''
axionLib.setTheta_i(theta_i,self.voidAx)
self.theta_i=theta_i
self.theta_osc=0
self.relic=0
self.a_peak=[]
self.T_peak=[]
self.theta_peak=[]
self.zeta_peak=[]
self.adiabatic_invariant=[]
self.rho_axion_peak=[]
self.a=[]
self.T=[]
self.theta=[]
self.zeta=[]
self.rho_axion=[]
self.dtheta=[]
self.dzeta=[]
+
+ def restart(self):
+ self.setTheta_i(self.theta_i)
+ axionLib.restart(self.voidAx)
+
def solveAxion(self):
'''
Solve the Axion eom.
Running this function we get:
T_osc [GeV] (self.T_osc), a_osc/a_i (self.a_osc), \theta_osc (self.theta_osc),
and \Omega h^2 (in self.relic).
This function also returns the it took to run it.
'''
time0=time()
axionLib.MAKE(self.voidAx)#run the solver
#get the results
points=(cdouble * 4)()
axionLib.getResults(points,self.voidAx)
self.T_osc=points[0]
self.a_osc=points[1]
self.theta_osc=points[2]
self.relic=points[3]
#return the time it took to run
return time()-time0
def getPeaks(self):
'''
This function stores the peaks of the oscillation
in the variables (with selfexplanatory names):
self.a_peak
self.T_peak
self.theta_peak
self.zeta_peak
self.adiabatic_invariant
self.rho_axion_peak
'''
peakSize=axionLib.getPeakSize(self.voidAx)
ArrP = cdouble * peakSize
self.a_peak=ArrP()
self.T_peak=ArrP()
self.theta_peak=ArrP()
self.zeta_peak=ArrP()
self.adiabatic_invariant=ArrP()
self.rho_axion_peak=ArrP()
axionLib.getPeaks(self.a_peak,self.T_peak,self.theta_peak,self.zeta_peak,self.rho_axion_peak,self.adiabatic_invariant,self.voidAx)
self.a_peak=np_array(list(self.a_peak))
self.theta_peak=np_array(list(self.theta_peak))
self.T_peak=np_array(list(self.T_peak))
self.zeta_peak=np_array(list(self.zeta_peak))
self.rho_axion_peak=np_array(list(self.rho_axion_peak))
self.adiabatic_invariant=np_array(list(self.adiabatic_invariant))
def getPoints(self):
'''
This function stores the points of integration
in the variables (with selfexplanatory names):
self.a
self.T
self.theta
self.zeta
self.rho_axion
'''
pointSize=axionLib.getPointSize(self.voidAx)
Arr = cdouble * pointSize
self.a=Arr()
self.T=Arr()
self.theta=Arr()
self.zeta=Arr()
self.rho_axion=Arr()
axionLib.getPoints(self.a,self.T,self.theta,self.zeta,self.rho_axion,self.voidAx)
self.a=np_array(list(self.a))
self.T=np_array(list(self.T))
self.theta=np_array(list(self.theta))
self.zeta=np_array(list(self.zeta))
self.rho_axion=np_array(list(self.rho_axion))
def getErrors(self):
'''
This function stores the local errors of theta and zeta in
self.dtheta and self.dzeta
'''
pointSize=axionLib.getPointSize(self.voidAx)
Arr = cdouble * pointSize
self.dtheta=Arr()
self.dzeta=Arr()
axionLib.getErrors(self.dtheta,self.dzeta,self.voidAx)
self.dtheta=np_array(list(self.dtheta))
self.dzeta=np_array(list(self.dzeta))
diff --git a/src/interfacePy/AxionMass/AxionMass.py b/src/interfacePy/AxionMass/AxionMass.py
index 6376ebf..ba8f7c8 100644
--- a/src/interfacePy/AxionMass/AxionMass.py
+++ b/src/interfacePy/AxionMass/AxionMass.py
@@ -1,110 +1,117 @@
from ctypes import CDLL, c_longdouble, c_double, c_void_p, c_char_p, CFUNCTYPE
from sys import path as sysPath
from os import path as osPath
sysPath.append(osPath.join(osPath.dirname(__file__), '../../'))
from misc_dir.path import _PATH_
from misc_dir.type import cdouble
void_p=c_void_p
char_p=c_char_p
#load the library
AxionMassLib = CDLL(_PATH_+'/lib/libma.so')
AxionMassLib.INIT_interpolation.argtypes=char_p, cdouble, cdouble
AxionMassLib.INIT_interpolation.restype=void_p
AxionMassLib.INIT_function.argtypes=CFUNCTYPE(cdouble, cdouble, cdouble),
AxionMassLib.INIT_function.restype=void_p
AxionMassLib.DEL.argtypes= void_p,
AxionMassLib.DEL.restype = None
AxionMassLib.ma2.argtypes=cdouble, cdouble, void_p
AxionMassLib.ma2.restype=cdouble
AxionMassLib.getTMin.argtypes=void_p,
AxionMassLib.getTMin.restype=cdouble
AxionMassLib.getTMax.argtypes=void_p,
AxionMassLib.getTMax.restype=cdouble
AxionMassLib.getChiMin.argtypes=void_p,
AxionMassLib.getChiMin.restype=cdouble
AxionMassLib.getChiMax.argtypes=void_p,
AxionMassLib.getChiMax.restype=cdouble
AxionMassLib.set_ma2_MAX.argtypes=CFUNCTYPE(cdouble, cdouble, cdouble),void_p
AxionMassLib.set_ma2_MAX.restype=None
AxionMassLib.set_ma2_MIN.argtypes=CFUNCTYPE(cdouble, cdouble, cdouble),void_p
AxionMassLib.set_ma2_MIN.restype=None
+AxionMassLib.set_ma2.argtypes=CFUNCTYPE(cdouble, cdouble, cdouble),void_p
+AxionMassLib.set_ma2.restype=None
+
class AxionMass:
'''AxionMass class: This is wrapper for the mimes::AxionMass class.
There member function available are:
pointer(): returns the void pointer of the underlying AxionMass instance. This is useful because
the AxionMass class needs a void pointer to a mimes::AxionMass instance.
ma2(T,fa): the axion mass squared (in GeV^2) at temperature T (in GeV) and fa (in GeV)
'''
def __init__(self,*args):
'''
The constructor of the class.
If only one agrument is given, it is assumed that it is the function of the axion mass squared.
If three arguments are given, the first should be a path of a datafile that contains T (in GeV)
and chi (in GeV^4), with increasing T. The second and third arguments are the minimum and maximum
interpolation temperatures.
'''
self.c_ma2=CFUNCTYPE(cdouble, cdouble, cdouble)(lambda T,fa:0)
self.c_ma2_MAX=CFUNCTYPE(cdouble, cdouble, cdouble)(lambda T,fa:0)
self.c_ma2_MIN=CFUNCTYPE(cdouble, cdouble, cdouble)(lambda T,fa:0)
if len(args) == 1:
self.c_ma2=CFUNCTYPE(cdouble, cdouble, cdouble)(args[0])
self.voidAxM=AxionMassLib.INIT_function( self.c_ma2 )
if len(args) == 3:
path, minT,maxT=args[0],args[1],args[2]
_file_=char_p(bytes(path, encoding='utf-8'))
self.voidAxM=AxionMassLib.INIT_interpolation(_file_,minT,maxT)
def pointer(self):
return self.voidAxM
def __del__(self):
AxionMassLib.DEL(self.voidAxM)
del self.voidAxM
del self.c_ma2
del self.c_ma2_MAX
del self.c_ma2_MIN
def getTMin(self):
return AxionMassLib.getTMin(self.voidAxM)
def getTMax(self):
return AxionMassLib.getTMax(self.voidAxM)
def getChiMin(self):
return AxionMassLib.getChiMin(self.voidAxM)
def getChiMax(self):
return AxionMassLib.getChiMax(self.voidAxM)
def set_ma2_MAX(self,ma2_MAX):
self.c_ma2_MAX= CFUNCTYPE(cdouble, cdouble, cdouble)(ma2_MAX)
AxionMassLib.set_ma2_MAX(self.c_ma2_MAX,self.voidAxM)
def set_ma2_MIN(self,ma2_MIN):
self.c_ma2_MIN= CFUNCTYPE(cdouble, cdouble, cdouble)(ma2_MIN)
AxionMassLib.set_ma2_MIN(self.c_ma2_MIN,self.voidAxM)
+ def set_ma2(self,ma2):
+ self.c_ma2= CFUNCTYPE(cdouble, cdouble, cdouble)(ma2)
+ AxionMassLib.set_ma2(self.c_ma2,self.voidAxM)
+
def ma2(self,T,fa):
return AxionMassLib.ma2(T,fa,self.pointer())
if __name__=="__main__":
pass
\ No newline at end of file