#  NSC++: Non-Standard Cosmologies in C++

Dimitrios Karamitros

###### Abstract

We introduce NSC++, a header-only C++ library that simulates the evolution of
the plasma and a decaying fluid in the early Universe. NSC++ can be used in
C++ programs or called directly from python scripts without significant
overhead. There is no special installation process or external dependencies.
Furthermore, there are example programs that can be modified to handle several
cases.

Keywords: Cosmology, simulation, high energy physics.

Program summary:

Program title: NSC++.

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

Programming language: C++ (C++17$+$) and python (3.7$+$).

Licensing provisions: MIT license.

Nature of problem: Solves equations that describe the evolution of the plasma
along with a fluid that increases the entropy of the plasma during the early
Universe.

Solution method: Embedded Runge-Kutta for the numerical integration of the
system of differential equations. The user can choose between explicit and
Rosenbrock methods. There are several of Butcher tableaux already implemented,
but the user can implement their own. The interpolations of the relativistic
degrees of freedom are accomplished by cubic spline interpolation.

Restrictions: The value of the pressure over the energy density and the energy
loss parameter of the fluid are assumed to be constant. It is also assumed
that there is no energy leaking from the plasma to the fluid.

###### Contents

  1. 1 Introduction
  2. 2 Evolution equations
    1. 2.1 Notation
  3. 3 NSC++ usage
    1. 3.1 Assumptions
    2. 3.2 Using NSC++
    3. 3.3 Using NSC++ in python
  4. 4 Example
    1. 4.1 C++
    2. 4.2 python
    3. 4.3 Result
  5. 5 Summary
  6. 6 Acknowledgements
  7. A Evolution class in C++
  8. B Evolution class in python
  9. C Quick guide to the user input

##  1 Introduction

The configuration of the components of the Universe today might have been
influenced by its nature at times that cannot be probed. Therefore, there is
the possibility that the Universe experienced a period of non-standard
evolution. A scenario that deviates from the standard cosmological case is
called non-standard cosmology (NSC). NSCs in general affect the evolution of
other components of the Universe, and their detailed examination can help us
probe the nature of the Universe at temperatures inaccessible by direct
observations (e.g. of CMB). For example, the study of dark matter within NSCs
can allow possibilities that are otherwise excluded. This means that if dark
matter is observed and turns out to be excluded under the standard
cosmological model, then its nature can reveal aspects of the Universe at
times relevant to the dark matter production (see, e.g. [1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14, 15, 16]). Several non-standard cosmological
scenarios have been studied in the literature (e.g. refs. [17, 18, 19, 20, 21,
22] with some reviews found in refs. [23, 24]), however there is no available
tool that will automate and standardise the simulation the evolution of the
Universe. NSC++ is an attempt to accomplish this by being an easy-to-use and
modifiable C++ library that can handle the most common NSCs. Since in most
cases the expansion mode of the Universe is just a background, NSC++ can be
used to simulate it and use the results to study the evolution of other
components of the Universe (e.g. dark matter or lepton number) under this
background. In particular, NSC++ can be used along with MiMeS [25], since it
can output the evolution of the temperature and Hubble parameter which can be
used as an input by MiMeS.

We note that there is another package, built specifically for the evolution of
cosmic relics in matter dominated Universe, EvoEMD [26]. NSC++ is different
because the evolution of the temperature of the plasma is computed by
manipulating the equation that describes the evolution of the entropy density,
while EvoEMD solves

$\dfrac{d\rho_{\rm R}}{dt}=-4H\rho_{\rm R}+C_{\rm EMD}\;,$ (1.1)

with $C_{\rm EMD}$ the additional matter contribution. This equation relies on
the approximation $\rho_{\rm R}\approx 3/4\,sT$ that does not hold in general.
Also, NSC++ allows the user to consider different equations of state for the
fluid, $\Phi$. As we explain in this article, NSC++ adheres to design
principles similar to MiMeS [25]; i.e. the general usage of both is almost
identical, NSC++ only needs a C++ compiler that supports the c++17 standard
with no other external dependences, and there will always be dedicated
versions of the underline numerical libraries [27, 28] guaranteed to work with
the current version of NSC++. Moreover, there are multiple Runge-Kutta (RK)
methods [27] included with the possibility to implement new ones by the user.
This can help the users inspect the numerical validity of the results.

This article is organised as follows. Section 2 is brief reminder of the
system that NSC++ simulates and a description of the internal notation it
uses. The next section describes the usage of NSC++ in general, from
downloading to compiling, along with a brief description of the various
classes and functions available to the user. In Section 4 there are complete
examples written in both C++ and python, and a brief mention of all the
results one can output. Section 5 summarises the article and provides future
directions. Appendices (A) and (B) describe in detail the implementation of
the classes available to the user, while in Appendix (C) there are tables that
summarise all the possible inputs and options.

##  2 Evolution equations

In order to model an NSC, we assume that the Universe at early times was
dominated by two components; the plasma and a fluid ($\Phi$) with an equation
of state

$p_{\Phi}=(c/3-1)\rho_{\Phi}\;,$ (2.1)

with $c$ a constant.

The plasma is assumed to be always in thermal equilibrium, with energy and
entropy densities defined as

$\displaystyle s_{\rm R}(T)$ $\displaystyle=\dfrac{2\pi^{2}}{45}\,h_{\rm
eff}{}(T)\,T^{3}$ (2.2) $\displaystyle\rho_{\rm R}(T)$
$\displaystyle=\dfrac{\pi^{2}}{30}\,g_{\rm eff}{}(T)\,T^{4}\;,$ (2.3)

where $h_{\rm eff}{}(T)$ and $g_{\rm eff}{}(T)$ count the effective
relativistic degrees of freedom (RDOF) of the plasma.

If $\Phi$ dominates for some temperature range, it must lose energy to the
plasma, since the contribution to the energy budget of the Universe of such
fluid is severely constraint at temperatures lower than
$T\sim\mathcal{O}(10)~{}{\rm MeV}$ [24]. We model this by introducing a
constant energy loss rate, $\Gamma_{\Phi}$. The equations that describe the
evolution of both components are

