DirtyGrid I: 3D Dust Radiative Transfer Modeling of Spectral Energy Distributions of Dusty Stellar Populations

Due to the ability of dust to absorb, scatter, and re-emit radiation, deriving the global SED of a stellar population mixed with interstellar dust is a radiative transfer problem. The three-dimensional (3D) radiative transfer problem can be represented by an integral-differential equation (Steinacker et al. 2013). However, in the general case, this equation has no analytic solution. We use DIRTY, a 3D Monte Carlo model to compute the radiative transfer of stellar, gas, and dust emission through a dusty medium. DIRTYv1 is described in detail in Gordon et al. (2001) and Misselt et al. (2001). The current version of the code (DIRTYv2) is based on the algorithms in DIRTYv1, but with the addition of a full dust grain model, a tighter integration of the dust radiative transfer and emission calculations, improved robustness of the input/output, and the enhancement of continuous photon absorption to accelerate the calculations. Here we provide a summary of DIRTY and detail the new features in the current version. Readers interested in the radiative transfer details should refer to the original DIRTY papers and Steinacker et al. (2013).

A single run of a DIRTY model requires three inputs: (1) the energy sources, (2) the physical properties of the absorptive medium (e.g., the dust grain physics), and (3) the relative distribution of both dust and stars: the model geometry. In the context of a DIRTY model, "geometry" refers to both the relative distribution of sources and dust throughout the model space (global geometry), as well as the clumpiness of the dust distribution (local geometry). By construction, in DIRTY all three ingredients may be arbitrarily specified. In general, the inputs are specified based on the physical problem at hand. While the local geometry can be of arbitrary complexity (e.g., independent dust densities at each point in the model space), it is convenient to specify the local distribution of dust ("clumpiness") by a two-phase medium characterized by the filling factor of high density regions embedded in a uniform medium and the density ratio between the high density regions and the uniform medium (Witt & Gordon 1996).

DIRTY has two main components: a radiative transfer routine and a dust emission routine. Given a distribution of emission sources within an absorptive medium, the purpose of the radiative transfer routine is to calculate the distribution of absorbed and escaped fluxes according to the physics of the interaction between photons and dust grains and the relative distribution between the sources and the dust. The same radiative transfer routine is used for both primary emission (from stars and gas) and secondary emission (from dust). The dust emission routine calculates the emission spectrum from the absorption spectrum according to the physics of dust heating and cooling. The fact that dust emission can be re-absorbed by dust during radiative transfer introduces a circular dependency between the dust absorption and emission spectra, so DIRTY takes an iterative approach to obtain a self-consistent solution. The iterative procedure is illustrated in Figure 2. This contrasts with the simpler way of modeling dust heating with a radiation field scaled from observations in the solar neighborhood (e.g., Draine et al. 2007).

Zoom In Zoom Out Reset image size

The DIRTY radiative transfer algorithm utilizes a weighted Monte Carlo approach with absorption-scattering split, forced first scattering, the peel-off technique, and continuous absorption, which are summarized below. We use a Cartesian local mean intensity storage grid to discretize the model geometry into regions of constant optical depth and radiation fields (see Section 2.4 for details). Steinacker et al. (2013) reviewed 3D dust radiative transfer techniques in general (not specific to DIRTY) and devoted a section to each of the performance enhancing techniques we list above. For the same level of accuracy achieved, these techniques allow us to reduce the amount of computation required, which is critical for a large scale project like the DirtyGrid. Alternatively, for the same computation time, one could use these techniques to enhance the accuracy of the results.

