Stellar Population Inference with Prospector

In Prospector, the generation of a basic spectrum and photometry is handled by stellar population sources. These source objects take a set of parameters and compute a mass-normalized rest-frame spectrum.

In this paper, we will model both simple and composite stellar populations, describing star clusters and galaxies, respectively. The spectra for these sources in Prospector are derived from the Flexible Stellar Population Synthesis (FSPS; Conroy et al. 2009, 2010) code. The FSPS code is accessed through the python-FSPS ⁹ (Foreman-Mackey et al. 2014) bindings. Given a set of stellar population parameters, python-FSPS is used to calculate the emergent galaxy spectrum by combining simple stellar populations (SSPs) with the SFH and absorption and emission by other galaxy components. If necessary, SSPs may be generated by FSPS during this computation; otherwise, the SSPs are cached by python-FSPS, and only the calculations required to produce the spectrum of the composite stellar population (as well as nebular emission, dust, and other effects) are made, significantly speeding up the calculation. This allows for the spectrum of a stellar population to be constructed efficiently on the fly; no parameter grids are necessary, except those used internally by FSPS.

The FSPS code includes detailed stellar populations with a variety of choices for stellar isochrones and stellar spectral libraries, self-consistent absorption and emission by dust, self-consistent nebular emission, and intergalactic medium (IGM) attenuation. More details about the ingredients of FSPS, especially those that have been added or modified since Conroy et al. (2010), can be found in Appendix A. Every parameter available in the FSPS code is a potential parameter of a Prospector model and can either be free or fixed. This includes parameters that affect the computation of the SSPs like the slope of the IMF, shifts in the temperatures or luminosities of AGB isochrones, and changing blue-straggler fraction, as well as more traditional composite stellar population parameters like stellar metallicity and dust attenuation (including the slope of the attenuation curve and the age dependence of the optical depth). Nebular emission parameters may also be varied. Varying parameters that affect the SSPs will trigger their regeneration, which is computationally expensive. As will be discussed in Section 2.2, for several parameters available in FSPS—including the redshift and spectral smoothing—the computations are handled by Prospector externally to FSPS.

Prospector has access to all of the parametric SFHs built into FSPS. These forms include exponentially declining SFRs, exponentially declining SFRs with a linear rise (delayed-exponential), constant, single-age bursts, and a flexible four-parameter delayed-exponential model that transitions to a linear ramp at late times described in Simha et al. (2014). In addition, Prospector includes the ability to model linear combinations of these SFHs. For example, it is possible to model SFR(t) as the sum of an exponentially declining SFH and a constant SFH, each with separate dust obscuration, to represent the stellar population of the bulge and disk, respectively. It is advised that any user-defined SFH should be tested by generating and fitting mock data to ensure that the information content of the data is well matched to the complexity of the model (Section 2.5).

Through the use of tabular SFHs in FSPS, the Prospector code can also be used to model and fit more flexible or complicated SFHs. As currently implemented in Prospector, these are the so-called "nonparametric" SFHs, high-dimensional piecewise constant (or step-function) SFHs, the fundamental parameters of which are stellar mass formed in each of N bins of lookback time of user-defined width. Additional forms (e.g., piecewise linear, with parameters giving the SFR at particular times) may be included in the future. Nonparametric SFHs are highly flexible and capable of modeling SFHs with sharp features or unusual shapes that are forbidden by usual parametric forms. It is difficult to constrain flexible SFHs with typical galaxy observations, so the chosen prior is very important (Sections 2.4.2, 3.2 Leja et al. 2019). Through careful parameter transformation, a number of different priors on the amplitudes are possible; Prospector is packaged with four families of such priors, but the user may adjust these or explore more complicated priors.

Finally, the modular nature of Prospector allows for sources other than the ones available through FSPS to be used, in principle. All that is required is a Python object that, given a parameter dictionary and a set of wavelengths and filters, can compute and return the mass-normalized rest-frame spectrum.

When modeling data, the models need to be "projected" into the data space in various ways. These effects are handled by code that operates on the spectrum returned by any of the sources described above. They include redshifting, normalization by the total formed stellar mass and distance dimming, and filter projections. If required, filter projections are handled by the pure Python code sedpy ¹⁰ (Johnson 2019), where the addition of user-defined filters is straightforward.

To model observed spectra, the normalized, redshifted spectrum must be broadened to match the observed line-spread function (LSF) including physical and instrumental broadening. Prospector includes code for fast spectral broadening in velocity or wavelength space using fast Fourier transforms (FFTs). The broadening kernels are currently limited to Gaussians, but more complex kernels may be added in the future. Broadening with FFTs uses the analytic Fourier transforms of the broadening kernels (e.g., Cappellari 2017). A detailed discussion of the Prospector treatment of broadening due to the line-of-sight velocity distribution, the instrumental resolution, and considering the library resolution is given in Appendix C. At the time of writing, the highest resolution models available in FSPS and hence in Prospector are based on the empirical library of stellar spectra MILES (Sánchez-Blázquez et al. 2006), with FWHM ≈ 2.5 Å from 3525 to 7400 Å. Modeling even moderate-resolution spectroscopic observations of rest-frame wavelengths outside of this region, or high-resolution rest-frame optical data, will require new stellar libraries. Alternatively, the data may be smoothed to a lower resolution, but this induces correlated noise that would need to be accounted for in the likelihood (see Appendix D). Finally, because the kinematics of the gas and stars may be different, the nebular lines may need a different broadening and redshift from the stellar continuum. Prospector has this functionality: see Appendix E for details.

Adjustments to the wavelength calibration can be modeled using polynomials, following Koleva et al. (2009). Additional modeling of sky emission and transmission is under development.

The treatment of spectrophotometric calibration and the uncertainty thereof is of critical importance when modeling spectra, and especially when combining spectroscopic and photometric information (Sections 2.4, 3.1.2, Appendix D). It is common when modeling spectroscopy to remove the sensitivity of the inferred physical parameters to the continuum shape, due to the difficulty of accurate spectrophotometry (e.g., Koleva et al. 2009; Cappellari 2017). Another way to think of this is that the information in the continuum shape of the spectroscopic data should not flow to the physical parameters of interest but rather to nuisance parameters describing the uncertain spectral calibration. Here we describe approaches to the modeling of spectrophotometric calibration that are available in Prospector. We assume throughout that the spectra and the photometry probe the same stellar populations, i.e., that there are no substantial color-dependent aperture corrections. The overall scaling of the spectroscopic calibration vector incorporates any gross aperture corrections, such that normalization-sensitive parameters (e.g., stellar mass) correspond to the population probed by the integrated photometry.

The most straightforward approach is to include a parameterized function of wavelength or pixel, typically a polynomial, that multiplies the model spectrum to produce the observed spectrum. One can then fit for the parameters of this calibration function simultaneously with the physical parameters (e.g., Lee et al. 2011; Becker et al. 2015). The calibration parameters can be marginalized over when considering physical parameter estimates. It is straightforward to include priors on these parameters to incorporate external knowledge about the spectrophotometric calibration accuracy. However, it can be technically difficult to make this simultaneous fit when the number of calibration function parameters is very large (e.g., for a very high-order polynomial). In Prospector, this difficulty can be avoided by optimizing the parameters of the calibration function at each model generation. In this case, the computational cost is small, but the resulting posterior distributions do not fully incorporate the uncertainties in the spectrophotometric calibration. It is possible—for certain prior distributions and calibration functions—to analytically marginalize over the parameters at each likelihood call, but this option is not yet available in Prospector. A different approach that avoids an explicit calibration function is to model calibration effects as covariant or correlated noise, which is discussed further in Appendix D.

