© The Authors 2025

1. Introduction

The ΛCDM model is regarded as the simplest theoretical framework, providing a remarkable fit to the bulk of cosmological observables, namely the cosmic microwave background (CMB) temperature and polarisation spectra (Planck Collaboration VI 2020; Aiola et al. 2020; Dutcher et al. 2021), baryon acoustic oscillation (BAO) measurements (Eisenstein et al. 2005; Beutler et al. 2011; Ross et al. 2015; Zhao et al. 2022; Adame et al. 2025), the distance-redshift relation with type Ia supernovae (Perlmutter et al. 1999; Scolnic et al. 2018; Brout et al. 2022) and weak gravitational lensing (Abbott et al. 2020, 2022; Dvornik et al. 2023; Dark Energy Survey and Kilo-Degree Survey Collaboration 2023), among others. However, the era of precision cosmology has revealed a number of observational inconsistencies that could point to a fundamental defect in the model, the most famous of which is the Hubble tension (see Di Valentino et al. 2021, for a review). Indeed, ever since the first release of the Planck cosmological analysis of the CMB (Planck Collaboration XVI 2014), a pervasive discrepancy between early- and late-time measurements of the cosmological parameter, H₀, has plagued the scientific community. The former, represented by CMB, BAO, and Big Bang nucleosynthesis (BBN) analyses, assume a Lambda cold dark matter (ΛCDM) model and agree upon low values of the Hubble constant (e.g. H₀ = 67.27 ± 0.60 km s⁻¹ Mpc⁻¹ as reported by Planck Collaboration VI 2020). The latter, characterised by ‘direct’ model-independent approaches, are represented (although not limited to) the cosmic distance ladder approach and systematically predict high values, such as the last determination of H₀ = 73.04 ± 1.05 km s⁻¹ Mpc⁻¹ by Riess et al. (2022).

In this troubled arena, the importance of emergent, independent and complementary cosmological probes should not be overlooked (Moresco et al. 2022), especially if they require no assumptions on the underlying cosmological model. This is the case of cosmic chronometers, which have come forth as a promising alternative to probe the expansion history of the Universe. In general terms, cosmic chronometers are astrophysical objects that can be statistically observed and dated at different redshifts, thus allowing us to reconstruct the time evolution of a large-scale homogeneous and isotropic Universe without assuming a specific cosmological model. This idea was first put forward by Jimenez & Loeb (2002), who suggested the use of passively evolving galaxies as cosmic chronometers, arguing that identifying the shift between the age distribution of two such populations at different epochs would yield an estimate of the Hubble parameter at an effective redshift between them. Although this exact approach was not used, the subsequent study by Simon et al. (2005) put the underlying idea to the test and provided measurements of the Hubble parameter at different redshifts using data from the Gemini Deep Deep Survey (GDDS; Abraham et al. 2004) via the single age differences of galaxies determined to be in passive evolution after a very short burst of star formation. Indeed, within the downsizing scenario (Cowie et al. 1996), the most massive elliptical galaxies are older, formed earlier and in a shorter burst, after which they have been in passive evolution (Kauffmann et al. 2003; Gallazzi et al. 2005; Thomas et al. 2005, 2010; Conroy et al. 2014). In other words, at a given redshift, the stellar populations of passive elliptical galaxies with the same stellar mass should be increasingly coeval, as one examines progressively more massive or α-enhanced galaxies. As a consequence, the cosmic chronometer approach relies on an accurate identification of massive and passively evolving elliptical galaxies at any given redshift (Borghi et al. 2022a).

The other essential aspect is the determination of the stellar age of galaxies. Although the cosmic chronometer approach is differential in nature (in the sense that the relevant quantity for the analysis is not the absolute ages but their redshift evolution), a robust dating of the stellar populations is still needed. Two main methods have been used to derive stellar ages from passive galaxy spectra, the first of which is full spectral fitting (FSF). Full spectral fitting techniques (Heavens et al. 2000; Cid Fernandes et al. 2005; Tojeiro et al. 2007; Chevallard & Charlot 2016; Wilkinson et al. 2017; Carnall et al. 2018) involve deriving the physical properties of a galaxy by fitting the entire observed spectrum with synthetic stellar population models, thus considering all available spectral features as well as the continuum. This approach has recently been applied in the cosmic chronometer literature (Simon et al. 2005; Zhang et al. 2014; Ratsimbazafy et al. 2017; Jiao et al. 2022; Tomasetti et al. 2023), producing 16 distinct data points in the Hubble diagram. The second method involves the modelling of specific spectral features that can correlate with age, metallicity, and star formation history. One example of this is the characteristic D4000_n break in passive galaxy spectra, which has been primarily employed by Moresco et al. (2011, 2012, 2016) owing to its simple tight correlation to stellar age at a fixed metallicity that can be calibrated using stellar population synthesis models. A total of 15 measurements of H(z) have so far been obtained for z < 2 using this method, with overall tighter constraints with respect to FSF approaches. However, several other spectral features have been traditionally employed to determine the age of passive elliptical galaxies given their sensitivity to stellar population properties.

Lick indices (Worthey et al. 1994; Worthey & Ottaviani 1997) are a set of 25 well-studied absorption features in the spectra of elliptical galaxies. They are defined by relatively wide bandpasses and two additional windows characterising pseudocontinua blue- and redwards of the central region. Their measurement in early-type galaxy (ETG) spectra signalled the presence of supersolar Mg/Fe ratios (Worthey et al. 1992; Davies et al. 1993; Carollo & Danziger 1994; Mehlert et al. 1998; Jørgensen 1999; Longhetti et al. 2000), which strongly suggested short star formation timescales in agreement with the downsizing scenario. This led to the development of stellar population models that predicted the strength of Lick absorption features as a function of single-burst stellar age, metallicity, and varying element abundances (Thomas et al. 2003, 2004, 2011; Tantalo & Chiosi 2004; Schiavon 2007; Chung et al. 2013). In addition, the alternative avenue of using models that predict the full spectral energy distribution (SED) at a moderate-high resolution (Schiavon et al. 2002; Bruzual & Charlot 2003; Vazdekis et al. 2010, 2015; Conroy & van Dokkum 2012; Conroy et al. 2018; Knowles et al. 2023) can also be exploited to predict Lick index strengths and derive stellar population parameters. However, the Lick index method has only been employed by Borghi et al. (2022b) in the cosmic chronometer context, providing a single measurement of the Hubble parameter at z ∼ 0.7.

The aim of this work is to provide independent measurements of the Hubble parameter at low redshifts (z ≲ 0.4) via the cosmic chronometer technique. These measurements are alternative to those based on FSF (Simon et al. 2005; Zhang et al. 2014) and the D4000_n break (Moresco et al. 2012), not only because of the method of estimating the stellar ages but also due to the nature of the approach. Indeed, instead of using finite differences to estimate the redshift derivative of cosmic time, we trace the continuous redshift evolution of the ages via a cosmographic model-independent approach. The cosmographic approach (Capozziello et al. 2011) consists in directly fitting H(z) (or derivative and/or integral functions) when H(z) is expressed as a Taylor series on kinematical parameters (deceleration q₀, jerk j₀, snap s₀). The expansion on the y−redshift holds well up to y ∼ 0.4 (Capozziello et al. 2011), which translates to z ∼ 0.7.

We provide a thorough description of the careful selection of the cosmic chronometer sample and employ two different stellar population synthesis (SPS) models, namely Thomas et al. (2011) and Knowles et al. (2023), to determine the age of the stellar populations from Lick index strength predictions. It should be noted that some models allow for a variety of internal configurations, and different models can significantly alter the final readout of H(z). For example, one can find in the literature the work by Moresco et al. (2020), which analyses the impact of the chosen initial mass function (IMF), stellar library, and overall SPS model on the final H(z) estimation. The results from that study highlighted that the global SPS model and the stellar library are the main drivers of disagreement on the final H readout, while the IMF played only a minor role. A comprehensive analysis of the various contributions to the systematic uncertainty budget (Moresco et al. 2022) is necessary in order to predict the current and future power of cosmic chronometers in providing a reliable estimation of the Hubble parameter. However, this work aims to produce a robust statistical measurement of H(z) through the ages derived with the aforementioned models and by accounting only for the systematic uncertainty contribution driven by data management.

The paper is structured as follows: Section 2 introduces the cosmic chronometer framework, outlines the cosmographic approach to the Hubble parameter, and provides a detailed description of the selection criteria for the cosmic chronometer sample. Section 3 details the methodology involved in the process of measuring and fitting Lick absorption indices to the models, including a careful description of different aspects, namely spectrum stacking, instrumental resolution and velocity dispersion corrections. The fit of the stellar ages to sample the Hubble parameter within the cosmographic framework is also explained. Section 4 presents the outcome of the cosmographic fit for both the TMJ and Knowles modelled ages. With the former, we present the estimation of the Hubble parameter as our central result, accompanied by a brief discussion of the sources of systematic uncertainty. Section 5 is included to discuss other implications of our results beyond the cosmographic analysis. We leveraged the low-redshift data to obtain an estimate of the scaling relations of the ETGs using the two SPS models. The final part of this section covers the analysis of the oscillations in the t − z relation and the role that index strengths play. Finally, Section 6 summarises the main conclusions.

2. Cosmic chronometers

In this section, we briefly introduce the theoretical framework of cosmic chronometers, focusing on the cosmographic approach that we follow throughout the paper. Then, the cosmic chronometer sample will be presented, with a detailed description of the selection criteria and its physical properties.

2.1. Theoretical framework

The cosmic chronometer approach typically relies on approximating the differential expression of the Hubble parameter given by

$$\begin{array}{c}\hfill H\equiv \frac{1}{a}\frac{\text{d}a}{\text{d}{t}_{\text{U}}}=-\frac{1}{1+z}\frac{\text{d}z}{\text{d}{t}_{\text{U}}},\end{array}$$(1)

by the ratio of the intervals in redshift (Δz) and age¹ (Δt) between two well-defined samples of cosmic chronometers. These two groups statistically trace the same population at different redshifts (Jimenez & Loeb 2002), i.e. (t₁, z₁) and (t₂, z₂). Indeed, since the average formation time of cosmic chronometers observed at different redshifts is defined to be equal, it factors out when computing the difference in cosmic time between two redshifts, that is, Δt_U = t_U(z₂)−t_U(z₁) = t₂ + t_U(z_form)−(t₁ + t_U(z_form)) = Δt. In other words, the time elapsed between two redshifts and the difference in the age of the stellar populations of cosmic chronometers observed at those redshifts are equivalent. This gives a local measurement of H at the mean redshift ⟨z⟩≡(z₁+z₂)/2 of the ensemble. Therefore,

$$\begin{array}{c}\hfill H(⟨z⟩)\approx -\frac{1}{1+⟨z⟩}\frac{\mathrm{\Delta }z}{\mathrm{\Delta }t}·\end{array}$$(2)

However, we decided to take a different path from that of Equation (2) and follow a cosmographic approach, valid for relatively low redshifts (Capozziello et al. 2011). As observed when the full results are presented in Section 4, the overall distribution of ages for each coeval population can feature fluctuations that would critically affect (i.e. reading H < 0) the point estimation of H obtained following (2). This issue is not new; indeed Tomasetti et al. (2023) found an anomaly at higher redshift by which they would measure a negative value for H if used (Appendix A: The z < 1.07 anomaly). Since we aim at drawing an estimation for H based on a wealth of data from a rich catalogue, it could be argued that these anomalies could be addressed statistically, with H < 0 measurements being outliers. Despite this, there are two more caveats that make us prefer the cosmographic analysis based on a functional fit. The first is the manual pairing of groups of galaxies (or stacks in our case), which implies a further arbitrary choice that shall be accompanied with a sensible analysis of its effect as a systematic uncertainty (Borghi et al. 2022b). The other is the fact that a continuous fit of t(z) allows us not to reduce the redshift window for which our estimation is solid. In other words, a point measurement of H, based on a pair of points, takes as fiducial redshift the average redshift of the pair. Instead, following our methodology, each point individually contributes with its own spectroscopic redshift to the estimation and so our original range of redshift is kept.

Working directly with t(z) requires integrating the inverse of the Hubble function. By doing so, one gets the evolution of cosmic time as a function of redshift,

$$\begin{array}{c}\hfill t_{\text{U}}(z)=t_{\text{U}}(0)-{\displaystyle {\int }_0^z}\frac{1}{1+z\prime }\frac{\text{d}z\prime }{H(z\prime )},\end{array}$$(3)

where t_U(0) is the current age of the Universe. In terms of the stellar ages of cosmic chronometers, we can rewrite (3) as

$$\begin{array}{c}\hfill t(z)=t_0-{\displaystyle {\int }_0^z}\frac{1}{1+z\prime }\frac{\text{d}z\prime }{H(z\prime )},\end{array}$$(4)

where t₀ ≡ t_U(0)−t_U(z_form) is the age of the cosmic chronometer population at z = 0. Therefore, the continuous redshift evolution of the ages can be exploited to derive the expansion history of the Universe via the parametrisation of the Hubble parameter. Following Capozziello et al. (2011), one can perform a series expansion that overcomes convergence issues at high redshift, with the advantage of being independent of the underlying cosmological model. Under the introduction of the y-redshift, defined as y ≡ z/(1 + z), the expansion to order k of H(y) reads

$$\begin{array}{c}\hfill H(y)=H_0·{\displaystyle \sum _{i=0}^k}({\mathcal{H}}_0^i·y^i)+\mathcal{O}(y^{k+1}),\end{array}$$(5)

where the terms written as ℋ₀ⁱ are dimensionless as the Hubble constant has been factored out so that the sum reflects just the kinematics of the expansion of the Universe. Given the redshifts of interest in this work, it is enough to use the expansion up to second order (k = 2). The relevant coefficients are ℋ₀⁰ ≡ 1, ℋ₀¹ ≡ ℋ₀¹(q₀) and ℋ₀² ≡ ℋ₀²(q₀, j₀), where q₀ and j₀ are the so-called deceleration and jerk parameters. With all this, we can rewrite the cosmographic expansion of H(z) as