$\displaystyle\dfrac{ds_{\rm R}(T)}{dt}$ $\displaystyle=-3\,H(T)\,s_{\rm
R}(T)+\dfrac{\Gamma_{\Phi}}{T}\,\rho_{\Phi}(T)$ (2.4)
$\displaystyle\dfrac{d\rho_{\Phi}(T)}{dt}$
$\displaystyle=-c\,H(T)\,\rho_{\Phi}(T)-\Gamma_{\Phi}\,\rho_{\Phi}(T)\;,$
(2.5)

with $H$ the Hubble parameter given by

$H(T)=\sqrt{\dfrac{8}{3\,m_{P}^{2}}\ \left[\rho_{\rm
R}(T)+\rho_{\Phi}(T)\right]}\;,$ (2.6)

where $m_{P}=1.22\times 10^{19}~{}{\rm GeV}$.

Following refs. [15, 6], we can parametrise $\Gamma_{\Phi}$ in terms of the
temperature that $\Phi$ would decay ($T_{\rm end}{}$) if the Universe was
radiation dominated, as

$\Gamma_{\Phi}\equiv H_{\rm R}(T=T_{\rm end}{})\;,$ (2.7)

where $H_{\rm R}$ the Hubble parameter for a radiation dominated Universe
obtained from eq. (2.6) with $\rho_{\Phi}=0$.

###  2.1 Notation

In order to solve this system of equations, we change the integration variable
to

$u=\log\frac{a}{a_{\rm i}{}}\;,$ (2.8)

with $a$ the scale factor of the Universe, and $a_{\rm i}{}$ its value at some
initial time. The initial condition can now be chosen at $u_{\rm i}{}=0$,
corresponding to $T=T_{\rm i}{}$ and an initial value for
$\rho_{\Phi}=\rho_{\Phi\,{\rm i}}$. The form of eq. (2.8) implies that the
evolution of the Universe depends only on its relative scale; i.e. if $a_{\rm
i}{}\neq 0$, its value does not matter.

We can reparametrise the temperature of the plasma and the energy density of
$\Phi$ as

$\displaystyle T(u)$ $\displaystyle=T_{\rm i}{}\,f_{\rm R}(u)\,e^{-u}$ (2.9)
$\displaystyle\rho_{\Phi}(u)$ $\displaystyle=\rho_{\Phi\,{\rm
i}}\,f_{\Phi}(u)\,e^{-cu}\;,$ (2.10)

with $f_{\rm R}(u)$ and $f_{\Phi}(u)$ functions of $u$ that obey the initial
conditions $f_{\rm R}(u=0)=f_{\Phi}(u=0)=1$ and the equations

$\displaystyle\dfrac{d\log f_{\rm R}}{du}$
$\displaystyle=1-1/\delta_{h}{}-\dfrac{1}{3}\dfrac{\Gamma_{\Phi}}{H\,\delta_{h}{}}\,\dfrac{\rho_{\Phi}}{T\,s}$
(2.11) $\displaystyle\dfrac{d\log f_{\Phi}}{du}$
$\displaystyle=-\dfrac{\Gamma_{\Phi}}{H}\;,$ (2.12)

with

$\delta_{h}{}=1+\dfrac{1}{3}\dfrac{d\log h_{\rm eff}{}}{d\log T}\;.$ (2.13)

Notice that for $\Gamma_{\Phi}=0$ and $\delta_{h}{}=1$, the functions
introduced are $f_{\rm R}=f_{\Phi}=1$. These are the equations that NSC++
actually solves, using as user input for $T_{\rm end}{}$, $c$, $T_{\rm i}{}$,
and $r=\rho_{\rm R\,i}/\rho_{\Phi\,{\rm i}}$. Then, it automatically
transforms them to the temperature of the plasma and the fluid energy density.
From these, the user can obtain all other relevant cosmological quantities
such as the Hubble parameter or the entropy density of the plasma.

NSC++ closely follows the notation of ref. [6], identifying (at most) four
points where the behaviour of the system changes; ${\rm E1}{}$, ${\rm D1}{}$,
${\rm E2}{}$, and ${\rm D2}{}$. The point ${\rm E1}{}$ (${\rm E2}{}$) is
defined when $\rho_{\Phi}=\rho_{\rm R}$ for the first (second) time. The point
${\rm D1}{}$ (${\rm D2}{}$) is defined as the point at which $\left[(H_{\rm
R,\,end}/H)\,(\rho_{\Phi}/\rho_{\rm R})\right]$ exceeds (drops below) $40\%$;
i.e. before ${\rm D1}{}$ and after ${\rm D2}{}$ the energy loss rate from the
fluid to the plasma is less than $10\%$ of the free dilution rate of
$\rho_{\rm R}$. For the sake of performance, NSC++ identifies these points
approximately during the integration of eqs. (2.11) and (2.12).

In general, there are two cases that exhibit qualitatively different
evolution, which are shown in Figure (1). In Figure (1(a)), we show the
evolution of the comoving energy densities of radiation and $\Phi$ as
functions of $u$ for $T_{\rm end}{}=10^{-2}~{}{\rm GeV}$, $c=3$, $T_{\rm
i}{}=10^{7}~{}{\rm GeV}$, and $r=10^{-2}$. In Figure (1(b)), we show the
evolution of these quantities for $T_{\rm end}{}=10^{4}~{}{\rm GeV}$, $c=6$,
$T_{\rm i}{}=10^{7}~{}{\rm GeV}$, and $r=10^{8}$. These two figures correspond
to the two qualitative typical behaviours we expect for $c<4$ and $c>4$. In
the former, $\Phi$ increases its energy contribution, until it decays away. In
the latter, $\rho_{\Phi}$ monotonically decreases. In this case, the Universe
will become radiation dominated regardless $\Gamma_{\Phi}$, but if
$\rho_{\Phi\,{\rm i}}$ and $\Gamma_{\Phi}$ are large enough the energy of the
plasma will increase. Notice that for $c>4$, There can only be one point of
equality between $\rho_{\rm R}$ and $\rho_{\Phi}$, which we call ${\rm E1}{}$
as this is the naming conversion in NSC++.

(a)

(b)

