Deriving Physical Properties from Broadband Photometry with Prospector: Description of the Model and a Demonstration of its Accuracy Using 129 Galaxies in the Local Universe

In order to extract emission line fluxes, we perform simultaneous fits to the [O iii] 4959, [O iii] 5007, Hβ, [N ii] 6549, [N ii] 6583, and Hα features. The simultaneous fits allow overlapping lines to be properly separated. Each emission and absorption feature is modeled with a Gaussian, with a variable amplitude and width. The center of each line is fixed relative to one another, but the overall redshift of the model is allowed to vary around the fiducial redshift provided by the Brown et al. (2014) catalog. For doublet emission lines ([O iii] and [N ii]), the relative amplitudes of the emission lines are fixed to their fundamental ratios, and the line widths are also fixed to be equal.

The best-fit Prospector-α spectrum (see Section 4.2) is adopted as the stellar continuum. A simple approach to measuring emission line strength from spectra in the case of mixed absorption and emission is to separately measure the underlying absorption in the model stellar continuum, and add that to the observed emission line strength. However, the observed spectra often show gas emission resolved separately from stellar absorption, especially when the width of the stellar absorption line profile is dominated by age-dependent stellar broadening. In these cases, the line emission fills the center of the line, but the wings of the absorption are not filled with emission. To account for this complex line profile, the underlying stellar continuum is subtracted from the observed spectra, and the emission line flux is measured from the residuals.

Before subtracting the stellar continuum model, it must first be smoothed to the appropriate resolution and normalized to the observed spectrum. The observed spectral resolution is controlled by instrumental broadening, atmospheric effects, and the intrinsic line of sight velocity dispersion. The appropriate smoothing resolution is determined by performing a first-pass fit to the observed spectrum with the Gaussian model for absorption and emission lines. The two emission lines with the highest equivalent width are selected, and the average of their line widths is used to smooth the model spectrum. This implicitly assumes the gas and stars have similar velocity dispersions. The normalization is done by fitting two separate linear components for the continuum: one between 4700 Å $\lt \,\lambda \,\lt \,5100$ Å, and one between 6350 Å $\lt \,\lambda \,\lt \,6680$ Å. This is done for both the model spectrum from the best-fit to the photometry and the observed spectrum, and the ratio between the linear components is used to renormalize the model spectrum separately in each wavelength range. Since the model is fit to the photometry, and the observed photometry and observed spectra have consistent normalizations by design (Brown et al. 2014), this normalization factor is typically <5%.

The Brown et al. (2014) spectral atlas does not include error estimates for the spectral fluxes. To estimate errors, the smoothed, normalized stellar continuum model is subtracted from the observed spectra, with known emission lines masked. A Gaussian is fit to the histogram of the residuals. The error on the spectra fluxes is taken to be the width of the Gaussian. This error is typically between 3%–15% of the flux. This provides a conservative estimate of the true instrumental noise, as it also includes inevitable mismatches between the model spectrum and the observed spectrum.

To derive errors on the line fluxes, a bootstrap analysis is performed: for each pixel, the flux is perturbed with a random value drawn from a Gaussian centered at zero, with a width equal to the errors derived above, then the line fluxes are re-fit. This is performed 100 times for each galaxy, and the final line fluxes are taken to be the 50th percentile of this distribution, with the +/− 1σ errors taken to be the 84th and 16th percentiles, respectively. For Hα and Hβ, an additional error due to uncertainty in the absorption correction is added, equal to the 1σ variance in the underlying model stellar absorption as sampled from the model posteriors.

D_n4000 is measured using the Balogh et al. (1999) definition: the ratio of the average flux density ${{\rm{F}}}_{\nu }$ in the narrow wavelength bands 3850–3950 Å and 4000–4100 Å. This is narrower and thus less sensitive to reddening effects than the original Bruzual (1983); Hamilton (1985) definition. This definition is used to measure D_n4000 in both the observed and the model spectra, and both are convolved to a common resolution of 200 km s⁻¹ before measurement. This is done to remove any potential systematics from different spectral resolutions between the model and the observed spectrum.

Hδ is measured separately from the other line features. This is done for two reasons: (1) the Hδ emission and absorption are often of comparable size, resulting in a complex line morphology with a weak net signal (often poorly fit by a Gaussian), and (2) because nearby CN features affect simple linear continuum fits, which introduces a dependence on the detailed elemental abundance ratios (Prochaska et al. 2007). The relatively low S/N of the optical spectra also make the measurement of the intrinsically weak Hδ feature challenging.

To measure Hδ, fixed wavelength indices are defined over which to measure the continuum and the strength of the Hδ feature (e.g., Worthey et al. 1994; Worthey & Ottaviani 1997; Nelan et al. 2005). Defining these indices is challenging, since the intrinsic Hδ absorption line profile is sensitive to the age of the composite stellar populations (Prochaska et al. 2007),. For example, K-dwarf stars have Hδ full-widths at half-maximum (FWHMs) of ∼1–2 Å, while A stars have FWHMs of ∼20 Å, with the wings of the absorption profile extending ∼50 Å from the center (Prochaska et al. 2007).

In order to accurately measure the strength of the Hδ feature for a variety of line profile shapes, we thus adopt a wide and narrow set of indices. The narrow set of indices measures Hδ as the average flux between 4095 Å < λ < 4114 Å, the blue continuum as the average flux between 4072 Å < λ < 4095 Å, and the red continuum as the average flux between 4114 Å < λ < 4130 Å. The broad set of indices measures Hδ as the average flux between 4082 Å < λ < 4124 Å, the blue continuum as the average flux between 4042 Å < λ < 4082 Å, and the red continuum as the average flux between 4124 Å < λ < 4151 Å.

In accordance with the age dependence of the line profile, the wide index is used to measure the strength of the Hδ feature for galaxies with a half-mass assembly time of $\lt 8\,\,\mathrm{Gyr}$, while the narrow index is adopted for galaxies with a half-mass assembly time of $\gt 8\,\,\mathrm{Gyr}$. The half-mass assembly time is measured directly from the model posteriors for each galaxy. This classification scheme qualitatively correlates the width of the Hδ feature in the observed spectra. The same index to measure the Hδ strength is used in both the model and the observations, and the model measurement is performed at the same 200 km s⁻¹ resolution. In the observed spectra, this method of measuring line fluxes inherently measures the sum of Hδ absorption and emission. Errors are derived with the same bootstrapping method described in Section 2.2.1.