The parameterized calibration function available by default in Prospector is a Chebyshev polynomial of the wavelength vector of user-defined order N_p . The scale of the calibration features that can be captured with this model is ∼$({\lambda }_{\max }-{\lambda }_{\min })/{N}_{p}$. Typically, optimization of the polynomial will serve most users' needs. When no photometry is available, the parameters of the calibration vector are unconstrained except by the space of possible continuum shapes allowed by the physical model. This will lead to unconstrained stellar mass and dust attenuation, as these parameters are typically constrained by the shape and normalization of the continuum.

2.4.1. Likelihood

The basic log-likelihood calculation in Prospector is effectively a χ² calculation for both the spectral and photometric data. The two log-likelihoods are then summed. This gives the following expression for the $\mathrm{ln}$-likelihood:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L }({d}_{s},{d}_{p}\,| \,\theta ,\phi ,\alpha ,\beta ) & = & \mathrm{ln}{ \mathcal L }({d}_{s}\,| \,\theta ,\phi ,\alpha )\\ & & +\,\mathrm{ln}{ \mathcal L }({d}_{p}\,| \,\theta ,\beta ),\end{array}\end{eqnarray} \tag{ 1 }$$

where d_s is the spectroscopic data, d_p is the photometric data, the parameters θ describe the physical model and generally include any of the python-FSPS parameters, the parameters ϕ describe spectroscopic instrumental effects and calibration, the parameters α describe the spectroscopic noise model, and the parameters β describe the photometric noise model.

By default, we assume known and independent Gaussian flux uncertainties. Under this assumption, the $\mathrm{ln}$-likelihood for the spectroscopic and photometric terms is simply χ². More complex noise models are available in Prospector. These may be necessary to model a variety of instrumental artifacts and are discussed further in Appendix D.

2.4.2. Priors

2.4.3. Optimization and Sampling

Given all these ingredients in Prospector, it is the responsibility of the user to define a model that is appropriate for their data and scientific question. In Prospector this amounts to choosing a set of parameters, deciding which are free and specifying prior distributions for those, and specifying the values of the remaining fixed parameters. It may also require defining any parameter transformations to convert from the sampling parameters to the native parameters used to generate spectra and SEDs. Several examples of this process are presented in Sections 3 and 4. A significant advantage of on-the-fly model generation is that changing the model definition—the free parameters, the priors, and the parameterization—does not require rebuilding large grids of models.

When constructing a model for inference in Prospector, it is often useful to start with a simple model and to add complexity as required by the data and scientific question. A first step is to make predictions with a simple model and some fiducial or very approximate parameter values and to compare this to some representative data. Once the model seems like it has reasonable components, the next step is to perform a fit and then, crucially, to examine residuals. This can reveal missing components in the model that were not immediately apparent, such as unmodeled weak nebular emission lines, resolution mismatch, or the presence of hot dust emission. Examining the posterior PDFs can identify priors that are too restrictive or misspecified. Finally, with a basic model in hand, one can explore letting parameters that were previously fixed be free, though with informative priors.

There is a natural desire to avoid complexity in the model by limiting the number of free parameters. However, one of the principal advantages of Bayesian inference and lack of grid is the ability to add extra parameters that may affect the SED and marginalize over them if they are not scientifically interesting. Adding such parameters—the values of which are not extremely well known a priori but that might affect the SED—will cause other parameter estimates to become less certain and thus guard against overinterpreting the data. If the parameter is not well constrained by the data—that is, if it does not affect the data in a unique and significant way—then this will be reflected in the posterior PDFs, which will take the form of the prior for that parameter. Furthermore, the constraints on other parameters will only be affected if they are somehow degenerate with the new parameter; in this case, not including the additional parameter would bias the uncertainties low on these other parameters. Thus, one should not shy away from including poorly constrained parameters that are known to be physically important or to affect the SED, though in such cases, it is important to put thought into the adopted prior distribution. The only drawback to adding parameters is that MCMC sampling can become inefficient as the number of parameters increases.

When the data are only weakly informative about the parameters, as can often occur in SED fitting, the form of the priors can play a strong role in the posterior estimation. It is often useful to place strongly informative priors on some parameters, as opposed to weakly informative priors such as uniform distributions. If one has external information about the plausible values that a parameter might take—for example, that the shape of the attenuation curve is most likely to be a certain value but that there is scatter—then it is best to use that in the form of a prior. This is preferred to the typical approach of fixing the parameter to the single most likely value: that effectively constitutes an extremely strong prior and will often result in overly confident posteriors.

However, the importance of priors can also raise concerns about the interpretation of the results—did you actually learn anything from fitting your data, or are the results dominated by the prior assumptions? When this is a concern, we recommend testing the sensitivity of the results to the form of the prior. This can be done most directly by running fits with a different but still plausible prior and seeing if and how the results change. It can also be useful to generate SEDs or other derived quantities by directly sampling the chosen prior distribution and inspecting the distribution of those quantities to see what values would be preferred or allowed by the model.

Since the pioneering work of Tinsley (1980), it has been common to fit observed SEDs to galaxy models where the SFH is represented by a relatively simple functional form characterized by a few parameters, SFR(t) = f(t, θ). The commonly used exponentially declining SFH arises in closed-box evolutionary models with a constant star formation efficiency. In addition to the exponentially declining and the related delayed-exponential SFHs, recent studies have proposed a variety of functional forms thought to be more or less well matched to the SFHs of real galaxies (Gladders et al. 2013; Simha et al. 2014; Diemer et al. 2017; Carnall et al. 2018). Parametric SFH models attempt to describe the potentially complex evolution of galaxies with only a few numbers. This is very desirable for limiting the size of SED grids and making any Bayesian inference with MCMC more tractable.

Here we present a fit to mock data with a parametric SFH model. This provides a comparison point for the extensive literature using parametric SED models, and the simplicity of the model allows us to focus on different aspects of the Prospector fitting. For all of the fits presented in this section, we use a delayed-exponential SFH, where the SFR ψ as a function of lookback time t_ℓ is given by

$$\begin{eqnarray}&&\psi ({t}_{{\ell }})\propto \displaystyle \frac{{t}_{\mathrm{age}}-{t}_{{\ell }}}{\tau }{e}^{-({t}_{\mathrm{age}}-{t}_{{\ell }})/\tau },\quad 0\lt {t}_{{\ell }}\lt {t}_{\mathrm{age}}\end{eqnarray} \tag{ 2 }$$

with parameters t_age and τ. While a number of treatments of dust attenuation are available through FSPS (see Appendix A), we model the dust attenuation with a simple dust screen affecting stars of all ages equally and specify a fixed wavelength dependence of the optical depth based on Kriek & Conroy (2013). The normalization of this attenuation curve—the optical depth at 5500 Å—is left as a free parameter. We also include as a free parameter a single metallicity for all stars. For simplicity, we tie the metallicity of the nebular emission model (Appendix A) to this stellar metallicity, but this is not required in general. The ionization parameter U of the nebular emission is left as a free parameter. Finally, the overall normalization of the SFH is left as a free parameter, corresponding to the total mass of stars formed or the integral of the SFH. We include dust emission in the model, based on the three-parameter models of Draine & Li (2007; see Appendix A for details), but we do not fit for these parameters, as we do not include data at infrared wavelengths that are sensitive to these parameters. Unless otherwise stated, values of the fixed parameters used to generate the mock data are the same as the parameters used to fit the mock data.

