# A census of star formation histories of massive galaxies at 0.6 < z < 1 from
spectro-photometric modeling using Bagpipes and Prospector

Yasha Kaushal Department of Physics and Astronomy and PITT PACC, University of
Pittsburgh, 3941 O’Hara St, Pittsburgh, PA 15213, USA  Angelos Nersesian
Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, B-9000
Gent, Belgium  Rachel Bezanson Department of Physics and Astronomy and PITT
PACC, University of Pittsburgh, 3941 O’Hara St, Pittsburgh, PA 15213, USA
Arjen van der Wel Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan
281 S9, B-9000 Gent, Belgium  Joel Leja Department of Astronomy &
Astrophysics, The Pennsylvania State University, University Park, PA 16802,
USA  Institute for Computational & Data Sciences, The Pennsylvania State
University, University Park, PA 16802, USA  Institute for Gravitation and the
Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
Adam Carnall Institute for Astronomy, School of Physics and Astronomy,
University of Edinburgh, Royal Observatory, Edinburgh, EH9 3HJ, UK  Anna
Gallazzi INAF-Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5,
I-50125 Firenze, Italy  Stefano Zibetti INAF-Osservatorio Astrofisico di
Arcetri, Largo Enrico Fermi 5, I-50125 Firenze, Italy  Gourav Khullar
Department of Physics and Astronomy and PITT PACC, University of Pittsburgh,
3941 O’Hara St, Pittsburgh, PA 15213, USA  Marijn Franx Leiden Observatory,
Leiden University, P.O.Box 9513, NL-2300 AA Leiden, The Netherlands  Adam
Muzzin Department of Physics and Astronomy, York University, 4700 Keele St.,
Toronto, Ontario, M3J 1P3, Canada  Anna de Graaff Leiden Observatory, Leiden
University, P.O.Box 9513, NL-2300 AA Leiden, The Netherlands  Max-Planck-
Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany  Camilla
Pacifici Space Telescope Science Institute, 3700 San Martin Drive, Baltimore,
MD 21218, USA  Katherine E. Whitaker Department of Astronomy, University of
Massachusetts, Amherst, MA 01003, USA  Cosmic Dawn Center (DAWN), Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200,
Denmark  Eric F. Bell Department of Astronomy, University of Michigan, Ann
Arbor, MI 48109, USA  Marco Martorano Sterrenkundig Observatorium,
Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium

###### Abstract

We present individual star-formation histories of $\sim 3000$ massive galaxies
(log($\mathrm{M_{*}/M_{\odot}}$) > 10.5) from the Large Early Galaxy
Astrophysics Census (LEGA-C) spectroscopic survey at a lookback time of
$\sim$7 billion years and quantify the population trends leveraging 20hr-deep
integrated spectra of these $\sim$ 1800 star-forming and $\sim$ 1200 quiescent
galaxies at 0.6 < $z$ < 1.0. Essentially all galaxies at this epoch contain
stars of age < 3 Gyr, in contrast with older massive galaxies today,
facilitating better recovery of previous generations of star formation at
cosmic noon and earlier. We conduct spectro-photometric analysis using
parametric and non-parametric Bayesian SPS modeling tools - Bagpipes and
Prospector to constrain the median star-formation histories of this mass-
complete sample and characterize population trends. A consistent picture
arises for the late-time stellar mass growth when quantified as $t_{50}$ and
$t_{90}$, corresponding to the age of the universe when galaxies formed 50%
and 90% of their total stellar mass, although the two methods disagree at the
earliest formation times (e.g. $t_{10}$). Our results reveal trends in both
stellar mass and stellar velocity dispersion as in the local universe - low-
mass galaxies with shallower potential wells grow their stellar masses later
in cosmic history compared to high-mass galaxies. Unlike local quiescent
galaxies, the median duration of late-time star-formation ($\tau_{SF,late}$ =
$t_{90}$ \- $t_{50}$) does not consistently depend on the stellar mass. This
census sets a benchmark for future deep spectro-photometric studies of the
more distant universe.

galaxies: evolution, galaxies: formation, galaxies: stellar content, galaxies:
star formation, galaxies: statistics

tablenum SIXtablenum

##  1 Introduction

Galaxies are complex systems of multiple stellar populations. Observational
constraints on the timescales of star-formation and hence stellar mass growth
and assembly enable us to understand the role of various physical mechanisms
and their environment in guiding the formation and cosmic evolution of
galaxies (Panter et al., 2007; Leitner, 2012). Understanding these processes
on various physical scales (radial and spatially resolved) and temporal scales
also helps shed light on some of the long-standing puzzles, such as when and
how galaxies stop forming stars (quenching) (Schawinski et al., 2014; Wild et
al., 2016; Schreiber et al., 2016; Maltby et al., 2018; Carnall et al., 2018)
and the transition from star-forming to quiescent that leaves a bimodal
population. (Bell et al., 2004; Muzzin et al., 2013a; Vulcani et al., 2014;
Smethurst et al., 2015; Taylor et al., 2015).

Our current knowledge of the formation and cosmic evolution of galaxies is
fundamentally limited by trade-offs between the quality and quantity of
observations. Building a coherent picture of galaxy evolution demands (a)
looking back in time for galaxies in the younger universe without compromising
the quality of observations, and (b) a statistically large representative
sample size to connect the dots between similar populations of galaxies at
multiple epochs. This is difficult owing to technological and observational
challenges as light from galaxies at increasing distances dims and shifts to
less accessible infrared regime.

As stellar sources predominantly emit in UV/optical to NIR wavelengths,
broadband multi-wavelength spectral energy distribution (SED) of a galaxy
contains information about its total stellar mass and dust reddening. Stellar
Population Synthesis (SPS) modeling is a widely used technique to derive
physical properties (e.g., age, mass, dust, metallicity) of a composite
stellar population (often a galaxy). This modeling is often limited to
broadband photometric data due to its relative availability and the speed and
simplicity of modeling a handful of data points compared to modeling a higher
resolution spectrum. However, broadband photometry itself may not be
sufficient to completely break the dust-age-metallicity degeneracy or to
constrain higher order moments of star formation histories (SFHs) (Leja et al.
(2017, 2019); Tacchella et al. (2022), A. Nersesian et al. submitted). The
numerous, old low mass stars have long lifetimes and slow spectral evolution,
yet leave weak imprints on observed data owing to their faint intrinsic
luminosities. Therefore, inferring SFHs from photometry only is strongly
susceptible to even small perturbations in data, leading to different recovery
of SFHs (see e.g. Ocvirk et al. (2006)). The recovered SFHs from mocks show
biases up to 0.2 dex in the mass-weighted formation times. This is even worse
with real-world broadband photometric data having more systematic
uncertainties involved (Leja et al., 2019; Carnall et al., 2019a). Photometry
alone is not always sufficient to differentiate amongst modeling assumptions
(Belli et al., 2018; Carnall et al., 2019a; Tacchella et al., 2022), thus
properties measured from SED-only fits will always suffer from systematic
uncertainties e.g. flux calibrations and different physical models obtained
from stellar templates that vary with uncertainties on the data.

To robustly recover and analyze the timescales of a galaxy’s stellar mass
growth, we need to decipher spectral signatures containing strong imprints of
stellar evolution. Continuum spectroscopy in the rest-frame optical regime
contains information about the nature of past star-formation (bursty, uniform,
rising, declining etc.) and metal enrichment within a galaxy in features like
G4300, Fe4383, Fe4531, $Mg_{2}$, Balmer lines and
4000-$\mathrm{\SIUnitSymbolAngstrom}$ break. These signatures evolve most
rapidly in young stellar populations, since more massive galaxies today are
generally older by several Gyr, studies of their SFHs should be more accurate
and precise at earlier times. To obtain better signal-to-noise ratios, many
analyses of high-redshift galaxy surveys involve mass matched stacking of
spectra, vanishing any information about any individual galaxy’s evolution
(Schiavon et al., 2006; Choi et al., 2014; Siudek et al., 2017; Cullen et al.,
2019). High signal-to-noise continuum spectroscopic data when combined with
deep broadband photometric data has the potential to produce much stronger
constraints on SFHs of galaxies (Gallazzi et al., 2008; Pacifici et al., 2012;
Thomas et al., 2017; Carnall et al., 2019b; Iyer, 2019; Tacchella et al.,
2022; Webb et al., 2020) and provide clues on the number of major star
formation episodes and the timescales of rejuvenation, starbursts, and
quiescence. Spectroscopy alone also suffers from systematics like instrumental
noise and outlier pixels and emphasizes the importance of spectro-photometric
modeling.

Figure 1: Stellar mass vs redshift (left), median signal to noise of LEGA-C
spectra (center) and stellar velocity dispersion (right) for all 3005 primary
LEGA-C Galaxies from the DR3 catalog with density distributions of star-
forming (blue) and quiescent (red) populations based on rest-frame UVJ colors.
Our effective sample selections are shown with dashdot lines (mass-selected)
and solid lines (velocity dispersion-selected).

It is worth noting that even spectro-photometric modeling can be
insufficiently sensitive to more slowly evolving old stellar populations, more
so at later epochs. However, it has been argued that photometry can only probe
the last $\sim$1 Gyr of a SFH, whereas adding spectroscopic data can help to
probe the SFH further back in time (Chaves-Montero & Hearin, 2020). Many
studies have been conducted to recover SFHs of local galaxies (Thomas et al.,
2005; Cid Fernandes, 2007; Panter et al., 2007; Tojeiro et al., 2009; McDermid
et al., 2015; Citro et al., 2016; Ibarra-Medel et al., 2016). However, this
approach of using fossil records at $z$ = 0 may not lead to reliable results,
especially for more massive galaxies for which the major episode of their
star-formation happened very early in the universe and these galaxies are now
left with predominantly old stellar populations that suffer from strong
outshining effects. This underscores the importance deep spectroscopic studies
with high signal-to-noise at large lookback times.

Recent advances in computational and sampling techniques (Skilling, 2006;
Feroz & Hobson, 2008; Feroz & Skilling, 2013; Feroz et al., 2019) have led to
the development of a variety of tools that provide fast full-Bayesian fitting
of models to spectro-photometric data, like MCSED (Bowman et al., 2020),
BEAGLE (Chevallard & Charlot, 2016), Prospector (Johnson et al., 2021; Johnson
& Leja, 2017), and Bagpipes (Carnall et al., 2018). These Bayesian methods
give a better handle on the priors assumed and are robust to classical problem
of over-fitting the data with complex models e.g. (Leja et al., 2019).
However, different SSP libraries and modeling assumptions in these tools can
influence the derived star-formation histories (Martins, 2021; Pacifici et
al., 2023).

Though real SFHs of galaxies are complex, a simplified way to model them is by
parametrizing the SFHs using a functional form. The most commonly used forms
are exponentially declining (Mortlock et al., 2017; McLure et al., 2018; Wu et
al., 2018a), delayed exponentially declining (Ciesla et al., 2017; Chevallard
et al., 2019), log-normal (Gladders et al., 2013; Abramson et al., 2015;
Diemer et al., 2017; Cohn, 2018) and double-power law (Ciesla et al., 2017;
Carnall et al., 2018). These analytic prescriptions have been shown to match
well with many SFHs from simulations (Simha et al., 2014; Diemer et al., 2017)
and are widely used as they minimize computational requirements. Increasing
complexity by adding bursts of star-formation can bring parametric SFHs even
closer to realistic scenarios. However, capturing events like rejuvenation and
sudden quenching can still be challenging for parametric models. A more
flexible method is non-parametric SFHs - which adopts a series of periods of
constant star-formation in fixed or flexible time bins (Chauke et al., 2018;
Cappellari, 2017; Leja et al., 2017; Iyer & Gawiser, 2017; Ocvirk et al.,
2006; Cid Fernandes et al., 2005). Non-parametric approaches have higher
flexibility with a wider range of possible solutions and therefore can allow a
broader range of priors on SFHs and decrease biases on recovered results with
more realistic episodes of star-formation in predefined time bins (Iyer &
Gawiser, 2017; Leja et al., 2019; Suess et al., 2022a, b).

Until recently, SPS modeling of statistically representative populations has
been limited to the analysis of photometry-only data sets. Numerous efforts
have been put toward recovering the SFHs of both star-forming and quiescent
populations from deep broadband photometric surveys spanning wide redshift
ranges (Dye, 2008; Wuyts et al., 2009; Pforr et al., 2012; Pacifici et al.,
2016b, a; Iyer & Gawiser, 2017; Iyer, 2019; Aufort et al., 2020; Olsen et al.,
2021) due to greater data availability and lower computational demands. In
contrast, fewer studies have tested modeling spectro-photometric data for high
redshift galaxies (Carnall et al., 2019b; Estrada-Carpenter et al., 2020;
Forrest et al., 2020; Johnson et al., 2021; Tacchella et al., 2022; Khullar et
al., 2022; Hamadouche et al., 2023) owing to the dearth of high signal-to-
noise continuum spectroscopy at significant lookback times. Large
spectroscopic surveys like DEEP2 (Newman et al., 2013), MOSDEF (Kriek et al.,
2015) and KBSS (Rudie et al., 2012; Steidel et al., 2014; Strom et al., 2017)
have been primarily sufficient to characterize emission line properties of
thousands of star-forming galaxies. On the other hand, smaller spectroscopic
studies of hundreds of quiescent galaxies have provided windows into the
stellar populations of quiescent systems out to z$\sim$2, but with a
significant bias towards the brightest, most massive subset.