Stellar mass, defined as the mass of existing stars and stellar remnants, is a free parameter in the Prospector-α model. Stellar emission powers all sources of radiation in our model (gas, dust, and stars); thus, the stellar mass sets the overall normalization of the SED. Since the normalization constant is defined in mass instead of luminosity, it often has covariances with the SFH parameters. A flat prior is adopted on the logarithm of the stellar mass, 5 < log(M/M_⊙) < 14, and a Chabrier (2003) IMF is used.

A nonparametric SFH prescription is adopted in the Prospector-α model. Using a nonparametric SFH has the great advantage of avoiding the poorly quantified systematics introduced by using a parameterized SFH (Conroy 2013): for example, exponentially declining τ models inherently couple the early-time to late-time SFH (Simha et al. 2014). These systematics are particularly challenging when determining proper errors on the derived SFH.

In the literature, nonparametric SFHs have primarily been used when performing full-spectrum fitting to high-resolution, high S/N data, e.g., STECMAP (Ocvirk et al. 2006) or VESPA (Tojeiro et al. 2007), or both low-resolution spectra and photometry (Kelson et al. 2014; Dressler et al. 2016). This is because, when performing classical minimization routines, non-parametric SFHs are very flexible and have too many degenerate parameters, leading to poorly behaved and noisy best-fit solutions. Bayesian statistics are more suitable for this problem, particularly when coupled with a Monte Carlo Markov chain (MCMC) sampler which can fully explore the posterior probability for poorly constrained parameters. Using a non-parametric SFH with Bayesian statistics will result in clean, accurate errors on SFH parameters compared to previous approaches, at the expense of longer computational time.

Fundamentally, the selection of the number, size, and location of SFH bins is a tradeoff between the required computational time and maximizing the amount of SFH information extracted from the data. This is challenging for a model which is fit to a large sample of galaxies with diverse SFHs, as the recoverability of the SFH given a set of data is highly dependent on the SFH of the galaxy itself (e.g., Tojeiro et al. 2007). To inform our selection of the SFH bins, we use a procedure which is described in detail in B. D. Johnson et al. (2017, in preparation). Here, it is briefly summarized.

The goal is to estimate the fractional uncertainty on the SFH bins given a set of data (i.e., broadband filters) and the associated noise parameters. The fractional uncertainty on each bin is used to select the bin spacing and size. This is done with the Cramer–Rao bound, which is the matrix inverse of the Fisher information matrix. The Fisher information matrix measures the fractional curvature of the likelihood function near the location of maximum likelihood: small curvature implies a broad posterior, and thus more uncertainty in the parameters. In practice, the SFH recoverability is estimated assuming only a single metallicity and no dust. These simplifying assumptions mean that this procedure will specify a lower limit on the uncertainties in the bin sizes.

As a simple estimate, two exponentially declining SFHs are used: a star-forming galaxy with τ = 1 Gyr observed after t = 1 Gyr, and a quiescent galaxy with τ = 1 Gyr observed after 10 Gyr. This procedure suggests that the six time bins shown in the bottom six panels of Figure 2 recover all SFH information for the chosen model galaxies with the Brown et al. (2014) photometry. The spectrum of each SFH bin is shown separately in Figure 4. The model SFR is constant within each bin.

Zoom In Zoom Out Reset image size

The most straightforward way to fit these bins to the data would be to fit directly for the mass in each bin. In practice, the fractional mass in each bin (i.e., the "shape" of the SFH) is used instead, and the mass normalization, M, is fit as a separate parameter. This separation helps MCMC convergence by restricting the SFH parameter space to a smaller, well-defined volume, and also provides cleaner posteriors. We fit for the weight in each bin, f_n, such that

$$\begin{eqnarray}&&{m}_{n}=\displaystyle \frac{{t}_{n}{f}_{n}}{{\displaystyle \sum }_{N}{t}_{n}{f}_{n}},\end{eqnarray} \tag{ 1 }$$

where m_n is the fractional mass in each bin, N is the number of bins, and t_n is the amount of time in each bin in years. The meaning of f_n is more clear in the following equation:

$$\begin{eqnarray}&&{\mathrm{SFR}}_{n}={\rm{M}}\displaystyle \frac{{m}_{n}}{{t}_{n}}={\rm{M}}\displaystyle \frac{{f}_{n}}{{\displaystyle \sum }_{N}{t}_{n}{f}_{n}},\end{eqnarray} \tag{ 2 }$$

where M is the total mass and SFR_n is the SFR in a bin. The weights f_n are thus proportional to the sSFR (and the SFR) within each time bin. This maps directly into the prediction of physical quantities such as the age of the stellar populations and the strength of spectral absorption and emission features. We require that ${f}_{n}\geqslant 0$, i.e., no bin can have a negative amount of stars formed.

In order to ensure a unique mapping between the weights f_n and the shape of the resulting SFH, the weights are constrained during the fitting process such that

$$\begin{eqnarray}&&\displaystyle \sum _{n}\ {f}_{n}=1.\end{eqnarray} \tag{ 3 }$$

In practice, given the sampling algorithm, applying this constraint directly would be highly inefficient. This is because in each sampling step, the sampler chooses a new set of parameters without taking constraints into account, then tests whether the parameters are within the constraints: for this constraint, the probability of randomly chosen parameters which fulfill this constraint is vanishingly small.

To alleviate this inefficiency, we instead allow $N-1$ bin fractions to be free model parameters with the constraint that

$$\begin{eqnarray}&&\displaystyle \sum _{n-1}\ {f}_{n}\leqslant 1.\end{eqnarray} \tag{ 4 }$$

The Nth bin fraction is then determined by

$$\begin{eqnarray}&&{f}_{N}=1-\displaystyle \sum _{1}^{N-1}{f}_{n}.\end{eqnarray} \tag{ 5 }$$

These constraints means that, when the data have little to no constraining power over the SFH, the fractional variables f_n will follow a flat Dirichlet distribution. We have verified numerically that this puts a Dirichlet prior on both the N-1 explicit fractional variables and the Nth implicit fractional variable. The marginal probability distribution for a single f_n, which corresponds to the effective Dirichlet prior on each f_n, is

