Wavelengths in SKIRT

Introduction

With few exceptions, the cross sections for interactions between matter and radiation depend heavily on the radiation's wavelength. Furthermore, some processes affect the wavelength of an interacting photon (e.g. Doppler shifts caused by relative motion) or drastically "convert" a given amount of radiative energy from one wavelength to another (e.g. absorption and re-emission by dust grains). As a result, wavelength plays an important role in a radiative transfer simulation code such as SKIRT.

In general, the process goes like this. When a photon packet is launched, it receives a certain amount of energy at a given wavelength. While the packet moves through the spatial domain of the simulation, some of its energy is absorbed by the transfer medium, and its wavelength may be adjusted by interactions with the medium. At each scattering event, a fraction of the packet's energy is peeled-off towards each instrument to record the observed flux. The energy absorbed by the medium is re-emitted at a later stage by launching new photon packets at different wavelengths.

This brief summary passes over a lot of details and complexities, but even so it hints at the many areas in the simulation where wavelength-related aspects need to be handled and configured. On the other hand, not every simulation requires this full general treatment. For example, when studying dust extinction in the optical wavelength range, there is no need for handling dust emission.

This concept note describes many of these issues in more detail:

Wavelength regimes

With respect to the treatment of wavelengths, SKIRT has two distinct simulation modes called panchromatic and oligochromatic. One of these modes must be selected early on in the simulation configuration process, and the choice affects many of the other options in the configuration.

Panchromatic simulation mode

Panchromatic simulation mode is the most powerful and commonly used mode in SKIRT. A panchromatic simulation operates with a continuous range of wavelengths that often spans ultraviolet, optical, infrared and submillimeter bands. Such a simulation supports the full SKIRT feature set, including for example thermal emission by dust grains and the effects of kinematics. If the wavelength range is limited to (part of) the ultraviolet and/or optical spectrum, it obviously doesn't make sense to include thermal emission, but the other features (such as kinematics) remain available.

In panchromatic mode, the user must specify a primary source wavelength range in which the SED assigned to each source must be considered, and configure wavelength grids for instruments, radiation field, and dust emission as discussed in the remainder of this text.

Oligochromatic simulation mode

Oligochromatic simulation mode is used for studying extinction effects (i.e. absorption and scattering) at a limited number of specific wavelengths (or just a single wavelength) in the ultraviolet and/or optical spectrum. The simulation operates at these discrete wavelength values only, and is therefore unable to calculate wavelength-integrated quantities such as dust temperatures or thermal emission, or to properly handle wavelength shifts such as those caused by kinematics. On the other hand, because of the small number of wavelengths in use, the total number of photon packets launched can usually be substantially smaller than for panchromatic simulations. In other words, one trades speed for a reduction of the supported feature set.

In oligochromatic mode, the user must specify the list of discrete wavelengths to be used in the simulation. The luminosity of each source is still derived from its associated SED or SED family, and luminosity normalization proceeds as usual (and may even happen at wavelengths that differ from the oligochromatic wavelengths). However, newly emitted photon packets are automatically assigned one of the distinct oligochromatic wavelengths. Furthermore, all instruments are forced to use a special wavelength grid with narrow bins around the oligochromatic wavelengths.

Photon packet properties

Wavelength and weight

The fundamental properties of a photon packet include its wavelength and its weight. The wavelength property specifies the wavelength (or equivalently, the frequency) of all photons in the packet. The weight property specifies the number of photons carried by the packet, or more precisely the number of photons per unit of time (because SKIRT solves the time-independent radiation transfer equation). Refer to the PhotonPacket class for more details.

At launch (from a primary or secondary source), a photon packet receives a wavelength sampled from the source spectrum and a luminosity, i.e. its share of the total luminosity of the source. The wavelength is stored as given. The luminosity is converted to a weight (number of photons) for storage in the photon packet.

Note
We use the term luminosity to indicate energy per unit time carried by a photon packet, although, strictly speaking, it should only be used to indicate energy per unit time emitted by a source.