$$\begin{array}{c}\hfill \frac{H(z)}{H_0}\approx 1+(1+q_0)\frac{z}{1+z}+(1+\frac{j_0}{2}+q_0-\frac{q_0^2}{2}){\left(\frac{z}{1+z}\right)}^2.\end{array}$$(6)

Introducing the above equation into (4) yields the cosmographic model for the redshift evolution of the stellar ages of cosmic chronometers. In particular, the ΛCDM model is recovered by setting q₀ = Ω_m/2 − Ω_Λ and j₀ = 1.

2.2. Cosmic chronometer sample

As described in the introduction, we aimed at selecting the most massive and passively evolving sample of ETGs. This is the group of red and metal-rich elliptical galaxies that have not undergone any major process of star formation since the end of the formation of the bulk of their stellar mass (Carretero et al. 2007; Sargent et al. 2015). This occurred at z ≳ 1 (Yildirim et al. 2017) and in a short burst (Cimatti et al. 2006), implying α-enhanced stellar populations. It is well established in the literature that through a cross-selection in morphology (Clemens et al. 2006; Shimasaku et al. 2001), colour (Moresco et al. 2016), emission lines (Moresco et al. 2011) and the ratio between the CaIIK and CaIIH absorption lines (Borghi et al. 2022a), a clean sample of these objects can be obtained.

Since this work targets low (z ≲ 0.4) redshift, the sources are extracted from the Sloan Digital Sky Survey (SDSS), in particular from the SDSS Legacy Survey corresponding to SDSS-I/II (York et al. 2000; Abazajian et al. 2009). In a first step, we submitted a query to the SDSS catalogue selecting all galaxies satisfying a series of criteria including morphology, velocity dispersion (a proxy for stellar mass) and other quality flags for their spectra. This process is performed via a SQL query and is described as follows:

1. Objects identified as a galaxy with no star formation activity, using the columns Class and subClass from the specObj view.

2. Velocity dispersion restricted to lie between 100 and 400 km s⁻¹, in order to select massive galaxies with stellar masses greater than ∼10^10.5 M_⊙, preventing contamination from rotation-dominated objects (Veale et al. 2017) that are typically star-forming (Di Teodoro et al. 2016).

3. Elliptical morphology, imposed through the restriction of the concentration or compactness parameter as indicated in Shimasaku et al. (2001). In other words, the ratio between the radii containing 50% and 90% of the Petrosian flux has to be less than 0.33.

4. Rest-frame wavelength window to measure spectral features. Indeed, requiring that all 25 Lick indices (and the CaII H and CaII K lines for reasons that are explained below) be observed would imply that the spectra should be available in rest-frame wavelengths between 3600 and 6500 Å. Although this would be a conservative cut, as the reddest indices (NaD, TiO₁ and TiO₂) are not used in the analysis, the window is actually extended to 6700 Å to be able to later exclude objects with emission in Balmer Hα (6526 Å). Since the SDSS observed window reaches up to ≈10 000 Å, this imposes the restriction that the redshift be below z ∼ 0.5. This is not a problem for the scope of this work, but it is important to keep in mind should one wish to extend the analysis to higher redshifts using SDSS data. In that case, one could drop the upper limit on the rest-frame window down to some value around 5100 Å (effectively including [O III] λ5007), or slightly higher if one wishes to measure the Fe5270 and Fe5335 indices, used to build the metallicity-sensitive [MgFe] and [MgFe′]. For reference, z ∼ 0.9 can be achieved by only dropping the requirement on the absence of emission in Hα.

We began from a parent sample, constructed by requiring that the objects be galaxies and satisfy the wavelength coverage condition, which amounts to a total of 669475 galaxies. This number was reduced to 506475 once the velocity dispersion constraint was imposed; then, both the star-forming flag and the concentration parameter conditions further diminished the number to 135554 galaxies, which constitute what we call the SQL query sample. It should be noted that the most restrictive criterion is the concentration parameter, which reduces the number of sources by two thirds.

Once the SQL query sample was defined, we proceeded to download the spectra to perform a further spectroscopic selection and remove objects that could show some remaining star-formation. To this end, we inspected emission lines associated with young stellar systems and ongoing star formation, namely [OII]λ3727, [OIII]λ5007, Hβ and Hα (Sobral et al. 2015; Suzuki et al. 2016; Khostovan et al. 2020). The SDSS provides the measurements of the equivalent width (EW) of emission lines as performed by Tremonti et al. (2004). For completeness, we decided to also measure the EW ourselves to cross-check with the SDSS values. We found good agreement and, for the purposes of this work, we decided to make a combined (stricter) cut requiring galaxies to meet the selection criteria for both measurements. In particular, we imposed that any of the EWs be larger than −2 Å, although values between −5 and −2 Å were allowed if the S/N of the line was smaller than 3, meaning that some noise can be mistakenly identified as an emission. For reference, previous studies with cosmic chronometers have agreed to exclude any galaxy with some of these lines showing emissions stronger than −5 Å (Moresco et al. 2011, 2012; Tomasetti et al. 2023; Borghi et al. 2022a). The criterion on the S/N varies slightly instead; for example, Moresco et al. (2016) chose a maximum threshold of 2 in S/N, while Tomasetti et al. (2023) relaxed it to a value of 3. Borghi et al. (2022a) did not impose any restrictions, although they found that the final sample rarely showed a S/N higher than 3.

To complement this criterion, we decided to put in practice one of the results from Borghi et al. (2022a) by making a cut in H/K, the ratio between the CaII H and the CaII K pseudo-Lick indices. Larger values of this quantity indicate the presence of star formation (Cincunegui et al. 2007; Moresco et al. 2018; Pimbblet et al. 2019; Werle et al. 2020) as CaII H deepens due to its blending with the Balmer Hε from young hot stars, while CaII K remains unaffected. Following Borghi et al. (2022a), we select sources that satisfy H/K < 1.2.

The final sample is made up of 90396 sources and extends up to z ∼ 0.38. Fig. 1 shows how the distribution of some absorption features (D4000_n and three Balmer indices) changes from the parent to the final sample. The typical bimodality separating early- from late-type galaxies is clearly observed in the D4000_n, Hδ_A and Hγ_A distributions of the parent sample (in grey), which is already broken with the query sample (in blue), clearly defining a much purer passive ETG sample, which is culminated in the final step after adding spectrometry (maroon). Note how the distributions for D4000_n and Hδ_A are equivalent to the ones observed by Borghi et al. (2022a).

Fig. 1. Distribution of D4000n and three Balmer-line indices over the parent (grey), query (blue), and final (maroon) galaxy samples.

Lastly, studies with ETGs typically rely on a photometric selection, including cuts in galaxy colour that separate between early and late types. These criteria have been known for a long time, as Strateva et al. (2001) already showed that the distribution of SDSS galaxies on the u − g/g − r plane presented a strong bimodality and proposed u − r = 2.22 as an optimal separator. More recently, cuts in colour have been applied for the selection of cosmic chronometer candidates. For instance, Moresco et al. (2016) selected early-type galaxies from BOSS (Baryon Oscillation Spectroscopic Survey) from the g − r/r − i plane following Eisenstein et al. (2011), Masters et al. (2011), White et al. (2011), and Maraston et al. (2013), while Borghi et al. (2022a) used (NUV − r)₀/(r − J)₀ (Ilbert et al. 2013) and Tomasetti et al. (2023) used (U − V)₀/(V − J)₀ colour² planes (McLure et al. 2018) to identify passive galaxies in the LEGA-C (van der Wel et al. 2016) and VANDALS (McLure et al. 2018) surveys, respectively.

Since criteria based on colours typically have non-negligible percentages of contamination or incompleteness, we did not apply a further cut based on colours. However, we did want to make sure that the final sample was physically reasonable in terms of its colour distribution. For the SDSS photometric bands at our disposal, Holden et al. (2012) defined a passive galaxy region on the (u − r)₀/(r − z)₀ plane by the polygon