Figure 1:  (a) The evolution of the comoving energy densities of radiation
(black) and $\Phi$ (dashed red) in the case of early matter domination ($c=3$)
with $T_{\rm end}{}=10^{-2}~{}{\rm GeV}$, $T_{\rm i}{}=10^{7}~{}{\rm GeV}$,
and $r=10^{-2}$. (b) The same for kination ($c=6$) with $T_{\rm
end}{}=10^{4}~{}{\rm GeV}$, $T_{\rm i}{}=10^{7}~{}{\rm GeV}$, and $r=10^{8}$.
The grey vertical lines show the approximate location of the points ${\rm
E_{1,2}}$ and ${\rm D_{1,2}}$.

##  3 NSC++ usage

NSC++ can be obtained directly form github. The safest choice is to download
the “stable” branch, which will always be the latest tested version. The
preferred method for downloading NSC++ is by running: 111The git utility can
be installed by following https://github.com/git-guides/install-git.

⬇

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

NSC++ relies on NaBBODES [27] and SimpleSplines [28], which are developed
independently. In order to get NSC++ with the latest version of these
libraries, one needs to run the following commands

⬇

1 git clone https://github.com/dkaramit/NSCpp.git

2 cd NSCpp

3 git submodule init

4 git submodule update \--remote

These commands download the master branch of NSC++, with the latest versions
of NaBBODES and SimpleSplines. All versions that are not labelled as stable
must be treated with caution, as they may not work properly.

By executing “bash configure.sh" in the root directory of NSC++, the library
will be available for usage within a C++ program. This bash script, writes
some paths in some files for convenience, formats the data for the RDOF in an
acceptable format (if a relative path is given in the NSCpp/Paths.mk file),
and creates some directories that may be useful later. In order to call NSC++
inside python, “make lib" must also be run, which will compile shared
libraries that needed for the python interface.

NSC++ already comes with several examples in C++ and python, which can be
found in the directory NSCpp/UserSpace. In order to compile the C++ examples,
one just has to run “make examples" in the root directory of NSC++ or make
inside the directory of each C++ program.

There is also an option to compile some test programs using make check, which
will create executables inside exec. These can then be executed, in order to
ensure that the code actually works; e.g. no segmentation fault occurs.

###  3.1 Assumptions

NSC++ is designed to be as model agnostic as possible. It can only be used for
constant $c$ and $\Gamma_{\Phi}$, and in cases where there are no plasma
energy leaks into the fluid. In cases that these assumptions do not hold,
knowledge of the exact nature of $\Phi$ is generally needed. It is also
assumed that the plasma does not fall out of thermal equilibrium between the
temperatures we integrate (defined later as Ti and TSTOP), and its energy and
entropy densities are always given by eqs. (2.2) and (2.3).

###  3.2 Using NSC++

Out-of-the-box, and after “bash configure.sh" is run successfully, NSC++ can
be used in a C++ program by including in the source file the following line:

⬇

1 #include "NSCpp.hpp"

Notice that if this file is not in the root directory of NSC++, we need to
compile it using the flag -Ipath-to-root, "path-to-root" the relative (or
absolute) path to the root directory of NSC++; e.g. if the .cpp is in the
NSCpp/UserSpace/Cpp/NSC directory, this flag should be -I../../../.

The class that deals with the system (2.11) and (2.12), is the
nsc::Evolution<LD,SOLVER,METHOD>. In order to begin, we declare an instance of
this class as

⬇

1 nsc::Evolution<LD,Solver,Method> BE;

Here, LD should be the numeric type to be used. Only floating point types are
allowed, and it is recommended to use double or long double. Moreover Solver
and Method depend on the type of RK method the user chooses. The available
choices are shown in Table 3.

The evolution of the Universe depends on the evolution of the plasma, which is
encoded in the class nsc::Cosmo<LD>. An instance of this class is declared as

⬇

1 nsc::Cosmo<LD> plasma(’path_to_data’, minT, maxT);

Here, ’path_to_data’ is the (relative of absolute) path to a file with
tabulated data of the relativistic degrees of freedom of the plasma . It is
assumed that the file has three columns with: $T$ in ${\rm GeV}$, $h_{\rm
eff}{}$, and $g_{\rm eff}{}$. The temperature given in this file has to
increase monotonically, otherwise the program exits with error core $1$.
222NSC++ already includes the standard model RDOF found in [29] in the
directory NSCpp/src/data. If a path is provided in the NSCpp/Paths.mk file
before “bash configure.sh" is run, there will be a global (constant) variable
called cosmo_PATH which can be used after NSCpp.hpp is included. The arguments
minT and maxT dictate the temperatures between which $g_{\rm eff}{}$ and
$h_{\rm eff}{}$ will be interpolated. For temperatures beyond these, both
$g_{\rm eff}{}$ and $h_{\rm eff}{}$ are assumed to be constant. A detailed
description of this class can be found in the Appendix of ref. [25].

The system can be solved by including:

⬇

1 BE.solveNSC(TEND, c, Ti, ratio, TSTOP, umax, &plasma,

2 {

3 .initial_step_size=initial_step_size, .minimum_step_size=minimum_step_size,

4 .maximum_step_size=maximum_step_size,
.absolute_tolerance=absolute_tolerance,

5 .relative_tolerance=relative_tolerance, .beta=beta,

6 .fac_max=fac_max, .fac_min=fac_min, .maximum_No_steps=maximum_No_steps

7 }

8 );

The various parameters are:

  1. 1.

TEND: the values of the $T_{\rm end}{}$ (in ${\rm GeV}$) as introduced in eq.
(2.4). 333Note that $\Gamma_{\Phi}=0$ corresponds to $T_{\rm end}{}=0$.

  2. 2.

c: the parameter that characterises the equation of state of the fluid,
defined as in eq. (2.1).

  3. 3.

Ti: the initial temperature (in ${\rm GeV}$) of integration. Has to be above
TEND.

  4. 4.

ratio: the value of $\rho_{\Phi}/\rho_{\rm R}$ at Ti.

  5. 5.

TSTOP: If the temperature drops below this, integration stops.

  6. 6.

umax : If $u>$umax the integration stops. It can be used to stop internation
if relevant bounds (e.g. [30]) are violated. It is also a failsafe parameter
that ensures that integration stops even in cases where the Universe expands
extremely fast.

  7. 7.

plasma : A pointer to an instance of the Cosmo class.

  8. 8.