During a photon packet’s life cycle, updates can occur to both its weight, e.g. because of biasing, and its wavelength, e.g. after being scattered by a moving medium. Because these updates can be fractional, both wavelength and weight are stored as floating-point values. Adjusting a photon packet's wavelength indirectly affects the luminosity represented by the packet, because the latter is directly proportional to the frequency and thus inversely proportional to the wavelength.

Kinematics

The wavelength of a photon packet is defined relative to the model coordinate system. In other words, a medium at rest relative to the model coordinate system sees this wavelength. Velocities of sources and media are also defined relative to the model coordinate system. The bulk velocity of the aggregate medium in a spatial grid cell is determined by averaging the bulk velocities of all media components weighed by density. Instruments are considered to be at rest relative to the model coordinate system.

When a photon packet is launched, its wavelength is Doppler shifted according to the component of the source velocity in the photon packet’s direction. When a photon packet interacts with a medium, the perceived wavelength is derived by Doppler shifting the packet’s wavelength according to the component of the medium velocity in the photon packet’s incoming direction. Specifically, registration of a photon packet’s contribution to the radiation field uses this perceived wavelength. After a scattering interaction, the photon packet’s wavelength is replaced by the perceived wavelength, Doppler shifted according to the component of the medium velocity in the photon packet’s outgoing direction.

Emitting photon packets

Distributing photon packets over sources

The source system (see the SourceSystem class) distributes the emitted photon packets over the primary sources in the simulation based on three types of information: the luminosity weight of each source, a user-configured weight for each source (defaulted to equal weight for all sources), and a user-configured bias factor (defaulted to 0.5) that shifts between distributing according to the above weights and simply allocating an equal number of packets to each source. The default values cause half of the photon packets to be distributed proportionally to the luminosity of the sources, and the other half is distributed equally over the sources. Changing the source weights allow a user to assign more importance to particular sources, and the bias factor can be adjusted to swing between the proportional and linear allocation schemes.

In turn, each source (see the Source class) determines how to further distribute the photon packets allocated to it. For example, a geometric source distributes photon packets according to its (built-in) spatial density distribution. An imported source has a scheme similar to the source system, using a user-configured composite bias factor to shift between a distribution according to the luminosity weights of each sub-source (particle or cell) and a linear equal-weight allocation.

The secondary source system (see the SecondarySourceSystem class) also uses the latter scheme to distribute photon packets among the cells of the spatial grid configured for the simulation.

Sampling wavelengths

Each radiation source provides a mechanism to sample a random wavelength from its spectral energy distribution (see the SED class and its subclasses) and assign a corresponding luminosity such that the bolometric luminosity of the source is properly distributed across the photon packets. All currently implemented SEDs employ the inversion method using the tabulated cumulative distribution to sample wavelengths. It would be possible, however, to employ a specialized sampling routine for analytically defined SEDs.

It is important that all aspects of the SED are properly sampled, for example wavelength ranges with low luminosities or narrow line features. To this end, the sampling procedure employs composite biasing to combine sampling from two distributions: the SED itself (favoring wavelengths that carry a lot of energy) and a custom bias distribution (favoring specific wavelength areas). Both the composite bias factor and the bias distribution are user-configurable. By default, the procedure has a bias factor of 0.5 and a bias distribution that is logarithmic in wavelength.

The default scheme ensures that the low-luminosity tails of a typical spectrum are properly sampled, while still favoring the higher-luminosity areas. Even narrow spectral features are properly sam- pled because half of the wavelengths are selected from the source SED at full spectral resolution. Lowering the bias fraction focuses more photon packets into high-luminosity areas because the composite distribution more closely follows the source spectrum. Vice versa, a bias fraction close to unity causes the source spectrum to be essentially ignored for the purpose of wavelength sampling.

The default bias distribution is usually appropriate for wavelength ranges spanning multiple decades, where one aims for a constant spectral resolution over the entire range (modulated with the source spectrum as per the bias fraction). For narrow wavelength ranges, perhaps corresponding to a particular spectrograph or spanning a given emission line, a linear distribution of the photon packet wavelengths might be more appropriate. To this end, SKIRT offers a built-in uniform wavelength bias distribution in addition to the logarithmic distribution. Users can also load a custom distribution from file for maximum flexibility. For example, one might want to strongly favor a wavelength interval of special interest even if the model's sources are not particularly luminous in that interval.

