CLASS MANUAL
Functions
thermodynamics.c File Reference
#include "thermodynamics.h"
#include "history.h"
#include "hyrectools.h"
#include "helium.h"
#include "wrap_hyrec.h"
+ Include dependency graph for thermodynamics.c:

Functions

int thermodynamics_at_z (struct background *pba, struct thermodynamics *pth, double z, enum interpolation_method inter_mode, int *last_index, double *pvecback, double *pvecthermo)
 
int thermodynamics_init (struct precision *ppr, struct background *pba, struct thermodynamics *pth)
 
int thermodynamics_free (struct thermodynamics *pth)
 
int thermodynamics_helium_from_bbn (struct precision *ppr, struct background *pba, struct thermodynamics *pth)
 
int thermodynamics_checks (struct precision *ppr, struct background *pba, struct thermodynamics *pth)
 
int thermodynamics_workspace_init (struct precision *ppr, struct background *pba, struct thermodynamics *pth, struct thermo_workspace *ptw)
 
int thermodynamics_indices (struct background *pba, struct thermodynamics *pth, struct thermo_workspace *ptw)
 
int thermodynamics_lists (struct precision *ppr, struct background *pba, struct thermodynamics *pth, struct thermo_workspace *ptw)
 
int thermodynamics_set_parameters_reionization (struct precision *ppr, struct background *pba, struct thermodynamics *pth, struct thermo_reionization_parameters *preio)
 
int thermodynamics_solve (struct precision *ppr, struct background *pba, struct thermodynamics *pth, struct thermo_workspace *ptw, double *pvecback)
 
int thermodynamics_calculate_remaining_quantities (struct precision *ppr, struct background *pba, struct thermodynamics *pth, double *pvecback)
 
int thermodynamics_output_summary (struct background *pba, struct thermodynamics *pth)
 
int thermodynamics_workspace_free (struct thermodynamics *pth, struct thermo_workspace *ptw)
 
int thermodynamics_vector_init (struct precision *ppr, struct background *pba, struct thermodynamics *pth, double mz, struct thermo_workspace *ptw)
 
int thermodynamics_reionization_evolve_with_tau (struct thermodynamics_parameters_and_workspace *ptpaw, double mz_ini, double mz_end, double *mz_output, int mz_size)
 
int thermodynamics_derivs (double mz, double *y, double *dy, void *parameters_and_workspace, ErrorMsg error_message)
 
int thermodynamics_timescale (double mz, void *thermo_parameters_and_workspace, double *timescale, ErrorMsg error_message)
 
int thermodynamics_sources (double mz, double *y, double *dy, int index_z, void *thermo_parameters_and_workspace, ErrorMsg error_message)
 
int thermodynamics_reionization_get_tau (struct precision *ppr, struct background *pba, struct thermodynamics *pth, struct thermo_workspace *ptw)
 
int thermodynamics_vector_free (struct thermo_vector *tv)
 
int thermodynamics_calculate_conformal_drag_time (struct background *pba, struct thermodynamics *pth, double *pvecback)
 
int thermodynamics_calculate_damping_scale (struct background *pba, struct thermodynamics *pth, double *pvecback)
 
int thermodynamics_calculate_opticals (struct precision *ppr, struct thermodynamics *pth)
 
int thermodynamics_calculate_idm_and_idr_quantities (struct precision *ppr, struct background *pba, struct thermodynamics *pth, double *pvecback)
 
int thermodynamics_calculate_recombination_quantities (struct precision *ppr, struct background *pba, struct thermodynamics *pth, double *pvecback)
 
int thermodynamics_calculate_drag_quantities (struct precision *ppr, struct background *pba, struct thermodynamics *pth, double *pvecback)
 
int thermodynamics_ionization_fractions (double z, double *y, struct background *pba, struct thermodynamics *pth, struct thermo_workspace *ptw, int current_ap)
 
int thermodynamics_reionization_function (double z, struct thermodynamics *pth, struct thermo_reionization_parameters *preio, double *x)
 
int thermodynamics_output_titles (struct background *pba, struct thermodynamics *pth, char titles[_MAXTITLESTRINGLENGTH_])
 
int thermodynamics_output_data (struct background *pba, struct thermodynamics *pth, int number_of_titles, double *data)
 
int thermodynamics_idm_quantities (struct background *pba, double z, double *y, double *dy, struct thermodynamics *pth, struct thermo_workspace *ptw, double *pvecback)
 
int thermodynamics_obtain_z_ini (struct precision *ppr, struct background *pba, struct thermodynamics *pth, struct thermo_workspace *ptw)
 
