New Observational $H(z)$ Data from Full-Spectrum Fitting of Cosmic Chronometers in the LEGA-C Survey

In this work, we perform a full-spectrum fitting of 350 massive and passive galaxies selected as cosmic chronometers from the LEGA-C ESO public survey to derive their stellar ages, metallicities, and star-formation histories. We extensively test our results by assessing their dependence on the possible contribution of dust, calibration of noise and signal, and the use of photometric data in addition to spectral information; we as well identify indicators of the correct convergence of the results, including the shape of the posterior distributions, the analysis of specific spectral features, and the correct reproduction of the observed spectrum. We derive a clear age-redshift trend compatible with the aging in a standard cosmological model, showing a clear downsizing pattern, with more massive galaxies being formed at higher redshift ($z_f\sim2.5$) with respect to lower massive ones ($z_f\sim2$). From these data, we measure the differential aging of this population of cosmic chronometers to derive a new measurement of the Hubble parameter, obtaining $H(z=0.8) = 113.1 \pm 15.1 (\mathrm{stat.}) ^{+29.1}_{-11.3} (\mathrm{syst.})\ \mathrm{ km\ s^{-1}\ Mpc^{-1}}$. This analysis allows us for the first time to compare the differential ages of cosmic chronometers measured on the same sample with two completely different methods, the full-spectrum fit (this work) and the analysis of Lick indices, known to correlate with the age and metallicity of the stellar populations \citep{Borghi2022a}. Albeit an understood offset in the absolute ages, the differential ages have proven to be extremely compatible between the two methods, despite the very different data, assumptions, and models considered, demonstrating the robustness of the method.