The Large Early Galaxy Astrophysics Census (LEGA-C) (van der Wel et al., 2021;
Straatman et al., 2018; van der Wel et al., 2016) provides a novel opportunity
to characterize the full population of Milky Way-mass and larger progenitors
at significant lookback time. This deep spectroscopic survey of $\sim$3000
star-forming and quiescent galaxies looking 6-8 billion years back in time
(0.6 $<$ z $<$ 1) includes deep imaging available for each galaxy from the
UltraVISTA survey (Muzzin et al., 2013b). At this redshift, most stars in
LEGA-C survey galaxies are $<$3 billion years in age, enabling more robust
characterization of their SFHs.

The primary goal of this paper is to measure the SFHs of the full LEGA-C DR3
sample by applying two commonly used Bayesian SPS modeling techniques (with
modeling choices optimized for each) on spectro-photometric data and
investigate the evolution of massive galaxies before z$\sim$0.8. We quantify
these SFHs in two widely used metrics, namely $t_{50}$ and $t_{90}$,
corresponding to the times when a galaxy formed 50% and 90% of its total
stellar mass respectively and study the population trends of these formation
times with stellar mass and stellar velocity dispersion.

The structure of the paper is as follows. Section 2 describes the LEGA-C
dataset and our approach to spectro-photometric modeling using Bagpipes and
Prospector, and some example demonstration of modeling results. Section 3
describes the SFHs of star-forming and quiescent galaxies and of the full
population as recovered from spectro-photometric fits. We show cumulative
median mass growth trends in stellar mass bins and further quantify population
trends of $t_{50}$ and $t_{90}$ with stellar mass and stellar velocity
dispersion. In Section 4 we expand on the interpretation of our results with
respect to previous low and high redshift studies and how this impacts our
understanding of the formation of both star-formation and quiescent systems.
In Section 5 we conclude our study and highlight some major takeaways from
this analysis. We also speculate on potential future works that could be a
successor of this study to help us better constrain the evolution of massive
galaxies and answer more broader questions. Throughout this paper we assume
$\Omega_{M}$ = 0.3, $\Omega_{\Lambda}$ = 0.7, $H_{0}$ = 70 km s-1 Mpc-1 and
all magnitudes are in quoted in the AB system.

##  2 Data and Modeling Methods

###  2.1 Data and Sample

The Large Early Galaxies Astrophysics Census (LEGA-C) (van der Wel et al.,
2021; Straatman et al., 2018; van der Wel et al., 2016) is a 130-night public
spectroscopic survey of $\sim$3000 Ks-band selected galaxies targeting
redshift range 0.6 < $z$ < 1.0 in the COSMOS field, looking 6-8 billion years
back in time. Each spectrum has an approximate observed spectral coverage of
6300$\mathrm{\SIUnitSymbolAngstrom}$ \- 8800$\mathrm{\SIUnitSymbolAngstrom}$
corresponding to 3600$\mathrm{\SIUnitSymbolAngstrom}$ \-
5200$\mathrm{\SIUnitSymbolAngstrom}$ rest-frame optical regime. The survey was
conducted on Very Large Telescope (VLT) using ViMOS (VIsible Multi Object
Spectrograph) (Le Fevre et al., 2000) and completed its third and final data
release (DR3) in August 2021 (van der Wel et al., 2021).

The full spectroscopic sample consists of 4081 galaxies - 3029 primary targets
and 1052 fillers. Targets were Ks-band selected from the UltraVISTA catalog
(Muzzin et al., 2013b) to include massive galaxies above $\sim$
$10^{10}M_{\odot}$ and to have Ks-band magnitudes brighter than a redshift
dependent limit of [20.7 - 7.5log((1+z)/1.8)]. This selection criterion is
independent of any derived quantities from the spectra. The full primary
sample can be re-weighted using sampling and volume corrections following van
der Wel et al. (2021) to be mass representative above
log($\mathrm{M_{*}/M_{\odot}}$) > 10.5, which we adopt as the mass threshold
for this study. Each galaxy was observed for $\sim$ 20 hours yielding S/N
$\approx$ 20 $$\mathrm{\SIUnitSymbolAngstrom}$^{-1}$ continuum with high
fidelity absorption and emission line features for both dusty, blue as well as
faint, red galaxies.

As the LEGA-C survey targets only a subset of the full photometric sample, we
account for missing galaxies with appropriate weighting to individual galaxies
to ensure spectroscopic completeness above the previously specified magnitude
limit in a full census of the properties of massive galaxies at that epoch.
Hence throughout this paper we apply a multiplicative factor Tcor \-
corresponding to volume and sample correction - from the DR3 catalog (van der
Wel et al., 2021) to individual object counts to make it representative of
that redshift. The photometric information is taken from UltraVISTA catalog
(Muzzin et al., 2013b). We use BvrizYJ bands for our primary spectro-
photometric analysis.

The stellar mass (log($\mathrm{M_{*}/M_{\odot}}$)), spectroscopic redshift
(z), signal-to-noise (S/N) and integrated stellar velocity dispersion
($\sigma_{*}$) density distributions for the primary sample are shown in Fig.
1. All values are taken from the LEGA-C DR3 catalog with (total) stellar mass
estimates from Prospector photometry-only stellar populations fits
(LOGM_MEDIAN), spectroscopic redshift from the LEGA-C spectra (Z_SPEC), the
median S/N per pixel of the LEGA-C spectrum (S/N) and stellar velocity
dispersion (SIGMA_STARS) estimated from Gaussian broadening of theoretical
single stellar population models as described in Bezanson et al. (2018). We
split the sample into star-forming and quiescent populations using Muzzin et
al. (2013c) rest-frame U-V and V-J color-color cuts, identifying quiescent
galaxies with U-V > (V-J) $\times$ 0.88 + 0.69. Note that throughout this
study, we adopt the LEGA-C DR3 catalog value of $\sigma_{*}$ and
$M_{*}/M_{\odot}$ to maintain a single set of labels for each individual
galaxy and to facilitate consistent comparisons with other studies. There are
1774 star-forming (medians: log($\mathrm{M_{*}/M_{\odot}}$) $\sim$ 10.8 and
$\sigma_{*}$ $\sim$ 140 km s-1) and 1231 quiescent galaxies (median:
log($\mathrm{M_{*}/M_{\odot}}$) $\sim$ 11.2, $\sigma_{*}$ $\sim$ 200 km s-1)
in the primary sample of 3005 objects shown in red and blue colors
respectively (obtained using PRIMARY flag in LEGA-C DR3 catalog, more details
in van der Wel et al. (2021)). Throughout this paper, we focus on objects
above the approximate mass completeness limit of the LEGA-C survey (10.5 <
log($\mathrm{M_{*}/M_{\odot}}$) $\leq 12$ and 100 < $\sigma_{*}$ $\leq 300$ km
s-1). For mass-limited sample we have an effective sample size of 2703 unique
galaxies (1459 star-forming and 1244 quiescents) and the velocity-dispersion
sample includes 2823 unique galaxies (1575 star-forming and 1248 quiescents).
This sample selection is shown with dashed and solid lines in Fig 1.

Table 1: Bagpipes modeling parameters and prior distribution functions  No | Property | Parameters | Symbol/Unit | Prior | Range | Remark
---|---|---|---|---|---|---
1 | SFH | alpha | $\alpha$ | logarithmic | (0.01,1000) | falling slope
2 | (double-power law) | beta | $\beta$ | logarithmic | (0.01,1000) | rising slope
3 |  | tau |  $\tau$/Gyr | uniform | (0.1,15) | peak of star-formation
4 |  | stellar metallicity | $Z_{*}/Z_{\odot}$ | uniform | (0.02,2.5) | From grid interpolation
5 |  | stellar mass | $M_{*}/M_{\odot}$ | logarithmic | (0,13) | priors from Carnall et al. (2019b)
6 | Burst-1 | age | Gyr | uniform | (0,13) | delta function burst
7 |  | stellar mass | $M_{*}/M_{\odot}$ | logarithmic | (0,10) | with 3 free parameters
8 |  | stellar metallicity | $Z_{*}/Z_{\odot}$ | uniform | (0,2.5) |
9 | Burst-2 | age | Gyr | uniform | (0,13) | same as above
10 |  | stellar mass | $M_{*}/M_{\odot}$ | logarithmic | (0,10) |
11 |  | stellar metallicity | $Z_{*}/Z_{\odot}$ | uniform | (0,2.5) |
12 | Dust | V-band attenuation | Av / mag | uniform | (0.,2.0) | Charlot & Fall (2000)
13 |  | slope of attenuation | n | Gaussian | (0.3,2.5) |  $\mu$ = 0.7, $\sigma$ = 0.3
14 | Stellar Velocity Dispersion | sigma |  $\sigma_{\star}$/(km/s) | logarithmic | (40,400) | free parameter
15 | Redshift | z | $z_{spec}$ | LEGA-C |  | fixed parameter
16 | Spectral White Noise | scaling factor | a | logarithmic | (0.1,10) | uncorrelated spectroscopic noise
17 | Calibration | 0th order | c0 | Gaussian | (0.9,1.1) |  $\mu=1$, $\sigma=0.05$
18 | Calibration | 1st order | c1 | Gaussian | (-0.5,0.5) |  $\mu=0$, $\sigma=0.1$
19 | Calibration | 2nd order | c2 | Gaussian | (-0.5,0.5) |  $\mu=0$, $\sigma=0.1$
Table 2: Bagpipes  Spectro-photometric modeling results

ID | Mask | LEGAC_ID | $z_{spec}$ | $\sigma_{*,cat}$ | $log(M_{*,cat}/M_{\odot})$ | $log(M_{*,fit}/M_{\odot})$ | SFR | t50 | t90
---|---|---|---|---|---|---|---|---|---
|  |  |  | (km/s) |  |  | ($M_{\odot}/yr$) | (Gyr) | (Gyr)
103041 | 17 | 1159 | 0.82 | 75.9 | 10.52 | 10.63$\pm_{0.07}^{0.1}$ | 5.91$\pm_{1.14}^{2.02}$ | 5.07$\pm_{0.18}^{0.24}$ | 6.01$\pm_{0.14}^{0.13}$
103061 | 16 | 1160 | 0.72 | 193.1 | 11.03 | 11.25$\pm_{0.03}^{0.03}$ | 0.00$\pm_{0.00}^{0.00}$ | 3.54$\pm_{0.03}^{0.24}$ | 5.72$\pm_{0.06}^{0.03}$
103155 | 19 | 1161 | 0.64 | 103.1 | 11.07 | 11.36$\pm_{0.06}^{0.07}$ | 17.26$\pm_{2.54}^{3.67}$ | 2.79$\pm_{0.03}^{0.03}$ | 6.29$\pm_{0.01}^{0.00}$
103179 | 16 | 1162 | 0.62 | 108.1 | 10.75 | 11.02$\pm_{0.04}^{0.04}$ | 0.00$\pm_{0.00}^{0.00}$ | 3.22$\pm_{0.29}^{0.41}$ | 5.18$\pm_{0.15}^{0.20}$
103274 | 14 | 1163 | 0.92 | 175.9 | 11.06 | 11.43$\pm_{0.04}^{0.04}$ | 6.12$\pm_{0.47}^{0.74}$ | 1.02$\pm_{0.12}^{0.30}$ | 3.56$\pm_{0.26}^{0.52}$

**footnotetext: We display a table snippet for formatting; the full data table
is available to download in the MRT format.

###  2.2 Spectro-photometric Modeling

In this work, we use two state-of-the-art Bayesian SPS modeling tools,
Bagpipes (Carnall et al., 2018) and Prospector (Johnson & Leja, 2017; Johnson
et al., 2021), to perform a full spectro-photometric fitting of the primary
LEGA-C galaxies. Details of the two methods used in this work are discussed in
sections 2.2.1 and 2.2.2. We use BvrizYJ bands in our primary spectro-
photometric analysis following the conclusions in van der Wel et al. (2021)’s
Appendix B. That study finds that using photometry longward of rest-frame
$\sim$0.8 $\mu$m disagrees with popular SPS models and leads to systematic
errors in fits to the SED of $\sim$20% of their sample, which in turn
propagated into the stellar mass, star formation rate, dust attenuation, and
other population parameters. We check the robustness of our results by
performing comparisons with two more parametric set of fits altering the
photometric data - (1) spectra +vriz bands, and (2) spectra + 20 bands (B to
24 micron, namely - B, g, IA484, IA527, V, IA624, r, IA679, IA738, i, z, y, J,
H, Ks, ch1, ch2, ch3, ch4, mips24) covering 15 Subaru bands (B to Ks), 4 IRAC
bands (ch1 to ch4) and 1 Spitzer MIPS24 band, ranging from UV to far-IR
wavelengths. Results of this comparison are shown in Figure A1 in our Appendix
A.

####  2.2.1 Bagpipes