Types of wavelength grids

Wavelength discretizations

As described in a previous section, photon packets can have an arbitrary wavelength value in some continuous range. However, the only practical option in other areas of the simulation, such as defining material properties or recording fluxes, is to discretize the wavelength range. These various discretizations are essentially uncoupled, which leads to a potentially large number of different wavelength grids, each specialized for a particular purpose:

• Wavelength-dependent properties of transfer media are tabulated in resource files on some private wavelength grid.
• The SEDs or SED families assigned to primary sources are tabulated on some private wavelength grid.
• A radiation field wavelength grid is required to register the radiation field during primary and secondary emission.
• A dust emission grid is required to store the dust emission spectrum calculated for each cell during secondary emission.
• Each instrument requires a wavelength grid for binning detected photon packets.

While some of these discretizations are defined by the input resource (e.g. material properties and SEDs), the other discretizations (most notably those for the instruments) must be configured by the user. To this end, SKIRT includes a number of wavelength grids as presented below.

A band object in SKIRT represents the transmission curve of a particular observational filter as a function of wavelength. Key operations offered by all band objects include obtaining the transmission at a given wavelength and calculating the mean specific luminosity for a given SED after convolution with the transmission curve. SKIRT offers a set of built-in band objects for standard filters, such as the Johnson filters, and for common observatories, such as GALEX, SDSS and Herschel. Refer to the documentation of the BroadBand class for a list of supported bands. Other bands can be loaded from file, given a tabulated transmission curve.

A band object can be used to normalize the luminosity of a source to a given mean specific luminosity for the band, which often corresponds more precisely to an observed quantity than specifying a specific luminosity at a particular wavelength. For more information, see the NormalizedSource and BandLuminosityNormalization classes.

More interestingly, perhaps, it is also possible to equip an instrument with a "wavelength grid" built from a list of (possibly overlapping) bands. In that case, each band represents a separate bin of the instrument. When a photon packet arrives, its contribution is multiplied by the transmission at the packet's wavelength for each band before being accumulated in the corresponding bin. This amounts to "on-the-fly" convolution of the detected flux with the transmission curve of each band.

A radiative transfer simulation is often performed with the aim of comparing its results with observations. In that case, using a band wavelength grid produces directly comparable output. The alternative is to run the simulation using a regular wavelength grid with fairly narrow bins, and perform the convolution after the fact. For proper results, the instrument wavelength grid must resolve all spectral features of the sources, including emission or absorption lines which may be Doppler shifted because of kinematic effects. This may require a large number of bins, with correspondingly large memory requirements.

Wavelength grid classes

The following diagram shows a portion of the class inheritance tree for wavelength grids (connections starting with a triangle denote inheritance).

The DisjointWavelengthGrid class represents wavelength grids with non-overlapping bins and constant transmission across each bin. There are currently two kinds of disjoint wavelength grids. The first kind has adjacent bins discretizing a wavelength range by filling it with consecutive bins that, together, cover all wavelengths in the range. The LogWavelengthGrid and NestedLogWavelengthGrid classes are examples of this kind of grid. The second kind has a number of distinct, non-adjacent wavelength bins each enveloping a characteristic wavelength with a bin width that is smaller than the separation between the characteristic wavelengths. The FileWavelengthGrid class can be configured to represent either adjacent or non-adjacent wavelength grids. For more information, refer to the documentation of these respective classes.

The BandWavelengthGrid class represents wavelength grids where each bin is defined by the transmission curve of a particular broadband. The intervals in which the transmission is nonzero are allowed to overlap, so the bins are not necessarily disjoint. The PredefinedBandWavelengthGrid class offers a list of predefined bands comprising the GALEX, SDSS, 2MASS, WISE and HERSCHEL broadbands. Each of these sets can be included or excluded through configuration flags. The ConfigurableBandWavelengthGrid class (not shown in the diagram above) allows configuring an arbitrary list of bands.

Configuring wavelength grids

Instrument wavelength grids

As mentioned above, each SKIRT instrument requires a wavelength grid for binning the observed fluxes contributed by detected photon packets. Each instrument can be configured with its own distinct wavelength grid. To avoid repetition in configurations with several instruments, a default wavelength grid can be configured for the instrument system. The default grid is adopted by instruments for which no individual grid has been specified.

