Nonparametric Star Formation History Reconstruction with Gaussian Processes. I. Counting Major Episodes of Star Formation

The main thrust of the Dense Basis method is to encode the maximum amount of information about the SFHs of galaxies using a minimal number of parameters. In this respect, the formalism used to describe SFHs in this work can be used as a module in existing sophisticated inference methods developed in public SED fitting codes like Bagpipes (Carnall et al. 2018), Beagle (Chevallard & Charlot 2016), CIGALE (Noll et al. 2009), and Prospector (Leja et al. 2017) to flexibly describe SFHs. As shown in our validation tests (Section 3), this description minimizes the bias in estimating SFHs at all lookback times, compared to existing parametric and nonparametric methods (Lee et al. 2010; Ciesla et al. 2017; Iyer & Gawiser 2017; Carnall et al. 2019; Leja et al. 2019a).

This subsection describes the methodology for creating SFHs using this formalism for a given number of parameters, and the following subsections handle the construction of the SED fitting machinery needed to infer the optimal number of SFH parameters corresponding to the amount of information available in individual galaxy SEDs.

We define an SFH by the tuple: (M_*, SFR, {t_X}), where M_* is the stellar mass, SFR is the star formation rate at the epoch of observation, and the set {t_X} contains the N "shape" parameters that describe the SFH. The shape parameters, {t_X}, are N lookback times at which the galaxy formed equally spaced quantiles of its total mass (Pacifici et al. 2016; Behroozi et al. 2018). For the first few values of N, we can write

$$\begin{eqnarray*}\begin{array}{rcl}N & = & 1\,\,P=\{{t}_{50}\}\\ N & = & 2\,\,P=\{{t}_{33},{t}_{67}\}\\ N & = & 3\,\,P=\{{t}_{25},{t}_{50},{t}_{75}\}\\ N & = & 4\,\,P=\{{t}_{20},{t}_{40},{t}_{60},{t}_{80}\}.\\ & ... & \end{array}\end{eqnarray*}$$

This can be seen graphically in Figure 1, where we show the {t_x} parameters with N = 2, 4, 9 (vertical dashed black lines) for a mock SFH (blue line), along with our construction of the SFH using these parameters (black solid lines). As expected, the shape of the mock SFH is better captured as N increases, with multiple episodes captured using four parameters. In practice, we find that it is possible to recover multiple episodes of star formation with as few as three {t_x} parameters. Together with the stellar mass and SFR, this tuple describes a set of integral constraints that describe the shape and overall normalization of the SFH. For a galaxy at redshift z, when the universe was t_z Gyr old, the constraints are:

$$\begin{eqnarray}\begin{array}{rcl}{M}_{* }({t}_{z}) & = & {\displaystyle \int }_{t=0}^{{t}_{z}}\mathrm{SFH}(t){f}_{\mathrm{ret}}(t-{t}_{z},Z){dt}\\ & & \times \,\left\{\displaystyle \frac{{{iM}}_{* ,\mathrm{tot}}({t}_{z})}{N+1}={\displaystyle \int }_{t=0}^{{t}_{x,i}}\mathrm{SFH}(t){dt}\right\}\forall i\in N\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{M}_{* }=\displaystyle \frac{{\int }_{t=0}^{{t}_{z}}\mathrm{SFH}(t){f}_{\mathrm{ret}}(t-{t}_{z},Z){dt}}{{\int }_{t=0}^{{t}_{z}}\mathrm{SFH}(t){dt}}{M}_{* ,\mathrm{tot}}.\end{eqnarray} \tag{ 2 }$$

The second line is a set of N equations, one for each parameter in the set {t_x} that requires that the galaxy form "x" fraction of its total mass by time t_x.

Zoom In Zoom Out Reset image size

This description of a galaxy's SFH already offers several advantages over methods found in the literature, a few of which are summarized below.

• 1.  
  Not being restricted to a particular functional form minimizes bias due to SFH parameterization (Iyer et al. 2018).
• 2.  
  Describing an SFH using {t_X} reduces the discrepancy in S/N per parameter, in comparison to methods that determine the SFR in bins of lookback time: for example, t₂₀ might be less well-constrained compared to t₈₀, but the overall signal depends on the shape of the SFH. For example, the parameterization will not try to constrain the SFR in the first year after the Big Bang unless enough stars were formed that early to provide a discernible signal in the SED. This can be compared to methods that adaptively choose time bins, e.g., VESPA (Tojeiro et al. 2007).
• 3.  
  This provides a novel framework for compressing the amount of information present in an SFH to a small set of numbers (given a way to reconstruct an SFH from a tuple), and hence for comparing SFHs across different simulations and observations on the same footing.
• 4.  
  The distribution of different {t_X} among galaxies at a given epoch within a simulation can be extremely useful in defining and checking the physical assumptions of the SFH priors during SED fitting.

Reconstructing an SFH from the tuple (M_*, SFR, {t_X}) can be done in multiple ways, but we seek to minimize the information lost in doing so, while remaining computationally inexpensive. As with all compression methods, we approach the true SFH as the number of parameters N in the set $\{{t}_{X}\}\to \infty$, but we would like to minimize the loss even with a relatively small number of parameters.

Reconstructing an SFH requires quantifying the integral constraints as points on a fractional mass (M_*,tot(t))—cosmic time (t) plane and drawing a piecewise smooth curve passing through these points. Rescaling the mass axis and differentiating this cumulative curve would then yield the SFH as a star formation rate at each lookback time. The simplest approximation would be to connect each point such that the resulting SFH is piecewise linear, but while this provides an acceptable solution, it is not very physical—in the sense that taking the derivative yields an SFR with jump discontinuities. While generalizing to polynomials for the interpolation causes problems with the derivative going negative in parts of the SFH, methods such as tensioned cubic splines and piecewise-cubic Hermite polynomial interpolation (De Boor et al. 1978; Butt & Brodlie 1993) provide more sophisticated solutions to this problem.

In this work, we use GPR (Rasmussen & Williams 2006; Leistedt & Hogg 2017) implemented through the george python package (Ambikasaran et al. 2015; Foreman-Mackey 2015) to create a smooth SFH along with uncertainties following a physically motivated covariance function (kernel) for a given SFH tuple. The Gaussian process framework uses a set of constraints, given by the set of Equation (1), along with a covariance function (or kernel) to estimate the probability of SFH(t) at a given time t. We use a Matern32 kernel (Seeger 2004) in the present application, given by:

$$\begin{eqnarray}&&k({r}^{2})=\left(1+\sqrt{3{r}^{2}}\right)\exp \left(-\sqrt{3{r}^{2}}\right)\end{eqnarray} \tag{ 3 }$$