The modeling of redshift is described and treated separately in the fits below. However, the redshift of the mock is set to 0.1, and the prior distribution on t_age does not allow solutions older than the age of the universe at this redshift. The values of the other mock parameters (and the priors used when modeling the SED) are given in Table 1.

Table 1. Summary of Parameters and Priors Used in Mock Demonstrations

Parameter	Mock Value	Prior Functions
Stellar Mass a [M⊙]	1010	LogUniform(107, 1013)
Age tage [Gyr]	12	Uniform(0.1, tuniv b )
τSF [Gyr]	4	LogUniform(0.01, 100)
z	0.1	${ \mathcal N }(0.1$, 0.001)
$\mathrm{log}({Z}_{\star }/{Z}_{\odot })$	−0.3	Uniform(−2.0, 0.2)
τ5500	0.3	Uniform(0, 4)
$\mathrm{log}\,U$	−2	Uniform(−4, −1)
σv b [km s−1]	200	Uniform(30, 400)

Notes. Uniform(x, y) indicates a uniform prior distribution from x to y, while ${ \mathcal N }(\mu ,\sigma )$ indicates a normal distribution for the prior, and LogUniform(x, y) indicates a distribution that is uniform in the log of the parameter in the parameter range (x, y)

^a This is the integral of the SFH. ^b The maximum age is set by the age of the universe at the model redshift in a WMAP9 cosmology. ^c Velocity dispersion is only included as a model parameter when fitting spectra

Download table as:  ASCIITypeset image

3.1.1. Inference from Photometry

In this section, we demonstrate inference of the parameters of the simple model described above from only the mock photometric data. The mock photometry is generated using the same model used for the inference, with parameters given in Table 1. We generate photometry in the filters of several large-area sky surveys, including GALEX FUV and NUV, SDSS ugriz, and the 2MASS JHK_s bands. We assume uncorrelated uncertainties of 5% in the photometry in each band and add a realization of the noise to the mock photometry.

For this demonstration, we assume, when fitting the photometry, that the redshift is constrained to the true value to within 0.001. That is, we include the redshift as a free parameter, but with a very narrow prior distribution as might be expected from low signal-to-noise ratio (S/N) grism spectra, yielding a model with seven free parameters. We then use Prospector to explore the posterior probability distribution of the free parameters for this model. For the inference, we use nested sampling, with default parameters. In Figure 1, we show the mock spectrum and the mock photometry, the latter perturbed by a realization of the noise. We also show the inferred values for the true photometry as shaded boxes with heights indicating the 16th–84th percentile range of the posterior probability. Uncertainty-normalized residuals for the highest probability posterior sample are also shown. Examining such residuals is a key step in evaluating the performance of any model and identifying any deficiencies that require additional model components or flexibility to adequately describe the data.

Zoom In Zoom Out Reset image size

In Figure 1, we also show projections of the posterior probability distribution for several of the free parameters (the redshift and the nebular ionization parameter are not constrained by the data, and the posterior distributions are indistinguishable from the prior). The input mock parameter values are also shown in black, and the adopted prior probabilities are indicated in the one-dimensional projections. The posterior PDFs for a number of the parameters are quite broad. Notably, there is a strong degeneracy between t_age and τ_SF such that while neither is well constrained, their ratio is reasonably well constrained. The ratio t_age/τ is closely related to the specific star formation rate. The well-known dust-age–metallicity degeneracy is mapped in complete detail.

It is worth pointing out that the peak of the marginalized posterior probability distributions are not generally centered on the true values. This is a common and expected occurrence, caused in part by the perturbation of the mock fluxes by a realization of the noise and in part by the complicated shape of the posterior in multidimensional parameter space combined with the adopted priors. This highlights the potential drawbacks of simply using the maximum a posteriori parameters.

3.1.2. Incorporating Information from Spectroscopy

We now turn to the inference of mock parameters from mock spectroscopic data, both alone and in combination with photometric data. We generate a mock spectrum from the model, smooth the mock spectroscopy by an additional velocity dispersion of σ_v = 200 km s⁻¹, and assume an instrumental dispersion of 2 Å per pixel from rest-frame 3800–7100 Å (corresponding to the region of the spectrum covered by the high-resolution empirical MILES stellar library). The applied physical velocity smoothing dominates over the intrinsic resolution of the MILES stellar library, but for clarity and simplicity, we also assume that the instrument has the same LSF as the MILES library redshifted to z = 0.1 (∼2.75 Å FWHM). See Appendix C for approaches to data where the instrumental or library broadening is a significant contributor to the final LSF of either the data or model. We assume 10% uncertainties in the flux of each pixel, and a realization of this noise is added to the data. To mimic realistic observations where the spectrophotometric calibration may not be known, and to highlight that we are not using information from the continuum shape, we divide the mock spectrum by a sixth-order polynomial fitted to regions of the spectrum free from strong absorption or emission lines.

To model the continuum-normalized mock spectrum, we include all free parameters used in the fit to the photometry (Table 1), and we also include the galaxy velocity dispersion as a free parameter. To account for any uncertainty in the spectroscopic calibration, we use the method of multiplying the model spectrum by the maximum-likelihood polynomial that describes the ratio between the spectroscopic data and the model before computing the likelihood of the data. We use a 12th order polynomial that is determined at each likelihood call. By including this polynomial multiplication, we do not use information from the spectral continuum shape; that information is supplied by the photometry if it is also being fit.

In Figure 2, we compare the posterior PDFs for several model parameters inferred from the continuum-normalized spectroscopy alone (in maroon contours) to those inferred from the photometry alone (blue contours, previous section). Given the strong degeneracy between t_age and τ when fitting to photometry alone, we instead show the posterior PDF for the ratio t_age/τ. In this comparison, we see that the posterior PDF for the dust attenuation is much broader when inferred from continuum-normalized spectroscopy than when inferred from the photometry, and is essentially the same as the prior. This is expected because all of the continuum shape information in the data is absorbed by the multiplicative polynomial. In other words, the dust attenuation is completely degenerate with the spectrophotometric calibration, which we have assumed to be completely unknown. The stellar mass is similarly unconstrained. The metallicity is tightly constrained by the spectroscopy; this is not surprising, as in addition to the metallicity information provided by absorption lines in the stellar continuum, we have tied the nebular emission model to the stellar metallicity and include the emission lines in the fitted data.

Zoom In Zoom Out Reset image size

Figure 2 shows that broadband UV through NIR photometry and continuum-normalized spectroscopy carry different information about the model parameters. It is natural then to combine the constraints from each. We do this by simultaneously modeling and fitting the photometry and spectroscopy using the combined likelihood described in Section 2.4. ¹² Because we have included a polynomial calibration or normalization vector when modeling the spectra, the information about the continuum shape and any associated parameter is entirely supplied by the photometric data. If we did not include this uncertainty in the spectroscopic calibration, then the much larger number of spectroscopic points with comparable S/N to the photometric data would overwhelm any information in the photometry, unless two models with equally likely spectra somehow had significantly different photometry. This latter situation would likely occur only if the SFH and dust model were very flexible, or the spectra were of much lower S/N than the photometry.