initial_stepsize (optional): Size of the initial step the solver takes.

  9. 9.

minimum_stepsize (optional): The step-size will be limited to be above number.

  10. 10.

maximum_stepsize (optional): The step-size will be limited to be below number.

  11. 11.

absolute_tolerance (optional): Absolute tolerance of the RK solver.

  12. 12.

relative_tolerance (optional): Relative tolerance of the RK solver.

  13. 13.

beta (optional): Aggressiveness of the adaptation process of the solver.
Generally, it should be around but less than 1.

  14. 14.

fac_max, fac_min (optional): The solver will keep step-size to from increasing
more than fac_max and less than fac_min at each iteration. This makes
adaptation of the step-size more stable.

  15. 15.

maximum_No_steps (optional): Maximum steps the solver can take. Quits if this
number is reached even if the temperature is still larger than TEND.

A detailed description of the effect of the various optional parameters on the
operation of the ODE solver can be found in the Appendices of ref. [25].
Notice that the optional parameters are passed using aggregation. 444This
means that all arguments inside {...} define a temporary an instance of a
struct that simply holds these arguments.

The solveNSC function returns true if the temperature reached TSTOP for
$u<=$umax and false otherwise. This can be used to test whether the
integration completed successfully. The evolution of $T$ and $\rho_{\Phi}$ is
stored in BE.T and BE.rhoPhi, with each element corresponding to integration
points $u$ (accessed as BE.u). The local errors of $T$ and $\rho_{\Phi}$
correspond to the variables BE.dT and BE.drhoPhi. One can also obtain the
points BE.TE1 ($T_{\rm E1}{}$), BE.TE2 ($T_{\rm E2}{}$), BE.TD1 ($T_{\rm
D1}{}$), BE.TD2 ($T_{\rm D2}{}$). The corresponding values of $u$ are BE.uE1,
BE.uE2, BE.uD1, BE.uD2.

###  3.3 Using NSC++ in python

The python interface modules are in the directory NSCpp/src/interfacePy. The
usage of the classes is made as similar as possible to the C++ case. However,
one should keep in mind that the various template arguments discussed in the
C++ case have to be chosen at compile-time. That is, for the python interface,
one needs to choose the numeric type and RK method to be used when the shared
libraries are compiled. This is accomplished by assigning the relevant
variables in NSCpp/Definitions.mk before running “bash configure.sh" and “make
lib". The various options are the same as in Section 4 of ref. [25], and
summarised in Table 4.

The two relevant classes are defined in the module interfacePy, and can be
imported in a python script as

⬇

1 from sys import path as sysPath

2 sysPath.append(’path_to_src’)

3 from interfacePy.Evolution import Evolution

4 from interfacePy.Cosmo import Cosmo

The string ’path_to_src’ must be the relative path to the NSCpp/src directory.
For example, if the script is located in NSCpp/UserSpace/Python, ’path_to_src’
should be ’../../src’.

We can define an Evolution instance as follows

⬇

1 BE=Evolution()

An instance of the Cosmo class is declared as

⬇

1 plasma=Cosmo(’path_to_data’, minT, maxT)

The arguments are the same as in C++.

The system, then, is solved by running:

⬇

1 BE.solveNSC(TEND, c, Ti, ratio, TSTOP, umax, plasma,

2 initial_step_size, minimum_step_size, maximum_step_size, absolute_tolerance,

3 relative_tolerance, beta, fac_max, fac_min, maximum_No_steps)

4 )

The arguments are identical to the C++ case, outlined in Table 1, but the
instance of Cosmo (plasma) is passed by value. A brief description of this
function can also be found by running BE.solveNSC after loading the module.
This function returns True (False) if the temperature reached TSTOP for
$u<=$umax ($u>$umax), and the time it took to execute it in seconds. In
contrast to the C++ implementation, this only gives us access to the points
where the behaviour changes; the corresponding variables are BE.TE1, BE.TE2,
BE.TD1, BE.TD2, BE.uE1, BE.uE2, BE.uD1, and BE.uD2. In order to get the
evolution of $T$ and $\rho_{\Phi}$, we need to run

⬇

1 BE.getPoints()

This will fill the numpy [31] arrays with the integration variable and the
solutions; BE.u, BE.T, and BE.rhoPhi. The following line:

⬇

1 BE.getErrors()

gives access to the local errors BE.dT and BE.drhoPhi.

It is essential to manually delete both BE and plasma if they are no longer
needed or before reassigning them. The reason is that there are underlying
C-pointers that manage the python interface for each instance of the classes.
If the instances are reassigned without deletion, these pointers can no longer
be accessed. This is a memory-leak, which can only be fixed by exiting the
script. That is, one must simply run

⬇

1 del BE, plasma

after BE and plasma served their purpose or before reassigning these
variables.

##  4 Example

In this section we will show an example code that solves the system of Figure
(1(a)) in both C++ and python.

###  4.1 C++

In C++, we need to include the header file NSCpp.hpp from the root directory
of NSC++. Then, we declare an instance of the Evolution class. The system is
solved by calling the solveNSC function with the inputs described in Table 1.
The following code can be used to solve the system (2.4) and (2.5). This code
will print the values of $T_{\rm E_{1,2}}$ and $T_{\rm D_{1,2}}$, and all the
integration points with their local errors.

⬇

1#include<iostream>

2#include<iomanip>

3

4//Include everything you need from NSC++

5#include"NSCpp.hpp"

6