$$\begin{eqnarray}&&\mathrm{PDF}={(1-x)}^{N-2}\displaystyle \frac{{\rm{\Gamma }}(N)}{{\rm{\Gamma }}(N-1)}.\end{eqnarray} \tag{ 6 }$$

which is a Beta distribution. So for the SFR in a single bin with N = 6 bins, the prior probability distribution is $p{({f}_{n})\propto (1-{f}_{n})}^{4}$ with $0\leqslant {f}_{n}\leqslant 1$.

The upshot of the above parameterization for the SFH is that there is a roughly Gaussian-shaped prior on the logarithm of the sSFR in each time bin, i.e., on log(SFR _bin/M_tot). The center of this prior is at log(1/t_univ) ∼ −10.1 yr⁻¹, which is a constant SFR. The FWHM of the Gaussian prior is roughly 1.5 dex. Thus, in the absence of strong evidence from the data (where "strong" is determined by the FWHM of the prior), Prospector-α will assign a constant SFR(t) to a galaxy.

A clear advantage of the nonparametric SFH formulation over libraries of SFHs from simulations is that the prior on sSFR(t) can be explicitly written down as has been done here. The SFH posterior for an individual galaxy is a combination of the prior and the constraining power of the data. By knowing the prior, it is straightforward to discern how much of the posterior comes directly from the data, and how much of it comes from assumptions that are built-in to every SFH formulation.

Stellar metallicity affects the optical-to-NIR flux ratio (see Figure 3), which is important in determining ages, dust attenuation, and masses (Bell & de Jong 2001; Mitchell et al. 2013). The stellar metallicity also helps to set the normalization and shape of the SED blueward of the Lyman limit, which is used to generate emission line fluxes.

When fitting galaxy SEDs, it is important to interpolate between metallicities rather than using a discrete set of metallicities, else the resulting stellar mass estimates will be biased (Mitchell et al. 2013). We refine this approach further and use a metallicity distribution function rather than interpolated metallicities. Simple interpolation between stellar spectra of different metallicities will introduce non-physical features in the interpolated spectra, because stellar SEDs have a nonlinear response to changes in metallicity. The metallicity distribution function avoids interpolation, producing more physically motivated SEDs.

Stellar metallicity does not vary with age in Prospector-α, and there is no age–metallicity relationship implemented. This is in contrast to expectations in galaxies, where an age–metallicity relationship is expected due to the build-up of heavy elements with the passing of cosmic time. In testing, we have found that implementing an age–metallicity relationship can impart a unique shape to the UV-NIR SED which cannot be replicated by a fixed stellar metallicity, though this depends on both the assumed SFH and metal enrichment history. For example, quiescent galaxies are less sensitive to an age–metallicity relationship, as their SFR timescales are shorter. The complexity of implementation combined with the poor constraining power of broadband photometry prevent inclusion of an age–metallicity relationship in Prospector-α. Ideally, implementation would be guided by performing a detailed comparison of stellar metallicities recovered by Prospector-α to robust spectroscopic metallicities from a sample of galaxies with a wide range in sSFR. This is a topic for future work.

The range in metallicity in our model is determined by the coverage of the Padova isochrones (Marigo & Girardi 2007; Marigo et al. 2008) in FSPS, which extend from $-2.00\lt \mathrm{log}({\rm{Z}}/{{\rm{Z}}}_{\odot })\lt 0.20$. The metallicity distribution function is implemented by summing simple stellar populations at discrete metallicity values, using a triangular weighting scheme. This is illustrated in Figure 3. The metallicity is specified in units of log(Z/Z${}_{\odot }$), where log(Z/Z${}_{\odot }$) defines the center of the triangular weighting scheme. A flat prior over −2.0 < log(${\rm{Z}}/{{\rm{Z}}}_{\odot })\lt 0.2$ is used in the Prospector-α model.

We use the two-component Charlot & Fall (2000) dust attenuation model, which postulates separate birth-cloud and diffuse dust screens. The birth-cloud dust simulates the embedding of young stars in molecular clouds and H ii regions, the net effect of which is extra attenuation toward young stars. We allow the normalization of the birth-cloud and diffuse dust components to vary separately in the fits, along with a power-law index describing the shape of the attenuation curve for the diffuse dust component. The variance of the SED with these parameters is shown in Figure 3.

It is difficult to achieve tight posteriors on dust attenuation parameters from UV-NIR broadband photometry alone, due to the well-known age–dust degeneracy (Papovich et al. 2001) (see Section 3.1.6 for further discussion). This degeneracy can be broken either by directly measuring the dust emission with MIR or FIR broadband filters (Burgarella et al. 2005), or with excellent photometric coverage of the age-sensitive D_n4000 break (MacArthur et al. 2004). In the Brown et al. (2014) sample, the MIR is sampled via WISE and Spitzer photometry, and for a limited subsample, the FIR is measured with Herschel photometry.

The birth-cloud component (${\hat{\tau }}_{\lambda ,1}$) in our model has an attenuation curve that scales with wavelength as ${\lambda }^{-1}$, and attenuates stellar emission only from stars formed in the last t = 10 Myr, the typical timescale for the disruption of a molecular cloud (Blitz & Shu 1980). It also attenuates nebular emission. It has the following functional form:

$$\begin{eqnarray}&&{\hat{\tau }}_{\lambda ,1}={\hat{\tau }}_{1}{(\lambda /5500\mathring{\rm A} )}^{-1.0}.\end{eqnarray} \tag{ 7 }$$

The diffuse component (${\hat{\tau }}_{\lambda ,2}$) has a variable attenuation curve, and attenuates all stellar and nebular emission from the galaxy. For the wavelength scaling, we use the following prescription from Noll et al. (2009):

$$\begin{eqnarray}&&{\hat{\tau }}_{\lambda ,2}=\displaystyle \frac{{\hat{\tau }}_{2}}{4.05}[{k}^{\prime }(\lambda )+D(\lambda )]{\left(\displaystyle \frac{\lambda }{{\lambda }_{V}}\right)}^{n}.\end{eqnarray} \tag{ 8 }$$

