Inference of the physical properties of stellar populations from observed photometry and spectroscopy is a key goal in the study of galaxy evolution. In recent years the quality and quantity of the available data have increased and there have been corresponding efforts to increase the realism of the stellar population models used to interpret these observations. Describing the observed galaxy spectral energy distributions in detail now requires physical models with a large number of highly correlated parameters. These models do not fit easily on grids and necessitate a full exploration of the available parameter space. We present PROSPECTOR a flexible code for inferring stellar population parameters from photometry and spectroscopy spanning UV through IR wavelengths. This code is based on forward modeling the data and Monte Carlo sampling the posterior parameter distribution enabling complex models and exploration of moderate dimensional parameter spaces. We describe the key ingredients of the code and discuss the general philosophy driving the design of these ingredients. We demonstrate some capabilities of the code on several data sets including mock and real data.
['Galaxy evolution', 'Spectral energy distribution', 'Astronomy data modeling', '594', '2129', '1859', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Instrumentation and Methods for Astrophysics']
Stellar Population Inference with Prospector
{'Stellar Population Inference with Prospector': 'Benjamin D. Johnson, 1 Joel Leja, 2, 3, 4 Charlie Conroy, 1 and Joshua S. Speagle 5, 6, 7, 1, ∗  \n- 2\n- 1 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA 3 Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA, USA 4 Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA 5 Department of Statistical Sciences, University of Toronto, Toronto, ON M5S 3G3, Canada', 'ABSTRACT': 'Inference of the physical properties of stellar populations from observed photometry and spectroscopy is a key goal in the study of galaxy evolution. In recent years the quality and quantity of the available data has increased, and there have been corresponding efforts to increase the realism of the stellar population models used to interpret these observations. Describing the observed galaxy spectral energy distributions in detail now requires physical models with a large number of highly correlated parameters. These models do not fit easily on grids and necessitate a full exploration of the available parameter space. We present Prospector , a flexible code for inferring stellar population parameters from photometry and spectroscopy spanning UV through IR wavelengths. This code is based on forward modeling the data and Monte Carlo sampling the posterior parameter distribution, enabling complex models and exploration of moderate dimensional parameter spaces. We describe the key ingredients of the code and discuss the general philosophy driving the design of these ingredients. We demonstrate some capabilities of the code on several datasets, including mock and real data. /github  \nKeywords: galaxies: evolution - galaxies: stellar content - galaxies: fundamental parameters methods: statistical', '1. INTRODUCTION': "Fitting galaxy spectral energy distributions (SEDs) is the primary method for converting galaxy observations into inferred physical properties. These physical properties can in turn be compared to predictions from theories of galaxy evolution. Accordingly, galaxy SED fitting has a rich history as one of the primary methods of hypothesis testing in extragalactic astronomy (see reviews by Faber 1977; Tinsley 1980; Walcher et al. 2011; Conroy 2013).  \nThere are three key components of SED fitting: a physical model, a set of observations such as photometry or spectroscopy, and a statistical inference framework which connects the model to the observations. The physical models are typically chosen to be stellar population synthesis codes, which combine inputs such as isochrones, stellar spectra, an initial mass function (IMF), a star formation and chemical evolution history, dust attenuation, and (sometimes) dust and nebular emission to create the SEDs of complex stellar populations. There are many publicly available stellar pop-  \nulation synthesis codes which treat these elements with varying degrees of complexity, including GALAXEV (Bruzual & Charlot 2003), P ' EGASE Fioc & RoccaVolmerange (1997); Le Borgne et al. (2004), Maraston (2005), FSPS (Conroy et al. 2009), and BPASS (Eldridge et al. 2017). In some cases, non-stellar effects such as nebular and dust emission can be included separately from the core stellar populations synthesis codes.  \nThe statistical inference framework performs the heavy task of matching complex physical models of galaxy emission to observations. This is usually accomplished through a likelihood function. Many techniques for SED-fitting in the literature thus far have been based on maximum-likelihood optimization, a process sometimes called inversion (Walcher et al. 2011). The optimization may be through non-linear or, for restricted models, linear methods, or a mix of both (Heavens et al. 2000; Cid Fernandes et al. 2005; Ocvirk et al. 2006; Tojeiro et al. 2007; Koleva et al. 2008). These techniques are popular in part because they are fast and simple. However, many of the likelihood spaces in SED-fitting are both highly non-Gaussian and ill-conditioned such that a small change in the input (e.g., noise in galaxy photometry) can lead to a large change in the output (e.g., star formation histories) (Ocvirk et al. 2006). In  \nthis case the pure maximum-likelihood estimation and associated inversion techniques cause large noise amplification, resulting in unstable maximum-likelihood solutions and causing severe difficulties in accurately assessing the uncertainties. While regularization can mitigate noise amplification, it complicates interpretation. Furthermore, it can be difficult to incorporate important nonlinear model components such as dust attenuation and spectral calibration in a maximum-likelihood framework.  \nA solution to these problems can be found in Bayesian forward-modeling techniques. These techniques are able to accurately extract complicated, large, and correlated parameter uncertainties from galaxy observations. Early Bayesian codes such as Kauffmann et al. 2003; Salim et al. 2007, CIGALE (Burgarella et al. 2005), and MAGPHYS (da Cunha et al. 2008), were successful in this effort. These codes often worked by comparing precomputed parameter grids or 'libraries' of SED models to the observations. While likelihood evaluation over a grid can be computationally very fast, grids grow in size exponentially with the number of parameters. As a result grids only allow a relatively small number of parameters or coarse spacing. Furthermore, grids and libraries are large on disk and slow to regenerate, making it difficult in practice to test the effect of the grid assumptions on the inferred parameters.  \nThe success of these early Bayesian forward-modeling codes motivated an expansion of the technique into gridless 'on-the-fly' modeling coupled with Markov chain Monte Carlo (MCMC) algorithms to efficiently explore the resulting parameter space. This approach allows priors and fit parameters to be changed quickly on an object-by-object basis, at the cost of longer per-object computational times. This approach was pioneered by GalMC (Acquaviva et al. 2011), followed later by codes such as BEAGLE (Chevallard & Charlot 2016) and BAGPIPES (Carnall et al. 2018). Here we present Prospector , a modular galaxy stellar populations inference code based on Bayesian forward modeling and Monte Carlo sampling of the parameter space using FSPS stellar populations. Prospector has been in development since 2014 and, in addition to extensive testing on nearby galaxies (Leja et al. 2017, 2018), has seen early use in modeling the SEDs of dwarf galaxies (Pandya et al. 2018; Greco et al. 2018), modeling spatially resolved spectroscopy (Patr'ıcio et al. 2018; Oyarz'un et al. 2019), in modeling the host galaxies of supernovae and gravitational wave events (Blanchard et al. 2017a,b; Nicholl et al. 2017), and in modeling gravitationally lensed high-redshift quiescent galaxies (Ebeling et al. 2018).  \nThe paper is structured as follows. In § 2 we describe the ingredients of the Prospector code including the physical galaxy model, modelling of observational effects, and the statistical framework. We also discuss some general considerations when building models from  \nthese ingredients. In § 3 we present several demonstrations of Prospector fits to mock data, including both photometry and spectroscopy. In § 4 we give examples of Prospector operating on real spectroscopic and photometric data. In § 5 we discuss lessons learned during the development of Prospector , and several avenues for future improvement. Finally, the appendices provide implementation details. All fits and figures in this paper were made using Prospector version 1.0.0.", '2. Prospector INGREDIENTS': 'The Prospector code works by forward modeling the observed spectrum and photometry of complex stellar populations given a combination of parameters describing that population and instrumental parameters, computing a likelihood and posterior probability for the model based on a noise model and prior distributions, and Monte Carlo sampling the posterior probability distribution. In this section we describe, in turn, each of the ingredients necessary for this process, ending with some advice for building models from these ingredients. Readers interested in demonstrations of the code applied to mock and real data can skip to the next sections.', '2.1. Stellar Population Sources': "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 restframe spectrum.  \nIn 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 1 (Foreman-Mackey et al. 2014) bindings. Given a set of stellar population parameters python-FSPS is used to calculate an spectrum integrated over stellar age and metallicity. If necessary, new simple stellar populations (SSPs) may be generated by FSPS during this computation; otherwise the SSPs are cached by pythonFSPS 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 .  \nThe 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 IGM attenuation. More details about the ingredients of FSPS , especially those that have been added or mod-  \nified 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 § 2.2, for several parameters available in FSPS - including the redshift and spectral smoothing - the computations are handled by Prospector externally to FSPS .  \nProspector has access to all of the parametric star formation histories (SFHs) built into FSPS . These forms include exponentially declining SFRs, exponentially declining SFRs with a linear rise (delayedexponential), constant, single age bursts, and a flexible 4-parameter delayed-exponential model which 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 ( § 2.5).  \nThrough 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 so-called 'nonparametric' SFHs, high dimensional piece-wise constant (or step-function) star formation histories, 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. piece-wise linear, with parameters giving the SFR at particular times) may be included in the future. Non-parametric SFHs are highly flexible and capable of modeling star formation histories with sharp features or unusual shapes which are forbidden by the usual parametric forms. It is difficult to constrain flexible star formation histories with typical galaxy observations, so the chosen prior is very important ( § 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.  \nFinally, the modular nature of Prospector allows for sources other than the ones available through FSPS  \nto 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.", '2.2. Redshift, Filters, LSFs': "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 2 (Johnson 2019) where the addition of user-defined filters is straightforward.  \nTo model observed spectra the normalized, redshifted spectrum must be broadened to match the observed linespread 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 the empirical MILES library, with FWHM ≈ 2 . 5 ˚ A from 3525-7400 ˚ A. 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 than the stellar continuum. Prospector has this functionality: see Appendix E for details.  \nAdjustments to the wavelength calibration can be modeled using polynomials, following Koleva et al. (2009). Additional modeling of sky emission and transmission is under development.", '2.3. Spectrophotometric Calibration': "The treatment of spectro-photometric calibration and the uncertainty thereon is of critical importance when modeling spectra, and especially when combining spectroscopic and photometric information ( § 2.4, § 3.1.2, Appendix D.) It is common when modeling spectroscopy to  \nremove the sensitivity of the inferred physical parameters to the continuum shape, due to the difficulty of accurate spectro-photometry (Koleva et al. 2009; Cappellari 2017, e.g.). 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 not 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.  \nThe 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 co-variant or correlated noise, which is discussed further in Appendix D.  \nThe 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 ∼ ( λ max -λ 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 con-  \nined by the shape and normalization of the continuum.", '2.4.1. Likelihood': 'The basic log-likelihood calculation in Prospector is effectively a χ 2 calculation for both the spectral and the photometric data. The two log-likelihoods are then summed. This gives the following expression for the lnlikelihood:  \nln L ( d s , d p | θ, φ, α, β ) = ln L ( d s | θ, φ, α ) (1) + ln L ( d p | θ, β ) ,  \nwhere 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.  \nBy default we assume known and independent Gaussian flux uncertainties. Under this assumption the lnlikelihood for the spectroscopic and photometric terms is simple χ 2 . 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': "Before each likelihood call, the prior probability of every free parameter is computed. A prior probability distribution must be specified for each of the parameters that is allowed to vary during a fit. Several simple distributions are provided including uniform, normal or Gaussian, clipped normal, log-normal, log-uniform, Student's-t, and Beta distributions. The parameters describing these distributions - such as the minimum and maximum for the uniform distribution, or the mean and dispersion for the normal distribution - must be specified as well. Users may provide their own prior distributions. It is also possible in principle to specify joint priors on several parameters (e.g. that one parameter be less than another parameter), though in practice this is more easily achieved through parameter transformations (see below).  \nIn the often highly degenerate space of SED modeling, priors can play a strong role in determining the shape of the posterior probability distribution (e.g. Carnall et al. 2019a; Leja et al. 2019). We therefore strongly encourage exploration of the sensitivity of the results from Prospector to choices of prior, as discussed further in § 2.5.  \nRelated to priors, Prospector offers methods for transforming variables. This works by labeling some parameters as dependent on other parameters via a given  \ntransformation function. This mechanism allows one to sample in transforms of the native parameters (e.g. in τ -1 instead of τ ) that may be more convenient, to tie parameters to other, varying, parameters (e.g. to tie the gas phase metallicity to the stellar metallicity), and to avoid joint priors such as limits on population age that are a function of the varying redshift of the model.", '2.4.3. Optimization and Sampling': "The posterior probability distribution is sampled using Monte Carlo techniques, and optimization can also be performed (though is not recommended for the reasons described in the Introduction, except as an initial guess for sampling algorithms.)  \nThe primary inference technique used by Prospector is nested sampling. Nested sampling (Skilling 2004) is an algorithm based on successive draws from the prior distribution at increasing values of the likelihood. It uses the change in the effective prior volume of iso-likelihood contours as a function of likelihood to estimate the Bayesian evidence, and the successive draws - along with their associated weights - can be used to estimate the posterior distribution of the parameters. Because of the way the prior volume is sampled, nested sampling is well-suited to posteriors that are multi-modal (e.g. for photometric redshifts). Furthermore, no optimization is required before sampling begins, and stopping criteria based on estimates of the remaining evidence can be defined. Nested sampling in Prospector uses the dynesty pure-Python code package (Speagle 2020), an implementation of algorithms described in Feroz et al. (2009); Higson et al. (2019).  \nWithin Prospector it is also possible to use the ensemble MCMC algorithm of Goodman & Weare (2010), as implemented in the emcee Python package (ForemanMackey et al. 2013). This algorithm is based on the use of many 'walkers' in parameter space, the distribution of which is used to construct parameter proposals for the next Monte Carlo step in the chain.  \nThe results of the Monte Carlo sampling are a list of parameter values, which represent (weighted) samples from the joint posterior probability distribution for all the parameters, given the data. These samples provide a natural and useful tool for quantifying parameter uncertainties (marginalized over other parameters of the model) and the degeneracies between parameters. Given samples from the posterior probability distribution, it is up to the user to decide what to report. This may include full set of samples or various projections of the posterior PDF. It is also possible to report moments or quantiles of the marginalized posterior PDF for each parameter, though care must be taken in the interpretation when the posteriors are very broad and/or nonGaussian. For complicated posterior distributions the set of marginalized one dimensional posterior medians may describe a model that is itself very unlikely. For these reasons it is often best to work with the entire  \nposterior distribution (as approximated by the Monte Carlo samples) for subsequent analyses. For visualizing the model in the space of the data the highest posterior probability sample 3 is also useful.", '2.5. Model Building': '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 § 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.  \nWhen 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 un-modeled 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 mis-specified. Finally, with a basic model in hand, one can explore letting parameters that were previously fixed be free, though with informative priors.  \nThere 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 over-interpreting 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.  \nFigure 1. An example of simple parameter inference from a mock UV through Near-IR photometric SED with Prospector . Lower left panels: Several projections of the posterior PDF for all 5 fitted parameters of a simple delayed exponential SFH model (blue contours and histograms) compared to the true input values (black points, black dashed lines). The one-dimensional projections also show the prior probability distributions (green dotted lines). Inset Top: Mock SED for a delayed exponential SFH. We show the true input mock spectrum, smoothed to R ∼ 500 (black line), the mock photometry with S/N= 20 (orange circles), and the posterior PDFs for the true photometry (blue rectangles). The height of each rectangle represents the 16th-84th percentile range of the posterior PDF for the true photometry. Inset Bottom: Uncertainty normalized residual ( χ ) between the mock photometry and the predicted photometry of the most probable Monte Carlo sample (orange).  \n<!-- image -->  \nFurthermore, 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  \nadding parameters is that MCMC sampling can become inefficient as the number of parameters increases.  \nWhen 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 informative priors on some parameters, by which we mean here priors that are not uniform. 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  \nto 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.  \nHowever, 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.', '3. DEMONSTRATIONS WITH MOCK DATA': "In this section we provide demonstrations of some of the capabilities of the Prospector code for fitting galaxy spectra and/or photometry using mock data. These demonstrations include the combination of photometric and spectroscopic information in a single fit, the inference of flexible, many-parameter SFHs, and the constraints provided on a complex model with different combinations of photometric data, including a single optical photometric point.  \nThroughout these demonstrations with mock data we use the MIST isochrones (Choi et al. 2016) and the MILES spectral library (S'anchez-Bl'azquez et al. 2006) within FSPS to construct both the mocks and the models with which they are fitted. Because we are fitting mock data with the same stellar population models used to generate that data, the results of this section allow us to demonstrate and test the Prospector methodology and examine the constraints that can placed in the limit of perfect stellar population models.  \nThe demonstrations in this section are not intended as exhaustive explorations of the information content of different kinds of data and the suitability of different kinds of models, but rather provide examples of the capabilities of the Prospector code for fitting different kinds or combinations of data with a variety of galaxy SED models.", '3.1. A Simple Parametric SFH 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 star formation history 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  \nhave 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.  \nHere 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 glyph[lscript] is given by  \nψ ( t glyph[lscript] ) ∝ t age -t glyph[lscript] τ e -( t age -t glyph[lscript] ) /τ , 0 < t glyph[lscript] < t age (2)  \nwith 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 ˚ A - 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 3-parameter models of Draine & Li (2007) (see Appendix A for details), but we do not fit for these parameters, as we do not include sufficiently red data to be 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.  \nThe 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.', '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  \nFigure 2. Comparison of parameter inference from photometry, continuum normalized spectroscopy, and the combination. Left panels: Several projections of the posterior PDF for selected parameters of a parametric SFH model. The parameters inferred from photometry only (blue shaded regions), continuum normalized spectroscopy only (red shaded regions), and photometry plus continuum normalized spectroscopy (orange shaded regions) are compared to the true input parameters (black points, black dashed lines). The prior probability distributions, scaled by an arbitrary factor, are shown in the one-dimensional projections (green dotted lines). Top right panel: The true spectrum, roughly continuum normalized as described in the text, is shown (black) along with the noisy mock spectrum (gray)  \n<!-- image -->  \nGALEX 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.  \nFor 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 free parameter, but with a very narrow prior distribution as might be expected from low S/N grism spectra, yielding a model with 7 free parameters. We then  \nuse 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.  \nTable 1. Summary of parameters and priors used in mock demonstrations  \nNote -Uniform( x , y ) indicates a uniform prior from x to y , while N ( µ, σ ) indicates a normal distribution for the prior and LogUniform( x , y ) indicates distribution that is uniform in the log of the parameter in the parameter range ( x , y )  \nExamining 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.  \nIn Figure 1 we also show projections of the posterior probability distribution for several of the free parameters (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 onedimensional 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-agemetallicity degeneracy is mapped in complete detail.  \nIt 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 = 200kms -1 and assume an instrumental dispersion of 2 ˚ A per pixel from rest-frame 3800 ˚ A-7100 ˚ A (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 line-spread function as the MILES library redshifted to z = 0 . 1 ( ∼ 2 . 75 ˚ A FWHM). See Appendix C for approaches to data where the instrumental or library broadening is a significant contributor to the final linespread function 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 6th-order polynomial fitted to regions of the spectrum free from strong absorption or emission lines.  \nTo 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 12thorder 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 also being fit.  \nIn 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 since all 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.  \nThe 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.  \nFigure 2 shows that broadband UV through NIR photometry and continuum normalized spectroscopy carry different information about 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 § 2.4 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 be likely occur only if the SFH and dust model where very flexible, or the spectra were of much lower S/N than the photometry.  \nIn 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 log(Z / Z glyph[circledot] ), which are partially degenerate with dust attenuation in the photometry alone, is increased by inclusion of the spectroscopic data.', '3.2. Inferring the SFHs of Illustris galaxies': "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 (Gallagher et al. 1984, e.g.). 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  \nFigure 3. Posterior probabilities for SFHs inferred from S/N ∼ 100 optical spectra with no nebular emission and perfectly known spectrophotometric calibration. The mock data is generated from an Illustris SFH (black). The shaded regions show the 16th-84th percentile credible interval for the inferred parametric SFH (blue) and non-parametric SFH (maroon). The non-parametric SFH uses 14 bins with a continuity prior and the average input SFR in the same bins is shown as dotted gray lines. The insets in each panel show the recent SFH ( < 1 Gyr) in more detail.  \n<!-- image -->  \nuncertainty on the SFH (Carnall et al. 2019a; Leja et al. 2019).  \nOne 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 semi-analytic 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.  \nAnother possibility is to use piece-wise constant or δ -function SFHs where the mass in each temporal bin be-  \nomes 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 piece-wise constant basis has the advantage of being readily interpreted.  \nProspector supports non-parametric models using a piece-wise constant basis ( § 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:  \nlog r j =log( s j /s j +1 ) (3) s 0 =1 s i = s 0 ∏ j<i j =0 r -1 j m i = s i ∆ t i M i = m i 10 logmass ∑ i m i ,  \nwhere ∆ t i is the width of of the i th bin of lookback time, in years. We place a Student's t -distribution prior on each of the log r j parameters, centered on 0 (i.e. constant SFR); the Student 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.  \nTo demonstrate the capabilities of Prospector to infer both parametric and non-parametric 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 star formation history. This selection provides the greatest challenge for the inference machinery, and are also weighted against 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  \nspectrum. 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 = 100kms -1 . As for the mocks in § 3.1 we assume an observed frame instrumental linespread function equivalent to the MILES spectral library at the assumed redshift. We assume a 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 non-parametric 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.  \nFor the parametric SFH we use the same delayed exponential model as in § 3.1. For the non-parametric SFH model we adopt 14 bins approximately evenly spaced in log t glyph[lscript] . 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.  \nIn Figure 3 we show the input Illustris SFHs and the posterior PDFs for the SFHs inferred with both a parametric and non-parametric 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 non-parametric 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 is remains evident, though at a lower level - the recovered non-parametric SFH will tend towards 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 non-parametric 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).", '3.3. A Complex Model: The Effect of More Data': '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 com-  \nFigure 4. Change in posterior PDF for the SED and parameters of a complex model as the number of photometric bands increases. The model includes a continuity SFH with 6 temporal bins, a complex dust attenuation model, dust emission parameters, and AGN contribution to the MIR. There are 16 parameters total. Nebular emission is also included. The model is fit to a varying number of mock photometric data points, indicating how the model becomes more tightly constrained with the addition of multi-wavelength data (animation). Left panel: The posterior prediction for the SED (blue shaded region, indicating the 16th-84th percentile range for the flux at every wavelength), resulting from a fit to photometry in the indicated bands (orange circles) and compared to the true input SED (black line). Right panels: The marginalized posterior (blue shaded regions) and prior (dotted green lines) probability distributions for selected parameters of the model, compared to the true input parameters (dashed black lines). Several of the parameters shown are derived from the fundamental fitted model parameters, see text for details. We show parameters of the stellar population (top row), the dust attenuation model (second row), the dust emission model (3rd row) and the AGN parameters (bottom row).  \nplex 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.  \nOne 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, since models are  \nonly generated in regions of high probability as needed, instead of over the whole parameter space.  \nThe model we explore is described by a nonparametric SFH with 6 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 towards 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 2 parameters (Nenkova et al. 2008a; Leja et al. 2018, A) and dust emission described by the 3 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  \nTable 2. Summary of parameters and priors used in the complex model  \nNote -Uniform( x , y ) indicates a uniform distribution from x to y , while N c ( µ, σ, 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 ).  \nprior distributions favor typical dust attenuation curves and disfavor large AGN contributions to the mid-IR.  \nWe generate broadband photometry in standard UVFIR 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.  \nIn 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 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 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 very high dimensional model.  \nAs 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.', '4. DEMONSTRATIONS WITH REAL DATA': 'We now turn to demonstrations of Prospector usage with real data.', '4.1. Globular Cluster Spectroscopy and Photometry': '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 simple stellar populations. 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 (BHB), 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 clus-  \nFigure 5. Example of a fit to globular cluster data, showing results for NGC1851. Top left panel: Comparison of the observed photometry (white circles and error bars) to the posterior probability distributions for the true photometry (blue rectangles; height gives the 16th-84th percentile range) and the true intrinsic spectrum (shaded red, showing 16th-84th percentile range). We also show the intrinsic spectrum for the most probable posterior sample (thin black line). All spectra are smoothed to R ∼ 500 for display purposes. Middle left panel: Uncertainty normalized residual between the observed spectrum and the model spectrum of the most probable posterior sample (thin red line). Note that these are normalized by the nominal uncertainties (median S/N= 320), and do not include the fitted noise inflation term. The uncertainty normalized residuals for the photometry are also shown (white circles). Bottom left panel: The inferred spectrophotometric calibration vector (red shaded region, 16th84th percentile), defined as the ratio between the observed spectrum and the true intrinsic spectrum. Right panels: Projection of the posterior PDF for selected parameters of the model, inferred from the combination of photometry and spectroscopy. The priors for the parameters shown are uniform.  \n<!-- image -->  \nter population, then they may be less useful as tests of stellar population models.  \nHere 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 longslit spectroscopy of 40 globular clusters. These spectra span the wavelength range 3400-6400 ˚ A, with a spectral resolution that is wavelength dependent but averages 3.2 ˚ A FWHM. The nominal S/N is often greater than 100. We take UBVRI photometry from Harris (1996), who provide total (integrated) magnitudes of clusters. We assign uncertainties of 10% to the photometric fluxes.  \nTo 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 log(Z / Z glyph[circledot] ) (assuming scaled solar abundances). We fix the distance to the value given by Harris (1996) since for integrated SEDs at non-cosmological distances the distance and mass are perfectly degenerate. We broaden the model library spectra to match the instrumental resolution as  \nreported 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.  \nIn 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 ∼ 2.  \nIn Figure 6 we show a comparison of the parameters that we infer for many GCs, compared to literature values as compiled by Koleva et al. (2008) and to values inferred from the same spectra using the UlySS code (Kol-  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 6. Comparison of globular cluster parameters inferred with Prospector to literature values and to values inferred with ULySS by Koleva et al. (2008). Literature ages are largely from CMD fits, as tabulated in Koleva et al. (2008). Clusters without reliable age determinations were assigned an age of 10 Gyr. Literature metallicities are from high resolution studies of individual stars, again as tabulated in Koleva et al. (2008). Literature extinction is computed from the E(B-V) values in Harris (1996). The color coding is by a measure of horizontal branch color (HBR) in each cluster. The bluest or hottest horizontal branches have HBR ∼ +1. Prior limits on metallicity and age are indicated as gray dotted lines. Top Row: Stellar metallicity. Middle Row: Globular cluster age. Bottom panel: Dust extinction.  \n<!-- image -->  \neva 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 un-modeled α -element enhancements. The median offset is 0.04 dex and the r.m.s. 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 ).  \nPrevious 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, since they are not included in the modeling, will tend to result 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 un-modeled α enhancement or other limitations in the isochrones. We note that differences between 10 and 14 Gyr are quite subtle in integrated light.  \nOverall, we find that Prospector performs as well as other stellar population fitting programs when applied to globular clusters.', '4.2. Photometric Redshift: GNz-11': '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 GNz11. 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 followup grism spectroscopy with HST confirms a redshift z = 11 . 09 (Oesch et al. 2016).  \nThe 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 restframe λ < 1216 ˚ A 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 5 temporal bins, with the continuity prior discussed in § 3.2. The exact width and extent of these bins are dependent on the redshift  \nFigure 7. Photometric redshift inference with Prospector . It is possible to include the redshift as free parameter in the model. Here we show the results of inference from the photometry of GNz-11. Top row: The fitted observed photometry (orange with black error bars) from Oesch et al. (2016) as well as the posterior PDFs for the true photometry (blue rectangles) and the predicted spectrum of the most probable sample from the Monte Carlo chain (red) which has a redshift of 10.35. Second row: Uncertainty normalized residual ( χ ) between the observed photometry and the predicted photometry for the most probably Monte Carlo sample. Third row: The posterior PDF for the redshift (blue shaded region). The prior is uniform over the range 2 < z < 13. Also indicated is the redshift inferred by (Oesch et al. 2016) from HST grism data (black vertical dashed line). The text gives the median of the redshift PDF, with uncertainties derived from the 16th and 84th percentiles. Fourth row: Marginalized posterior PDFs for selected model parameters (blue histograms) as well as the prior probabilities for each parameter (green dotted lines) Fifth row: Inferred SFH. The SFH is modeled with 5 bins in lookback time, which change in width depending on the model redshift.  \n<!-- image -->  \nof 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 dispersion of 0.3. The ideal redshift prior would be the product of the redshiftdependent differential cosmological volume element and the integral of the redshift-dependent luminosity function above some limiting magnitude. However, for clarity and due to a 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.  \nThe 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 . 311 . 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 glyph[circledot] , low dust attenuation, a weak constraint on the stellar metallicity, and no constraint on the IGM opacity scaling factor beyond the imposed prior distribution.', '4.3. SDSS Post-Starburst Spectrum': "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 poststarburst galaxy to demonstrate both the flexible SFH model and the flexible emission line modeling available within Prospector .  \nWe fit SDSS J133751.29+354755.7 from the poststarburst catalog of Goto (2007), specifically targeted for its high signal-to-noise 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.  \nA 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  \nis 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 non-parametric 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.  \nThe results of this fit are shown in Figure 8. Prospector is able to fit the data with high precision; the bestfit model produces excellent χ 2 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 are assigned any significant weight. As expected, a significant post-starburst event is detected within the previous ∼ 300 Myr, suggesting a factor-oftwenty 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.  \nThe spectral fit includes a zoom-in on the H α and [N II] doublet. The [NII]/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, AGN, 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 [NII]/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.", '5. DISCUSSION': "5.1. Lessons Learned and Practical Advice  \n<!-- image -->  \nFigure 8. A simultaneous Prospector fit to the spectroscopy and photometry of an SDSS post-starburst galaxy (Goto 2007). Top row: The observed ugriz photometry (black), compared to the model posterior prediction for the photometry (red) and the model spectrum of the most probable sample. Second row: The observed SDSS spectrum (black), compared to the model posterior prediction for the spectrum (red). The inset panel zooms in on the H α and [NII] emission lines, showing the continuum model (grey), emission lines predicted from the CLOUDY grid (blue), and emission lines from the nebular marginalization model (red). The lower panel shows the normalization vector between the spectroscopy and photometry. Third row: Marginalized 1D posterior distributions for selected model parameters, along with their median and 1 σ confidence intervals. Fourth row: The recovered star formation history. The factor-of-ten decrease in the star formation rate at ∼ 300 Myr is the post-starburst event.  \n<!-- image -->  \nThe development of Prospector began in 2014, and it has undergone continual development, refinement, and refactoring since. 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 makes 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.  \nPriors 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.  \nPriors 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 through 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.  \nWhile 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 populations 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.", '5.2. Future directions': '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 populations 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 pre-computed (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.  \nWhile 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 HII regions, radio synchrotron emission, and X-ray through optical emission from AGN. Additionally, Lymanα 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.  \nWe 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 selfconsistent nebular emission models may prove useful.  \nPerforming 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 highdimensional 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 ex-  \nample corresponding to the correlation function between SFRs in different time periods (e.g. Kelson et al. 2016; Tacchella et al. 2020) may be implemented.  \nMultiple spectra described by the same physical parameters but substantially different instrumental parameters (e.g. spectral resolution or calibration) can not 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.  \nIn 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 remaining challenges. These efforts are necessary to fully exploit the information content of rich datasets that will be provided by future facilities, such as James Webb Space Telescope .', 'ACKNOWLEDGMENTS': "We are grateful for discussions early in the project in with Dan Weisz, Dan Foreman-Mackey, and David Hogg, as well as conversations with Adam Carnall and Sandro Tacchella. We are also grateful for early testing of the code and suggestions by John Moustakas, Antara Basu-Zych, Dylan Nelson, Imad Pasha, Tom Zick, Song Huang, and Johnny Greco. BDJ and CC acknowledge support from the Packard Foundation and NSF grants AST-1313280 and AST-1524161. This research made extensive use of NASA's Astrophysics Data System Bibliographic Services. Computations in this paper were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University.  \nFunding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.  \nSoftware: dynesty (v1.0.0 Speagle 2020), emcee (v2.2.1 Foreman-Mackey et al. 2013), FSPS (v3.1 Conroy et al. 2009, 2010), python-FSPS (v0.3.0 Foreman-Mackey et al. 2014), Astropy (v4.0 Astropy Collaboration et al. 2013, 2018), NumPy (v1.17.2 Harris et al. 2020), SciPy (v1.3.1 Virtanen et al. 2020), matplotlib (v3.1.1 Hunter 2007), IPython (Perez & Granger 2007), f2py (Peterson 2009)", 'A. FSPS UPDATES': "There have been many updates and improvements to the FSPS code since 2010. In this section we provide a brief overview of the main features now available in the code.  \nIsochrone libraries are available from the MIST (Choi et al. 2016), PARSEC (Bressan et al. 2012), BaSTI (Cordier et al. 2007), Padova (Marigo et al. 2008), Geneva (Meynet et al. 1994), and BPASS (Eldridge et al. 2017) models. The BPASS models are provided as simple stellar populations and are therefore treated differently than the other isochrone sets: for example, the isochrone-based parameters that affect the strength of the TP-AGB cannot be applied to the BPASS models, nor can the IMF be changed.  \nAs in Conroy & Gunn (2010), the available base spectral libraries include BaSeL and MILES. Those libraries are supplemented with the TP-AGB library from Lan¸con & Wood (2000), the post-AGB library from Rauch (2003), and the WMBASIC hot star library from Eldridge et al. (2017).  \nThere are a variety of options for the dust absorption model including the Milky Way extinction curve of Cardelli et al. (1989) with variable slope and UV bump strength, empirical attenuation curves from Calzetti et al. (2000) and Kriek & Conroy (2013), the two-component dust model of Charlot & Fall (2000) with separate normalization of each component and adjustable attenuation curve slopes, and the radiative transfer models from Witt & Gordon (2000). 'Diffuse' dust emission is provided by the Draine & Li (2007) templates. Circumstellar dust around AGB stars is implemented via the Villaume et al. (2015) templates, while dust associated with AGN is modeled using a simplified implementation of the Nenkova et al. (2008b) templates.  \nNebular emission associated with HII regions is included using Cloudy lookup tables computed by Byler et al. (2017). Attenuation by the intergalactic medium (IGM) follows the prescription described in Madau (1995), with an optional overall scaling factor to account for line-of-sight variations.", 'B. SFH INTEGRALS AND PARAMETRIC FORMS': 'The fundamental equation for computing the spectra of composite stellar populations is  \nF λ = ∫ t univ 0 dt ∫ Z max Z min dZ ψ ( t, Z ) s λ ( t, Z ) e -τ λ ( t,Z ) (B1)  \nwhere t is the lookback time , ψ ( t, Z ) is the SFH, s λ ( t, Z ) is the spectrum of a unit mass of stars of age t and metallicity Z , and τ λ ( t, Z ) is the effective dust opacity towards stars of age t and metallicity Z .  \nIn practice, this is approximated by a sum over discrete simple stellar populations (SSPs). These SSPs are the s λ ( t, Z ) for J specific values of t (indexed as t j ) and K specific values Z (indexed as Z k ). They are calculated using isochrones and stellar spectral libraries. The sum is then  \nF λ = J ∑ j =1 K ∑ k =1 m j,k s λ,j,k e -τ λ,j,k , (B2)  \ns λ,j,k = s λ ( t j , Z k ) , (B3)  \nwhere the m j,k are the formed masses or weights of each SSP.  \nIt is tempting to use trapezoidal integration of the SFR at t j , Z k to calculate these weights. However, if the SFR has a discontinuous derivative or is highlynonlinear between the SSP points, this can lead to substantial errors. Increasing the effective temporal resolution of the SSPs via interpolation can mitigate (but not solve) this problem, at the expense of the increased memory needed to store the interpolated SSPs, and the increased computational time needed to sum them with the appropriate weights.  \nA more robust and scalable solution is to integrate the product of the interpolation weights and the SFH. That is, if a spectrum at time t j < t < t j +1 can be written as  \ns λ ( t ) = a j ( t ) s λ ( t j ) + b j +1 ( t ) s λ ( t j +1 ) , (B4) a j ( t ) = log( t j +1 ) -log( t ) log( t j +1 ) -log( t j ) , b j +1 ( t ) = log t -log( t j ) log( t j +1 ) -log( t j ) ,  \nwhere we have chosen a linear interpolation in log t (because changes in the SSP fluxes are generally linear in log t , not in linear time), then we wish to calculate  \nm j = ∫ t j +1 t j dt Ψ( t ) a j ( t ) + ∫ t j t j -1 dt Ψ( t ) b j ( t ) . (B5)  \nThis gives the contribution to the weights of the SFR in the bins that are older and younger than t j , and for parametric SFHs this can be solved analytically for each j . This is what is used in the FSPS code for the available  \nparametric SFHs, as of v3.0. Combinations of parameterized SFHs can be easily considered by calculating the above quantity for each component and summing the m j values with the appropriate scaling to retain the proper mass fractions in each component. Similarly, for tabular SFHs we make the assumption that the SFR is a linear function between the tabulated time points, and solve the equation above for each t j and segment separately, accumulating the weights m j from each segment.', 'C. SMOOTHING': "When comparing model galaxy spectra to observed spectra there are (at least) three spectral smoothing resolutions (or smoothing kernels) that must be considered. These are: 1) the line-of-sight velocity distribution (LOSVD) of stars in the restframe of the galaxy (L g ); 2) the observed frame line-spread function of the instrument (L i ); and 3) the restframe line-spread function of the model spectral library (L glyph[lscript] ) . Along with any imposed smoothing of the model spectrum, these combine to create the line-spread functions of the data (L d ) and the model (L m ) as follows:  \nL d ( λ o ) = L i ( λ o ) ⊗ L g ( λ r ) , L m ( λ o ) = L glyph[lscript] ( λ r ) ⊗ L k ( λ o ) , λ o =(1 + z ) λ r , (C6)  \nwhere λ o is the observed frame wavelength, λ r is the restframe wavelength, L k is the kernel corresponding to extra smoothing applied while modeling the data, and ⊗ denotes convolution. While the LOSVD of the galaxy can be handled by straightforward and fast smoothing in velocity space (e.g. Cappellari 2017), the instrument and library resolution are both potentially complicated functions of observed-frame and rest-frame wavelength respectively. This can be difficult to implement in a speedy convolution algorithm.  \nWhen L i = L glyph[lscript] or both are narrow compared to L g they can be neglected and the best match between the data and model is obtained when L k ∼ L g . Varying L k to determine the best match can be done using fast FFT based smoothing in velocity space. However, when the LOSVD of the galaxy is comparable to or only slightly broader than the instrumental or library LSF (at any observed wavelength) then these latter effects must be taken into account. The possible wavelength dependence of the instrumental and library resolutions makes this difficult, as standard FFT-based algorithms do not admit smoothing kernels that are wavelength dependent. The usual solution is to use slower brute-force convolution algorithms to smooth the library spectra once , before any other calculations are made or any fitting is done, such that L ' glyph[lscript] ( λ r ) ∼ L i ( λ o ). This is possible provided the library everywhere has higher resolution than the instrument, that L i is well known, and that the redshift of the model does not vary significantly.  \nIf any of these conditions are not met - for example if one wishes to infer and marginalize over uncertainties in the instrumental LSF - then another approach must be used. A solution is to include, in each generation of the model, a convolution by the kernel L k defined by the difference between the rest-frame library resolution and the observed frame instrumental and physical resolution. This potentially wavelength dependent difference kernel may depend on the modeled systemic velocity. Because the convolution is done at each model generation step, this scheme also allows for simultaneous fitting of the LSF.  \nIn Prospector we allow for both options: either smoothing the library once, before fitting, by the difference between the observed-frame instrumental and rest-frame library resolution or smoothing at each model generation step by this difference, calculated on the fly for the model systemic velocity, velocity dispersion, and instrumental LSF parameters. We assume that both the instrumental and library resolution can be approximated by Gaussians at each wavelength, so that the difference kernel is also represented by a wavelength dependent Gaussian.  \nIt is desirable to have algorithms that can quickly and accurately produce spectra smoothed by a wavelength dependent kernel, so that instrumental smoothing can be done during the fitting process in a flexible way. Since the difficulty with FFT-based algorithms in this case is that the kernel is varying as a function of wavelength, the solution is to resample the model spectrum onto a space in which the kernel is not varying (analogous to resampling a spectrum from a regular wavelength grid to a regular velocity grid when smoothing by a LOSVD.)  \nAssuming the line-spread function can be described by a Gaussian, this can be accomplished by defining a coordinate system from the cumulative distribution function of the inverse of the dispersion as a function of wavelength  \nx ( λ ) ≡ K ∫ λ λ min dλ ' σ ( λ ' ) , (C7)  \nwhere λ is the wavelength, σ ( λ ) is the width of the Gaussian LSF at any wavelength, and K is a normalization constant. A regular sampling in x will provide an equal number of samples per σ , even as σ changes with wavelength. Thus, convolving a spectrum resampled onto an even grid in x with a single Gaussian (in x ) will produce a spectrum with a smoothing in λ that is λ dependent.  \nThe integral that defines x can be done analytically or, for more flexibility, numerically. The constant K and the final spacing in x should be chosen so that the interpolated spectrum adequately samples at least the output line-spread function, and preferably the input LSF, for some definition of adequate. Such an algorithm is implemented within Prospector . Further improvements to the smoothing functionality in Prospector may include a larger variety of smoothing kernels (including  \nhigher order LOSVD moments and kernels appropriate for pixelization.)", 'D. ALTERNATIVE NOISE MODELS': "In addition to the default assumption of known Gaussian independent noise described in § 2.4, we have also implemented more complex likelihood calculations. These allow for properties of the noise itself to be inferred, or for the presence of correlated noise to be modeled. This is accomplished though a flexible noise covariance matrix construction. The ln-likelihood for the more complex noise model is  \nln L ( d | Θ , γ ) = -1 2 [ ∆ T C -1 ( γ ) ∆ + ln || 2 π C ( γ ) || ] (D8) ∆= d -µ (Θ) , (D9)  \nwhere d is the data vector of length n , µ is the mean physical and instrumental model which depends on the parameters Θ, ∆ is the residual, and C is the variancecovariance matrix which depends on the noise model parameters γ .  \nThe parameters Θ generally include the physical parameters θ . They may also include instrumental effects such as wavelength-scale distortions or the instrumental resolution as parameters φ .  \nIn Prospector the covariance matrix is specified as a sum of weighted kernels applied to given metrics , where the kernels, weight vectors, and metrics can be specified at run time. More precisely, we use the following formulation of C  \nC i,j = ∑ glyph[lscript] w i,glyph[lscript] w j,glyph[lscript] K glyph[lscript] ( x i , x j , γ glyph[lscript] ) , (D10)  \nwhere w is the weight vector of length n , x is the metric (also of length n ), K glyph[lscript] is the glyph[lscript] th kernel, and γ glyph[lscript] are the noise parameters for that kernel.  \nThe kernels K can in principle be anything that results in a positive semi-definite matrix when multiplied by the weight vectors w . The following are available by default in Prospector :  \n- · uncorrelated: K ( x i , x j | a ) = aδ i,j . This is a simple diagonal matrix parameterized by an amplitude a for the diagonal elements. All other elements are 0.\n- · photo-cal: K ( x i , x j | a ) = a for i, j corresponding to user-specified photometric bands. This adds correlated noise to these bands with amplitude a . All other matrix elements are 0.\n- · exp-squared: K ( x i , x j | a, L ) = ae -( x i -x j ) 2 / (2 L 2 ) . This is the exponential squared kernel, parameterized by an amplitude a and scale length L .  \nA simple use for this more complicated noise model is to incorporate uncertainty in the noise estimates themselves. This is accomplished by allowing the matrix C  \nto depend on parameters γ that are then inferred and can be marginalized over. For example, to allow for the possibility that the supplied uncertainties are underestimated by a constant but unknown factor b a diagonal kernel K 1 ( x i , x j ) = b 2 δ i,j can be combined with weights w i, 1 given by the nominal supplied uncertainties, where δ i,j is the Kroneckerδ . Increasing b will decrease the first term in the likelihood of Equation D8, but this is eventually balanced by an increase in the normalization given by the determinant of C , preventing the uncertainties from becoming arbitrarily high. Inferring such an inflation in the uncertainties can be useful for providing more accurate final parameter uncertainties. An additive noise component, i.e. a noise floor, can also be incorporated with an additional diagonal kernel K 2 ( x i , x j ) = aδ i,j and weights w i, 2 = 1.  \nThe covariance matrix can be used to account for correlated noise in the data. Correlated noise can be induced during extraction of 2-D spectra, smoothing of the data, or interpolation to a particular wavelength grid (e.g. Bolton & Schlegel 2010). When the pixel covariances are known from error propagation during extraction, they can be supplied in matrix form as the kernel K , with weights w i = 1.  \nSky subtraction residuals are likely to be correlated over at least several pixels and may also be treated as additive covariant noise (e.g. Carnall et al. 2019b). The scale length and amplitude of these correlations may be supplied or inferred.  \nFinally, the use of weight vectors also allows calibration to be modeled as covariant or correlated noise. Here the covariant noise may be due calibration uncertainties or to inaccuracies in the spectral models. Since calibration is multiplicative with respect to flux, modeling calibration noise can be accomplished by using the predicted model flux vector µ (Θ) as the weight vector w . This approach was used for photometry, where calibration uncertainties from different instruments or observatories are often correlated, by Gordon et al. (2014). It can be extended to spectroscopic data. This approach has the benefit of being more flexible than than polynomials or other parameterized functions, which will introduce systematic effects if they are not a good representation of the actual spectrophotometric calibration function. It also has the benefit that a model where the data can be described well purely by parameters unrelated to the calibration (e.g. absorption lines) will be higher likelihood than one where correlated noise is necessary. This mitigates the problem that physical information in the data may be lost to the spectrophotometric calibration model if the latter is overly flexible. However, this approach requires a realistic model for the correlations induced by mis-calibration, which can be complex, multi-scale and must be separated from correlations induced solely by data processing (e.g. by interpolation or spectral extraction.) This approach may also incur high computational cost.  \nFinally, we also include a rudimentary outlier model in Prospector . This model approximately accounts for the presence of data which deviates from the model much more than should be allowed by the noise, due for example to cosmic rays or strong, un-masked sky subtraction residuals. The outliers are modeled as a small fraction of data points that are drawn from a normal distribution centered on the model fluxes but with much larger than the nominal supplied uncertainties. By analytically marginalizing over the exact location of the outliers (e.g. Hogg et al. 2010), this can be expressed as  \nln p ' ( d | Θ , β ) = (1 -f out ) ln p ( d | Θ , β ) (D11) +f out ∑ i d i -µ i (Θ) 2 ( s out σ i ) 2 ,  \nwhere f out is the fraction of outlier data points, and s out is the inflation of the nominal noise for outlier points, which defaults to 50. Both of these parameters may be supplied or inferred. At the present time this outlier model cannot be used with correlated noise models.", 'E. EMISSION LINE MARGINALIZATION': 'As described in § A the emission line model included with FSPS is based on Cloudy (Ferland et al. 2013; Byler et al. 2017). This treatment is sufficient in many contexts, but deviations of emission line ratios and luminosities from the two-parameter model implemented there may occur for a variety of reasons. These include the presence of AGN, shocks, or other non-stellar ionizing sources, abundance variations, and complicated geometries. Furthermore, the kinematics of the ionized gas may be different than that of the stars. It is therefore desirable to include some additional flexibility in the emission line modeling.  \nThis is a challenge; there are 128 emission lines included in FSPS , and it is not practical to sample directly for the amplitude of each line. With some reasonable assumptions, it is possible to introduce analytic Gaussian emission line modeling into the Prospector likelihood which do not require explicitly sampling for emission line amplitudes. When the prior on the amplitudes is also Gaussian, the posterior probability for the amplitude is also Gaussian and can be analytically marginalized in each likelihood call. The emission line width and central wavelength offset is sampled for as well; the offset and width are assumed to be the same for each emission line to ease the sampling, though this assumption can be relaxed by the user if necessary.  \nFor clarity we define here a Gaussian distribution N of the variable x with mean µ and covariance matrix Σ as  \nN ( x | µ , Σ ) = 1 √ || 2 π Σ || e -1 2 ( x -µ ) T Σ -1 ( x -µ ) . (E12)', 'E.1. Analytically Marginalizing over emission lines': '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  \nL ( D | θ ) = ∫ d α L ( D | θ, α ) p ( α ) (E13)  \nwhere 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 loglikelihood, conditional on both θ and α , as  \nln L ( D | θ, α ) = ( ∆ -Gα ) T Σ -1 ( ∆ -Gα ) , (E14) ∆ = { D λ -M λ } , G = { G λ , j } , G λ,j ≡N ( λ | (1 + z ) λ 0 ,j , σ v c (1 + z ) λ 0 ,j )  \nwhere Σ 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 θ ),  \nThe 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, ˆ α . 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 ˆ α and has a width described by the covariance matrix Σ ˆ α Our expression for the likelihood function is then  \nL ( D | θ, α ) = L ( D | θ, ˆ α ) N ( α | ˆ α , Σ ˆ α ) N ( ˆ α | ˆ α , Σ ˆ α ) . (E15)  \nThe denominator is a non-unity value required to make the equality hold for α = ˆ α . When the prior is uniform ( p ( α ) is constant) the integral in the marginalized likelihood of equation E13 - ignoring a constant factor that is independent of θ and α - reduces to  \nL ( D | θ ) = L ( D | θ, ˆ α ) N ( ˆ α | ˆ α , Σ ˆ α ) . (E16)  \nWe 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 ˆ α as  \nˆ α = ( G T Σ -1 G ) -1 G T Σ -1 ∆ . (E17)  \nWe are also interested in the widths or uncertainty on the MLE value, described by the covariance matrix Σ ˆ α .  \nThis is obtained from the inverse second derivative of the likelihood with respect to α , resulting in  \nΣ ˆ α = ( G T Σ -1 G ) -1 . (E18)  \nNotably, 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.  \nUsing 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.', 'E.2. Gaussian priors on the emission line amplitudes': 'To incorporate priors we can multiply the factorized likelihood in equation E15 by the prior. This is straightforward if the prior is another Gaussian with mean ˘ α and covariance Σ ˘ α : that is, p ( α ) = N ( α | ˘ α , Σ ˘ α ). The product of two Gaussians is another Gaussian; using well known formulae for the products of normal distributions we obtain  \nL ( D | θ, α ) p ( α ) = L ( D | θ, ˆ α ) N ( ˆ α | ˘ α , Σ ˆ α +Σ ˘ α ) N ( ˆ α | ˆ α , Σ ˆ α ) N ( α | ¯ α , Σ ¯ α ) (E19)  \nwhere ¯ α and Σ ¯ α are the mean and covariance of the new Gaussian. In fact this new Gaussian does not need to be evaluated, since the integral in the marginalized likelihood of equation E13 reduces to  \nL ( D | θ ) = L ( D | θ, ˆ α ) N ( ˆ α | ˘ α , Σ ˆ α +Σ ˘ α ) N ( ˆ α | ˆ α , Σ ˆ α ) . (E20)  \nWe do not explore non-Gaussian priors, as these are considerably more difficult to implement.', 'E.3. Practical usage': '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 β .  \nWhen applying the Gaussian priors, the center of the prior for each line is taken as the prediction from the  \nCloudy implementation. This means that the lines are assumed to be powered by normal star formation, but  \ndeviations are permitted. The width of the prior is specified by the user, in units of (true/predicted) luminosity.', 'REFERENCES': "Villaume, A., Conroy, C., & Johnson, B. D. 2015, ApJ, 806, 82  \n- Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261  \nVogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518  \n- Walcher, J., Groves, B., Budav'ari, T., & Dale, D. 2011, Ap&SS, 331, 1\n- Weisz, D. R., Dalcanton, J. J., Williams, B. F., et al. 2011, ApJ, 739, 5  \nWilkins, S. M., Gonzalez-Perez, V., Lacey, C. G., & Baugh, C. M. 2012, MNRAS, 427, 1490 Witt, A. N., & Gordon, K. D. 2000, ApJ, 528, 799 Yan, R., Newman, J. A., Faber, S. M., et al. 2006, ApJ, 648, 281  \nYang, Y., Tremonti, C. A., Zabludoff, A. I., & Zaritsky, D. 2006, ApJL, 646, L33"}