INTRODUCTION
Since the discovery of the accelerating expansion of the Universe (Riess et al. 1998;Perlmutter et al. 1999), the cosmological community has been working to understand the mechanism of this expansion. Modern cosmology postulates that dark energy, an unknown form of energy with negative pressure, is driving the accelerated expansion of the late Universe and that the gravitational effect of Cold Dark Matter (CDM) shapes the large-scale structure of the Universe, a model dubbed ΛCDM. Numerous cosmological probes and observations, including the Cosmic Microwave Background (CMB, eg. Smoot et al. 1992;Bennett et al. 2003;Planck Collaboration et al. 2014;Swetz et al. 2011;Carlstrom et al. 2011;Planck Collaboration et al. 2020), Baryon Acoustic Oscillations (BAO, eg. Percival et al. 2001;Cole et al. 2005;Eisenstein et al. 2005), Type Ia supernovae (SNe, eg. Sullivan et al. 2011;Suzuki et al. 2012;Betoule et al. 2014;Scolnic et al. 2018), weak gravitational lensing (Bartelmann & Schneider 2001), and cluster counts (Allen et al. 2011), have been proposed and extensively studied to determine the Universe's large-scale structure and evolution. After more than twenty years of unremitting efforts, we are now in the golden age of precision cosmology, with measurements and constraints on cosmological parameters reaching the percent level.
The Hubble constant H 0 has long been a critical observable of observational cosmology (Freedman & Madore 2010), and its value is directly related to our current estimate of the Universe's age. However, the two probes that represent the most precise level of measurement today, SNe and CMB, produce significant discrepancies beyond the order of 4σ (for an extensive review, see Di Valentino et al. 2021). The increase of observational evidence supporting this discrepancy between observations of the early and late Universe, has undoubtedly kicked off a crisis in modern cosmology (Verde et al. 2019;Davis 2019;Riess 2020;Abdalla et al. 2022). At the moment, there are suggested theories to try to explain it, even if not definitive (Di Valentino et al. 2021). Alternative cosmological probes (Moresco et al. 2022) can play an important role in obtaining additional independent, high-precision measurements to assess the current Hubble tension's reliability. It also becomes evident that a single probe is not adequate to constrain the properties and evolution of the Universe accurately and completely. The Hubble parameter H(z) is the physical quantity that most directly describes the history of the Universe's expansion, and its measurement has advanced significantly over the last decade or so. We can not only reconstruct these H(z) measurements to extrapolate H 0 at zero-redshift, shedding light on a different path to explore the crisis, but also enhance our ability to understand the nature of dark energy, which dominates the late Universe just covered by the observable range of H(z).
Existing observational H(z) data (referred to as OHD, see e.g. Zhang et al. 2010;Ma & Zhang 2011) are mainly based on two probes: the differential age method and the radial BAO size method (Benítez et al. 2009). The former can be obtained with Cosmic Chronometers (CC, Jimenez & Loeb 2002;Moresco et al. 2018Moresco et al. , 2020) by measuring the differential age-redshift relation of massive and passive galaxies throughout the Universe. Any systematic offset introduced by the galaxy age measurement method will be canceled out when deriving the differential age. A total of 32 H(z) measurements have been obtained (Jimenez et al. 2003;Simon et al. 2005;Stern et al. 2010;Zhang et al. 2014;Moresco 2015;Moresco et al. 2016;Ratsimbazafy et al. 2017;Borghi et al. 2022b) and are currently widely used to test cosmological models. These measurements are regarded as cosmological model-independent since the principle is not dependent on the choice of cosmological models. The second method is based on the inverse proportionality between H(z) and the differential radial (comoving) distance, which can be traced by measuring the radial size of BAO features at different red-shifts. This method, however, requires knowledge of the comoving BAO scale (r BAO ), which is derived from the CMB measurements. This fact makes this probe not fully cosmology-independent, since typically in the derivation of the sound horizon scale from CMB a cosmology model is assumed. Additionally, gravitational waves can be used as standard sirens (Schutz 1986;Holz & Hughes 2005;Abbott et al. 2017), to study H(z), with promising perspectives for the next decade (Farr et al. 2019). Finally, the phenomenal growth of Fast Radio Burst (FRB) observations also expands the H(z) measurement possibilities (Wu et al. 2020).
Selecting a pure passive sample and measuring the age difference between galaxies are the two bases of the CC method. Various strategies have been proposed to distinguish 'passive' from 'star-forming' galaxies, including morphological selections of spheroidal systems (following Hubble 1936), cuts on color-color diagrams (e.g. UVJ, Williams et al. 2009;NUVrJ, Ilbert et al. 2013) or on a color-mass diagram (e.g. Peng et al. 2010), and Spectral Energy Distribution (SED) fitting (e.g. Pozzetti et al. 2010a). Combining multiple criteria and maximizing the overlap of complementary information (photometric and spectroscopic) result in a significantly more effective method of selecting a pure sample (Moresco et al. 2013(Moresco et al. , 2018Borghi et al. 2022a). While spectral line analysis enables us to obtain extremely precise redshift values, the situation is much more complicated in determining the age, which can not be directly observed, but can be estimated using photometry (SED), single spectral regions (e.g., D4000, Lick indices), or the entire spectrum features (full-spectral fitting). However, each of these methods may suffer various systematics caused by parameter degeneracies. Moresco (2011) explains that SED-fitting, which is commonly used to derive galaxies' ages, is incapable of fully breaking the age-metallicity degeneracy; also, the age degenerates to τ in the delayed exponential Star-Forming History (SFHs). Moresco et al. (2011) proposed an innovative method that consists in not using the age but rather a direct spectroscopic observable (the 4000Å break) to measure H(z), making the decoupling of systematic and statistical error easier. Most recently, Borghi et al. (2022b) obtained for the first time a new H(z) measurement using the Lick indices method. Their analysis takes advantage of the high signal-to-noise (S/N) spectroscopic data of the LEGA-C (DR2) survey (Straatman et al. 2018) of galaxies at z ∼ 0.7. However, the Lick indices method does not allow to flexibly study the galaxies' star formation histories, which are useful to better exclude possible biases in the derived age-redshift relation and, therefore, on H(z).
For the purpose of optimizing the set of the spectral absorption features for Lick indices fitting, Borghi et al. (2022b) only make use of a sub-sample (140 galaxies) of the selected passive sample, while the full-spectrum fitting is not subject to this issue.
In this paper, we perform full-spectrum fitting to derive the ages and star-formation histories of passive galaxies in LEGA-C DR2, then use them as cosmic chronometers to obtain a new H(z) measurement. The dataset is introduced in Section 2, and the fundamental principles and details of the full-spectrum fitting in Section 3. In Section 4 we present the results on physical parameters and the strategies to improve the performance of their estimation. In Section 5 we detail the procedure for applying the CC method, presenting and discussing the final H(z) measurement results. The conclusions are presented in Section 6.

DATA
In this section, we describe the spectroscopic and photometric data used in this analysis, and the selection criteria adopted to select the sample of cosmic chronometers.
Spectroscopic Data -The spectroscopic data are taken from the Large Early Galaxy Astrophysics Census (LEGA-C), an ESO 130-night public survey of ∼ 3200 K s -band selected galaxies conducted with VLT/VIMOS (Le Fèvre et al. 2003) on the Very Large Telescope. The 20-hour long integrations produce continuum spectra with an average S/N ∼ 21.8 per pixel (0.6Å) for massive galaxies (M 10 11 M ). The second data release (LEGA-C DR2, Straatman et al. 2018) includes 1988 spectra in the redshift range 0.6 z 1.0 covering the observed wavelength range ∼ 6300−8800Å, with an effective spectral resolution of R ∼ 3500. We add to the LEGA-C dataset the spectral indices measurements from Borghi et al. (2022a), providing a catalog of Lick indices measurement including also the recent CaII H/K diagnostic (a useful tracer of recent episodes of star formation, Moresco et al. 2018).
Photometric Data -One of the advantage of the LEGA-C sources is that, being observed in the COS-MOS field, a wealth of multi-wavelength photometric observations are available (e.g., Muzzin et al. 2013;Laigle et al. 2016;Weaver et al. 2022). In this work, following Straatman et al. (2018) we adopt the Ultra Deep Survey with the VISTA telescope (UltraVISTA) photometric catalog from (Muzzin et al. 2013). We use a total of 21 photometric bands, namely , IB427, IB464,  IA484, IB505, IA527, IB574, IA624, IA679, IB709,  IA738, IA767, IB827 ,u,V,zp,Y,J,H,Ks,ch1,ch2. For a given filter x, we compute the total flux f x,tot by applying the equation where f Ks,tot is the total K s band flux from SExtractor's flux auto that has been corrected using the growth curve of the point-spread function (PSF) stars and f Ks is the specific Ks-band flux (see Muzzin et al. 2013). The SEDs of the catalogues are in good agreement, but differences in calibration and measurement precision may affect the age-z relation. We will further explore the use of different photometric calalogues in a follow-up analysis.  Borghi et al. (2022a) selected a pure sample of 350 passive galaxies in LEGA-C DR2, minimizing any residual contamination from star-forming outliers. The distribution of some key parameters describing this population is shown in Figure 1. This passive sample has a median redshift of z = 0.735 with two peaks around z ∼ 0.745 and z ∼ 0.839. The median values of the σ and log 10 (M /M ) distributions increases from 165.7 km s −1 (10.72) to 205.7 km s −1 (10.95) respectively. Its specific star formation rate (sSFR) distribution has a median logarithmic value of log 10 (sSFR) = −12.2 yr −1 , which is ∼ 1 dex lower than what is typically used to define a galaxy 'passive' (e.g., Pozzetti et al. 2010b).

METHOD
Full Spectrum Fitting -To perform the full-spectrum fitting we use the Bagpipes code developed by Carnall et al. (2018). Bagpipes models the observed spectrum of a galaxy f obs into f H (Θ) based on an hypothesis H of the physics involved described by parameters Θ. The posterior distribution P (Θ | f λ , H) obtained from the Bayes theorem is sampled with the nested sampling algorithm Mu-liNest (Buchner et al. 2014). Here H includes the modelling of the star formation rate SFR(t Ui ), the simple stellar-population SSP(t i , λ, Z j ), and the neutral and ionized interstellar medium (ISM) radiative transmission function T 0 (t i , λ) and T + (t i , λ), which are used to  Figure 1. Distributions of four key parameters for the LEGA-C DR2 parent sample (blue) and the 350 passive galaxies analyzed in this work (yellow). The redshift (z), stellar velocity dispersion (σ ) and stellar mass (M ) are taken from LEGA-C DR2, while the specific star formation rate (sSFR) are from COSMOS2015. The arrows mark out the median values.
simulate the luminosity function of a galaxy, where t U , t, λ, Z are the cosmic time, the age of the stellar population, the wavelength of spectral line and stellar metallicity, respectively, and the subscripts i and j denote summations for all the age bins and SFH components, respectively. H also includes the modelling of the intergalactic medium (IGM) radiative transfer to finally simulate the observed flux, where λ obs = (1 + z obs )λ, D L (z obs ) is the luminosity distance and T IGM is the transmission function of the IGM. The nebular emission lines and continuum come from pre-computed CLOUDY (Ferland et al. 2017) grids with only one free parameter, the logarithmic ionization parameter (log 10 (U )). We apply the Charlot & Fall (2000) model out of the four choices of dust attenuation models (see detail descriptions in Carnall et al. 2018) that Bagpipes provides. The likelihood function can be written in a logarithmic form as, where σ is the observation error of the fluxes, and here sums over all i-th wavelength pixels.
In addition to the spectrum, Bagpipes allows the inclusion of photometric data points in the fit, thus enabling modelling a galaxy SED on a wide wavelength range, from far-ultraviolet to microwave regimes. Another significant advantage is that it is possible to adaptively test different SFH choices (e.g., single burst, constant, exponentially declining, as well as a combination of them). Models within the Bagpipes code are resampled in an age grid ∆t i based on the SSP model generated using the 2016 version of the Bruzual & Charlot (2003) (BC16) models. Bagpipes is structured around three core classes, which are galaxy for loading observational data, model galaxy for generating model galaxy spectra and fit for fitting models to observational data.
The code is open source and publicly avaliable 1 . To extract parameters values and associated uncertainties from the posterior distributions, we adopt the median and 16-84 th percentiles, respectively.
Star Forming History Choice -Most of the CC analyses, including Borghi et al. (2022b), assume single burst star formation histories as an ideal simplification of the real SFH. This model assumes that the total mass of a galaxy suddenly formed at a specific cosmic time, which is characterized by a delta function SFR(t U ) ∝ δ(t U ). A more realistic SFH model is necessary to test the robustness of the age-z relations and the H(z) obtained. In this work, we extend this analysis by testing two other wellestablished SFH models, namely the double power law (DPL) and the delayed exponentially declining (DED) model based on the CC sample that Borghi et al. (2022a) compiled. The DPL model separates the rising and declining phases of the SFH using two separate power-law slopes, SFR(t U ) ∝ (t U /τ ) α + t U /τ ) −β −1 , where α is the falling slope, β is the rising slope and τ is related to  Figure 2. Full spectrum fitting results for an example galaxy (ID=215424) obtained with spectroscopic data alone (top panel) and adding photometry (second and third panel). The observational spectra from LEGA-C DR2 (Straatman et al. 2018) and photometric data points from the UltraVISTA catalogue of Muzzin et al. (2013) are shown in light blue and blue, respectively. The best-fit BAGPIPES spectra and photometric points from fitting spectroscopy and spectroscopy+photometry are shown in black and orange, respectively. The bottom panel shows the best-fit BAGPIPES SFR (as a function of lookback time) correspondingly, the solid curves are the median posteriors while the shade regions are the 1σ confidence regions, the horizontal axis is the lookback time since t(z obs ).
(but not the same as) the peak time. The DED model assumes that the star-formation starts at some time T 0 and increases gradually to its peak, after which it declines exponentially with some timescale τ , As detailed in this Section 3, we fit the selected sample separately using the DPL and DED model and get compatible median reduced chi-squared values for the spectrum of χ 2 ν = 1.96 and χ 2 ν = 1.92. According to the principle of the Ockham's Razor, we choose the one fewer free parameter model -DED (Equation 6) for our following analysis.
Removing Cosmological Prior -The age of a stellar population, t, is defined as the look-back time between its observed redshift and the beginning of its star- In fact, the cosmic time (t U ) at a given redshift is not a direct observable, we can only calculate its value based on a cosmological model. Bagpipes use ΛCDM as its default cosmological model, with the default parameters Ω M = 0.3, Ω Λ = 0.7 and H 0 = 70 km s −1 Mpc −1 . In principle, when the age of a galaxy exceeds the age of the Universe, it is reasonable to consider it non-physical.
To address this issue, Bagpipes assumes that the star formation rate SFR(t Ui ) = 0 when the retrieved age is larger than the estimated age of the Universe at the given redshift. While this assumption is typically neglected in galaxy evolution studies, being of relative interest for the results, imposing an upper limit on the retrieved age based on a cosmological model is to be strictly avoided in our analysis. In particular, such a prior could induce cosmological biases in the age estimates and circular arguments in the derivation of the Hubble parameter, since if the ages of the oldest objects are set to the age of the Universe of the reference cosmological model the method would artificially provide, by definition, the reference cosmology. Fortunately, we can avoid the above situation by releasing the upper limit for t U in Bagpipes to a value that our sample galaxies cannot exceed, such as 20 Gyr at all redshifts. This modification changes the upper boundary of the age sampling, without affecting the sampling grid. Galaxies' age can be easily retrieved by subtracting the formation time with the new upper limit we set, and the age-z slop will not be affected by the cosmological assumptions. Besides, we notice that the cosmological model is also used when calculating the D L in Equation 4. The D L doesn't interact with t U in the rest of the code and because of the negligible dependence of dH(z) ∼ O 2 (dD L ), it is acceptable to ignore the issue of D L affected by the choice of cosmological model (see figure 1 in Jimenez & Loeb 2002). In conclusion, with these modifications we "erase" the effect of the cosmological prior on the galaxy age estimation in the original Bagpipes code.
Adding Photometric Data -The spectroscopy covers a relatively narrow wavelength range of the galaxies' entire spectrum compared with the photometry. Fitting spectroscopic data alone, due to the lack of enough information, is incapable of fully modelling the line features and breaking the degeneracy between parameters, especially the age-τ degeneracy in our analysis, as well as the CaII H and K lines that are essential for diagnostic of passive galaxies. Adding photometric data will improve the performance of fitting by providing additional information. We employ the Ks selected UltraVISTA photometries in our analysis as detailed in Section 2.
In Figure 2 we show the full spectrum fitting results from an example galaxy (ID 215424) obtained with spectroscopic data alone and adding photometry. We observe that in the posterior the CaII K line is less deep than the H line for the spectroscopic (only) fitting, contradicting the observational data. On the contrary, this feature is well reproduced after including the photometric data in the fit. This same behavior is observed, in general, for the entire sample. In particular the median percentage difference between the observed and reconstructed H/K (see also Table 1 and Section 4.1 for a more extended discussion) is significantly reduced from the value of 11.93 ± 6.76 (%) to 6.46 ± 4.34 (%). From the histograms in Figure 3, we observe both long tails of derived ages and τ distributions for the spectroscopic (only) fitting, indicating the existence of the age-τ degeneracy, while the tails are significantly suppressed after adding the photometric fitting.
Exploring Parameter Space -In our analysis, we fully explore all the possible parameter space consist of the physical properties of the galaxy, the parameters of modelling the noise and the dust attenuation, together with the calibration of the fluxes. We consider two diagnostics, the reduced chi-square (χ 2 ν ) and the percentage difference of H/K (see topical description in Section 4.1) to quantify the agreement between the posterior and observed spectrum and to verify the improvement of the fit when adding more information on the parameters.
Fitting redshift is not necessary, for each galaxy we set its z to the value of the high precision spectroscopic redshift derived from numerous high S/N absorption features observed the LEGA-C DR2 spectra. We take Uniform priors in reasonably wide parameter spaces for the age (t), star formation timescale (τ ), stellar metallicity (Z), and the log of the mass formed log 10 (M form ) as described in Table 1. The stellar velocity dispersion σ is a direct observable derived from the line-broadening of the spectrum, we test the effect of using a wide uninformative prior σ ∼ U[0, 400] and narrower Gaussian prior of σ ∼ G[σ ; err σ ], set by the measurements provided in the LEGA-C DR2 catalog. We also test the impact of varying the observed flux calibration according to a second-order Chebyshev polynomial with coefficients' prior set to P 0 ∼ G[1, 0.25], P 1 ∼ G[0, 0.25] and P 2 ∼ G[0, 0.25] for each of the orders respectively. We also test the white noise model as an additional component, by adopting a uniform prior for the logarithmic white noise scaling parameter log 10 (S noise ) ∼ U[0, 10].
Finally, even if we expect it to be negligible for this sample, we additionally model the dust effect by assuming the Calzetti et al. (2000) attenuation curve and a Gaussian prior on the absolute attenuation at 5500Å, A V ∼ G[1, 0.25], a multiplicative factor of η = 2 for stars in birth clouds, and an attenuation power-law slope of n ∼ G[0.7, 0.3] in the range n ∈ [0.3, 2.5]. For further details on each model component please check Carnall et al. (2018). When doing the analyses, we do not expect the nebular and dust emission modellings to have significant impacts on our results, since our galaxies were selected to have negligible or no emission lines contribution (Borghi et al. 2022b).
We compare the results obtained from using or not the aforementioned components and find all the χ 2 ν are compatible (on the second digits). The modelling of H/K (see topical description in Section 4.1) are significantly improved (on the first digits of ∆ H/K ) when not using the calibration or modelling the noise. Since the LEGA-C DR2 spectra are flux-calibrated using Ul-traVISTA's photometric SEDs, we prefer not to impact on the data with an additional calibration. But the other option is also acceptable and we get compatible H(z) measurements for them during our following step analyses. We decide to use as a baseline the set of parameters that better reproduces the CaII H/K feature (the last row of Table 1, since it has been proven to be a very powerful and important diagnostic to trace the purity of CC samples (see Moresco et al. 2018;Borghi et al. 2022a;Moresco et al. 2022) and having a significantly different H/K in the reconstructed spectrum would mean not to correctly reproduce the behaviour of our data, possibly resulting in biases in the results; moreover, we notice that all the models present a compatible χ 2 ν , but this model has also the advantage to avoid additional calibration, since the LEGA-C DR2 spectra are flux-calibrated using UltraVISTA's photometric data.

