◀ Back

Cartesian outputs

While running the THESAN simulations we output high time cadence Cartesian outputs of the full particle data onto a regular Cartesian grid using a first order Particle-In-Cell approach. THESAN-1 uses a 10243 grid while the medium resolution runs use 5123 grids, setting the cell sizes to ∼93 and ∼186 ckpc, respectively. For THESAN-1, we also provide a coarse-grained version of the Cartesian output on a 5123 grid, to enable a faithful comparison with the medium-resolution runs. While most of the quantities like the gas and stellar density, ionization fraction, and temperature are binned on a regular spatial grid, observables like the 21 cm emission and Lyman-alpha luminosities are binned in redshift space which accounts for Doppler shifting due to the peculiar velocities of the galaxies (i.e. individual gas cells). This redshift-space binning is carried out assuming the observed direction is aligned with the +z axis. The Cartesian outputs are written every ∼2.8 Myr amounting to about 400 snapshots over the duration of the simulation. These lower resolution higher time-cadence outputs are ideal for studying the large scale properties of the reionization process.

We now provide the full list of Cartesian output quantities. The main simulation attributes are provided in the Header HDF5 group:
Header attributes ▾
Attribute Dimensions Units Description
BoxSize 1 \({\rm ckpc} / h\) Spatial extent of the periodic box (in comoving units).
HubbleParam 1 \(100\,{\rm km/s/Mpc}\) Hubble constant (little \(h\) in standard units).
NumFiles 1 - Number of subfiles per output.
NumPixels 1 - Number of pixels (cells) in each Cartesian dimension.
Omega0 1 - The cosmological density parameter for matter.
OmegaBaryon 1 - The cosmological density parameter for baryonic matter.
OmegaLambda 1 - The cosmological density parameter for the cosmological constant.
Redshift 1 - The redshift corresponding to the current output.
Time 1 - The scale factor \(a=1/(1+z)\) corresponding to the current output.
UnitLength_in_cm 1 - Code length units in cm.
UnitMass_in_g 1 - Code mass units in g.
UnitVelocity_in_cm_per_s 1 - Code velocity units in cm/s.
The following HDF5 datasets are on 3D spatial grids:
Datasets on spatial grid ▾
Attribute Dimensions Units Description
Density N \((10^{10}\,{\rm M}_{\odot} / h) / ({\rm ckpc} / h)^3\) Comoving gas mass density (calculated as gas mass/volume).
DensityStars N \((10^{10}\,{\rm M}_{\odot} / h) / ({\rm ckpc} / h)^3\) Comoving stellar mass density (calculated as stellar mass/volume).
DensityDust N \((10^{10}\,{\rm M}_{\odot} / h) / ({\rm ckpc} / h)^3\) Comoving dust mass density (calculated as dust mass/volume).
DensityMetals N \((10^{10}\,{\rm M}_{\odot} / h) / ({\rm ckpc} / h)^3\) Comoving metal mass density (calculated as metal mass/volume).
HII_Fraction N - Ionized hydrogen fraction \(n_{\rm HI} / n_{\rm H}\) (mass-weighted).
HeIII_Fraction N - Doubly ionized helium fraction \(n_{\rm HeIII} / n_{\rm H}\) (mass-weighted).
Temperature N \({\rm K}\) Gas temperature (mass-weighted).
IonEnergy N \((10^{10}\,{\rm M}_{\odot} / h) \cdot ({\rm km/s})^2\) Ionizing radiation energy, \(E_{\rm ion}\).
IonFlux N,3 \((10^{10}\,{\rm M}_{\odot} / h) \cdot ({\rm km/s})^3\) Ionizing radiation flux, \({\bf F} = (F_x,F_y,F_z)\).

Note: The IonEnergy and IonFlux fields are not densities. They are also affected by the reduced speed of light \(\tilde{c} = 0.2\,c\). For example, after converting to cgs units the photo-ionization integral would be \(\Gamma_{\rm HI} = \tilde{c} \sigma_{\rm HI} {\tt IonEnergy} / {\tt MeanPhotonEnergy} / {\tt VoxelVolume}\).

