Importing snapshots

# Introduction

Analytical or quasi-analytical models often fall short in describing the geometric complexity of the real world. More realistic numerical models can be produced by computer programs that simulate the (magneto-)hydrodynamical evolution of an astrophysical system over time. Such codes usually output a series of snapshots. The data provided in each snapshot define the state of the simulated system at a particular time in its evolution, including the spatial distribution and other relevant properties of the dynamical components, such as for example the stars and the interstellar medium in a galaxy.

Because of the nonlocal and nonlinear behavior of radiation transport, predicting the observable properties of these hydrodynamical models requires a full, three-dimensional radiative transfer simulation in all but the most trivial cases. Assuming that the radiation crossing times are much shorter than the dynamical time scales, and that radiation pressure is negligible, the radiative transfer and hydrodynamical simulations can be completely uncoupled. In other words, the radiative transfer simulation can then be performed as a post-processing step based solely on the information contained in a particular snapshot. The post-processing procedure extracts the spatial distribution of the radiation sources (e.g., stars) and the obscuring medium (e.g., dust), plus the corresponding properties, from the hydrodynamical snapshot data, and then performs a radiative transfer simulation on the resulting setup.

The radiation sources and the obscuring medium are obviously handled differently during a SKIRT simulation, implying different needs for the corresponding data structures. The most frequent demand on the data structure representing the radiation sources is to randomly sample a photon packet launch position (from the spatial luminosity distribution) and wavelength (from the spectral energy distribution at that position). This operation usually occurs millions or even billions of times during a simulation, and thus should execute fairly quickly. On the other hand, the key operation for the obscuring medium is to determine its density at a given location, or its mass in a given volume. These functions are needed while constructing the spatial grid, i.e. for deciding where to subdivide cells in adaptive grids, and for calculating the medium mass in each cell of the final spatial grid. Although this is an intensive process, it is limited to the setup phase of the simulation, because, once constructed, the spatial grid includes all information required for tracing photon packets through its domain.

Hydrodynamical simulation codes historically employ one of two schemes: a Lagrangian formulation based on moving particles (smoothed particle hydrodynamics or SPH), for example Gadget and Gasoline; or a Eulerian approach based on a non-moving spatial grid, often an adaptive grid with multiple refinement levels (AMR), for example RAMSES, Enzo, and AMR-VAC. Recent codes including AREPO, TESS, and Shadowfax introduce a new scheme that employs a moving mesh based on a Voronoi tessellation of the spatial domain. This new scheme is claimed to combine the best features of SPH and the traditional Eulerian approach, and it is becoming increasingly popular.

In this note we discuss SKIRT's mechanisms for importing radiation sources and obscuring media from snapshots produced by each of these three schemes. Note that each source component and each medium component in SKIRT has its own distinct import mechanism. It is perfectly possible to import, say, radiation sources from particles and medium densities from a mesh-based grid. Such a combination is not uncommon. Several codes, for example, use an N-body solver for the pressure-less dark matter and stellar particles while employing a mesh-based hydrodynamical solver for the gas medium in the same simulation.

The discussion is organized in sections as follows:

# Snapshots

## Snapshot file formats

All snapshot data must be provided to SKIRT for import in a predefined column text format. In other words, an external procedure (usually a Python script) is needed to extract the appropriate information from the native snapshot data format and output the result as an ASCII column text file. While this seems harsh, it is almost never an issue in practice. In most cases, one needs to perform some additional processing anyway, such as locating the proper data files, performing a coordinate transformation, limiting the imported data to a given aperture, calculating indirectly derived properties, etc. Furthermore, the alternative would be to implement a variety of data formats and point SKIRT to the relevant data through a potentially very complex configuration mechanism that would never be as flexible as the Python language.

There are three supported snapshot paradigmns (particle, adaptive mesh, and Voronoi mesh), each with their own format for describing the spatial geometry. In addition to describing the spatial geometry, the input file must contain the properties required by the component (e.g. source, medium or geometry) being imported. The number and type of such properties depends on the component's type and on its configured options. This results in the following formats:

• Particle (SPH): the input file has a line for each particle, listing its position and smoothing length plus the component-specific properties.
• Adaptive mesh (AMR): the input file has lines defining the nonleaf/leaf structure of the mesh, plus a line for each leaf cell with its component-specific properties.
• Voronoi mesh: the input file has a line for each generating site, listing its position plus the component-specific properties.