In the weighted Monte Carlo approach, instead of following individual photons as discrete particles, we assign weights to photon packets and keep track of the changes in the weights due to interaction events. We trace the trajectory of a photon packet and determine the distance traveled according a random number sampled from a probability distribution based on the optical depth in the path. At each interaction site, we split the weight of the photon packet into an absorbed portion and a scattered fraction according to the dust albedo (the absorption-scattering split). The direction of scattering is determined from a scattering phase function. Instead of depositing the entire absorbed portion at the interaction site, we distribute the energy along the path. This is known as continuous absorption and is new to DIRTYv2. In addition, we allow the photon packet to contribute to the escaped flux at each interaction site to improve the spatial sampling (the peel-off technique). In low optical depth models, a large number of photon packets may exit the model geometry without any interactions. To avoid "wasting" photon packets, we force the first scattering so that all of them contribute to the absorbed and scattered fluxes, and we correct for the bias in probability by adjusting the weights.

For a fixed number of interaction events, continuous absorption greatly improves the sampling of the radiation field while keeping the sampling of scattering the same. However, this calculation is relatively expensive. As a result, for a fixed amount of computation, this technique improves the determination of the local radiation field but worsens the sampling of scattering. In other words, it shifts the computational balance to focus more on populating the radiation field. For the parameter space covered by the DirtyGrid, this is beneficial because usually the scattered fluxes are better populated than the radiation field, and stochastic heating of dust grains is very sensitive to noise in the radiation field. Conversely, continuous absorption may not be suitable for models of extremely high optical depth or models that only calculate the equilibrium emission of dust grains.

We repeat the radiative transfer for a number of photon packets until the Monte Carlo uncertainty in the direct and scattered fluxes is below a desired limit (0.5%–1.0%, depending on the geometry), or until the maximum limit of 10 million photon packets per wavelength has been reached. We then repeat the process for all the wavelengths in a non-uniform wavelength grid with 174 points. This wavelength grid is composed of a uniform grid in log wavelength between 0.0912 and 1000 μm, with extra points placed to sample where there are features in the dust absorption, scattering, or emission (e.g., the mid-infrared for the PAH emission features and silicate absorption features). Including more points would not change the solution, as the relevant dust features are well sampled, while including fewer points would result in an inaccurate solution.

Since DirtyGrid is intended to provide a grid of integrated galactic SEDs over a large range of properties, the sources for the radiative transfer model are stellar population models. Stellar population models calculate the integrated SEDs of stellar populations optionally including the emission due to ionized gas. These SEDs determine the intensity and wavelength dependence of primary emission in the context of radiative transfer. Examples of stellar population models in the literature include PEGASE.2 (Fioc & Rocca-Volmerange 1997,1999), Starburst99 (Leitherer et al. 1999), GALAXEV (Bruzual & Charlot 2003), and FSPS (Conroy et al. 2009). We select PEGASE.2 (hereafter PEGASE) for the DirtyGrid due to its ability to model both young and old stars well, its flexibility, and our previous experience in using it with DIRTY. For example, Gordon et al. (1999) have used PEGASE with DIRTY to model the dusty starburst nucleus of M33. We expect that using a different stellar evolutionary synthesis model would produce somewhat different output SEDs, but the changes are expected to be fairly small when compared to the overall shape and intensity of the SEDs. We plan to investigate using different stellar evolutionary models in future work. In the DirtyGrid, we vary four of the parameters of the stellar population model: the age of the stellar population and its formation history, the metallicity, and the SED scaling factor. Figure 3 shows an example of PEGASE SEDs.

Zoom In Zoom Out Reset image size

PEGASE uses the isochrone synthesis method, which integrates an IMF, evolutionary tracks, and a library of stellar spectra, to produce the integrated SEDs of a stellar population. In a single star formation event, the IMF describes the relative abundance of newly formed stars of different masses. We used the Salpeter (1955) IMF from 0.1 to 120 M_⊙. This IMF was used, as it is a common one to assume for SFR studies and one of the applications of the DirtyGrid is to investigate the various SFR recipes. Kroupa (2001) has found variation in the IMFs between globular clusters and the Galactic-field, but it is uncertain how or whether the IMF change in different environments and redshifts. Adopting a cutoff in the low mass range of the IMF can lead to a factor of ∼1.5 change in the SFR determined from empirical SFR indicators (Kennicutt & Evans 2012).