The free parameters in this equation are ${\hat{\tau }}_{\lambda ,2}$, which controls the normalization of the diffuse dust, and n the diffuse dust attenuation index. $k^{\prime} (\lambda$) is the (fixed) Calzetti et al. (2000) attenuation curve. D(λ) is a Lorentzian-like Drude profile describing the UV dust bump. We tie the strength of the UV dust absorption bump to the best-fit diffuse dust attenuation index, following the results of Kriek & Conroy (2013). We adopt a flat prior over $0\lt {\hat{\tau }}_{\lambda ,2}\lt 4.0$, a flat prior over $0\lt {\hat{\tau }}_{\lambda ,1}\lt 4.0$, and a flat prior over $-2.2\lt {\text{}}n\lt 0.4$. The upper limit on ${\text{}}n$ is chosen to disallow a flat attenuation curve, which would cause ${\hat{\tau }}_{\lambda ,2}$ to be nearly fully degenerate with the normalization of the SED.

The similar effects of ${\hat{\tau }}_{\lambda ,1}$ and ${\hat{\tau }}_{\lambda ,2}$ on the stellar SED means that they are often degenerate (though this is SFH-dependent). It is important to distinguish between these parameters to properly predict emission line equivalent widths. Previous work suggests that the total optical depth toward nebular emission lines is ∼twice that of the stellar component (Calzetti et al. 1994; Price et al. 2014), though the exact value varies with SFR and galaxy inclination (Wild et al. 2011). Translating this result into the Charlot & Fall (2000) model implies that ${\hat{\tau }}_{\lambda ,1}\sim {\hat{\tau }}_{\lambda ,2}$, as ${\hat{\tau }}_{\lambda ,1}$ applies to the entire galaxy, whereas star-forming regions only experience attenuation from ${\hat{\tau }}_{\lambda ,2}$. We adopt a joint prior on the ratio of the two: $0\lt ({\hat{\tau }}_{\lambda ,1}/{\hat{\tau }}_{\lambda ,2})\lt 2.0$, which allows some reasonable variation around the fiducial results in the literature.

To model the dust emission from galaxies, we assume energy balance, where all starlight attenuated by the dust is re-emitted in the IR (da Cunha et al. 2008). This assumption means that the MIR and FIR photometry are additional constraints both on the total amount of dust attenuation, and on the dust-free stellar SED. Simultaneous modeling of the dust and stellar emission is key to accurate SFRs, particularly for galaxies with lower sSFRs where dust-heating via old stellar populations becomes more important (e.g., Utomo et al. 2014). While a useful tool, energy balance cannot be used to determine L_IR from the UV-MIR SED alone without making some assumptions about the shape of the IR SED (see Appendix C for further discussion).

We use the Draine & Li (2007) dust emission templates to describe the shape of the IR SED. This model is based on the silicate-graphite-PAH model of interstellar dust (Mathis et al. 1977; Draine & Lee 1984). The size distribution of carbonaceous and silicate grains is calibrated to be consistent with the wavelength-dependent extinction in the Milky Way. We note that in our model, the shape of the dust emission is modeled independently from the shape of the dust attenuation curve. This is a reasonable implementation because the shape of the attenuation (not extinction) curve is largely set by the geometry and column depth of the dust, rather than the composition of the dust (Chevallard et al. 2013). In the future, it may be useful to tie the strength of the 2175 Å extinction bump to the strength of the PAH emission, as it has long been suggested that these features both arise from the same dust grains (Stecher & Donn 1965; Draine 1989).

The Draine & Li (2007) model has three free parameters controlling the shape of the IR SED: U_min, ${\gamma }_{{\rm{e}}}$, and Q_PAH. The effect of these three parameters on the resulting SED is demonstrated in Figure 3.

U_min and ${\gamma }_{{\rm{e}}}$ together control the shape and location of the thermal dust emission bump in the IR SED. Physically, they describe the fraction of dust mass exposed to the distribution of starlight intensities U, which is parameterized in the Draine & Li (2007) model by a combination of a delta function and a power-law distribution over U_min $\lt \,U\lt {U}_{\max }$:

$$\begin{eqnarray}\displaystyle \frac{{{dM}}_{\mathrm{dust}}}{{dU}} & = & (1-{\gamma }_{{\rm{e}}}){M}_{\mathrm{dust}}\delta (U-{U}_{\min })\\ & & +{\gamma }_{{\rm{e}}}{M}_{\mathrm{dust}}\displaystyle \frac{(\alpha -1)}{({U}_{\min }^{1-\alpha }-{U}_{\max }^{1-\alpha })}{U}^{-\alpha },\end{eqnarray} \tag{ 9 }$$

with ${{dM}}_{\mathrm{dust}}$ is the mass of dust heated by starlight intensities between U and U + dU, ${M}_{\mathrm{dust}}$ is the total dust mass, and α is a power-law index. Draine & Li (2007) finds that MIR photometry of galaxies in the SINGS survey (Kennicutt et al. 2003) is well-reproduced with $\alpha =2$ and ${U}_{\max }={10}^{6}$, which we adopt here. ${M}_{\mathrm{dust}}$ is determined by the normalization of the SED.

Specifically, U_min represents the minimum starlight intensity to which the dust mass is exposed, and ${\gamma }_{{\rm{e}}}$ represents the fraction of dust mass which is exposed to this minimum starlight intensity. The minimum starlight intensity roughly represents the intensity experienced by dust in the general diffuse interstellar medium (ISM). The inclusion of this delta function at the minimum starlight intensity is the major difference between the Draine & Li (2007) and Dale & Helou (2002) models for ${{dM}}_{\mathrm{dust}}/{dU}$.

The final free parameter, Q_PAH, describes the fraction of total dust mass that is in PAHs. This parameter is particularly important because a substantial fraction or a majority of the MIR emission comes from strong PAH emission features.

The PAH ionized fractions and absorption cross-sections in Draine & Li (2007) are updated from the Li & Draine (2001) model by using observations of SINGS galaxies with 5–38 μm Spitzer IRS spectra (Smith et al. 2007). Draine & Li (2007) provide IR spectral templates for 0.1 < Q_PAH < 4.58. We note that with the excellent MIR coverage of the Brown et al. (2014) spectral atlas, Q_PAH can be very well-determined photometrically (see Appendix A for further details).