Bagpipes (Bayesian Analysis of Galaxies for Physical Inference and Parameter
EStimation) is a SPS modeling package built on the updated BC03 (Bruzual &
Charlot, 2003) spectral
library111www.bruzual.org/~gbruzual/bc03/Updated_version_2016 with the 2016
version of the MILES library of empirical spectra that includes 2.5
$\mathrm{\SIUnitSymbolAngstrom}$ resolution in
3525$\mathrm{\SIUnitSymbolAngstrom}$ \- 7500$\mathrm{\SIUnitSymbolAngstrom}$
wavelength range (Falcón-Barroso et al., 2011). It is built on a Kroupa (2001)
IMF assumption and utilizes a Multi-Nest nested sampling algorithm (Feroz et
al., 2019) to produce posterior distributions of physical parameters. We
perform parametric full spectral SPS modeling of LEGA-C spectra with BvrizYJ
UltraVISTA broadband photometry using a double-power law SFH (Ciesla et al.,
2017; Carnall et al., 2018), parameterized as:

$\text{SFH(t)}\propto\left[\left(\frac{t}{\tau}\right)^{\alpha}+\left(\frac{t}{\tau}\right)^{-\beta}\right]^{-1}$
(1)

Table 3: Prospector modeling parameters and prior distribution functions  No | Parameter | Description | Prior
---|---|---|---
1 | $\sigma_{*}$ | Velocity dispersion | fixed to LEGA-C values
2 | $\log\left(M/\mathrm{M}_{\odot}\right)$ | total stellar mass formed | uniform in log space: min=7, max=12
3 | $\log\left(Z/\mathrm{Z}_{\odot}\right)$ | stellar metallicity | uniform in log space: min=-1.98, max=0.4
4 | SFR ratios | Ratios of adjacent SFRs | Student’s t-distribution ($\sigma=0.3$, $\nu=2$)
5 | $z_{spec}$ | redshift | prior: LEGA-C spectroscopic redshift $\pm 0.005$
6 | $\hat{\tau}_{\lambda,2}$ | diffuse dust optical depth | clipped normal: min=0, max=4, mean=0.3, $\sigma=1$
7 | $\hat{\tau}_{\lambda,1}$ | birth-cloud dust optical depth | clipped normal in ($\hat{\tau}_{\lambda,1}/\hat{\tau}_{\lambda,2}$):
|  |  | min=0, max=4, mean=0.3, $\sigma=1$
8 | $n$ |  Kriek & Conroy (2013) dust law slope | uniform: min=-1, max=0.4
9 | $\log\left(Z_{\mathrm{gas}}/\mathrm{Z}_{\odot}\right)$ | Gas-phase metallicity | uniform: min=-2, max=0.5
10 | $\log U$ | Ionization parameter | uniform: min=-4, max=-1
11 | $\sigma_{\mathrm{eline}}$ | Emission line amplitude | uniform: min=30, max=300
12 | $\nu$ | Spectral white noise | uniform: min=1, max=3
13 | $cp$ | Photometric calibration | uniform: min=$10^{-5}$, max=0.5
14 | $cs$ | Spectroscopic calibration | uniform: min=$10^{-5}$, max$=0.5$
Table 4: Prospector  Spectro-photometric modeling results

ID | Mask | LEGAC_ID | $z_{spec}$ | $\sigma_{*,cat}$ | $log(M_{*,cat}/M_{\odot})$ | $log(M_{*,fit}/M_{\odot})$ | SFR | t50 | t90
---|---|---|---|---|---|---|---|---|---
|  |  |  | (km/s) |  |  | ($M_{\odot}/yr$) | (Gyr) | (Gyr)
103041 | 17 | 1159 | 0.82 | 75.9 | 10.52 | 10.64$\pm_{0.06}^{0.05}$ | 1.27$\pm_{0.47}^{8.80}$ | 3.72$\pm_{0.09}^{0.11}$ | 5.94$\pm_{0.43}^{0.03}$
103061 | 16 | 1160 | 0.72 | 193.1 | 11.03 | 11.37$\pm_{0.02}^{0.02}$ | 0.15$\pm_{0.12}^{2.50}$ | 2.40$\pm_{0.01}^{0.01}$ | 4.66$\pm_{0.19}^{0.36}$
103155 | 19 | 1161 | 0.64 | 103.1 | 11.07 | 11.19$\pm_{0.02}^{0.03}$ | 14.34$\pm_{2.01}^{15.61}$ | 4.32$\pm_{0.66}^{0.52}$ | 6.70$\pm_{0.02}^{0.34}$
103179 | 16 | 1162 | 0.62 | 108.1 | 10.75 | 11.04$\pm_{0.02}^{0.02}$ | 0.19$\pm_{0.13}^{2.05}$ | 2.60$\pm_{0.17}^{0.11}$ | 5.27$\pm_{0.26}^{0.27}$
103274 | 14 | 1163 | 0.92 | 175.9 | 11.06 | 11.35$\pm_{0.04}^{0.04}$ | 3.82$\pm_{1.80}^{7.41}$ | 2.36$\pm_{0.11}^{0.01}$ | 4.83$\pm_{0.36}^{0.04}$

**footnotetext: We display a table snippet for formatting; the full data table
is available to download in the MRT format.

This model has three free parameters describing the rising ($\beta$), falling
($\alpha$) and peak ($\tau$) of star formation, whereas other widely used
options (e.g. exponential, delayed-tau and log-normal) have two or less free
parameters. Hence, by construction, the double-power law SFH has more
flexibility. Tau models are shown to fail to recover mock SFHs (Carnall et
al., 2019a) and simulation SFHs (Pacifici et al., 2012). A double-power law
SFH is chosen owing to its ability to recover the redshift evolution of cosmic
SFRD (Behroozi et al., 2013; Gladders et al., 2013; Madau & Dickinson, 2014)
and owing to the agreement with simulation results (Pacifici et al., 2016b;
Diemer et al., 2017). The rising and falling slopes essentially do not change
beyond the prior ranges chosen (become either flat or vertical) owing to the
analytical functional form (Eq.1) and cover full variations of possible SFHs
within this limit. The stellar metallicity is assumed to be the same for all
stars born. This value is linearly interpolated on a grid of SSP models and is
allowed as a free parameter varying from (0.02,2.5)$Z_{\odot}$ uniformly in
linear space. The $Z_{\odot}$ value is assumed to be 0.02 in BC03 models. We
test the impact of choosing linear and logarithmic priors on the stellar
metallicity on the derived posterior values in our spectro-photometric fits
for a subset of galaxies and find no strong deviations. The last free SFH
parameter is the stellar mass formed in the entire lifetime of a galaxy until
the point of observation (without mass return to the ISM); we allow
$\mathrm{log(M_{*}/M_{\odot})}$ to vary logarithmically from (0, 13). Two
bursts on top of a double power law are included to account for any abrupt
variation in star formation activity. Each burst is given flexibility of age,
stellar mass and stellar metallicity and hence three free parameters. We use
Charlot & Fall (2000) dust model with two free parameters - V-band attenuation
and the slope of attenuation. We adopt a second order spectral calibration and
white uncorrelated noise for spectral pixels, for which a detailed description
can be found in Carnall et al. (2019b). Dust emission models from Draine & Li
(2007) are implemented with fixed $Q_{pah}$ = 2, $U_{min}$ = 1 and
$\gamma_{e}$ = 0.01. While testing multiple parameter options and analyzing
the full posterior distributions of the output stellar population properties,
we found that nebular emission line modeling in Bagpipes biased the stellar
metallicities of star-forming galaxies to high values
($\mathrm{log(Z_{*}/Z_{\odot})}$ > 0.35). Additionally, we were concerned that
the limiting ionizing radiation of young stars could inappropriately describe
emission lines produced by either evolved stars or AGN (and bias the SFR
estimates), especially in such a diverse and massive sample of galaxies
(Carnall et al., 2019b). Thus, we choose to mask the emission lines from the
fit. Stellar velocity dispersion is another free parameter modeled with a
variable Gaussian kernel in velocity space. Although the DR3 LEGA-C spectra
used in this study are flux calibrated using the UltraVISTA photometry (van
der Wel et al., 2021), we include an additional polynomial function of
wavelength to address any higher order spectro-photometric calibration
uncertainties. We use a second order Chebyshev polynomial function with
Gaussian priors - 0th order centered around unity and 1st and 2nd orders
centered around zero (see §4.3.1 in Carnall et al. (2019b) for more details).
This modeling requires an average of $\sim$70 CPU hours per galaxy. Further
description of all modeling parameters and prior distributions is shown in
Table 1. Modeling results are included in Table 2. We notice a subset (total
126 in the full sample / 4.2% of the full population) of quiescent population
that prefers very similar best-fitting SFHs with ages $t_{mw}\sim 3$Gyr and
3.160 $<$ $t_{lw}$ $<$ 3.163. We investigated these objects and found them
well constrained in parameter space, spanning a range in empirical properties
like spectral indices, UV-VJ colors and redshifts.

####  2.2.2 Prospector

Prospector is a Bayesian SPS modeling tool that allows a non-parametric
modeling of the SFH of a galaxy in piece-wise constant SFR time bins. (Johnson
et al., 2021). It deploys the Flexible Stellar Population Synthesis (FSPS)
package with the MILES stellar library and MIST isochrones to model stellar
properties (Conroy et al., 2009). We use a Chabrier (2003) IMF and the Kriek &
Conroy (2013) dust law with nebular continuum and line emissions modeled with
CLOUDY (Ferland et al., 2013). Note that this choice differs from the Bagpipes
modeling; we found that modeling nebular emission with Cloudy grids within
Prospector produced well behaved posterior distributions. We also tested the
impact of including/excluding physical modeling of emission lines in
Prospector on a a subset of total 300 galaxies (150 quiescent and 150 star-
forming) with significantly-detected emission lines and high signal-to-noise
spectra (OII EW > 4 and SN > 12). In general, the recovered SFHs agree for
each galaxy within uncertainties. In a small subset of quiescent galaxies
(N=11), non-physical modeling of the emission lines results in maximally old
stellar populations that are formed in dramatic, but uncertain, bursts of star
formation in the earliest time bin. Given the overall agreement between the
two sets of models and slightly more extended SFHs in the aforementioned
subset, we include the physical line ratio modeling in our Prospector fitting.
We use dynesty (Speagle, 2020) nested sampling option for posterior sampling,
similar to Bagpipes. This non-parametric approach is capable of recovering
complex SFHs and capturing abrupt star formation processes like sudden
quenching and rejuvenation events. On the flip side, fitting both the galaxy
spectra and SED is highly computationally expensive and requires about
$\sim$100 CPU hours per galaxy to fit an 8 fixed bin SFH model with 14-free
parameters.

In this work, we adopt a continuity prior piece-wise constant SFH with
Student’s-t distribution that fits the change in log(SFR(t)) in adjacent time
bins while weighing against abrupt changes in SFR(t). This prior has also been
shown to robustly reproduce mock and more importantly, simulated SFHs (Lower
et al., 2020). The pioneering work of Ocvirk et al. (2006) show that a maximum
of 8 episodes of SFH can be independently recovered from an optical spectra of
resolution R=10,000, S/N=100 and wavelength coverage $\lambda$ =
4000$\mathrm{\SIUnitSymbolAngstrom}$ \- 6800$\mathrm{\SIUnitSymbolAngstrom}$,
with the distinguishability of simple stellar populations proportional to the
separation in logarithmic time. Hence we use an 8 time bin SFH model (5
logarithmically spaced) in our analysis. The 8 bins of constant SFRs are
distributed as follows (in lookback time) -

0 < t < 30 Myr

30 Myr < t < 100 Myr

100 Myr < $t_{log1,2,3,4}$ < 0.85 $t_{univ}$ (5 log bins)

0.85 $t_{univ}$ < t < $t_{univ}$

The two most recent fixed bins capture signatures of any recent abrupt star
formation, one earliest fixed bin corresponding to oldest stellar populations
spans the first 15% of cosmic time, with 5 logarithmically-spaced bins in
between. Redshift is set to the LEGA-C spectroscopic redshift with allowed
$\pm$ 0.005 variation and stellar velocity dispersion is fixed to the LEGA-C
DR3 catalog values. A full description of the free and fixed parameters of the
Prospector model and their adopted priors are included in Table 3 and modeling
results are reported in Table 4.

One might be concerned that priors on the SFHs could drive differences in the
inferred SFHs that are derived from the two software packages. We test this by
drawing from the prior distributions in each set of models. Figure 2 depicts
the prior probability density from Bagpipes (orange) and Prospector (green)
when 1000 random samples are drawn from the SFH parametrization described in
Table 1 and Table 3 (more details in §2.2.1 and §2.2.2). For ease of
comparison, we assign a floor SFR value of 0.001 $M_{\odot}\;yr^{-1}$ to
arbitrarily small SFR values. The median distributions of prior SFHs are shown
in solid lines whereas 16th and 84th percentiles are shown with dotted lines.
The median values follow closely for both the codes, except at the earliest
times the analytic function requires SFR(t=0) = 0 whereas non-parametric
models assign non-zero star-formation in the earliest bin. This figure
suggests no strong biases in SFRs with lookback time from the priors adopted
in the two codes.

Figure 2: Prior density distribution of SFH obtained by drawing 1000 samples
from parametrization described in Table 1 (Bagpipes) and Table 3 (Prospector).
Thick lines show the medians (solid) and dashed 1-$\sigma$ (16th - 84th
percentile) scatter. The median values follow closely for both the codes,
except at the earliest times the analytic function requires SFR(t=0) = 0 and
non-parametric models assign non-zero star-formation in the earliest bin.

As shown in this section, we find the expected differences in the SFHs for
individual galaxies as derived by the two modeling methods. We explore the
impact of these choices on the full population of massive galaxies in the
remainder of the paper.