In Figure 2, we show the posterior probability distribution functions of the selected parameters obtained from a simultaneous fit to the photometry and the spectroscopy (gold contours). For all parameters, these are narrower than the posterior PDFs obtained from either photometry or spectroscopy alone. Interestingly, while the continuum-normalized spectroscopy places no constraint on the dust attenuation by itself, its combination with photometry results in even tighter constraints on the attenuation than from the photometry alone. This is because the precision on other parameters, such as $\mathrm{log}(Z/{Z}_{\odot })$, which are partially degenerate with dust attenuation in the photometry alone, is increased by the inclusion of the spectroscopic data.

The assumption of a particular functional form for the SFH of a galaxy amounts to an extremely strong prior on that SFH. However, it is unclear how well real galaxies are represented by these parametric models (e.g., Gallagher et al. 1984). Deviations from the assumed forms affect the accuracy of derived parameters (e.g., Lee et al. 2009; Pforr et al. 2012; Wilkins et al. 2012; Pacifici et al. 2015; Lower et al. 2020), as well as estimates of the uncertainty on the SFH (Carnall et al. 2019a; Leja et al. 2019).

One way to address the complex SFHs in real galaxies is to increase the flexibility of the model SFHs. Flexibility has been introduced previously by adding additional components (i.e., recent bursts), by choosing more flexible forms for the parametric SFH (e.g., Salim et al. 2007; Carnall et al. 2018), or using a library of SFHs based on semianalytic or hydrodynamic simulations (e.g., da Cunha et al. 2008; Pacifici et al. 2013; Iyer et al. 2019). These approaches can require large library sizes to densely sample the prior SFH parameter space and also incorporate the simulation as an SFH prior.

Another possibility is to use piecewise constant or δ-function SFHs where the mass in each temporal bin becomes a parameter of the model. These so-called nonparametric—actually many-parameter—models have become the standard in SFH inference from (semi)resolved stellar color–magnitude diagrams (e.g., Dolphin 1997; Harris & Zaritsky 2001; Weisz et al. 2011; Cook et al. 2019) and have been used in maximum-likelihood "inversion" codes for unresolved spectra (e.g., Cappellari & Emsellem 2004; Cid Fernandes et al. 2005; Ocvirk et al. 2006; Tojeiro et al. 2007; Koleva et al. 2009). In principle, any basis function can be used (e.g., Iyer & Gawiser 2017) though the piecewise constant basis has the advantage of being readily interpreted.

Prospector supports nonparametric models using a piecewise constant basis (Section 2.1). In contrast to maximum-likelihood techniques, Prospector allows for the exploration of the full posterior PDF for the component amplitudes. Inference with this kind of nonparametric SFH was explored in detail by Leja et al. (2019), where a number of different priors on the amplitudes of each bin were considered. A prior based on the ratio of the amplitudes in adjacent bins has some desirable qualities, and we adopt such a "continuity" prior in this section. In this model, the sampled parameters are the log of the ratio of the SFRs in adjacent time bins and the log of the total stellar mass. The mass M_i formed in each time bin—the native parameters used to generate the spectrum—depend on these sampled parameters through a transform:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\,{r}_{j} & = & \mathrm{log}({s}_{j}/{s}_{j+1})\\ {s}_{0} & = & 1\\ {s}_{i} & = & \displaystyle \frac{{s}_{0}}{{\prod }_{j=0}^{j\lt i}{r}_{j}^{-1}}\\ {m}_{i} & = & {s}_{i}\,{\rm{\Delta }}{t}_{i}\\ {M}_{i} & = & \,{m}_{i}\,\displaystyle \frac{{10}^{\mathrm{logmass}}}{{\displaystyle \sum }_{i}{m}_{i}},\end{array}\end{eqnarray} \tag{ 3 }$$

where Δt_i is the width of the ith bin of the lookback time, in years. We place a Student's t-distribution prior on each of the $\mathrm{log}\,{r}_{j}$ parameters, centered on 0 (i.e., constant SFR); the Student's t distribution is similar to a Gaussian prior but with more probability in the tails. The effect of this prior is that in the absence of a constraint from the data the bins default to a constant SFR.

To demonstrate the capabilities of Prospector to infer both parametric and nonparametric SFHs, we construct mock data based on several galaxies from the Illustris cosmological hydrodynamic simulation of galaxy formation (Vogelsberger et al. 2014; Diemer et al. 2017). The galaxies were selected to have sharp variations in both the recent and the extended SFH. This selection provides the greatest challenge for the inference machinery and are also disfavored by the chosen prior, ensuring that any agreement is truly driven by the data. These SFHs have a temporal resolution of 0.136 Gyr (i.e., t_univ/100), and we linearly interpolate between the provided SFRs to build the mock spectrum. We assume a single stellar metallicity and a simple foreground dust screen. The mock spectra are generated for a redshift z = 0.1 and smoothed by a velocity dispersion of σ_v = 100 km s⁻¹. As for the mocks in Section 3.1, we assume an observed-frame instrumental LSF equivalent to the MILES spectral library at the assumed redshift. We assume an S/N of 100 for each pixel. For simplicity, we do not include nebular emission. We infer the SFH from the mock spectra with both parametric and nonparametric SFH models. For both models, we simultaneously infer the SFH, the stellar metallicity, the dust attenuation optical depth for a simple foreground screen, the stellar velocity dispersion, and the redshift. To focus the analysis on the physical information provided by ideal data we assume that the spectrophotometric calibration is perfect.

For the parametric SFH, we use the same delayed-exponential model as in Section 3.1. For the nonparametric SFH model, we adopt 14 bins approximately evenly spaced in $\mathrm{log}\,{t}_{{\ell }}$. This choice of the number of bins and their exact widths is driven by experimentation, the difficulty in sampling a very high-dimensional parameter space when using an extremely large number of bins, and considering the degree of difference between the spectra in adjacent bins.

In Figure 3, we show the input Illustris SFHs and the posterior PDFs for the SFHs inferred with both a parametric and nonparametric model. In these examples, we find that the parametric model will tend to reproduce the shape of the very recent SFH—which often accounts for the stars that dominate the optical light—even at the expense of reproducing the shape of the SFH at earlier epochs if that is not allowed by prior. In contrast, the flexibility of the nonparametric SFH model is better able to recover the shape of the recent and older SFH. However, the influence of the adopted prior on SFR ratios remains evident, though at a lower level—the recovered nonparametric SFH will tend toward approximately constant SFR when the data are not constraining, and obviously, changes in the SFR within a given temporal bin are unrecoverable. We also find that the constraints on the SFH when using a parametric model are too strong when compared to the nonparametric SFH and less able to encompass the true range of SFHs consistent with the data (including the true input mock SFH). Finally, these differences in the inferred SFHs also affect inferred stellar masses, mass-weighted stellar ages, and star formation rates. For more discussion on these points, we refer the reader to Carnall et al. (2019a), Leja et al. (2019), and Lower et al. (2020).

Zoom In Zoom Out Reset image size