In all cases, columns must be in the fixed order prescribed by the SKIRT component controlling import of the file (see Snapshot classes, Sources, and Media). Values are either in the default units prescribed by that SKIRT component, or in the units defined through header lines in the data file being imported (see the documentation of the TextInFile class). In the latter case, any of the ski file units corresponding to the imported quantity can be used (see Supported units).

## Snapshot classes

The heart of the SKIRT snapshot import mechanism is implemented by the classes in the following hierarchy:

Snapshot is an abstract base class for representing snapshot data imported from a column text file. It can be used to define the properties of primary sources or transfer media. This base class offers facilities to handle the input options, such as determining which columns are to be imported from the text file, and defines a common interface for all snapshot types, for example to obtain the density at a given position. The actual implementation is provided in a derived class for each snapshot type (ParticleSnapshot, AdaptiveMeshSnapshot, and VoronoiMeshSnapshot).

The Snapshot classes are not exposed to the user through the configuration file. Instead, they are invoked under the hood by the various components that are configured as part of the simulation. Three distinct class hierarchies actually make use of the Snapshot classes:

Source                  Medium                  Geometry
ImportedSource          ImportedMedium          ImportedGeometry
ParticleSource          ParticleMedium          ParticleGeometry       --> ParticleSnapshot
VoronoiMeshSource       VoronoiMeshMedium       VoronoiMeshGeometry    --> VoronoiMeshSnapshot


Each of these class hierarchies implements the appropriate layer of functionality on top of the Snapshot classes. For example, sources can launch photo packets, media can quickly find the material density at a given location, and geometries offer a density distribution normalized to unity.

# Sources

A SKIRT simulation always includes one or more components describing the radiation sources in the model (see the SourceSystem and Source classes). Each source component defines (a) the spatial density distribution of the sources and (b) the spectral energy distribution (SED) at each location. In the context of post-processing hydrodynamical simulations, there are three distinct ways of accomplishing this:

• The imported snapshot defines the spatial density distribution and provides additional properties for each smoothed particle or grid cell from which an appropriate SED can be determined. For this purpose, SKIRT includes several parameterized SED families. For example, the Bruzual-Charlot family (the BruzualCharlotSEDFamily class) represents young and evolved stellar populations. The Mappings III SED family (the MappingsSEDFamily class) represents star-forming regions. New families can be added with minimal programming effort, provided the data for the SED templates is available in tabulated form.
• The imported snapshot defines just the spatial density distribution, and a particular SED (constant throughout the spatial domain) is configured from the SED components built into SKIRT (subclasses of the SED class).
• The radiation sources are completely defined within SKIRT by configuring a built-in geometry (subclasses of the Geometry class) and SED (subclasses of the SED class). In case a hydrodynamical simulation does not include any sources, this option can be used to complete the model by adding, for example, a central source to a dust disk, or background radiation to the dust in a molecular cloud.

One can even combine these options by using multiple source components, for example, for adding an active galactic nucleus to a spiral galaxy snapshot.

## Particle sources

A particle snapshot consists of a set of smoothed particles (or rather anchor points in a co-moving grid), each characterized by a suite of properties. Formally, the spatial luminosity distribution in a particle snapshot is defined by a list of $$N$$ particles and an assumed smoothing kernel $$W(r,h)$$, with each smoothed particle characterized by a position $$\mathbf{r}_j$$, a smoothing length $$h_j$$, and a luminosity contribution $$L_j$$. The total luminosity density at an arbitrary position $$\mathbf{r}$$ is then given by

$\Sigma(\mathbf{r}) = \sum_{j=1}^N L_j \, W(|\mathbf{r}-\mathbf{r}_j|,h_j)$

The ParticleSource class allows configuring one of several (spherically symmetric) smoothing kernels, including the frequently-used cubic spline kernel (the CubicSplineSmoothingKernel class) and a scaled Gaussian smoothing kernel (the ScaledGaussianSmoothingKernel class). If needed, implementing a new smoothing kernel should be fairly straightforward.

For launching photon packets, SKIRT needs to generate random positions sampled from the radiation source's density distribution. This is rather straightforward thanks to the composition method. The first step is the choice of a random particle, based on a discrete distribution where each particle is weighted by its relative luminosity contribution. The second step is generating a random position according to the luminosity distribution of the chosen particle, using the configured smoothing kernel.

### Spatial distribution and SED from snapshot