int thermodynamics_idm_initial_temperature (struct background *pba, struct thermodynamics *pth, double z_ini, struct thermo_diffeq_workspace *ptdw)
 

Detailed Description

Documented thermodynamics module

  • Julien Lesgourgues, 6.09.2010
  • Restructured by Nils Schoeneberg and Matteo Lucca, 27.02.2019
  • Evolver implementation by Daniel Meinert, spring 2019

Deals with the thermodynamical evolution. This module has two purposes:

  • at the beginning, to initialize the thermodynamics, i.e. to integrate the thermodynamical equations, and store all thermodynamical quantities as a function of redshift inside an interpolation table.
  • to provide a routine which allow other modules to evaluate any thermodynamical quantities at a given redshift value (by interpolating within the interpolation table).

The most important differential equations to compute the free electron fraction x at each step are provided either by the HyRec 2020 or RecFastCLASS code, located in the external/ directory. The thermodynamics module integrates these equations using the generic integrator (which can be set to ndf15, rkck4, etc.) The HyRec and RecFastCLASS algorithms are used and called in the same way by this module.

In summary, the following functions can be called from other modules:

  1. thermodynamics_init at the beginning (but after background_init)
  2. thermodynamics_at_z at any later time
  3. thermodynamics_free at the end, when no more calls to thermodynamics_at_z are needed

Function Documentation

◆ thermodynamics_at_z()

int thermodynamics_at_z ( struct background *  pba,
struct thermodynamics *  pth,
double  z,
enum interpolation_method  inter_mode,
int *  last_index,
double *  pvecback,
double *  pvecthermo 
)

Thermodynamics quantities at given redshift z. Evaluates all thermodynamics quantities at a given value of the redshift by reading the pre-computed table and interpolating.

Parameters
pbaInput: pointer to background structure
pthInput: pointer to the thermodynamics structure (containing pre-computed table)
zInput: redshift
inter_modeInput: interpolation mode (normal or growing_closeby)
last_indexInput/Output: index of the previous/current point in the interpolation array (input only for closeby mode, output for both)
pvecbackInput: vector of background quantities (used only in case z>z_initial for getting ddkappa and dddkappa; in that case, should be already allocated (!) and filled (!), with format short_info or larger; in other cases, will be ignored)
pvecthermoOutput: vector of thermodynamics quantities (assumed to be already allocated)
Returns
the error status

Summary:

  • define local variables
  • interpolate in table with array_interpolate_spline (normal mode) or array_interpolate_spline_growing_closeby (closeby mode)

◆ thermodynamics_init()

int thermodynamics_init ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth 
)

Initialize the thermodynamics structure, and in particular the thermodynamics interpolation table.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput/Output: pointer to initialized thermodynamics structure
Returns
the error status

Summary:

  • define local variables
  • update the user about which recombination code is being run
  • set flag for varying constants
  • compute and check primordial Helium mass fraction rho_He/(rho_H+rho_He)
  • infer primordial helium-to-hydrogen nucleon ratio n_He/n_H It is calculated via n_He/n_H = rho_He/(m_He/m_H * rho_H) = YHe * rho_b / (m_He/m_H * (1-YHe) rho_b) = YHe / (m_He/m_H * (1-YHe))
  • infer number of hydrogen nuclei today in m**-3
  • test whether all parameters are in the correct regime
  • allocate and assign all temporary structures and indices
  • initialize injection struct (not temporary)
  • assign reionisation parameters
  • solve recombination and reionization and store values of $ z, x_e, d \kappa / d \tau, T_b, c_b^2 $
  • the differential equation system is now completely solved
  • fill missing columns (quantities not computed during the differential evolution but related)
  • write information on thermal history in standard output
  • free workspace and local variables

◆ thermodynamics_free()

int thermodynamics_free ( struct thermodynamics *  pth)

Free all memory space allocated by thermodynamics_init.

Parameters
pthInput/Output: pointer to thermodynamics structure (to be freed)
Returns
the error status

◆ thermodynamics_helium_from_bbn()

int thermodynamics_helium_from_bbn ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth 
)

Infer the primordial helium mass fraction from standard BBN calculations, as a function of the baryon density and expansion rate during BBN.

This module is simpler then the one used in arXiv:0712.2826 because it neglects the impact of a possible significant chemical potentials for electron neutrinos. The full code with xi_nu_e could be introduced here later.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput/Output: pointer to initialized thermodynamics structure
Returns
the error status

Summary:

Define local variables

  • Infer effective number of neutrinos at the time of BBN
  • We randomly choose 0.1 MeV to be the temperature of BBN
  • compute Delta N_eff as defined in bbn file, i.e. $ \Delta N_{eff}=0$ means $ N_{eff}=3.046$. Note that even if 3.044 is a better default value, we must keep 3.046 here as long as the BBN file we are using has been computed assuming 3.046.
  • spline in one dimension (along deltaN)
  • interpolate in one dimension (along deltaN)
  • spline in remaining dimension (along omegab)
  • interpolate in remaining dimension (along omegab)
  • Take into account impact of varying alpha on helium fraction
  • deallocate arrays

◆ thermodynamics_checks()

int thermodynamics_checks ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth 
)

Check the thermodynamics structure parameters for bounds and critical values.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput: pointer to initialized thermodynamics structure
Returns
the error status

Summary:

  • check BBN Y_He fracion
  • tests in order to prevent divisions by zero
  • test initial condition for recombination

◆ thermodynamics_workspace_init()

int thermodynamics_workspace_init ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
struct thermo_workspace *  ptw 
)

Initialize the thermodynamics workspace.

The workspace contains the arrays used for solving differential equations (dubbed thermo_diffeq_workspace), and storing all approximations, reionization parameters, heating parameters.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput: pointer to the thermodynamics structure
ptwInput/Output: pointer to thermodynamics workspace
Returns
the error status

Summary:

Define local variables

  • number of z values
  • relevant cosmological parameters
  • relevant constants
  • Allocate and initialize differential equation workspace
  • define approximations
  • store all ending redshifts for each approximation
  • Rescale these redshifts in case of varying fundamental constants
  • store smoothing deltas for transitions at the beginning of each aproximation
  • Allocate reionisation parameter workspace

◆ thermodynamics_indices()

int thermodynamics_indices ( struct background *  pba,
struct thermodynamics *  pth,
struct thermo_workspace *  ptw 
)

Assign value to each relevant index in vectors of thermodynamical quantities, and the reionization parameters

Parameters
pbaInput: pointer to background structure
pthInput/Output: pointer to thermodynamics structure
ptwInput/Output: pointer to thermo workspace
Returns
the error status

Summary:

  • define local variables
  • initialization of all indices and flags in thermodynamics structure
  • initialization of all indices of parameters of reionization function

◆ thermodynamics_lists()

int thermodynamics_lists ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
struct thermo_workspace *  ptw 
)

Initialize the lists (of redshift, tau, etc.) of the thermodynamics struct

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput/Output: pointer to thermodynamics structure
ptwInput: pointer to thermo workspace
Returns
the error status

Summary:

Define local variables

  • allocate tables
  • define time sampling
  • store initial value of conformal time in the structure

◆ thermodynamics_set_parameters_reionization()

int thermodynamics_set_parameters_reionization ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
struct thermo_reionization_parameters *  preio 
)

This routine initializes reionization_parameters for the chosen scheme of reionization function.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput: pointer to the thermodynamics structure
preioInput/Output: pointer to the reionization parameters structure
Returns
the error status

Summary:

Define local variables

  • allocate the vector of parameters defining the function $ X_e(z) $
  • (a) no reionization
  • (b) if reionization implemented like in CAMB, or half tanh like in 1209.0247
  • --> set values of these parameters, excepted those depending on the reionization redshift
  • --> if reionization redshift given as an input, initialize the remaining values
  • --> if reionization optical depth given as an input, find reionization redshift by bisection and initialize the remaining values
  • (c) if reionization implemented with reio_bins_tanh scheme
  • (d) if reionization implemented with reio_many_tanh scheme
  • (e) if reionization implemented with reio_inter scheme

◆ thermodynamics_solve()

int thermodynamics_solve ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
struct thermo_workspace *  ptw,
double *  pvecback 
)

Integrate thermodynamics with your favorite recombination code. The default options are HyRec and RecFastCLASS.

Integrate thermodynamics with HyRec or Recfast, allocate and fill part of the thermodynamics interpolation table (the rest is filled in thermodynamics_calculate_remaining_quantitie).

Version modified by Daniel Meinert and Nils Schoeneberg to use the ndf15 evolver or any other evolver inherent to CLASS, modified again by Nils Schoeneberg to use wrappers.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput/Output: pointer to thermodynamics structure where results are stored
ptwInput: pointer to thermo_workspace structure used to communicate with generic evolver
pvecbackInput: pointer to an allocated (but empty) vector of background variables
Returns
the error status

Integrate thermodynamics with your favorite recombination code. The default options are HyRec and Recfast.

