In-depth usage
ECHO21 can be used to generate the thermal and ionization history of the intergalactic medium and hence, the cosmological global 21-cm signal. Addionally, one can use this code to study a simple analytical model of reionization and compute the CMB optical depth.
Single realization of 21-cm signal
There are just two steps to use ECHO21:
Give your choice of parameters
Run the solver
Thus, you first set up your cosmological and astrophysical parameters. These have to be supplied as a dictionary. Then specify your star formation rate density (SFRD) model, again as a dictionary. Once this is done, we can move on to creating a pipeline object. Finally, the function run_simulation() in the pipeline object runs the code and produces the outputs. The following script (say my_echo_script.py) helps you get started.
from echo21 import echopipeline
#Step-1 Set you parameter choices
cosmo = {'Ho':67.4,'Om_m':0.315,'Om_b':0.049,'sig8':0.811,'ns':0.965,'Tcmbo':2.725,'Yp':0.245}
astro = {'fLy':1,'sLy':2.64,'fX':1,'wX':1.5,'fesc':0.01}
#and choose your SFRD model type by defining a dictionary
sfrd = {'type':'phy','hmf':'press74','mdef':'fof','Tmin_vir':1e4}
#Step-2 Create an object and run
myobj = echopipeline.pipeline(cosmo=cosmo,astro=astro,sfrd=sfrd,path='/path/where/you/want/your/outputs/')
myobj.run_simulation()
#That's it.
Running the above script will generate an output folder in the path you gave in the path argument. Suppose you ran the script at 3:00:00 PM on 26th February 2025, then the output folder will have the name output_20250226-150000. To understand output structure, see Output structure below.
Please see the paper for an understanding of parameters fLy, sLy, etc. In brief, fLy, sLy, fX, wX, and fesc, are the Ly \(\alpha\) emissivity normalisation, power-law index of Ly \(\alpha\) SED, X-ray emissivity normalisation, power-law index of X-ray SED, and escape fraction of ionizing photons, respectively. Tmin_vir is the minimum virial temperature of the star-forming haloes.
Charting a parameter space of IGM: running ECHO21 in parallel mode
Now suppose you want to generate many 21-cm signals, neutral hydrogen fraction, and optical depths for different astrophysical (or cosmological) parameter. For this you simply have to provide your choice of parameters as a dictionary of lists or arrays.
Grid or no grid? there is an additional consideration. You can choose to generate the results for all possible combinations of the parameters (i.e., a grid) or only for the combinations of parameters at the same index in the lists/arrays. For example, if you have two parameters, each with three values, then in the grid case you will have 9 models corresponding to all possible combinations of the parameters, but in the no-grid case you will have only 3 models corresponding to the combinations of parameters at the same index. The choice is yours. By default, ECHO21 generates signals for all possible combinations of parameters (i.e., grid). To turn off this feature set grid_on=False when defining your pipeline object.
The following example shows you how. Replace astro = {'fLy':1,'sLy':2.64,'fX':1,'wX':1.5,'fesc':0.01} in my_echo_script.py by
astro = {'fLy':np.logspace(-2,2,5),'sLy':[-1,0,1],'fX':np.logspace(-2,2,5),'wX':[0,1,2],'fesc':[0.01,0.1,1]}
For minimum virial temperature you should modify the SFRD dictionary. So now instead of sfrd = {'type':'phy','hmf':'press74','mdef':'fof','Tmin_vir':1e4} you should have something like this
sfrd = {'type':'phy','hmf':'press74','mdef':'fof','Tmin_vir':np.logspace(2,6,5)}
Thus, the complete script to generate a large space of \(T_{21}\), \(x_{\mathrm{HI}}\), and \(\tau_{\mathrm{e}}\) with varying astrophysical parameters now looks like
import numpy as np
from echo21 import echopipeline
#Step-1 Set you parameter choices
cosmo = {'Ho':67.4,'Om_m':0.315,'Om_b':0.049,'sig8':0.811,'ns':0.965,'Tcmbo':2.725,'Yp':0.245}
astro = {'fLy':np.logspace(-2,2,5),'sLy':[-1,0,1],'fX':np.logspace(-2,2,5),'wX':[0,1,2],'fesc':[0.01,0.1,1]}
#and choose your SFRD model type by defining a dictionary
sfrd = {'type':'phy','hmf':'press74','mdef':'fof','Tmin_vir':np.logspace(2,6,5)}
#Step-2 Create an object and run
myobj = echopipeline.pipeline(cosmo=cosmo,astro=astro,sfrd=sfrd, grid_on=True, path='/path/where/you/want/your/outputs/')
myobj.run_simulation()
Now a total of \(5\times3\times5\times3\times3\times5=3375\) models will be generated corresponding to 5 values of \(f_{\mathrm{Ly}}\), 3 values of \(s_{\mathrm{Ly}}\), 5 values of \(f_{\mathrm{X}}\), 3 values of \(w_{\mathrm{X}}\), 3 values of \(f_{\mathrm{esc}}\), and 5 values of \(T_{\mathrm{vir}}\). (In the paper, I have used \(s\) for sLy and \(w\) for wX.)
Similarly, you can change the cosmo parameter in the above script to generate a large space of \(T_{21}\), \(x_{\mathrm{HI}}\), and \(\tau_{\mathrm{e}}\) with varying cosmological parameters. Further, ECHO21 is not limited to varying either astrophysical or cosmological parameters; both can be simultaneously varied.
In the above example we set grid_on to True. If you choose grid_on=False, then the parameters you wish to vary should have the same number of values. In this case, the code will generate models for combinations of parameters at the same index in the lists/arrays. As an exmaple
import numpy as np
from echo21 import echopipeline
#Step-1 Set you parameter choices
cosmo = {'Ho':67.4,'Om_m':0.315,'Om_b':0.049,'sig8':0.811,'ns':0.965,'Tcmbo':2.725,'Yp':0.245}
astro = {'fLy':np.logspace(-2,2,3),'sLy':[-1,0,1],'fX':np.logspace(-2,2,3),'wX':[0,1,2],'fesc':[0.01,0.1,1]}
#and choose your SFRD model type by defining a dictionary
sfrd = {'type':'phy','hmf':'press74','mdef':'fof','Tmin_vir':np.logspace(2,6,3)}
#Step-2 Create an object and run
myobj = echopipeline.pipeline(cosmo=cosmo,astro=astro,sfrd=sfrd, grid_on=False, path='/path/where/you/want/your/outputs/')
myobj.run_simulation()
The above script will generate 3 models only.
You can run the above script on your local PC as usual but with more than one CPU, as ECHO21 uses a master-worker CPU distribution. Thus, if you provide N CPUs, one CPU will act as the master CPU and remaining N-1 will act as worker CPUs. In general, generating a large number of models on a single CPU can be time consuming. To save time, you should utilize the parallel feature of ECHO21 and run the script my_echo_script.py as (say on four CPUs)
mpirun -np 4 python my_echo_script.py
Using a similar strategy you can now generate thousands of models in a few minutes with an appropriate choice of HPC resources.
Choosing a different HMF
Until now we have been using the Press & Schechter (1974) HMF. In ECHO21 you can choose a different HMF also. Suppose you want to generate a signal for Sheth & Tormen (1999) HMF. Then set 'sheth99' for the hmf keyword in the SFRD dictionary. For some HMFs you will have to change your definition of halo mass which is done by the keyword mdef. For example both Press & Schechter (1974) and Sheth & Tormen (1999) are based on the friends-of-friends definition (which is why we set 'fof' for mdef), but Tinker et al. (2008) is based on an integer multiple of mean matter density of the Universe. So you can give, say, '200m' for mdef. For a complete list of available HMFs see the COLOSSUS page.
Below is an example syntax for SFRD dictionary using Tinker et al. (2008) HMF.
sfrd = {'type':'phy','hmf':'tinker08','mdef':'200m','Tmin_vir':1e4}
Choosing a different SFRD model
Until now we have been working with physically-motivated star formation rate density (SFRD) models, which is why we had 'phy' for type in the SFRD dictionary. ECHO21 offers two additional models of SFRD – semi-empirical model and an empirically-motivated SFRD model. Let us first look at the semi-empirical model. The dictionary looks mostly the same as for the physically-motivated case, except now we use 'semi-emp' for type. Further, for this case now you also have an additional free parameter, t_star (default value 0.5). The dictionary now looks like
sfrd = {'type':'semi-emp','hmf':'press74','mdef':'fof','Tmin_vir':1e4, 't_star':0.5}
Let us now implement an empirically-motivated SFRD model. For this you need to set your SFRD type as 'emp' and choose the \(a\) parameter.
sfrd = {'type':'emp','a':0.257}
See section 2.2 from our paper for more details on SFRD.
Choosing the redshifts at which you want to evaluate the physical quantities
Before anything I want to clarify that I always work with \(1+z\) and NOT \(z\). So wherever, I write redshift I talk about \(1+z\). To avoid confusion I have used the capital letter zed (‘Z’) to represent \(1+z\).
Moving on to the main content of this section, when you do not specify the redshift range the code will evaluate the quantities at default redshifts. This default has 2300 values defined by the array Z_default given below.
import numpy as np
Z_cd = np.concatenate((1/np.linspace(1/60,1/5.05,200),1/np.linspace(1/5,1,100)))
Z_default = np.concatenate((np.linspace(1501,60,2001)[:-1],Z_cd))
When you run the code for a single set of parameters or vary cosmological parameters (irrespective of astrophysical ones) then the code will output the signal at redshits defined by Z_default by default. When you vary only astrophysical parameters then the code will output the signals at cosmic dawn redshifts defined by Z_cd.
How to give redshift values of your choice? Simple, just give your choice through the argument Z_eval when defining the pipeline object. For example, if you want to generate signal between \(1+z=30\) and \(1+z=10\) with 100 evenly spaced values then you should do the following
myZs = np.linspace(30,10,100)
myobj = echopipeline.pipeline(cosmo=cosmo,astro=astro,sfrd=sfrd,path='/path/where/you/want/your/outputs/',Z_eval=myZs)
Note: you don’t have to worry about giving redshifts in decreasing order. Whichever order you give, ECHO21 will always generate outputs for decreasing redshifts. When you are varying the astrophysical parameters only, the highest value of \(1+z\) should not be above 60.
Output structure
When you run ECHO21 for a single parameter the output folder will contain 9 files. These are redshifts (\(1+z\), not \(z\)), CMB temperature (Tcmb.npy), gas temperature (Tk.npy), spin temperature (Ts.npy), bulk IGM electron fraction (xe.npy), volume-filling factor (Q.npy), 21-cm signal (T21.npy), a text file glob_sig_<timestamp>.txt, and the class object echopipeline.pipeline as pipe.pkl. All .npy files are 1D arrays. They are evaluated at redshifts in the .npy file one_plus_z.npy. The CMB, gas, and spin temperatures are in units of kelvin and 21-cm signal is in units of milli kelvin (mK).
The text file contains all the basic information regarding your simulation such as the timestamp, execution time, cosmological & astrophysical parameters you provided. This file mentions the strongest 21-cm signal and the corresponding redshift, and the total CMB optical depth.
Simulation with a multi-valued parameters outputs differ slightly. The code will generate a summary file, pipeline object and an HDF5 file. The HDF5 file called, echo_output.h5 contains parameters as a pandas dataframe (params), and 1+redshift (Z), global signal (T21), neutral hydrogen fraction (xHI), and CMB optical depth (tau) as arrays. If \(N_{\mathrm{p}}, N_{z}\), and \(N_{\mathrm{mod}}\) are the number of parameters that are multi-valued, number of redshift points, and number of models, respectively then params is a table of shape \(N_{\mathrm{mod}} \times N_{\mathrm{p}}\), Z is 1D numpy array of length \(N_{z}\), T21 & xHI are 2D arrays of shape \(N_{\mathrm{mod}} \times N_{z}\), and tau is a 1D array of length \(N_{\mathrm{mod}}\).
Use the following lines of code to load the output:
import pandas as pd
with pd.HDFStore("echo_output.h5") as store:
params = store["params"]
one_plus_z = store["Z"].values
T21 = store["T21"].values
xHI = store["xHI"].values
tau = store["tau"].values
Thus, each row in params table corresponds to a unique model. Corresponding to this model (say at \(r^{\mathrm{th}}\) row) is the 21-cm signal at \(r^{\mathrm{th}}\) row, i.e., T21[r,:] (similarly for neutral hydrogen fraction and CMB optical depth). (Note that single-valued parameters do not appear in the table.) To see what values you gave, check the summary text file.