7int main(){

8 //Use cosmo_PATH to interpolate $h_{\rm eff}{}$ and $g_{\rm eff}{}$ from T=0
to T=mP.

9 nsc::Cosmo<long double> plasma(cosmo_PATH, 0, nsc::Cosmo<long double>::mP);

10

11 //Declare Evolution instance using the Rosenbrock method RODASPR2

12 nsc::Evolution<long double,1,RODASPR2<long double>> BE;

13

14 //Declare parameters

15 long double TEND=1e-2, c=3, Ti=1e7, ratio=1e-2, TSTOP=1e-4, umax=200;

16

17 bool check=BE.solveNSC(TEND, c, Ti, ratio, TSTOP, umax, &plasma,

18 {

19 .initial_step_size=1e-2, .minimum_step_size=1e-8, .maximum_step_size=1e-2,

20 .absolute_tolerance=1e-11, .relative_tolerance=1e-11, .beta=0.9,

21 .fac_max=1.2, .fac_min=0.8, .maximum_No_steps=10000000

22 });

23

24 if(check){

25 // If the solver returns true, the solver probably worked.

26 std::cout<<std::setprecision(5);

27 std::cout<<BE.TE1<<"\t"<<BE.TE2<<"\t"<<BE.TD1<<"\t"<<BE.TD2<<"\n";

28

29 //print the results and the errors

30 for(size_t i=0; i<BE.pointSize; ++i ){

31 std::cout<<std::left<<std::setw(15)<<BE.u[i]<<std::setw(15);

32 std::cout<<std::left<<BE.T[i]<<std::setw(15)<<BE.dT[i]<<std::setw(15);

33 std::cout<<std::left<<BE.rhoPhi[i]<<std::setw(15)<<BE.drhoPhi[i]<<"\n";

34 }

35 }else{

36 // If the solver returns false, you may need a larger umax.

37 std::cerr<<"Something␣went␣wrong.␣Try␣using␣larger␣value␣for␣umax\n";

38 exit(1);

39 }

40

41 return 0;

42}

Notice that if solveNSC returns false, there is an error message printed, and
the code exits with error code $1$.

In this program we have used long double as the numeric type in all
declarations, with the RODASPR2 Rosenbrock method [32] for the RK solver.
Assuming that this code is written in a file with name example.cpp, we can
compile it with

⬇

1 g++ -std=c++17 -O3 -lm -I’path_to_root’ example.cpp -o example.run

or

⬇

1 clang -lstdc++ -std=c++17 -O3 -lm -I’path_to_root’ example.cpp -o
example.run

with ’path_to_root’ the root directory of NSC++. This will create the
executable example.run. Note that the variable cosmo_PATH is defined in
NSCpp/src/misc_dir/path.hpp, created by the script NSCpp/configure.sh if a
relative path is given in NSCpp/PATHS.mk. As This variable is the absolute
path to the RDOF file, the executable can be copied to any other place in the
same system.

###  4.2 python

The python interface of NSC++ can handle same system as in the previous
subsection. This is shown in the following code:

⬇

1# append the path to the src directory

2from sys import path as sysPath

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

4

5#load the NSC module

6from interfacePy.Evolution import Evolution

7

8# load Cosmo and the Planck mass

9from interfacePy.Cosmo import Cosmo,mP

10

11#This gives you access to the path of the rdof file.

12from misc_dir.path import cosmo_PATH

13

14#Instance of the Cosmo class. Interpolate from T=0 to T=mP

15plasma=Cosmo(cosmo_PATH,0,mP)

16

17# Evolution instance

18BE=Evolution()

19

20# solve the system

21check,time=BE.solveNSC(TEND=1e-2, c=3, Ti=1e7, ratio=1e-2, umax=500,
TSTOP=1e-4, plasma=plasma,

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

23 absolute_tolerance=1e-11, relative_tolerance=1e-11, beta=0.9, fac_max=1.2,

24 fac_min=0.8,maximum_No_steps=10000000)

25

26if check:

27 # get points

28 BE.getPoints()

29 # get errors

30 BE.getErrors()

31

32print(BE.TE1,BE.TE2,BE.TD1,BE.TD2)

33

34for i,u in enumerate(BE.u):

35 print(u,BE.T[i],BE.dT[i],BE.rhoPhi[i],BE.drhoPhi[i])

36else:

37 print("Something␣went␣wrong.␣Try␣using␣larger␣value␣for␣umax")

38 exit(1)

39

40#run the destructors

41del BE

42del plasma

This script is assumed to be located in a subdirectory of the root directory
of NSC++, e.g. inside NSCpp/UserSpace. If it is in another directory, then
only the third line should be changed to sysPath.append(’path_to_root’), with
’path_to_root’ the path to the root directory of NSC++.

###  4.3 Result

In both C++ and python, the results we obtain are identical, and evolution of
the comoving energy densities is given by Figure (1(a)).

(a)

(b)

Figure 2:  (a) The local integration error that corresponds to the temperature
of the plasma (black) and $\rho_{\Phi}$ (dashed red) in the case of early
matter domination. (b) A histogram of the number of states taken by the solver
at various $u$. The parameters used are as in Figure 1(a).

The relative errors associated with the integration of eqs. (2.11) and (2.12)
are shown in Figure (2(a)). As we can see, the local errors appear to be
negligible compared to the values of the corresponding quantities, and they
appear to need controlling only as long as the decays of $\Phi$ are active.
Once the fluid has decayed away, its relative error seems to increase, but
this is only because $\rho_{\Phi}$ almost vanishes. In Figure (2(b)), we show
a histogram of the number of steps taken by the RK solver. In accordance to
the increase in the errors, the number of steps increases after ${\rm E1}{}$,
which means that the step size changes in order to regulate the integration
error.

We also show the relative local errors and histogram corresponding to Figure
(1(b)), which are obtained by changing the parameters of the script to those
shown in that figure.

(a)

(b)

Figure 3: Same as Figure 2 for the parameters of Figure 1(b).

In Figures (3(a)) and (3(b)) we can see that the errors are still small, and
after ${\rm D1}{}$ the step size is automatically adjusted in order to keep it
from increasing.

Figures of errors and histogram such as Figures (2) and (3) can be used to
identify cases that the solver behaves badly. In general, these plots do not
provide the whole picture. It should be considered good practice to use
different methods and choices of the various RK parameters (e.g. lower of
higher values of relative_tolerance and absolute_tolerance) in order to make
sure that the results are consistent. Ideally, these results should be also
compared against approximate solutions to eqs. (2.11) and (2.12) such as the
ones explained in detail in ref. [6].

##  5 Summary

We have introduced NSC++, a simple library that simulates the evolution of the
Universe assuming the appearance of a decaying fluid. NSC++ provides several
different methods to solve the radiation-fluid system and also gives access to
the local integration errors, which allows the user to try different
approaches and decide whether the results can be trusted.

In Section (2) we have shown the equations that NSC++ solve and briefly
described the general form of the solutions we expect. We also showed
explicitly the internal notation NSC++ uses, which should enable others modify
the code in order to meet their needs.