Summary:

  • define local variables
  • choose evolver
  • define the fields of the 'thermodynamics parameter and workspace' structure
  • define time sampling: create a local array of minus z values called mz (from mz=-zinitial growing towards mz=0)
  • define intervals for each approximation scheme
  • loop over intervals over which approximation scheme is uniform. For each interval:
  • --> (a) fix current approximation scheme.
  • --> (b) define the vector of quantities to be integrated over. If the current interval starts from the initial time zinitial, fill the vector with initial conditions. If it starts from an approximation switching point, redistribute correctly the values from the previous to the new vector. For both RECFAST and HYREC, the vector consists of Tmat, x_H, x_He, + others for exotic models
  • --> (c1) If we have the optical depth tau_reio as input the last evolver step (reionization approximation) is done separately in a loop, to find the approximate redshift of reionization given the input of tau_reio, using a bisection method. This is similar to the general CLASS shooting method, but doing this step here is more advantageous since we only need to do repeatedly the last approximation step of the integration, instead of the full background and thermodynamics module

--> (c2) otherwise, just integrate quantities over the current interval.

  • Compute reionization optical depth, if not supplied as input parameter
  • free quantities allocated at the beginning of the routine

◆ thermodynamics_calculate_remaining_quantities()

int thermodynamics_calculate_remaining_quantities ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
double *  pvecback 
)

Calculate those thermodynamics quantities which are not inside of the thermodynamics table already.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput/Output: pointer to initialized thermodynamics structure
pvecbackInput: pointer to some allocated pvecback
Returns
the error status

Summary:

  • fill tables of second derivatives with respect to z (in view of spline interpolation)

◆ thermodynamics_output_summary()

int thermodynamics_output_summary ( struct background *  pba,
struct thermodynamics *  pth 
)

In verbose mode, print basic information on the thermal history

Parameters
pbaInput: pointer to background structure
pthInput/Output: pointer to initialized thermodynamics structure
Returns
the error status

Summary:

Define local variables

◆ thermodynamics_workspace_free()

int thermodynamics_workspace_free ( struct thermodynamics *  pth,
struct thermo_workspace *  ptw 
)

Free the thermo_workspace structure (with the exception of the thermo_vector '->ptv' field, which is freed separately in thermo_vector_free).

Parameters
pthInput: pointer to initialized thermodynamics structure
ptwInput: pointer to perturbations_workspace structure to be freed
Returns
the error status

◆ thermodynamics_vector_init()

int thermodynamics_vector_init ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
double  mz,
struct thermo_workspace *  ptw 
)

Initialize the field '->ptv' of a thermo_diffeq_workspace structure, which is a thermo_vector structure. This structure contains indices and values of all quantities which need to be integrated with respect to time (and only them: quantities fixed analytically or obeying constraint equations are NOT included in this vector).

The routine sets and allocates the vector y, dy and used_in_output with the right size depending on the current approximation scheme stored in the workspace. Moreover the initial conditions for each approximation scheme are calculated and set correctly.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput: pointer to the thermodynamics structure
mzInput: negative redshift
ptwInput/Output: pointer to thermodynamics workspace
Returns
the error status

Summary:

Define local variables

◆ thermodynamics_reionization_evolve_with_tau()

int thermodynamics_reionization_evolve_with_tau ( struct thermodynamics_parameters_and_workspace *  ptpaw,
double  mz_ini,
double  mz_end,
double *  mz_output,
int  mz_size 
)

If the input for reionization is tau_reio, thermodynamics_solve() calls this function instead of the evolver for dealing with the last era (the reionization era).

Instead of computing the evolution of quantities during reionization for a fixed z_reio, as the evolver would do, this function finds z_reio by bisection. First we make an initial guess for z_reio with reionization_z_start_max and then find a z_reio which leads to the given tau_reio (in the range of tolerance reionization_optical_depth_tol).

Parameters
ptpawInput: pointer to parameters and workspace
mz_iniInput: initial redshift
mz_endInput: ending redshift
mz_outputInput: pointer to redshift array at which output should be written
mz_sizeInput: number of redshift values in this array
Returns
the error status

Summary:

Define local variables

  • Remame fields to avoid heavy notations
  • Choose evolver
  • ptvs will be a pointer towards the same thermo vector that was used in the previous approximation schemes; it contains values that will serve here to set initial conditions.
  • ptv is a pointer towards a whole new thermo vector used for the calculations in the bisection, that we must allocate and initialize
  • Initialize the values of the temporary vector
  • Evolve quantities through reionization assuming upper value of z_reio
  • Restore initial conditions
  • Evolve quantities through reionization assuming lower value of z_reio
  • Restore initial conditions
  • Evolve quantities through reionization, trying intermediate values of z_reio by bisection
  • Store the ionization redshift in the thermodynamics structure
  • Free tempeoraty thermo vector