Figure 3: Spectro-photometric modeling (left) of four representative quiescent
galaxies using LEGA-C spectra and UltraVISTA photometry (grey) along with best
fit models from Bagpipes (orange) and Prospector (green). Lower sub-panels
show $\chi$ values (data-model/error) with grey bands showing regions masked
in parametric modeling. Reduced chi-squared values are quoted in inset boxes.
The right panels show median SFHs (solid lines) and 16th - 84th percentile
distributions of the posteriors (shaded regions). The vertical dashed lines
from left to right correspond to 90%, 50% and 10% formation time stamps for
each tool (color coded accordingly). The right most column shows the posterior
distributions of stellar mass (without mass loss), dust Av, and stellar
metallicity. Both modeling techniques yield good fits to the data and overall
the SFHs agree reasonably well in late times ($t_{90}$), with some divergence
at early times ($t_{10}$).

Figure 4: Four representative star-forming galaxies with modeling outputs in
the same format as described in Figure 3. In general, parametric and non-
parametric SFHs for individual star-forming galaxies agree reasonably well at
all times.

###  2.3 Modeling Results and Examples

First, we emphasize that all SFHs are allowed to start from the Big Bang,
however the analytic SFHs (Bagpipes) naturally exhibit more flexibility in
onset time (with negligible star-formation in early times for some cases) and
slope than the early bins in piece-wise constant non-parametric SFHs
(Prospector). Although in principle the first bin in the latter models could
exhibit negligible star formation, the fits prefer at least some non-zero
average SFR at the earliest times. Also, since bins represent SFRs averaged
over an extended period of time, they are more likely to be non-zero. This
discrepancy is partially driven by differences in modeling choices and also
likely reflects a lack of constraining power in the data at the earliest
times, even with the high S/N of the LEGA-C dataset. Fundamentally such
information is in the prior dominated regime; in this paper we quantify the
ultimate impact of popular choices on aggregate SFHs of the populations.

Fig. 3 and Fig. 4 show representative examples of spectro-photometric modeling
and recovered SFHs of four quiescent and four star-forming galaxies
respectively. In the left column, each panel shows an observed LEGA-C spectrum
and UltraVISTA photometry (grey), along with the best-fitting models from
Bagpipes (orange) and Prospector (green). The insets show the full BvrizYJ
photometric SEDs and models. We note that the models fit both the
spectroscopic and photometric data very well. The residuals are quantified in
$\chi$ values in the bottom panel ($\chi$ = (observed flux - model
flux)/error). Grey bands indicate regions that are masked to avoid emission
lines in the Bagpipes modeling. The right large panels show the corresponding
SFHs. Note that we chose to logarithmically scale the lookback time
(horizontal axis) to highlight the most robustly measured epochs at late
times.

For quiescent galaxies (Fig. 3), we note that both models provide excellent
fits to the observed data and median SFHs agree reasonably well, more so in
late times than early times, when quantified in the formation time metrics
$t_{10}$, $t_{50}$ and $t_{90}$ ($\Delta t_{10}|_{median}$ = 2.05 Gyr, $\Delta
t_{50}|_{median}$ = 1.47 Gyr and $\Delta t_{90}|_{median}$ = 0.22 Gyr with
uncertainties (calculated from posteriors) of order of 0.27/0.13/0.20 Gyr in
each delta metric respectively. Note that these are representative of the full
LEGA-C quiescent population statistics ($\Delta t_{10}|_{median}$ = 2.41 Gyr,
$\Delta t_{50}|_{median}$ = 1.38 Gyr and $\Delta t_{90}|_{median}$ = 0.23
Gyr). One subtle difference emerges at the earliest times when Prospector
consistently assigns finite star formation in the first time bin leading to
older stellar populations, whereas it is not necessarily the case for
parametric SFHs in Bagpipes. For star-forming galaxies (Fig. 4), we see better
agreement in median SFRs at later times compared to quiescent examples above,
with differences in timescales within the uncertainties. ($\Delta
t_{10}|_{median}$ = 0.49 Gyr, $\Delta t_{50}|_{median}$ = 0.28 Gyr and $\Delta
t_{90}|_{median}$ = 0.14 Gyr with uncertainties of order of 1.25/0.84/0.22 Gyr
in each delta metric respectively). Note that these are representative of the
full LEGA-C star-forming population statistics ($\Delta t_{10}|_{median}$ =
1.30 Gyr, $\Delta t_{50}|_{median}$ = 0.34 Gyr and $\Delta t_{90}|_{median}$ =
0.17 Gyr). Further mass and sigma trends of the full LEGA-C sample would be
analyzed in the following sections.

##  3 The Star Formation Histories of Massive Galaxies

Figure 5: Individual (left) and population (center and right) SFHs recovered
from Bagpipes and Prospector for two example samples of massive (11.5 $<$
log($\mathrm{M_{*}/M_{\odot}}$) $<$ 12) quiescent galaxies (top row (a)) and
star-forming galaxies (bottom row (b)). The middle panel shows individual SFHs
in the sSFR and the right panel shows 16th, 50th and 84th percentile
population distributions. Note that Bagpipes generates relatively higher
population scatter due to variation in the onset and truncation of the double
power law parametrization, in contrast with the early onset of star formation
inferred from non-parametric models.

In this section, we combine the median star-formation histories derived for
each galaxy in the LEGA-C sample to characterize the overall growth histories
of massive galaxies. We first interpolate all SFHs on a uniform age grid from
(0.01, 8) Gyr and impose a minimum star-formation rate (floor) of 0.001
$M_{\odot}\;yr^{-1}$ to arbitrarily small SFR values. For each galaxy, we take
the median SFH from the posterior distributions and combine these medians to
calculate population medians in bins of - (1) stellar mass and (2) stellar
velocity dispersions from LEGA-C DR3 catalog. For ease of presentation, we
first focus on three coarse stellar mass bins representative of LEGA-C primary
sample - a) 10.5 < log($\mathrm{M_{*}/M_{\odot}}$) < 11 b) 11 <
log($\mathrm{M_{*}/M_{\odot}}$) < 11.5 c) 11.5 <
log($\mathrm{M_{*}/M_{\odot}}$) < 12 and four coarse stellar velocity
dispersion bins - a) 100 < $\mathrm{\sigma_{*}}$(km/s) < 150 b) 150 <
$\mathrm{\sigma_{*}}$(km/s) < 200 c) 200 < $\mathrm{\sigma_{*}}$(km/s) < 250
and d) 250 < $\mathrm{\sigma_{*}}$(km/s) < 300 and later increase the
resolution to finer bins to characterize population trends. Note that the mass
and sigma values used to bin galaxies are taken from the LEGA-C DR3 catalog
(LOGM_MEDIAN,SIGMA_STARS)and are mass loss corrected, described in detail in
§2. These catalog values have mean offsets of +0.05 dex and -0.02 dex with
Bagpipes and Prospector modeling results respectively. In each mass/velocity
dispersion bin, each individual SFH is weighted by $T_{cor}$ (which includes
both a volume and sample correction factor, for details see van der Wel et al.
(2021) Appendix A) in the calculations of population median and scatter to
account for the LEGA-C survey targeting strategy.

To demonstrate our approach to combining the posteriors from individual fits,
we show the individual median and population trends in an example subset of
LEGA-C data in Figure 5. The top row shows individual median SFHs from the two
methods for the most massive quiescent galaxies in the stellar mass bin 11.5 <
log($\mathrm{M_{*}/M_{\odot}}$) < 12\. Bagpipes SFHs are shown in orange and
Prospector SFHs are shown in green. The middle and right panels show these
SFHs in specific star-formation rate (log (sSFR)) as a function of time. sSFR
at an epoch is obtained by dividing SFR at that epoch by the total stellar
mass formed up until that epoch excluding mass loss. The center panel includes
median SFHs for individual galaxies and the right panel collapses these to
show only the population distributions. The bottom row shows these trends for
most massive star-forming population. Note that Bagpipes has more variety in
onset, duration and quenching of star-formation for individual galaxies, which
are often imprinted on the population trends even though the individual
parametric SFHs are smooth. In addition to the 126 galaxies with questionable
fits that are discussed in Section §2, there is a significant number of
quiescent galaxies in Bagpipes with $t_{50}$ formation times between 3 Gyr to
4 Gyr, as evident in the SFHs in Figure 5 (similarly for other mass bins, see
Appendix Figure A2).

Figure 6: Median specific star-formation rates (log($\mathrm{sSFR}$)) with
lookback time for populations of star-forming (left), quiescent (center) and
all (right) galaxies sorted in bins of their stellar mass. Colors indicate
stellar mass bins - blue < green < red in increasing mass. Bagpipes/
Prospector results are shown in dotted-lines/solid-lines with hashed/shaded
regions representing 16th-84th percentile population scatter. Bottom most row
shows median trends of all mass bins with thin shaded regions representing
standard error on the medians via bootstrapping. Transition region from star
forming to quiescent is shown in grey band. Star-forming population shows
overall better agreement between two methods than quiescent population. Note
that Bagpipes has more bursty pathways to reach quiescence. Figure 7: The
cumulative stellar mass fractions versus time for star-forming (left),
quiescent (center) and all (right) galaxies sorted by stellar mass. The color
and plotting scheme is the same as in Fig 4. This represents a complete census
as individual galaxies are weighted by correction factors to reflect the
LEGA-C targeting. Lower mass galaxies grow more slowly than their massive
counterparts. Non-parametric SFHs (Prospector, solid) start earlier and are
more extended, especially for quiescent galaxies.

Figure 6 expands the right panel in Figure 5 to show the $\mathrm{sSFR}$
evolution in three stellar mass bins (rows) for: star-forming galaxies (left),
quiescent galaxies (center), and the full population (right). Each panel shows
the population medians (dashed lines/solid lines) and 16th - 84th percentile
population scatter (hashed/filled regions) for Bagpipes and Prospector
respectively. The bottom row combines the median trends of all mass bins with
thin shaded regions showing 1-$\sigma$ error on the medians calculated from
bootstrapping. We adopt the transition and quiescence boundaries as (1/3
$t_{H}$) and (1/20 $t_{H}$) from Tacchella et al. (2022), where $t_{H}$ is the
Hubble time at $z_{median}=0.8$. Although individual galaxies show a variety
of SFHs, the population median trends are largely independent of stellar mass,
with slight trends at early times that are dominated by the demographics of
the LEGA-C sample and modeling priors. In part, differences in the onset of
star formation can be partially due to slight differences in the mass-
dependent redshift distributions of the LEGA-C sample. However, as discussed
in §2, the modeling priors also introduce subtle differences. Bagpipes
consistently measures slightly higher median sSFRs, especially within the
quiescent population and shows higher population scatter compared to
Prospector. Median sSFR trends within the star-forming population show better
agreement between the two codes. To quantify this, for each mass bin, we
calculate the median of the differences between the curves. In Figure 6, from
top to bottom - median $\Delta sSFR(t)|_{SF}$ = 0.08, 0.12, 0.31 $yr^{-1}$
(left column), median $\Delta sSFR(t)|_{Q}$ = 0.22, 0.25, 0.29 $yr^{-1}$
(middle column) and median $\Delta sSFR(t)|_{ALL}$ = 0.11, 0.18, 0.30
$yr^{-1}$ (right column). For star-forming galaxies, both codes infer that
massive galaxies fall to lower sSFRs at the point of observation than the less
massive ones. Some differences are seen in the median trends of the most
massive galaxies (red) from the two codes that could be attributed to small
sample size. For quiescent galaxies, we see a large diversity of Bagpipes sSFR
tracks. We do not find that the population of massive quiescent galaxies shut
off before lower mass counterparts. Note that this is different from local
universe findings (Gallazzi et al., 2006; Nelan et al., 2005).

Figure 7 focuses on the median trends in the cumulative stellar mass growth
and adopts the same plotting conventions as Figure 6, except the horizontal
axis now shows the age of the Universe, starting at the Big Bang. The hashed
and filled regions show the 16th - 84th percentile distributions of individual
median SFHs from Bagpipes and Prospector respectively. Horizontal grey dotted
lines correspond to 10%, 50% and 90% formation thresholds. These timescales,
which we hereby define as $t_{10}$, $t_{50}$, $t_{90}$, are widely used
metrics to quantify the timescales of stellar mass growth in a galaxy
including all progenitors up to the point of observation (Behroozi et al.,
2019; Pacifici et al., 2016a; Weisz et al., 2011, 2008). These values are
quantified for the full primary sample in Table 3 and Table 4. As we find in
the previous two figures, parametric models in general show higher population
scatter than non-parametric ones due to galaxy to galaxy variations. The
bottom row compiles all median trends to compare the rate of stellar mass
growth across different masses. This cumulative view highlights both mass
trends within the populations and discrepancies at early times due to modeling
degeneracies and prior assumptions. At any given time, low mass galaxies have
grown less than massive galaxies. Significant discrepancies in $t_{10}$, and
to some degree in $t_{50}$, suggest that SFHs are prior dominated at large
lookback times. Perhaps similar deep spectroscopic observations of galaxies at
even earlier times (e.g. with JWST) will help to break modeling degeneracies
and resolve this issue, but at this epoch, the measurements of the earliest
SFHs of massive galaxies will remain prior dominated, even with high-S/N
spectroscopic data.