$$\begin{array}{c}\hfill \{\begin{array}{cc}\hfill {(\text{u}-\text{r})}_0& >2.26\hfill \\ \hfill {(\text{r}-\text{z})}_0& <0.75\hfill \\ \hfill {(\text{u}-\text{r})}_0& >0.76+2.5{(\text{r}-\text{z})}_0,\hfill \end{array}\end{array}$$(7)

which is plotted with a dashed black line in Fig. 2. We also depict, for the 0.05 < z < 0.06 redshift bin, the rest-frame colours of all galaxies in the parent sample (grey) and all the subsequent subsamples, namely adding a velocity dispersion cut (cyan), the star-forming flag and concentration parameter (blue) and the spectroscopic selection (maroon) defining the final sample. K-corrections have been applied to the observed magnitudes following Chilingarian et al. (2010) and O’Mill et al. (2011). The behaviour is qualitatively similar for all redshifts in the sample, although with different incompleteness fractions. In all cases, these numbers never exceed the expected 18% of missing passive galaxies (Holden et al. 2012). We can safely conclude that the selection procedure is physically solid.

Fig. 2. Distribution of sources in the (u − r)0 vs. (r − z)0 colour-colour plane. The parent sample is shown in grey, while the increasingly restrictive subsamples (adding velocity dispersion, query and spectroscopic information) are shown in light blue, dark blue and maroon, respectively. The dashed black line defines the passive galaxy region (Holden et al. 2012).

3. Data treatment and methodology

The cosmic chronometer framework involves deriving stellar ages from the spectra of the final galaxy sample detailed in the previous section. To this end, Lick absorption features were measured and fitted to the predictions of stellar population models for Lick index strengths. However, the process is not straightforward and needs to be carefully explained. In particular, high S/N spectra are needed in order for the index measurements to be reliable. This is accomplished by stacking the flux densities of single galaxies of similar nature (Ferreras et al. 2019; Sánchez et al. 2020; Parikh et al. 2021). Once these ‘raw’ measurements are performed, corrections are applied to account for the astrophysical and observational effects that are tangled in the spectra and affect the real value of the indices, namely instrumental resolution and velocity dispersion. Indices can then be fit with two flux-calibrated models, the Lick-index model of Thomas et al. (2011) (referred to as TMJ for the rest of the paper) and the single stellar population (SSP) model based on the sMILES library Knowles et al. (2023) (referred to as Knowles) via a Markov chain Monte Carlo (MCMC) algorithm. An estimation of the luminosity-weighted ages, metallicities and α-enhancement of the populations is obtained. Finally, the age-redshift relations are fitted within a cosmographic approach to derive a model-independent sampling of the Hubble parameter up to z ∼ 0.4.

3.1. Stacking

Lick index measurements demand a high enough S/N of the spectra, which is not guaranteed by individual galaxies and can compromise the reliability of the estimation of age, metallicity, and α-enhancement. Single galaxy spectra are thus shifted to the rest frame, normalised, and then stacked. In this process, we followed Borghi et al. (2022a) to derive the composite spectra and their associated uncertainties. We chose a central stable window (rest frame) between 4200 and 4400 Å and normalised all spectra as $\stackrel{\sim }{F}(\lambda )\equiv F(\lambda )/⟨F([4200,4400])⟩$, where ⟨…⟩ denotes the median. Since the SDSS spectra are not sampled on the same wavelength array, they need to be resampled into a common grid and interpolated, for which we set a wavelength step of 0.1 Å. The flux of the stacked spectrum is computed through the median of the individual normalised spectra at every wavelength and the uncertainties are given by the normalised median absolute deviation as $\delta F_{\text{stack}}(\lambda )=1.4826·⟨F_{i\in \text{stack}}(\lambda )-⟨F_{i\in \text{stack}}⟩(\lambda )⟩/\sqrt{n}$, where n is the number of objects in the stack and F_i is each individual galaxy spectrum (Hoaglin et al. 1983).

To determine which sources to stack together, we required that the final (stacked) spectra all have a similar S/N, a quantity that clearly depends on the number of sources included in each stack and their redshift. In other words, we wished to estimate the function N(z, ⟨S/N⟩) that would tells us the number of sources N needed at a given redshift z to obtain a median ⟨S/N⟩. For this reason, we built an auxiliary set of stacks through a 2D grid in velocity dispersion and redshift. We set step sizes of 20 km s⁻¹ (from 100 to 400 km/s) and 0.01 (from 0.05 to 0.40), respectively. This is shown in the left panel of Fig. 3, where each pixel (auxiliary stack) is coloured according to the number of galaxies it contains. To sample N(z, ⟨S/N⟩). For each auxiliary stack, we measured its S/N and kept track of the number of single galaxies included within. The result is plotted in the right panel of Fig. 3, where each point corresponds to a different pixel. As expected, stacks need more galaxies at a higher redshift to attain a given ⟨S/N⟩ and, at high redshift, it is uncommon to observe pixels with ≳320 sources.

Fig. 3. Properties of the ancillary stacks. Left panel: Distribution of the number of galaxies included in the ancillary stacks defined by a linear constant spacing in both σ and z. Right panel: Relation between the S/N and the number of galaxies in each ancillary stack as a function of redshift.

Finally, we performed a linear fit of log(⟨S/N⟩) as a function of log(N) in bins of z to sample N(z) at a fixed ⟨S/N⟩, which we set to 300 as the reference value. At this point we obtained a set of values of log(N) that correspond to a set of values of z, with ⟨S/N⟩ fixed to 300. The final N(z) function is an interpolation of this set of data. The final stacks were then constructed taking into account the effect of velocity dispersion, for which we defined six groups logarithmically spaced from 10² to 10^2.6 km s⁻¹. This choice roughly corresponds to setting the bin limits at [100, 125, 160, 200, 250, 320, 400] km s⁻¹. Then, the procedure to assemble the stacks in each group is the following:

1. We ordered galaxies in ascending redshift.

2. We ‘opened’ a new stack and included, one by one, galaxies in order of redshift (loop explained in steps 3–6).

3. After each inclusion, we computed the average redshift of the stack, ⟨z⟩_loc, and keep track of the number of galaxies included, n_loc.

4. We computed, using the N(z) function we built from the auxiliary stacks, the number of sources that will be necessary to obtain ⟨S/N⟩∼300 at redshift ⟨z⟩_loc. This is, N_ref ≡ N(⟨z⟩_loc).

5. When N_ref > n_loc, we continued the loop with the following galaxy (step 3). When N_ref < n_loc we ‘closed’ the stack (step 6).

6. When the stack was ‘closed’, we proceeded to measure the spectral indices. Then, a new stack was ‘opened’, starting with the very next galaxy in redshift order. We repeated this methodology up to the last galaxy in the sequence.

The final stacks can be found in Figs. 6 and 8, distributed in the plane t − z according to the ages we estimated using the SPS models that we present in Section 3.4. Unless explicitly stated otherwise, any distribution of ages presented in the remainder of the article is computed on stacks built trying to match a median S/N of 300.

From the final stacks, all those that included an internal distribution of sources that exceeds in width an average dispersion of Δz = 0.01 (σ_z = 0.005) were discarded, in order to guarantee that in each stack only the most similar individual galaxies would be included. We observe that some of those data points that cover vast redshift windows end up showing a behaviour typical of outliers (Figs. 6 and 8). This final cut also makes it impossible to fit the highest velocity dispersion group (320 < σ [km s⁻¹] < 400), as we are left with no valid points. In any case, this group could be physically problematic as these velocity dispersions would be associated (Cappellari et al. 2006; Zahid et al. 2016) to extremely massive galaxies of close to 10¹² M_⊙ in stellar mass.

3.2. Lick index measurements

The original set of Lick indices (Worthey et al. 1994; Trager et al. 1998) includes 21 features identifying atomic and molecular absorptions typical of the light of SSPs of passively ageing galaxies, name-ordered in increasing wavelength from CN₁ to TiO₂. Then, the Balmer −δ and −γ indices were also included (Worthey & Ottaviani 1997) given their sensitivity to stellar ages, and they were given two definitions: a wide band labelled with _A and a narrower one labelled with _F. All indices, whether atomic or molecular, are defined by a central band (λ_cb and λ_cr), over which an integral was performed, and two additional bands bluewards (λ_bb and λ_br) and redwards (λ_rb and λ_rr) of the central region. These were used to define the pseudocontinuum, estimated as a straight line joining the median flux over each of the two bands, that is,

$$\begin{array}{c}\hfill F_{\text{p-c}}(\lambda )=\frac{F_{\text{r}}(\lambda -{\lambda }_{\text{b}})-F_{\text{b}}(\lambda -{\lambda }_{\text{r}})}{{\lambda }_{\text{r}}-{\lambda }_{\text{b}}},\end{array}$$(8)

where λ_b, λ_r, F_b and F_r are the central wavelengths and the median fluxes of the pseudocontinuum regions. Then, indices are computed as

$$\begin{array}{c}\hfill I_{\text{at.}}={\displaystyle {\int }_{{\lambda }_{\text{cb}}}^{{\lambda }_{\text{cr}}}}(1-\frac{F(\lambda )}{F_{\text{p-c}}(\lambda )})\text{d}\lambda ,\end{array}$$(9)

if they are atomic, or

$$\begin{array}{c}\hfill I_{\text{mol.}}=-2.5log\left(\frac{1}{{\lambda }_{\text{cr}}-{\lambda }_{\text{cb}}}{\displaystyle {\int }_{{\lambda }_{\text{cb}}}^{{\lambda }_{\text{cr}}}}\frac{F(\lambda )}{F_{\text{p-c}}(\lambda )}\text{d}\lambda \right),\end{array}$$(10)

if they are molecular. Atomic indices have units of wavelength and are defined in the same way as EWs, while molecular indices are dimensionless. With these definitions, any other atomic or molecular absorption (or emission) can also be defined. Indeed, we used the atomic definition for the computation of the emission in [O II] λ3727, [O III] λ5007 and Hα for the spectroscopic part of the selection of sources, while Hβ is already included in the original group of Lick indices. All index measurements were performed using the qcrpyLick code developed by Borghi et al. (2022a). By default, the code also predicts other features beyond the set of 25 Lick indices; namely the two CaII (both K and H) the D4000 (Hamilton 1985) and the D4000_n (Balogh et al. 1999). The errors associated with the indices are computed through a Monte Carlo resampling of the spectra using the uncertainties of the fluxes.

3.3. Spectral resolution and velocity dispersion

The index predictions from the models have an intrinsic resolution equal to that of the stellar library on which they are based. For this work, they were both built from the MILES library, which has an rms wavelength resolution of λ_IR ≈ 1.06 Å. For a fair comparison, index measurements should be at the same resolution as the models. However, SDSS spectra have an effective resolution arising from the convolution of the intrinsic SDSS instrumental resolution, and then galaxy spectra suffer from the intrinsic effect of velocity dispersion of the stars in the galaxy. Therefore, in principle, the data and/or model spectra need to be corrected so that resolutions are matched. In particular, given a flux density F_R0 at a resolution R₀, the downgraded spectrum at a resolution R is given by the convolution

$$\begin{array}{c}\hfill F_{\text{R}}(\lambda )=\frac{1}{\sqrt{2\pi {\sigma }_{\text{match}}^2}}{\displaystyle \int \text{d k}}{F}_{{\text{R}}_0}(\text{k})exp(-\frac{1}{2}{\left(\frac{\lambda -\text{k}}{{\sigma }_{\text{match}}}\right)}^2),\end{array}$$(11)

where σ_match is the Gaussian filter rms that needs to be applied in order to get the desired resolution. The value of this quantity can be obtained as

$$\begin{array}{c}\hfill {\sigma }_{\text{R}}^2={\sigma }_{{\text{R}}_0}^2+{\sigma }_{\text{match}}^2={\sigma }_{\text{VD}}^2+{\sigma }_{\text{IR}}^2+{\sigma }_{\text{match}}^2,\end{array}$$(12)

where σ_R is the aimed final rms resolution. We identified the present resolution, σ_R0, of a (data) spectrum with the quadratic sum of both the instrumental resolution and velocity dispersion resolution before the degradation was performed. Velocity dispersion resolution is linear in wavelength and takes the form σ_VD = σ_VD(σ, λ) = λ ⋅ (σ/c), where σ is the velocity dispersion. In particular, if we wanted to take a spectrum with a certain velocity dispersion σ₀ and instrumental resolution σ_IR0 to a larger velocity dispersion σ_f and instrumental resolution σ_IRf, one should perform the convolution (11) with a filter

$$\begin{array}{cc}\hfill {\sigma }_{\text{match}}& =\sqrt{{\sigma }_{{\text{IR}}_f}^2+{\sigma }_{\text{VD}}^2({\sigma }_f,\lambda )-{\sigma }_{{\text{IR}}_0}^2-{\sigma }_{\text{VD}}^2({\sigma }_0,\lambda )},\hfill \\ \hfill & =\sqrt{{\sigma }_{{\text{IR}}_f}^2-{\sigma }_{{\text{IR}}_0}^2+{\lambda }^2\frac{{\sigma }_f^2-{\sigma }_0^2}{c^2}}·\hfill \end{array}$$(13)

Fig. 4 shows the instrumental resolution of SDSS for sources at different redshifts (for the particular case of 200 < σ [km/s]< 250), which is a non-trivial function of wavelength, characterised by a sharp step-function-like jump at around 5500 Å in observed wavelength. In particular, the rms instrumental resolution of the data increases from ∼1.0 Å on the blue side to ≳1.35 Å on the red side. It should be noted that, at very low redshifts, this implies that the instrumental resolution of the data is similar to that of the MILES stellar library within the range of interest for Lick indices. While this would mean that no instrumental resolution correction would be needed for these sources, the situation changes with redshift and the effect becomes apparent, as the reddest indices progressively enter the lower resolution regime. Although the effect of this correction is only minor, it is appropriate to take it into account. Therefore, all science (data) spectra will be degraded to a common instrumental resolution of 1.5 Å; that is, with σ_match = (1.5² − σ_IR0²(λ))^1/2.

Fig. 4. RMS wavelength resolution (σR) of SDSS stacked spectra as a function of rest-frame wavelength for different redshifts.

Nevertheless, the procedure to then fit Lick indices with the SPS models depends on the model used. When fitting with the Knowles model, we leveraged the provided model spectra by degrading them to match the velocity dispersion and instrumental resolution of the object we fitted in each case. In other words, we created an index-strength model in which velocity dispersion is just a configuration parameter. The transition from the (appropriately smoothed) spectra to Lick index predictions is always performed with the pyLick code (Borghi et al. 2022a). On the other hand, when using the TMJ model, an intermediate correction of the data needs to be performed. This model is built as a direct prediction of index strengths at MILES (1.06 Å) resolution, for which we need to derive a correction (covering both a degraded instrumental resolution of 1.5 Å and velocity dispersion) that is applied to the indices once they have been measured. The goal is to find a function C_I(σ) that takes an index I measured from a spectrum of a galaxy with velocity dispersion σ and intrinsic instrumental resolution σ_IR to the value I₀ it would be expected to have if measured at a MILES instrumental resolution and without the effect of velocity dispersion, mirroring the output of the TMJ model. Depending on whether the index is atomic or molecular, the correction is written as I₀ = I(σ)⋅C_I(σ) or I₀ = I(σ)+C_I(σ), respectively. In other words,

$$\begin{array}{c}\hfill C_I(\sigma ,{\sigma }_{\text{IR}})=\{\begin{array}{cc}I(0,{\sigma }_{\text{MILES}})/I(\sigma ,{\sigma }_{\text{IR}}),& \text{if atomic,}\\ I(0,{\sigma }_{\text{MILES}})-I(\sigma ,{\sigma }_{\text{IR}}),& \text{if molecular,}\\ \end{array}\end{array}$$(14)

where σ_IR = 1.5 Å. In order to compute the correction, several methods exist in the literature. Following Carson & Nichol (2010), we took the SDSS spectra³ of giant stars in M67 (NGC2682) and performed smoothings of progressively higher velocity dispersions. Then, the change in the indices was averaged over all stars, and the sampling of C_I(σ) was parametrised by a polynomial of order 3 for all indices. It should be noted that I₀ can be measured directly from the stars even though they have the SDSS resolution imprinted in their spectra, since, similarly to what happens with low redshift galaxy spectra, the entire set of Lick indices falls below the jump in resolution.

We highlight that the corrections C(σ) carry an associated uncertainty δC, which we calculated during their derivation and took into account. These uncertainties are small compared to the statistical uncertainty of the data (index measurements) and are treated altogether as statistical uncertainty, despite the systematic nature of the former. For further information, we include Appendix A, where we show the correction functions C(σ) we used and give the parameters of the 3rd order polynomials. We also show for reference the difference in stellar age estimated for the set of data due to the use of the two methodologies presented to address the velocity dispersion effect. We note that only the Knowles-modelled ages can be used for this particular study, as taking the TMJ model to different velocity dispersions is not possible.

3.4. Stellar population models

The measured indices were fitted via an MCMC algorithm to both the TMJ and Knowles models. The underlying assumption is that the spectra can be approximated as a coeval stellar population born in a single burst, which is the rationale behind the use of cosmic chronometers. While recent studies fitting Lick indices have typically employed direct Lick index models such as TMJ, we decided to also test the sMILES-based SSP model to compare the results. As mentioned in the introduction, we advise that the choice of SPS model can alter significantly the results (Moresco et al. 2020) in matters of both age and H estimations. For what concerns the specificities of the models included, TMJ uses the evolutionary synthesis code of Maraston (2005), the (Salpeter 1955) IMF, the MILES empirical library, and (Cassisi et al. 1997) stellar tracks. The Knowles model includes a larger semi-empirical library based on MILES, known as sMILES (Knowles et al. 2021). This library is built to correctly reproduce the MILES stellar spectra at the typical metallicities and α-abundances of the stars in our neighbourhood, but using theoretical tools to go beyond those restrictions in the parameter space. They use the BaSTI isochrones (Pietrinferni et al. 2004, 2006, 2021) and the model is available for five different IMFs, from which we chose the Chabrier IMF (Chabrier 2003). These models can be interpreted as functions such as

$$\begin{array}{c}\hfill I_i\equiv I_i(t,[\text{Z}/\text{H}],[\alpha /\text{Fe}]),\end{array}$$(15)

where the subindex _i denotes a given Lick index, which are a function of age (t), metallicity ([Z/H]), and abundance of α elements ([α/Fe]). Each ‘base’ model is sampled on a specific grid of parameters, which yields a discreet function that needs to be interpolated for the fit. Once this is done, the age, metallicity, and α-enhancement are found via an MCMC algorithm. In particular, we used emcee (Foreman-Mackey et al. 2013), the Python implementation of the Goodman & Weare affine invariant MCMC ensemble sampler (Goodman & Weare 2010). The log-likelihood is written as

$$\begin{array}{c}\hfill log\mathcal{L}=-{\displaystyle \sum _i}log\sqrt{2\pi {\sigma }_{I_i}^2}-\frac{1}{2}{\displaystyle \sum _i}{\left(\frac{{\widehat{I}}_i-I_i}{{\sigma }_{{\widehat{I}}_i}}\right)}^2,\end{array}$$(16)

where ${\sigma }_{{\widehat{I}}_i}$ are the uncertainties associated with the Lick index strength measurements, ${\widehat{I}}_i$, and I_i are the model indices from Eq. (15). Once a particular model for the fit is selected, one should pay attention to two key aspects: the set of indices chosen for the fit, and the grid onto which the base model I_i(t, [Z/H],[α/Fe]) is interpolated. Regarding the former, the recommendations provided by the authors should be followed not to bias the results. For the TMJ model, following their study using globular clusters without individual element variations, the use of certain indices is discouraged since they are not well reproduced by the data. In essence, our set comprises all of the indices of the baseline used by Johansson et al. (2012), that is, Hδ_A, Hδ_F, Fe4383, Fe4531, Hβ, Mgb, Fe5270, Fe5335, Fe5406; plus G4300, Hγ_A, Hγ_F, and the Mg₂, since they are well reproduced by globular cluster data in Thomas et al. (2011). In turn, for the Knowles model, we discarded indices known to be poorly modelled, too broad, or sensitive to the IMF, following Knowles et al. (2023). The set of Lick indices is thus Hδ_A, Hδ_F, G4300, Hγ_A, Hγ_F, Fe4383, Ca4455, Fe4531, Hβ, Mgb, Fe5270, Fe5335, Fe5406, Fe5709, Fe5782. We emphasise at this point that using a fixed set of indices across the entire dataset is strongly recommended, as Borghi et al. (2022b) found age estimates could vary by up to ∼2 Gyr depending on the indices used.

In general, Balmer indices are sensitive to age, iron indices to metallicity and magnesium indices to α-enhancement. Knowles et al. (2023) showed that the plane Hβ₀-[MgFe] could separate, almost orthogonally, ages and metallicities, where Hβ₀ is an optimised definition of Hβ. We note, however, that for our selection of very passively evolving ETGs, the Balmer-β feature alone produced extremely old ages, as Vazdekis et al. (2015) had already shown. Instead, the combination of the three Balmer lines (β, γ, δ) returns much more reasonable ages. When it comes to metallicity and α-enhancement, Knowles et al. (2023) showed that the plane Mg_b-Fe5335 or -Fe5270 separated well between different alpha enhancement groups. The three of them are included and recommended for their use in stellar population fits by themselves and Thomas et al. (2011). Moreover, a combination of these three indices makes the [MgFe] index, which is sensitive to metallicity but insensitive to α-enhancement.

Regarding the grid of the models, the trade-off between the computational cost and the precision when using an overly fine interpolating grid should be addressed. The base grid for the TMJ model features ages distributed in the interval [0.1,  15.0] Gyr, linearly spaced with steps of 0.1 and 0.2 Gyr between 0.1 and 1.0 Gyr and with steps of 1.0 Gyr until 15.0 Gyr. Metallicity ([Z/H]) is sampled in six values, namely { − 2.25, −1.35, −0.33, 0.00, 0.35, 0.67}. Finally, α-enhancement ([α/Fe]) is sampled in four values: { − 0.3, 0.0, 0.3, 0.5}. Regarding the Knowles model, the base grid features ages distributed in the interval [0.03, 14.00] Gyr, with a varying step size that increases with age, namely 0.01, 0.05, 0.1, 0.25 and 0.5 Gyr. Metallicity is sampled over ten values not evenly spaced from −1.79 to 0.26, meaning it does not reach strong metallicities beyond 0.26; in contrast to TMJ, which goes up to 0.67. α-enhancement is equally the least thoroughly sampled parameter, with five values: { − 0.2, 0.0, 0.2, 0.4, 0.6}.

Some of these values are typically not attained for massive passive galaxies (Parikh et al. 2021; Saracco et al. 2023). For example, metallicity is expected to range from slightly subsolar to supersolar metallicities, but not beyond [Z/H]∼0.3 (Clemens et al. 2006, 2009; Thomas et al. 2010; Maiolino & Mannucci 2019). Moreover, we find that allowing for metallicities that go beyond the limit given by the Knowles model could bias a portion of the sample into smaller ages due to the age-metallicity relation and the poor sensitivity of the models at high stellar ages, which is not resolved by employing a thinner grid. For this reason, in addition to the computational cost, we decided to restrain the parameter limits of both models to t ∈[0.1,  14.0] Gyr, [Z/H]∈[−0.35,  0.26] and [α/Fe]∈[0.0,  0.4] and interpolate into a finer grid with step sizes of Δt = 0.1 Gyr, Δ[Z/H] = 0.01 and Δ[α/Fe] = 0.01. This choice of parameter limits is also supported by Fig. 5, where a third finer version of the base models (with 50 equally spaced points for each parameter in the original extent of the grids) was run on the galaxy stacks (shown in orange) as a preliminary test. As expected for massive ETGs, metallicities are either solar or supersolar and there is a significant abundance of α-elements with respect to iron. For comparison, in blue, we show the result for the fit to the spectra of single galaxies. The distributions in this case are much more scattered and do not reflect the expected behaviour for massive passively evolving galaxies. Indeed, subsolar metallicities and iron enhancement with respect to α-elements are found for a significant fraction of galaxies, especially when using the Knowles model. Additionally, extremely high values for the ages are found in a large number of galaxies, with an additional seemingly young population of ∼3 Gyr obtained through the TMJ model. This underscores the importance of using high S/N spectra in these studies, as results behave as expected when the fits are performed with the stacked data.

Fig. 5. Age, metallicity, and α-enhancement distributions for the TMJ (left panels) and Knowles (right panels) models, interpolated to a grid of 50 points in each parameter within the original ranges. Results are shown for the Lick index fitted to both single galaxies (in blue) and stacks (in yellow) of the final cosmic chronometer sample.

3.5. Cosmographic fit

The last step is to fit the stellar age-redshift relation in different velocity dispersion bins to derive the cosmographic parameters (H₀, q₀, j₀). This is performed using (4) and (6), which allow us to sample the Hubble parameter up to z ≲ 0.35 without assuming a cosmological model. As anticipated at the end of Section 3.1, we divided the sample into six logarithmic dispersion bins; namely [100,125), [125, 160), [160, 200), [200, 250), [250,320), and [320, 400) km/ s⁻¹. In principle, this choice implies that six sets of age-redshift data should be fit jointly. Since every set characterises a different cosmic chronometer population defined by its mass, six distinct t₀ parameters need to be defined (i.e. t_{0, 100 − 125}, t_{0, 125 − 160}, t_{0, 160 − 200}, t_{0, 200 − 250}, t_{0, 250 − 320}, t_{0, 320 − 400}), which stand for the average stellar age of each ensemble at z = 0. With all this, the log-likelihood reads

$$\begin{array}{c}\hfill log\mathcal{L}=-{\displaystyle \sum _{i,v}}\frac{1}{2}{\left(\frac{{\widehat{t}}_{i,v}-t(z_{i,v};t_{0,v},H_0,q_0,j_0)}{{\sigma }_{{\text{t}}_{i,v}}}\right)}^2,\end{array}$$(17)

where i runs through the stacks in each velocity dispersion group, represented by v. For all parameters, we assigned uniform prior probability distributions. The limits for these distributions were set between 0 and 20 Gyr for the zero-ages (t_0, v), 20 and 120 km s⁻¹ Mpc⁻¹ for the Hubble constant, −10 and 10 for the deceleration parameter and −100 and 100 for the jerk parameter. The choices for the cosmographic parameters were made in order to guarantee that the posterior distributions were constrained in the main cases (TMJ, Knowles), even if this allows the tails of the q₀ and j₀ distributions to attain values far from the conventional ones found in the literature. As shown below, it is safer to do this due to the low sensitivity of the model on these kinematical parameters, particularly the jerk.

Additionally, we required the derivative of the Hubble parameter, dH/dz, to be positive, as suggested by the global trend observed from model-independent measurements of H(z) (Simon et al. 2005; Stern et al. 2010; Zhang et al. 2014; Ratsimbazafy et al. 2017; Jiao et al. 2022; Tomasetti et al. 2023; Moresco et al. 2012, 2016; Borghi et al. 2022a). Physically, this translates into allowing the expansion of the Universe to accelerate but not too much, with a limit on q(z) >  − 1 at any redshift. In particular, the deceleration parameter today, q₀, has to follow this condition. This could be interpreted, under the introduction of a w₀CDM model, as the condition of reducing the amount of phantom dark energy states available. Specifically, the mathematical condition is w₀ > 1/(Ω_m − 1), which allows for w₀ < −1 only when Ω_m ≫ 0. In particular, w₀ would present a lower limit in ≈ − 1.43 for Ω_m = 0.3. Under the assumption of the w₀CDM model, we can combine our estimations for q₀ and j₀ (Capozziello et al. 2011) to get the estimations of Ω_m and w₀ in those models, following

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_m(q_0,j_0)=\frac{2(j_0-q_0-2q_0^2)}{1+2j_0-6q_0}\end{array}$$(18)

and

$$\begin{array}{c}\hfill w_0(q_0,j_0)=\frac{1+2j_0-6q_0}{-3+6q_0}·\end{array}$$(19)

Note that even if we apply these transformations, we can obtain values for Ω_m that exit the range [0, 1]. With these transformations, we can have a hint of the blind preference for values of Ω_m and w₀ aligned with the results from model-dependent probes.

4. Results

The present section covers the cosmographic analysis of the data, after the fit using the two SPS models. Using TMJ-modelled ages, we find a credible interval for H₀ aligned with other probes in the literature. Further discussion is provided for the other parameters included in the fit (q₀, j₀), as well as for the cosmographic readout of the Knowles-modelled ages. Through the introduction of a w₀CDM model, we can get some insight into the implications of our estimation for the parameter of the equation of state of dark energy, w₀. Finally, we examine the posterior sampling of our cosmographic fit in the Hubble diagram and compare our results with the point observations in the literature.

The right panels of Figs. 6 and 8 show the ages derived using the TMJ and Knowles models, respectively, where each data point corresponds to a different stack. We found that more massive galaxies are older at every redshift, and, at fixed velocity dispersion, the age decreases with redshift. This is in agreement with the cosmic chronometer ground principle: for a given velocity dispersion, we are statistically observing the same population of passive galaxies across cosmic time, and thus the age evolution of the same galaxy population can be traced. The stratification in velocity dispersion groups is clear and follows the downsizing scenario framework for the formation of ETGs. This aspect will be treated in detail in Section 5.1.

Fig. 6. Joint cosmographic fit of t(z) with ages estimated using the TMJ model. On the left, the corner-plot for the posterior probability distributions of H0, q0, and j0 is shown, displaying 1σ and 2σ contour levels. On the right, the t − z dispersion is shown. The empty circles in strong colour represent the data points we used for the fit, while the faded crosses represent the points, resulting from of a set of galaxies with average dispersion in redshift greater than 0.01. The shaded regions in the same colour as the data points are the 1σ of the sampling of the posterior distribution coming from the parameters on the left. The zero-ages of the velocity dispersion groups are not shown here, but we refer to the corner plot in Fig. B.2; however, they are used for the sampling and to set a zero-age for the Planck Collaboration VI (2020) evolution of t(z), shown as dashed black lines. The dashed q0 = −1 lines in the corner plots for the planes q0 − H0 and j0 − q0 represent the natural limit imposed by dH/dz > 0. This is the baseline t − z relation referred to in this work.

Before proceeding to the cosmographic fit of the ages, some points need to be addressed. First, the uncertainties in the ages derived from the Knowles model are systematically much smaller by a factor of ∼6 − 7. To investigate the potential reasons behind this, we examined the differences in procedure when deriving ages with the two models. The first difference is the set of Lick indices to be fit, as explained in Section 3.4. However, after fitting a common set of indices for the two models (defined by their overlap), no relevant changes are observed in ages or their uncertainties. The second difference lies in the treatment of the velocity dispersion effect. While for the TMJ model, we had no other option than applying a correction–which introduces an extra systematic uncertainty associated with the computation of the correction–for the Knowles model, we could directly fit the data with the model indices, in which velocity dispersion is applied a priori to the synthetic spectra. However, using the TMJ methodology with the Knowles model is possible, allowing us to rule out the treatment of the velocity dispersion as the reason behind the disagreement on the value of uncertainties. In particular, we tried fitting the Knowles model at 0 velocity dispersion with the indices (data) corrected using the velocity dispersion corrections, analogously to the exact procedure of the TMJ model. The results were negligible changes in the ages, and the uncertainties remained the same size (further details are given in Appendix A). Therefore, the only explanation is the nature of the models themselves; that is, the Knowles model seems to be much more precise. However, this will not necessarily play in our favour, due to the intrinsic scatter in the data.

The second point to address is the redshift dispersion within the stacks. As depicted by the horizontal error bars in the right panels of the aforementioned Figures 6 and 8, there is a (small) number of stacks (plotted with bold faded crosses) containing galaxies whose redshift covers a relatively wide range. Due to the stacking procedure, this typically occurs at the highest redshifts for a given velocity dispersion bin or for stacks with 320 < σ [km s⁻¹] < 400. Since these results are less reliable, we discard them from the start.

Fig. 7. Posterior probability distribution for H0, Ωm, and w0. The sampled parameters were q0 and j0, then converted to w0CDM parameters through Equations (18) and (19). The sampling was run under the ΛCDM-model condition of 0 < Ωm(q0, j0) < 1. The brown hatched region in the Ωm − w0 plane represents the region allowed by dH/dz > 0 ⇔ w0 > 1/(Ωm − 1).

The last point to cover is related to a fluctuating pattern found in the ages with respect to redshift. Indeed, as seen more clearly for the TMJ model (Fig. 6), the ages are found to oscillate at various redshifts, deviating from the expected smooth decline with redshift. This issue is particularly enhanced at 0.16 ≲ z ≲ 0.19. This one particular ‘oscillation’ is observed across all three velocity dispersion bins for which there are data at these redshifts. While we do not expect a large impact on the cosmographic results from the fit due to its symmetry with respect to the global trend, we found this issue of great relevance and therefore decided to discuss it in depth in Section 5.2.

Once the distribution of ages has been individually assessed, we can turn to the cosmographic fit. In the right panel of Fig. 6, apart from the t(z) scatter of the data, we also show the 1σ sampling of the posterior distribution for t(z), as obtained from the MCMC algorithm with the standard priors (0 < t_0,i < 20 Gyr, 20 < H₀ < 120 km s⁻¹ Mpc⁻¹, −10 < q₀ < 10, −100 < j₀ < 100; as previously discussed in Section 3.5). For comparative purposes, we also included the ΛCDM derived t(z) following the results from Planck Collaboration VI (2020) as dashed black lines, taking as t_0, v the peak of each posterior distribution on these parameters. Notice that the exclusion of widely z-extended data points eliminates all three belonging to the 320 < σ [km s⁻¹] < 400 group, making it impossible for us to estimate t_{0, 320 − 400}.

The corner plot corresponding to this sampling is shown in the left panel of the figure. The zero-redshift ages (t_0, v) of the individual cosmic chronometer samples are not plotted for the sake of a clear view of the other posteriors, but we anticipate they are extremely well constrained–at the exception of, obviously, t_{0, 320 − 400}-, as seen in Appendix C. We also note here that the estimations that we give throughout this discussion for the fit parameters include only statistical uncertainties. Systematic effects in both our methodology and from the SPS models will be addressed later (discussion around Fig. 9).

Fig. 8. t − z relations for the ages estimated using the Knowles model. The structure of this figure is the same as that in Fig. 6. The grey contours represent the case with fv factors included, to account for the possible underestimation of age uncertainties. The t − z samplings on the right panel refer to this case. The contours in purple represent the case where these factors were not taken into account. The sampling of the associated posterior for t − z is not represented on the right panel for this case.

Fig. 9. Reconstruction of the Hubble parameter. Our sampling is shown in shades of blue, corresponding to the posterior distribution obtained for H0, q0, and j0 from Fig. 6. The width of the tendency represents the range between the 16th and 84th percentiles of our posterior sampling. The opacity of the colour shade represents the reliability of the H(z) measurement, which we based on the amount of sources (individual galaxies) at each z. The dashed black line represents the H(z) evolution according to Planck Collaboration VI (2020). The plot is extended up to z = 0.9 in order to show the result at the lowest redshift obtained through the Lick indices fit (dark red), by Borghi et al. (2022b). Other measurements from the available literature are represented in yellow (obtained through FSF) and orange (obtained through the measurement of D4000 of its narrow counterpart D4000n).

The Hubble constant, H₀, is constrained to a mean value of $70.0_{-7.6}^{+4.1}\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$, which includes both the Planck Collaboration VI (2020) and Riess et al. (2022) estimations within the credible interval. The estimation of q₀ is relevant, clearly rejecting values below −1, having a seemingly strong cut-off at that value. For a direct view, we drew this limit in the corner plots presenting the q₀ − H₀ and j₀ − q₀ planes. As we explain further on, this is due to the limitation on q(z) (and thereby on q₀) imposed by dH/dz > 0. To inspect this in further detail, we examine the 16th and 84th percentiles and find that the modal value (maximum a posterior probability), mod(q₀) =  − 0.95, is to the left of the percentile 16, P₁₆(q₀) =  − 0.85, while the percentile 84 lies at P₈₄(q₀) = 0.62. It is clear that the posterior is skewed towards q₀ = −1 from its right. Indeed, when we inspect the MCMC chain thoroughly, we find that 99.91% of the chain chose values greater than or equal to −1, with only 0.09% of samples remaining below, possibly due to the remainder in the approximation of H(z). As previously stated in Section 3.5, we interpret this as proof that the second order Taylor expansion of H(z) is a good approximation in the redshift window we treat. Knowing this, we do not expect the region q₀ < −1 to be explored in any case by the MCMC chain, but we kept the lower limit of the prior in −10 as a consistency check. In order to assess how much this restriction on the derivative of H means at a statistical level, we explored an extended parameter (0 < H₀ [km s⁻¹ Mpc⁻¹] < 300, −20 < q₀ < 20, −100 < j₀ < 300) space after lifting the condition. We found that marginals were much wider, in some cases beyond the limits set as standard, and a clear anti-correlation between H₀ and q₀ that effectively extends to q₀ < −1 (posterior shown in Appendix B). The minimum of the χ² in this case occurs at H₀ ∼ 220 km s⁻¹ Mpc⁻¹, but the typical values in this and all the regions explored clearly outside ΛCDM (H₀ ∼ 70 km s⁻¹ Mpc⁻¹, q₀ ∼ −0.5, j₀ ∼ 1) are similar to those found in the well-delimited posteriors in Fig. 6 (the baseline).

Finally, we assert that the cosmographic fit is highly insensitive to the jerk parameter, j₀, as the window needed in the prior to obtain a constrained posterior covers a range equivalent to more than ten times the value required by ΛCDM, j₀ = 1. Our estimation for the parameter is $j_0=0.9_{-3.7}^{+5.5}$, and the related corner plots show no degeneration of the jerk with either of the other two cosmographic parameters. However, we found, in the unrestricted fit without the limitation on the derivative of H, a seemingly functional relation between the preferred combinations of H₀ and j₀ that accepts 50 < H₀ [km s⁻¹ Mpc⁻¹] < 100 only when −50 < j₀ < 50.

Our assumption that the oscillations in redshift do not affect the global fit is supported a posteriori by an additional run of the cosmographic fit in which we excluded the data points with redshifts in the range 0.15 < z < 0.20. The marginal posterior distributions of the parameters are similar to those observed in Fig. 6, with the median and minimum χ² in the last 1000 steps of the MCMC chain being roughly half the value found in the original case, the same order of magnitude. The marginals and corner plots for this case can be found in Appendix B.

So far, results remain model-independent in the framework of cosmology, but with the transformations in Equations (18) and (19) we can take a look at the implications for a w₀CDM model. To this end, we decided to run an extra MCMC imposing 0 < Ω_m(q₀, j₀) < 1. We do so because, if we depended solely on the results for the sampling shown in Fig. 6, the regions Ω_m > 1 and Ω_m < 0 would be populated, adding spurious probability to w₀ and H₀ coming from the sampling there. The probability density regions of the figure show the sampling when 0 < Ω_m < 1 is enforced, and the brown hatched area represents the region allowed by w₀ > 1/(Ω_m − 1). This condition on the Ω_m − w₀ plane is the result of having imposed dH/dz > 0. Due to the effect of the smoothing of the contour levels, it might seem that there is a non-negligible amount of sampling at w₀ < 1/(Ω_m − 1); however, when the steps of the MCMC in that region are counted, we find only 0.19% of the total amount. This result is similar to the 0.09% of samples that had q₀ < −1 in the paragraph above, and it can be attributed to the remainder in the approximation of H(z) for its second order Taylor expansion. The newly constrained H₀ keeps almost the exact same posterior distribution as before, constrained now to $H_0=70.7_{-5.1}^{+4.3}\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$. Ω_m is mostly unconstrained, although central values seem to be preferred. Finally, we have a readout of w₀, with a posterior probability defined by $w_0=-1.{01}_{-1.75}^{+0.03}$. This result is compatible with ΛCDM, although the fit seems to prefer the combinations of w₀ and Ω_m close to the limit w₀ = 1/(Ω_m − 1), tightly limiting the right-hand tail of the posterior of w₀.

We now turn to the Knowles model. We apply the same conditions we did for the sampling in Fig. 6. The direct application of the MCMC algorithm to the data yields the results shown in purple in the left panel of Fig. 8. The posteriors constrain the cosmographic parameters to $H_0=41.7_{-0.6}^{+0.6}\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$, $q_0=-0.{99}_{-0.01}^{+0.08}$ and $j_0=55.4_{-2.4}^{+2.6}$. The result for q₀ follows what we already observed in the TMJ case. The other two marginals, however, are well off and extremely thin, distancing from literature values by several σ. Such precise results for parameters so strongly deviating from their expected value prompt a careful reconsideration of all possible sources of error. First, we note that the narrowness of the posterior distributions is primarily driven by the size of the observational errors. Then, after analysing the last thousand iterations of the MCMC chain, we find a discrepancy of three orders of magnitude in the minimum χ² compared to the fit in Fig. 6.

This suggests that the ages inferred using the Knowles model may not be useful to gain sensitivity on the cosmographic parameters in their current form. We conclude that the nature of the Knowles model (Appendix A), affects the performance of the subsequent cosmographic fit by enlarging the scatter of the ages while reducing the value of their uncertainties. The overall effect is the restriction of the posteriors to clearly off values and worsening the overall χ² values of the sampling. If there was an underestimation of the ages that we cannot control, it can be addressed by introducing a set of f factors to the likelihood. We do so following the documentation of emcee, by substituting

$$\begin{array}{c}\hfill {\sigma }_{i,v}^2\to {\sigma }_{i,v}^2+t_{i,v}^2(z_{i,v};H_0,q_0,j_0)·f_v^2.\end{array}$$(20)

In essence, an additional parameter f_v is included in the MCMC for each velocity dispersion bin, with a view to capture the lacking variance at least to a certain degree. They are assigned wide uniform priors satisfying −10 < log f_v < 1. The corresponding corner plot for the cosmographic parameters is shown in grey in the left panel of Fig. 8, with the Hubble constant constrained to $38.2_{-5.1}^{+12.9}$ km s⁻¹ Mpc⁻¹. The result is still incompatible with other probes at 2σ, and the central value is lower than before. However, the width of the posterior allows for a 3σ agreement with H₀ ∼ 70 km s⁻¹ Mpc⁻¹. The posterior of q₀ is again skewed towards −1 (mod(q₀)∼ − 0.9, P₁₆(q₀) =  − 0.17), although this time the positive tail of the sampling extends to the upper limit of the prior. If we observe the corner plot relative to the plane q₀ − H₀ we can see a degeneracy that prefers higher H₀ values (still low in comparison to other probes) for smaller q₀ values. Finally, the model is still insensitive to j₀ as in the TMJ case, but we find a constraint ($j_0=5.8_{-5.1}^{+42.6}$), thanks to the introduction of f_v factors, that is compatible with j₀ = 1 within 1σ. We notice that this time the posterior of j₀ is right-skewed.

In the right panel of the figure, we omitted the 1σ sampling of t(z) without f_v factors. We did so to simplify the figure and avoid unnecessary visual noise, but they are tight and leave most of the data points out by several σ. For the rest, Fig. 8 contains the same elements as Fig. 6.

It is clear that f_v factors somewhat aligned the results with the expected values coming from other probes, even if they do not capture the full complexity of the underestimated variance. The power of these factors is clear, as the average value of the minimum χ² over the last thousand steps of the MCMC chain is now of the order of magnitude obtained for the TMJ fit. As a final comment on the Knowles model relative to the t − z relation, it could happen that the modelled ages are less sensitive to the index-oscillations observed in the TMJ case, as they are not as easily perceived on the right panel of Fig. 8. This idea could be tested further if the intrinsic scatter of the data was reduced, and if true, would constitute a strong feature of the Knowles model while the fluctuating index-strength issue is understood. So far, we believe that caution is advisable when using the Knowles model for cosmographic purposes, as it seems to present some irregularities regarding the optimal tracking of the evolution of ETG ages across redshift. Some of the difference between models in the evolution of the differential age with redshift can be attributed to the effect of model variation when it comes to the estimation of the Hubble parameter (Moresco et al. 2020). For what is known in the literature, the intra-model variation of the IMF is expected to be only a minor issue, but other age diagnostics or index selection could be inspected in future work.

4.1. Assessment of systematic uncertainties

So far, the uncertainties that have been included in the estimation of cosmographic parameters are mainly of statistical nature. Only the uncertainty associated with the velocity dispersion correction functions, which is added to the index measurements, is of a systematic nature. However, as discussed in Appendix A, the contribution to the error budget is small. In the end, it is treated as part of the statistical uncertainties, as the stellar model fit makes it impossible to get them disentangled. Apart from these, there are other sources of systematic uncertainty that have been thoroughly treated in the literature. We distinguish between those that are external to this work (e.g. SPS models, their configuration: IMF, isochrones) and those that arise from our methodology. We do not treat the former, performing the H(z) estimation using only the TMJ ages, as we conclude that the Knowles model needs further inspection before deriving conclusions from it. For the rest of this section, all computations refer to ages estimated using the TMJ model.

Following the stream of our methodology, we can cite the following three sources of systematic uncertainty: the S/N level of the stacked spectra, the set of Lick indices used for the stellar parameter fit, the choice of the final reliable dataset for the cosmographic fit (effect of outliers and internal homogeneity of the stacks). We run a test on the S/N level and the redshift extension of the data points, to assess their contribution to the systematic error budget. Addressing the impact of the set of Lick indices chosen for the fit would require a complicated analysis in which several combinations of indices, selected to be sensitive enough to variations on the SPS model parameters, are included. Beyond the feasibility of the study, this work does not aim to test the performance of index combinations in the TMJ model, which was already done by Johansson et al. (2012). Indeed, we chose our indices as the baseline described in that work, plus another four indices that Thomas et al. (2011) proved to be well calibrated with globular clusters. The most recent study by Borghi et al. (2022b) showed that the contribution of the set of indices to the systematic error budget in H(z) would be of the order of that related to the binning (in our case the stacking). Finally, the exclusion of data points that averaged galaxies over a large redshift interval is addressed by performing the cosmographic fit both releasing and tightening the threshold for exclusion in σ_z.

In order to build the stacks we followed the methodology described in Sect. 3.1, which leaves the selection of the S/N level as the sole arbitrary choice. We assess its contribution to the systematic uncertainty by creating 11 groups of stacks, built up to have S/N levels from 200 to 400 in steps of 20. The central group at S/N ∼ 300 is the one that we used as reference throughout the paper. We obtain the posterior distributions from the cosmographic fit in each case, compute the median, over the last 1000 steps of the MCMC chains, of the standard deviation between the 11 groups at each redshift. This is

$$\begin{array}{c}\hfill {\sigma }_{H,\text{syst}}(z)={\text{med}}_k\left(\sqrt{{\text{Var}}_i\left(H_i(z;H_{0k},q_{0k},j_{0k})\right)}\right),\end{array}$$(21)

where the index _k represents the sample from the MCMC chain and the index _i represents each of the 11 stacking groups. Equation (21) is also used to evaluate the contribution to the systematic uncertainty introduced by the choice of the σ_z threshold when excluding data points prior to the cosmographic fit. For this, we consider threshold values of σ_z = {∅,0.05, 0.02, 0.01, 0.005, 0.002}, where we use ∅ to denote the absence of any threshold (i.e. the inclusion of all stacks). The contribution of the two sources of systematic uncertainty to the total is shown in Fig. 10, in the following section, as it is presented in units of the statistical uncertainty of the posterior sampling of H(z).

Fig. 10. Systematic uncertainties arising from the methodology directly accounted for in this work. They are presented in units of 1σ of the statistical uncertainty σstat at each redshift. The individual contributions of the level of S/N selected in the stacking procedure and the limit to the dispersion of sources in redshift in each stack are presented in yellow and dark orange, respectively. In dark blue, we present the overall systematic uncertainty.

4.2. The sampling of H(z)

As a final result, we present in Fig. 9 the sampling of H(z), using the posteriors from the baseline, in which we fit the TMJ-modelled ages imposing dH/dz > 0; that is, the MCMC run shown and discussed in Fig. 6. At this point, we introduce the systematic uncertainties that we discussed in the prior section. We do so following the typical quadratic sum: σ_total² = σ_stat.² + σ_syst.². In Fig. 10, each of the contributions taken into account is presented in units of 1σ of the statistical uncertainty from the sampling of the MCMC run. These are quite stable and slightly greater than the level of statistical uncertainty. The one associated with the S/N level chosen in the stacking procedure (yellow) dominates at low redshift, while the one related to the exclusion of data points (dark orange) is slightly greater at higher redshifts. Overall, the systematic uncertainty is stable and of the order of twice the statistical uncertainty.

In Figure 9, the final credible interval for H(z) is presented in blue, where we chose to plot it in different shades according to the number of sources (individual galaxies) we had at each redshift. The light blue represents the contribution of all sources of uncertainty that have been studied in this work, while the darker blue represents only the 1σ statistical uncertainty interval, computed as the region between the 16th and 84th percentiles of the sampling of the MCMC fit (parameters sampling in Fig. 6). As a reference, we also included, as dashed light blue lines, the limits of the 1σ credible interval in case the main sources of systematic uncertainty not computed in this work (Lick index selection, SPS models) were collectively the same size (Borghi et al. 2022b) as the systematic uncertainty taken into account.

For the sake of comparing the potential of this work with literature, we include in yellow, orange, and dark red the local model-independent measurements of H(z) available to date up to z = 1. The colour of the observations depends on the methodology used for the analysis in each case. Full spectral fitting (Simon et al. 2005; Zhang et al. 2014; Ratsimbazafy et al. 2017; Jiao et al. 2022), in yellow, shows greater uncertainties than ages derived using the age-sensitive D4000_n spectral feature (Moresco et al. 2012, 2016; Loubser et al. 2025), in orange. Our results extend along a redshift range (z ≲ 0.4) mostly populated by FSF measurements, and we see that our estimation is competitive in comparison. Within the same range, four/five local measurements through the D4000_n feature show a relatively greater precision than ours. Finally, we present the ΛCDM derived H(z) from Planck Collaboration VI (2020), which sits within the 1σ of our estimation. Naturally, our constraint power is stronger H(z) when more sources are available, and it loses precision towards z ≳ 0.3. This pushes us to consider the extension of this methodology to higher redshift by the use of BOSS and eBOSS (extended BOSS) data, using a Taylor expansion of H(z) of appropriate order.

5. Further discussion

This section is included to fully exploit the results shown in the previous section. In the first place, we perform the linear fit of the scatter of log(t), [Z/H], and [α/Fe] with log(σ), known as scaling relations. These functions try to give a quantitative guideline for the status of the currently well-established downsizing scenario, by which more massive ETGs tend to have formed earlier in cosmic history, be more metallic, and α-enhanced. We obtain these relations for both the SPS models we consider in this work, using the low redshift (z < 0.1) end of our sample of sources. The second part of this section includes the discussion around the oscillations observed in the t − z relation in Fig. 6. In particular, we focus on the coordinated fluctuation, in all central and statistically significant velocity dispersion groups, at 0.16 < z < 0.19. We unveil the tight correlation of this oscillation with a similar pattern on the index strength of some indices (Hγ_F and Fe4383) in that redshift range. A deeper inspection of the global trends of index strengths with redshift unveils an apparently overlooked behaviour that was already present in previous measurements of these Lick indices (i.e. Tremonti et al. 2004).

5.1. Stellar archaeology and scaling relations

We perform the study on the scaling relations by selecting all galaxies with redshifts 0.05 < z < 0.10, 18 955 objects, which makes up around a 21% of the whole sample. Given that the default stacks described in Section 3.1 are divided into six velocity dispersion bins (but we aim to derive a smooth relation from σ ∼ 100 to σ ∼ 350 km s⁻¹), we construct another set of stacks for the purpose of this subsection. We stacked the galaxies in the subsample in groups of 400 (ordered by velocity dispersion), since that number should ensure a high enough S/N at this redshift according to Fig. 3. Lick indices are then measured and fit using both the TMJ and Knowles models.

Fig. 11 shows the stellar population parameter trends in terms of velocity dispersion for the low-redshift sample. In particular, the results using the TMJ model are presented in the left panel, shown as black data points. The corresponding log-log linear fit is shown as a solid black line and, for comparison, the trends found by Thomas et al. (2010) and Johansson et al. (2012)⁴ are plotted with dotted red and blue lines, respectively. As a test of robustness, the coloured rectangles represent the results for the default stacks used for the cosmic chronometer analysis (Figs. 6 and 8) that had an average redshift below 0.1. The running of all three parameters behaves as expected with respect to velocity dispersion: the more massive the (passive) elliptical galaxy, the older, the more metal-rich, and the more α-enhanced it is. However, there is a caveat regarding the metallicities for log σ ≳ 2.3: the posterior distributions saturate at the high end of the prior, indicating a preference for values higher than 0.26, the limit set as explained in Section 3.4. While the Z − σ fit omitted these data points, it should be noticed that relaxing this upper bound and allowing metallicities to extend up to 0.67 breaks the smooth increase in age with velocity dispersion. Indeed, due to age-metallicity degeneracy, the values of [Z/H] over 0.3 cause the ages to decrease and produce an unstable trend towards the massive end for certain cases. As a consequence, we consider it safer to retain the hard prior on the maximum value of the metallicity.

Fig. 11. Left panel: Age, metallicity, and α-enhancement trends with velocity dispersion at low redshift (0.05 < z < 0.10). The black data points represent the results from stacks of 400 galaxies. The dashed red and blue tendencies represent the Thomas et al. (2010) and Johansson et al. (2012) results (central values). The best log-log fit is plotted as a solid black line. We also depict the results from the ‘default’ stacks as coloured rectangles representing the different velocity dispersion groups from dark blue (100 < σ [km s−1] < 125) to orange (250 < σ [km s−1] < 320). Note that the colour code of the velocity dispersion groups corresponds to the one used in Figs. 6 and 8. Right panel: Same as the left panel but with the Knowles model. Red triangles are the results by Knowles et al. (2023) for stacks of massive ETGs, which were previously constructed by La Barbera et al. (2013) using the same model.

Our results differ quantitatively from those in the literature. In particular, we predict systematically lower ages than Johansson et al. (2012), with an approximately constant offset of ∼2 Gyr. We also observe lower ages compared to the Thomas et al. (2010) results in lower velocity dispersion data. Our results also indicate that galaxies tend to be more metal-rich and α-enhanced. There are a number of reasons for this quantitative discrepancy. First, the Lick index model used in Thomas et al. (2010) is based on the Thomas et al. (2003) and Thomas et al. (2004) models, the precursors of the TMJ model (used both in Johansson et al. 2012 and in this work) that were not flux-calibrated. Second, both Thomas et al. (2010) and Johansson et al. (2012) performed an analysis of the spectra of single galaxies from the SDSS that were only morphologically selected to be ETGs. On the other hand, we work with high ⟨S/N⟩ stacked spectra of galaxies and aim at a selection of the most passively evolving sources. Third, the choice of Lick indices and the fitting procedure are different, since we fix a priori a set of indices for all galaxies based on the quality of their calibration with respect to globular cluster data as indicated in Thomas et al. (2011). However, the qualitative physical picture remains unchanged and in line with the downsizing scenario.

The right panel of Fig. 11 shows the corresponding results using the Knowles model, depicted in black. For comparison, we plot the results from Knowles et al. (2023) as red triangles, obtained via the same model and Lick index choice we used, taking the set of ETG stacks from La Barbera et al. (2013). The age, metallicity and α-enhancement for the default stacks used in the cosmic chronometer analysis are also shown as coloured rectangles. The agreement among the three sets of results is remarkable and serves as a further consistency check that the stacks have been constructed in a reliable manner. With regard to the galactic parameters, the qualitative trends remain unchanged and are consistent with a downsizing scenario. However, the Knowles model predicts systematically higher ages than TMJ, with an approximately constant offset of ∼2 Gyr across all velocity dispersions. This translates directly into much smaller metallicities that reach even subsolar values for σ ≲ 160 km s⁻¹. Lastly, the sample is clearly α-enhanced, although the values are slightly smaller than those of the TMJ model across all masses.

The average ages, metallicities and α-enhancements obtained via a log-log linear fit of the above data can be regarded as the mean values of these quantities at the population level of massive quiescent galaxies. Even if the overall shape of the distributions for the single galaxy estimation is underwhelming, the mean of this distribution uncannily aligns with the results obtained using the more reliable stacking procedure. This is clearly seen in the upper-left panel of Fig. C.1, where the ages derived from Lick index fitting of single galaxies (plotted as density contours) are compared to the best fit for the stacking analysis with the TMJ model, depicted as a solid black line). While the median of single-galaxy ages (plotted as a thick blue line) seems to align with the values obtained from the stacks, a very large number of single massive galaxies are assigned very low or extremely high ages, as portrayed by the high density (darker blue contour) group (lobe) of individual galaxies in between 2 and 5 Gyr in the panel. This is also true for the metallicity and α-enhancement trends, which are shown in the other two panels of the figure. This is a further hint that the S/N of the spectra needs to be high enough to derive stellar population parameters in a reliable way.

However, we can exploit the distribution of the single-galaxy stellar population parameters to derive a mean dispersion at the population level for the age, metallicity and α-enhancement of massive quiescent galaxies at low redshift, that is, σ_{log t}, σ_[Z/H] and σ_[α/Fe]. Indeed, for every velocity dispersion bin, the standard deviation of the stellar population parameters with respect to their mean value can be computed. Via a Bootstrap resampling procedure, we can assign an uncertainty to the standard deviation and perform a simple linear fit with respect to velocity dispersion, as shown in the lower panels of Fig. C.1. As expected, the intrinsic dispersion in the ages (and thus in the formation times) decreases with velocity dispersion, in agreement with the popular view that the more massive the local quiescent galaxies, the more synchronous the formation of their progenitors was. The procedure is then repeated for the Knowles model, yielding the same qualitative behaviour, as shown in the right panels. However, the scatter in the ages for this model drops abruptly for large velocity dispersions due to the very high ages that it systematically assigns to the most massive galaxies. Through this approach, we provide an analytical recipe for the age, metallicity and α-enhancement distribution of local quiescent galaxies and, in particular, for the formation time of their progenitors. These scaling relations can be used as inputs, for example, in studies relying on the formation time distribution of the progenitors of local quiescent galaxies (Bosi et al. 2025). The results are summarised in the left column of Table 1.

5.2. An oscillating pattern in the index-redshift relations

When analysing the scatter of ages with redshift in Fig. 6, we revealed some synchronous oscillations at various redshift intervals. These effects seem to break the global tendency in a manner that natural statistical fluctuations cannot explain. The effect is most clearly seen for the TMJ model and around z ∼ 0.17, although other small maxima could be identified at z ∼ 0.08 and z ∼ 0.13. After ensuring that the fluctuations were independent of the S/N of the stacks, we decided to look into the index measurements, aiming to assess whether there were specific indices responsible for the oscillatory pattern.

Fig. 12 shows the age difference with respect to a Planck cosmology for all stacks in the three central velocity dispersion groups⁵ at redshifts 0.16 < z < 0.19 (where the main oscillation is found) as a function of several Lick index strengths; namely Hγ_F, Fe4383, Hβ and Fe5335. It is clear that Hγ_F and Fe4383 portray a strong correlation with respect to the age deviation from Planck that is not present in the other two indices. It should be noted that this does not mean that Hβ is not sensitive to stellar ages; in fact, it correlates very well when plotted over the whole sample. What the figure tells us is that certain indices are more responsible than others for the main fluctuation in the data irrespective of their intrinsic correlation with stellar ages.

Fig. 12. Difference between our estimated ages and the corresponding Planck Collaboration VI (2020) trends (as shown in the form of dashed black lines in Figure 6) with respect to the index strength of HγF (left panel) and Fe4383 (right panel). For the sake of clarity, only the intermediate velocity dispersion groups, with enough data in the interest region 0.15 < z < 0.21, are represented.

As a consequence, we decided to look more closely at the index measurements and found an oscillatory behaviour with respect to redshift for a number of them. A paradigmatic example of this finding is shown in Fig. 13, where the coloured dots are our measurements of the Hγ_F index for every galaxy stack in the four central velocity dispersion groups (again, the colour code for velocity dispersion is the same as in previous figures). As depicted in the bottom subpanel, a synchronous oscillation of ∼10% with respect to the smoothed global trend (depicted as a solid black line) is clearly visible for three velocity dispersion bins around z ∼ 0.17, where the main fluctuation feature is found for the ages. To determine if this could be related to the stacking procedure, we decided to measure the index strength for all single galaxies used in the stacks individually. The mean index value for every redshift is shown as a solid line, with the same colour code. As it is clearly observed, the single-galaxy trend follows perfectly that of the stacks, proving that the oscillatory phenomenon comes from a different nature. Once the stacking was discarded as the reason behind this effect, we could only turn to some flaw in the measurement of the index strengths. Fortunately, the SDSS archive offers measurements of Lick indices as well in their galSpecIndx table (Tremonti et al. 2004; Brinchmann et al. 2004). The mean for every redshift is shown with thin dashed lines, again with the same colours. The agreement is extremely perfect for z < 0.25, to the point that the solid and dashed lines overlap almost at every point with our measurements coming from both single galaxies and stacks.

Fig. 13. HγF trend with redshift. The measurements from the stacks are shown as data points, coloured by the corresponding velocity dispersion bin according to the colour code used in previous figures. The solid black lines represent the global trend via smoothed moving means, while the solid coloured lines are our measurements for single galaxies. The results from the SDSS database are depicted with coloured dashed lines.

Since the trait is particularly enhanced at 0.16 ≲ z ≲ 0.19, we can inspect our methodology to check for any particular divergence at that redshift. We have a strong non-linear behaviour in the SDSS instrumental resolution with wavelength at ∼5500 Å (rest frame) that is extended down to ∼4500 Å at z ∼ 0.3, as can be observed in Figure 4. One could suspect of the treatment of the instrumental resolution as a driver of (or element that enhances) the fluctuation. At the aforementioned redshift window, all spectral features falling within 4810 Å and 5350 Å could be affected by the discontinuity in the instrumental resolution of SDSS. The indices that are included in that wavelength range, partially or totally, are Hβ, Fe5015, Mg₁, Mg₂, Mg_b, Fe5270 and Fe5335. From these, Hβ and the bluest iron index were inspected in Fig. 12, but showed no correlation in the age-index strength plane. The ones that present the clearest correlation in that plane, Hγ_F and Fe4383, are safely found in the region where σ_IR is still ∼1.06 Å, so for the moment we can say that the (non-linear wavelength-dependent) instrumental resolution of SDSS is not the cause we are looking for.

The tests above show that there is an oscillatory trend in the Hγ_F index with respect to redshift that does not depend on galaxy mass or the S/N of the spectrum. It is neither related to the stacking procedure nor the measurement itself, pointing to something rooted inside the individual galaxy spectra. Moreover, when we extend this analysis to other indices, we see that this strange behaviour is not only typical of Hγ_F, nor is it in that index that the trait is most enhanced. The pattern (although with different strengths and at different redshifts) is found as well, for instance, in the Fe4383, Fe4531, Hδ_A and Hβ index strengths. All of these are included in the fitting set, yet many others not included likewise show oscillating patterns or fluctuations (Appendix D includes the oscillations found in several other indices apart from Hγ_F, shown in Fig. 13), as is the case for D4000_n and C₂4668. The latter is the most striking, with oscillations enhanced, periodic, and stable for a long redshift range, and more relevant than the index-velocity dispersion relation itself. The discontinuity feature suffers instead a significant drop at z ∼ 0.10 (Fig. D.3), breaking the expected stable and monotonous decrease in the feature with redshift. We call attention to the fact that such behaviour could undermine the good work of the D4000_n estimator at the referred redshift. Indeed, D4000_n estimations of H(z) at z < 0.4 are relatively scarce (Fig. 9), and Moresco et al. (2011) had already limited the work with this feature at z > 0.15, where it seems to gain stability.

Since the factors within our control do not seem to explain the observed oscillations, we believe that the only possibilities should be related to the SDSS data treatment. In particular, either the issue reflects a problem with the interpolation in the flux calibration procedure or is associated with atmospheric emission. As explained in Stoughton et al. (2002), these effects are taken into account in the SDSS data reduction pipeline by modelling hot stars and subtracting the empty sky. Regarding the former, the instrument response function is estimated using the observed spectra of hot young stars, whose spectra are well understood and calibrated. These spectra feature very deep Balmer lines, and modelling the continuum of these stars requires interpolating the flux at their wavelength region. This can introduce a bias in the response function at the rest-frame position of Balmer lines (since calibration stars are local) and, as a consequence, in the measurement of Lick absorption features that fall at the rest-frame wavelength of these Balmer lines. We take the Hγ_F index as an example. As shown in Fig. 13, this index displays a strong oscillation between z ∼ 0.16 and z ∼ 0.19. Taking into account the definition of both the central and pseudocontinuum regions of this index, a quick computation shows that the rest-frame wavelength corresponding to this index for these redshifts lies well below 3800 Å, far from the Hδ Balmer index, the only one that could potentially be responsible for a bias in the measurement. While it is true that other smaller fluctuations in Hγ_F, such as the minimum around z = 0.08, might fall at the right rest-frame wavelength for Hδ, we checked that the interpolation issue could not account for the vast majority of oscillations in all indices. The second possibility would be related to telluric absorption or emission lines not correctly accounted for in the SDSS data reduction pipeline. Trying the wealth of existing atmospheric lines that could be responsible for the oscillating patterns is beyond the scope of this work, but preliminary tests seem to show that they cannot contribute to the bulk of this anomaly. Further work will address these findings.

6. Conclusions

In this study, we employed a cosmographic approach in the context of cosmic chronometers to derive a continuous measurement of the Hubble parameter up to z ∼ 0.4 via the fitting of Lick index absorption lines from massive and passively evolving galaxies. To this end, we carefully constructed a sample of such objects from the SDSS Legacy data – resulting in a sample of 90 396 objects meeting the expected physical requirements for these galaxies – in agreement with other samples in the literature. The spectra are reliably stacked in bins of both velocity dispersion and redshift to ensure a similar and high enough S/N. After carefully considering instrumental resolution and velocity dispersion effects, Lick absorption line indices were measured and fitted with two stellar population models, namely TMJ (Thomas et al. 2011) and Knowles (Knowles et al. 2023), in order to derive stellar population parameters (t, [Z/H], [α/Fe]). We emphasise the importance of using stacked spectra for the analysis, as demonstrated by the stellar population parameter distributions obtained from the fit of single galaxy Lick indices, which feature a much less reliable and stable behaviour.

Results cover the cosmographic fit of the stellar ages, t(z; H₀, q₀, j₀), through the integration of the inverse of the Hubble parameter, expanded in the Taylor series of z/(1 + z) up to order 2. The credible intervals given for the cosmographic parameters {H₀, q₀, j₀} include only the statistical uncertainty associated with the scatter of the data and error propagation throughout the data treatment process. A brief summary of the systematic uncertainties due to our manipulation of data is included and applied directly to the final H(z) estimation. Namely, we tested the effect of the choice of the S/N level for the stacked spectra and the exclusion of stacks that average galaxies extended in redshift with σ_z greater than the limit we imposed (i.e. the baseline is σ_z < 0.005). We observed that these account for an impact similar to twice the statistical uncertainty. Other sources of systematic uncertainty related to the SPS models are not analysed in this work. However, we included an additional margin in the estimation of H(z) under the assumption that these could represent a similar value to the those we accounted for. This is based on recent studies that assessed the effect of the SPS model and the Lick index combination on the H(z) readout.

The TMJ model provides ΛCDM-aligned posteriors for our set of cosmographic parameters when the condition dH/dz > 0 is included. This yields a model-independent measurement of the Hubble constant, $H_0=70.0_{-7.6}^{+4.1}$ (stat.) km s⁻¹ Mpc⁻¹, that covers both late-time and early-time probes results. The deceleration parameter is skewed towards the limit imposed by the condition on the derivative of H, which translates into q(z) >  − 1; in particular q₀ > −1. The removal of this condition reveals the compatibility of the TMJ-modelled ages with combinations of cosmographic parameters in q₀ < −1 and j₀ > 50, which in those cases requires H₀ > 100. We gave a cosmological readout of our results by introducing a w₀CDM model in which our estimated q₀ and j₀ (baseline) can be translated into Ω_m and w₀. By requiring that 0 < Ω_m(q₀, j₀) < 1, we find that the condition dH/dz > 0 only removes some phantom solutions from the possible results, while keeping H₀ virtually unchanged. The associated dark energy equation of the state parameter, w₀, is compatible with a cosmological constant, with some but not all phantom states allowed.

Our final result consists of drawing H(z) using the MCMC chains from the baseline case (TMJ model with dH/dz > 0). The agreement with previous punctual measurements is good, and we conclude that our methodology is competitive with the analyses based on FSF (whose uncertainties are somewhat greater than ours) and covers a region (z < 0.15) where D4000_n measurements from SDSS proved to be problematic. We associate the latter trait with the unveiled pattern of oscillations and fluctuations in the index-redshift relations for several Lick indices and other features, in particular with the non-monotonicity in D4000_n and D4000.

Regarding the Knowles model, the distribution of stellar ages on the t − z plane posed some problems for the cosmographic fit. Uncertainties in the derived ages were smaller by a factor of ∼7 in comparison to the TMJ-modelled ages and were therefore much smaller than the intrinsic scatter in the data, which was also slightly higher than in the TMJ case. These issues were partially mitigated by the introduction of additional f-factor parameters in the likelihood. In our assessment, the overall age distribution remains less robust than that produced by the TMJ model. Hence, we advise some caution when using the Knowles-modelled ages for cosmographic applications, at least until the model has been more thoroughly examined. This includes evaluating its various assumptions, such as the choice of IMF, which is expected to contribute only marginally (Moresco et al. 2020); exploring alternative spectral diagnostics beyond those used in our fit; or even considering different methodologies other than spectral-feature fitting, taking advantage of the model being presented as full synthetic spectra for the given combinations of age, metallicity, and α-enhancement.

As a by-product, we gave relations for the evolution of stellar population parameters with respect to velocity dispersion for the stellar population models. As expected, these scaling relations are in agreement with the downsizing scenario, although with quantitative differences with respect to results in the literature. Beyond these, we also offered a quantitative description for the dispersion relations of these parameters.

We conclude our study by presenting the finding of an unexpected oscillating pattern in the distribution of index strengths with respect to redshift for a number of Lick indices, particularly those used in the cosmographic fits. These were found when trying to understand the fluctuating behaviour of the age-redshift dispersion. We observe that some index strength oscillations were more correlated with age oscillations than others, namely the Balmer Hγ feature (both _A and _F definitions) and some iron indices (in particular Fe4383), which are the main drivers of the oscillation around z ∼ 0.17. We investigated their potential origin and discarded a statistical nature due to the synchronised behaviour of the oscillations across velocity dispersion bins. A methodological origin was also ruled out, as fluctuations are found both in the stacked and single-galaxy spectra and are also present in the Lick index measurements provided by the SDSS (Tremonti et al. 2004; Brinchmann et al. 2004). Furthermore, flux calibration issues related to the interpolation of the continuum in the SDSS data reduction pipeline do not seem to account for such strong features, nor does a preliminary analysis of telluric absorption lines. To our knowledge, such oscillations have never been reported before and will be further studied in future work. In particular, the presence or absence of such oscillations in spectroscopic data from the GAMA survey will help inform us about the role of the SDSS spectrograph and data reduction in their origin.

The big picture for future works is that SPS models will remain a critical bottleneck and source of systematic uncertainty. Deep spectroscopic surveys (e.g. DESI and Euclid), as well as fully exploiting BOSS and eBOSS data, may provide the necessary leverage at z > 0.5. At higher redshifts, the parameters of the cosmographic fit may be more sensitive to kinematical parameters, and a custom H(z) expansion in the vicinity of a non-zero redshift is even possible, giving us the counterparts of our set of cosmographic parameters at higher redshift.

Acknowledgments

This work was partially funded from the projects: “Data Science methods for MultiMessenger Astrophysics & Multi-Survey Cosmology” funded by the Italian Ministry of University and Research, Programmazione triennale 2021/2023 (DM n.2503 dd. 9 December 2019), Programma Congiunto Scuole; EU H2020-MSCA-ITN-2019 n. 860744 BiD4BESt: Big Data applications for black hole Evolution STudies; Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union – NextGenerationEU – and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2 within the activities of Spoke 3 (Astrophysics and Cosmos Observations); European Union – NextGenerationEU under the PRIN MUR 2022 project n. 20224JR28W “Charting unexplored avenues in Dark Matter”; INAF Large Grant 2022 funding scheme with the project “MeerKAT and LOFAR Team up: a Unique Radio Window on Galaxy/AGN co-Evolution”; INAF GO-GTO Normal 2023 funding scheme with the project “Serendipitous H-ATLAS-fields Observations of Radio Extragalactic Sources (SHORES)”. We further want to mention the helpful contribution and counsel of Alexandre Vazdekis, Nicola Borghi, Michele Moresco and Christy A. Tremonti.

References

Appendix A: Correction of the velocity dispersion effect

The process to transform the raw measurements of Lick indices (or in general spectral features) to data ready for the cosmographic fit was outlined in Sec. 3.3. In this Appendix we provide with the functional form of the corrections that need to be applied in order to take data at σ [km s⁻¹] velocity dispersion and instrumental spectral resolution of 1.5Å (root mean squared wavelength resolution) to zero velocity dispersion and σ_IR = 1.06Å.

Figure A.3 shows the correction functions for all 25 basic Lick indices. The choice of colour depends on whether the feature has been used in the fit with each of the models. All the corrections are computed as 3rd order polynomials in σ [km s⁻¹], with central values for the parameters a_i, C_I(σ, σ_IR) = ∑_ia_iσⁱ, shown in table A.1. For the atomic indices the correction takes the form C(σ, σ_IR) = I(0, σ_MILES)/I(σ, σ_IR), while for the molecular indices (CN₁, CN₂, Mg₁, Mg₂, TiO₁, TiO₂) it is computed as C(σ, σ_IR) = I(0, σ_MILES)−I(σ, σ_IR). σ_IR = 1.5Å and σ_MILES = 1.06Å. The uncertainty associated with the correction δC is small in all cases, due to the little scatter (over the sample of SDSS giant stars) and good approximation of the average computed I₀/I(σ) to a third order polynomial. Both Hδ features show a thicker line, representing the 1σ exceeding the thickness of the line indicating the fiducial value.

In this paper, we presented two different ways to address the correction of the velocity dispersion effect imprinted in galaxy spectra. The common way is the aforementioned one, through the correction functions, and applied a posteriori. The other one, in which we presented on the Knowles model, consists in correcting the SPS model to make velocity dispersion become a fine tuning configuration. However, this can only be applied to Lick index models that come from direct integration over spectra, which is not the case of the TMJ model. Even if, as discussed in the text (Sec. 4), the ages derived with the Knowles model have to be taken with caution, we analyse the effect that these two different approaches have on the age distribution. To visualise the difference we present in Fig. A.1 the rate Δt/t_corr between the age difference derived using both methodologies Δt ≡ t_in-model − t_corr, where t_in-model is the age derived using the model adapted to different velocity dispersions and t_corr is the age derived using the corrections.

Fig. A.1. Effect on the modelled age (Knowles model) of varying the methodology to account for the velocity dispersion effect. Colours are chosen to represent the velocity dispersion groups in Fig. 8.

Fig. A.2. Comparison of the uncertainties derived from Knowles-modelled age resulting from the two possible treatments of the velocity dispersion effect.

Fig. A.3. Velocity dispersion correction functions C(σ, σIR) for the 25 Lick indices. Yellow lines were chosen to represent the indices that were used for the fit with the two SPS models, while blue and red indicate those that were used only with one of them (Knowles and TMJ, respectively). Grey tendencies represent the corrections for the indices that were not used with any of the models.

The global distribution of Δt/t_corr is Gaussian and with a mean close to 0. Indeed, P₁₆(Δt/t_corr) =  − 0.9%, Med(Δt/t_corr) = 0.2%, and P₈₄(Δt/t_corr) = 1.9%. We conclude that even if the ages from the two methodologies seem highly compatible, some systematic uncertainty can be due to the treatment of the velocity dispersion effect. The uncertainty of the correction functions is already included as a part of the statistical uncertainty in the moment the galactic parameters fit is performed over the set of Lick index measurements. In general, this cannot be accounted for when the model is modified to include velocity dispersion as a fine-tuning parameter. Interested in assessing the impact of the different treatments of the velocity dispersion effect onto the uncertainties derived for the age estimation, we performed a fit of the corrected Lick indices with the Knowles model at zero velocity dispersion. The comparison of the difference in the uncertainties in age coming from both fits, σ_{t, in-model} − σ_t, corr, is presented in Fig. A.2.

The distributions for all the velocity dispersion groups, on the right panel in Fig. A.2, shows a perfect compatibility of the uncertainties in the two cases. We see that the choice of methodology to treat the velocity dispersion effect is irrelevant to the age uncertainties, which supports the idea that the nature of each model is the root of the differences observed between Figs. 6 and 8.

Appendix B: Cosmographic fits

In Sec. 4 we commented on the results of several cosmographic fits. Of those obtained from the TMJ-modelled ages, we showed the posterior distributions for the baseline (Fig. 6) and the w₀CDM readout of our results (Fig. 7). Now we show some of the remaining corner plots that we did not present in the main text. In Fig. B.1 we unveil the posterior distributions for the case without the dH/dz > 0 limitation (top) and the case in which we omit the data from the region of the strong oscillation in 0.15 < z < 0.20 (bottom). In the former case, we shadowed the region with q₀ < −1, indicating all the parameter space that is limited by the condition on the derivative of H in the baseline fit. For the later, we highlight how the exclusion of the data affected by the greater oscillation in the index-redshift trend does not particularly improve the overall fit.

Fig. B.1. Sampling of the parameter space relative to the cosmographic parameters. On top, the general case in which the dH/dz > 0 condition was removed; for clarity, not only is q0 = −1 shown but also the entire q0 < −1 region. On the bottom, the baseline case is shown with the exclusion of the data in the range 0.15 < z < 0.20.

Finally, we show an example of the full sampling of the parameter space for the baseline scenario, including the ages at zero redshift for each velocity dispersion group, t_0, v, in Fig. B.2. As commented in the main text, the posterior distributions for the present-time ages are well defined within the prior we set, with typical 1σ of the order of 0.3 Gyr.

Fig. B.2. Full posterior probability distributions of the parameters included in the fit from Fig. 6.

Appendix C: Scaling relations

This appendix is included to show the distribution of single galaxies in the stellar parameter - velocity dispersion planes, as in Fig. C.1, where in blue we presented the results for the TMJ model and in light red for the Knowles model. In black we plot the tendency obtained from the custom stacks that are mapped in Fig. 11, showing the good agreement between the moving mean of the single galaxies (darker blue / red line in each case) with the fit coming from the stacks. In the lower panels, the bootstrap methodology explained in the main text give us a relation for the dispersion of the parameters and the velocity dispersion.

Fig. C.1. Density contours of the distribution age (upper panels), metallicity (middle panels), and α-enhancement (lower panels) for single galaxies. The thick blue line represents the median, while the black line is the best log-log fit from the stacking procedure. The bottom subpanels of each panel show the scatter in the parameter distribution. Results for the TMJ and Knowles models are shown in blue (left panels) and pink (right panels).

Appendix D: Index strength oscillations

This appendix is included to show the behaviour of several indices with redshift. We included the remaining Balmer indices (Fig. D.1), as they present oscillations that represent deviations from the average value of the index of up to a 10%. However, the Hδ_F index does not exhibit a clear coherence between the different velocity dispersion groups. Two iron indices (Fe4383, Fe4531) are presented (Fig. D.2) as the oscillations are clearly synchronous among the velocity dispersion groups. Even if the value of the fluctuation is relatively small relative to the index strength, the expected smooth behaviour is broken clearly.

In Fig. D.3 we present the D4000 and D4000_n spectral features, showing a clear non-monotonous behaviour. Particularly, from z ∼ 0.10 to z ∼ 0.15, the D4000_n grows with redshift, opposing the expected soft decreasing trend. Then, in Fig. D.4, the C₂4668 and TiO₁ lines are shown to exhibit the case of an extreme oscillatory behaviour. In all cases the SDSS measurements (thin dashed lines), stacked data (dots) and single galaxies average (thick line) agree, backing our thesis that these are not spurious fluctuations arising from our treatment of the data.

Fig. D.1. Index-redshift relation for Balmer lines. On top HδA (left) and HδF (right). Below HγA (left) and Hβ. The structure of each panel is equivalent to Fig. 13.

Fig. D.2. Index-redshift relation for Iron lines Fe4383 (left) and Fe4531 (right). The structure of each panel is equivalent to Fig. 13.

Fig. D.3. Index-redshift relation for the discontinuity features, D4000 (left) and D4000n (right). The structure of each panel is equivalent toFig. 13.

Fig. D.4. Index-redshift relation for indices not included in the stellar parameter fits. C24668 and TiO1 are shown on the left and right, respectively. The structure of each panel is equivalent to Fig. 13.

All Tables

All Figures

	Fig. 1. Distribution of D4000n and three Balmer-line indices over the parent (grey), query (blue), and final (maroon) galaxy samples.
In the text

	Fig. 2. Distribution of sources in the (u − r)0 vs. (r − z)0 colour-colour plane. The parent sample is shown in grey, while the increasingly restrictive subsamples (adding velocity dispersion, query and spectroscopic information) are shown in light blue, dark blue and maroon, respectively. The dashed black line defines the passive galaxy region (Holden et al. 2012).
In the text

	Fig. 3. Properties of the ancillary stacks. Left panel: Distribution of the number of galaxies included in the ancillary stacks defined by a linear constant spacing in both σ and z. Right panel: Relation between the S/N and the number of galaxies in each ancillary stack as a function of redshift.
In the text

	Fig. 4. RMS wavelength resolution (σR) of SDSS stacked spectra as a function of rest-frame wavelength for different redshifts.
In the text

	Fig. 5. Age, metallicity, and α-enhancement distributions for the TMJ (left panels) and Knowles (right panels) models, interpolated to a grid of 50 points in each parameter within the original ranges. Results are shown for the Lick index fitted to both single galaxies (in blue) and stacks (in yellow) of the final cosmic chronometer sample.
In the text

Fig. 6. Joint cosmographic fit of t(z) with ages estimated using the TMJ model. On the left, the corner-plot for the posterior probability distributions of H0, q0, and j0 is shown, displaying 1σ and 2σ contour levels. On the right, the t − z dispersion is shown. The empty circles in strong colour represent the data points we used for the fit, while the faded crosses represent the points, resulting from of a set of galaxies with average dispersion in redshift greater than 0.01. The shaded regions in the same colour as the data points are the 1σ of the sampling of the posterior distribution coming from the parameters on the left. The zero-ages of the velocity dispersion groups are not shown here, but we refer to the corner plot in Fig. B.2; however, they are used for the sampling and to set a zero-age for the Planck Collaboration VI (2020) evolution of t(z), shown as dashed black lines. The dashed q0 = −1 lines in the corner plots for the planes q0 − H0 and j0 − q0 represent the natural limit imposed by dH/dz > 0. This is the baseline t − z relation referred to in this work.

In the text

	Fig. 7. Posterior probability distribution for H0, Ωm, and w0. The sampled parameters were q0 and j0, then converted to w0CDM parameters through Equations (18) and (19). The sampling was run under the ΛCDM-model condition of 0 < Ωm(q0, j0) < 1. The brown hatched region in the Ωm − w0 plane represents the region allowed by dH/dz > 0 ⇔ w0 > 1/(Ωm − 1).
In the text

Fig. 8. t − z relations for the ages estimated using the Knowles model. The structure of this figure is the same as that in Fig. 6. The grey contours represent the case with fv factors included, to account for the possible underestimation of age uncertainties. The t − z samplings on the right panel refer to this case. The contours in purple represent the case where these factors were not taken into account. The sampling of the associated posterior for t − z is not represented on the right panel for this case.

In the text

Fig. 9. Reconstruction of the Hubble parameter. Our sampling is shown in shades of blue, corresponding to the posterior distribution obtained for H0, q0, and j0 from Fig. 6. The width of the tendency represents the range between the 16th and 84th percentiles of our posterior sampling. The opacity of the colour shade represents the reliability of the H(z) measurement, which we based on the amount of sources (individual galaxies) at each z. The dashed black line represents the H(z) evolution according to Planck Collaboration VI (2020). The plot is extended up to z = 0.9 in order to show the result at the lowest redshift obtained through the Lick indices fit (dark red), by Borghi et al. (2022b). Other measurements from the available literature are represented in yellow (obtained through FSF) and orange (obtained through the measurement of D4000 of its narrow counterpart D4000n).

In the text

Fig. 10. Systematic uncertainties arising from the methodology directly accounted for in this work. They are presented in units of 1σ of the statistical uncertainty σstat at each redshift. The individual contributions of the level of S/N selected in the stacking procedure and the limit to the dispersion of sources in redshift in each stack are presented in yellow and dark orange, respectively. In dark blue, we present the overall systematic uncertainty.

In the text

Fig. 11. Left panel: Age, metallicity, and α-enhancement trends with velocity dispersion at low redshift (0.05 < z < 0.10). The black data points represent the results from stacks of 400 galaxies. The dashed red and blue tendencies represent the Thomas et al. (2010) and Johansson et al. (2012) results (central values). The best log-log fit is plotted as a solid black line. We also depict the results from the ‘default’ stacks as coloured rectangles representing the different velocity dispersion groups from dark blue (100 < σ [km s−1] < 125) to orange (250 < σ [km s−1] < 320). Note that the colour code of the velocity dispersion groups corresponds to the one used in Figs. 6 and 8. Right panel: Same as the left panel but with the Knowles model. Red triangles are the results by Knowles et al. (2023) for stacks of massive ETGs, which were previously constructed by La Barbera et al. (2013) using the same model.

In the text

	Fig. 12. Difference between our estimated ages and the corresponding Planck Collaboration VI (2020) trends (as shown in the form of dashed black lines in Figure 6) with respect to the index strength of HγF (left panel) and Fe4383 (right panel). For the sake of clarity, only the intermediate velocity dispersion groups, with enough data in the interest region 0.15 < z < 0.21, are represented.
In the text

Fig. 13. HγF trend with redshift. The measurements from the stacks are shown as data points, coloured by the corresponding velocity dispersion bin according to the colour code used in previous figures. The solid black lines represent the global trend via smoothed moving means, while the solid coloured lines are our measurements for single galaxies. The results from the SDSS database are depicted with coloured dashed lines.

In the text

	Fig. A.1. Effect on the modelled age (Knowles model) of varying the methodology to account for the velocity dispersion effect. Colours are chosen to represent the velocity dispersion groups in Fig. 8.
In the text

	Fig. A.2. Comparison of the uncertainties derived from Knowles-modelled age resulting from the two possible treatments of the velocity dispersion effect.
In the text

	Fig. A.3. Velocity dispersion correction functions C(σ, σIR) for the 25 Lick indices. Yellow lines were chosen to represent the indices that were used for the fit with the two SPS models, while blue and red indicate those that were used only with one of them (Knowles and TMJ, respectively). Grey tendencies represent the corrections for the indices that were not used with any of the models.
In the text

	Fig. B.1. Sampling of the parameter space relative to the cosmographic parameters. On top, the general case in which the dH/dz > 0 condition was removed; for clarity, not only is q0 = −1 shown but also the entire q0 < −1 region. On the bottom, the baseline case is shown with the exclusion of the data in the range 0.15 < z < 0.20.
In the text

	Fig. B.2. Full posterior probability distributions of the parameters included in the fit from Fig. 6.
In the text

Fig. C.1. Density contours of the distribution age (upper panels), metallicity (middle panels), and α-enhancement (lower panels) for single galaxies. The thick blue line represents the median, while the black line is the best log-log fit from the stacking procedure. The bottom subpanels of each panel show the scatter in the parameter distribution. Results for the TMJ and Knowles models are shown in blue (left panels) and pink (right panels).

In the text

	Fig. D.1. Index-redshift relation for Balmer lines. On top HδA (left) and HδF (right). Below HγA (left) and Hβ. The structure of each panel is equivalent to Fig. 13.
In the text

	Fig. D.2. Index-redshift relation for Iron lines Fe4383 (left) and Fe4531 (right). The structure of each panel is equivalent to Fig. 13.
In the text

	Fig. D.3. Index-redshift relation for the discontinuity features, D4000 (left) and D4000n (right). The structure of each panel is equivalent toFig. 13.
In the text

	Fig. D.4. Index-redshift relation for indices not included in the stellar parameter fits. C24668 and TiO1 are shown on the left and right, respectively. The structure of each panel is equivalent to Fig. 13.
In the text