NSC++ can be downloaded and configured following the steps in Section (3).
This section also broadly describes the various classes and functions that can
be used by the user in both C++ and python, with more details given in the
Appendix. Complete code examples in C++ as well as python are shown in Section
(4).

NSC++ has limitations as it is assumed that there is only one additional
component ($\Phi$) with an equation of state of the given by eq. (2.1) and the
evolution of its energy density described by eq. (2.5). NSC++ will keep
evolving and new versions will be able to handle non-minimal cases such as the
ones studied in refs. [9, 33, 34, 35]). Another useful extension that will be
included in a future version is a module that compares against observational
constraints given an inflation model.

##  6 Acknowledgements

This work is supported by the Lancaster–Manchester–Sheffield Consortium for
Fundamental Physics, under STFC research grant ST/T001038/1.

## Appendix

In the appendix, we show a detailed description of the Evolution class in both
C++ and python. A detailed description of the Cosmo class as well as a basic
introduction of the RK methods used can be found in ref. [25].

##  Appendix A Evolution class in C++

The Evolution class uses the default constructors, so instances are declared
trivially as shown in Section 3. It has only one public member function,
solveNSC, with signature

⬇

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

2 bool Evolution<LD,Solver,Method>::solveNSC(const LD &TEND, const LD &c,
const LD &Ti,

3 const LD &ratio, const LD &TSTOP, const LD &umax, Cosmo<LD> *plasma,

4 const solverArgs<LD> &args={})

All the template and required input arguments are summarised in Table 1. This
function returns false if the integration stops above TSTOP and true
otherwise, as a way of providing a simple consistency check. The optional
argument, args, is an instance of a struct defined as

⬇

1 template<class LD>

2 struct solverArgs{

3 LD initial_step_size=1e-2, minimum_step_size=1e-8, maximum_step_size=1e-2;

4 LD absolute_tolerance=1e-8, relative_tolerance=1e-8;

5 LD beta=0.9, fac_max=1.2, fac_min=0.8;

6 unsigned int maximum_No_steps=10000000;

7 };

This helps us pass to solveNSC only the arguments we wish to change. For
example, we can choose absolute_tolerance=1e-11 and beta=0.8 by adding
{.absolute_tolerance=1e-11, .beta=0.8} in the place of args. This is called
aggregation in C++, and basically allows us to use named arguments. The caveat
is that we need to put the arguments in the order thy appear in the definition
of solverArgs. For example, using {.beta=0.8, .absolute_tolerance=1e-11} will
result in a compilation error.

The class Evolution has several member variables that store the results of the
solution of the system (2.11) and (2.12):

  1. 1.

std::vector<LD> u: integration steps of $u$.

  2. 2.

std::vector<LD> T: values of $T$ (in ${\rm GeV}$) that corresponds to every
step in u.

  3. 3.

std::vector<LD> rhoPhi: values of $\rho_{\Phi}$ (in ${\rm GeV}^{4}$) that
corresponds to every step in u.

  4. 4.

std::vector<LD> dT: local integration error of $T$ (in ${\rm GeV}$) at every
u.

  5. 5.

std::vector<LD> drhoPhi: local integration error of $\rho_{\Phi}$ (in ${\rm
GeV}^{4}$) at every u.

  6. 6.

LD TE1: value of $T_{\rm E1}{}$ (in ${\rm GeV}$). It is initialised at 0. If
$\rho_{\Phi}$ never reaches $\rho_{\rm R}$, This variable will be equal to
$T_{\rm i}{}$.

  7. 7.

LD TE2: value of $T_{\rm E2}{}$ (in ${\rm GeV}$). It is initialised at 0. If
$\rho_{\Phi}$ never reaches $\rho_{\rm R}$ or $c>4$, this variable will be
equal to $T_{\rm i}{}$.

  8. 8.

LD TD1: value of $T_{\rm D1}{}$ (in ${\rm GeV}$). It is initialised at 0. If
the energy injection rate never becomes more than $10\%$ of the free dilution
rate of $\rho_{\rm R}$, this variable obtains the value of $T_{\rm i}{}$.

  9. 9.

LD TD2: value of $T_{\rm D2}{}$ (in ${\rm GeV}$). It is initialised at 0. If
the energy injection rate never becomes more than $10\%$ of the free dilution
rate of $\rho_{\rm R}$, this variable obtains the value of $T_{\rm i}{}$.

  10. 10.

LD uE1: value of $u_{\rm E1}{}$. It is initialised at 0. If $\rho_{\Phi}$
never reaches $\rho_{\rm R}$, This variable will be equal to $0$.

  11. 11.

LD uE2: value of $u_{\rm E2}{}$. It is initialised at 0. If $\rho_{\Phi}$
never reaches $\rho_{\rm R}$ or $c>4$, this variable will be equal to $0$.

  12. 12.

LD uD1: value of $u_{\rm D1}{}$. It is initialised at 0. If the energy
injection rate never becomes more than $10\%$ of the free dilution rate of
$\rho_{\rm R}$, this variable obtains the value of $0$.

  13. 13.

LD uD2: value of $u_{\rm D2}{}$. It is initialised at 0. If the energy
injection rate never becomes more than $10\%$ of the free dilution rate of
$\rho_{\rm R}$, this variable obtains the value of $0$.

##  Appendix B Evolution class in python

The Evolution class is very similar to the C++ case. It uses the default
constructor, with instances declared as shown in Section 3. The solveNSC, has
signature

⬇

1 solveNSC(TEND, c, Ti, ratio, TSTOP, umax, plasma,

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

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

4 maximum_No_steps=10000000)

All the input arguments are summarised in Table 1. The difference between this
and the C++ version is that the plasma argument is an instance of the python
Cosmo class. This function returns a tuple. The first element of the returned
tuple is False if the integration stops above TSTOP and True otherwise, while
the second is the time it took to solve the system (in seconds).

The member variables are the same as in C++, with numpy arrays used in place
of the C++ vectors. However, solving the system using solveNSC does not fill
the arrays with the integration steps. The member function

⬇

1 getPoints()

Fills the arrays u, T, and rhoPhi with the integration steps.

In order to fill the arrays with the local errors we make use of the member
function

⬇

1 getErrors()