The evolutionary tracks describe the trajectory of stars in the Hertzsprung-Russell diagram. PEGASE mainly uses tracks from the Padova group, but they also supplemented the Padova tracks in a number of special evolutionary phases. Most importantly, they use the equations from Groenewegen & de Jong (1993) to create pseudo-tracks for the thermally pulsing asymptotic giant branch (TP-AGB) stars, which are short-lived but luminous.

In general, the SFR of a galaxy is a function of time. The typical timescale for gas depletion is about 3 Gyr (Pflamm-Altenburg & Kroupa 2009), so it is reasonable to expect the SFR to decrease over time. The time-dependent profile of star formation can be modeled as an exponentially decaying burst (Searle et al. 1973; Conti et al. 2003). However, given limited computational resources, we choose to model only the following two extreme cases: instantaneous burst and constant star formation, which correspond to time constants of infinity and zero, respectively. These two cases bracket the possible real world scenarios. For example, the fraction of young stars in an exponentially decaying star formation model is lower than that of a constant star formation model, and higher than an instantaneous burst model, for any finite decaying time constants. Our choice reduces the "star formation history" dimension of our parameter grid to only two points. With further increase in the capability of supercomputers and more data to constrain the parameter space, we may include exponentially decaying star formation into a future version of DirtyGrid, potentially increasing the number of points in the "star formation history" dimension by an order of magnitude.

Closely related to the star formation history is the age of the stellar population. The age of the stellar population reflects the time elapsed since the "instantaneous" burst of start formation in the case of a burst model, or the time over which the model has been forming stars at a constant rate. To cover a range of plausible ages for the stellar populations, the DirtyGrid includes both burst and constant stellar populations with ages between 1 Myr and 13 Gyr.

The integrated SED of a stellar population depends on the initial metallicity of the stellar population. Populations with different metallicity trace different evolutionary tracks and populate different isochrones. Hence population synthesis models, including PEGASE, also specify the metallicity of the stellar population. To cover a realistic range of integrated galactic metallicities, in DirtyGrid we compute SEDs for stellar populations spanning a range of metallicities between Z of 0.0001 and 0.1, where solar metallicity is ∼0.013.

To control the intensity of stellar emission, we scale the input SED by the stellar surface density. For instantaneous star formation (ISR) models, the scaling factor is the stellar mass surface density ${\sigma }_{M}$. For constant star formation (CSR) models, the scaling factor is the SFR surface density ${\sigma }_{\mathrm{SFR}}$. The total mass of stars formed in CSR models is the product of ${\sigma }_{\mathrm{SFR}}$ and the age of the stellar population. In both cases, the final stellar mass depends on the evolution of the stellar population, which is an output of PEGASE, and due to the finite life time of stars is always smaller than the total mass of stars formed. We choose the range of stellar mass or SFR by comparing the model outputs to SINGS (Kennicutt et al. 2003) and LVL galaxies (Dale et al. 2009) and star formation regions in very nearby galaxies. The range we use covers the range derived empirically from these observations.

DIRTY does not calculate the radiative transfer through gas, which would add significant complexity to the code. The effect of gas, which is particular important for young stellar populations due to the significant UV radiation, is approximated by a gaseous nebula that is local to the stars and optically thick in the Lyman lines. Under extremely intense UV radiation, the absorption rate of Lyman continuum photons by gas is limited by the number of neutral hydrogen atoms, so some of these photons should be absorbed by dust instead. In contrast, in our models these photons are always absorbed by gas. The difference is negligible in typical H ii regions (Walcher et al. 2011). We add the gas continuum and recombination lines emission, as calculated by PEGASE, to the stellar emission.

The combined spectrum at a given age, star formation history, metallicity, and intensity defines the radiation sources used as input for the DIRTY radiative transfer model.