Figure 8: $t_{90}$ (brown) and $t_{50}$ (teal) formation times/redshifts of
star-forming, quiescent and total galaxy populations versus stellar mass.
Error bars/shaded regions indicate the 16th - 84th percentile population
scatter. The bottom row represents the difference in median values of the
Bagpipes and Prospector formation times. The far right panel shows the
vertical axis in equivalent mass-weighted ages. As the redshift of individual
galaxies is different, this is a first order age estimate assuming all
galaxies are observed at median LEGA-C redshift of z $\sim$ 0.8. Late
formation times ($t_{90}$) are consistent between models, with greater
disagreement at earlier times ($t_{50}$), especially for quiescent galaxies.
The strong mass trends within the full population (right panel) are primarily
driven by star forming galaxies that dominate the low-mass end. Figure 9:
Similar to Fig. 8, now binned in stellar velocity dispersion. We find a
striking correlation between the formation timescales and stellar velocity
dispersion in full population (right). This trend is in part driven by
demographics; the quiescent fraction depends more strongly on stellar velocity
dispersion than stellar mass (e.g., Taylor et al., 2022).

Figures 8 and 9 show $t_{50}$ (teal) and $t_{90}$ (brown) the formation time
scales for the populations of star-forming (left), quiescent (center) and all
galaxies (right) versus stellar mass (6) and stellar velocity dispersion (8).
We exclude $t_{10}$ because of the dramatic disagreement between the two sets
of models at early times, suggesting that its measurements are completely
prior dominated. Bagpipes median trends and population scatter are shown with
dotted lines and error bars and Prospector as solid lines with shaded bands.
The bottom row shows the difference between the median values derived from the
two methods with (very small) error bars indicating the standard error in the
medians. The measured late formation times ($t_{90}$) agree well across all
stellar masses, with a slight offset in early formation times ($t_{50}$)
especially of quiescent systems ($t_{50}$ \- $\Delta_{median,SF}$ < 0.9 Gyr,
$\Delta_{median,Q}$ < 1.8 Gyr, $\Delta_{median,ALL}$ < 0.9 Gyr and $t_{90}$ \-
$\Delta_{median,SF}$ < 0.4 Gyr, $\Delta_{median,Q}$ < 0.3 Gyr,
$\Delta_{median,ALL}$ < 0.6 Gyr). This consistent offset in $t_{50}$ for
quiescent systems could be partially driven by preference for more flexible
onset of star-formation in Bagpipes compared to Prospector. We find clear
trends in both $t_{50}$ and $t_{90}$ with stellar mass amongst star-forming
galaxies. In contrast, the median formation times are mostly independent of
mass amongst the quiescent galaxies. The overall correlation with mass is
reflected by the full population, in which the younger star-forming galaxies
dominate at low masses and are rare at the massive end. We note that similar
correlations with stellar mass remain within finer redshift bins ($\Delta
z=0.1$). Correlations of formation times with stellar velocity dispersion
qualitatively similar to those with stellar mass (Figure 9). Again, the median
$t_{50}$ and $t_{90}$ formation times are essentially uniform for the
quiescent population. However, because the fraction of star-forming galaxies
plummets at high stellar velocity dispersion (especially above
$\sigma_{\star}\gtrsim 200$ km s-1) (Taylor et al., 2022), the formation
timescales of the full population exhibits an even stronger correlation with
stellar velocity dispersion (Fig. 9,right panel).

Finally, we investigate the duration of late time star formation ($t_{90}$ \-
$t_{50}$) from Prospector (green) and Bagpipes (orange) in Figure 10. Some
studies have used ($t_{80}$ \- $t_{20}$) and ($t_{90}$ \- $t_{10}$) timescales
(Pacifici et al., 2016a; Tacchella et al., 2022) to quantify full duration of
star-formation, but we choose to avoid large uncertainties at earlier times
and quantify late-time star-formation as $\tau_{SF,late}\equiv t_{90}$ \-
$t_{50}$. Figure 10 shows these timescales versus stellar mass (top row) and
stellar velocity dispersion (bottom row) with error bars/shaded regions
representing the population scatter, following the plotting conventions of
Figures 8 and 9. Here we see that while agreement is reasonably good between
e.g., median $t_{50}$ in the two models, systematic offsets in the duration of
star-formation remain. Non-parametric models from Prospector exhibit
consistently more extended star formation histories than those derived with
double-power law models. This discrepancy at the earliest times has been well-
documented (e.g., Carnall et al., 2019a; Leja et al., 2019). Notably, the
median $\tau_{SF,late}$ from Bagpipes is dramatically ($\sim$2 Gyr) more rapid
for quiescent galaxies, with an enormous population scatter, reflecting the
larger variety of SFHs (see e.g., Fig. 5). Mock recovery testing of
$\tau_{SF,late}$ using spectro-photometric data could reveal if this timescale
is recoverable. Leja et al. (2019) and Carnall et al. (2019b) have shown that
with mock photometric data, Prospector was able to recover late-time star-
formation better than Bagpipes. Note that low $\tau_{SF,late}$ values in
Bagpipes are partially driven by flexibility in onset times as well as in the
late time SFH parameter (e.g. falling slope ‘alpha’) to adapt steep/shallow
values to give a range of $t_{90}$ \- $t_{50}$ values. In fact, the population
scatter on Bagpipes star formation duration detracts from the informative
nature of the median trends. However, significant offset between
$\tau_{SF,late}$ of star-forming and quiescent galaxies does imprint a trend
in the full population median value with stellar velocity dispersion (bottom
row, right-most panel), where demographics are more cleanly separated than in
stellar mass (e.g., Franx et al., 2008; Belli et al., 2014; Taylor et al.,
2022, Fig. 9). The median duration of star-formation derived by Prospector
does not vary significantly within the star-forming population, but quiescent
galaxies exhibit a slight decrease with stellar mass and stellar velocity
dispersion ($\lesssim 0.5$ Gyr across the full range). This is consistent with
weak correlations between alpha enhancement and stellar velocity dispersion
for a subset of quiescent galaxies from the same LEGA-C dataset Beverage et
al. (2023).

##  4 Discussion

Figure 10: The distribution of duration of star-formation at late times
($\tau_{SF,late}\equiv t_{90}$ \- $t_{50}$) in stellar mass (top row) and
stellar velocity dispersion (bottom row) bins for population of star-forming
(left), quiescent (center) and all (right) as recovered from Bagpipes and
Prospector. Green solid lines and orange solid points represent 50th
percentile values (median) and broad green shaded regions and orange error
bars represent 16th-84th percentile population scatter in that bin. The
standard error on the medians is represented in slightly darker thin shaded
regions calculated from bootstrapping. Non-parametric models are consistently
more extended than the double-power law SFHs. We find no strong correlation
between $\tau_{SF,late}$ and stellar mass or stellar velocity dispersion. This
lack of correlation is in contrast with trends found amongst analogous massive
elliptical galaxies today (e.g., Graves et al., 2009; Greene et al., 2022).
The significant offset between $\tau_{SF,late}$ of star-forming and quiescent
galaxies does imprint a trend in the full population median value with stellar
velocity dispersion (bottom row, right-most panel), where demographics are
more cleanly separated than in stellar mass.

Star-formation activity has been demonstrated to depend both on stellar mass
(Calvi et al., 2018; Pacifici et al., 2016a; Kauffmann et al., 2004) and on
stellar velocity dispersion (Franx et al., 2008). It has been long debated
whether stellar velocity dispersion is a more fundamental property than
stellar mass (Sharma et al., 2021; Zahid et al., 2018; Wake et al., 2012), as
it is directly connected to the total gravitational potential well (including
effects of dark matter halo and supermassive black hole at galactic center) in
which a galaxy resides (Elahi et al., 2018; Dutton et al., 2010).

This census of star formation histories builds on extensive studies based on
complete photometric (e.g., Olsen et al., 2021; Aufort et al., 2020; Iyer,
2019; Iyer & Gawiser, 2017; Pacifici et al., 2016b, a; Dye, 2008) and smaller
and/or less complete, and more observationally expensive, spectroscopic
dataset (e.g., Gallazzi et al., 2008; Sánchez-Blázquez et al., 2011; Carnall
et al., 2019b, 2022; Tacchella et al., 2022). For the first time at
significant lookback time, we present a comprehensive view of the star-
formation histories of the full population of massive galaxies (not just
quiescent or star forming galaxies separately). This builds on a preliminary
study by Chauke et al. (2018) based on early LEGA-C data (607 galaxies) with
simpler modeling assumptions: adopting fixed solar metallicity and piece-wise
constant star-formation rates using CSP templates from FSPS package (Conroy et
al., 2009). This current paper expands in both scope and detail to yield a
comprehensive analysis of the full LEGA-C sample.

Figure 11: Comparison of median formation times $t_{50}$ (teal) and $t_{90}$
(brown) measured forward in the direction of the age of the universe in bins
of dynamical masses. These are derived for the progenitor populations from the
LEGA-C sample and compared to those inferred from local elliptical galaxies
from ATLAS3D (after correcting for stellar mass growth between $z\sim 0.8$ and
$z\sim 0$) in McDermid et al. (2015). The left panel shows the comparison with
the full LEGA-C population (star-forming and quiescent progenitors) and the
right panel compares only quiescent galaxies to the local ellipticals. Error
bars and shaded regions represent 1-$\sigma$ dispersion in the population
medians. SPS modeling of local ellipticals finds ubiquitously earlier
($t_{50}$ $\sim$ 2 Gyr) star formation (diamonds), which is even more extreme
than even non-parametric SFHs for LEGA-C progenitors. The strong trends in
later star formation ($t_{90}$) suggest the importance of star forming
progenitors joining the population at $M_{dyn}\sim 10^{11}M_{\odot}$ since
$z\sim 0.8$.

Our results are in qualitative agreement with previous analyses, and
quantitative differences can generally be attributed to modeling systematics.
Previous studies have found similar correlations between formation times
$t_{50}$ and stellar mass and stellar velocity dispersion e.g. for 10.5 <
log($\mathrm{M_{*}/M_{\odot}}$) < 12 linear relation goes from (forward in
time) 1.8 Gyr to 0.9 Gyr for quiescent sample in Tacchella et al. (2022), 1.9
Gyr to 1.4 Gyr for quiescent and 3 Gyr to 2.5 Gyr for star-forming in
(Ferreras et al., 2019), $\sim$ 1 Gyr for quiescent sample and $\sim$ 4 Gyr to
2 Gyr for star forming sample across full mass range in (Chauke et al., 2018),
which is similar to what we find in our analysis in Fig 8 and Fig 9.
Interestingly, both of our modeling methods recover significant population
scatter in the late time formation timescales of star forming galaxies,
however the median $\tau_{SF,late}$ is almost independent of stellar mass.
This result is in contrast with previous studies. One interesting comparison
is with Estrada-Carpenter et al. (2020), in which $\tau_{SF,late}$ is defined
as quenching timescale ($t_{Q}$). Using Prospector to fit non-parametric SFHs
to a much smaller sample of $\sim$30 quiescent systems at 0.7 < $z$ < 1,
finding median formation redshift $z50$ roughly corresponding to 3.3 Gyr - 1.5
Gyr formation times for 10.5 < log($\mathrm{M_{*}/M_{\odot}}$) < 12 which is
consistent within our measured population scatter (Fig.8: Bagpipes \- 4.2 Gyr
- 3.0 Gyr, Prospector \- 2.8 Gyr - 2.3 Gyr). However, they find a broader
range of quenching timescales varying from 0.4 Gyr to 2.2 Gyr which is
different from our $\tau_{SF,late}$ derived using Prospector (Fig.10: 1.7 Gyr
to 2.7 Gyr) but similar to what we find in Bagpipes (Fig.10: 0.2 Gyr to 2.3
Gyr). We can only speculate that this subtle difference could be driven by
priors, sample size or differences in data characteristics (e.g., the use of
low-resolution grism spectra). Beverage et al. (2023) study 135 LEGA-C massive
quiescent galaxies based on elemental abundance patterns and found a slightly
stronger correlation between formation times and stellar velocity dispersion.
The highest dispersion galaxies (> 250 km/s ) formed the earliest around 3.1
Gyr and are most metal rich whereas low dispersion galaxies (< 150 Km/s)
formed around 4.4 Gyr. This is qualitatively similar to our findings in Fig.9,
middle panel.

Interestingly, although the qualitative trends agree well with cluster
galaxies at a similar epoch, there may be some hint of environmental effects
on SFHs. Khullar et al. (2022) find $t_{50}$ of massive quiescent cluster
galaxies to be inversely correlated with stellar mass at 0.3 < $z$ < 1.4 with
values within 2 Gyr. This is $\sim$ 1 Gyr/2 Gyr earlier than our
Prospector/Bagpipes estimates. Again, this difference could be driven by
modeling systematics (e.g. that study adopts a delayed tau SFH) or a physical
manifestation of environmental processes that combine to quench galaxies
earlier or the generally speed up galaxy evolution in clusters. Even within
the LEGA-C dataset age sensitive indices (Dn4000 and H$\delta$) of quiescent
galaxies depend on the environment, not just in cluster versus field; galaxies
in over-dense regions are older and formed earlier than those in less extreme
environments (Sobral et al., 2022). We note that direct comparisons between
different studies can be challenging; Webb et al. (2020) perform a cluster
analysis at similar redshifts and find little difference between the ages of
galaxies in the field versus cluster environments at fixed mass.