We first consider the case where both the spatial distribution and the SED at each location are determined from the SPH snapshot. The ParticleSource class expects a text file with a single line for each particle. The first four columns specify the $$x$$, $$y$$ and $$z$$ coordinates and the smoothing length $$h$$ for the particle. The number and interpretation of the subsequent columns depends on the SED family configured for this ParticleSource instance. For example, for the Bruzual-Charlot SED family (BruzualCharlotSEDFamily), the remaining three columns provide the properties of the stellar population represented by the particle. The first extra column specifies the initial mass of the stellar population, i.e. the particle's mass $$M_\mathrm{init}$$ at the time it was born. The second extra column specifies the metallicity $$Z$$ of the stellar population, and the third extra column represents the age $$t$$ of the stellar population.

The ParticleSource class can also be configured to take into account the Doppler shift caused by the velocity of the radiation sources relative to the model coordinate system. If this option is enabled, the input text file must provide three additional columns (after the smoothing length and before the additional properties for the SED family), specifying the $$v_\mathrm{x}$$, $$v_\mathrm{y}$$ and $$v_\mathrm{z}$$ components of the particle's velocity. Given this velocity $$\mathbf{v}$$, the redshift $$z$$ experienced by a photon packet launched from the particle in direction $$\mathbf{k}$$ (given as a unit vector) is determined by

$z = -\frac{\mathbf{v}\cdot\mathbf{k}}{c}$

where $$c$$ is the speed of light in vacuum. As a function of the rest wavelength $$\lambda_0$$, the redshifted wavelength $$\lambda_z$$ is given by $$\lambda_z=(1+z)\,\lambda_0$$. Since the redshift value depends on the angle between the particle and photon packet trajectories, the emission is no longer isotropic.

The default units for these quantities can be overridden by including columns header info in the file as described in the TextInFile class header. For example, here is an input file describing just a single source particle placed at the origin of the model coordinate system:

# Star particle for testing SED family
# Column 1: position x (pc)
# Column 2: position y (pc)
# Column 3: position z (pc)
# Column 4: smoothing length (pc)
# Column 5: mass (Msun)
# Column 6: metallicity (1)
# Column 7: age (Gyr)
0 0 0 1   10  0.04  5


### Spatial distribution from snapshot with built-in SED

If the hydrodynamical simulation does not track the relevant properties for defining the local emission spectrum, we can still import the spatial luminosity distribution from the snapshot and assign a built-in SED that is constant across the spatial domain. In this case, we use the GeometricSource class, and we configure it with an instance of the ParticleGeometry class to import the snapshot. Similar to the ParticleSource class described earlier, the ParticleGeometry class expects a text file with a single line for each particle. The first four columns specify the $$x$$, $$y$$ and $$z$$ coordinates and the smoothing length $$h$$ for the particle. The fifth column specifies the luminosity contribution of this particle, in arbitrary units since the ParticleGeometry class normalizes the total luminosity to unity. Because of this normalization, we need to separately supply the actual luminosity of the source through one of the regular options offered by the GeometricSource class.

Simulations with an a adaptive mesh paradigm use a hierarchical grid that is refined to a variable level depending on the resolution requirements in various regions of the model. While the resolution requirements and thus the grid may change as the system evolves, each particular snapshot corresponds to a unique and well-defined hierarchical grid. Thus, an adaptive mesh snapshot includes some description of the structure of the hierarchical grid, implicitly or explicitly defining the spatial extent of each grid cell, and for each cell the values of a suite of physical properties.

While it is possible to use this technique with spherical or cylindrical coordinate systems, perhaps to benefit from certain quasi-symmetries in the model, SKIRT only supports adaptive mesh snapshots using Cartesian coordinates. In a Cartesian hierarchical grid, the cuboidal spatial domain is recursively subdivided into smaller cuboids according to some scheme until the desired resolution in each region is reached.

### Snapshot data format

The AdaptiveMeshSource class expects a text file in a format that describes the topological structure of the Cartesian hierarchical grid and lists the relevant physical properties in each cell. The size of the spatial domain must be defined separately in the SKIRT configuration file, i.e. this information is not part of the imported data.