ANALYSIS
In this section, we present our analysis of the physical parameters derived. We start by defining three diagnostic criteria to assess the reliability of our results, and flag the constraints not properly converged. We will explore their use to improve the robustness of the derived parameter, if needed, concluding by presenting our baseline results on which will be based the cosmological analysis.

Breaking The Degeneracies
In Figure 3 we present the distributions of the derived parameters (z, σ , stellar age, and τ ) obtained from our analysis of both the spectroscopic data alone and from the fit of the spectroscopic and photometric data combined. We start by observing that the pure spectroscopic case presents significant tails in the age and τ distributions, with ages larger than the age of the Universe and τ up to the maximum value allowed by the prior. We will therefore define a set of diagnostic criteria to check the accuracy of our results, apply them to our constraints, and verify if and how they can impact the nonphysical large values just discussed.
Improperly Converged Constraints -The values estimated from the posterior distributions obtained from Bagpipes may not be fully reliable depending on several issues. As an example, this is the case when the   posterior distribution exhibits multi-modal peaks, or is very skewed toward the edge of the domain allowed by the priors. Considering the large number of galaxies and combinations of parameters explored, we develop two efficient automatic algorithms to assist us in recognizing these issues and subsequently flagging the corresponding galaxies.
We identify the multi-modality by counting the convex inflection points of the one-dimensional posterior dis-tribution functions (PDFs) after applying a Gaussian kernel smoothing function with a bandwidth with size 12% of the full posterior range. This technique is robust against spurious detection of close-by peaks in distributions that are not significantly multi-modal. We validate this method through visual inspection, verifying that PDFs with zero inflection points are actually very flat and uninformative, and those with more than one are significantly multi-modal. Therefore, we keep those PDFs with only one convex inflection point. Throughout the paper, we refer to this specific flag as unimodal.
To determine whether a posterior is very skewed toward the parameter space boundary, the most intuitive method is to examine whether the estimated upper and lower bounds exceed the parameter space. Here, we adopt a more flexible approach based on the skewness computed as the adjusted Fisher-Pearson standardized moment coefficient of the PDF: if the absolute value of the skewness is > 1, the chain is considered highly skewed or asymmetric, indicating that it is converging toward the edge of the parameter space. We validate this algorithm by visually inspecting the onedimensional posterior PDFs of a fully fitted catalogue. Throughout the paper, we refer to this specific flag as not truncated.
In our Bayesian analysis, there are numerous free parameters involved, including the model components of noise and dust as discussed in Section 3 for which in principle one may apply these automatic inspection techniques. In this work, we consider only four key parameters associated with galaxies' physical properties, namely, age, τ , Z, and M . The final unimodal (or not truncated) flag is then taken from the intersection (AND logic) of the individual flags obtained from all four parameters.
We find that the not truncated flag removes more galaxies in the lower age regime, by increasing the median of the age from t = 4.1 ± 1.5 to 5.6 ± 2.4 Gyr and the value of τ from τ = 0.5±0.3 to 0.7±0.3 Gyr respectively, while the unimodal selection negligibly affects the shape of the posterior distribution of the derived parameters, as well as their median values. Though potentially improving its reliability, the not truncated and unimodal flags shrink the sample by approximately 20% and 50%-60%.
Poorly modelled CaII H/K -The CaII K and H lines, centered respectively at 3934 and 3969Å rest frame (see Figure 4), are two prominent features in galaxy spectra. In galaxies dominated by an old stellar population, it is usually found that the K line is deeper than the H line, this being opposite in presence of young star-forming components (see figure 5 in Moresco et al. 2018). In our analysis, we adopt the definition of H/K introduced by Fanfani (2019) that consists in measuring the ratio of two pseudo-Lick indices CaII K and H, i.e. H/K = I H /I K . This technique is less sensitive to potential bias introduced by noise peaks in the spectrum with respect to using the H and K flux minima, i.e. |H/K| min = F min (H)/F min (K).
This diagnostic has been used to test the presence or absence of a contaminant population, therefore describ-ing the purity of the selected sample (Moresco et al. 2018;Borghi et al. 2022a). In particular, Borghi et al. (2022a) found that H/K < 1.2 is safely equivalent to |H/K| min > 1, and H/K < 1.1 well reproduces other selection criteria including NUVrJ (Ilbert et al. 2013) and sSFR/yr < −11 cut. They tested that the current sample of passive galaxies has a typical value of H/K = 0.96 ± 0.08, validating the purity of the selection. We, as well, do not observe any correlation between galaxies' properties (especially age) and the H/K in both the observed and posterior spectra, indicating that no significant contribution from young stellar component is present in our passive galaxies sample. This also excludes the possible presence of galaxies that experienced recent rejuvenation events as found in Chauke et al. (2019).
It is, therefore, plausible to consider that the posterior spectra which do not adequately reproduce this feature are not completely appropriate fits, and should be excluded. We calculate the H/K ratios in the observational LEGA-C DR2 spectra and in our posterior spectra using PyLick 2 . Then, we associate a flag with removing galaxies from the sample according to the following criterion: if the discrepancy between the observed and the inferred value is greater than 10%, we consider that the feature has not been well reproduced, resulting in a potential deviation in the age estimation of the galaxy. We define the flag corresponding to this criterion H/K well modelled. Applying this flag, we observe a considerable reduction in the long-tail shapes of the age and τ distributions in Figure 3, which clearly demonstrates the importance and validity of this criterion.
Inconsistency between Photometric and Spectroscopic Data -The calibration of spectroscopic and photometric data is a complex process, and systematic differences between data obtained using different calibration pipelines are possible. To ensure a proper combination of the different sets of data, we must compare the photometric points (when available) in the wavelength region covered by the observed spectrum to the flux of the spectrum at the effective position of the corresponding photometric filters(see description in Section 2); if the difference between the two is significant, it means that there is an inconsistency between the photometric and spectroscopic points, and therefore these data cannot be fit jointly. We define the median absolute pull (MAP) between the photometric fluxes and the spectral fluxes in the form of where the spectral fluxes are averaged over extremely narrow wavelength windows (10 pixels ∼ 6Å) centered at the effective wavelength of the filters. We conservatively keep only the galaxies with a MAP ≤ 1. Such flagging is stringent enough to rule out most of the largest potential deviations. Throughout the paper, we refer to this specific flag as better SEDs. Applying this flag, we observe the median of the age and σ are equivalent to the values when applying the H/K well modelled flag, even though not significant, suggesting the inconsistency has a negative effect on simultaneously fitting spectroscopic and photometric data.
Inspecting Figure 3, we notice how in the case when only the spectrum is fitted, the various flags defined are crucial to significantly help reduce the tails of nonphysically large ages (age > 7 Gyr), and as well the largest values of τ (τ > 1 Gyr), that were galaxies severely affected by the age-τ degeneracy fitted with a larger value of age and τ . This effect is in particular evident by comparing the median values of each distribution as a function of the different flags applied, as also reported in the figure. At the same time, analyzing the distribution of the fit obtained from the joint analysis of spectrum and photometry, we notice how the impact of the flagging is significantly smaller, and that the addition of the photometric bands per-se helps in reducing the degeneracy between parameters and obtaining well-converged fits, with an extremely negligible fraction of points at high ages and τ . It remains also interesting to notice that despite wide and uninformative priors on age and τ , the derived best fits confirm the fact that these objects have been selected extremely accurately, and that they are old objects formed over relatively small timescales.
In conclusion, we decide to keep as our baseline for the analysis the fit that includes spectroscopic and photometric data without including any flagging, maximizing in this way both the accuracy of the results and the final statistics. In Section 5.2 we will discuss and quantify the impact of applying these flags on our cosmological result. Each galaxy is color-coded by its stellar velocity dispersion (σ ). The first two and the last panel span the entire parameter space explored (for a more detailed description, see Section 3. Figure 4 presents the final constraints to the stellar population properties, namely the best-fit stellar age, metallicity (Z), mass (M * ), and star formation timescale (τ ), by assuming a delayed exponentially declining star formation history (see the last row of Table 1 for more details on the model ingredients and adopted priors). We use the light-weighted properties instead of the mass-weighted ones, considering that the latter are more sensitive to the choice of the star-forming history parametrization (Conroy 2013). Each galaxy is colorcoded by its stellar velocity dispersion. Interestingly, even if do not bind the upper age value with a cosmological prior, we find that this population of galaxies qualitatively follows the descending trend predicted from the ΛCDM model assuming the Planck Collaboration et al. (2020) parameters, with less than 10% of galaxies higher than the reference ΛCDM boundary. The median age is t = 3.60 ± 0.82 Gyr. For comparison purposes only, when assuming a baseline ΛCDM model, this value corresponds to a typical formation time of ∼ 3 Gyr after the Big Bang, or z f ∼ 2.5, in agreement with a wealth of literature data (e.g. Gallazzi et al. 2014;Carnall et al. 2019;Belli et al. 2019;Carnall et al. 2022;Tacchella et al. 2022). In Section 5 we further study the age(z) relation and its trends with σ .

Galaxy Properties
For the first time, we are able to quantitatively study the star-formation timescale τ -redshift relation for this sample of cosmic chronometers. Even if we adopt a wide Uniform prior on τ ∼ U[0, 2] Gyr, we find typical τ = 0.36 ± 0.13 Gyr, with about 80% of the galaxies having τ < 0.5 Gyr. Most importantly, from this analysis we find no significant dependence on z. This is an additional confirmation that the current sample of passive galaxies is very homogeneous in its physical properties over different z. In addition, we do not find a statistically significant trend of τ with σ . This is not in contradiction with the idea that more massive galaxies formed in shorter timescales (as expected from the downsizing-scenario; Cowie et al. 1996). On the contrary, we are selecting the very massive and passive envelope of objects, so that we expect shortest τ and no correlation with mass. This may also explain the reason why our SFHs are shorter the those derived by Chauke et al. (2018) on LEGA-C quiescent galaxies. But for a definitive answer, the impact of different SFH assumptions must be further assessed. Similar values of τ for quiescent galaxies at intermediate redshifts were also obtained by Pacifici et al. (2016), Carnall et al. (2019), and Tacchella et al. (2022) We find stellar metallicities with slightly undersolar values, Z/Z = 0.84 ± 0.41, in agreement with Borghi et al. (2022a) who find typical Z/Z ∼ 1.1 using Lick indices. This result further strengthens the idea that it exists a population of massive and passive galaxies which, at least up to z ∼ 0.8, does not evolve significantly in its metal content and has values similar to those of their local counterparts (see also Thomas et al. 2011;Gallazzi et al. 2014;Onodera et al. 2015;Estrada-Carpenter et al. 2019). However, we observe an evolution toward smaller Z with increasing redshift as observed in Beverage et al. (2021) for 68 massive quiescent LEGA-C galaxies, or in Carnall et al. (2022) with VANDELS galaxies at z ∼ 1.2. We traced that this effect is due to a degeneracy between metallicity and dust in the fit, because it completely disappears when we remove that parameter from the fit. We notice also that, instead, the differential age measured is very stable, since the Hubble parameter derived in that configuration varies only by 1.1% with respect to our baseline, well below the currently estimated error. We discuss this point in Appendix A.
The stellar masses derived in the analysis correlate with the observed stellar velocity dispersion. This is a well-established result, usually interpreted with the idea that galaxies with a larger gravitational potential well are capable of retaining more gas and therefore forming more stars.
In conclusion, our sample of cosmic chronometers shows ages and star-formation timescales supporting the scenario that they must have formed at early epochs and with very short star formation events quickly exhausting their gas reservoir and then evolved passively.

Binning Parameters
To apply the CC approach, we need to derive from the age(z) relation obtained in our analysis the differential age evolution ∆t in a given redshift bin ∆z. Since this measurement involves the estimation of a derivative, it is typically convenient, in case of noisy data, to increase the S/N of the data by averaging different values, having a consequently more robust estimate of the differential age. This same approach has been adopted in most of the CC studies, see e.g. ; Moresco (2015); Moresco et al. (2016); Borghi et al. (2022b).
Following the previous works by Moresco et al. (2016);Borghi et al. (2022b), we decide to average our data not only as a function of redshift, but also as a function of the velocity dispersion; this last step is in particular important on the one side since it allows us to detect possible trends of the physical parameters as a function of σ (i.e. as a function of the stellar mass), but at the same time because, as highlighted by Thomas et al. (2011) stellar populations of different stellar mass correspond to population formed at different times and over different timescales. Performing an analysis at almost constant velocity dispersion (or stellar mass) ensures the homogeneity of the tracers compared and as a consequence a non-biased determination of the Hubble parameter. For more details, see Moresco et al. (2022).
We note here that before binning and averaging our data, we further excluded 15 objects from our sample since they had a redshift significantly different from the bulk of the population (see Figure 3); we, therefore, imposed a cut 0.6 ≤ z ≤ 1.0, ending up with 335 galaxies.
From now on, all the results will be referred to this sample. Several different choices of binning can be adopted, including the type and number of bins, as well as the type of average statistics adopted. In particular, we can choose to divide the data into N × M bins of redshift and velocity dispersion, to have bins of fixed widths or divided into equi-populated quantiles, and to estimate, within each bin, the averaged quantities with different methods (mean, median, weighted mean). We must consider a trade-off between the benefit of avoiding uneven data distribution by using quantile bins and the benefit of improved population separation by using fixed bins. Additionally, we must use a sufficient number of bins to achieve both good statistics and an optimal sampling of the age-z trend. If the number is too large, the statistics of each bin deteriorate and create large oscillations; on the other hand, if the number is too small, the evolution trend is smoothed out.
We consider median statistics instead of arithmetic (or weighted) average because it is less sensitive to outliers (including badly constrained galaxies, see Section 4.1. In particular, we decided to avoid the use of weighted average because we observed a positive correlation between the estimated ages and their uncertainties, with younger galaxies having smaller uncertainties. Therefore, by using a weighted average, we would have biased our result toward smaller ages. We use the sampling error of the median in the form of σ = MAD/ √ N , where MAD is the median absolute deviation and N is the number of galaxies in that bin. We test all the possible choices for binning, including combinations of quantile and fixed binning for z and σ , combinations by taking reasonable values from {4, 2, 1} and z in {2, 1} for N zbin and N σ bin respectively. We do not use additional bins to ensure that each bin contains a sufficient number of galaxies. We have approximately an average of 40 galaxies per bin in the case of 4×2 bins that significantly decreases to 20 when the H/K selection is applied. After a careful comparison of the various options, we decide to consider as our baseline binning 4 fixed z bins combined with 2 quantile σ bins and median statistics, since it provides the best trade-off between a large enough redshift range probed, a separation in velocity dispersion allowing us to study the mass effect, and the same number of bins to compare our results to Borghi et al. (2022b). We discuss the rationale of this choice in further detail, as well as quantify its impact on the results, in Section 5.2.
In Figure 5, we compare our baseline binned age(z) to the one obtained by Borghi et al. (2022b). As a first point, we observe an offset between the absolute ages estimated in the two methods of about 0.61 ± 0.05 Gyr. This difference can be actually explained and interpreted by taking into account the different SFH adopted in the two analyses. In Borghi et al. (2022b), the theoretical models for the Lick indices where available only for SSP, while in this analysis, we assumed a more complex and realistic SFH, which is also one of the improvements of this analysis with respect to the previous one. The net effect is the bias in age observed.
It is however striking to observe the accuracy with which both derived ages evolve as a function of redshift. Qualitatively, they both follow extremely well the cosmological lines reported as a reference in the figure, demonstrating that despite the difference in the method, the difference in the assumed SFH, the slightly different threshold in σ adopted, the different number of objects, the differential ages agrees extremely well. This will be demonstrated quantitatively in the following section, but it is important to stress here that this is the first time that two different methods to derive differential ages on a common sample of pure massive and passive galaxies (cosmic chronometers) have been performed. The agreement found between the trends of ∆t is, therefore, an additional piece of evidence supporting the robustness of the CC method as a cosmological probe.
Additionally, we observe a distinct mass-downsizing pattern for which more massive galaxies (σ > 208 kms −1 ) exhibit a higher redshift of formation (z f ∼ 2.5) with respect to less massive ones (z f ∼ 2). Our z f estimates are approximately 0.5 higher than those in Borghi et al. (2022b), still due to the use of a more extensive and realistic SFH (DED model), that allows the star formation to start earlier and persist over a longer period of time with respect to the SSP. Our results also achieve higher z and smaller errors of binned ages due to the use of a larger sample of galaxies (335) with respect to Borghi et al. (2022b) (140), which reduces statistical errors, particularly in the third z bin. With the same number of z bins, our result also gives a larger range of redshift, allowing us to probe to Universe up to a slightly higher z.

Measuring H(z)
We compute H(z) by using binned ages and redshifts described in the previous section, To minimize the impact of fluctuations in the data, we do not use consecutive bins to calculate ∆z and ∆t, but a difference approach based on non-adjacent bins, estimating the difference, in a given σ bin, between the i th and the (i+N/2) th point, where N is the number of redshift bins defined. As presented in Section 5.1 we use two types of bin type: fixed and quantiles (i.e. flexible to ensure an equal number of objects in each bin). This strategy requires making an even number of redshift bins to avoid covariance between results caused by the multiple uses of the same data point. This approach allows us also to estimate the difference between points where the expected age evolution is larger than the associated error, making the estimate of ∆t less noisy and more robust. This differential approach of the cosmic chronometer method plays a crucial role in minimizing the rejuvenation in the star-formation history (Moresco et al. 2022).
We perform this evaluation in each σ regime, obtaining (N bin,σ × N bin,z ) /2 Hubble parameter measurements. Finally, these values are combined to get a single and more accurate estimate of H(z) with an inverse-variance weighted average (as also done, e.g. in Moresco et al. 2016;Borghi et al. 2022b).
The different choices of how to bin and select our data, as well as the assumed SFH, are potential sources of systematical uncertainties for the final H(z). In our analysis, we do not account for other systematics introduced by other assumptions of Stellar Population Synthesis (SPS) model, as we discuss in Section 6 (for a detailed treatment, see Moresco et al. 2020Moresco et al. , 2022. In summary, starting from our baseline result, we estimate the impact on the cosmological results by adopting: different choices of binning (Section 5.1), different flagging methods (Section 4.1), or a different SFH assumption that decouples the rising and declining slopes of the SFH (i.e., double power-law or DPL, see Section 3). This will allow us to estimate the systematic errors due to these effects to be associated to our measurement. The results are shown in Figure 7.
First of all, we observe that the configurations contributing the most to a systematic difference in H(z) are those in which all galaxies are averaged together in a redshift bin independently of their σ . The larger shift underlines, even more, the need to perform the analysis of CC carefully selecting the sample in bins of velocity dispersion (or stellar mass), otherwise, the assumption of having a homogeneous sample of chronometers is dropped, and mixing different galaxy populations will exacerbate the progenitor bias (van Dokkum et al. 2000). As shown in Figure 7, we end up smoothing the evolutionary trend, obtaining a higher H(z) and a larger scatter, resulting in larger statistical uncertainties. We also note that by using (equi-populated) quantile bins in z we obtain higher H(z) values with respect to fixed bins. This behavior may be explained by the uneven redshift distribution of our galaxies (see Figure 1). In particular, because there are fewer high-redshift galaxies, the high-redshift quantile bins span a much wider interval, thus flattening the age(z) relation and, ultimately, increasing H(z) and its associated uncertainty. On the contrary, the σ distribution is approximately Gaussian, which makes the result of quantile and fixed σ bins not significantly different.
As for the number of z bins, the results show that a smaller or larger number of bins produce a higher H(z), which is reasonable because, due to the redshift distribution, a smaller number of bins does not allow to correctly map the slope of the age redshift relation, since the larger dz and lower statistics at high z would artificially flatten the median age(z). Besides, taking more bins, the data will be noise dominated.
Even-though flagging posteriors slightly changes the value of H(z), all the related results are compatible with our baseline result as shown in Figure 7.
The various choices of the SFH could further contribute to systematic uncertainties in the measurement of H(z). To address this point, we fit our sample with a different more flexible SFH model, DPL (see section 3), commonly used in other BAGPIPES analyses . To evaluate the difference brought on by a change in the SFH, we solely alter the SFH assumed and keep the other fitting configuration unchanged. In Figure 6, we make a direct comparison between the ages obtained using the two models, showing that the ages estimated with the two models are compatible. To quantitatively assess the impact of choosing a different SFH on our result, we estimate the Hubble parameter with the DPL SFH. With all other configurations of analyses unchanged, we obtain a H(z = 0.80) = 122.0 ± 21.1 km s −1 Mpc −1 using the results by fitting the DPL model, showing a 7.9% difference comparing to our baseline H(z). We take this as an estimation of the systematic uncertainty caused by choice of the SFH model.
To assess the systematic error, we, therefore, consider our baseline result and quantify how much the results are perturbed by three sources of systematics, namely varying the binning scheme, the applied quality flags, and by assuming a different double power law SFH. We estimate the median difference between our baseline H(z) and the measurements obtained from each source, taking In this equation, we need to account for the difference caused by redshift evolution, ∆H model = H model (z base )− H model (z j ), where the choice of the assumed cosmological model negligibly affects in the case of minuscule redshift difference. We compute the total systematic uncertainty by summing each contribution in quadrature and calculating the upper and lower σ syst.tot. separately.
In summary, we obtain a new cosmology-independent measurement of H(z) = 113.1 ± 15.1(stat.) +29.1 −11.3 (syst.) km s −1 Mpc −1 at z = 0.80. Our measurement is consistent with the H(z = 0.75) = 98.8 ± 33.6 km s −1 Mpc −1 that Borghi et al. (2022b) obtained. We notice that the statistical error decreases from 24.8km s −1 Mpc −1 to 15.1 km s −1 Mpc −1 , i.e. by a factor of approximately 0.6, equivalent to the inverse of 335/140 due to the increasing number of cosmic chronometers used for our H(z) measurement. The biggest contribution to this uncertainty is given by the binning scheme, and in particular when H(z) is computed without separat- ing the galaxies into two σ subsamples. By excluding this contribution the upper systematical error decrease to 26.1 km s −1 Mpc −1 . In Figure 8, we compare our final result with all the currently available H(z), finding a noticeable consistency.

CONCLUSIONS
In this paper, we analyze a sample of 350 massive and passive galaxies mostly (95%) at 0.6 z 1.0 extracted from LEGA-C DR2 in Borghi et al. (2022a), deriving their physical properties from a full-spectral fitting analysis. Given that the ultimate goal of this work is apply the CC method with this dataset, it was optimized and thoroughly tested to minimize the possible contamination of the sample by young, star-forming objects. The derived age(z) relation is then used to constrain the differential ages dt, and to provide a new estimate of the Hubble parameter H(z). Our analysis will also allow us to compare for the first time on the same dataset the 4 fi x e d x 1 q u a n t .
6 fi x e d x 1 q u a n t .
2 fi x e d x 2 q u a n t .
6 fi x e d x 2 q u a n t .
4 fi x e d x 2 fi x e d 4 q u a n t . x 2 fi x e d 4 q u a n t . x 2 q u a n t . . Final H(z) value obtained from full spectral fitting of 335 cosmic chronometers in LEGA-C DR2 with statistical (red) and total (black) uncertainty. The violet point is the value obtained by Borghi et al. (2022b) via Lick indices analysis of 140 galaxies of this dataset. The gray points are all the other OHDs available from the literature (Jimenez et al. 2003;Simon et al. 2005;Stern et al. 2010;Zhang et al. 2014;Moresco 2015;Moresco et al. 2016;Ratsimbazafy et al. 2017). For illustrative purposes only, we include the H(z) prediction assuming a ΛCDM model with Planck Collaboration et al. (2020) parameters. differential ages dt derived with two different and independent methods, namely Lick indices and full-spectral fitting, and to investigate and validate the capability of the CC approach to robustly derive H(z).
Here, we utilize the public code Bagpipes  to derive stellar ages, metallicities, and SFH fitting both the spectroscopic data alone and the spectroscopic and photometric data jointly for all the individual galaxies in our sample. Our main results are summarized as follows.
• We first extend Bagpipes by removing the cosmological prior on the derived ages, to avoid possible biases on our cosmological results due to the assumption of a fiducial cosmology. We also adopt flat uninformative priors on all the derived quantities, namely the stellar age, metallicity, and SFH. We then explore the dependence of our results on the SFH assumed, the parameters included in the fit, and the priors considered. We opt to use a delayed exponentially declined SFH, improving with respect to the analysis of Borghi et al. (2022a) where SSP where assumed, but minimizing the number of free parameters for the functional form of the SFH, since we verified that the gain in the quality of the fit was marginal with other choices.
• We find the results obtained from the fit to the spectroscopy alone to be less accurate and more scattered than the results obtained by fitting the combination of spectroscopy and photometry, with in particular larger tails toward higher ages and τ . We defined a set of indicators of quality of the fit and convergence criteria, based on the inspection of the posterior distribution and on the accuracy with which the best-fit was reproducing specific observational features in the spectrum, namely the CaII H/K known to be correlated with potential episodes of recent star formation. We demonstrate that the nonphysical scatter in the derived parameter obtained using only spectroscopic data can be lifted by applying masks defined on these indicators.
• We observe that the inclusion of photometric data (21 bands in this analysis) allows the fit to converge correctly even without applying the convergence criteria previously discussed, reducing the degeneracy between parameters. As a consequence, in this framework we are able to maximize the final number of objects with a correct fit, and we decide to consider this as the baseline of our analysis.
• We find that the measured age(z) relation are well compatible with a cosmological aging as a function of redshift even with assumed a flat prior on age∈ U(0, 20) Gyr. Our results present also a clear downsizing trend when divided into two bins of velocity dispersion, with galaxies with σ > 108 km s −1 with a formation redshift z f ∼ 2.5 and the ones with σ < 108 km s −1 with a formation redshift z f ∼ 2. Even though we consider a significantly wider prior in our analysis, these galaxies show very short star formation timescales with a median value of τ = 0.36 ± 0.13, • The average measured stellar metallicity is Z/Z = 0.84 ± 0.41 with a small hint of evolution as a function of redshift. We prove that this tension or fitted evolution is due to the degeneracy between the metallicity and dust in the fit, and demonstrate that it has an almost negligible effect on the differential age.
• We compare for the first time the stellar ages derived from two very different and independent methods applied on the same sample, the fullspectrum fitting and the Lick indices analysis. We find that the absolute ages derived present an offset of 0.61 ± 0.05 Gyr, which can be understood by recalling that the Lick index models assume an idealized SSP SFH that slightly bias the absolute values toward younger ages. The agreement, however, between the differential ages is striking, and it is important to underline here that the CC approach is based on the measurement of dz/dt, that perfectly agrees within errors in the two measurements. This result is in particular important because it demonstrates the robustness of the method and the stability of the dt measurement, confirming that it can be derived with significantly less biases than absolute ages.
• From the analysis of the binned age(z) relation, we derive a new H(z) measurement H(z = 0.80) = 113.1 ± 15.1(stat.). We verify that our result is fully compatible with the one by Borghi et al. (2022b), even if at a slightly larger redshift and with a slightly smaller statistical error due to the different number of objects in the final sample used by the two methods (in Borghi et al. 2022b, to ensure homogeneity in the analysis we decided to consider in the final sample only the spectra for which the same number of spectral features where observable). we also test that our result is consistent with other literature OHDs, and as well with the prediction of ΛCDM model assuming the Planck Collaboration et al. (2020) cosmological parameters.
• We assess the systematics involved in the results by varying the methods with which the binned age(z) relations are derived, changing the number of bins, the method to estimate the average value, the assumed SFH model in the fit, and testing the application of the masks described above. We estimate a systematic error of +29.1 −11.3(syst.) km s −1 Mpc −1 , mainly dominated by the large variation in the results obtained when the sample is not divided into two σ (or stellar mass) bins, suggesting that an analysis in specific ranges of masses is fundamental to ensure the homogeneity in formation of the CCs considered.
We underline that since the sample used in this analysis and in the one of Borghi et al. (2022b) are drawn from the same parent sample, it should be avoided to use them in combination, since the measurements will be highly covariant. We also underline that in the current analysis we decided not to explore the further dependence of our result on other assumptions of the SPS models. Regarding the dependence on different SFH, we verified that within Bagpipes the SFH choice is somehow limited (we could have chosen amongst SSP, DED, and DPL, but the SSP yields discrete pattern on galaxies' ages, weakening the reliability of the differential ages obtained, which however can be improved by a re-run of the SPS on more refined grids.), therefore it would not have allowed a full estimate of this effect. Moreover, the fit obtained with the DPL SFH provided similar results to the ones obtained here, but with a higher number of free parameters, not justifying in our case the choice of that SFH. In the end, we chose the DED model as our baseline, and the DPL for evaluating the SFH choice caused systematic uncertainty.
Further investigation is crucial to go beyond analytic forms of SFHs. We notice that non-parametric SFHs (e.g., Leja et al. 2019) could be more flexible options than these analytic approximations in describing the full diversity of SFH shapes, and they are becoming the recently popular alternatives to parametric ones. Therefore, we acknowledge that to assess systematic effects they could be interesting alternatives to exploit. We explored this possibility using the latest update of Bagpipes, which includes the possibility of using a nonparametric SFH proposed by (Leja et al. 2019), considering a model with a continuity prior on ∆ log(SFR) between adjacent time bins (7 bins in total). The results we obtained are extremely encouraging, and point toward the fact that the SFH uncertainty is currently not dominant in our analysis. We found that the Hubble parameter estimated considering non-parametric SFH is compatible within 0.27σ with respect to the one we derive in our analysis. A larger statistical uncertainty is also expected in this case, since a more flexible SFH with a larger number of free parameters is considered. However, further checks and verifications that go beyond the scope of this paper are needed before including this result as further systematic uncertainty in our analysis, and we defer it to a following paper.
Dry mergers of less massive galaxies hosting younger stellar populations may bias the light-weighted age estimation toward younger ages. However, this effect is expected to be sub-dominant since this analysis is based on differential ages and spans a limited redshift interval ). While the sample selection criteria in section 2 are aimed to minimize contamination of the young component, a residual, even if minor, evolution could bias our measurement, in Moresco et al. (2018) a recipe was provided to include in the covariance matrix of the H(z) measurement an error contribution due to this effect, that we quantify to be negligible in this sample.
The rejuvenation in the star-formation history could introduce bias in the measurement of the differential ages. This has been addressed in detail in ; Moresco et al. (2018), where it has been quantified the impact of these effects on the Hubble parameter (see in particular App. A.1 of , and Sect. 5 and 6 of Moresco et al. (2018)). The main point to stress here is that the selection process in section 2 and the differential approach of the cosmic chronometer method in section 5 play a crucial role in minimizing this effect. Moreover, the systematic uncertainty originating from the binning method takes into account the potential progenitor bias, as we explain in section 5.2.
To include in the current measurement a proper full systematic covariance matrix, we therefore suggest the reader to follow the procedure described in details in Moresco et al. (2022) and Moresco et al. (2020), including the missing statistical effects not already included in this analysis.
In conclusion, this work provides a further important piece of evidence supporting the robustness of the CC method as a cosmological probe, showing the potential of the full-spectral fitting approach as another different method to derive the relative ages of massive and passive galaxies. It is interesting to notice here that our sample is almost a factor 2 larger than the final sample used by Borghi et al. (2022b), and the statistical error between the two measurement scales as expected as √ N . This is very promising in the view of several current (SDSS BOSS Data Release 16, Ahumada et al. 2020) and incoming spectroscopic surveys (such as Euclid Laureijs et al. 2011 andWang et al. 2019), that will significantly improve the census of massive and passive galaxies, especially at z > 1. . Median binned metallicity-redshift relations for our baseline results based on a full-spectrum fitting that contains a dust model (bottom, Carnall et al. 2018) and for the fitting under the same configuration, but with the dust model switched off (upper). The blue and red points represent the lower and higher σ bins, divided using the median value of each sample as the threshold. For each bin, the vertical error bars are the errors associated to the median ages, while horizontal bars denote the bin width.