Instruments can use any of the wavelength grids discussed in the previous section. For disjoint grids the binning process is straightforward. In case a BandWavelengthGrid has been specified, an arriving photon packet is registered for each band in the wavelength grid after multiplying its luminosity contribution by the band’s transmission factor corresponding to the packet’s wavelength.

One useful technique is to configure two instruments for the same line of sight: an SEDInstrument for recording a high-resolution spectrum, and a FrameInstrument for recording a small set of broadband images. Each instrument performs different binning on the same set of arriving photon packets: the SED instrument spatially integrates fluxes into narrow spectral bins, while the frame instrument does a spectral convolution for the fluxes in every spatial pixel. This leads to a very efficient use of the photon packets being traced through the system. To achieve a similar result with previous SKIRT versions, one would configure a single instrument that records a data cube with both the required spatial and spectral resolution, and perform the two-way binning after the fact. Apart from requiring an extra processing step, the three-dimensional data structures in this approach can become very large.

Internal wavelength grids

For a panchromatic simulation including dust emission, SKIRT requires the user to configure two wavelength grids that affect its internal operation: the radiation field wavelength grid and the dust emission wavelength grid. Only disjoint wavelength grids with adjacent bins can be used for these purposes. Future implementations may allow determining these grids automatically using some heuristic. In the meantime, this section offers some guidance.

To focus the discussion, we consider a panchromatic SKIRT simulation that produces synthetic observations from UV to submm wavelengths for a typical hydrodynamcially simulated galaxy with a stellar mass of $$1.75\times 10^{10}~\mathrm{M}_\odot$$ represented by more than 125000 stellar particles, and a dust mass of $$3.8\times 10^{7}~\mathrm{M}_\odot$$ derived from over 20000 gas particles. These numbers are sufficiently large for the radiative transfer simulation to be nontrivial, without being a limiting factor for running several tests with varying configuration parameters. To properly resolve the spatial structure of the input model, we configure a spatial octree grid with just over one million cells.

The radiation field wavelength grid defines the bins used to record the energy deposited by photon packets in each spatial cell. Generally, we expect the radiation field at shorter wavelengths to dominate the dust heating process, with the longer wavelengths having a minimal effect. We can also presume that the precise wavelength of an incoming photon packet might not be so important, as long as its energy is properly categorized, possibly allowing fairly wide wavelength bins. Convergence tests indicate that, for the simulation described above, it is sufficient to configure a radiation field wavelength grid with just 40 points from 0.02 to 10 micron, distributed evenly in log space (using a LogWavelengthGrid). This is an important result, because the memory requirements for a SKIRT simulation critically depend on this number of bins.

The dust emission wavelength grid controls the resolution of the dust emission spectrum calculated for each spatial cell. Especially when taking into account the stochastic heating of small dust grains, this spectrum contains many narrow infrared features. It is thus desirable to configure a grid that can properly resolve these features. Memory usage is not an issue because the emission spectrum is stored just once per execution thread. The performance impact is very limited as well because sampling these emission spectra is not the bottleneck of the calculation. For the simulation described above, we configure an emission wavelength grid with a resolution of 100 bins per dex in the overall range from 0.2 to 2000 micron, and 200 narrower bins in the range from 3 to 25 micron, for a total of 508 bins (using a NestedLogWavelengthGrid).

Note
Before adopting these or similar wavelength grids for other SKIRT simulations, appropriate convergence tests should be performed.

Memory usage

There are obviously countless areas in the SKIRT code that consume memory. In many cases, however, the overall memory requirements are dominated by just a few components, namely the radiation field storage and the instrument data cubes.

Consider a typical panchromatic dust continuum simulation as the one described in the previous subsection. We assume that the radiation field wavelength grid has $$N_{\lambda,\mathrm{rf}}$$ points and that the spatial grid has $$N_\mathrm{cell}$$ points. The data structure for storing the radiation field then has a size in bytes of

$S_\mathrm{rf} = 8 N_{\lambda,\mathrm{rf}} N_\mathrm{cell}$