The dusty medium through which source photons propagate in a radiative transfer simulation has a profound effect on the output spectrum of the simulation. The optical depth, albedo, and scattering phase function of the dust determines the probability that a photon will interact, as well as the amount of the source energy that is absorbed or scattered. Once absorbed, the radiation is thermalized and re-emitted according to the emission cross sections of the dust population, producing a dust emission spectrum.

To integrate the dust population into the radiative transfer simulation, a dust grain model is required. A dust model specifies the composition, size distribution, and micro-physics (e.g., the optical properties) of the dust grain population. Since the DirtyGrid is designed to provide a grid of SEDs representative of the integrated SED of galaxies over a wide range of galactic parameters, the dust inputs must reflect the range of dust characteristics observed in galaxies. Therefore, with DirtyGrid, we adopt the dust grain models described in Weingartner & Draine (2001), who used an analytical expression to parametrize the size distributions of "astronomical" carbonaceous and silicate grains. We summarize the Weingartner & Draine (2001) models we use and refer the interested reader to that paper for the details. They fit the extinction curves of sightlines in the MW, Large Magellanic Cloud (LMC), and Small Magellanic Cloud (SMC) to produce sets of best-fit parameters of the size distributions. For the MW, we pick the R_V = 3.1, ${b}_{C}=6.0\times {10}^{5}$ model; for LMC, we choose the ${b}_{C}=1.0\times {10}^{5}$ model of the "LMC 2" environment (Meaburn 1980; Misselt et al. 1999, the supergiant shell on the southeast side of 30 Dor); for SMC we take the model of the "SMC bar" environment (Gordon & Clayton 1998), which is the only model presented for the SMC. The grains are spherical, and their radii range from 3.5 Å to 2.5 μm. The spherical grain approximation enables the use of Mie theory to compute the optical scattering and absorption cross sections. The silicate grains use the same properties for all grain sizes, unlike the carbonaceous grains. The smallest carbonaceous grains have PAH properties, the largest grains have graphite properties, and intermediate sizes have properties that are a mixture of the two. By including MW, LMC, and SMC dust grain models in the DirtyGrid, our SEDs represent the full range of dust properties that have been observed in galaxies. Figure 4 gives the extinction curves for the three different models of dust grains used in the DirtyGrid.

Zoom In Zoom Out Reset image size

The radiative transfer routine in DIRTY calculates the total energy absorbed by dust in each of the grid cells in our local mean intensity storage grid. DIRTY then launches the dust emission routine to calculate the emission spectrum in each grid cell independently (and sequentially). For each species of dust grain (carbonaceous and silicate), starting at the smallest grain size, we calculate the stochastic heating of dust grains. The calculation of stochastic heating is computationally expensive, but various authors (e.g., Leger & Puget 1984) have shown that stochastic heating is needed to explain the near/mid-IR emission in excess of that expected from dust in pure thermal equilibrium in many astronomical systems. Using the method of Guhathakurta & Draine (1989), we populate a matrix with the probabilities of transitions between different internal energy states, each of which has an associated temperature. The solution of the transition matrix gives us the probability of finding the dust grain in each of the states, P(T). The emission at this particular grain size is an integration of the modified blackbody function (taking into account the effects of the changing emissivity across wavelengths) weighted by P(T). As the calculation moves to grains of increasing sizes, we reach a regime in which the grains are in equilibrium with the radiation field and P(T) becomes strongly peaked. At this point we switch over to the simpler equilibrium heating calculation, in which grains emit in a modified blackbody function with a single temperature. We find the equilibrium temperature, and therefore the equilibrium emission as a function of wavelengths, by simply balancing the energy absorbed and emitted by dust grains of this particular size. In the end, we sum over the emission from all grain sizes and species to obtain the total dust emission spectrum in the current grid cell. A rigorous discussion on dust modeling in DIRTY is available in Misselt et al. (2001). Recently, Camps et al. (2015) compared the dust emission algorithms of six radiative transfer codes, including our latest code, and found the all six codes consistent within 10%.