where the covariance function contains a scale length hyperparameter (r² ∝ (Δt)² that sets how much the SFR(t) can vary from the SFR(t+Δt) separated by a time interval Δt. The hyperparameter in the kernels essentially encodes the amount of stochasticity in the SFHs, and is described further in Section 3.5. In practice, this can be thought of as setting the tension in a string that passes through all the constraints in cumulative mass–cosmic time space, and therefore affects the shape and amount of ringing that can happen between two quantiles. When using a spline to reconstruct the SFH, this ringing can sometimes cause the SFR to be negative. Our choice of physically motivated kernel minimizes this behavior and limits it to extreme SFH tuples (e.g., t₂₅ = 0.1 Myr, t₅₀ = 10 Gyr, t₇₅ = 10.1 Gyr, i.e., the SFH formed the first and third quarters of its mass in bursts of 0.1 Myr, separated by a long period of relative quiescence), which account for <3% across our basis. In these cases, we set the relevant portion (i.e., SFR(t) < 0) to 0 and check that the overall error in stellar mass due to this is <5%. By doing so, we ensure that SFR ≥ 0 at all times across our entire basis. Setting the SFR to 0 in cases like this is physically motivated—stretches with SFR = 0 might be found in nature as well, as seen in SFH reconstructions of local dwarfs (Weisz et al. 2011) or in high-resolution simulations (Wright et al. 2019), for galaxies with very bursty SFHs. Examples of approximating SFHs using this formalism are shown in Figure 1. The advantage of this method is that both parameters (M_*, SFR, {t_x}) and uncertainties on these parameters can be passed as arguments while constructing the SFH posterior. This is tested using SFHs from the Santa Cruz semi-analytic models (SAM) (Somerville et al. 2008; Porter et al. 2014) and the MUFASA Davé et al. (2016) simulation to minimize loss during the reconstruction of SFHs; in both cases, the error in approximating the true SFH decreases as we increase the number of parameters. The Santa Cruz semi-analytic model uses dark matter halo merger trees from the Bolshoi-Planck simulation (Klypin et al. 2011) in combination with analytic recipes for the radiative cooling of gas, collapse of cold gas to form rotationally supported disks, star formation, feedback, and chemical enrichment. The MUFASA hydrodynamical simulation uses the GIZMO (Hopkins 2015) meshless code, including prescriptions for H₂-based star formation, chemical evolution, and kinetic outflows following scalings from the FIRE (Hopkins et al. 2014) simulation, as well as an evolving halo-mass–based quenching mechanism. Both simulations reproduce a wide range of present-day observables, including z ∼ 0 mass functions and scaling relations, with further details in Somerville et al. (2008), Porter et al. (2014), and Davé et al. (2016).

While we are primarily interested in the SFH, determining it from the galaxy's SED requires us to account for several other factors, such as the chemical enrichment, stellar initial mass function (IMF), dust attenuation model, and absorption by the IGM. We then formulate the SFH estimation as an inference problem given by:

$$\begin{eqnarray}&&\begin{array}{l}P(\mathrm{SFH},{A}_{V},z...| {F}_{\nu ,j}^{\mathrm{obs}})\\ =\,\displaystyle \frac{P({F}_{\nu ,j}^{\mathrm{obs}}| \mathrm{SFH},{A}_{V},Z,z...)P(\mathrm{SFH},{A}_{V},Z,z...)}{P({F}_{\nu ,j}^{\mathrm{obs}})}.\end{array}\end{eqnarray} \tag{ 4 }$$

The term $P({F}_{\nu ,j}^{\mathrm{obs}}| \mathrm{SFH},{A}_{V},Z,z)$ is the likelihood, given by ${ \mathcal L }\propto \exp (-{\chi }^{2}/2)$, where

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{j=1}^{{N}_{\mathrm{filters}}}{\left(\displaystyle \frac{{F}_{\nu ,j}^{\mathrm{obs}}-{F}_{\nu ,j}^{\mathrm{model}}(\mathrm{SFH},{A}_{V},Z,z)}{{\sigma }_{j}}\right)}^{2}.\end{eqnarray} \tag{ 5 }$$

The term $P(\mathrm{SFH},{A}_{V},Z,z...)$ denotes the prior distribution for the model. If we assume uncorrelated priors for all the parameters, this can be written as $P(\mathrm{SFH})P({A}_{V})P(Z)P(z)...$. Here, ${F}_{\nu ,j}^{\mathrm{obs}}$ is the observed photometry being fit, in the jth photometric filter. SFH denotes the SFH tuple (M_*, SFR, {t_X}), A_V is the dust model, Z the stellar metallicity, and z is the redshift. In addition to this, we need to consider the stellar population synthesis models, stellar IMF, absorption by the intergalactic medium, and a self-consistent implementation of nebular emission lines using CLOUDY through the Flexible Stellar Population Synthesis (FSPS) code. We adopt the Calzetti attenuation law for the dust attenuation (Calzetti 2001), with one free parameter (because the Calzetti law couples the birth-cloud attenuation to that from older stars), one for stellar metallicity, and a Chabrier IMF (Chabrier 2003) with no free parameters. We define our SFH parameterization in Section 2.1 with N+2 parameters, where N is the number of SFH percentiles, given by the set {t_x}. We use FSPS (Conroy et al. 2009; Conroy & Gunn 2010) to generate spectra corresponding to the Basel stellar tracks and the Padova isochrones. With N SFH quantiles, the model then has N+5 (M_*, SFR, A_V, Z, z) free parameters that need to be determined from the data. We construct the method in such a way that N itself is a variable that is tuned to extract the maximum amount of information present in a galaxy's SED. It is important to keep in mind that these modeling choices can impose an implicit prior on our SED fits. The effects of testing our model assumptions and priors are further explored in Section 3, where we find that our models and priors are suitably robust, considering the S/N and wavelength coverage of our data set.

We implement the posterior computation numerically, using a brute-force Bayesian approach similar to Pacifici et al. (2012, 2016) and Da Cunha et al. (2008), with a large pregrid of model SEDs constructed through random draws from the prior distributions corresponding to each free parameter in Equation (5). To ensure that the pregrid samples the priors finely enough and is effectively dense in SED space, we perform fits to a sample of 1000 galaxies while varying the size of our pregrid. Using this, we estimate the optimal size of the pregrid as the point where the improvement in median χ² for the sample as a function of pregrid size is negligible, leading to a pregrid with ∼900,000 SEDs.

To construct the pregrid, we draw random values from our prior distributions for stellar metallicity, dust attenuation, and SFH parameters (M_*, SFR, {t_x}). For metallicity, we adopt a flat prior on $\mathrm{log}Z/{Z}_{\odot }$ (Pacifici et al. 2016; Carnall et al. 2018; Leja et al. 2019a), an exponential prior on dust attenuation (Iyer & Gawiser 2017), and a Dirichlet prior for the lookback times {t_x} that specify the shape of the SFH. The Dirichlet prior is a generalization of the Beta distribution to N variables, such that a random draw yields N random numbers x_i that satisfy ${\sum }_{i}{x}_{i}=1$. More details about this prior can be found in Leja et al. (2017, 2019a). In practice, for a galaxy SFH at redshift z, we generate the set {t_x} by performing a random draw multiplied by the age of the universe at that redshift, giving the set of lookback times at which a galaxy formed various quantiles of its stellar mass. The Dirichlet prior has a single tunable parameter α that specifies how correlated the values are. In our case, values of this parameter α < 1 result in values that can be arbitrarily close, leading to extremely spiky SFHs because galaxies have to assemble a significant fraction of their mass in a very short period of time, while α > 1 leads to smoother SFHs with more evenly spaced values that nevertheless have considerable diversity. In practice, we use a value of α = 5, which leads to a distribution of parameters that is similar to what we find in SAM and MUFASA. This can be seen in Figure 2. For each galaxy in our pregrid, we turn the corresponding (M_*, SFR, {t_x}) tuple into an SFH via our Gaussian process routine, assuming fiducial 1% uncertainties. Using this SFH, in combination with the values for metallicity, dust, and redshift, we generate a spectrum through the FSPS stellar population synthesis routine. Multiplying this spectrum with the appropriate filter transmission curves gives the SED corresponding to each tuple in the pregrid.

SFH uncertainties are computed after the fit is performed, using the posterior for each parameter describing the tuple (M_*, SFR, {t_x}). For each galaxy, 100 self-consistent random draws are performed from the posterior with the covariances between parameters taken into account, and corresponding SFH realizations are constructed using the Gaussian process routine. For a given realization, if the set of parameters already exists in the pregrid, the corresponding precomputed SFH is used to decrease the computational cost. Figure 3 shows 20 draws from the SFH posterior for that galaxy using this approach, depicted as thin black lines in the SFH inset panel. The 68% confidence intervals are then constructed by taking the 16th to 84th percentiles of the SFR distribution at each point in lookback time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Gaussian process–based SFH (GP-SFH) formalism allows us to gain independence from having to make a choice of functional form for the shape of the SFH, without having to bin SFHs in lookback time. This results in smooth, effectively nonparametric SFHs that minimize the bias and scatter in SFH reconstruction, as seen in Section 3. However, because we would like the method to be truly nonparametric, we require that the number of {t_x} parameters be optimized for the amount of information present in a given noisy SED. To implement this in practice, we generate pregrids for SFHs with the set {t_x} ranging from three to nine parameters, because it is difficult to specify the shape of complex SFHs with less than three parameters and it is impractical to recover more than nine from broadband SEDs. We then fit each observed SED with all seven pregrids, and obtain the most appropriate number of parameters using an appropriate model selection criterion. Ideally, the Bayesian evidence would be used for this model selection step (Liddle 2007). However, in practice, the numerical computation of the evidence is expensive due to the need for a nested sampler, and cannot be completed for the number of galaxies typically found in large photometric surveys. Having tested the properties of the likelihood surfaces using this method, however, we find that while they may be multimodal, they generally do not contain pathological features that necessitate this numerically expensive procedure. In light of this, we perform our model selection using an approximation of the evidence, given by the Bayesian Information Criterion (BIC) (Schwarz et al. 1978; Liddle 2007), defined as

$$\begin{eqnarray}&&\mathrm{BIC}=k\mathrm{ln}(n)-2\mathrm{ln}{{ \mathcal L }}_{\max }\end{eqnarray} \tag{ 6 }$$

where n is the number of photometric data points, k is the number of parameters in the SFH tuple, and ${ \mathcal L }$ is the maximized value of the likelihood function given by Equation (5). The latter term in this equation is a measure of how well the model describes the data, and the former term is a penalty for an increased number of parameters. The results from the two model selection criteria are equivalent when all the parameters are independent (Kass & Wasserman 1995; Szydłowski et al. 2015). However, if parameters are strongly covariant, the volume will reflect this (because the intrinsic dimensionality is lower than the number of dimensions describing the likelihood) while the BIC does not. This can lead to the optimal number of parameters estimated with the BIC ≤ the optimal number of parameters estimated using the evidence. Our validation tests show that using the BIC does not cause any significant issues in fitting an ensemble of SEDs to recover their SFHs. By finding the minima in the BIC, we find the number of SFH parameters that can be robustly extracted from the SED being fit. This leads to a truly nonparametric description of the SFH, based on the amount of information about the different stellar populations encoded in the galaxy's SED. An example of this for a single galaxy is shown in the bottom right plot of Figure 2. Iyer & Gawiser (2017) used a similar nonparametric estimate of the number of SFH components using a F-test, but it was computationally expensive to implement as the number of SFH components grew. Tojeiro et al. (2007) and Dye (2008) also used nonparametric methods to estimate the optimal number of time bins during SFH reconstruction, although this is prone to bin edge effects, as described in Leja et al. (2019a). Our Gaussian process SFH reconstruction method can be summarized as follows:

• 1.  
  We introduce a formalism based on Bayesian statistics to describe the SFH, which we will express not as a functional form but rather as N-tuples, plus M_* and SFR at the time of observation. To generate a curve SFR(t) from the tuples, the GPR needs a covariance function or kernel, which determines how SFR distributions at t and t + Δt are related, i.e., how smooth the SFHs are.
• 2.  
  The Bayesian likelihood is computed for each of a sufficiently large number of model SFHs realized using GPR, having known values of the n-tuples, M_*, and SFR-today (the pregrid) drawn from our chosen prior. By estimating the likelihood of each basis N-tuple, we can compute posteriors for the N-tuples (and M_* and SFR-today) along with their credible intervals and covariances.
• 3.  
  The GPR is now used to interpolate between these (noisy) estimates of the n-tuples to obtain draws from the posterior distribution. Each of these draws is a model SFH. The distribution of these draws at each point in time gives the uncertainty in the SFH at that time.
• 4.  
  This procedure is repeated as we vary the number of parameters in the N-tuples, using a Bayesian model selection criterion to evaluate the optimal number of parameters.

Unlike simple SFH parameterizations commonly used in SED fitting, the Gaussian process–based SFHs can have multiple maxima, even with four or fewer t_X parameters, as seen in Figure 1. Hence, it is possible to analyze them using a peak-finding algorithm to quantify the number of major episodes of star formation in a galaxy's past. For this particular analysis, we use the best-fit SFH for each galaxy, because the median SFH for each galaxy is biased toward smooth SFHs. This is an effect due to our choice of kernel, which prefers smooth solutions when the {t_x} parameters are uncertain, as seen in the left column of Figure 1. With tighter constraints on {t_x} from higher-S/N data, this problem is alleviated. As a result, while the number of major episodes (N_ep) estimated using this method for individual galaxies is susceptible to noise, the overall distribution is seen to be recovered without any significant bias, as shown in Section 3.

For each SFH, we quantify the number of episodes as follows: we first find the number of peaks in an SFH as the set of points that satisfy dSFH/dt = 0 and ${d}^{2}\mathrm{SFH}/{{dt}}^{2}\lt 0$. To separate multiple peaks within an overall episode of star formation from different episodes, we impose a peak-prominence criterion by requiring that

$$\begin{eqnarray}&&(\mathrm{log}\,{\mathrm{SFR}}_{\mathrm{peak}}-\mathrm{log}\,{\mathrm{SFR}}_{\min ,\mathrm{local}})\gt 1.5+\displaystyle \frac{1.5}{4}\mathrm{log}\displaystyle \frac{{M}_{* }}{{10}^{8}{M}_{\odot }}\end{eqnarray} \tag{ 7 }$$

where SFR_min,local is the minima between two peaks in the SFH. This condition is shown (in Section 3) to minimize the type-1 (overestimating N_ep) and type-2 (underestimating N_ep) errors in our validation sample. It arises because the sensitivity to star formation drops approximately logarithmically with time (Ocvirk et al. 2006). Because more massive galaxies tend to have older stellar populations, we found that a mass-independent peak-prominence criterion caused a mass-dependent bias in our estimates of the number of episodes. We correct for this effect by requiring a more (less) stringent dip in the SFH for a massive (low-mass) galaxy, and while this does not improve the result for every galaxy, it accurately recovers the distribution, which is the quantity that we are interested in.

To quantify possible biases in our SFH reconstruction, we plot the difference between the true and recovered SFH for our sample of validation galaxies in Figure 4. The two panels show the linear and logarithmic differences between the true and the recovered SFHs, which respectively represent the difference and ratio of the true versus predicted SFR(t). This difference is smaller than 0.3 dex for most of the lookback time. While the reconstructed SFHs cannot recover every short episode of star formation, as seen in the errors for individual galaxies (thin black lines), the overall SFH for the ensemble of galaxies is unbiased out to nearly ∼8 Gyr. In the ${\rm{\Delta }}\mathrm{log}\,\mathrm{SFR}$ error plot, the bias blows up as we approach the Big Bang. This is because the SFHs in the mocks can abruptly fall to 0 close to the Big Bang while the Gaussian process SFHs smoothly decline to 0 at t = 0, leading to a one-sided error. However, the arithmetic ΔSFR plot shows that this is, in fact, a very small difference, exaggerated by the fact that log ${\mathrm{SFR}}_{\mathrm{true}}\to -\infty$ as SFR_true → 0 close to the Big Bang.

Zoom In Zoom Out Reset image size

In addition to the SFHs, Figure 5 shows the results of estimating traditional SED fit parameters—stellar mass, star formation rate averaged over the last 100 Myr, dust attenuation, and stellar metallicity for the validation catalog. The stellar masses are constrained better than the traditional scatter around the mean of ∼0.14 dex (Mobasher et al. 2015). Star formation rates have a larger scatter (0.29 dex) due to degeneracies with dust attenuation (0.13 dex) and metallicity (0.32 dex). At low stellar masses, we are prone to overestimating the overall mass due to our choices of SFH prior dominating the fit for low-S/N SEDs. Restricting the fits to galaxies with H < 25 is found to largely exclude these low-S/N objects, leading to more robust estimates of these parameters.

Zoom In Zoom Out Reset image size

As described in Section 2, SFH uncertainties are computed at each point in lookback time for various percentiles, ranging from the 10th (55–65th percentile bounds) to the 95th (2.5–97.5th) percentile credible intervals of SFHs constructed from 100 draws from the SFH posterior (M_*, SFR, {t_x}) using the Gaussian process routine. To check whether these represent true 68% uncertainties, we consider the uncertainties estimated for each of the 10,000 galaxies in the mock catalog that we fit and compute the fraction of the true SFH that lies within the uncertainties for each galaxy. This is similar to the concept of a "coverage probability" (Levasseur et al. 2017) used to evaluate the accuracy and precision of uncertainty estimates. The distribution of this fraction is given in the left panel of Figure 6, which shows the fraction of the truth that lies, on average, within a given credible interval. For an ideal method, we expect this to lie along the 1:1 line; our method is extremely close to this across the entire evaluated range. This indicates that our uncertainties are robust, expanding on the tests in Iyer & Gawiser (2017), which solely checked the 68% uncertainties. Figure 6 shows two examples: one where the uncertainties are representative of the truth, and another where both the reconstruction and the uncertainties miss a relatively short (<0.5 Gyr) episode of star formation that follows a relatively quiescent period, due to the smoother reconstructed SFH producing a comparable SED given the photometric noise.

Zoom In Zoom Out Reset image size

For each galaxy in the mock sample, we compute the number of major episodes of star formation (N_ep) using the peak-finding algorithm described in Section 2. For the true SFHs, we compute the number of episodes requiring a dip of $\mathrm{log}\,{\mathrm{SFR}}_{\mathrm{peak}}-\mathrm{log}\,{\mathrm{SFR}}_{\min ,\mathrm{local}}\geqslant 0.3$ dex, to model the scenario of a galaxy on the star-forming sequence dropping below the sequence before rejoining it. The distribution of galaxies with multiple strong episodes of star formation for the mock sample is then given by the thick dashed black lines in Figure 7. We then estimate the number of episodes for the reconstructed SFHs with the same criterion, and find a flatter trend than observed in the mocks. This is a combination of two effects. (1) The S/N contributed to the overall SED by a stellar population of a certain age decreases roughly logarithmically with lookback time (Ocvirk et al. 2006), leading to poorer constraints on intermediate and older stellar populations. (2) There is a mass-dependent correlation toward older stellar populations as galaxies grow more massive. As a result of this, galaxies that formed most of their mass prior to an extended period of quiescence can be erroneously classified as having multiple episodes, because when the SFR at the time of observation is low, it is easier to have a fluctuation of ≥0.3 dex with a small variation in SFR. We find that introducing a stellar mass-dependent threshold to determine the number of peaks given by Equation (7) accounts for these mass-dependent effects on the distribution, allowing us to better reproduce the true trends in N_ep with stellar mass and redshift, as seen in Figure 7. The thin lines also show the number of overestimates and underestimates, such that N_ep,true = N_ep,rec – underestimates + overestimates. Because we are using the best-fit SFHs for this computation, we find that our results are quite sensitive to noise for individual galaxies. However, the distributions are estimated robustly, across a range of stellar masses and redshifts. Additionally, the differing behavior of the observational sample from our mock catalog is indicative of the fact that our mass-dependent threshold does not impose an artificial trend on the distribution. Reconciling the differences between the two distributions would be an interesting topic for further study.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To generate a curve in SFR versus lookback time space, we use GPR to perform interpolation in cumulative mass ($\int \mathrm{SFR}(t){dt}$) as a function of time, conditioned on constraints given by the N-tuples (SFR, M_*, {t_X}), as described in Section 2. The choice of kernel can make a significant difference to this interpolation, because it determines the extent to which the cumulative mass (and hence the SFR) at two times (t and t') are correlated (Rasmussen & Williams 2006; Foreman-Mackey 2015). In the current implementation, we choose to focus on kernels that are one-dimensional and stationary, which means that the kernel only depends on the separation between two times and not the times themselves. This choice of kernel is physically motivated in the sense that SFRs over short timescales are thought to be strongly correlated, with this correlation dropping off as the time separation increases, similar to the models explored by Kelson (2014) and Caplar & Tacchella (2019). To determine the most appropriate kernel from those widely used in the literature, we consider four stationary kernels implemented in the george Python package: the squared-exponential or Radial Basis kernel, the Matern32 and Matern52 kernels from the Matern class of kernels (Seeger 2004), and the Rational Quadratic kernel, which can be thought of as adding together many Radial Basis kernels with differing length scales. We test all of these different kernel functions using the SFHs from our mock catalog. An example of the reconstruction for a single SFH using different kernel choices is shown in Figure 8. In general, we find that the choice of kernel only makes a percent-level difference in computing the median SFH, although the uncertainties can be significantly different. We choose the Matern32 kernel because it gives the best reconstructions and most accurate uncertainties, as explored in Section 3.3. The Matern32 kernel as implemented in the george python package is given by the equation

$$\begin{eqnarray}&&k({r}^{2})=\left(1+\sqrt{3{r}^{2}}\right)\exp \left(-\sqrt{3{r}^{2}}\right),\end{eqnarray} \tag{ 8 }$$

where r² is related to the separation between two times. Because our kernel is one-dimensional and stationary (independent of choice of origin), this simplifies to r² ∝ (t − t')², where r² is related to the separation by a scale length. This scale length is determined directly within the george package by minimizing the log-likelihood of the data under the Gaussian process model. While the scale length varies with the number of SFH percentiles (N), because there is more information to condition the likelihood with for greater N, it is fixed across the pregrid for a given choice of N. In Figure 9, we test the effect of increasing and reducing the scale length used in the Gaussian process kernel. This causes the covariance matrix to become more (less) coupled to SFRs separated by a given time interval. While this does not have a major effect on the reconstructed SFH, it can cause a decrease (increase) in the uncertainties estimated using random draws from the SFH posterior. Taken in conjunction with our validation of the SFH credible intervals in Section 3.3, we find that our choice of kernel and hyperparameters robustly reconstructs SFHs with accurate uncertainties. In addition to the basic kernels described in this section, it is possible to design more physically motivated kernels by computing the covariance function of SFHs from cosmological simulations. However, because this is subject to a wide variety of systematics depending on resolution, simulation size, and differing physical prescriptions, such an analysis requires a detailed comparison of different simulations. Although this out of the scope of the current paper, it is being actively investigated for a follow-up paper.

Zoom In Zoom Out Reset image size

We test the effectiveness of the Gaussian process–based SFH reconstruction method in compressing data about the SFH and encoding the information in a small number of parameters. Because the SFH percentiles store information about the shape of the SFHs, the accuracy of describing an arbitrary SFH should decrease as we increase the number of parameters (N). We see this to be the case in Figure 10, where the error in describing an ensemble of randomly drawn SFHs from our mock catalog decreases with increasing N. This analysis is done purely in SFH space as a test of the loss of information in going from the full SFH to the GP-SFH tuple. In addition to this, we compare the method with existing nonparametric techniques of approximating an arbitrary function using: 1) polynomial coefficients, and 2) binning the SFH linearly in time such that the SFR in any given bin is a constant. The latter approach and its variants—such as binning in logarithmic lookback time or using adaptive bins—are seen across the literature in works that explore nonparametric SFH reconstruction (Heavens et al. 2000; Tojeiro et al. 2007; Da Cunha et al. 2008; Dye 2008; Leja et al. 2017). While these are generally effective, they have the disadvantages of introducing unphysical discontinuities at bin edges. Repeating the test with these two methods, we find that their error approaches the GP-SFH error for large N, with the GP-SFH performing better for any given N. Tuning the number of parameters during SED fitting therefore optimally estimates the amount of SFH information from each SED being fit.

Zoom In Zoom Out Reset image size

We now apply the Dense Basis SED fitting method as described in Section 2 to the sample of CANDELS galaxies at 0.5 < z < 3.0, to reconstruct the SFH of each galaxy with uncertainties. Figure 12 shows the distribution of the number of SFH parameters that individual galaxies are best fit with in each of the five fields, with a slight trend toward increasing amounts of SFH information recovered with more photometric bands.

Zoom In Zoom Out Reset image size

In Figure 13, we show the distributions of reconstructed SFHs in the five CANDELS fields and how they evolve with mass and redshift. The median and interquartile range are computed at each point in time as in Pacifici et al. (2016), who performed a similar calculation for quiescent galaxies using a basis of SFHs derived from a semi-analytic model (De Lucia & Blaizot 2007). In the current analysis, we have used a redshift bin width of δz = ±0.1, and a mass bin width of δM_* = ±0.25 dex, such that the M_* ∼ 10¹⁰M_⊙ bin includes all galaxies with 10^9.75 < M_* < 10^10.25M_⊙. Dashed black lines in each panel show the mean SFR (≡M_*/t_univ), assuming constant SFR for that redshift and mass bin, and the f_samp at the top left corner of each panel shows the fraction of all galaxies in our sample at that redshift that fall in a particular mass bin. Because there can be a few galaxies that fall outside the plotted mass range, this may not sum to unity.

Zoom In Zoom Out Reset image size

We see that the median SFH across the five CANDELS fields are remarkably similar, as expected. There are discrepancies in a few of the bins, most notably at $2\lt z\lt 3$ and M_* ∼ 10¹⁰M_⊙ for UDS, which could be the result of correlated photometric noise or smearing effects in the photo-z that was used for the calculation. In the highest redshift bin, the portion of the SED with wavelengths greater than rest-frame 1.6 μm is not well-sampled, leading to poorer constraints on the SFHs in the the last row.

In general, SFHs tend to rise at high redshifts and low stellar masses, similar to those from cosmological simulations, with a turnover and subsequent decline as we go to higher masses and lower redshifts. In agreement with Pacifici et al. (2012) and Behroozi et al. (2018), massive galaxies tend to peak earlier in their SFHs, and galaxies on average tend to move toward quiescence at lower masses as the universe grows older, with the threshold changing from nearly 10^11.5M_⊙ at z ∼ 3 to 10¹⁰M_⊙ at z ∼ 0.5, without the need for any implicit assumptions about the SFHs, such as the well-motivated assumption of Behroozi et al. (2018) that earlier-forming halos get lower SFRs. The additional advantage of the Dense Basis method is that, in addition to the average SFHs, the individual SFHs of galaxies reconstructed from the observations allow us to explore the additional factors that drive the diversity of SFHs at a given mass and epoch. This allows us to extend our analysis beyond the physics encoded in the stellar mass–halo mass relation, which gives a constraint upon the first-order behavior of SFHs.

Most studies of galaxy SFHs focus on the overall rise and fall of an ensemble of SFHs (Leitner 2012; Pacifici et al. 2016; Ciesla et al. 2017; Leja et al. 2017), which has led to a well-constrained understanding of the overall behavior seen in Figure 13. However, with smooth, nonparametric SFHs, it is now possible to ask questions about the second-order statistics of an ensemble of SFHs, analyzing the departures from this overall behavior in the form of periods of relative quiescence between episodes of star formation.

In Figure 14, we show the fraction of galaxies with different numbers of major episodes of star formation in a galaxy's past at redshifts 0.5 < z < 3.0. Because the number of episodes is a discrete quantity, Poisson noise dominates the formal uncertainties in an individual galaxy's number of episodes while calculating functions of the sample, as discussed in Section 3. The different fields (colored lines), are in good agreement with each other and the median of the full sample (solid black line).

Zoom In Zoom Out Reset image size

We find that, at low redshifts, the fraction of galaxies with multiple major episodes of star formation decreases as we go up in stellar mass above 10^10.5M_⊙, in agreement with Iyer & Gawiser (2017). In addition to this, we find a slight decrease in the overall fraction of galaxies with multiple episodes with increasing redshift at any given mass, with the notable exception of M_* ≈ 10^10.5M_⊙, which does not show a noticeable evolution with redshift. Although we have accounted for S/N variations, the decrease at lower masses could be at least partially due to having insufficient S/N to resolve multiple episodes of star formation as we go to lower masses and higher redshifts. A few explanations are possible for the behavior at high masses: AGN feedback quenching galaxies (Weinberger et al. 2018) could lead to SFHs that form most of their stars at by z ∼ 3, which could look like a single early episode of star formation without sufficient S/N in the SEDs to resolve the older populations at z ∼ 1. This is made more probable by the fact that, while most galaxies with multiple episodes are found to lie on the SFR–M_* correlation, the greatest number of galaxies with low SFRs at the time of observation and multiple episodes occurs at masses close to 10^10.5M_⊙. Another reason could be the central limit theorem (Kelson 2014): massive galaxies that grow primarily through mergers (Brinchmann et al. 2004; Bundy et al. 2005; Pérez-González et al. 2008) at early times could be composed of multiple progenitors. As the number of progenitors grows with mass, by the central limit theorem, their SFHs should look smoother than those for less-massive galaxies. Low-mass galaxies, by comparison, should have more stochastic SFHs because they are growing most of their mass in situ, which would be in close agreement with the findings of Guo et al. (2016), Matthee & Schaye (2019), Emami et al. (2018), Shivaei et al. (2015), and Broussard et al. (2019).

If a galaxy is found to have multiple strong episodes of star formation in its lifetime, an interesting question would be whether the galaxy was actively forming stars and the star formation was temporarily suppressed by a quenching attempt (short time interval between peaks), as opposed to a galaxy that was on its way to quiescence but restarted star formation due to an inflow of pristine gas or merger. This is especially interesting within the context of rejuvenation of galaxy SFRs (Fang et al. 2012), because it sets timescales for how long a galaxy spends off the star-forming sequence when it makes such an excursion. To quantify this, we measure the time interval between multiple peaks for the subsample of galaxies with N_ep > 1. We plot this as a function of redshift and mass in Figure 15, finding that although the separation does not vary strongly with mass, it does show a strong trend with redshift. However, upon normalizing by the age of the universe at different redshifts, this trend is significantly decreased, leaving us with a roughly constant timescale across which galaxy rejuvenation occurs, given by ${t}_{{\rm{\Delta }}\mathrm{peak}}\sim {0.42}_{-0.10}^{+0.15}{t}_{\mathrm{univ}}\,\mathrm{Gyr}$, where t_univ is the age of the universe at the redshift of observation. This is similar to the result by Abramson et al. (2016), which found the transit time through the green valley (i.e., half the time between two peaks for a rejuvenating SFH) to be ∼0.2t_univ Gyr, roughly independent of redshift (see also Pandya et al. (2017), which studies the transition timescales in massive CANDELS galaxies using a statistical analysis of their number densities). This is also related to the results of Pacifici et al. (2016), who that found that the width of the SFH for quiescent galaxies is roughly constant across stellar mass and redshift when the age of the universe is factored out, as well as those of Muzzin et al. (2014), who found that the post-starburst spectra of galaxies at z ∼ 1 are fit well with a quenching timescale of ${0.4}_{-0.3}^{+0.4}\,\mathrm{Gyr}$. Fang et al. (2012) identify a subsample of galaxies at z ∼ 0.1 that could linger in the green valley for ${ \mathcal O }$(Gyr). The astute reader may wonder how significant it is to find that two episodes are typically separated by roughly half the age of the universe at the time of observation, as this also corresponds to the median period of a generic sine wave possessing two peaks without that interval. Given the present data quality, it is difficult to test this further by investigating the separation between the earliest two star formation episodes in galaxies whose SFHs show 3–5 major episodes of star formation, but such a test should be done with higher-S/N spectrophotometry. At present, we can compare the fit for separation between episodes against the predictions of galaxy formation models, finding similar trends albeit with a slightly smaller value of ≈0.3 ± 0.15t_univ Gyr. The consistency of our result with a variety of similar results across a range of redshifts summarized above is an additional reassuring check. In Behroozi et al. (2018), galaxy rejuvenation is a generic feature of a population, with the timescale depending on the mode by which a halo is accreting mass, whether through mergers, accretion from another halo, or infall. In this scenario, rejuvenation occurs more often when the quenched population evolves more slowly than the halo dynamical time, during which it can switch between modes of accretion, increasing or decreasing the SFR as a result. This results in estimates of the fraction of galaxies with past rejuvenation as a function of mass, decreasing as redshift increases and showing a trend in stellar mass that is consistent with our results. Our result seems to indicate that the rejuvenation timescales remain relatively constant over cosmic time. It is important to note that the scatter in this quantity is quite large, and the evolution in the quenched fraction happens most rapidly at redshifts 0 < z < 0.5 (Muzzin et al. 2013; Behroozi et al. 2018; Donnari et al. 2019; Hahn et al. 2019), so extrapolation to that regime needs to be tested with further data.

Zoom In Zoom Out Reset image size

In keeping with the finding of Pacifici et al. (2016) that the widths of the SFHs of passive galaxies are roughly constant upon factoring out the age of the universe, and our similar finding for the time interval between two peaks of star formation, we consider galaxy SFHs binned in t₅₀. We bin galaxies in four bins, 0.1t_univ < t₅₀ < 0.3t_univ, 0.3t_univ < t₅₀ < 0.5t_univ, 0.5t_univ < t₅₀ < 0.7t_univ, and 0.7t_univ < t₅₀ < 0.9t_univ, at different redshifts. We show the results in Figure 16. As in Figure 13, each panel lists the fraction of the sample in a particular bin. However, in this case, the fractions are no longer tracing the mass function of galaxies. Instead, the four bins in time serve as proxies for galaxy SFHs at different stages of their lifetimes. This enables us to identify different populations of galaxies, including starbursting galaxies, late bloomers (Dressler et al. 2016), star-forming galaxies, post-starburst or green-valley galaxies (Fang et al. 2012), and quiescent galaxies, using different redshift-dependent t₅₀ cuts, either independently or in combination with other factors like t₂₅, t₇₅, size, and morphology. To further interpret these SFHs, Figure 17 displays the positions of the galaxies within each panel in Figure 16 on the SFR–M_* plane. We find that the left (right) columns in both figures select star-forming (quiescent) galaxies that lie on (off) the star-forming sequence. At intermediate values of t₅₀, the populations are a combination of star-forming and quiescent galaxies, with the selection gradually shifting from star-forming to quiescent galaxies. The intermediate t₅₀ panels also show an excess of galaxies with multiple strong episodes of star formation (N_ep > 1). While this is an intuitive result, because galaxies that have assembled most of their mass recently or those that have long since shut off their star formation are not very likely to contain multiple episodes of star formation, it has important implications for analyses that assume galaxies evolve along smooth SFH trajectories (Leitner 2012) or are described by simple parametric forms (Dressler & Abramson 2014; Ciesla et al. 2017; Lee et al. 2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The morphologies of galaxies are seen to strongly correlate with their stellar masses and redshifts (Conselice 2014), as well as sSFR (Whitaker et al. 2015). While this is a combined effect of the different processes that regulate star formation within galaxies, including mergers, gas accretion through inflows, and stellar and AGN feedback (Boselli & Gavazzi 2006; Hopkins et al. 2008, 2014; Anglés-Alcázar et al. 2017; Weinberger et al. 2018), it is difficult to observationally disentangle the relative strengths of these effects. However, the different timescales that these processes act on enable us to discriminate between the relative effects of these processes, if we can observationally constrain the timescales on which morphological transformation occurs across a population of galaxies—as was done for groups in, e.g., Kovač et al. (2010). The reconstructed SFHs of galaxies provide a direct probe of these timescales by correlating the morphologies of galaxies with their SFHs as compared to indirect measurements of timescales through the frequencies of different morphologies, which are subject to a variety of systematics and selection effects.

We use the CANDELS wide morphology catalogs by Kartaltepe et al. (2015) to study the SFHs of galaxies with different morphological features at 0.5 < z < 1.0. We limit our redshift range to avoid the effects of small number statistics of classifications as we go to higher redshifts. The morphology catalogs contain visual classifications for over 50,000 objects spanning 0 < z < 4 with f160w < 24.5, which gives a large overlap with our sample. The catalogs contain flags for main morphology class (disk, spheroid, peculiar/irregular, point source/compact, and unclassifiable), a class for mergers and other interactions, and structure flags for bars, tidal features, spiral arms, and more. For each class and flag, the catalog reports the fraction of classifiers who were confident about the existence of that feature.

We use this to analyze the SFHs of six classes of galaxies, described as follows.

• 1.  
  Disk: (f_disk > 0.9) AND (f_sph, f_irreg < 0.1). This includes the set of all galaxies classified as disky.
• 2.  
  Disk-dominated galaxies: (f_{disk dom} > 0.9). Disks with a central bulge where the disk dominates the structure.
• 3.  
  Bulge-dominated galaxies: (f_{bulge dom} > 0.9). Disks with a central bulge where the bulge dominates the structure.
• 4.  
  Spheroid: (f_sphk > 0.9) AND (f_disk, f_irreg < 0.1). This includes the set of all galaxies classified as spheroidal.
• 5.  
  Galaxies with spiral arms: (f_arms > 0.9).
• 6.  
  Mergers and interactions: galaxies that either appear to have undergone a merger, as evidenced by tidal features, structures such as loops, or highly irregular outer isophotes (f_merger > 0.9), or are interacting with a companion galaxy within the segmentation map from SExtractor (f_int1 > 0.9).

The results of this analysis are shown in Figure 18. The Figure contains six sets of panels, one for each subsample of galaxies. The top panel shows where the galaxies lie on the SFR–M_* correlation, and the bottom panel shows the median SFH for the subsample of these galaxies at M_* ∈ [10¹⁰, 10^10.5] M_⊙. This is useful to test feedback-driven models of quenching that posit a correlation between bulge–total ratios and SFH shape (Zolotov et al. 2015; Belfiore et al. 2016; Tacchella et al. 2016; Abramson et al. 2018). While the SFHs of our pure disk population at M_* ∼ 10^10.25 M_⊙ seem to actively form stars throughout their lifetime, with the maximal peak in their SFHs' maximum SFR close to the time of observation, the galaxies containing a disk and bulge component seem to show a downward trend in their median SFHs, with disk-dominated galaxies peaking earlier on average than pure disks without a bulge, followed by a decline in SFR. This trend continues to bulge-dominated galaxies and spheroids, showing an evolution in timescales that can be tested with simulations implementing different models for quiescence. For each population, we consider the time since the SFH of each galaxy peaked (t_fall) and use this distribution to quantify the timescale on which galaxy SFHs began their decline. For each morphologically distinct population, we find the median and 68% of the t_fall distribution for galaxies at 0.5 < z < 1.0 and 10¹⁰ < M_* < 10^10.5M_⊙, given by the following.

• 1.  
  Mergers: ${t}_{\mathrm{fall}}:{0.00}_{+0.39}^{-0.00}\,\mathrm{Gyr}$.
• 2.  
  Galaxies with spiral arms: ${t}_{\mathrm{fall}}:{0.60}_{+1.54}^{-0.54}\,\mathrm{Gyr}$.
• 3.  
  Disks: ${t}_{\mathrm{fall}}={0.81}_{+2.56}^{-0.80}\,\mathrm{Gyr}$.
• 4.  
  Disk-dominated galaxies: ${t}_{\mathrm{fall}}:{0.70}_{+2.73}^{-0.38}\,\mathrm{Gyr}$.
• 5.  
  Bulge-dominated galaxies:${t}_{\mathrm{fall}}:{2.15}_{+3.07}^{-1.55}\,\mathrm{Gyr}$.
• 6.  
  Spheroids: ${t}_{\mathrm{fall}}:{2.50}_{+2.25}^{-1.60}\,\mathrm{Gyr}$.

Zoom In Zoom Out Reset image size

In the absence of a bulge component, we also see that galaxies with spiral arms inhabit the high-stellar-mass, high-SFR portion of the SFR–M_* plane, continuing to actively form stars until they lose rotational support or start forming bulges. We also see that mergers and interactions show a noticeable increase in recent SFR, over timescales within the last ∼0.5 Gyr of their SFHs. The timescales for these morphological transformations can be further constrained by determining the resolved SFHs of individual galaxies using IFU surveys like SDSS-IV MaNGA and CALIFA (Delgado et al. 2014; Belfiore et al. 2016).

Because nonparametric methods make no explicit assumption about the form of SFHs, they are only as good as the data being used for SED fitting. In this regard, there are three main avenues for data-driven improvement: better wavelength resolution, better wavelength coverage, and better S/N. Spectroscopy contains more information about stellar populations of different ages and metallicity, as compared to broadband photometry, but often suffers from wavelength-dependent flux calibration issues that need to be accounted for prior to fitting. Panchromatic SEDs allow us to test models of dust attenuation and re-emission to better constrain dust effects while estimating the SFHs of galaxies.

The SFH reconstructions we obtain are also subject to several modeling uncertainties: Stellar Population Synthesis models can introduce systematics into the mapping between physical parameters and observed photometry (Conroy & Gunn 2010; Han & Han 2019). Differences in dust models can introduce systematics into the measurement of recent SFR, which would then propagate into differences in the SFH. C. Pacifici et al. (2019, in preparation) compares the results from 14 different SED fitting codes applied to the same sample of CANDELS/GOODS-S galaxies at z ∼ 1. This allows us to examine the effects of intercode variability and model assumptions for dust, IMF, SFH, and SPS models, including the effects of binary populations. Additionally, Han & Han (2019) tested multiple SPS models, SFH assumptions and dust models using a comprehensive Bayesian formalism that allowed them to estimate the Bayesian evidence in comparing different models.

Finally, the choices of prior assumed during SED fitting are extremely influential in the estimates of physical parameters and their covariances. While we have tried to be agnostic about the priors in this work, it is important to note that an informative prior could be especially useful while fitting noisy, low-S/N data with limited wavelength coverage. Predictive checks could be put in place to ensure that the priors do not introduce significant biases into the estimates, or cause artificially tight correlations due to regression to the mean. These informative priors could be developed by studying the distributions of physical quantities at a particular epoch from a small subset of high-S/N observations, scaling relations, and mass functions, as well as semi-analytic or empirical models that encode the physics that lead to these observables, explicitly quantifying the covariance between star formation, chemical attenuation, and dust enrichment and destruction histories.

The SFHs of galaxies allow us to probe the behavior of mass functions and scaling relations out to higher redshifts than is currently possible. At low to intermediate redshifts, this can be used as a consistency check, to ensure that the reconstructed SFHs are not biased due to noise or prior assumptions. At high redshifts, this can be a powerful tool to increase observational statistics and push measurements out to redshifts higher than those directly accessible through observations. Because we can only measure the SFH summed over all progenitors, it is important to consider the effects of mergers while propagating galaxies backward in time for periods longer than the typical merger timescale (Mantha et al. 2017; Duncan et al. 2019) at a given redshift and stellar mass.

Iyer et al. (2018) propagated galaxies backward in time along their SFHs in the form of trajectories in SFR–M_* space, to probe the high-redshift low-stellar-mass regime of the SFR–M_* correlation. They found that the projected correlation at intermediate redshifts matches the observed distribution well, and extended it by nearly two orders of magnitude, out to z ∼ 6, where observations are extremely faint. C. Pacifici et al. (2019, in preparation) implements a validation test by reconstructing the stellar mass function using galaxies at lower redshifts and comparing them to the stellar mass function obtained through direct fits. Leja et al. (2019b) use SFHs to probe the cosmic SFRD using galaxies at low redshifts, finding that the apparent mismatch between the mass functions and star formation rate functions is alleviated by using nonparametric SFHs.

The smooth, nonparametric SFHs obtained with the improved Dense Basis method offer a window into the pasts of different galaxy populations. While interesting in itself, this insight has the potential to be combined with a variety of ancillary data to probe a wide range of previously inaccessible quantities, some of which we briefly describe below.

• 1.  
  Higher-S/N observations or spectrophotometric data would help obtain better SFH constraints, allowing us to better constrain the number of major episodes of star formation and the timescales on which rejuvenation, starbursts, and quiescence occur at different epochs. While current observations studying the separation between multiple major episodes or transition timescales (this work; Abramson et al. 2016; Pandya et al. 2017) show a flat trend with mass and a linear one with the age of the universe, it is an interesting problem to understand the physical mechanisms responsible for this trend and the dispersion of ∼0.25t_univ Gyr using simulations.
• 2.  
  Spatially resolved SFHs computed using IFU data from surveys like SDSS-IV MaNGA and CALIFA allow us to better understand the correlation between the SFH and morphology and discriminate between inside-out versus outside-in scenarios for galaxy growth and quenching (Goddard et al. 2016), better examine the connection between the physical properties of individual regions within galaxies and their SFHs (Rowlands et al. 2018), and test scaling relations at different regimes (Hsieh et al. 2017). Care needs to be exercised in interpreting these results because we only see where the stellar populations are today. Additional kinematic information would help alleviate this problem to a certain extent.
• 3.  
  Correlating SFHs with environment, size, kinematics, central density, and morphology in addition to stellar mass, SFR, and redshift could help build a unified picture of how galaxies evolve, with the SFHs providing a link between galaxy populations of different types and their earlier progenitors. Although this approach is similar to empirical models (Behroozi et al. 2018; Moster et al. 2018), it has the advantage of much richer observational constraints from the individual SFHs of galaxies. A comparison of SFH distributions between simulations and observations would allow us to qualify additional factors that are not directly accessible.
• 4.  
  The Gaussian process–based parameterization can also be used as a general compression method for compressing and storing PDFs from all kinds of codes, similar to Malz et al. (2018).

While we have performed an extensive range of validation tests (Section 3) to ensure that all the quantities reported in this work are robust, there are some caveats to keep in mind while extending the SED fitting to different data sets or the analysis beyond what is performed here.

• 1.  
  SFHs are not mass accretion histories. The SFH is a record of when the stars present in a galaxy at the time of observation were formed, as opposed to the mass accretion history, which is a record of when those stars entered the galaxy. These two quantities are the same for stars formed in situ, but differ when the stars were brought in through mergers. This needs to be taken into account in certain kinds of analysis; for example, when using a mass- and redshift-dependent merger fraction to correct for mass functions calculated by propagating galaxies backward in time along their SFHs.
• 2.  
  Lack of sensitivity to the shortest timescales. The smooth SFHs reconstructed using SED fitting in this work cannot capture starbursts that can happen on extremely short timescales of $\sim { \mathcal O }(10)\,\mathrm{Myr}$. While fitting galaxies from the semi-analytic model that contain such starbursts, we generally find that the overall stellar mass is well-recovered, but the starburst is smeared out over larger timescales, depending on when it occurred. This needs to be accounted for in the uncertainty budget; for example, while calculating the scatter along the SFR–M_* correlation using galaxies propagated backward in time along their SFR–M_* trajectories.
• 3.  
  Non-uniform sensitivity to variations in SFH. SED fitting is more sensitive to recent star formation than it is to star formation older than a few Gyr. As we go back in time, our SFH reconstruction transitions from being likelihood-dominated to being prior-dominated; while we show that this does not cause biases in our SFH reconstruction in Section 3, it does mean that we are less sensitive to sharp variations in SFH at large lookback times than when closer to the time of observation.
• 4.  
  Correlation with chemical enrichment histories: While we have considered the problem of estimating the SFHs of galaxies in this work, in practice they are highly correlated with the chemical enrichment histories of galaxies. Although metallicity is poorly constrained in our current observations, when working with higher-S/N data or spectra, the analysis should include a joint model for SFR(t) and Z(t), which can be achieved through joint priors on the metallicity given by Z(M_*, {t_x}), informed using simulations.