We adopt a flat prior of $0.1\,\lt$ U_min $\lt \,15$ and $0.0\,\lt$ ${\gamma }_{{\rm{e}}}$ $\lt \,0.15$. For Q_PAH, we extrapolate the Spitzer IRS spectra out past Q_PAH = 4.58, and adopt a flat prior over $0.1\,\lt$ Q_PAH $\lt \,7.0$. The extrapolation is linear in both flux and Q_PAH. In brief, the adopted priors are consistent with both the SINGS sample (Draine et al. 2007) and the Brown et al. (2014) galaxies with Herschel photometry, and adopting more permissive priors would bias the L_IR, SFR, and dust attenuation in galaxies without Herschel photometry. We discuss the origin and effects of these priors further in Appendix C.

We generate nebular continuum and line emission using the CLOUDY (Ferland et al. 1998, 2013) implementation within FSPS. This is described in detail by Byler et al. (2016), and summarized briefly here.

The CLOUDY implementation in FSPS is a grid in the following parameters: (1) simple stellar population (SSP) age, (2) SSP and gas-phase metallicity, and (3) the ionization parameter, U, which is the ratio of ionizing photons to the total hydrogen density. Before the SSPs within FSPS are combined into composite stellar populations, the number of ionizing photons in each SSP is calculated and removed from the SSP SEDs. The energy removed from the SSPs is then rerouted into line luminosities and nebular continua. This is only done for SSPs with ages $\lt 10\,\,\mathrm{Myr}$, as the ionizing contribution from older SSPs is very small. This also implicitly assumes an ionizing photon escape fraction of zero.

Within the Prospector-α model, the gas-phase metallicity is set equal to the model stellar metallicity, and the ionization parameter U is fixed to 0.01 in all fits. We note that the ionization parameter has very little effect on the Balmer emission line luminosities (e.g., Hα or Hβ), though it has a strong effect on the luminosity of forbidden lines ([O ii], [O iii], etc.). In principle, it could be fit as well, though broadband photometry alone provides very little constraining power. For the Brown et al. (2014) sample, the contribution to the broadband photometry from emission lines is relatively low (≲5%) due to the relatively low sSFRs, so U is fixed for simplicity. In the future, it may be useful to allow U to vary, perhaps with a joint prior between sSFR and U, as suggested by recent work which finds a strong relationship between the two (e.g., Dickey et al. 2016). This would be particularly relevant for galaxies with high sSFR, where the nebular contribution to the broadband flux can become significant (for typical star-forming galaxies at z ∼ 1–2, it is 5%–10% or greater in certain rest-frame optical bands; Pacifici et al. 2015).

In the Prospector-α model, after production, emission lines fluxes experience attenuation from both the birth-cloud and diffuse dust component (see Section 3.1.3). The emission line equivalent widths are unaffected by attenuation from the diffuse dust component, as it attenuates the stellar emission equally. However, the emission line equivalent widths are sensitive to the birth-cloud attenuation, since it only attenuates emission from young stars and nebular emission. Both nebular continuum and line emission are added to the full SED model, and their effects are included in the broadband photometry.

Here we describe the typical posteriors for an application of the Prospector-α model to the Brown et al. (2014) photometry. Figure 5 shows a "corner"⁷ plot for NGC 0628, a mildly star-forming galaxy chosen to illustrate several important parameter degeneracies.

Zoom In Zoom Out Reset image size

The panels along the diagonal show the posterior probability functions for each model parameter. The median of each posterior is printed above the diagonals, with the ±1σ credible intervals also indicated. Most model parameters have Gaussian-like posteriors, and the exceptions in the SFH and dust parameters are a result of explicit or implicit priors in the Prospector-α model. The fractional SFH parameters largely have $p{({f}_{n})\propto (1-{f}_{n})}^{4}$, as a result of the implicit prior set by the SFH implementation described in Section 3.1.1. This is not true for the youngest SFH time bin, where the evidence from the photometry is strong enough to reshape the posterior into a Gaussian. The birth-cloud dust parameter posterior, dust1, is truncated where dust1 = dust2, as a result of the dust priors described in Section 3.1.3. Finally, the dust_index parameter, which controls the shape of the dust attenuation curve, is against the prior limit described in Section 3.1.3.

The panels in each column show the correlation between the posteriors for each model parameter. These illustrate several degeneracies which are important in understanding the limited precision for certain parameters, including the dust–metallicity–SFR degeneracy (dust2-logzsol-sfr_fraction1, respectively), and the amount of dust and the shape of the attenuation curve (dust2 and dust_index). The parameters controlling the shape of the dust emission (duste_gamma, duste_qpah, and duste_umin) are also mildly degenerate with one another, due to the limited sampling of the IR SED. These degeneracies can be accurately explored by Prospector due to the on-the-fly model generation and the MCMC sampling.

We construct the Prospector-α galaxy model within Prospector.⁸ Prospector is a Python package which serves as a framework to bring together a galaxy model with a noise model in a posterior probability function. That probability function is then passed to algorithms which either maximize or sample the posterior probability function. Prospector has a highly flexible, modular design: input observations can be photometric, spectroscopic, or both; the sampling/minimizing procedure can be changed to suit the problem at hand; and the likelihood function, physical model, and priors are cleanly separated from one another.

We adopt a simple chi-squared noise model, i.e., we model the photometric noise with Gaussian independent errors of known dispersion. We note that Prospector has a number of more sophisticated noise models available, including estimation of covariant noise between sets of photometric filters, and estimating an overall jitter term from the available photometric data.

We use a two-step process to find the posterior probability distribution function for our model. The first step is a simple minimization scheme intended to find the best-fit model parameters. The second step is an MCMC sampler, which is initialized in the region of the best-fit parameters from the previous step.

The minimization step uses the Powell minimization scheme from the scipy Python package. We take advantage of multiprocessing to perform this process. For N available processors, N sets of initial parameter values are generated randomly within the constraints of the priors, and then each processor performs an individual Powell minimization scheme for a specified number of iterations, or to a specified likelihood tolerance, whichever is reached first. The best-fitting model parameters from the Powell routine are used as the initial conditions for the next step.