##  Appendix C Quick guide to the user input

In this section, we present tables with all the inputs the various functions
need.

User input for solving eqs. (2.11) and (2.12).
---
TEND |  The value of $T_{\rm end}{}$ in ${\rm GeV}$.
c |  Value of $c$ as defined through eq. (2.1).
Ti |  The value of $T_{\rm i}{}$ in ${\rm GeV}$.
ratio |  The value of $\rho_{\Phi}/\rho_{\rm R}$ at $T_{\rm i}{}$.
plasma |  A pointer to an instance of the Cosmo class.
TSTOP |  Once $T<$TSTOP, integration stops.
umax |  If $u>$umax the integration stops.
initial_stepsize |  Initial step-size of the solver. Default value: $10^{-2}$.
minimum_stepsize |  Lower limit of the step-size. Default value: $10^{-8}$.
maximum_stepsize |  Upper limit of the step-size. Default value: $10^{-2}$.
absolute_tolerance | Absolute tolerance of the RK solver. Default value: $10^{-8}$.
relative_tolerance | Relative tolerance of the RK solver. Default value: $10^{-8}$.
beta | Aggressiveness of the adaptation strategy. Default value: $0.9$.
fac_max, fac_min | The step-size does not change more than fac_max and less than fac_min within a trial step . Default values: $1.2$ and $0.8$, respectively.
maximum_No_steps | If integration needs more than maximum_No_steps integration stops. Default value: $10^{7}$.
Table 1: Table of the arguments of the
nsc::Evolution<LD,Solver,Method>::solveNSC method. User input for
interpolating RDOF.
---
path | Path to a file with $T$ (${\rm GeV}$), $h_{\rm eff}{}$, and $g_{\rm eff}{}$. The temperature has to monotonically increase.
minT | The RDOF will be interpolated starting from the closest temperature that satisfies $T_{\rm min}\leq$minT. Below this, both $h_{\rm eff}{}$ and $g_{\rm eff}{}$ will be assumed to be their value at $T_{\rm min}$.
maxT | The RDOF will be interpolated up to the closest temperature that satisfies $T_{\rm max}\geq$maxT. Above this, both $h_{\rm eff}{}$ and $g_{\rm eff}{}$ will be assumed to be their value at $T_{\rm max}$.
Table 2: Table of the arguments of the nsc::Cosmo<LD> class constructor.
Template arguments.
---
LD | This template argument is the numeric type that NSC++ will use. The preferred choice is long double. However, in many cases double can be used. Notice that the instance of the Cosmo class must use identical numeric type to the one that used by the instance of Evolution, in order to be able to run the solveNSC method.
Solver | This is the second template argument of the nsc::Evolution<LD,Solver,Method> class. The available choices are $1$ for Rosenbrock method, and $2$ for explicit RK method.
Method | The third template argument of the Evolution class. Its value depends on the choice of Solver; For Solver==1, Method can be RODASPR2<LD> (fourth order) [32], ROS34PW2<LD> (third order) [36], GRK4A<LD> or GRK4T<LD> (fourth order) [37]. For Solver==2, there are a few options, but the only choice that works in most cases is DormandPrince<LD> (fourth order) [38]. The Method classes need a template argument, LD, which must be the same as the first template argument of the nsc::Evolution<LD,Solver,Method> class. If one defines their own Butcher table, they would have to follow their definitions and assumptions.
Table 3: Template arguments of the various NSC++ classes. User compile-time
options. Variables in the various Definitions.mk files.
---
rootDir | The relative path of root directory of NSC++. Relevant only when compiling using make. Available in all Definitions.mk.
LONG | long for long double or empty for double. This is defines a macro in the source files of the various C++ examples. Available in Definitions.mk inside the various subdirectories of NSCpp/UserSpace/Cpp.
LONGpy | long or empty. Same as LONG, applies to the python modules. Available in NSCpp/Definitions.mk.
SOLVER | In order to use a Rosenbrock method SOLVER=$1$. For explicit RK method, SOLVER=$2$. This defines a macro that is passes as the second template argument of nsc::Evolution<LD,Solver,Method>. The corresponding variable in NSCpp/Definitions.mk applies to the python modules. The variable in NSCpp/UserSpace/Cpp/NSC/Definitions.mk applies to the example in the same directory.
METHOD | Depending on the solver, this variable should name one of its available methods. For SOLVER=$1$, METHOD=RODASPR2(fourth order) or ROS34PW2(third order). For SOLVER=$2$, METHOD=DormandPrince (seventh order). There is a macro (METHOD) used by the shared library NSCpp/lib/libNSC.so. The corresponding variable in NSCpp/Definitions.mk applies to the python modules. The variable in NSCpp/UserSpace/Cpp/NSC/Definitions.mk applies to the example in that directory.
Compiler options
CC | The preferred C++ compiler (g++ by default). Corresponding variable in all Definitions.mk files.
OPT | Available options are OPT=O1, O2, O3 (default). This variable defines the optimization level of the compiler. The variable can be changed in all Definitions.mk files. In the root directory of NSC++, the optimization level applies to the python modules (i.e. the shared libraries), while in the subdirectories of NSCpp/UserSpace/Cpp it only applies to example inside them.
Table 4: User compile-time input and options. These are available in the
various Definitions.mk files, which are used when compiling using make.