Next, we investigate whether our results from LEGA-C at $z\sim 0.8$ differ
significantly from our understanding of this time period as recovered from
previous modeling of their presumed descendent galaxies in the local Universe.
We start by comparing to SPS modeling of local elliptical galaxies from
McDermid et al. (2015). Figure 11 shows a comparison of the median formation
times $t_{50}$ and $t_{90}$ as a function of the Virial mass for all galaxies
in our sample at < $z$ > = 0.8 from Bagpipes (dashed line) and Prospector
(solid-line) relative to quiescent systems at z $\sim$ 0 (diamonds) from
McDermid et al. (2015), after subtracting the minimal cumulative stellar mass
growth between $z\sim 0$ and $z\sim 0.6$. The error bars and shaded regions
represent 1-$\sigma$ dispersion in the population medians. The vertical axis
is forward in time with t=0 corresponding to the Big Bang. Note that these
elliptical galaxies at $z\sim 0$ could either be quiescent or star-forming at
our redshift, hence for a first order comparison with LEGA-C, we compare to
the full (left) and quiescent only (right) population. Virial masses have been
calibrated for the majority of the LEGA-C sample (van der Wel et al., 2022)
using spatially resolved stellar kinematics (van Houdt et al., 2021). For this
simplistic calculation we adopt a median offset of 0.33 dex between the
stellar and dynamical mass for each LEGA-C individual object. McDermid et al.
(2015) find that most massive galaxies (log$\mathrm{(M_{*}/M_{\odot})}$ > 11)
have already accumulated > 95% of their stellar mass by $z$ $\sim$ 0.6,
whereas low mass galaxies have formed $\sim$90% of their stellar mass. They
also find ubiquitously earlier ($t_{50}=2Gyr$) stellar mass growth (diamonds
with dotted line), than the view from LEGA-C using either modeling framework.
This reveals a natural uncertainty of stellar population modeling for the
oldest local stellar systems. We see similar dynamical mass trends in the left
panel for z$\sim$0.8 galaxies (both star-forming and quiescent) and local
ellipticals especially toward the low mass end. This correlation is weaker
when only quiescent galaxies from the LEGA-C sample are included. This
emphasizes the importance of star-forming progenitors joining the local
elliptical population.

Thomas et al. (2005, 2010) provide clear evidence from a $z$ $\sim$ 0 sample
of galaxies that the duration and timing of galaxy formation depends strongly
on the stellar mass, arguing that this ultimately reflects halo assembly.
However, even at $z\sim 0.8$, we do not find these trends, even when star-
forming progenitors are included. Some of these differences could be
attributed to modeling systematics e.g. use of Gaussian SFH prescription to
describe population average, chemical enrichment models/abundance patterns,
various environmental dependence and most importantly aperture effects -
probing the central few kpc (not the full galaxy) with a 3” diameter fiber
that makes results more mass dependent for older ages and higher
metallicities, whereas slit aperture effects are not that extreme.
Furthermore, the view from maximally old galaxies in local Universe will
always be complicated by the relative challenge of modeling slowly evolving
old stellar populations. However, it is hard to escape the conclusion that
suggests that the trends found in local ellipticals reflect the transformation
and assembly/merger history of massive galaxies and not just their in-situ
SFHs.

Figure 12: Formation times versus stellar mass for massive quiescent galaxies
across redshifts. The grey lines and bands indicate relations at lower and
higher redshifts respectively (Choi et al., 2014; Carnall et al., 2019b). The
highest redshift measurement from JWST is shown with black star(Carnall et
al., 2023a). Our results are shown in orange (parametric SFHs) and green (non-
parametric SFHs). Shaded color bands represent 16th - 84th percentile
population scatter in the median values. Scaling relations from cosmological
simulations (dashed/dotted lines, Davé et al., 2019; Nelson et al., 2019)
agree well with our Prospector models. The evolution in this relation for
quiescent galaxies emphasizes the importance of including star-forming
progenitors and ex-situ contributions when comparing the stellar populations
of galaxies across cosmic time.

Finally, we compare our measured formation time trends for quiescent galaxies
with similar studies at other epochs. Figure 12 shows the comparison of median
formation times ($t_{50}$) of massive quiescent systems from our study (orange
and green bands indicating populations scatter) with other empirical
(spectroscopic) studies and simulations across redshifts. Choi et al. (2014)
results (grey solid lines) are derived from full-spectral fitting based on
alpha abundance measurements and SSP-equivalent ages. The grey band indicates
mean formation times recovered from the VANDELS survey (1 < $z$ < 1.3) in
Carnall et al. (2019b) using Bagpipes. The single JWST spectroscopic
observation around $z$ = 4.65 from Carnall et al. (2023b) is shown with a
star-symbol. We note that the incredibly small uncertainty on formation time
reflects the incredible data quality and limits imposed by the age of the
Universe at that early time. Trends found in cosmological simulations at $z$ =
1 are shown with dashed (SIMBA (Davé et al., 2019)) and dotted
(IllustrisTNG-100 (Nelson et al., 2019)) lines. For details on the simulation
selection, see Carnall et al. (2019b). Although modeling differences could
introduce significant systematic offsets with respect to our fits (e.g.
extended SFH versus SSPs, scaled solar metallicity versus alpha enrichment),
it is interesting to speculate whether some of the observed offset in $t_{50}$
between epochs is physical, perhaps driven by late additions to the quiescent
population. This progenitor bias could exaggerate the evolution of these
scaling relations without changing the stellar populations of any individual
galaxies. Superficially the older LEGA-C ages from Prospector are roughly
consistent with the trend found for quiescent galaxies in VANDELS at higher z,
however our Bagpipes modeling adopts much more similar priors and provide the
more consistent comparison. Thus, the observed offset between the grey and
orange bands suggests that the relation evolves even between $z\sim 1.15$ and
$z\sim 0.8$. Choi et al. (2014) modeled stacked optical spectra of massive
quiescent galaxies from 0.1 < $z$ < 0.7 from SDSS and AGES Survey. Although
our spectro-photometric modeling of individual galaxies differs from their
analysis, the similar offset to later formation times in lower-redshift
quiescent populations is consistent with new additions shifting the scaling
relations, perhaps more efficiently at lower mass. There is one caveat in
comparing with local elliptical galaxies: galaxy mergers are also major
avenues of mass growth and therefore will also contribute to the spatially-
integrated ages of galaxies (e.g. $t_{50}$ and $t_{90}$). With SPS modeling,
we measure the formation times of stars from all progenitors of a galaxy at
$t_{obs}$, which can be quite different from assembly times ($t_{a}$). Hill et
al. (2017b) demonstrate that this can be a significant effect, by $z$$\sim$0.1
assembly times and light-weighted formation times can differ by up to $\sim
2.5$ Gyr (see also e.g., Hill et al., 2017a; Muzzin, 2017). We hence note that
progenitor-descendant linking is further complicated by uncertain importance
of merging.

Note that one of the key findings of this study is the impact of priors of SPS
modeling on early formation time recovery at z $\sim$ 0.8, even with high S/N
data. A part of this is also driven by the outshining of old populations. A
natural conclusion, therefore, is that the earliest SFHs can be better
understood by modeling the light from galaxies at progressively earlier times,
when this becomes less of a concern. Future spectroscopic studies of higher
redshifts can precisely pinpoint the emergence of massive quiescent systems
and can help settle uncertainties at the earliest times. Recent results from
JWST photometry-only SED analyses suggest that massive galaxies may indeed
form efficiently and rapidly within 1 Gyr of the Big Bang (e.g. Labbé et al.
(2023)). However, given uncertainties in the emission line contributions to
photometry (e.g. McKinney et al. (2023)), spectroscopic data are needed more
than ever. JWST/NIRSpec has been shown to provide the sensitivity required to
constrain the SFHs and and quenching timescales of galaxies at higher
redshifts and lower masses: e.g., the Carnall et al. (2023b)’s extreme massive
quiescent galaxy at $z$ = 4.65 included in Figure 12. Other recent studies
include Looser et al. (2023) and Nanayakkara et al. (2023). Larger, more
complete spectroscopic samples from JWST will be crucial to map the early
evolution of today’s massive quiescent galaxy population. These samples will
hopefully shed more light on the modeling degeneracies at the earliest times.

##  5 Conclusions

In this paper, we analyze the median SFHs of massive
(log($\mathrm{M_{*}/M_{\odot}}$) > 10.5) galaxies at $0.6<z<1.0$ from the
LEGA-C survey by quantifying their formation times ($t_{50}$ and $t_{90}$) and
investigate the population trends in this census with stellar masses and
stellar velocity dispersions. From our spectro-photometric analysis, we
conclude that:

  1. 1.

The two modeling methods (parametric and non-parametric) yield consistent
late-time star-formation histories for star-forming and quiescent galaxies.
Non-parametric SFHs consistently prefer earlier stellar mass formation,
especially for quiescent systems. Potential reasons could be - stellar
population outshining driving this inference into a systematics-limited regime
emphasizing the dominance of modeling priors at the earliest times or a host
of different modeling assumption e.g. SPS models, treatment of emission lines,
dust, parametrization of SFHs, etc. (Figure 2,5).

  2. 2.

Although individual galaxies show a variety of SFH pathways, the median sSFR
evolution is similar for all mass bins. Bagpipes SFHs shows population scatter
at any time. Both codes infer that massive star-forming galaxies exhibit
falling sSFRs with redshift. Neither model finds a mass-dependence in the
median time at which quiescent galaxies cross a sSFR threshold (Figure 6).

  3. 3.

From median trends, we find that for both non-parametric and parametric SFH
modeling, quiescent galaxies formed 90% of their stars well within $\sim$4 Gyr
and $\sim$5.5 Gyr age of the universe for massive and less massive end
respectively. Our analysis reflects the established impact of these priors at
earlier times, yielding differences in mass-weighted ages $t_{50}$ $\sim$1.5
Gyr) within that population, with SFHs that further diverge at earlier times.
(Figures 7,8).

  4. 4.

Lower mass galaxies are slightly younger than massive galaxies (Figure 7),
however this trend is relatively weak. This is perhaps due to the fact that
the full LEGA-C sample only probes a limited dynamic range in stellar mass.
This picture is also consistent with simulations (Iyer et al., 2020; Torrey et
al., 2018; Shamshiri et al., 2015). The formation times of star-forming
galaxies correlate with stellar mass and stellar velocity dispersion, but
similar trends do not exist not amongst the quiescent galaxies. For combined
population, the formation times depend strongly on the $M_{*}$ and
$\sigma_{*}$ as quiescent fraction increases with each property. (Figure 8,9).

  5. 5.

For the star-forming population, both codes recover consistent median
formation times ($|t_{Bagpipes}-t_{Prospector}|<300$ Myr); the most massive
systems formed within $t_{50}$/$t_{90}$ $\sim$3 Gyr/5.5 Gyr and the less
massive ones within $\sim$5 Gyr/7 Gyr (Figure 8, or similarly for high/low
$\sigma_{\star}$, Fig. 9). Even for these galaxies, SFRs were on the decline
by the time of observation, corroborating a number of previous studies (e.g.,
Feulner et al., 2005; Ilbert et al., 2015).

  6. 6.

The late time duration of star formation ($\tau_{SF,late}$ = $t_{90}-t_{50}$)
does not exhibit a significant correlation with $M_{\star}$ or
$\sigma_{\star}$ for quiescent or star-forming galaxy populations. At face
value, this is in contrast with expectations from the local universe (e.g.
Thomas et al. (2005)). We posit that this suggests the importance of
transformation, quenching and/or ex-situ evolution even at late times ($z$ <
1) (Figure 10).

Remaining open questions include connecting quenching timescales (defined from
recovered SFHs), environment, and processes driving quenching in quiescent
systems to stellar mass and stellar velocity dispersion. A number of studies
have begun to answer those questions using photometric data (Mao et al., 2022)
and age sensitive spectral indices (Sobral et al., 2022; Wu et al., 2018b),
but not much effort has been put towards understanding this evolution from the
SFH of individual galaxies. Another interesting comparison would be to compare
our SFHs to those derived from simulation using PCA analysis (e.g., Sparre et
al., 2015).

This SFH census of a representative sample of $\sim 3000$ massive galaxies
using two state-of-the-art Bayesian modeling tools will set the benchmark at
$z\sim 1$. It is a "cosmic mid-point" in time, providing a stepping stone to
future spectroscopic studies of high redshift galaxies (e.g., from
JWST/NIRSpec, VLT/MOONS (Maiolino et al., 2020) or Subaru/PFS (Greene et al.,
2022)). By including two independent sets of models, we hope that this
intermediate redshift study would be a link to other similar studies and help
connect our understanding of evolution of stellar mass growth and assembly in
galaxies from cosmic noon and beyond to present day universe.

##  6 Acknowledgements

RB and YK gratefully acknowledge funding for project KA2019-105551 provided by
the Robert C. Smith Fund and the Betsy R. Clark Fund of The Pittsburgh
Foundation. RB acknowledges support from the Research Corporation for
Scientific Advancement (RCSA) Cottrell Scholar Award ID No: 27587. RB and YK
also acknowledge support from NSF grant AST-2144314. AG acknowledges support
from INAF-Minigrant-2022 "LEGA-C" 1.05.12.04.01. This research was supported
in part by the University of Pittsburgh Center for Research Computing,
RRID:SCR_022735, through the resources provided. Specifically, this work used
the H2P cluster, which is supported by NSF award number OAC-2117681. We
specifically acknowledge the assistance of Prof. Kim Wong to setup the
framework at the cluster. YK thanks Brett Andrews, Alan Pearl, David Setton
and Biprateep Dey for their constructive feedback and suggestions that greatly
improved this work.