In this section, we explore how the inference of a complex SED model is affected by increasing the number of photometric bands used to constrain it. The inferences we have shown to this point have been of relatively simple models, with few parameters. However, more complex models may be required to adequately describe the increasing quality of observations. Furthermore, making robust inferences from even a small number of low-quality data—and accurately representing the uncertainties on the underlying parameters—requires modeling all the physical properties that might affect that data. For example, there is evidence that the effective attenuation curve in galaxies is not fixed, so any quantities inferred assuming a fixed attenuation curve will tend to have underestimated uncertainties.

One reason to make a number of simplifying assumptions about the relevant physical ingredients is that generating and storing grids for high-dimensional parameter spaces is computationally challenging. One of the most important features of on-the-fly model generation with Monte Carlo sampling of the posterior is the ability to infer a larger number of physical parameters than can be considered in grid-based approaches, because models are only generated in regions of high probability as needed, instead of over the whole parameter space.

The model we explore is described by a nonparametric SFH with six temporal bins, a single stellar metallicity for all stars, a dust screen attenuation curve affecting all stars and described by a free power-law index and a V-band optical depth, an extra attenuation toward young stars described by a power-law attenuation curve with fixed power-law index, nebular emission described by the gas-phase metallicity and ionization parameter, AGN dust torus emission described by two parameters (Nenkova et al. 2008a; Leja et al. 2018, A) and dust emission described by the three parameters of the Draine & Li (2007) model. In total, this model has 16 free parameters, listed in Table 2. We assume a known redshift, z = 0.1. We adopt informative priors on many of these parameters, particularly the dust attenuation parameters and the AGN parameters. These prior distributions favor typical dust attenuation curves and disfavor large AGN contributions to the mid-IR.

Table 2. Summary of Parameters and Priors Used in the Complex Model

Parameter	Mock Value	Prior Functions
SFH
$\mathrm{log}{M}_{\star }$ a [M⊙]	10	Uniform(7, 13)
$\mathrm{log}{r}_{i}$ b	0	StudentT(0, 0.3, 2)
$\mathrm{log}({Z}_{\star }/{Z}_{\odot })$	−0.3	Uniform( − 2.0, 0.2)
Dust Attenuation
τ5500,diffuse	0.5	${{ \mathcal N }}_{c}(0.3,1.0,0,4)$
τyoung/τdiffuse	1	${{ \mathcal N }}_{c}(1,0.3,0,1.5)$
Γdust	0	${{ \mathcal N }}_{c}(0,0.5,-1,0.4)$
Nebular Emission
$\mathrm{log}\,{U}_{\mathrm{neb}}$	−2	Uniform( − 4, − 1)
Dust Emission c
${U}_{\min ,\mathrm{dust}}$	2	Uniform(0.1, 25)
QPAH	1	Uniform(0.5, 7)
γdust	10−3	LogUniform(10−3, 0.5)
AGN d
fAGN	0.05	LogUniform(10−5, 3)
τAGN	20	LogUniform(5, 150)

Notes. Uniform(x, y) indicates a uniform distribution from x to y, while ${{ \mathcal N }}_{c}(\mu ,\sigma ,x,y)$ indicates a normal distribution, clipped to the range (x, y), and LogUniform(x, y) indicates distribution that is uniform in the log of the parameter in the parameter range (x, y).

^a This is the integral of the SFH. ^b r_i is the ratio of the SFR in temporal bin i to that in the adjacent bin. There are five such parameters that describe the SFH. ^c Parameters for the Draine & Li (2007) dust emission model; see Appendix A. ^d Parameters for the Nenkova et al. (2008a) AGN torus emission model; see Appendix A.

Download table as:  ASCIITypeset image

We generate broadband photometry in standard UV-FIR filters for a mock SED generated for a particular value of the model parameters. We then infer the model parameters using increasingly large subsets of the photometry, beginning with a single SDSS r-band data point. The other filter sets add SDSS optical bands, then 2MASS near-IR bands, then GALEX UV bands, then WISE mid-IR bands, and finally Herschel far-IR bands.

In Figure 4, we show the resulting posterior distributions for parameters of the model and for the SED. To summarize the inferred SFHs, we derive and present the recent SFR and the mass-weighted stellar age of the population. We see that for a small number of photometric bands, the posteriors of nearly every parameter are determined by the prior distribution. Nevertheless, even with a single photometric point, we are able to obtain a (poor) constraint on the total stellar mass. For this particular mock SED, the posterior PDF is shifted to a slightly larger stellar mass than the true value. This is largely because the true dust attenuation—which is entirely unconstrained except by the prior—is at the lower end of the prior distribution. Therefore, the prior (and posterior) distribution for the SED is redder and has a larger mass to light ratio than the true SED. This exercise demonstrates how through the use of informative priors and Bayesian reasoning even a single data point can be informative in the context of a very high-dimensional model.

As the number of photometric bands increases, the posterior constraints on many parameters become stronger and, more importantly, distinct from the prior distribution. Adding even a single additional band leads to a much stronger constraint on the dust attenuation and hence mass-to-light ratio than a single band (e.g., Bell & de Jong 2001). With the full SED, nearly every parameter has a posterior PDF different from the prior.

Globular clusters have long been used as a strong test of stellar population models (e.g., Bruzual & Charlot 2003; Koleva et al. 2008). They are generally single-age, single-metallicity populations, making them real-world examples of SSPs. Many clusters in the Milky Way have high-quality ages based on the main-sequence turnoff and metallicities determined from high-resolution spectroscopy of individual stars. However, globular clusters present unique challenges—at low metallicities, they have very prominent blue horizontal branches (BHBs), which may be a consequence of the unusual light-element abundance variations that have been detected in nearly all clusters (Carretta et al. 2009). If these unusual features are confined to the cluster population, then they may be less useful as tests of stellar population models.

Here we fit the integrated spectra and photometry of a sample of Milky Way globular clusters with Prospector. The sample is drawn from the spectral atlas of Schiavon et al. (2005), who conducted drift-scan long-slit spectroscopy of 40 globular clusters. These spectra span the wavelength range 3400–6400 Å, with a spectral resolution that is wavelength dependent but averages 3.2 Å FWHM. The nominal S/N is often greater than 100. We take UBVRI photometry from Harris (1996), who provides the total (integrated) magnitudes of clusters. We assign uncertainties of 10% to the photometric fluxes.

To model this data, we treat each cluster as a single-age stellar population. The free parameters of this population are its mass, age, and stellar metallicity $\mathrm{log}(Z/{Z}_{\odot })$ (assuming scaled solar abundances). We fix the distance to the value given by Harris (1996) because for integrated SEDs at noncosmological distances, the distance and mass are perfectly degenerate. We broaden the model library spectra to match the instrumental resolution as reported by Schiavon et al. (2005) and fit for an additional stellar velocity dispersion σ_v as well as a systemic redshift. Extinction is modeled with the attenuation curve of Cardelli et al. (1989), with the normalization of this curve left as a free parameter. We also multiply the model spectrum by a maximum-likelihood 15th order Chebyshev polynomial at each likelihood call to account for imperfections in the spectrophotometric calibration.

In Figure 5, we show an example of the spectroscopic and photometric data, the residual of the posterior prediction for the SED from the data, the posterior prediction for the calibration vector, and the inferred posterior probability distribution functions for the key physical parameters. The photometry is typically well fit by the models, and the spectra can be fit simultaneously when accounting for offsets from the photometry that vary smoothly with wavelength by ∼10% and inflation of the nominal noise estimates by factors of ∼2.