## References

  * [1] J. McDonald, WIMP Densities in Decaying Particle Dominated Cosmology, Phys. Rev. D 43 (1991) 1063–1068.

  * [2] F. D’Eramo, N. Fernandez, and S. Profumo, When the Universe Expands Too Fast: Relentless Dark Matter, JCAP 05 (2017) 012, [arXiv:1703.04793].
  * [3] K. Redmond and A. L. Erickcek, New Constraints on Dark Matter Production during Kination, Phys. Rev. D 96 (2017), no. 4 043511, [arXiv:1704.01056].
  * [4] F. D’Eramo, N. Fernandez, and S. Profumo, Dark Matter Freeze-in Production in Fast-Expanding Universes, JCAP 02 (2018) 046, [arXiv:1712.07453].
  * [5] N. Bernal, J. Rubio, and H. Veermäe, Boosting Ultraviolet Freeze-in in NO Models, JCAP 06 (2020) 047, [arXiv:2004.13706].
  * [6] P. Arias, D. Karamitros, and L. Roszkowski, Frozen-in fermionic singlet dark matter in non-standard cosmology with a decaying fluid, JCAP 05 (2021) 041, [arXiv:2012.07202].
  * [7] P. Arias, N. Bernal, D. Karamitros, C. Maldonado, L. Roszkowski, and M. Venegas, New opportunities for axion dark matter searches in nonstandard cosmological models, JCAP 11 (2021) 003, [arXiv:2107.13588].
  * [8] B. Barman, P. Ghosh, F. S. Queiroz, and A. K. Saha, Scalar multiplet dark matter in a fast expanding Universe: Resurrection of the desert region, Phys. Rev. D 104 (2021), no. 1 015040, [arXiv:2101.10175].
  * [9] K. R. Dienes, L. Heurtier, F. Huang, D. Kim, T. M. P. Tait, and B. Thomas, Stasis in an expanding universe: A recipe for stable mixed-component cosmological eras, Phys. Rev. D 105 (2022), no. 2 023530, [arXiv:2111.04753].
  * [10] A. Banerjee and D. Chowdhury, Fingerprints of freeze-in dark matter in an early matter-dominated era, SciPost Phys. 13 (2022), no. 2 022, [arXiv:2204.03670].
  * [11] E. Hardy, Higgs portal dark matter in non-standard cosmological histories, JHEP 06 (2018) 043, [arXiv:1804.06783].
  * [12] N. Bernal, C. Cosme, T. Tenkanen, and V. Vaskonen, Scalar singlet dark matter in non-standard cosmologies, Eur. Phys. J. C 79 (2019), no. 1 30, [arXiv:1806.11122].
  * [13] P. Arias, N. Bernal, A. Herrera, and C. Maldonado, Reconstructing Non-standard Cosmologies with Dark Matter, JCAP 10 (2019) 047, [arXiv:1906.04183].
  * [14] R. Allahverdi and J. K. Osiński, Freeze-in Production of Dark Matter Prior to Early Matter Domination, Phys. Rev. D 101 (2020), no. 6 063503, [arXiv:1909.01457].
  * [15] N. Bernal, F. Elahi, C. Maldonado, and J. Unwin, Ultraviolet Freeze-in and Non-Standard Cosmologies, JCAP 11 (2019) 026, [arXiv:1909.07992].
  * [16] C. Cosme, M. Dutra, T. Ma, Y. Wu, and L. Yang, Neutrino portal to FIMP dark matter with an early matter era, JHEP 03 (2021) 026, [arXiv:2003.01723].
  * [17] A. Vilenkin and L. H. Ford, Gravitational Effects upon Cosmological Phase Transitions, Phys. Rev. D 26 (1982) 1231.
  * [18] G. D. Coughlan, W. Fischler, E. W. Kolb, S. Raby, and G. G. Ross, Cosmological Problems for the Polonyi Potential, Phys. Lett. B 131 (1983) 59–64.
  * [19] B. Ratra and P. J. E. Peebles, Cosmological Consequences of a Rolling Homogeneous Scalar Field, Phys. Rev. D 37 (1988) 3406.
  * [20] G. F. Giudice, E. W. Kolb, and A. Riotto, Largest temperature of the radiation era and its cosmological implications, Phys. Rev. D 64 (2001) 023508, [hep-ph/0005123].
  * [21] C. L. Gardner, Quintessence and the transition to an accelerating universe, Nucl. Phys. B 707 (2005) 278–300, [astro-ph/0407604].
  * [22] I. Dalianis and Y. Watanabe, Probing the BSM physics with CMB precision cosmology: an application to supersymmetry, JHEP 02 (2018) 118, [arXiv:1801.05736].
  * [23] S. Tsujikawa, Quintessence: A Review, Class. Quant. Grav. 30 (2013) 214003, [arXiv:1304.1961].
  * [24] R. Allahverdi et al., The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe, arXiv:2006.16182.
  * [25] D. Karamitros, MiMeS: Misalignment mechanism solver, Comput. Phys. Commun. 275 (2022) 108311, [arXiv:2110.12253].
  * [26] M. Dutra and Y. Wu, EvoEMD: cosmic evolution with an early matter-dominated era, arXiv:2111.15665.
  * [27] D. Karamitros, NaBBODES: Not a Black Box Ordinary Differential Equation Solver in C++, 2019.
  * [28] D. Karamitros, SimpleSplines: A header-only library for linear and cubic spline interpolation in C++, 2021\.
  * [29] K. Saikawa and S. Shirai, Precise WIMP Dark Matter Abundance and Standard Model Thermodynamics, JCAP 08 (2020) 011, [arXiv:2005.03544].
  * [30] Planck Collaboration, Y. Akrami et al., Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641 (2020) A10, [arXiv:1807.06211].
  * [31] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585 (Sept., 2020) 357–362.
  * [32] J. Rang and L. Angermann, New rosenbrock w-methods of order 3 for partial differential algebraic equations of index 1, BIT Numerical Mathematics 45 (2005) 761–787.
  * [33] A. Ahmed, B. Grzadkowski, and A. Socha, Implications of time-dependent inflaton decay on reheating and dark matter production, Phys. Lett. B 831 (2022) 137201, [arXiv:2111.06065].
  * [34] B. Barman, N. Bernal, Y. Xu, and O. Zapata, Ultraviolet freeze-in with a time-dependent inflaton decay, JCAP 07 (2022), no. 07 019, [arXiv:2202.12906].
  * [35] P. Arias, N. Bernal, J. K. Osiński, and L. Roszkowski, Dark Matter Axions in the Early Universe with a Period of Increasing Temperature, arXiv:2207.07677.
  * [36] J. Rang, Improved traditional rosenbrock–wanner methods for stiff odes and daes, Journal of Computational and Applied Mathematics 286 (2015) 128–144.
  * [37] P. Rentrop and P. Kaps, Generalized Runge-Kutta Methods of Order Four with Stepsize Control for Stiff Ordinary Differential Equations., Numerische Mathematik 33 (1979) 55–68.
  * [38] J. Dormand and P. Prince, A family of embedded runge-kutta formulae, Journal of Computational and Applied Mathematics 6 (1980), no. 1 19–26.