In addition to the physics that describes how photons will interact with dust and how the dust will respond to the input of energy, the total amount of dust in the model must be specified. Rather than directly specifying the mass of dust, in the DirtyGrid, the optical depth at V band (${\tau }_{v}$) is specified as an input parameter to DIRTY. ${\tau }_{v}$ is measured from the center of the model to the edge, averaged over all $4\pi$ steradians of solid angle. For normal disk galaxies, Holwerda et al. (2007) has found that ${\tau }_{v}$ is on the order of unity. Dale et al. (2006) also found that the average attenuation is ${A}_{V}\approx 1.0$ (i.e., ${\tau }_{V}\approx 0.92$) for a large portion of SINGS galaxies. In the DirtyGrid, the ${\tau }_{v}$ dimension in the parameter space span over 2 dex of optical depths from 0.1 to 10, which we believe is sufficient to cover the majority of observed galaxies.

In the DIRTY model, the global geometry describes the relative distributions of dust and stars, while the local geometry determines the clumpiness of dust. Geometry has a profound effect on radiative transfer, both on the attenuation (e.g., by providing paths of different optical depths) and the emission (e.g., through the non-trivial spatial and temperature distribution of dust). While an extinction curve can be accurately measured by comparing the spectrum of a reddened and an unreddened star of the same spectral type (and has been done for sightlines in the MW, LMC, and SMC), in the complex geometries generally associated with galactic environments, only "attenuation" curves can be measured. A dust attenuation curve is defined as the wavelength dependence of attenuation by dust, including the effects of absorption, scattering out of the beam, and scattering into the beam, and is usually measured for regions of galaxies or whole galaxies. In contrast to the traditional extinction curve, determination of attenuation curves is more challenging because it is difficult to independently constrain the underlying stellar populations and radiative transfer effects (anisotropic scattering, differential extinction between sources, dust exposed to different radiation field intensities, and dust clumpiness). Baes & Dejonghe (2001) found that scattering approximations such as isotropic scattering can be a major source of error in radiative transfer. Even in the absence of scattering, Disney et al. (1989) found that the choice of model geometry significantly affects attenuation.

Commonly used global geometries for radiative transfer modeling of galactic environments include plane-parallel slabs, spheres, double-exponentials, and arbitrary distributions (see Calzetti 2001 and references therein). Plane-parallel slabs are simplest to model, but Calzetti et al. (1994) found that such a model does not agree with the UV-optical spectra of starburst galaxies, and they empirically calibrated a single attenuation curve of starburst galaxies as a third-degree polynomial in $1/\lambda$. To account for that fact that young stars are likely to be embedded in their birth cloud and therefore more attenuated then old stars, Silva et al. (1998) and Charlot & Fall (2000) created models that provide additional opacity around the young stars.

In the DirtyGrid, we use the Cloudy, Dusty, and Shell global geometries, as defined by Witt & Gordon (2000), that are described below. These geometries are spherically symmetric and efficient for radiative transfer. DIRTY supports arbitrary stellar and dust geometries, and we implement the dust geometries as spatially resolved blocks in our local mean intensity storage grid, which has a resolution of 30 × 30 × 30. As we embedded these spherical geometries in a cubical grid, the grid cells outside the sphere are empty and explicitly not traversed during the radiative transfer calculation. A visualization of the cross section of the geometries is shown in Figure 5. In the Cloudy geometry, stars are uniformly distributed in the sphere, and a dusty core is located within 0.69 times of the system radius. In the Dusty geometry, stars and dust are uniformly distributed in the sphere. In the Shell geometry, stars extend only to 0.3 times of the system radius and are surrounded by a shell of dust from 0.3 to 1.0 times of the system radius. Witt & Gordon (2000) found that a Shell geometry (with clumpiness) with ${\tau }_{V}\approx 1.5$ can reproduce the Calzetti Attenuation Law (Calzetti et al. 1994; Calzetti 1997).

Zoom In Zoom Out Reset image size