This work made extensive use of the publicly available tools Bagpipes (Carnall
et al., 2018), Prospector (Johnson et al., 2021; Johnson & Leja, 2017), Python
programming language (van Rossum, 1995), astropy (Astropy Collaboration et
al., 2013), matplotlib (Hunter, 2007), numpy (Harris et al., 2020), scipy
(Jones et al., 2001).

## References

  * Abramson et al. (2015) Abramson, L. E., Gladders, M. D., Dressler, A., et al. 2015, ApJ, 801, L12, doi: 10.1088/2041-8205/801/1/L12
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Aufort et al. (2020) Aufort, G., Ciesla, L., Pudlo, P., & Buat, V. 2020, A&A, 635, A136, doi: 10.1051/0004-6361/201936788
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Bell et al. (2004) Bell, E. F., Wolf, C., Meisenheimer, K., et al. 2004, ApJ, 608, 752, doi: 10.1086/420778
  * Belli et al. (2018) Belli, S., Contursi, A., & Davies, R. I. 2018, MNRAS, 478, 2097, doi: 10.1093/mnras/sty1236
  * Belli et al. (2014) Belli, S., Newman, A. B., & Ellis, R. S. 2014, ApJ, 783, 117, doi: 10.1088/0004-637X/783/2/117
  * Beverage et al. (2023) Beverage, A. G., Kriek, M., Conroy, C., et al. 2023, ApJ, 948, 140, doi: 10.3847/1538-4357/acc176
  * Bezanson et al. (2018) Bezanson, R., van der Wel, A., Straatman, C., et al. 2018, ApJ, 868, L36, doi: 10.3847/2041-8213/aaf16b
  * Bowman et al. (2020) Bowman, W. P., Zeimann, G. R., Nagaraj, G., et al. 2020, ApJ, 899, 7, doi: 10.3847/1538-4357/ab9f3c
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x
  * Calvi et al. (2018) Calvi, R., Vulcani, B., Poggianti, B. M., et al. 2018, MNRAS, 481, 3456, doi: 10.1093/mnras/sty2476
  * Cappellari (2017) Cappellari, M. 2017, MNRAS, 466, 798, doi: 10.1093/mnras/stw3020
  * Carnall et al. (2019a) Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019a, ApJ, 873, 44, doi: 10.3847/1538-4357/ab04a2
  * Carnall et al. (2018) Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379, doi: 10.1093/mnras/sty2169
  * Carnall et al. (2019b) Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2019b, MNRAS, 490, 417, doi: 10.1093/mnras/stz2544
  * Carnall et al. (2022) —. 2022, ApJ, 929, 131, doi: 10.3847/1538-4357/ac5b62
  * Carnall et al. (2023a) —. 2023a, Nature, 619, 716, doi: 10.1038/s41586-023-06158-6
  * Carnall et al. (2023b) Carnall, A. C., McLeod, D. J., McLure, R. J., et al. 2023b, MNRAS, 520, 3974, doi: 10.1093/mnras/stad369
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Charlot & Fall (2000) Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718, doi: 10.1086/309250
  * Chauke et al. (2018) Chauke, P., van der Wel, A., Pacifici, C., et al. 2018, ApJ, 861, 13, doi: 10.3847/1538-4357/aac324
  * Chaves-Montero & Hearin (2020) Chaves-Montero, J., & Hearin, A. 2020, MNRAS, 495, 2088, doi: 10.1093/mnras/staa1230
  * Chevallard & Charlot (2016) Chevallard, J., & Charlot, S. 2016, MNRAS, 462, 1415, doi: 10.1093/mnras/stw1756
  * Chevallard et al. (2019) Chevallard, J., Curtis-Lake, E., Charlot, S., et al. 2019, MNRAS, 483, 2621, doi: 10.1093/mnras/sty2426
  * Choi et al. (2014) Choi, J., Conroy, C., Moustakas, J., et al. 2014, ApJ, 792, 95, doi: 10.1088/0004-637X/792/2/95
  * Cid Fernandes (2007) Cid Fernandes, R. 2007, in Stellar Populations as Building Blocks of Galaxies, ed. A. Vazdekis & R. Peletier, Vol. 241, 461–469, doi: 10.1017/S1743921307008794
  * Cid Fernandes et al. (2005) Cid Fernandes, R., Mateus, A., Sodré, L., Stasińska, G., & Gomes, J. M. 2005, MNRAS, 358, 363, doi: 10.1111/j.1365-2966.2005.08752.x
  * Ciesla et al. (2017) Ciesla, L., Elbaz, D., & Fensch, J. 2017, A&A, 608, A41, doi: 10.1051/0004-6361/201731036
  * Citro et al. (2016) Citro, A., Pozzetti, L., Moresco, M., & Cimatti, A. 2016, A&A, 592, A19, doi: 10.1051/0004-6361/201527772
  * Cohn (2018) Cohn, J. D. 2018, MNRAS, 478, 2291, doi: 10.1093/mnras/sty1148
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486
  * Cullen et al. (2019) Cullen, F., McLure, R. J., Dunlop, J. S., et al. 2019, MNRAS, 487, 2038, doi: 10.1093/mnras/stz1402
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827, doi: 10.1093/mnras/stz937
  * Diemer et al. (2017) Diemer, B., Sparre, M., Abramson, L. E., & Torrey, P. 2017, ApJ, 839, 26, doi: 10.3847/1538-4357/aa68e5
  * Draine & Li (2007) Draine, B. T., & Li, A. 2007, ApJ, 657, 810, doi: 10.1086/511055
  * Dutton et al. (2010) Dutton, A. A., Conroy, C., van den Bosch, F. C., Prada, F., & More, S. 2010, MNRAS, 407, 2, doi: 10.1111/j.1365-2966.2010.16911.x
  * Dye (2008) Dye, S. 2008, MNRAS, 389, 1293, doi: 10.1111/j.1365-2966.2008.13639.x
  * Elahi et al. (2018) Elahi, P. J., Power, C., Lagos, C. d. P., Poulton, R., & Robotham, A. S. G. 2018, MNRAS, 477, 616, doi: 10.1093/mnras/sty590
  * Estrada-Carpenter et al. (2020) Estrada-Carpenter, V., Papovich, C., Momcheva, I., et al. 2020, ApJ, 898, 171, doi: 10.3847/1538-4357/aba004
  * Falcón-Barroso et al. (2011) Falcón-Barroso, J., Sánchez-Blázquez, P., Vazdekis, A., et al. 2011, A&A, 532, A95, doi: 10.1051/0004-6361/201116842
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mexicana Astron. Astrofis., 49, 137.  https://arxiv.org/abs/1302.4485
  * Feroz & Hobson (2008) Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449, doi: 10.1111/j.1365-2966.2007.12353.x
  * Feroz et al. (2019) Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, The Open Journal of Astrophysics, 2, 10, doi: 10.21105/astro.1306.2144
  * Feroz & Skilling (2013) Feroz, F., & Skilling, J. 2013, in American Institute of Physics Conference Series, Vol. 1553, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 32nd International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. U. von Toussaint, 106–113, doi: 10.1063/1.4819989
  * Ferreras et al. (2019) Ferreras, I., Pasquali, A., Pirzkal, N., et al. 2019, MNRAS, 486, 1358, doi: 10.1093/mnras/stz849
  * Feulner et al. (2005) Feulner, G., Goranova, Y., Drory, N., Hopp, U., & Bender, R. 2005, MNRAS, 358, L1, doi: 10.1111/j.1745-3933.2005.00012.x
  * Forrest et al. (2020) Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020, ApJ, 903, 47, doi: 10.3847/1538-4357/abb819
  * Franx et al. (2008) Franx, M., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2008, ApJ, 688, 770, doi: 10.1086/592431
  * Gallazzi et al. (2008) Gallazzi, A., Brinchmann, J., Charlot, S., & White, S. D. M. 2008, MNRAS, 383, 1439, doi: 10.1111/j.1365-2966.2007.12632.x
  * Gallazzi et al. (2006) Gallazzi, A., Charlot, S., Brinchmann, J., & White, S. D. M. 2006, MNRAS, 370, 1106, doi: 10.1111/j.1365-2966.2006.10548.x
  * Gladders et al. (2013) Gladders, M. D., Oemler, A., Dressler, A., et al. 2013, ApJ, 770, 64, doi: 10.1088/0004-637X/770/1/64
  * Graves et al. (2009) Graves, G. J., Faber, S. M., & Schiavon, R. P. 2009, ApJ, 698, 1590, doi: 10.1088/0004-637X/698/2/1590
  * Greene et al. (2022) Greene, J., Bezanson, R., Ouchi, M., Silverman, J., & the PFS Galaxy Evolution Working Group. 2022, arXiv e-prints, arXiv:2206.14908, doi: 10.48550/arXiv.2206.14908
  * Hamadouche et al. (2023) Hamadouche, M. L., Carnall, A. C., McLure, R. J., et al. 2023, MNRAS, 521, 5400, doi: 10.1093/mnras/stad773
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hill et al. (2017a) Hill, A. R., Muzzin, A., Franx, M., & Marchesini, D. 2017a, ApJ, 849, L26, doi: 10.3847/2041-8213/aa951a
  * Hill et al. (2017b) Hill, A. R., Muzzin, A., Franx, M., et al. 2017b, ApJ, 837, 147, doi: 10.3847/1538-4357/aa61fe
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ibarra-Medel et al. (2016) Ibarra-Medel, H. J., Sánchez, S. F., Avila-Reese, V., et al. 2016, MNRAS, 463, 2799, doi: 10.1093/mnras/stw2126
  * Ilbert et al. (2015) Ilbert, O., Arnouts, S., Le Floc’h, E., et al. 2015, A&A, 579, A2, doi: 10.1051/0004-6361/201425176
  * Iyer (2019) Iyer, K. 2019, in American Astronomical Society Meeting Abstracts, Vol. 233, American Astronomical Society Meeting Abstracts #233, 429.05
  * Iyer & Gawiser (2017) Iyer, K., & Gawiser, E. 2017, ApJ, 838, 127, doi: 10.3847/1538-4357/aa63f0
  * Iyer et al. (2020) Iyer, K. G., Tacchella, S., Genel, S., et al. 2020, MNRAS, 498, 430, doi: 10.1093/mnras/staa2150
  * Johnson & Leja (2017) Johnson, B., & Leja, J. 2017, Bd-J/Prospector: Initial Release, v0.1, Zenodo, Zenodo, doi: 10.5281/zenodo.1116491
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22, doi: 10.3847/1538-4365/abef67
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python.  http://www.scipy.org/
  * Kauffmann et al. (2004) Kauffmann, G., White, S. D. M., Heckman, T. M., et al. 2004, MNRAS, 353, 713, doi: 10.1111/j.1365-2966.2004.08117.x
  * Khullar et al. (2022) Khullar, G., Bayliss, M. B., Gladders, M. D., et al. 2022, ApJ, 934, 177, doi: 10.3847/1538-4357/ac7c0c
  * Kriek & Conroy (2013) Kriek, M., & Conroy, C. 2013, ApJ, 775, L16, doi: 10.1088/2041-8205/775/1/L16
  * Kriek et al. (2015) Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15, doi: 10.1088/0067-0049/218/2/15
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Labbé et al. (2023) Labbé, I., van Dokkum, P., Nelson, E., et al. 2023, Nature, 616, 266, doi: 10.1038/s41586-023-05786-2
  * Le Fevre et al. (2000) Le Fevre, O., Saisse, M., Mancini, D., et al. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4008, Optical and IR Telescope Instrumentation and Detectors, ed. M. Iye & A. F. Moorwood, 546–557, doi: 10.1117/12.395513
  * Leitner (2012) Leitner, S. N. 2012, ApJ, 745, 149, doi: 10.1088/0004-637X/745/2/149
  * Leja et al. (2019) Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3, doi: 10.3847/1538-4357/ab133c
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170, doi: 10.3847/1538-4357/aa5ffe
  * Looser et al. (2023) Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2023, arXiv e-prints, arXiv:2302.14155, doi: 10.48550/arXiv.2302.14155
  * Lower et al. (2020) Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33, doi: 10.3847/1538-4357/abbfa7
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Maiolino et al. (2020) Maiolino, R., Cirasuolo, M., Afonso, J., et al. 2020, The Messenger, 180, 24, doi: 10.18727/0722-6691/5197
  * Maltby et al. (2018) Maltby, D. T., Almaini, O., Wild, V., et al. 2018, MNRAS, 480, 381, doi: 10.1093/mnras/sty1794
  * Mao et al. (2022) Mao, Z., Kodama, T., Pérez-Martínez, J. M., et al. 2022, A&A, 666, A141, doi: 10.1051/0004-6361/202243733
  * Martins (2021) Martins, L. P. 2021, in Galaxy Evolution and Feedback across Different Environments, ed. T. Storchi Bergmann, W. Forman, R. Overzier, & R. Riffel, Vol. 359, 386–390, doi: 10.1017/S1743921320001647
  * McDermid et al. (2015) McDermid, R. M., Alatalo, K., Blitz, L., et al. 2015, MNRAS, 448, 3484, doi: 10.1093/mnras/stv105
  * McKinney et al. (2023) McKinney, J., Finnerty, L., Casey, C. M., et al. 2023, ApJ, 946, L39, doi: 10.3847/2041-8213/acc322
  * McLure et al. (2018) McLure, R. J., Pentericci, L., Cimatti, A., et al. 2018, MNRAS, 479, 25, doi: 10.1093/mnras/sty1213
  * Mortlock et al. (2017) Mortlock, A., McLure, R. J., Bowler, R. A. A., et al. 2017, MNRAS, 465, 672, doi: 10.1093/mnras/stw2728
  * Muzzin (2017) Muzzin, A. 2017, in Early stages of Galaxy Cluster Formation, 24, doi: 10.5281/zenodo.833363
  * Muzzin et al. (2013a) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013a, ApJS, 206, 8, doi: 10.1088/0067-0049/206/1/8
  * Muzzin et al. (2013b) —. 2013b, ApJS, 206, 8, doi: 10.1088/0067-0049/206/1/8
  * Muzzin et al. (2013c) —. 2013c, ApJ, 777, 18, doi: 10.1088/0004-637X/777/1/18
  * Nanayakkara et al. (2023) Nanayakkara, T., Glazebrook, K., Jacobs, C., et al. 2023, ApJ, 947, L26, doi: 10.3847/2041-8213/acbfb9
  * Nelan et al. (2005) Nelan, J. E., Smith, R. J., Hudson, M. J., et al. 2005, ApJ, 632, 137, doi: 10.1086/431962
  * Nelson et al. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x
  * Newman et al. (2013) Newman, J. A., Cooper, M. C., Davis, M., et al. 2013, ApJS, 208, 5, doi: 10.1088/0067-0049/208/1/5
  * Ocvirk et al. (2006) Ocvirk, P., Pichon, C., Lançon, A., & Thiébaut, E. 2006, MNRAS, 365, 46, doi: 10.1111/j.1365-2966.2005.09182.x
  * Olsen et al. (2021) Olsen, C., Gawiser, E., Iyer, K., et al. 2021, ApJ, 913, 45, doi: 10.3847/1538-4357/abf3c2
  * Pacifici et al. (2012) Pacifici, C., Charlot, S., Blaizot, J., & Brinchmann, J. 2012, MNRAS, 421, 2002, doi: 10.1111/j.1365-2966.2012.20431.x
  * Pacifici et al. (2016a) Pacifici, C., Oh, S., Oh, K., Lee, J., & Yi, S. K. 2016a, ApJ, 824, 45, doi: 10.3847/0004-637X/824/1/45
  * Pacifici et al. (2016b) Pacifici, C., Kassin, S. A., Weiner, B. J., et al. 2016b, ApJ, 832, 79, doi: 10.3847/0004-637X/832/1/79
  * Pacifici et al. (2023) Pacifici, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, 141, doi: 10.3847/1538-4357/acacff
  * Panter et al. (2007) Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2007, MNRAS, 378, 1550, doi: 10.1111/j.1365-2966.2007.11909.x
  * Pforr et al. (2012) Pforr, J., Maraston, C., & Tonini, C. 2012, MNRAS, 422, 3285, doi: 10.1111/j.1365-2966.2012.20848.x
  * Rudie et al. (2012) Rudie, G. C., Steidel, C. C., Trainor, R. F., et al. 2012, ApJ, 750, 67, doi: 10.1088/0004-637X/750/1/67
  * Sánchez-Blázquez et al. (2011) Sánchez-Blázquez, P., Ocvirk, P., Gibson, B. K., Pérez, I., & Peletier, R. F. 2011, MNRAS, 415, 709, doi: 10.1111/j.1365-2966.2011.18749.x
  * Schawinski et al. (2014) Schawinski, K., Urry, C. M., Simmons, B. D., et al. 2014, MNRAS, 440, 889, doi: 10.1093/mnras/stu327
  * Schiavon et al. (2006) Schiavon, R. P., Faber, S. M., Konidaris, N., et al. 2006, ApJ, 651, L93, doi: 10.1086/509074
  * Schreiber et al. (2016) Schreiber, C., Elbaz, D., Pannella, M., et al. 2016, A&A, 589, A35, doi: 10.1051/0004-6361/201527200
  * Shamshiri et al. (2015) Shamshiri, S., Thomas, P. A., Henriques, B. M., et al. 2015, MNRAS, 451, 2681, doi: 10.1093/mnras/stv883
  * Sharma et al. (2021) Sharma, S., Hayden, M. R., Bland-Hawthorn, J., et al. 2021, MNRAS, 506, 1761, doi: 10.1093/mnras/stab1086
  * Simha et al. (2014) Simha, V., Weinberg, D. H., Conroy, C., et al. 2014, arXiv e-prints, arXiv:1404.0402.  https://arxiv.org/abs/1404.0402
  * Siudek et al. (2017) Siudek, M., Małek, K., Scodeggio, M., et al. 2017, A&A, 597, A107, doi: 10.1051/0004-6361/201628951
  * Skilling (2006) Skilling, J. 2006, Bayesian Analysis, 1, 833 , doi: 10.1214/06-BA127
  * Smethurst et al. (2015) Smethurst, R. J., Lintott, C. J., Simmons, B. D., et al. 2015, MNRAS, 450, 435, doi: 10.1093/mnras/stv161
  * Sobral et al. (2022) Sobral, D., van der Wel, A., Bezanson, R., et al. 2022, ApJ, 926, 117, doi: 10.3847/1538-4357/ac4419
  * Sparre et al. (2015) Sparre, M., Hayward, C. C., Springel, V., et al. 2015, MNRAS, 447, 3548, doi: 10.1093/mnras/stu2713
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Steidel et al. (2014) Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165, doi: 10.1088/0004-637X/795/2/165
  * Straatman et al. (2018) Straatman, C. M. S., van der Wel, A., Bezanson, R., et al. 2018, ApJS, 239, 27, doi: 10.3847/1538-4365/aae37a
  * Strom et al. (2017) Strom, A. L., Steidel, C. C., Rudie, G. C., et al. 2017, ApJ, 836, 164, doi: 10.3847/1538-4357/836/2/164
  * Suess et al. (2022a) Suess, K. A., Leja, J., Johnson, B. D., et al. 2022a, ApJ, 935, 146, doi: 10.3847/1538-4357/ac82b0
  * Suess et al. (2022b) Suess, K. A., Kriek, M., Bezanson, R., et al. 2022b, ApJ, 926, 89, doi: 10.3847/1538-4357/ac404a
  * Tacchella et al. (2022) Tacchella, S., Conroy, C., Faber, S. M., et al. 2022, ApJ, 926, 134, doi: 10.3847/1538-4357/ac449b
  * Taylor et al. (2015) Taylor, E. N., Hopkins, A. M., Baldry, I. K., et al. 2015, MNRAS, 446, 2144, doi: 10.1093/mnras/stu1900
  * Taylor et al. (2022) Taylor, L., Bezanson, R., van der Wel, A., et al. 2022, ApJ, 939, 90, doi: 10.3847/1538-4357/ac9796
  * Thomas et al. (2005) Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673, doi: 10.1086/426932
  * Thomas et al. (2010) Thomas, D., Maraston, C., Schawinski, K., Sarzi, M., & Silk, J. 2010, MNRAS, 404, 1775, doi: 10.1111/j.1365-2966.2010.16427.x
  * Thomas et al. (2017) Thomas, R., Le Fèvre, O., Scodeggio, M., et al. 2017, A&A, 602, A35, doi: 10.1051/0004-6361/201628141
  * Tojeiro et al. (2009) Tojeiro, R., Wilkins, S., Heavens, A. F., Panter, B., & Jimenez, R. 2009, ApJS, 185, 1, doi: 10.1088/0067-0049/185/1/1
  * Torrey et al. (2018) Torrey, P., Vogelsberger, M., Hernquist, L., et al. 2018, MNRAS, 477, L16, doi: 10.1093/mnrasl/sly031
  * van der Wel et al. (2016) van der Wel, A., Noeske, K., Bezanson, R., et al. 2016, ApJS, 223, 29, doi: 10.3847/0067-0049/223/2/29
  * van der Wel et al. (2021) van der Wel, A., Bezanson, R., D’Eugenio, F., et al. 2021, ApJS, 256, 44, doi: 10.3847/1538-4365/ac1356
  * van der Wel et al. (2022) van der Wel, A., van Houdt, J., Bezanson, R., et al. 2022, ApJ, 936, 9, doi: 10.3847/1538-4357/ac83c5
  * van Houdt et al. (2021) van Houdt, J., van der Wel, A., Bezanson, R., et al. 2021, ApJ, 923, 11, doi: 10.3847/1538-4357/ac1f29
  * van Rossum (1995) van Rossum, G. 1995, CWI Technical Report, CS-R9526
  * Vulcani et al. (2014) Vulcani, B., Bamford, S. P., Häußler, B., et al. 2014, MNRAS, 441, 1340, doi: 10.1093/mnras/stu632
  * Wake et al. (2012) Wake, D. A., van Dokkum, P. G., & Franx, M. 2012, ApJ, 751, L44, doi: 10.1088/2041-8205/751/2/L44
  * Webb et al. (2020) Webb, K., Balogh, M. L., Leja, J., et al. 2020, MNRAS, 498, 5317, doi: 10.1093/mnras/staa2752
  * Weisz et al. (2008) Weisz, D. R., Skillman, E. D., Cannon, J. M., et al. 2008, ApJ, 689, 160, doi: 10.1086/592323
  * Weisz et al. (2011) Weisz, D. R., Dalcanton, J. J., Williams, B. F., et al. 2011, ApJ, 739, 5, doi: 10.1088/0004-637X/739/1/5
  * Wild et al. (2016) Wild, V., Almaini, O., Dunlop, J., et al. 2016, MNRAS, 463, 832, doi: 10.1093/mnras/stw1996
  * Wu et al. (2018a) Wu, P.-F., van der Wel, A., Bezanson, R., et al. 2018a, ApJ, 868, 37, doi: 10.3847/1538-4357/aae822
  * Wu et al. (2018b) Wu, P.-F., van der Wel, A., Gallazzi, A., et al. 2018b, ApJ, 855, 85, doi: 10.3847/1538-4357/aab0a6
  * Wuyts et al. (2009) Wuyts, S., Franx, M., Cox, T. J., et al. 2009, ApJ, 696, 348, doi: 10.1088/0004-637X/696/1/348
  * Zahid et al. (2018) Zahid, H. J., Sohn, J., & Geller, M. J. 2018, ApJ, 859, 96, doi: 10.3847/1538-4357/aabe31