◆ thermodynamics_derivs()

int thermodynamics_derivs ( double  mz,
double *  y,
double *  dy,
void *  parameters_and_workspace,
ErrorMsg  error_message 
)

Subroutine evaluating the derivative of thermodynamical quantities with respect to negative redshift mz=-z.

Automatically recognizes the current approximation interval and computes the derivatives for this interval of the vector y, which contains (Tmat, x_H, x_He) + others for exotic models.

Derivatives are obtained either by calling either HyRec 2020 (Lee and Ali-Haimoud 2020, 2007.14114) or RecFastCLASS (that is, RecFast version 1.5, modified by Daniel Meinert and Nils Schoeneberg for better precision and smoothness at early times). See credits and licences in the wrappers (in external/...)

This is one of the few functions in the code which are passed to the generic_evolver routine. Since generic_evolver should work with functions passed from various modules, the format of the arguments is a bit special:

  • fixed parameters and workspaces are passed through a generic pointer. Here, this pointer contains the precision, background and thermodynamics structures, plus a background vector, but generic_evolver doesn't know its precise structure.
  • the error management is a bit special: errors are not written as usual to pth->error_message, but to a generic error_message passed in the list of arguments.
Parameters
mzInput: negative redshift mz = -z
yInput: vector of variable to integrate
dyOutput: its derivative (already allocated)
parameters_and_workspaceInput: pointer to fixed parameters (e.g. indices) and workspace (already allocated)
error_messageOutput: error message

Summary:

Define local variables

  • Rename structure fields (just to avoid heavy notations)
  • Get background/thermo quantities in this point

Set Tmat from the evolver (it is always evolved) and store it in the workspace.

  • The input vector y contains thermodynamic variables like (Tmat, x_H,x_He). The goal of this function is: 1) Depending on the chosen code and current approximation, to use either analytical approximations or the vector y to calculate x_e; 2) To compute re-ionization effects on x_e; The output of this function is stored in the workspace ptdw
  • If needed, calculate heating effects (i.e. any possible energy deposition rates affecting the evolution equations for x and Tmat)
  • Derivative of the ionization fractions

--> use Recfast or HyRec to get the derivatives d(x_H)/dz and d(x_He)/dz, and store the result directly in the vector dy. This gives the derivative of the ionization fractions from recombination only (not from reionization). Of course, the full treatment would involve the actual evolution equations for x_H and x_He during reionization, but these are not yet fully implemented.

< Pass value of fsR = relative alpha (fine-structure)

< Pass value of meR = relative m_e (effective eletron mass)

  • Derivative of the matter temperature (relevant for both Recfast and HyRec cases)
  • Calculate quantities for interacting dark matter
  • Derivative of the Dark Matter temperature
  • If we have extreme heatings, recombination does not fully happen and/or re-ionization happens before a redshift of reionization_z_start_max (default = 50). We want to catch this unphysical regime, because it would lead to further errors (and/or unphysical calculations) within our recombination codes
  • invert all derivatives (because the evolver evolves with -z, not with +z)

◆ thermodynamics_timescale()

int thermodynamics_timescale ( double  mz,
void *  thermo_parameters_and_workspace,
double *  timescale,
ErrorMsg  error_message 
)

This function is relevant for the rk evolver, not ndf15. It estimates a timescale 'delta z' over which quantitites vary. The rk evolver divides large intervals in steps given by this timescale multiplied by ppr->thermo_integration_stepsize.

This is one of the few functions in the code which is passed to the generic_evolver routine. Since generic_evolver should work with functions passed from various modules, the format of the arguments is a bit special:

  • fixed parameters and workspaces are passed through a generic pointer. generic_evolver doesn't know the content of this pointer.
  • the error management is a bit special: errors are not written as usual to pth->error_message, but to a generic error_message passed in the list of arguments.
Parameters
mzInput: minus the redshift
thermo_parameters_and_workspaceInput: pointer to parameters and workspace
timescaleOutput: pointer to the timescale
error_messageOutput: possible errors are written here
Returns
the error status

◆ thermodynamics_sources()

int thermodynamics_sources ( double  mz,
double *  y,
double *  dy,
int  index_z,
void *  thermo_parameters_and_workspace,
ErrorMsg  error_message 
)

This function is passed to the generic evolver and is called whenever we want to store values for a given mz.

The ionization fraction is either computed within a call to thermodynamics_derivs(). Moreover there is an automatic smoothing enabled which smoothes out the the ionization_fraction after each approximation switch. This is also the place where HyRec is asked to evolve x(z) using its internal system of differential equations over the next range [z_i, z_i+1], and to store the result in a temporary table.