In addition to global geometry effects, dust clumpiness is also important. For a given mass of dust, clumpiness reduces the overall opacity and therefore the efficiency of dust in converting UV-optical radiation into IR. Popescu et al. (2000) found that while a single exponential diffuse dust disk can reproduce the optical and NIR emission, a clumpy distribution of dust that is spatially correlated with stars is needed to account for the FIR and sub-mm emission. Duval et al. (2014) found that a clumpy and dusty ISM appears more transparent to radiation (both line and continuum) compared to an equivalent homogeneous ISM of equal dust optical depth. Witt & Gordon (1996) found that a minimum of two phases of non-zero density is necessary to model the interstellar dust medium in the MW. In our clumpy models, we use a two-phase medium with a filling factor of 15% and a density ratio of 100 to implement clumpiness; Witt & Gordon (1996) found that these values reproduce the cloud mass spectrum and fractal dimension of the diffuse clouds in the MW. For comparison, we also run homogeneous models, which is equivalent to setting the density ratio to unity.

While spherical geometries are efficient for computation, they may not represent the realistic geometries of whole galaxies. For example, Baes et al. (2003) modeled spiral galaxies with double exponential disks, which more closely resembles spiral galaxies than our spherical models. Realistic geometries are suitable for studying individual objects, as they can provide spatial images and in principal reproduce detailed features better. However, they are not suitable for DirtyGrid because the radiative transfer in such geometries would take many orders of magnitude longer to compute, and the broadband SED may not contain sufficient information to constrain the additional parameters. Instead, we use simplified geometries to span the possible range in galactic environments. For the same amount of dust, the Shell geometry provides the highest efficiency in attenuation and is similar to the birth cloud of young stars, while the Cloudy geometry gives the lowest attenuation and is similar to the environments in elliptical galaxies or the bulge of spiral galaxies. In essence, we model pieces of galaxies, and in the approximate limit where the global UV–IR SEDs of galaxies are dominated by a few representative types of non-interacting pieces, we can use a combination of models to fit the global SEDs. Within a single stellar population, we do not segregate the stars by age and place young stars in more embedded dust. However, we can approximate the segregation by combining the outputs of two models of different stellar ages, with the younger model having a higher optical depth than the older model.

In the DirtyGrid, the stellar mass and the radius of the spherical geometry (the model size) are degenerate parameters. We use the same ratio of the clump size to the model size for all models. In each grid cell, the dust emission spectrum is calculated from the dust absorption spectrum. At each wavelength, the amount of energy absorbed depends on the radiation intensity. As we increase the model size (and therefore the grid cell size) by a factor of x, the radiation intensity in each cell drops by a factor of x². A corresponding increase in the stellar mass by a factor of x² restores the radiation intensity to the original level. The output SED will have the same shape but scaled up by a factor of x². Effectively, stellar mass per radius squared, essentially the mass surface density, determines the SED shape. Ivezic & Elitzur (1997) have a rigorous derivation of the scaling behavior in radiation transfer. In practice, we run all models with a 10 kpc radius and scale the output as needed.

Prior to running the first stage sampling, we ran a pilot grid with high resolution in each of the three continuous variables (stellar age, stellar mass, and amount of dust) to study the behavior of the SEDs as a function of these variables. We found that the band fluxes vary smoothly as a function of optical depth. While the precise details of the output spectrum are sensitive to the parameters, the band integrated fluxes (e.g., using the Johnson UBVRI band pass filters) are far less sensitive. This allows us to sample optical depths sparsely and then interpolate over intermediate optical depths to get band fluxes at the intermediate values. Note that the first stage sampling does not have to be precise, because the second stage sampling fills in where the sampling is insufficient in the first stage.