The hierarchical grid structure is organized into a tree. Each tree node represents a cuboidal portion of the domain, called its extent. The root node's extent is the complete domain. A nonleaf node distributes its extent over its child nodes using a regular linear grid. The number of subdivisions is defined separately for each node and can differ for each spatial dimension. An octree, for example, would subdivide each nonleaf node into $$2\times2\times2$$ child nodes. Typical AMR schemes have much larger subdivision counts that sometimes vary within the grid; for example, the root node may be subdivided differently. A leaf node represents a cell that is not subdivided any further and that holds a data value for each field; the field values are considered to be constant over the leaf node's extent. Collectively the leaf nodes form a partition of the domain, i.e. their extents cover the complete domain without overlapping one another.

The leaf nodes or cells in this three-dimensional data structure can be arranged in a linear sequence using so-called Morton ordering. This ordering is obtained by performing a depth-first traversal of the tree, where each nonleaf node visits its children in the order x-first, then y, then z. The process is illustrated in the following figure for a two-dimensional structure. The root node is subdivided into $$4\times3$$ subnodes; some of the nodes near the center are subdivided into $$2\times2$$ subnodes, and two of the cells in the upper row are subdivided into $$1\times2$$ subnodes. The solid line connects the leaf nodes in Morton order, starting at the lower left node with Morton index 0 and performing a depth-first traversal with each nonleaf node visiting its children in the order first x (horizontal) and then y (vertical).

The overall makeup of the snapshot text file reflects this structure and ordering. Below is a text representation of the 2D grid shown in the above figure. Because SKIRT deals with 3D grids, we've added a third dimension that is never subdivided. The lines starting with an exclamation mark indicate a subdivided node, e.g., the root node is subdivided in $$4\times3\times1$$ cells. The other lines specify the field values for a particular cell. In this example, the first field value specifies the serial number or Morton index of the cell, and the second value specifies the cell's "density" reflected by the gray level in the above figure. In an actual snapshot, these illustrative values would be replaced by the relevant physical quantities.

# Snapshot text data file
! 4 3 1
0 0.4
1 0.4
2 0.4
3 0.4
4 0.4
! 2 2 1
5 0.6
6 0.6
7 0.6
! 2 2 1
8 0.8
9 0.8
10 0.8
11 0.8
! 2 2 1
12 0.6
13 0.6
14 0.6
15 0.6
16 0.4
17 0.4
! 1 2 1
18 0.5
19 0.5
! 1 2 1
20 0.5
21 0.5
22 0.4


More formally, each line in the text file describes a particular tree node (nonleaf or leaf), and the lines are given in Morton order. Specifically, each line in the file can be of one of the following types:

• Comment: lines having a number sign or hash (#) as the first non-whitespace character, lines containing only whitespace and empty lines are ignored (and do not count in the Morton order).
• Nonleaf: a nonleaf line has an exclamation mark (!) as the first non-whitespace character, followed by optional whitespace and then three whitespace-separated positive integer numbers $$N_x,N_y,N_z$$. These three numbers specify the number of child nodes carried by this node in each spatial direction. The child nodes are on a regular linear grid as described above.
• Leaf: a leaf node contains one or more whitespace-separated floating point numbers reflecting the physical quantities associated with the leaf node, depending on the user configuration settings for the snapshot. The default units for these quantities can be overridden by including columns header info in the file as described in the TextInFile class header.

Note that there is no need to include the cell positions in the lines for the leaf nodes because the geometry follows from the Morton order.

### Spatial distribution and SED from snapshot

As before, we first consider the case where both the spatial distribution and the SED at each location are determined from the snapshot. The AdaptiveMeshSource class expects a snapshot data file in the format described in the previous section, supplemented with the extent of the spatial domain in each direction (specified separately as part of the configuration), because that information is not stored in the file.

As with other imported sources (see Spatial distribution and SED from snapshot), the AdaptiveMeshSource class allows configuring an SED family from which a specific template is selected for each leaf node. To this end, each leaf-node line in the import file specifies the parameters for the configured SED family, possibly preceded by the bulk velocity components.

Given this information, SKIRT calculates the luminosity contribution of each leaf node or cell. This is used to determine the probability that a given cell will be selected as the launch site for a new photon packet. Generating a random launch position from a uniform distribution within the selected cuboidal cell is of course trivial.

### Spatial distribution from snapshot with built-in SED

Similar to the procedure described for SPH snapshots in Spatial distribution from snapshot with built-in SED, we can import the spatial luminosity distribution from a mesh-based snapshot and assign a built-in SED that is constant across the spatial domain. In this case, we configure a GeometricSource with an AdaptiveMeshGeometry to import the snapshot, separately specifying the extent of the spatial domain as before. The values in the import file now specify the luminosity contribution for each cell.

The AdaptiveMeshGeometry class offers various options for compatibility with medium density distributions (see Loading snapshots). In the context of source luminosities, the massType option is the only relevant one:

• for an "integrated mass" type, the value in the import file indicates the luminosity contribution of the cell;
• for an "average density" type, the value in the import file is multiplied by the volume of the cell to determine the cell's luminosity contribution.

In both cases the units of the specified values are irrelevant because the AdaptiveMeshGeometry class normalizes the total luminosity to unity. Because of this normalization, we need to separately supply the actual luminosity of the source through one of the regular options offered by the GeometricSource class.

## Voronoi mesh sources

Due to the nature of a Voronoi tessellation, the geometry of the grid is completely defined by the positions of the generating sites. It is hence not necessary for a snapshot to store all the vertices and edges of each of the cells. Similar to our approach for other snapshot types, SKIRT reads the properties for a list of Voronoi cells from a text file in a simple format. SKIRT constructs the Voronoi grid from the positions of the generating sites, and assigns the corresponding physical values to each cell, assuming that the values are constant across the cell's extent.

### Launching photon packets

The generation of random positions for launching photon packets is similar to the procedure for hierarchical grids. In the first step, we pick a random cell using a discrete distribution where each cell is weighted by its relative luminosity contribution. The second step, generating a random position from within the chosen cell, is significantly more complex than in the case of a cuboidal cell. To the best of our knowledge, there are no dedicated techniques to generate a random point from a Voronoi cell. There are two possible options.

The first option is to partition the cell into a set of tetrahedra, subsequently select a random tetrahedron from a discrete distribution where every tetrahedron is weighted by its relative volume, and finally generate a random position from the selected tetrahedron. Specific algorithms are available for both the tetrahedrization of convex polyhedra and the generation of random positions from a tetrahedron.

The second option, which is more straightforward and which we adopted in SKIRT, is to use the rejection technique. As the reference distribution we use a uniform density in a cuboidal volume, defined as the 3D bounding box of the cell. During construction of the Voronoi grid, for each cell, we store the cell bounding box (easily obtained from the cell vertices because Voronoi cells are convex polyhedra) and a list of all neighboring cells. After generating a candidate random position uniformly from the cell bounding box, we decide whether the position is actually inside the Voronoi cell by checking that it is closer to the cell's generating site than to the generating sites of all its neighbors.

Our tests have shown that, depending on the distribution of the generating Voronoi sites, the average ratio of the volume of the bounding box of a Voronoi cell over the actual cell volume is about 3 to 4. This ratio represents the average rejection rate for the random position generation.

The procedures for loading a spatial and/or spectral distribution from a Voronoi-tessellation-based snapshot into SKIRT are very similar to what was described for the other snapshot types in Particle sources and Adaptive mesh sources.

The VoronoiMeshSource class expects a snapshot data file in a straightforward text column format, containing one line per cell. The first three columns provide the $$x$$, $$y$$, and $$z$$ coordinates of the generating site for the cell. The subsequent columns specify the relevant physical quantities for determining the SED assigned to the cell depending on the configured SED family.

The VoronoiMeshGeometry class, when configured as part of a GeometricSource, imports just the spatial luminosity distribution from a Voronoi-based snapshot, allowing to separately assign a built-in SED that is constant across the spatial domain. The file format is similar as described above. The luminosity contribution of each Voronoi cell is now given by the corresponding value in the import file, possibly multiplied by the volume of the cell as described in Spatial distribution from snapshot with built-in SED. Again, because the VoronoiMeshGeometry class normalizes the total luminosity to unity, we need to separately supply the actual luminosity of the source through one of the standard normalization options.

# Media

Apart from radiation sources, any nontrivial SKIRT simulation also includes a medium system with one or more medium components describing the transfer medium in the model (see the MediumSystem and Medium classes). Each medium component defines (a) the spatial density distribution of the medium and (b) the relevant material properties of the medium at each location. In addition, the medium system also configures a spatial grid that will be used to discretize the spatial domain for the purposes of the radiative transfer simulation (see the SpatialGrid class).

In the context of post-processing hydrodynamical simulations, there are several ways of accomplishing this. Let us first examine the most common situation where a medium component is assigned material properties that do not vary across the spatial domain (although different medium components may be assigned different material properties). In this case, SKIRT offers three distinct ways to define a component's medium density distribution:

• An imported snapshot fully defines the medium density distribution including the total mass, either directly or using a simple formula, for example to derive a dust density from a gas density. With this option, the medium component is configured using one of the ImportedMedium subclasses (ParticleMedium, AdaptiveMeshMedium, VoronoiMeshMedium). These classes can import velocity information for each particle or cell in the imported snapshot.
• An imported snapshot defines the medium density distribution normalized to unity, either directly or using a simple formula, and the total mass is specified separately in the configuration file. With this option, the medium component is configured using the GeometricMedium class and one of the ImportedGeometry subclasses (ParticleGeometry, AdaptiveMeshGeometry, VoronoiMeshGeometry). These classes cannot import spatially varying velocity information.
• The medium density distribution is defined by configuring a built-in geometry (subclasses of the Geometry class). This option can be useful to combine medium distributions imported from a snapshot with built-in, semi-analytical distributions.

In any case, each medium component must be assigned a specific material mixture, i.e. one of the subclasses of the MaterialMix class, which defines the optical and calorimetric properties of the medium, considered to be constant across the spatial domain.

Let us now examine the situation where it is desired to associate multiple sets of material properties with a single medium component depending on the location of the medium. This feature is supported only when configuring one of the ImportedMedium subclasses (i.e. the first option in the list above). Specifically, these classes allow associating a material mix family with an imported medium component (see the MaterialMixFamily class). In this case, the appropriate material mix from the family is selected by one or more extra parameters specified in the imported file for each particle or cell.

Currently, the only implemented material mix family is the SelectDustMixFamily class, which allows an imported medium component to select one of the user-configured dust mixtures for each particle or cell. Future material mix families may provide more sophisticated options, for example varying the form of the dust grain size distribution depending on parameters in the import file. However, there are important limitations. Because each material mix is represented by a distinct C++ object, a material mix family must by necessity consist of a discrete number of family members. And also, because a material mix object often consumes a significant amount of memory for holding material properties and precalculated values, the number of material mixes that can be offered by a material mix family is limited by the available memory resources.

Finally, a spatial grid must be properly configured to provide adequate spatial resolution in areas with high medium densities or large density gradients. When importing a hydrodynamical simulation snapshot, the best option is to select an adaptive grid, such as a hierarchical octree or binary tree grid (PolicyTreeSpatialGrid) or a Voronoi grid (VoronoiMeshSpatialGrid). To help evaluate and configure the quality of a spatial grid, it is instructive to review the convergence information and the density cuts along the coordinate planes generated by SKIRT upon request (see the SpatialGridConvergenceProbe and DefaultMediaDensityCutsProbe classes). The theoretical density maps shows the medium density as it is defined through the medium components, while the gridded density shows the medium density after discretization by the spatial grid.

## Particle media

### Calculating densities

The mechanisms for loading the medium density from a particle snapshot are similar to those discussed in Particle sources. There are, however, some differences and additional considerations. The total medium density at an arbitrary position is obtained through

$\rho(\mathbf{r}) = \sum_{j=1}^N M_j \, W(|\mathbf{r}-\mathbf{r}_j|,h_j)$

where $$M_j$$ now indicates the mass contribution for each particle. This value may simply be specified in the snapshot as a single value, or it can be derived from multiple values using a simple formula depending on the options configured for the ParticleMedium class.

The formula is provided mostly for compatibility with previous SKIRT versions. It implements a basic heuristic to calculate the dust mass corresponding to the gas mass in the input data, assuming that the amount of dust is proportional to the metal fraction in the gas, except in areas where the gas is too hot to form dust. In other words, dropping the particle index, we have

$M_\mathrm{dust} = \begin{cases} f_\mathrm{dust}\,Z\,M_\mathrm{gas} & \mathrm{if}\; T<T_\mathrm{max} \\ 0 & \mathrm{otherwise} \end{cases}$

where $$M_\mathrm{gas}$$, $$Z$$, and $$T$$ are the gas mass, metallicity, and temperature given by the gas particle's properties in the snapshot, and $$f_\mathrm{dust}$$ and $$T_\mathrm{max}$$ are constant parameters specified when configuring SKIRT.

As for sources, the ParticleMedium class allows configuring one of several (spherically symmetric) smoothing kernels, including the frequently-used cubic spline kernel (the CubicSplineSmoothingKernel class) and a scaled Gaussian smoothing kernel (the ScaledGaussianSmoothingKernel class). Corresponding to the kernels used in SPH simulations, these kernels have a finite support so that only a relatively small number of terms in the sum above have a non-zero contribution. SKIRT uses this property to accelerate the calculation of the medium density at a given spatial position even for a potentially large number of smoothed particles. Essentially, the procedure involves placing an appropriate grid over the spatial domain and constructing a list of overlapping smoothed particles for each cell in that grid. For more information, see the ParticleSnapshot and SmoothedParticleGrid classes.

The ParticleMedium class expects a text file in a format similar to what was described in Spatial distribution and SED from snapshot for radiation sources. Again, the first four columns specify the $$x$$, $$y$$ and $$z$$ coordinates and the smoothing length $$h$$ for each particle. The subsequent columns now specify the medium mass $$M$$, and optionally the metallicity $$Z$$ and temperature $$T$$ (in K) depending on the configured import options.

The ParticleMedium class can also be configured to take into account the Doppler shift caused by the velocity of the medium relative to the model coordinate system when interacting with photon packets. If this option is enabled, the input text file must provide three additional columns specifying the $$v_\mathrm{x}$$, $$v_\mathrm{y}$$ and $$v_\mathrm{z}$$ components of the particle's velocity.

Alternatively, the configuration can include a GeometricMedium with a ParticleGeometry. The expected text file has the same format as described above (without support for velocities). However, the ParticleGeometry class normalizes the total mass of the density distribution to unity, which means that we need to supply the actual total mass through one of the standard normalization mechanisms offered by the GeometricMedium class.

## Mesh-based media

### Calculating densities

Determining the medium density at an arbitrary position in a mesh-based distribution comes down to identifying the cell that contains that position and returning the density associated with this cell. For a hierarchical grid this entails a simple search that starts at the root node and recursively descends into the child node that happens to contain the given position until a leaf node has been reached. With the cuboidal cells in a Cartesian adaptive mesh this is rather straightforward.

In the case of a Voronoi mesh, however, the cell identification is not as simple. Due to the nature of Voronoi tessellations, locating the appropriate cell is essentially a nearest neighbor search. Rather than looping over all possible cells (or equivalently, all generating sites), SKIRT implements a dual level search grid to accelerate the procedure. For more information, see the VoronoiMeshSnapshot class.

The file formats expected by the AdaptiveMeshMedium and VoronoiMeshMedium classes are very similar to those described for radiation sources respectively in Snapshot data format and Loading snapshots. The properties determining the SED for radiation sources are in this case replaced by the properties defining the medium mass contributions and optionally the velocity components for each cell as described in Loading snapshots for particles.

Alternatively, we can configure an AdaptiveMeshGeometry or VoronoiMeshGeometry in a GeometricMedium. These classes normalize the total mass of the density distribution to unity, so that we need to supply the actual total mass through one of the standard normalization mechanisms offered by the GeometricMedium class.

To support various use cases, all of these classes offer a massType option with four possible values, affecting the interpretation of the imported value in the first column:

• Mass (integrated mass density): specifies the total mass in the cell; appropriate for materials such as dust.
• MassDensity (mass density): value is multiplied by the cell volume to obtain the integrated mass in the cell.
• Number (integrated number density): specifies the number of atoms or electrons in the cell; appropriate for gas.
• NumberDensity (number density): value is multiplied by the cell volume to obtain the integrated number in the cell.

### Configuring the spatial grid

For mesh-based hydrodynamical snapshots, SKIRT offers the option to discretize the spatial domain for the radiative transfer simulation itself using the imported grid (see the AdaptiveMeshSpatialGrid VoronoiMeshSpatialGrid classes). In this case, rather than constructing a spatial grid using some configured scheme based on sampling the density in the imported distribution, SKIRT directly adopts the grid on which the snapshot has been defined.

While this seems a logical thing to do, it is not always the best choice. The resolution requirements for the radiative transfer treatment may differ from those for the hydrodynamical simulation. We might not need the same overall resolution, and fine-grained cells may be needed in different areas of the domain. Also, sometimes the radiative transfer spatial grid can be made smaller than the full snapshot domain, cutting off the outer areas that do not contain significant amounts of material anyway. And finally, SKIRT's octree spatial grid (the PolicyTreeSpatialGrid class) is highly optimized for tracing photon packets, so the actual photon shooting simulation phase might be a lot faster.