This is one of the few functions in the code which is passed to the generic_evolver routine. Since generic_evolver should work with functions passed from various modules, the format of the arguments is a bit special:

  • fixed parameters and workspaces are passed through a generic pointer. generic_evolver doesn't know the content of this pointer.
  • the error management is a bit special: errors are not written as usual to pth->error_message, but to a generic error_message passed in the list of arguments.

All quantities are computed by a simple call to thermodynamics_derivs, which computes all necessary quantities and stores them in the ptdw thermo_diffeq_workspace structure

Parameters
mzInput: negative redshift, belonging to array mz_output
yInput: vector of evolved thermodynamical quantities
dyInput: derivatives of this vector w.r.t redshift
index_zInput: index in the array mz_output
thermo_parameters_and_workspaceInput/Output: in input, all parameters needed by thermodynamics_derivs; in output, recombination table
error_messageOutput: error message
Returns
the error status

Summary:

Define local variables

  • Rename structure fields (just to avoid heavy notations)
  • Recalculate all quantities at this current redshift: we need at least pvecback, ptdw->x_reio, dy[ptv->index_ti_D_Tmat]
  • In the recfast case, we manually smooth the results a bit
  • Store the results in the table. Results are obtained in order of decreasing z, and stored in order of growing z

◆ thermodynamics_reionization_get_tau()

int thermodynamics_reionization_get_tau ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
struct thermo_workspace *  ptw 
)

Get the optical depth of reionization tau_reio for a given thermodynamical history.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput: pointer to the thermodynamics structure
ptwInput: pointer to thermodynamics workspace
Returns
the error status

Summary:

Define local variables

We are searching now for the start of reionization. This will be the time at which the optical depth tau_reio will be computed.

Note that the value reionization_parameters[index_reio_start] is only the start of the reionization function added manually, but not necessarily the total start of reionization. Reionization could be longer/shifted by energy injection.

The actual the definition of tau_reio is not unique and unambiguous. We defined it here to be the optical depth up to the time at which there is a global minimum in the free electron fraction. We search for this time by iterating over the thermodynamics table, in order to find the corresponding index_reio_start.

  • --> spline $ d \tau / dz $ with respect to z in view of integrating for optical depth between 0 and the just found starting index
  • --> integrate for optical depth

◆ thermodynamics_vector_free()

int thermodynamics_vector_free ( struct thermo_vector *  tv)

Free the thermo_vector structure, which is the '->ptv' field of the thermodynamics_differential_workspace ptdw structure

Parameters
tvInput: pointer to thermo_vector structure to be freed
Returns
the error status

◆ thermodynamics_calculate_conformal_drag_time()

int thermodynamics_calculate_conformal_drag_time ( struct background *  pba,
struct thermodynamics *  pth,
double *  pvecback 
)

Compute the baryon drag conformal time tau_d = [int_{tau_today}^{tau} dtau -dkappa_d/dtau]

Parameters
pbaInput: pointer to background structure
pthInput/Output: pointer to initialized thermodynamics structure
pvecbackInput: Initialized vector of background quantities
Returns
the error status

Summary:

Define local variables

  • compute minus the baryon drag interaction rate time, -dkappa_d/dtau = -[1/R * kappa'], with R = 3 rho_b / 4 rho_gamma, stored temporarily in column ddkappa
  • compute second derivative of this rate, -[1/R * kappa']'', stored temporarily in column dddkappa
  • compute tau_d = [int_{tau_today}^{tau} dtau -dkappa_d/dtau]

◆ thermodynamics_calculate_damping_scale()

int thermodynamics_calculate_damping_scale ( struct background *  pba,
struct thermodynamics *  pth,
double *  pvecback 
)

Compute the damping scale r_d = 2pi/k_d = 2pi * [int_{tau_ini}^{tau} dtau (1/kappa') 1/6 (R^2+16/15(1+R))/(1+R)^2]^1/2 = 2pi * [int_{tau_ini}^{tau} dtau (1/kappa') 1/6 (R^2/(1+R)+16/15)/(1+R)]^1/2

which is like in CosmoTherm (CT), but slightly different from Wayne Hu (WH)'s thesis eq. (5.59): The factor 16/15 in CT is 4/5 in WH, but 16/15 is taking more effects into account

Parameters
pbaInput: pointer to background structure
pthInput/Output: pointer to initialized thermodynamics structure
pvecbackInput: Initialized vector of background quantities
Returns
the error status

Summary:

Define local variables

◆ thermodynamics_calculate_opticals()

int thermodynamics_calculate_opticals ( struct precision *  ppr,
struct thermodynamics *  pth 
)

Calculate quantities relating to optical phenomena like kappa' and exp(-kappa) and the visibility function, optical depth (including dark matter photon interactions if necessary)

Parameters
pprInput: pointer to precision structure
pthInput/Output: pointer to thermodynamics structure
Returns
the error status

Summary:

Define local quantities

  • --> second derivative with respect to tau of dkappa (in view of of spline interpolation)
  • --> first derivative with respect to tau of dkappa (using spline interpolation)
  • --> compute -kappa = [int_{tau_today}^{tau} dtau dkappa/dtau], store temporarily in column "g"

if there is idm_g, calculate its optical contributions

  • --> second derivative with respect to tau of dmu (in view of of spline interpolation)
  • --> first derivative with respect to tau of dmu (using spline interpolation)
  • --> compute -mu = [int_{tau_today}^{tau} dtau dmu/dtau], store temporarily in column "exp_mu_idm_g"
  • --> compute visibility: $ g= (d \kappa/d \tau) e^{- \kappa} $
  • —> compute g
  • —> compute g' (the plus sign of the second term is correct, see def of -kappa in thermodynamics module!)
  • —> compute g''
  • —> store g
  • —> compute variation rate
  • smooth the rate (details of smoothing unimportant: only the order of magnitude of the rate matters)
  • --> derivatives of baryon sound speed (only computed if some non-minimal tight-coupling schemes is requested)
  • —> second derivative with respect to tau of cb2
  • —> first derivative with respect to tau of cb2 (using spline interpolation)

◆ thermodynamics_calculate_idm_and_idr_quantities()

int thermodynamics_calculate_idm_and_idr_quantities ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
double *  pvecback 
)

Calculate dark optical depths, idm_b interation rate, and other quantities relevant for idm and idr.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput/Output: pointer to initialized thermo structure
pvecbackInput: Initialized vector of background quantities
Returns
the error status

Summary:

  • Define local variables
  • Calculate optical functions tau_idm_dr, tau_idr, g_idm_dr
  • Find interacting dark radiation free-streaming time
  • Search after the above free-streaming time for idr free-streaming as well
  • --> second derivative with respect to tau of R_idm_b (in view of spline interpolation)
  • --> first derivative with respect to tau of R_idm_b (using spline interpolation)

◆ thermodynamics_calculate_recombination_quantities()

int thermodynamics_calculate_recombination_quantities ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
double *  pvecback 
)

Calculate various quantities at the time of recombination, as well as the time tau_cut at which visibility gets negligible and one can assume pure free-streaming.

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput/Output: pointer to initialized thermodynamics structure
pvecbackInput: pointer to some allocated pvecback
Returns
the error status

Summary:

Define local variables

  • find maximum of g
  • find conformal recombination time using background_tau_of_z
  • find damping scale at recombination (using linear interpolation)
  • find time (always after recombination) at which tau_c/tau falls below some threshold, defining tau_free_streaming
  • find time above which visibility falls below a given fraction of its maximum
  • find z_star (when optical depth kappa crosses one, using linear interpolation) and sound horizon at that time

◆ thermodynamics_calculate_drag_quantities()

int thermodynamics_calculate_drag_quantities ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
double *  pvecback 
)

Calculate various quantities at the time of ending of baryon drag (It is precisely where tau_d crosses one)

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput/Output: pointer to initialized thermodynamics structure
pvecbackInput: pointer to some allocated pvecback
Returns
the error status

Summary:

Define local variables

  • find baryon drag time (when tau_d crosses one, using linear interpolation) and sound horizon at that time

◆ thermodynamics_ionization_fractions()

int thermodynamics_ionization_fractions ( double  z,
double *  y,
struct background *  pba,
struct thermodynamics *  pth,
struct thermo_workspace *  ptw,
int  current_ap 
)

Compute ionization fractions with the RecFast of HyRec algorithm.

Compute ionization fractions using either the vector y or, in some approximation schemes, some analytic approximations. The output of this function is located in the ptw->ptdw workspace. We need to assign:

  • in the RecFast scheme only: – ptdw->x_H, ptdw-> x_He (all neglecting reionisation, which is accounted for later on);
  • in both schemes: – ptdw->x_noreio (neglecting reionisation); – ptdw->x_reio (if reionisation is going on; obtained by calling thermodynamics_reionization_function() and adding something to ptdw->x_noreio)
Parameters
zInput: redshift
yInput: vector of quantities to integrate with evolver
pthInput: pointer to thermodynamics structure
pbaInput: pointer to background structure
ptwInput/Output: pointer to thermo workspace. Contains output for x, ...
current_apInput: index of current approximation scheme
Returns
the error status