Figure 6 shows that the optical V band fluxes vary smoothly over two dex of optical depth (from ${\tau }_{v}$ = 0.1 to 10.0). To quantitatively measure the error interpolation introduces, we ran the model with 25 different optical depths and selected stellar ages and geometries with all other parameters fixed (solar metallicity, instantaneous star formation, and clumpy MW dust). The 25 different values are logarithmically evenly spaced between ${\tau }_{v}$ = 0.1 and 10. Then we take 1 model from every 4 models, resulting in 7 (including the ones on the edge) models sparsely sampling the optical depth. From these 7 models we calculate the band fluxes, and use linear interpolation (in log space) to predict the band fluxes at the intermediate optical depths, and finally compare the interpolated fluxes and actual fluxes. For the limited parameter space we looked at, they generally agree within 1%. The behavior for stellar mass is also smooth, and we choose to use only 15 samples to span over 7 dex of stellar mass values, including the edge points.

Zoom In Zoom Out Reset image size

The behavior of stellar ages, however, requires many more samples to resolve. For the purpose of the first stage sampling in stellar age, we study the band integrated fluxes of the unreddened spectra of stellar populations, which do not require radiative transfer. In Figure 7, we plot the V band fluxes against the full range of stellar ages, sampled at 47 points. Between 10 Myr and 13 Gyr, the age samples are roughly evenly spaced in log space (12 per dex), but below 10 Myr we included all integer multiples of Myr because PEGASE does not output results at fractional Myr ages. Having a constant amount of young stars and a steadily increasing amount of old stars, the constant star formation models give a much smoother curve. On the other hand, instantaneous star formation models exhibit relatively sharp changes in V band fluxes as they evolve. This is due to evolved stars entering rapidly changing phases or the death of luminous stars. In order to achieve <10% accuracy in band fluxes, we find that we need 14 samples of ages for constant star formation models, but the number for instantaneous star formation models is much higher and depends on the metallicity.

Zoom In Zoom Out Reset image size

To pick the minimum number of samples needed to reach 10% error due to age interpolation in burst models, we did an exhaustive search for the optimal set for each metallicity. From 10 Myr to 10 Gyr, we divided the range of age into 6 equal segments (in log space). For example, the first segment goes from 10 to 32 Myr. In each segment, we have five possible sample locations that are uniformly distributed. At each location, the sample is either present or absent. For each of the ${2}^{5}=32$ combinations, we compute the band integrated fluxes and interpolate the fluxes at the missing samples. We compare the interpolated fluxes with the correct fluxes to get the fractional error. With a error target of <10%, we pick the minimum number of samples needed in each segment and arrive at a list of selected ages. We run the age selection algorithm for different star formation types and metallicity separately. Figure 8 shows our age selection for models with instantaneous star formation from 10 Myr to 10 Gyr. The number of age samples for burst models ranges from 22 to 28.

Zoom In Zoom Out Reset image size

We improve upon the first stage sampling by filling in models at locations with insufficient interpolation accuracy. Because of possible nonlinearity in the SED landscape, we cannot determine the error in the interpolated fluxes at locations of the parameter space that we have not sampled. Even in rapidly changing areas, the interpolated fluxes could be accurate if the changes are smooth. As an estimate of the interpolation accuracy, we again compared the actual band fluxes of each model with the interpolated band fluxes using the nearby points. Above a threshold of error (15%), we fill in more samples around such models. To choose the second stage samples, we apply the following procedure for each of the three continuous dimensions (denoted as x):

• 1.  
  Loop over each samples in the grid, except the ones with the maximum or minimum allowed values of x.
• 2.  
  While ignoring the current sample, interpolate the broadband fluxes between the model that has a higher value of x and the model that has a lower value of x.
• 3.  
  Compare the broadband fluxes of the current sample with the interpolated fluxes.
• 4.  
  If the error in any of the bands is greater than 15%, include the two neighboring grid points of the current sample (in the x dimension) in the second stage sampling.

We calculate the error of the interpolated fluxes, assuming that the actual model has the true values. The procedure identified 225K locations to fill in second stage samples. In the final DirtyGrid, we interpolate over the parameter space using the combined set of first and second stage samples.