Zoom In Zoom Out Reset image size

In Figure 6, we show a comparison of the parameters that we infer for many globular clusters, compared to literature values as compiled by Koleva et al. (2008) and to values inferred from the same spectra using the UlySS code (Koleva et al. 2009). The literature values are based on stellar CMD fitting and metallicity measurements of individual stars. We find that the inferred metallicities from combined spectroscopy and photometry with Prospector are in good agreement with literature values, even in the presence of unmodeled α-element enhancements. The median offset is 0.04 dex and the rms scatter is ∼0.2 dex, including outliers at very low metallicity where the stellar libraries in our models are incomplete. Good agreement is also achieved between the inferred and literature values for the reddening (A_V ).

Zoom In Zoom Out Reset image size

Previous work has shown the integrated-light spectral fitting of globular cluster data returns generally unreliable ages (Koleva et al. 2008; Conroy et al. 2018). We find similar results with Prospector. The primary issue is the impact of hot BHB stars which, because they are not included in the modeling, will tend to result in artificially young inferred ages (see Koleva et al. 2009 for a discussion of this effect). Indeed, in Figure 6, one sees that the clusters with the most significant BHB populations (as inferred from the HBR parameter) have the youngest ages from Prospector. There are also some cases where integrated-light measurements (both from Prospector and Koleva et al. 2008) result in maximally old ages. The origin of this shortcoming is unclear but could be due to other model systematics such as unmodeled α enhancement or other limitations in the isochrones. We note that differences between 10 and 14 Gyr are quite subtle in integrated light.

Overall, we find that Prospector performs as well as other stellar population fitting programs when applied to globular clusters.

Inferring galaxy redshifts from photometric data is one of the core applications of stellar population models. Here, we demonstrate the use of Prospector for this task, using the high-redshift galaxy candidate GNz-11. This object was selected as a high-redshift galaxy candidate on the basis of photometry with Hubble Space Telescope and Spitzer in the GOODS-North field, and follow-up grism spectroscopy with HST confirms a redshift z = 11.09 (Oesch et al. 2016).

The model used to fit the data includes the galaxy redshift as a free parameter. The other free parameters are the total stellar mass formed, the stellar and nebular metallicity, the nebular ionization parameter, and three parameters describing the dust attenuation. Due to the potential for the photometry to probe the rest-frame λ < 1216 Å spectrum at high redshift, we include a redshift-dependent IGM attenuation following Madau (1995). This includes a free parameter that scales the total IGM opacity, intended to account for line-of-sight variations in the total opacity (see Appendix A for details). The SFH is modeled with five temporal bins, with the continuity prior discussed in Section 3.2. The exact width and extent of these bins are dependent on the redshift of a given model, bounded such that no stars are older than the age of the universe. The priors are similar to those given in Table 2. For the multiplicative scaling of the IGM attenuation curve, we adopt a Gaussian prior distribution centered on 1, with a dispersion of 0.3. The ideal redshift prior would be the product of the redshift-dependent differential cosmological volume element and the integral of the redshift-dependent luminosity function above some limiting magnitude. However, for clarity and due to the lack of knowledge about the luminosity function at these epochs, we adopt a uniform prior over the range 1 < z < 13. There are 13 free parameters in the model.

The results of fitting this model to the photometric data are shown in Figure 7. The photometric data are well fit, and the photometric residuals from the sample with the highest posterior probability are consistent with the stated errors. The posterior probability distribution for the redshift shows a broad peak at $z={10.8}_{-0.4}^{+0.4}$, and the highest posterior probability sample has z = 10.35, though many samples with nearly equal probability are distributed from z ∼ 10.3–11.4. This is in excellent agreement with the redshift of z = 11.09 ± 0.1 inferred from HST grism data. We do not find significant additional or secondary peaks in the posterior distribution of redshift using this model, though a much weaker peak is present at z ∼ 2. Posterior distributions for the SFH and several additional model parameters suggest a rising SFH, a formed stellar mass of ∼10^9.6 M_⊙, low dust attenuation, a weak constraint on the stellar metallicity, and no constraint on the IGM opacity scaling factor beyond the imposed prior distribution.

Zoom In Zoom Out Reset image size

We fit a post-starburst galaxy from the SDSS IV (Gunn et al. 2006; Blanton et al. 2017), chosen to highlight Prospector's capability to perform a simultaneous fit to spectra and photometry. A post-starburst galaxy is characterized by strong hydrogen absorption lines, indicating it has decreased its star formation activity sharply in the past ∼500 Myr. We select a post-starburst galaxy to demonstrate both the flexible SFH model and the flexible emission-line modeling available within Prospector.

We fit SDSS J133751.29+354755.7 from the post-starburst catalog of Goto (2007), specifically targeted for its high-S/N spectrum and substantive emission-line infilling. We fit the ugriz photometry and optical spectrum (Smee et al. 2013) from SDSS Data Release 16 (Ahumada et al. 2020). A noise floor of 5% is enforced on the photometry consistent with uncertainties in the underlying stellar models.

A 19 parameter physical model is used. The total mass, stellar metallicity, and velocity dispersion are set free. The redshift is allowed to vary within 0.01 around the SDSS catalog value of z = 0.073. The dust model is the two-parameter Charlot & Fall (2000) model, plus a flexible attenuation curve (Noll et al. 2009) with the UV bump tied to the slope of the curve (Kriek & Conroy 2013). An eight-bin nonparametric SFH is used, with the time bins spaced logarithmically and a continuity prior applied (Leja et al. 2019). Nebular emission is implemented using the Byler et al. (2017) grid, with metallicity, ionization parameter, and nebular velocity dispersion set free. Nebular line marginalization is also turned on, as described in Appendix E; in brief, the amplitude of each emission line is marginalized at each model call to fit the data. The prior on the emission-line fluxes is taken as a Gaussian, with both the center and the standard deviation taken to be the expected line luminosities from the Byler et al. (2017) physical model. A polynomial vector is used to normalize out the continuum during each model call, following Section 2.3. Finally, the pixel outlier model and jitter term from Appendix D are both adopted.

The results of this fit are shown in Figure 8. Prospector is able to fit the data with high precision; the best-fit model produces excellent χ² results for both the photometry and the spectrum. Though they are included in the fit, neither the spectral noise inflation nor the pixel outlier model is assigned any significant weight. As expected, a significant post-starburst event is detected within the previous ∼300 Myr, suggesting a factor of 20 decrease in the ongoing rate of star formation. The mass-weighted age is 1.11 Gyr, indicating a relatively short formation phase for this galaxy, and the dust attenuation is significant at A_V = 0.74.

Zoom In Zoom Out Reset image size