We use the emcee package (Foreman-Mackey et al. 2013) for the sampling step. emcee is a Python-based MCMC sampler based on an affine-invariant algorithm. In brief, emcee uses an ensemble of walkers to explore probability space. Each walker will attempt a specified number of steps in parameter space. A step is attempted by first proposing a new position, and then calculating the likelihood of the new position based on the data. The walker has a chance to move to the new position based on the ratio of the likelihood of its current position to that of the new position. After a specified number of steps, the algorithm is completed and the position and likelihood of each walker as a function of step number is written out. These positions and likelihoods are used to generate the posterior probability distribution function for each model parameter.

To initialize emcee, we generate walker positions in a Gaussian ball around the best-fitting model parameters from the Powell minimization phase. The dispersion of the Gaussian ball is customized for each parameter. The dispersions are chosen to represent the typical uncertainty in each parameter. If the first minimization process returns a set of best-fit parameters close to the global maximum, then convergence is swift and the probability space around the global maximum will be well-sampled. If the initial conditions are far from the global maximum, the walkers will attempt to reach and explore the global maximum before the process is completed. In some cases the process completes before the walkers have reached an equilibrium around the global maximum (before they have burned in), and the resulting parameter probability distributions will not be useful. We visually inspect plots of the walker location in each parameter as a function of iteration step for each galaxy to ensure they are burned in. If the MCMC process is not burned in, the fit is repeated with more iterations until equilibrium is achieved. Future versions of Prospector will have quantitative criteria to assess convergence.

It is critical that we use an MCMC algorithm to explore our model posterior, as a standard grid approach with 13 parameters and sufficient grid spacing would take up a prohibitive amount of memory. One way to limit the amount of grid space to explore is to use generating functions to sample reasonable combinations of parameters, and create static template libraries from this sampling of parameter space. However, this method necessarily results in opaque priors. On-the-fly stellar populations generation with FSPS combined with MCMC sampling allows Prospector to explore model parameter space free from the complex priors from pre-generated template libraries, at the expense of increased computational time.

Balmer emission lines are produced when hydrogen recombines after being ionized by high-energy photons. The conversion from number of ionizing photons to the intrinsic Hα luminosity is effectively a constant. The primary stellar source of ionizing photons is young, massive stars that live for a very short period of time (<10 Myr; e.g., Kennicutt 1998; Conroy 2013 and references therein). Testing model predictions of the Hα luminosity is thus, in large part, a test of the recent SFR.

In addition to SFR, the dust attenuation and the stellar metallicity affect the emission line luminosity. In the Prospector-α galaxy model, the dust attenuation of the line emission is determined by a combination of three parameters: the optical depth of the diffuse dust, the optical depth of the birth-cloud dust, and the shape of the diffuse dust attenuation curve. The birth-cloud dust and diffuse dust are largely degenerate in fits to the broadband photometry (see Figure 3): the ratio of the two is constrained in the model priors to mitigate uncertainty.

Stellar metallicity affects the Hα and Hβ luminosities by changing the production rate of ionizing photons at a fixed SFR: at lower metallicities, there is an increased number of ionizing photons, resulting in a higher Hα luminosity. To illustrate the effect of stellar metallicity on the Hα luminosities, we use the Kennicutt (1998) conversion between recent SFR and predicted Hα luminosity:

$$\begin{eqnarray}&&\mathrm{SFR}({M}_{\odot }\,{\mathrm{yr}}^{-1})=\displaystyle \frac{L({{\rm{H}}}_{\alpha })}{1.26\times {10}^{41}\ \mathrm{erg}\ {{\rm{s}}}^{-1}}.\end{eqnarray} \tag{ 10 }$$

This is derived using a Salpeter (1955) IMF and assumes a fixed, solar metallicity. We correct our SFRs from a Chabrier to a Salpeter scale by multiplying by a factor of 1.7, which accounts for the additional low-mass stars in a Salpeter IMF.

We plug the SFR inferred from the photometric fit into Equation 10 to give another estimate of Hα luminosity (in addition to the CLOUDY estimates), and attenuate the resulting Hα luminosities with the Prospector-α model dust attenuation.

In the left panel of Figure 14, we compare the analytical, fixed-metallicity Hα luminosities from Equation 10 to the observations. There is an additional scatter of 0.05 dex over the same comparison with the CLOUDY Hα predictions (Figure 9).

Zoom In Zoom Out Reset image size

The difference between the two estimates is because the CLOUDY predictions integrate the ionizing flux directly from the stellar model, thus including the effect of a variable stellar metallicity on the number of ionizing photons. We demonstrate this directly by comparing the ratio of the analytical Hα luminosities to the CLOUDY predictions in the right panel of Figure 14, as a function of model stellar metallicity. The ratio of the two is a nearly monotonic function of stellar metallicity.

Because the gas-phase and stellar metallicities are set equal in the Prospector-α model (Section 3.1.5), the right panel of Figure 14 includes the effect of changing the gas-phase abundance within CLOUDY. However, this effect is negligible compared to the effect of stellar metallicity.

The increased predictive power of the CLOUDY Hα predictions, compared to the Hα predictions from Equation 10, is indirect evidence that the Prospector-α metallicities for star-forming galaxies derived from the broadband photometry are accurate.

We have shown here that neglecting the effect of stellar metallicities can dominate the scatter between model and observed Hα. We thus caution against interpreting observational differences between Hα SFRs and other SFR indicators as variation in SFRs over a short timescale, without first taking stellar metallicity into account.

In this section we demonstrate how the model Hα luminosity is sensitive to the shape and normalization of the observed SED, and argue that full-SED fitting is the best way to predict Hα from the photometry.

Typical mock SEDs at three different Hα luminosities are generated by drawing model parameters that affect Hα luminosity (dust, metallicity, and SFH) randomly within their priors, while other model parameters are held fixed. The typical SEDs at each Hα luminosity are shown in the left panel of Figure 15.

Zoom In Zoom Out Reset image size

The primary signature of increasing Hα luminosity in the galaxy SED is an increase in the total IR and UV luminosities, L_IR + L_UV. This is because L_UV and L_IR are directly correlated with the recent SFR (Kennicutt 1998).

Measuring L_IR + L_UV alone does not uniquely determine the Hα luminosity. The right panel of Figure 15 shows individual model SEDs created at a fixed L_IR + L_UV. The SEDs are color-coded by the resulting observed Hα luminosity, which spans a range of ∼0.7 dex. This variation in Hα luminosity comes from variation in the dust attenuation at Hα wavelengths, which is visible in the right panel of Figure 15 as a gradual increase in L_IR/L_UV with decreasing Hα luminosity.