The following HDF5 datasets are on redshift-distorted spatial grids:
Datasets on spectral grid ▾
Attribute Dimensions Units Description
DensityHI N \((10^{10}\,{\rm M}_{\odot} / h) / ({\rm ckpc} / h)^3\) Comoving HI mass density (calculated as \(0.76 \times x_{\rm HI} \times\) gas mass/volume) for the 21 cm emission.
IonLuminosityStars N \((10^{10}\,{\rm M}_{\odot} / {\rm kpc}) \cdot ({\rm km/s})^3\) Ionizing luminosity emitted by stars.
IonLuminosityAGN N \((10^{10}\,{\rm M}_{\odot} / {\rm kpc}) \cdot ({\rm km/s})^3\) Ionizing luminosity emitted by AGN.
LyaLuminosityRec N \((10^{10}\,{\rm M}_{\odot} / {\rm kpc}) \cdot ({\rm km/s})^3\) Lyman-alpha luminosity from resolved recombination emission.
LyaLuminosityCol N \((10^{10}\,{\rm M}_{\odot} / {\rm kpc}) \cdot ({\rm km/s})^3\) Lyman-alpha luminosity from collisional excitation emission.
StarFormationRate N \({\rm M}_{\odot} / {\rm yr}\) Instantaneous star formation rate of this voxel.

Note: The IonLuminosity, LyaLuminosity, and StarFormationRate fields are not densities.

Similar to snapshots and group catalogs the Cartesian data are split into multiple files named as: output/cartesian_NNN/cartesian_NNN.CCC.hdf5, where NNN is the left-zero-padded Cartesian output number, and CCC is the left-zero-padded chunk file number. The three-dimensional grid is flattened into a C-ordered array (with the z-axis increasing fastest) and then split into chunks each saved in a different HDF5 file. For convenience, we provide example python code to read the Density field:
Example code ▾
import numpy as np
import h5py

def read_render(snap=300, ren_dir='Thesan-1/output'):
    # Read header info
    with h5py.File(f'{ren_dir}/cartesian_{snap:03d}/cartesian_{snap:03d}.000.hdf5', 'r') as f:
        g = f['Header'] # Access the HDF5 Header group
        NumFiles = g.attrs['NumFiles'] # Number of subfiles per output
        NumPixels = g.attrs['NumPixels'] # Number of cells in each Cartesian dimension
        NumPixels3 = NumPixels**3 # Total number of Cartesian cells
        BoxSize = g.attrs['BoxSize'] # Spatial extent of the periodic box [ckpc/h]
        HubbleParam = g.attrs['HubbleParam'] # Hubble constant [100 km/s/Mpc]
        Time = g.attrs['Time'] # Scale factor a=1/(1+z) of the current output
        Redshift = g.attrs['Redshift'] # Redshift of the current output
        Omega0 = g.attrs['Omega0'] # Cosmological matter density [rho_crit_0]
        OmegaBaryon = g.attrs['OmegaBaryon'] # Cosmological baryonic matter density [rho_crit_0]
        OmegaLambda = g.attrs['OmegaLambda'] # Cosmological constant density [rho_crit_0]
        UnitLength_in_cm = g.attrs['UnitLength_in_cm'] # Code length [cm]
        UnitMass_in_g = g.attrs['UnitMass_in_g'] # Code mass [g]
        UnitVelocity_in_cm_per_s = g.attrs['UnitVelocity_in_cm_per_s'] # Code velocity [cm/s]

    # Define various convenience factors
    a,z,h = Time, Redshift, HubbleParam
    length_to_cgs = a * UnitLength_in_cm / h
    volume_to_cgs = length_to_cgs**3
    mass_to_cgs = UnitMass_in_g / h
    density_to_cgs = mass_to_cgs / volume_to_cgs
    velocity_to_cgs = np.sqrt(a) * UnitVelocity_in_cm_per_s
    cell_volume = (length_to_cgs * BoxSize / float(NumPixels))**3
    print(f'snap = {snap}  |  BoxSize = {BoxSize} ckpc/h  |  NumPixels = {NumPixels}')

    # Read full data
    offset = 0
    Density = np.zeros(NumPixels3, dtype=np.float32)
    for sub in range(NumFiles):
        filename = f'{ren_dir}/cartesian_{snap:03d}/cartesian_{snap:03d}.{sub:03d}.hdf5'
        with h5py.File(filename, 'r') as f:
            n_buf = len(f['Density'])
            next_offset = offset + n_buf
            Density[offset:next_offset] = f['Density'][:]
            offset = next_offset # Keep track of file offset
    Density = Density.reshape([NumPixels,NumPixels,NumPixels]) # 3D density grid