Figure A1: Comparing 10$\%$ (top row), 50$\%$ (middle row) and 90$\%$ (bottom
row) formation times for star-forming (left), quiescent (center) and all
(right) galaxies from different flavors of spectro-photometric runs using
Bagpipes and Prospector. Our results are robust to the photometric information
provided in the fits and are defined by the high-resolution spectral
information. Error bars represent 16th - 84th percentile ranges.

Figure A2: Extended version of Figure 5 showing recovered SFHs of lower
stellar mass galaxies for Bagpipes and Prospector in SFR ($M_{\odot}/yr$) and
log(sSFR) ($yr^{-1}$) plane. Top two rows show galaxies in lowest stellar mass
bin 10.5 < log($\mathrm{M_{*}/M_{\odot}}$) < 11 for quiescent (a) and star-
forming (b) sample. Similarly, bottom two row shows the same for galaxies in
11 < log($\mathrm{M_{*}/M_{\odot}}$) < 11.5 stellar mass bin. Note that the
vertical axis limits of SFHs are different for the two populations to zoom
into the features that are prominent in sSFR plane.

Our primary spectro-photometric analysis was conducted on deep LEGA-C spectra
with BvrizYJ broadband photometric data. To investigate whether this modeling
choice impacts our final results, we fit the full LEGA-C primary sample with 2
more flavors in Bagpipes - (1) spectra + vriz bands and (2) spectra + 20
photometric bands (from B to 24 microns as described in Section §2.1). We
compare the median $t_{10}$, $t_{10}$ and $t_{90}$ formation times from these
flavors with those presented in this paper from Bagpipes and Prospector using
BvrizYJ bands only in Figure A1. Error bars represent the population scatter
within the samples of that mass bin. It is evident that the derived formation
times do not depend strongly on the photometric information provided to the
modeling.

The population scatter is higher for earlier formation times compared to later
ones with Prospector $t_{10}$ dominated by the priors on the very first time
bin. For quiescent population, we do see highest variation in lowest mass bin
10.5 < $\mathrm{log(M/M_{\odot}})$ < 11, primarily driven by low number of
samples ($\sim$290) in that mass range.

Figure A2 is an extended version of Figure 5 and shows SFHs of quiescent and
star-forming population combined in the other two low stellar mass bins e.g.
10.5 < log($\mathrm{M_{*}/M_{\odot}}$) < 11 and 11 <
log($\mathrm{M_{*}/M_{\odot}}$) < 11.5. All symbols and labels are same as
Figure 5. Top two rows show lowest mass bin for quiescent and star-forming
population whereas bottom two rows show middle mass bins. Note that the
vertical axis limits of SFHs in $M_{\odot}/yr$ unit are different for the two
populations to zoom into the features that are prominent in sSFR plane.