The spectral fit includes a zoom-in on the Hα and [N ii] doublet. The [N ii]/Hα ratio is very high, as is typical of post-starburst objects (Yan et al. 2006; Yang et al. 2006; French et al. 2015). This is consistent with SDSS J133751.29+354755.7 being a LINER; in such objects, the nebular emission lines are thought to be partially or fully powered by evolved stars, AGNs, or shocks (e.g., Kewley et al. 2006). This is in contrast to the built-in assumption in FSPS that emission lines are powered by young stars only. This difference has the potential to cause two distinct issues: first, additional young stars would have to be added in order to fully power this emission, and second, the observed [N ii]/Hα ratio is challenging and likely impossible to produce with star formation alone. Indeed, the prediction from the CLOUDY grid is indicated in blue and is a poor fit to the data. The nebular marginalization technique, indicated in red, solves both of these issues by allowing flexible extensions from the standard CLOUDY grid. It is suggested to allow nebular marginalization in spectral fits whenever there is a chance that the emission lines are powered by sources that do not lie on the standard CLOUDY grid.

The development of Prospector began in 2014, and it has undergone continual development, refinement, and refactoring since then. This development is usually guided by areas of practical focus, as the diversity and complexity of questions and data to which SED fitting is applied make general solutions nearly impossible. Because of this, Prospector has tended to become more modular with time, in order to maintain flexibility while adding new capabilities. A benefit of this modularity is that changing or extending the code to tackle new and more detailed data is relatively straightforward. In a sense, Prospector is a toolbox for SED fitting. The tradeoff for flexibility is often computational speed; some tasks that for particular problems could be optimized with a different code structure are instead less efficiently implemented in a modular structure.

Priors are an essential ingredient in SED fitting. The SEDs of single-age stellar populations change subtly with both age and metallicity, making SFH inference from an integrated SED particularly difficult. The SEDs clearly contain information; the question is how to extract the available information in the face of substantial degeneracies. The solution is to impose priors on the SFH, either explicitly or through a particular SFH parameterization, that minimize the impact of these degeneracies where the data are not informative yet still allow the data to "speak." Nevertheless, it is important to map the remaining degeneracies and to consider the effect of the adopted prior.

Priors on dust emission parameters can also play an important role when the far-IR SED is incompletely sampled. The inferred total infrared luminosity depends sensitively on these parameters in this case, and the imposed energy balance between dust emission and dust attenuation can affect many other parameters (e.g., Appendix C of Leja et al. 2017). It is important to adopt priors that reflect prior knowledge about the distribution of these parameters, for example, derived from subsamples with more complete IR coverage. Through the Bayesian framework, this prior knowledge can be propagated through to the posterior PDFs.

While many of the uncertainties and degeneracies that can be explored through a Bayesian framework are important to consider, significant systematic uncertainties in the stellar population synthesis models remain. The modularity of FSPS grants the ability to explore some of these as model parameters, for example, the temperature distribution of horizontal branch stars or the luminosity of the TP-AGB phase. However, many other uncertainties are "baked in" to the stellar population codes: these include binary evolution, stellar rotation, ionizing fluxes and, more generally, the blue/UV spectra of low-metallicity populations, other isochrone details, and the detailed physics of nebular emission. Results that depend on the details of these calculations should be treated with caution.

As described above, Prospector is under continual development, and we expect to add new capabilities in the future. Of special emphasis is computational speed. The computational time required is set by the product of two factors: the speed at which stellar population models are constructed and the number of models required to adequately sample the posterior distribution (which increases with the dimensionality of the model). In the former category, we are exploring the use of emulators (e.g., Alsing et al. 2020) to produce models quickly, with the disadvantage that the emulator must be trained on a precomputed (and therefore largely inflexible) grid. In the latter category, we are considering new Monte Carlo techniques based on gradients. Such techniques will also allow larger model dimensionality.

While the underlying FSPS stellar population synthesis models used by Prospector include most physical components of a galaxy—stellar emission, nebular emission due to stellar photoionization, dust attenuation, dust emission, and the emission of AGN dust torii—additional components may be added. These include ISM absorption, nebular emission from shocked ISM gas, free–free emission from H ii regions, radio synchrotron emission, and X-ray through optical emission from AGNs. Additionally, Lyα radiative transfer models may prove important for modeling high-redshift galaxy spectra in detail. For some of these components, straightforward ab initio physical models exist. For others, empirical relationships may have to be adopted.

We will also pursue improved treatment of the existing components. For the stellar populations, this includes the addition of α-element enhancement to the underlying SPS models in both the isochrones and the spectra. We also plan to expand the capability to model age-dependent stellar metallicity and dust attenuation, as well as additional SFH forms. Finally, more flexible and self-consistent nebular emission models may prove useful.

Performing optimized SFH inference is a perennial challenge. The SFH (=SFR(t, Z)) is fundamentally a density; inference of densities from samples or integrals thereof is a notoriously difficult, and perhaps even ill-posed, problem. A variety of solutions to the problem have been proposed. These can all be considered various choices about the priors that are placed on the form or properties of the SFH. Ideally, we would combine a prior with a very large number of time bins, or anchor points, without consideration of the uniqueness of the spectra for each bin. Then, the data would be able to dictate which bins (if any) were well constrained while the prior would take over where the data do not provide information. However, sampling in such high-dimensional spaces is difficult, though solutions based on scalable gradient-based sampling approaches (Hamiltonian Monte Carlo) are under development. Additionally, more complex priors on the shape of the SFH—for example, corresponding to the correlation function between SFRs in different time periods (e.g., Kelson et al. 2016; Tacchella et al. 2020)—may be implemented.

Multiple spectra described by the same physical parameters but substantially different instrumental parameters (e.g., spectral resolution or calibration) cannot yet be fit simultaneously in Prospector, though this capability will be added in the future. This can occur due to multiple exposures of the same object taken under different observing conditions, or to spectra from different instruments or spectroscopic orders, or both.

In summary, as the available data increase in quality and diversity, SED models and inference techniques have had to improve as well, while growing in complexity. Significant progress has been made, leading to the ability to model and make inferences from many different kinds and combinations of galaxy observations. Efforts are underway to tackle the remaining challenges. These efforts are necessary to fully exploit the information content of rich data sets that will be provided by future facilities, such as the James Webb Space Telescope.

First, we write the likelihood of the data, conditioned on the parameters θ, marginalized over several emission-line amplitudes α_j , and including the prior on these amplitudes, as

$$\begin{eqnarray}&&{ \mathcal L }(D\,| \,\theta )=\int \,d{\boldsymbol{\alpha }}\,{ \mathcal L }(D\,| \,\theta ,{\boldsymbol{\alpha }})\,p({\boldsymbol{\alpha }}),\end{eqnarray} \tag{ E2 }$$

where D is the data, θ are the parameters of the model not including the line amplitudes, and α is the vector of all line amplitudes {α_j }. We can then write the log-likelihood, conditional on both θ and α , as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L }(D\,| \,\theta ,{\boldsymbol{\alpha }}) & = & {\left({\boldsymbol{\Delta }}-{\boldsymbol{G}}{\boldsymbol{\alpha }}\right)}^{{\mathsf{T}}}\,{{\rm{\Sigma }}}^{-1}\,({\boldsymbol{\Delta }}-{\boldsymbol{G}}\,{\boldsymbol{\alpha }}),\\ {\boldsymbol{\Delta }} & = & \{{D}_{\lambda }-{M}_{\lambda }\},\\ {\boldsymbol{G}} & = & \{{G}_{\lambda },j\},\\ {G}_{\lambda ,j} & \equiv & { \mathcal N }(\lambda \,| \,(1+z)\,{\lambda }_{0,j},\displaystyle \frac{{\sigma }_{v}}{c}\,(1+z)\,{\lambda }_{0,j}),\end{array}\end{eqnarray} \tag{ E3 }$$