Summary:

Define local variables

  • Calculate x_noreio from Recfast/Hyrec in each approximation regime. Store the results in the workspace.
  • --> For credits, see external/wrap_recfast.c
  • --> first regime: H and Helium fully ionized
  • --> second regime: first Helium recombination (analytic approximation)
  • --> third regime: first Helium recombination finished, H and Helium fully ionized
  • --> fourth regime: second Helium recombination starts (analytic approximation)
  • --> fifth regime: Hydrogen recombination starts (analytic approximation) while Helium recombination continues (full equation)
  • --> sixth regime: full Hydrogen and Helium equations
  • --> seventh regime: calculate x_noreio during reionization (i.e. x before taking reionisation into account)
  • If z is during reionization, also calculate the reionized x

◆ thermodynamics_reionization_function()

int thermodynamics_reionization_function ( double  z,
struct thermodynamics *  pth,
struct thermo_reionization_parameters *  preio,
double *  x 
)

This subroutine contains the reionization function $ X_e(z) $ (one for each scheme) and gives x for a given z.

Parameters
zInput: redshift
pthInput: pointer to thermodynamics structure, to know which scheme is used
preioInput: pointer to reionization parameters of the function $ X_e(z) $
xOutput: $ X_e(z) $

Summary:

  • define local variables
  • no reionization means nothing to be added to xe_before
  • implementation of ionization function similar to the one in CAMB
  • --> case z > z_reio_start
  • --> case z < z_reio_start: hydrogen contribution (tanh of complicated argument)
  • --> case z < z_reio_start: helium contribution (tanh of simpler argument)
  • implementation of half-tangent like in 1209.0247
  • --> case z > z_reio_start
  • --> case z < z_reio_start: hydrogen contribution (tanh of complicated argument)
  • implementation of binned ionization function similar to astro-ph/0606552
  • --> case z > z_reio_start
  • implementation of many tanh jumps
  • --> case z > z_reio_start
  • implementation of reio_inter
  • --> case z > z_reio_start

◆ thermodynamics_output_titles()

int thermodynamics_output_titles ( struct background *  pba,
struct thermodynamics *  pth,
char  titles[_MAXTITLESTRINGLENGTH_] 
)

Function for formatting the titles to be output

Parameters
pbaInput: pointer to background structure
pthInput: pointer to the thermodynamics structure
titlesInput: titles string containing all titles
Returns
the error status

◆ thermodynamics_output_data()

int thermodynamics_output_data ( struct background *  pba,
struct thermodynamics *  pth,
int  number_of_titles,
double *  data 
)

Output the data for the output into files

Parameters
pbaInput: pointer to background structure
pthInput: pointer to the thermodynamics structure
number_of_titlesInput: number of titles
dataInput: pointer to data file
Returns
the error status

◆ thermodynamics_idm_quantities()

int thermodynamics_idm_quantities ( struct background *  pba,
double  z,
double *  y,
double *  dy,
struct thermodynamics *  pth,
struct thermo_workspace *  ptw,
double *  pvecback 
)

This routine computes the quantities connected to interacting dark matter with photons, baryons & dark radiation (idm), and interacting dark radiation (idr)

Parameters
pbaInput: pointer to background structure
zInput: redshift
yInput: vector of evolver quantities
dyInput: derivative of this vector
pthInput: pointer to thermodynamics structure
ptwInput/Output: pointer to thermo workspace
pvecbackInput: vector of background quantities
Returns
the error status

Summary:

Define local variables

  • First deal with any required dark radiation
  • Now deal with any required dark matter (and its interactions)

◆ thermodynamics_obtain_z_ini()

int thermodynamics_obtain_z_ini ( struct precision *  ppr,
struct background *  pba,
struct thermodynamics *  pth,
struct thermo_workspace *  ptw 
)

Check if the initial integration time and spacing needs adjusting, for example for interacting dark matter

Parameters
pprInput: pointer to precision structure
pbaInput: pointer to background structure
pthInput: pointer to thermo structure
ptwInput/Output: pointer to thermo workspace
Returns
the error status

◆ thermodynamics_idm_initial_temperature()

int thermodynamics_idm_initial_temperature ( struct background *  pba,
struct thermodynamics *  pth,
double  z_ini,
struct thermo_diffeq_workspace *  ptdw 
)

Set the correct initial temperature for the idm species

Parameters
pbaInput: pointer to background structure
pthInput: pointer to thermo structure
z_iniInput: z_ini
ptdwInput/Output: pointer to thermo differential equation workspace
Returns
the error status
  • define local parameters
  • obtain necessary background information