Getting Started

The best way to get started with all of the finer details of how to use the code and how it works is to review the notebook. Since the code is still in early development, all of the examples show the features of a voxel in the framework of this code. The three-dimensional uses of this code have thus-far been limited to models of the Milky Way as they appear in Yanitski (2023). Some of the results are shown in Three-dimensional PDR Models, but potential users are encouraged to contact yanitski@ph1.uni-koeln.de to iterate on their use-case.

The Single-Voxel model

Voxel representation

A diagram representing how the voxels work. The region in the red box is the ISM we want to model, with the clumps shown as scatter points. The blue boxes are the voxels. In order for the probabilistic calculation to be accurate, the clumps are randomly positioned in the voxel.

This is the base functionality of kosmatau3d. A voxel, or volumetric pixel, is one cell in a three-dimensional spatial grid. It contains one or more ensembles of KOSMA-\(\tau\) clumps (see Theory for more information about these), which each contribute to the emissivity and absorption in the voxel. It performs a probabilistic calculation of the emission by considering all possible combinations of clumps in the line-of-sight (see the schematic below).

Below is a minimal example of how to initialise and set the properties in a voxel. Note that there are plenty of other parameters that can be adjusted to modify the physics of the ensembles contained within the voxel, and most of them are explained in kosmatau3d.models.voxel.Voxel and in the jupyter notebook.

Initialise voxel object and set ISM properties

>>> from kosmatau3d import models
>>> vox = models.Voxel()
>>> vox.set_properties(
...     voxel_size=1,
...     clump_mass_number=[3],
...     clump_mass_range=[[0,2]],
...     ensemble_mass=1e1,
...     ensemble_density=1e5,
...     fuv=1e5,
... )

This will create a voxel with side length \(\ell_\mathrm{vox}=1\,\mathrm{pc}\) (so the volume is \(1\, \mathrm{pc}^3\)) containing an ensemble with three distinct clump masses: \(1\), \(10\), and \(100\, M_\odot\). The number of each of these clumps in the ensemble is calculated from the clump-mass distribution as explained in Theory. Since the mass of the ensemble in this example is lower than the maximum mass of the clumps in the ensemble, this example will only contain a fraction of the largest clump. This showcases a key property of voxels in kosmatau3d: the scale-invariance of the calculations.

Now that a voxel has been initialised with an ensemble, it is possible to obtain intrinsic properties such as the fractional abundance by,

>>> vox.get_abundances(species=["H2", "H", "C+", "CO"])

where the abundances are normalised by the total hydrogen abundance (\(n_\mathrm{H} = n_\mathrm{H^0} + 2\, n_\mathrm{H_2}\)). Extrinsic properties such as column density are also available at this point, but for the emission it is necessary to first execute the computation:

>>> vox.calculate_emission()

This will compute the emissivity and absorption for the line and continuum emission from the base KOSMA-\(\tau\) models. One can then access the result using,

>>> vox.get_species_emissivity()

Read the documentation of kosmatau3d.models.voxel.Voxel to find out the different arguments for this method. In the current procedure for calculating the emission, we use the peak intensity and optical depth as well as post-processing the HI \(21\, \mathrm{cm}\) line emission by making an isothermal approximation. for that reason we have split the accessing methods for the continuum and line emission into get_dust_... and get_species_..., respectively, and there is a kwarg hi that can be set to True to get the HI line emission. The emission values that are available are,

• emissivity \(\epsilon_\nu\) in \(\frac{K}{pc}\)

• absorption \(\kappa_\nu\) in \(\frac{1}{pc}\)

• intensity \(I_\nu\) in \(K\)

• optical depth \(\tau_\nu\) (dimensionless)

The intensity and optical depth require integrating over the length-scale of the voxel, so they should not be used in three-dimensional models.

Three-dimensional PDR Models

Galactic model

One Galactic model, where the voxels are coloured according to the [CII] \(158\, \mu\mathrm{m}\) intensity integrated over the spectrum (\(-350\) to \(350\, \mathrm{km\, s^{-1}}\)).

The most complex functionality of kosmatau3d, and the reason for its development since its conception (Andree-Labsch et al. 2017). This has so-far been used to model the Orion Bar and the Milky Way, though more models will soon be developed. The figure below depicts the [CII] \(158\, \mu\mathrm{m}\) integrated intensity in each voxel of one of the galactic models, which is then used to compute the synthetic observation.

The benefit of using kosmatau3d voxels for the model is two-fold: it uniquely accounts for the inhomogeneity and shadowing in PDRs and has a velocity component to its emission. Because of this, it is important to first get the velocity information for each voxel in the model correct before integrating the radiative transfer equation. The procedure is described in detail in Ch. IV of Yanitski (2023), and it results in a longitude-velocity diagram like below. We focus on galactic latitude \(b\! =\! 0\) to avoid the complications regarding partially-filled voxels.

We use the kosmatau3d.models.model.Model() to initialise all voxels, compute their radiative properties, and save the relevant data in FITS files. It shares many kwargs with set_properties(), though some are renamed to fit the context of a three-dimensional model. A minimal working example to create a galactic model with voxel size \(\ell_\mathrm{vox}=400\,\mathrm{pc}\) is,

Initialise galactic model and compute emission

>>> from kosmatau3d import models
>>> kwargs = { ... }
>>> galaxy = models.Model(resolution=400,
...     history_path='.',
...     folder='temp',
...     **kwargs)
>>> galaxy.calculate_model()

[CII] \(158\, \mu\mathrm{m}\)

The synthetic emission resulting from the model above. Note the large-scale velocity structure of the Milky Way is replicated.

Here kwargs can be used to specify any of the model parameters. A distinct difference in making the full model is that the kwargs are given when initialising the object instance rather than through a separate method. While all of the physical and emissive properties are calculated at this stage, a synthetic observation requires the kosmatau.radiative_transfer:

Compute synthetic observation

>>> import numpy as np
>>> from models.radiative_transfer import calculateObservation
>>> calculateObservation(directory='temp/',
...     slRange=[(-np.pi, np.pi), (-np.pi/2, np.pi/2)],
...     nsl=[180, 90])

This will result in a synthetic datacube of the region for all of the included transitions (by default all of them) and a subset of the dust continuum (where 22 wavelengths are used; enough to span the FIR emission). From the synthetic intensity datacube, it is possible to get the position-velocity diagram as below.

It should be noted, though, that the procedure described thus-far in this section is for one model, but for scientific modelling it is likely useful to analyse the sythetic emission from a grid of models to constrain some parameters. There is a convenient method to do this with the github repository. From the root directory of the repo, we can run a grid of models using,

$ mkdir ../kt3_models
$ python terminal_scripts/run_model_grid.py -f ../kt3_models -m 0

By default, this will run a grid of three models of varying resolution (specifically \(400\,\mathrm{pc}\), \(200\,\mathrm{pc}\), and \(100\,\mathrm{pc}\)), though you may notice that it takes a long time to finish. For that reason, it might be better to set -m 8 for example to multiprocess the radiative transfer calculation. At the moment, calculateModel() does not have the ability to utilise multiprocessing.

References

Andree-Labsch, S., Ossenkopf-Okada, V., & Röllig, M. 2017, Astronomy & Astrophysics, 598, A2

Yanitski 2023, The Milky Way with kosmatau3d: Modelling the Galactic cooling lines using clumpy PDRs, PhD thesis, (Universität zu Köln)