where Σ is the variance–covariance matrix of the data, M_λ is the model prediction without emission lines, Δ is a vector of residuals from the line-free model, G is a normalized Gaussian profile with a given width σ_v and center (specified as part of θ).

The key is that by using α_j only for the amplitude of the lines (i.e., not for the shape or center), the conditional log-likelihood is linear with respect to α . The likelihood can then be factored into two components. The first component is the likelihood of the model conditional on the maximum-likelihood estimate (MLE) of the line amplitudes, $\hat{{\boldsymbol{\alpha }}}$. This scales the second component, which is a multidimensional Gaussian describing the dependence of the likelihood on the line amplitudes α . This Gaussian is centered on $\hat{{\boldsymbol{\alpha }}}$ and has a width described by the covariance matrix ${{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}}$. Our expression for the likelihood function is then

$$\begin{eqnarray}&&{ \mathcal L }(D\,| \,\theta ,{\boldsymbol{\alpha }})={ \mathcal L }(D\,| \,\theta ,\hat{{\boldsymbol{\alpha }}})\,\displaystyle \frac{{ \mathcal N }({\boldsymbol{\alpha }}\,| \,\hat{{\boldsymbol{\alpha }}},{{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}})}{{ \mathcal N }(\hat{{\boldsymbol{\alpha }}}\,| \,\hat{{\boldsymbol{\alpha }}},{{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}})}.\end{eqnarray} \tag{ E4 }$$

The denominator is a non-unity value required to make the equality hold for ${\boldsymbol{\alpha }}=\hat{{\boldsymbol{\alpha }}}$. When the prior is uniform (p( α ) is constant), the integral in the marginalized likelihood of Equation (E2)—ignoring a constant factor that is independent of θ and α —reduces to

$$\begin{eqnarray}&&{ \mathcal L }(D\,| \,\theta )=\displaystyle \frac{{ \mathcal L }(D\,| \,\theta ,\hat{{\boldsymbol{\alpha }}})}{{ \mathcal N }(\hat{{\boldsymbol{\alpha }}}\,| \,\hat{{\boldsymbol{\alpha }}},{{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}})}.\end{eqnarray} \tag{ E5 }$$

We now derive the MLE for α . By taking the derivative of the likelihood with respect to α and setting the result to 0, we find the MLE value $\hat{{\boldsymbol{\alpha }}}$ as

$$\begin{eqnarray}&&\hat{{\boldsymbol{\alpha }}}={\left({{\boldsymbol{G}}}^{{\mathsf{T}}}{{\rm{\Sigma }}}^{-1}{\boldsymbol{G}}\right)}^{-1}\,{{\boldsymbol{G}}}^{{\mathsf{T}}}\,{{\rm{\Sigma }}}^{-1}\,{\boldsymbol{\Delta }}.\end{eqnarray} \tag{ E6 }$$

We are also interested in the widths or uncertainty on the MLE value, described by the covariance matrix ${{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}}$. This is obtained from the inverse second derivative of the likelihood with respect to α , resulting in

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}}={\left({{\boldsymbol{G}}}^{{\mathsf{T}}}{{\rm{\Sigma }}}^{-1}{\boldsymbol{G}}\right)}^{-1}.\end{eqnarray} \tag{ E7 }$$

Notably, by doing this calculation in vector space, we explicitly model possible covariances between the line amplitudes; this is important in cases where the emission lines overlap.

Using the MLE solutions for the emission-line amplitudes is correct for a uniform prior on α . In some cases, an informative prior on α is desirable. This is discussed in the next section.

To incorporate priors we can multiply the factorized likelihood in Equation (E4) by the prior. This is straightforward if the prior is another Gaussian with mean $\breve{{\boldsymbol{\alpha }}}$ and covariance ${{\rm{\Sigma }}}_{\breve{{\boldsymbol{\alpha }}}}$: that is, $p({\boldsymbol{\alpha }})={ \mathcal N }({\boldsymbol{\alpha }}\,| \,\breve{{\boldsymbol{\alpha }}},{{\rm{\Sigma }}}_{\breve{{\boldsymbol{\alpha }}}})$. The product of two Gaussians is another Gaussian; using well-known formulae for the products of normal distributions, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L }(D\,| \,\theta ,{\boldsymbol{\alpha }})p({\boldsymbol{\alpha }}) & = & { \mathcal L }(D\,| \,\theta ,\hat{{\boldsymbol{\alpha }}})\,\displaystyle \frac{{ \mathcal N }(\hat{{\boldsymbol{\alpha }}}\,| \,\breve{{\boldsymbol{\alpha }}},{{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}}+{{\rm{\Sigma }}}_{\breve{{\boldsymbol{\alpha }}}})}{{ \mathcal N }(\hat{{\boldsymbol{\alpha }}}\,| \,\hat{{\boldsymbol{\alpha }}},{{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}})}\\ & & \times { \mathcal N }({\boldsymbol{\alpha }}\,| \,\bar{{\boldsymbol{\alpha }}},{{\rm{\Sigma }}}_{\bar{{\boldsymbol{\alpha }}}}),\end{array}\end{eqnarray} \tag{ E8 }$$

where $\bar{{\boldsymbol{\alpha }}}$ and ${{\rm{\Sigma }}}_{\bar{{\boldsymbol{\alpha }}}}$ are the mean and covariance of the new Gaussian. In fact, this new Gaussian does not need to be evaluated, as the integral in the marginalized likelihood of Equation (E2) reduces to

$$\begin{eqnarray}&&{ \mathcal L }(D\,| \,\theta )={ \mathcal L }(D\,| \,\theta ,\hat{{\boldsymbol{\alpha }}})\,\displaystyle \frac{{ \mathcal N }(\hat{{\boldsymbol{\alpha }}}\,| \,\breve{{\boldsymbol{\alpha }}},{{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}}+{{\rm{\Sigma }}}_{\breve{{\boldsymbol{\alpha }}}})}{{ \mathcal N }(\hat{{\boldsymbol{\alpha }}}\,| \,\hat{{\boldsymbol{\alpha }}},{{\rm{\Sigma }}}_{\hat{{\boldsymbol{\alpha }}}})}.\end{eqnarray} \tag{ E9 }$$

We do not explore non-Gaussian priors, as these are considerably more difficult to implement.

The Cloudy implementation in FSPS includes 128 emission lines; any or all of them may be optimized or marginalized over using this scheme. In practice, it is suggested to only optimize or marginalize over emission lines included in an observed spectrum. The marginalization calculation is fast and adds little computational overhead in most cases. Note that neither optimization nor marginalization respects emission-line physics. This means that, for example, emission-line doublets may be assigned forbidden ratios if the data prefer it. The added emission lines may also have negative luminosities. In objects with weak or no nebular emission, the user may wish to turn off the optimization or limit it to strong emission features to prevent it from instead erroneously modeling differences between the model and observed continuum. This is particularly true in cases where a weak or nonexistent emission line lies on top of an informative absorption feature, such as Hα or Hβ.

When applying the Gaussian priors, the center of the prior for each line is taken as the prediction from the Cloudy implementation. This means that the lines are assumed to be powered by normal star formation, but deviations are permitted. The width of the prior is specified by the user, in units of (true/predicted) luminosity.