The upshot of Figure 15 is that accurate predictions of the Hα luminosity require: (1) estimates of the energy budget available for Hα photon production (L_UV + L_IR), and (2) the dust attenuation toward nebular regions at Hα wavelengths.

With only UV-MIR broadband photometry, it is possible to obtain a crude estimate of the total energy budget available for Hα photon production from measurements of L_UV and the UV slope β by using the L_UV/L_IR (known as "IRX")–β relationship⁹ (Meurer et al. 1999; Hao et al. 2011). This method will also give the nebular dust attenuation, if one assumes a shape for the dust attenuation curve and a conversion from stellar to nebular attenuation. However, the IRX–β relationship is sensitive to variations in the intrinsic UV slope from variation in the recent SFH (Kong et al. 2004; Boquien et al. 2012); furthermore, we show in Section 5.2.1 that the dust attenuation curve shows strong galaxy-to-galaxy variation.

Hao et al. (2011) perform the above conversion, from observed UV slope and luminosities to the dust-free UV luminosity, finding an RMS scatter of 0.23 dex between the observed and predicted Hα. However, they use the spectroscopic Balmer decrement to correct the Hα luminosity, rather than attempting to estimate ${A}_{{\rm{H}}\alpha }$ directly. They find that the scatter between ${A}_{{\rm{H}}\alpha }$ and A_FUV without FIR fluxes to be on the order of ∼0.6 magnitudes or ∼0.25 dex (see their Figure 15): this means that using UV-MIR photometric data only would increase the scatter further.

We conclude that the full UV-IR SED modeling approach that we present in this work is the most robust way to estimate Hα luminosities from the broadband photometry, as it requires no assumptions about the intrinsic stellar populations or the shape of the attenuation curve. We note that having direct L_IR measurements from Herschel reduces the uncertainty in both approaches, though it may be more helpful for the IRX–β approach, as L_IR recovery from UV-MIR imaging alone after applying the appropriate priors on the shape of the IR SED is excellent (see Appendix C).

In this section we investigate the scatter between the predicted and observed Balmer emission line luminosities in Figure 9, and to what extent this scatter is consistent with the uncertainty on the model predictions.

We first compare the residuals in Hα to the residuals in Hβ in Figure 16. Since the intrinsic flux ratio between the Balmer lines is nearly invariant, the correlations between emission line residuals are informative. Correlated residuals are due to errors in normalization of the Balmer lines: intrinsic emission line strength (errors in SFR or stellar metallicity), normalization of the dust attenuation, or errors in the underlying continuum absorption. Uncorrelated errors are due to errors in the reddening curve, which is a function of ${\hat{\tau }}_{\lambda ,1}$, ${\hat{\tau }}_{\lambda ,2}$, and n in the Prospector-α model.

Zoom In Zoom Out Reset image size

Figure 16 makes it clear that the dominant source of scatter is from the normalization of the Balmer lines, and errors in the reddening curve are a secondary contribution. This can also be calculated directly from the observed scatter in the Balmer decrements in Figure 11: reddening errors add a scatter of ∼ 0.08 dex to the Hα measurements, compared to an overall scatter of 0.18 dex.

We also investigate whether the observed scatter in emission line luminosities is consistent with predictions from the Prospector-α model. Each fit to the photometry results in some posterior probability distribution for the emission line fluxes. The observational errors on emission line fluxes are much smaller than the typical width of this posterior probability function.¹⁰ Thus, we can use the observations to test the model posteriors.

In Figure 17, we show the sum of the model posteriors for the Hα luminosity. The posteriors are stacked such that the median of each posterior is located at the origin, and the global offset is removed (the offset is shown in the corner of Figure 9). We also show the location of the observations for each of these fits, again stacked relative to the median of each posterior. If the model posteriors accurately describe all relevant emission line production and attenuation processes in these galaxies, then the two histograms would match each other to within Poisson errors.

Zoom In Zoom Out Reset image size

The agreement between the two histograms is good, suggesting that the posterior probability distributions are accurate. This means that the scatter in Figure 16 is largely consistent with the errors. We quantify this statement by calculating the fraction of observations that fall within the model 1σ range (i.e., the 16th and 84th percentiles of the posterior) and 2σ range (the 2.5th and 97.5th percentiles of the posterior). These are shown in the upper-left corner of Figure 17. This test shows that the model error bars are accurate to within 20%–30% for Hα and Hβ.

A relationship between the dust attenuation optical depth and the shape of the attenuation curve in galaxies has been predicted both analytically (Bruzual et al. 1988) and numerically, with radiative transfer codes (Witt et al. 1992; Witt & Gordon 2000; Gordon et al. 2001; Chevallard et al. 2013). The wavelength-dependent scattering properties of dust and the mixed star-dust geometry both contribute to this effect. In brief, at low optical depths, red light is scattered more isotropically and escapes the galaxy, while blue light is more forward-scattered and is subjected to more absorption. This results in a net steepening of the attenuation curve. At high optical depths, in a mixed star-dust geometry, the observed radiation primarily comes from stars at an optical depth less than unity. Redder photons travel more physical distance than bluer photons before absorption or scattering; the net effect is that larger optical depths result in a grayer attenuation curve.

In Figure 18, we plot the optical depth of the diffuse dust versus the slope of the attenuation curve at 5500 Å (d logτ/dλ), compared to a compilation of dusty radiative transfer models from Chevallard et al. (2013). The qualitative agreement is remarkable: as the observed dust optical depth increases, the observed attenuation curves become grayer. There is clear evidence for a systematic relationship between the dust optical depth and the shape of the attenuation curve, although the models predict a somewhat grayer attenuation curve than is observed at large optical depths. We caution that there is a correlation between the optical depth of the dust and the attenuation curve for individual galaxies which could partially mimic this relationship; however, it is not large enough to explain the observed sample trend.

Zoom In Zoom Out Reset image size