If the simulation supports dust self-absorption, it stores the absorbed energy from stellar and dust emission separately, so that the data structure has twice this size.

We further assume a frame instrument recording in $$N_{\lambda,\mathrm{ins}}$$ wavelength bins and $$N_\mathrm{x}$$ by $$N_\mathrm{y}$$ image pixels. The data structure for storing the detected fluxes has a size in bytes of

$S_\mathrm{ins} = 8 N_{\lambda,\mathrm{ins}} N_\mathrm{x} N_\mathrm{y}$

If the instrument is requested to keep track of individual flux components (such as direct and scattered light), or if there are similar instruments at other viewing angles, the required memory becomes a multiple of this size.

Previous SKIRT versions used a global wavelength grid for all purposes throughout the simulation. In other words, in the above equations, $$N_{\lambda,\mathrm{rf}}=N_{\lambda,\mathrm{ins}}=N_\lambda$$, and $$N_\lambda$$ must be sufficiently large to accomodate the requirements of all areas in the code. With one million spatial cells, 500 wavelengths, and a single 750 x 750 pixel frame instrument, the total size of the discussed data structures is 5.8 GB.

Because wavelength grids are now uncoupled, memory requirements can be dramatically reduced. Assume that we configure a radiation field wavelength grid with 40 points as the previous section, so that $$N_{\lambda,\mathrm{rf}}=40$$, and that we require image frames for 20 broadbands, so that $$N_{\lambda,\mathrm{ins}}=20$$, the total size of the data structures now becomes 0.4 GB.

To be fair, it should be noted that some studies will need a spatially and spectrally resolved instrument data cube in a given wavelength range, for example to evaluate the effects of kinematics, or to simulate integral-field spectroscopy observations. This affects $$N_{\lambda,\mathrm{ins}}$$ without changing $$N_{\lambda,\mathrm{rf}}$$ in the equations above. When we introduce other media types such as hydrogen gas, the wavelength resolution $$N_{\lambda,\mathrm{rf}}$$ of the radiation field storage will need to increase as well.

Discrete panchromatic simulation mode

In the regular – and highly recommended – panchromatic simulation mode, SKIRT samples all wavelengths in a given wavelength range. In contract, previous SKIRT versions and many other codes only sample a single characteristic wavelength within the bin. Depending on the form of the involved spectra, these two methods may (and often will) produce different results. While the current SKIRT method is designed to more closely reflect the physical reality, these differences can be very annoying when comparing SKIRT results with those of other radiative transfer simulation codes as part of, for example, a benchmark effort.

SKIRT therefore includes the DiscreteWavelengthDistribution class, which creates the ability to mimick codes that emit photon packets at discrete wavelengths rather than across a continuous range, while still using panchromatic simulation mode. Configuring such a simulation, however, involves some advanced options and requires great care, as described in more detail below.

When a DiscreteWavelengthDistribution instance is used as the wavelength bias distribution for a source with a composite bias factor of one, the source will emit photon packets only at the characteristic wavelengths of the configured grid, and the photon packets will be distributed with equal probability among those wavelengths. In this situation, it makes little sense to configure wavelength grids with different bins for detecting photon packets in other areas of the simulation. Doing so might cause unexpected results. For example, some bins might receive no photon packets at all and thus incorrectly report zero influx. It is therefore best to configure the same wavelength grid throughout the simulation.

To configure a simulation that uses discrete wavelengths in panchromatic mode, follow these guidelines:

• Choose a wavelength grid with adjacent bins (i.e. representing a continuous range), such as a LogWavelengthGrid or a FileWavelengthGrid with relativeHalfWidth set to zero.
• Configure all primary and secondary sources (including, e.g. dust emission) with a wavelength bias factor of one, and a wavelength bias distribution of type DiscreteWavelengthDistribution using the chosen wavelength grid.
• Explicitly configure all wavelength grids in the simulation to be the same chosen wavelength grid. This includes the radiation field wavelength grid, the dust emission wavelength grid, and the wavelength grids for all instruments and probes.
• Avoid using kinematics in the simulation. Small Doppler shifts from the emitted discrete wavelengths will stay within the same wavelength bin anyway. Larger shifts will cause photon packets to move between wavelength bins in a discontinuous way.