The dust attenuation curve is often a fixed parameter when fitting the photometric SEDs of galaxies, even for complex SED fitters with many free parameters. Yet there is substantial evidence that the attenuation curve varies: both theoretically, as discussed above, and observationally, with serious differences between the Milky Way (Cardelli et al. 1989), Small Magellanic Cloud (SMC; Prevot et al. 1984; Bouchet et al. 1985), and Large Magellanic Cloud (LMC; Fitzpatrick 1986), as well as in high redshift galaxies (Buat et al. 2012; Reddy et al. 2015; Salmon et al. 2016).

The assumption of a fixed attenuation curve can have profound effects on the results of SED fits. For example, assuming a fixed Calzetti et al. (2000) attenuation curve will underpredict the 1500 Å luminosity dust correction by ∼2 at low optical depths, and overpredict the dust corrections by a factor of 2–5 at high optical depths (Salmon et al. 2016). Scatter in the relationship between L_IR/L_UV and UV slope could be reduced substantially by taking into account the variation in the slope of the attenuation curve (Buat et al. 2012). Assuming the shape of the attenuation curve also affects the bulk properties of the galaxy population, Marchesini et al. (2009) show that the low-mass slope of the stellar mass function steepens and the normalization decreases by 20%–50%.

Given that the shape of the dust attenuation curve is well-recovered in photometric mock tests (Appendix A), the excellent accuracy in which the model fits predict the reddening from the Balmer decrements (Figure 11), and the biases introduced by a fixed attenuation curve in SED fits, it is wise to use a variable dust law when fitting galaxy SEDs.

In this section we investigate the PAH mass fraction, parameterized in the Draine & Li (2007) dust emission model as Q_PAH. Many studies use broadband photometry in the MIR to estimate total IR luminosities, either by using fixed IR templates (e.g., Daddi et al. 2007; Wuyts et al. 2008; Buat et al. 2009; Whitaker et al. 2014) or by fitting IR models to a limited wavelength range (e.g., Chang et al. 2015). The PAH mass fraction is of particular interest because PAH emission can dominate the integrated flux at MIR wavelengths. Thus, the natural variation of the PAH mass fraction in galaxies has serious consequences for estimating SFRs from photometry: unless Q_PAH can be constrained with available photometry, systematic errors in SFR estimates will be introduced by assuming a fixed IR template.

In Figure 19, we investigate the recovery of Q_PAH from photometric fits to the full SED. We show the Spitzer IRS spectra stacked in bins of median Q_PAH, which come from the posteriors of the fits to the broadband photometry. Only galaxies with log(M) $\gt \,10$ are included, to minimize confounding trends in stellar mass. The bin sizes are chosen so that the number of galaxies is the same for each bin (N = 20). Individual spectra are normalized such that the total flux observed in the Spitzer IRS spectrum is equal to unity, prior to stacking: this means that each galaxy contributes an equal amount of information to the stacks. The data are smoothed for presentation purposes. The following PAH features dominate the total PAH luminosity, and are marked with dashed lines: the 6.2 μm (∼10% of L_PAH), 7.7 μm (∼30% of L_PAH), 8.3 μm + 8.6 μm (∼15% of L_PAH), 11.3 μm (∼15% of L_PAH), 12.6 μm (∼10% of L_PAH, but blended with [Ne ii]), and 17 μm (∼10% of L_PAH) (Smith et al. 2007). Strong nebular emission lines are masked in gray bands.

Zoom In Zoom Out Reset image size

It is clear from Figure 19 that the EW of the PAH features increases as the Q_PAH from the photometric fits increases. This is an excellent validation of the photometric Q_PAH values derived from fitting the full SED. This result implies that PAH emission can be cleanly distinguished from other sources of MIR emission (e.g., circumstellar dust around AGB stars, or AGN; Lange et al. 2016) with full-SED modeling.

Measuring total IR luminosity in PAH features for individual galaxies, rather than stacks, would provide additional compelling spectral evidence for the photometric Q_PAH values. However, this measurement is challenging: it is difficult to define the continuum in the MIR region, resulting in biased EWs by factors of ∼2–3 (Smith et al. 2007). Furthermore, individual PAH EWs are sensitive to the intensity of the radiation field (particularly the presence of an AGN) and the galaxy metallicity (Smith et al. 2007; O'Dowd et al. 2009), which makes the interpretation of PAH EWs complex. A full validation of the IR SED models would involve a simultaneous fit to the observed Spitzer IRS, Akari spectra, and the photometry, with the Draine & Li (2007) models and a model for the IR emission of the AGN. We defer this more careful treatment to a later investigation.

Figure 20 shows the relationship between stellar mass and Q_PAH. There is a clear link with stellar mass, such that Q_PAH is very low below a stellar mass of log(M/M${}_{\odot }$) ∼ 8.8. Assuming a tight relationship between gas-phase metallicity and stellar mass, this is consistent with previous observations which show the PAH mass fraction is very low for galaxies below a gas-phase metallicity of log(O/H) + 12 ∼ 8.1, with considerable natural variation above this threshold (Roche et al. 1991; Engelbracht et al. 2005; Hunt et al. 2005; Draine et al. 2007). Full-SED modeling of the broadband photometry reproduces these trends in the PAH mass fraction. We note that the MIR is particularly well-sampled for these galaxies, with Spitzer/IRAC at 3–8 μm, WISE at 3–22 μm, and Spitzer/MIPS at 24 μm, and that Q_PAH is very well-recovered in mock tests; see Appendix A.

Zoom In Zoom Out Reset image size

It has been demonstrated in this section that Q_PAH, and trends in Q_PAH with stellar mass, can be reproduced with a combination of MIR photometry and full-SED modeling. This, along with encouraging mock test results in Appendix A, suggests that Q_PAH is accurately recovered by fits to the photometry. It is likely that full-SED modeling with MIR photometry can tighten the 8 μm/L_IR relationship (Elbaz et al. 2010, 2011), and thus refine SFR estimates, by including constraints from optical–IR energy conservation. It has recently been shown at z = 2 that L_IR can be accurately estimated with the FSPS stellar populations combined with the Draine & Li (2007) dust emission model (Shivaei et al. 2016a). This is particularly pertinent in light of recent results from Shivaei et al. (2016b), which demonstrate that systematic trends in PAH emission with mass bias global SFR estimates by ∼30% at $z\sim 2$. In future work, we will explore the use of Prospector to estimate accurate L_IR, 8 μm/L_IR, and SFRs of $z\gt 1$ galaxies.