© The Authors 2025

1. Introduction

One of the most intensely debated topics of modern astrophysics is the formation and quenching of massive quiescent galaxies. Observations and simulations suggest that quiescent galaxies can undergo different formation paths, which implies a variety of star formation histories (SFHs; e.g., Bevacqua et al. 2024) and quenching mechanisms (e.g., Tacchella et al. 2022), depending on the properties of the gas from which they formed and the subsequent mass growth and assembly, which can vary significantly depending on the local cosmological conditions and environment (e.g., De Lucia et al. 2006).

Many processes can cause the halt of star formation in galaxies (e.g., Man & Belli 2018 and references therein), and there is no current agreement on the dominant quenching mechanism, despite the efforts in the past decades (e.g., (Thomas et al. 2005; McDermid et al. 2015; Newman et al. 2018; Rowlands et al. 2018; Belli et al. 2019; Bluck et al. 2022; Piotrowska et al. 2022; De Lucia et al. 2024)). Quenching mechanisms can be broadly distinguished into rapid mechanisms (<1 Gyr), for instance, ejective active galactic nucleus (AGN) feedback and ram-pressure stripping, and slow mechanisms (>1 Gyr), for example, virial shocks and thermal AGN feedback. Simulations (e.g., Rodriguez-Gomez et al. 2019) and observations (e.g., Tacchella et al. 2022) both found a broad range of quenching timescales for galaxies in the local Universe, which implies that nearby quiescent galaxies experienced different quenching mechanisms. Alternatively, galaxies can simply run out of gas as a consequence of highly efficient star formation (e.g., Vazdekis et al. 1996, 1997; Salvador-Rusiñol et al. 2020).

At high redshifts, the processes that cause the galaxy quenching are typically shorter because the Universe is younger. For instance, massive quiescent galaxies found at z≳3 (e.g., Forrest et al. 2020a, b; Saracco et al. 2020; Glazebrook et al. 2024) require very high star formation efficiencies and short quenching times. As a consequence of an efficient starburst followed by rapid quenching, it is expected that high-redshift quiescent galaxies resemble local post-starburst galaxies (PSBs), namely very young galaxies with negligible ongoing star formation. Observations seem to support this. For instance, using a sample of nine galaxies at z∼3, D’Eugenio et al. (2020) showed that the typical near-UV/optical spectrum of high-z quiescent galaxies is characterized by a strong Balmer break that is larger than the 4000 Å   break and corresponds to ages of a few hundred million years, typical of PSBs. The ignition of the starburst and the main cause of the rapid quenching of these galaxies are still unclear, however.

As a result of instrumental limitations, only a limited number of (massive) quiescent galaxies have been observed at high redshift so far. The advent of the James Webb Space Telescope (JWST) has opened the possibility of performing large galaxy surveys at high redshift, which allows us to study the properties of a large number of high-redshift quiescent galaxies at different masses. Recent observations with JWST have revealed a large number of quiescent galaxies at very high redshift, up to z = 7.3 (Looser et al. 2024; Weibel et al. 2025), that span a wide range of masses (Carnall et al. 2023, 2024; Labbé et al. 2023; Looser et al. 2025, 2024; de Graaff et al. 2025a, b; Slob et al. 2024; Wang et al. 2024; Glazebrook et al. 2024; Alberts et al. 2024). Several low-mass (<10¹⁰ M_⊙) quiescent galaxies at high redshift exhibit bursty SFHs, that is, they are characterized by episodic events of short star formation followed by periods of prolonged quiescence. This suggests that they are just temporary (or mini-) quenched (Looser et al. 2025, 2024; Strait et al. 2023; Dome et al. 2024). For these galaxies, quenching via environmental interactions or ejective stellar feedback is currently favored over other mechanisms, such as supernova feedback (Looser et al. 2025; Sandles et al. 2023; Strait et al. 2023; Gelli et al. 2023, 2024; Carniani et al. 2024a; Zhang et al. 2024; Cutler et al. 2024; Hamadouche et al. 2024). For massive galaxies, many studies have reported a high incidence of strong AGN activity in recently quenched galaxies with PSB-like spectra associated with massive, high-velocity multiphase gas outflows (Belli et al. 2024; Davies et al. 2024; Park et al. 2024; Scholtz et al. 2024; Wu 2025; D’Eugenio et al. 2024; see also Tremonti et al. 2007; Baron et al. 2017; Förster Schreiber et al. 2019; Maltby et al. 2019), suggesting that ejective AGN feedback plays a major role in quenching massive galaxies at high redshift, which agrees with some simulations (Pontzen et al. 2017; Davis et al. 2019; Zheng et al. 2020). The population of high-redshift quiescent galaxies generally exhibits a rich variety of SFHs (e.g., Carnall et al. 2024; Park et al. 2024; Nanayakkara et al. 2025), however, which agrees with studies at lower redshift (e.g., Thomas et al. 2005; Panter et al. 2007; McDermid et al. 2015; Pacifici et al. 2016; Tacchella et al. 2022; Bevacqua et al. 2024).

A major issue related to estimating the properties of these high-redshift galaxies, including their mass and SFH, is the effect of different stellar population models on the inferred properties. It is well known that stellar population models struggle to reproduce the near-infrared (NIR) spectrum of galaxies (e.g., Eminian et al. 2008; Riffel et al. 2015, 2019; Eftekhari et al. 2021). For quiescent galaxies, the characterization of the NIR spectrum is crucial because a significant fraction of the stellar light is emitted in this wavelength region. In the past decade, much effort has been made in this regard, and current models can now rely on better and more complete stellar libraries, as well as on the characterization of reliable NIR spectral indices (Riffel et al. 2019; François et al. 2019; Morelli et al. 2020; Eftekhari et al. 2021). Despite the improvements, the (few) studies that were based on local galaxies showed that models still struggle to reproduce the NIR spectral indices (e.g., Eftekhari et al. 2022; La Barbera et al. 2024).

One of the main sources of uncertainty for the characterization of the NIR spectrum of galaxies is the contribution of thermally pulsating (TP-)AGB stars. TP-AGB stars are intermediate-mass stars (1−10 M_⊙) whose luminosity is maintained by a double-shell burning regime that determines the pulsation. Based on the fuel-consumption theorem (Renzini & Voli 1981), it is estimated that TP-AGB stars dominate the galaxy spectrum at ages between 0.2 and 2 Gyr and peak at about 1 Gyr, when they are expected to contribute up to ∼80% in the K band and up to ∼40% to the bolometric luminosity (Maraston 1998, 2005; Marigo et al. 2008). Including them in the model recipes then significantly affects the inferred galaxy properties (Pastorelli et al. 2020). For example, Maraston et al. (2006) showed that models with a heavy contribution of TP-AGBs predict systematically younger ages and lower masses by up to ∼60% compared to models with little contribution.

After decades, the disagreement persists for several reasons. Theoretically, the modeling of this evolutionary phase is hampered by our poor understanding of many key processes of stellar evolution (e.g., mass loss, dredge-up, or varying opacities). As a result, different models that adopt different recipes lead to different results. Observationally, TP-AGB stars are rare objects because their lifetimes are short, and only small samples of TP-AGB stars are observed in a few local globular clusters and the Magellanic clouds. Further, it has been shown that estimates of the TP-AGB contribution to the integrated light of galaxies derived from these samples are biased toward high values (Girardi et al. 2013). Additionally, TP-AGB stars are expected to contribute little to the optical spectrum of galaxies, but they are prominent in the NIR, where observations before the JWST era were severely limited.

The best galaxy candidates for probing the contribution of TP-AGB stars are high-redshift PSBs because they are young systems without ongoing star formation and no strong contamination from cosmological evolution. Observations have provided conflicting results so far, however, which led to a vivid debate of the actual contribution of TP-AGB stars. For instance, Kriek et al. (2010) found a low contribution of TP-AGB stars in PSBs observed at 0.7<z<2, while Capozzi et al. (2016) found a strong contribution for a sample of passive galaxies at 1.3<z<2.7, which agreed with recent results by Lu et al. (2025), who used prism observations with JWST of three PSBs. At lower redshifts, Conroy & Gunn (2010) found that models with a low contribution of TP-AGB stars reproduce the NIR colors of nearby PSBs well, while models with a heavy contribution do not. This result was further corroborated by Zibetti et al. (2013) based on results with spectroscopic data. On the other hand, MacArthur et al. (2010) showed that the bulges of local spiral galaxies require a strong contribution. In between these results, Melnick & De Propris (2014) showed that models with a moderate contribution of TP-AGB stars can also reproduce the NIR photometry of local PSBs. In brief, the results in the literature are mixed.

One problem with previous results in the literature is that they are based on broadband photometry, optical spectroscopy, or low-resolution NIR spectra, while good-quality NIR spectroscopic data are required for an accurate investigation of the spectral features. For instance, Riffel et al. (2015) showed that although nearby galaxies can exhibit strong features characteristic of TP-AGB stars (e.g., the CO and CN molecular bands), well-calibrated models based on empirical stellar libraries reproduce the observations better than models with a significant contribution of TP-AGB stars (see also Dahmer-Hahn et al. 2018). More recently, Eftekhari et al. (2022) showed that the mismatch between the observed CO indices in local massive early-type galaxies (ETGs) and the predictions of advanced stellar population models is strong, although local massive quiescent galaxies might not be an optimal sample for probing the TP-AGB phase because they are very old while the contribution of TP-AGB stars peaks at young ages. The JWST now enables us to obtain spectra of high-redshift young quiescent galaxies with a sufficiently high resolution to probe the contribution of TP-AGB stars directly based on the spectral indices.

We perform a spectral analysis of a quiescent PSB observed at z∼2. In section 2 we describe the JWST data and the libraries we used to analyze it. In section 3 we estimate the physical and stellar population properties and reconstruct the SFH of the galaxy from the rest-frame optical and IR spectra. In section 4 we perform a detailed analysis of the optical and NIR spectral indices. In Sect. 5 we discuss our results, and we summarize them in Section 6. We adopt a flat ΛCDM cosmology with H₀ = 70 km s⁻¹ and Ω_m = 0.3.

2. Data and templates

2.1. Observations and data reduction issues

We selected our target galaxy from the JWST Advanced Deep Extragalactic Survey (JADES) (Eisenstein et al. 2023; Bunker et al. 2024). For the spectroscopic data, we considered the data reduction of the third data release (D’Eugenio et al. 2025, hereafter DR3) of JADES. The spectroscopic data consists of observations with the two medium resolution (R ∼ 1000) gratings, G140M/F070LP and G395M/F290LP. For the photometric data, we considered the public JADES catalog (Hainline et al. 2024). To select our target, we initially considered only galaxies having spectroscopic observations with both the prism and gratings observed at z>1. Then, we used the EAZY code (Brammer et al. 2008) to fit the photometry and estimate the rest-frame U−V and V−J colors, from which we selected quiescent galaxies using Eq. (4) of Williams et al. (2009). Finally, we considered galaxies having signal-to-noise ratio (S/N) > 15 in all spectra. The final sample consists of two galaxies having IDs 138717 (hereafter, J-138717) and 199773. In this paper, we only study the former; the latter is presented in D’Eugenio et al. (2025). The photometric data of J-138717 consist of three NIRCam observations with filters F182M, F210M, and F444W, and four HST observations with filters F606W, F775W, F814W, F850LP. We considered the Kron photometry convolved by the point-spread function to the F444W spatial resolution (see Hainline et al. 2024 for details). The photometric catalog also provides the photometric redshift, z_phot = 2.16, and the half-light radius, ≈0.13 arcsec, corresponding to about 1.1 kpc. The total integration time for the G140M, G395M, and the prism spectra are 9.2 hr, 35.6 hr, and 40.5 hr, resulting in a median S/N of about 19 Å⁻¹, 54 Å⁻¹, and 59 Å⁻¹, respectively. We show the two grating spectra, the prism spectrum, and the composite imaging of J-138717 in Figure 1.

Fig. 1. Spectrum of the JWST of J-138717. The rest-frame wavelengths are calculated with respect to the spectroscopic redshift, z = 1.8845 (see sect. 3.1). The combined RGB imaging of the galaxy is shown in the upper right corner. The observed spectra are shown in the upper panels. The black lines show the medium-resolution G140M/F070LP and G395M/F290LP spectra, while the overplotted orange line shows the prism spectrum (we note that the prism spectrum also covers the rest-frame interval ∼6500−10 000 Å; see, e.g., Fig. A.1). The two lower panels show the comparison of the flux of the prism and grating spectra (black line). In both panels, the magenta line marks a difference of zero, and the gray lines and shadow regions indicate the median and the average errors, respectively. The vertical red line in the left panels marks λ = 1.27 μm, above which the spectrum is potentially affected by second-order contamination. In all panels, the vertical dashed lines correspond to different absorption and emission lines.

In Sect. 5 of D’Eugenio et al. (2025), the authors show that the JADES spectra suffer from flux calibration issues. In the example shown in their Fig. 5, they report a 10–15% mismatch in the flux for G140M, with respect to the prism spectrum, at wavelengths above ∼1.3 μm. This is a known issue for JWST NIRspec, however, because the F070LP filter of the G140M grating is affected by second-order contamination at wavelengths longer than λ = 1.27 μm. For the G395M grating, they show evidence of a systematic difference in the flux all over the spectral range. To address these issues, we also compare the prism and grating spectra of J-138717 in the lower panels of Fig. 1. Differently from the galaxy shown in Fig. 5 of D’Eugenio et al. (2025), the G140M and prism spectra of J-138717 are consistent with one another along almost all the spectrum, and in particular at λ>1.27 μm, indicating that, for J-138717, the second order contamination is not significant. On the other hand, at λ<0.9 μm the prism and grating spectra differ significantly, but this difference is likely due to the very low S/N of the G140M spectrum. The flux of the G395M spectrum is systematically offset by ∼12%, on average, with respect to the prism, which is comparable to the offset reported by D’Eugenio et al. (2025). This offset is roughly constant along most of the spectrum, while there is some dependence on the wavelength below ∼3.4 μm and above ∼4.3 μm (∼1.2 μm and ∼1.8 μm rest-frame, respectively).

To further address the flux calibration issues, we considered the independent data reduction of Brammer (2022) (updated the 2024 version; hereafter, B24). In Appendix A we discuss in detail the differences between the two spectra (see Fig. A.1 for a direct comparison of the spectra of two reductions). In brief, the G140M spectra of the DR3 and B24 are consistent with each other over the whole wavelength range, although the G140M spectrum of B24 is cut at 1.27 μm (∼4400 Å rest-frame), so we cannot compare it to the DR3 at longer wavelengths. The G395M spectra of the two reductions are consistent up to ∼4 μm (∼1.4 μm rest-frame), while they differ significantly at longer wavelengths. Finally, the prism spectrum shows a clear difference in the slope of the spectrum affecting the full wavelength range. Considering the nature of the problem, it is difficult to determine which is the most reliable data reduction. We refer to the spectra reduced by the DR3, but we discuss the consistency of the results with those obtained using the spectra reduced by B24 when needed.

2.2. Templates

Throughout this paper, we make use of several libraries. For the stellar templates, we considered the MILES (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011) and IRTF + IRTF-extended (Rayner et al. 2009; Villaume et al. 2017) libraries. The MILES library consists of 985 stars observed in the wavelength range ∼3500−7500 Å, with an average resolution R ∼ 2000. The IRTF and IRTF-extended libraries consist of a total of 385 stars with different wavelength ranges. For this paper, we use a version of the library where all stars have been rebinned to a common wavelength range ∼0.8−5.2 μm having an average resolution of R ∼ 2000 from 0.8 to 2.5 μm and R ∼ 2500 from 2.5 to 5.2 μm (see Röck et al. 2017 and La Barbera et al. 2024 for details); in the following, we refer to this version of the library simply as IRTF. As for the libraries of simple stellar population (SSP) models, we considered the following.

• (i)

  The EMILES models (Vazdekis et al. 2016). These models cover the wavelength range 1680−50 000 Å and are based on fully empirical stellar libraries, namely the NGSL library (Gregg et al. 2006) in the ultraviolet, the MILES (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011), the Indo-US (Valdes et al. 2004), and CaT (Cenarro et al. 2001) libraries in the optical/NIR, and the IRTF library (Villaume et al. 2017) in the NIR/IR, with a spectral resolution of 2.51 Å full-width at half maximum (FWHM) in the optical and 60 km s⁻¹ in the NIR (see Vazdekis et al. 2016 for details). For the present work, we make use of an updated version of the models that differs from that described in Vazdekis et al. (2016) only in the spectral range ∼0.82−2.5 μm, where 205 additional stellar spectra from the extended IRTF library have been included, on top of the previous ones, while the computation of the models follows the same prescriptions and ingredients as for the original EMILES models (Vazdekis et al., in prep.). They have been recently used in La Barbera et al. (2024).

• (ii)

  The X-Shooter Library of SSP models (Verro et al. 2022a, hereafter XSL) is entirely based on the X-Shooter stellar library (Verro et al. 2022b), covering the wavelength range 0.3−2.5 μm, with a spectral resolution R ∼ 10 000.

• (iii)

  The Conroy et al. (2018) models (hereafter CvD18) models are the updated version of the Conroy & van Dokkum (2012) models. They cover the wavelength range 0.35−2.5 μm and are based on the MILES and extended IRTF stellar library (Villaume et al. 2017) in the optical and NIR, respectively, having a fixed instrumental resolution of σ_inst = 100 km s⁻¹. One key feature of these models is that they include models with non-solar values in the abundance ratios for different elements.

• (iv)

  The Maraston (2005) models (hereafter M05) are semi-empirical models having a prescription for evolved stars determined by the fuel consumption theorem of Renzini & Voli (1981). They are based on the BaSel library (Westera et al. 2002) and the Lançon & Wood (2000) empirical library of TP-AGB stars, covering a wide range of wavelengths, 91−1.6×10⁶ Å, and having a relatively low resolution, R ∼ 100. In the rest-frame wavelength range covered by J-138717, 0.35−1.4 μm, this corresponds to a FWHM varying from 10 Å to 50 Å.

• (v)

  The Maraston (2013, described in Noël et al. 2013; hereafter, M13) models are an updated version of the M05 models, adopting the same libraries and having the same instrumental characteristics. The only two differences are: (i) they assume a different fuel consumption, resulting in a heavy contribution of TP-AGB stars for the M05 models and a mild contribution for the M13 models; (ii) the isochrones used to compute the M13 models have been re-calibrated in the age range 0.2−1 Gyr by fitting the Color-Magnitude diagram of the Magellanic Cloud star clusters.

We note that all these models are computed for different isochrones and IMF types, and are available for large ranges of ages and metallicities. Here, we considered the BaSTI isochrones (Pietrinferni et al. 2004) for EMILES, the PARSEC isochrones (Bressan et al. 2012; Marigo et al. 2013; Chen et al. 2015) for XSL, the MESA Isochrones and stellar tracks (MIST isochrones) (Choi et al. 2016; Dotter 2016) for CvD18, and the Padova isochrones (Girardi et al. 2000) for M05 and M13. As for the IMF, we adopted the Kroupa (2001) IMF for all the libraries, except for the EMILES models for which we used a bimodal IMF with logarithmic slope Γ_b = 1.3, which resembles the Kroupa IMF. As for the maximum age of the models we considered the age of the models closest to the age of the Universe at the redshift of the J-138717 (∼3.4 Gyr). We also constrain the metallicity between [M/H] = −1.80 and 0.35 dex. We verified that including models with older ages and lower and higher metallicities has a negligible impact on the results, while increasing the computational time. In Table 1 we summarize the relevant stellar population properties of the libraries we adopted.

3. Physical and stellar population parameter estimate

3.1. Integrated stellar velocity dispersion

We estimated the stellar velocity and velocity dispersion (σ_*) by fitting the grating spectra G140M and G395M, separately, using pPXF (Cappellari & Emsellem 2004; Cappellari 2017, 2023, updated to version 9). To account for residuals in the background subtraction during the data reduction, we included additive polynomials in the fit. To choose the degree of the polynomials, we initially performed a fit of the spectrum using additive polynomials of degree 10. From this fit, we computed the standard deviation of the residuals (σ_std) and masked all the spectral pixels deviating more than 3σ_std from the best fit. Then, we iteratively fit the spectrum with additive polynomials at increasingly higher degrees, starting from 2. For each iteration (degree), we computed the best-fit velocity and velocity dispersion and estimated the relative errors (see below). As a result, we found that the best-fitting velocity and velocity dispersions are consistent within the errors for the order of the polynomials larger than 10 and 12 for the G140M and G395M spectra, respectively. For this reason, as our fiducial estimate for the velocity dispersion and errors of J-138717 we considered the best-fitting solutions found using additive polynomials of degree 10 and 12 for G140M and G395M, respectively. To assess the errors, for each iteration we performed a wild bootstrapping of the residuals by resampling the best-fit spectrum with the noise Gaussianly distributed from the residuals. For each iteration, we performed 100 realizations and considered the standard deviation as the relative error.

As templates to run pPXF, we used both the stellar libraries and SSP models. We used the MILES stars to fit the G140M spectrum, corresponding to the optical rest-frame, and the IRTF library to fit the G395M spectrum, corresponding to the NIR rest-frame spectrum, because the MILES (IRTF) stars do not cover the NIR (optical) rest-frame wavelength range. As SSP models, we considered the EMILES, XSL, and CvD18 models. We did not consider the M05 and M13 models because they have lower resolutions than the gratings. We note that the IRTF, XSL, and CvD18 libraries are affected by tellurics at ∼1.3−1.4 μm (Verro et al. 2022a; Villaume et al. 2017), so we masked this rest-frame wavelength region when performing the fit for the G395M spectrum with these templates.

All estimates are consistent within the errors with average velocity dispersion σ_* = 198 km s⁻¹ and standard deviation of 10 km s⁻¹, which we adopt as our reference value and associated uncertainty. The systemic velocities are also consistent and correspond to a spectroscopic redshift z_spec≃1.8845 (lower than z_phot = 2.16) that we manually verified by checking the wavelengths of the most prominent emission and absorption lines.

3.2. Stellar population properties

To estimate the stellar age and metallicity of J-138717, we performed full-spectral fitting using pPXF on the G140M, G395M, and prism spectra separately. For the G140M spectrum, we fit the rest-frame wavelength range λ∼3600−6500 Å, masking out the region around the NaD line (λ∼5900) that is not well reproduced by stellar population models, as its strength is determined by neutral gas absorption; we further discuss NaD in Sect. 4.1. We excluded from the fit the NUV region because of the low S/N (≲5/Å, on average). For the G395M spectrum, we fit the rest-frame wavelength range λ∼1−1.4 μm, since at longer wavelengths the flux calibration is uncertain, as discussed in section 2.1. For consistency with the two gratings, we restricted the fit of the prism to the rest-frame spectral region ∼0.36−1.4 μm. For the prism, we had to mask the region 5000−6000 Å (rest-frame), because it cannot be fit simultaneously with the other wavelength ranges by any model library (or even other fitting codes). We note that we can fit the 5000−6000 Å region when considering the G140M spectrum. Since the G140M spectrum is consistent with that of the prism in that wavelength region, while there is a systematic offset with respect to the G395M spectrum (see Fig. 1), we attribute the mismatch to flux calibration issues of the observed spectrum in the NIR. We further discuss this problem in Section 5.3.

For the G140M spectrum, we fit a two-parameters (A_V, δ) attenuation law for the dust, as modeled in Cappellari (2023) (Eq. (23)) and multiplicative polynomials of degree 8. For the G395M spectrum, we only fit multiplicative polynomials of degree 14; here, we do not include a dust attenuation law in the fit as the NIR rest-frame spectrum is not significantly affected by the dust, and we verified that including it does not change significantly our results. Analogously to the additive polynomials described in Section 3.1, we chose the orders of the multiplicative polynomials by iteratively increasing the degree during the fit until it converged to a constant, consistent solution. For the prism, we fit the dust attenuation but we do not include multiplicative polynomials due to the relatively low number of spectral pixels. We verified that including them does not affect significantly the results. For the G140M and prism spectra, we simultaneously fit the stars and the gas emission lines of the [OII] doublet (λ = 3726,3728 Å), [NeIII] (λ = 3869 Å) [OIII] (λ = 5007 Å), and the Balmer lines. For the prism spectrum, we also fit the [NII] line (λ = 6584 Å) that is blended with Hα, as we discuss in the next sections. When fitting the G140M and G395M spectra we constrained the stellar velocity dispersion to be within 3 standard deviations from our reference value (i.e., within σ_* = 198±30 km s⁻¹; Sect. 3.1). On the other hand, we let pPXF fit the kinematics of the prism, since the instrumental resolution is much lower than σ_*.

We then estimated the best-fitting models as follows. For each spectrum, we performed a first fit and computed the standard deviation of the residuals (σ_std). Then, we performed a second fit and masked all the spectral pixels deviating >3σ_std. We consider the results of this second fit as our best fit. The best fits of the three spectra are shown in Fig. 2.

Fig. 2. Best-fit of the prism (upper panel), G140M (middle panel), and G395M (lower panel) of the EMILES models. The three spectra were fit separately. In each panel, the black line is the rest-frame observed spectrum. The red line is the best-fit to the stellar continuum. The orange line is the best-fit to the emission lines. The gray shaded regions are excluded from the fit. The green diamonds are the residuals and the blue horizontal line is their median value.

We calculate the light-weighted ages and metallicities from the best-fits as:

$${\log }_{10}\mathrm{Age}=\frac{{\sum }_i{\mathit{w}}_i{{\log }_{10}\mathrm{Age}}_i}{{\sum }_i{\mathit{w}}_i}$$(1)

$$[\mathrm{M}/\mathrm{H}]=\frac{{\sum }_i{\mathit{w}}_i{[\mathrm{M}/\mathrm{H}]}_i}{{\sum }_i{\mathit{w}}_i}$$(2)

where w_i is the weight of the i-th input template and the sums are performed over all the templates. To estimate the errors of the age, metallicity, and dust attenuation parameters, we performed a wild bootstrapping of the residuals. Specifically, we performed 100 realizations of each spectrum by shuffling the residuals over windows ∼300 Å wide, to account for wavelength correlated errors, and then re-distributing (i.e., adding) the shuffled residuals to the best-fitting spectrum Gaussianly. For each realization, we estimated the age and metallicity and took the standard deviations of all the realizations as the errors. Using the EMILES models, from the fits of the optical (G140M), NIR (G395M), and prism spectrum we estimate an age of 0.91±0.21, 0.95±0.25, and 0.89±0.11 Gyr, respectively, and a metallicity of −0.28±0.17, −0.38±0.15, and −0.35±0.08 dex. The attenuation parameters estimated from the optical and prism spectrum are A_V = 0.65±0.18,0.63±0.05 mag and δ=−0.82±0.29,−0.90±0.12, respectively. In Fig. 3 we show the best fits to the emission lines. We estimate F(Hα) = 4.7±0.5×10⁻¹⁹ erg cm² s⁻¹, F([NII]) = 1.5±0.1×10⁻¹⁸ erg cm² s⁻¹, F([OII]) = 5.2±0.5×10⁻¹⁹ erg cm² s⁻¹, F([NeIII]) = 1.7±0.5×10⁻¹⁹ erg cm² s⁻¹, and F([OIII]) = 4.0±0.2×10⁻¹⁹erg cm² s⁻¹; except Hα, the Balmer lines are consistent with having zero flux in emission.

Fig. 3. Zoom-in of the best fits for the emission lines. The lines and the symbols are the same as in Fig. 2. The magenta line in the lower right panel show the best-fit templates for the Hα and [NII] lines separately.

In Appendix A.2, we show the results of the fits with the other SSP libraries. Overall, the stellar population parameters are consistent when using different libraries or fitting different spectral ranges. In Appendix A.2, we also show how different the stellar population properties are when fitting the B24 prism spectrum. In general, the latter provides older ages and higher metallicities, while the attenuation parameters are consistent. This difference is expected since the B24 prism spectrum is flatter than that of the DR3. On the other hand, we find fully consistent results when fitting the G140M and G395M spectra for all SSP libraries.

Since the stellar continuum is attenuated by dust, possibly affecting the age and metallicity estimates, we re-derive these parameters using the Bayesian full-index fitting (FIF) method (Martín-Navarro et al. 2019, 2021). The FIF method is a hybrid approach between the full-spectral fitting and the line-strength analysis. The basic idea is to use as much spectral information as possible, like in the full-spectral fitting, but to limit the fit to a number of carefully selected spectral indices, based on their sensitivity to the stellar population parameters to be determined. For our science case, the key advantage of the FIF is that it relies on the spectral indices, thus minimizing the source of uncertainties affecting the stellar continuum such as the dust absorption or flux calibration issues.

The FIF method performs a pixel-by-pixel fit of all the spectral pixels within the bandpasses of the spectral features we fit. We used the emcee Bayesian Monte Carlo Markov chain sampler (Foreman-Mackey et al. 2013) to maximize the following likelihood function:

$$ln(\text{O}|\text{S})=-\frac{1}{2}{\displaystyle \sum _n}\left[\frac{(O_n-M_n{)}^2}{{\sigma }_n^2}-ln\left(\frac{1}{{\sigma }_n^2}\right)\right]\hspace{0.5em},$$(3)

where the sum extends over all the spectral pixels within the bandpasses we fit, with O_n and M_n being the observed and the model fluxes, respectively, of the n-th spectral pixel, while σ_n is the uncertainty on O_n. As for the spectral indices, we considered Hβ, Hγ, and Hδ, which are primarily sensitive to the galaxy's age, and the Fe4383, Fe4531, Fe5270, and Fe5335, which are sensitive to the metallicity. We also include the Mgb, which is sensitive to both metallicity and [α/Fe] abundance. To estimate the stellar population properties we adopted the EMILES models, and we verified that the XSL library provides consistent results. Given that the FIF is based on lower statistics with respect to the full-spectral fitting and does not rely on the stellar continuum, we did not fit the M05 and M13 models because the resolution of the observed spectrum must be reduced for them, and we did not fit the CvD18 models either because they are not computed at ages lower than 1 Gyr, which would thus affect the metallicity.

The FIF best-fitting age is $0.{81}_{-0.07}^{+0.05}$ Gyr, consistent with the results of the full-spectral fitting, with a metallicity of $-0.{19}_{-0.07}^{+0.11}$ dex, which is consistent with the estimate of the full-spectral fitting of the optical spectrum, albeit ∼0.1 dex higher (∼0.2 with respect to the fit of the NIR and prism). In Fig. 4 we show how the observed indices compare to the best-fit model (i.e., with the estimated median values of age and metallicity). The Balmer lines, which are more sensitive to the age, are well fit. The iron lines Fe4531 and Fe5270 are also well reproduced, while Fe4383, Fe5335, and Mgb are underestimated. We verified that when including other spectral indices sensitive to the metallicity (CN1, CN2, Ca4227, CN3883) the metallicity is slightly lower ($-0.{23}_{-0.11}^{+0.07}$), yet still systematically higher than the one calculated from the full-spectral fitting.

Fig. 4. Comparison of the observed indices (black dots show the spectral pixels) with the predictions of the FIF best-fit models (magenta line), with an age of 0.81 Gyr and a metallicity of −0.19 dex. The blue and red regions correspond to the wavelengths where the pseudo-continua are measured and the white regions are the wavelengths where the index is measured. In the bottom panels the residuals are plotted as green diamonds. The gray lines and shaded regions are the median of the residuals and the average observed errors on the indices, respectively, while the black dashed line marks the 0.

3.3. Mass, SFR, and SFH

To estimate the stellar mass of J-138717, we fit the photometry provided by the JADES photometric catalog using the convolved Kron magnitudes. To fit the data, we used the C++ implementation of the FAST code (Kriek et al. 2009)¹ fed with the Bruzual & Charlot (2003) models using a Kroupa IMF and [M/H] = –0.4 dex. Further, we included a Kriek & Conroy (2013) dust attenuation law in the fit. We estimate a stellar mass M_* = 3.5±0.2×10¹⁰ M_⊙². The FAST best fit is shown in Fig. 5. FAST also estimates an age of 0.9 Gyr and attenuation parameter A_V = 0.8 mag, similar to the values obtained from the full-spectral fitting.

Fig. 5. FAST best fit to the photometric data (red circles). Specifically, we used the convolved Kron magnitudes in the JWST filters F182M, F210M, and F444W, and the HST filters F606W, F775W, F814W, F850LP, as listed in the JADES photometric catalog.

We estimate the SFR from the Hα and [OII] emission lines as follows. The fluxes reported in the previous section correspond to luminosities L_Hα = 1.2±0.2×10⁴⁰ erg s⁻¹ and L_[OII] = 1.3±0.1×10⁴⁰ erg s⁻¹, given that $L=F\times 4\pi d_L^2$, where d_L = 1.44×10⁵ Mpc is the luminosity distance at the redshift of the galaxy. We correct the observed luminosities for the dust attenuation as L_corr=L_obs×e^τ_λ, with τ_λ=A_λ/1.086 and A_λ=k_λA_V/R_V, where A_V = 0.65 mag, as estimated by the fit of the optical spectrum, k_λ is the reddening curve for the gas, for which we assume a Calzetti et al. (2000) law with R_V = 4.05. Thus, τ_Hα≃1.045 and τ_[OII]≃0.593. We then calculate the SFR from Hα using Equation (2) of Kennicutt (1998), SFR(Hα) = 7.9×10⁻⁴² L_Hα,corr = 0.16±0.03 M_⊙ yr⁻¹, and from [OII] using Equation (4) of Kewley et al. (2004), SFR([OII]) = 6.58±1.65×10⁻⁴² L_[OII],corr = 0.24±0.09 M_⊙ yr⁻¹. Given the stellar mass estimated above, we derive the specific star formation rate (sSFR) as SFR(Hα)/M_* = 5±1×10⁻¹² yr⁻¹ and SFR([OII])/M_* = 7±2×10⁻¹² yr⁻¹, both about 1 dex below the limit defined by Gallazzi et al. (2014), that is, an sSFR < 1/5t_U(z)≃6×10⁻¹¹, where t_U≃3.4 Gyr is the age of the Universe at the galaxy's redshift. This confirms the quiescent nature of J-138717. We note that emission lines usually suffer more dust attenuation than the stars. For local starburst galaxies, Calzetti et al. (2000) estimated that the color excess for stars is a factor 0.44 smaller than the gas (Equation (3)). This would imply an A_V a factor ∼2.3 stronger for the gas than stars, such that we could underestimate the SFR by a factor ∼2−4 (depending on the line considered). Even in this case, however, we would measure an sSFR lower than 1/5t_U(z). We also note that the resolution of the prism did not allow us to resolve the blending between the Hα and [NII] emission lines, so it is possible that the SFR derived from Hα is underestimated. For this reason, we re-fit the prism spectrum without [NII], that is, assuming that all the emission lines are due to Hα. We measure F(Hα) = 1.1±0.5×10⁻¹⁸ erg cm²s⁻¹, resulting in a SFR(Hα) = 0.37±0.15 and then sSFR = 11±5×10⁻¹² yr⁻¹, still consistent with the previous estimates. On the other hand, we note that, if J-138717 is a LINER (see Sect. 5.1), then we are overestimating SFR(Hα), as the Hα emission would be caused by an AGN or post-AGB stars, and not by the young stars.

We reconstruct the SFH of J-138717 in a non-parametric way, by performing the fit of the G140M spectrum with pPXF using the same setup described in the previous section but allowing for regularization (Cappellari 2017, 2023). The SFH is shown in Fig. 6. From the cumulative sum of the mass-weights, shown in the right panel, we estimate the lookback times at which the galaxy formed 20%,50%,and 90% of its mass, $t_{20}^{\mathrm{lbt}}=0.8\pm 0.1$ Gyr, $t_{50}^{\mathrm{lbt}}=0.6\pm 0.1$ Gyr, and $t_{90}^{\mathrm{lbt}}=0.4\pm 0.1$ Gyr, corresponding to redshifts z₂₀≃2.6, z₅₀≃2.3, and z₉₀≃2.1, respectively. Given Δt=t₉₀−t₂₀≃0.4 Gyr, the average past SFR with which J-138717 formed its stellar mass is $\langle \text{SFR}\rangle ={\text{M}}_{*}/\mathrm{\Delta }t\simeq 88$ M_⊙ yr⁻¹. We verified that fitting the NIR and prism spectra provides qualitatively the same SFH, but the lookback times are systematically older (for the NIR we get $t_{20}^{\mathrm{lbt}}=1.2\pm 0.2$ Gyr, $t_{50}^{\mathrm{lbt}}=0.8\pm 0.1$ Gyr, and $t_{90}^{\mathrm{lbt}}=0.5\pm 0.1$ Gyr) and the two Δt are almost a factor 2 longer, resulting in a halved $\langle \text{SFR}\rangle$. This is in agreement with the older (∼0.1 Gyr) age found when performing the full-spectral fitting of the NIR compared to the optical (see Sect. 3.2). This difference could arise because the Balmer lines are more sensitive to the most recent stellar burst (i.e., youngest population), while the NIR range is more sensitive to the older stellar component, thus hinting at the presence of a non-negligible fraction of an older (≳1 Gyr) population. We caution that models in the NIR are uncertain at such young ages, as we show in the next section and further discuss in Section 5.

Fig. 6. Star formation history inferred from the regularized fit of the optical spectrum. Left panel: Map of the weights assigned by pPXF to each input template of given age and metallicity. Right panel: Cumulative sum of the weights. The green, blue, and red line correspond to $t_{20}^{\mathrm{lbt}},t_{50}^{\mathrm{lbt}},\text{and}{t}_{90}^{\mathrm{lbt}}$, respectively.

4. Analysis of the spectral indices

In this section, we compare the observations of the spectral indices with the predictions of the different models libraries, using the definitions of Trager et al. (1998) and Eftekhari et al. (2021) for the optical and NIR indices, respectively. Considering the data reduction issues discussed in Sect. 2.1 (see also Appendix A.1), we study only those indices whose measurements are consistent when measured from the spectra of the DR3 and B24 data reductions. We also consider the optical spectral indices at rest-frame wavelengths λ>4400 Å, although the G140M spectrum of B24 does not reach longer wavelengths, so we cannot evaluate the consistency of the measurements. Although consistent, some spectral indices are significantly different when measured from DR3 or B24. For this reason, we verified that the qualitative results of the following analysis do not change when considering the indices measured from B24. We found that the only index providing qualitatively different results is AlI1.67 (see Section 4.3).

To compare observations with model predictions we have to match the different resolutions. Among the libraries we considered, the M05 and M13 models have the lowest resolution, which is also lower than that of the gratings, varying from 10 Å−100 Å FWHM in the rest-frame wavelength range of J-138717. For this reason, we convolved the observed grating spectra and the EMILES, XSL, and CvD18 models by a Gaussian corresponding to the instrumental resolution of the M05 (and M13) models at the wavelengths corresponding to each spectral index.

As the reference age we consider the one estimated from the full-spectral fitting of the optical spectrum using the EMILES models (sect. 3.2), although for a proper comparison we should use the predictions from the different models. In practice, the qualitative results would not change for most of the libraries. The only exception is for the M13 models, which predict a systematically younger age (∼0.6 Gyr; see Fig. B.1) compared to all the other libraries. In many cases, M13 models predict values of the spectral indices closer to observations if an age of 0.6 Gyr is assumed. We highlight this consistency in the following sections, when needed.

4.1. Probing the C, Na, and α enhancement

In Fig. 7 we compare the observed C₂ spectral index at 4668 Å with predictions of the different models. None of the models can reproduce the observations. Since this C₂ index is particularly sensitive to the C abundance, the observed mismatch could indicate that J-138717 is C-enhanced. To probe this, using the CvD18 model at 1 Gyr and solar metallicity, we calculate the difference in C₂ for a C-enhancement of [C/Fe] = +0.15 dex, shown in Fig. 7 with a brown arrow. Assuming C-enhancement significantly reduces the mismatch, but observations are still inconsistent with model predictions. This could imply a C-enhancement much higher than [C/Fe] = +0.15. This is in contrast with the recent results of Beverage et al. (2025), who found a carbon deficiency in quiescent galaxies at comparable redshifts (although they also report a galaxy with a C-enhancement of +0.3 dex). We note, however, that C₂4668 is sensitive to the age (see the steep rise for the EMILES models in Fig. 7), and the possible presence of a population of stars older than ∼1 Gyr, as the SFH indicates, could account for part of the mismatch. Finally, we verified that changing the IMF or [α/Fe] has a negligible effect on this index at such a young age. As an example, in Fig. 7 we show the effect of increasing the [α/Fe] from solar-scaled to +0.4 dex³.

Fig. 7. C24668 spectral index as a function of age. The yellow star is the measured value from the G140M spectrum of J-138717. The green, red, gray, blue, and purple lines are the predictions of the M05, M13, XSL, CvD18, and EMILES models, respectively. Dark to light colors indicate higher to lower metallicities. Both observations and models have been smoothed to the instrumental resolution of the BaSeL stellar library. The brown (olive) line in the central (right) panel indicates by how much the tracks would shift towards higher index values due to a C-enhancement(α-enhancement) of 0.15 dex (0.4 dex), as predicted by the CvD18 (MILES) models.

Similarly to C₂4668, models largely underpredict the NaD spectral index. We verified that the strength of the NaD index cannot be recovered by any model even assuming a strong Na-enhancement, a steep IMF or α-enhancement, and can only be explained by the presence of significant dust absorption, as we infer from the fits (Sect. 3.2). In particular, we measure an equivalent width of 6.2±0.5 Å for NaD, consistent with the measurements of Belli et al. (2024) for a PSB observed at z = 2.45. Having EW(NaD) = 6.5±0.5 Å, the authors estimate a hydrogen column density N_H⁴ of the order of ∼10²¹ (see their Equation (2)). Following Rachford et al. (2009), we estimate N_H/E(B−V) = 5.8×10²¹ H cm⁻²mag⁻¹; assuming R_V=A_V/E(B−V) = 3.1 and substituting one gets A_V∼0.7 mag. Although this estimate is quite uncertain and relies on a number of assumptions, the result is consistent with the estimate of the full-spectral fitting and confirms the presence of a significant dust component.

As shown in Fig. 4, J-138717 has a significant Mgb absorption that is not well matched by the EMILES models and could be indicative of a super-solar [α/Fe] abundance. To probe α-enhancement, we performed the FIF, as described in Section 3.2, fitting [α/Fe] with the MILES models having [α/Fe] = +0.0 and +0.4. From this fit, we estimate [α/Fe] = 0.15±0.10 dex, hinting at α-enhancement. Using a more traditional approach (e.g., Thomas et al. 2003; Bevacqua et al. 2023), we constructed grids of Mgb vs. Fe3, with Fe3 = $\frac{1}{3}\left(\mathrm{Fe}4383+\mathrm{Fe}5270+\mathrm{Fe}5335\right)$ by interpolating the MILES models at different [α/Fe] and metallicity, having a fixed age of 0.8 Gyr. We then estimate [α/Fe] by minimizing χ², defined as

$${\chi }^2={\displaystyle \sum _j}{\left(\frac{{\mathrm{I}}_{\mathrm{obs},\mathrm{j}}-{\mathrm{I}}_{\mathrm{mod},\mathrm{j}}}{{\sigma }_{\mathrm{I},\mathrm{j}}}\right)}^2\hspace{0.5em}\hspace{0.5em}\hspace{0.5em},$$(4)

where j runs over the indices Mgb and Fe3, I_obs,j are the observed indices, I_mod,j are the indices of the model grid, and σ_I,j is the measured error. With this method, we estimate [α/Fe] = 0.19±0.13. Both the FIF and model grid suggest that J138717 might be α-enhanced, but the uncertainties are large and a higher S/N is required to safely constrain [α/Fe].

In summary, the analysis of C₂4668 suggests that J-138717 is strongly C-enhanced ([C/Fe]> 0.15 dex), but the actual enhancement is questioned by the possible presence of an old (>1 Gyr) stellar component affecting the index. The analysis of NaD does not allow us to estimate the [Na/Fe] abundance, since the strength of the index is dominated by the dust. Assuming that the measured EW(NaD) is determined by the dust absorption, we roughly estimate A_V∼0.7 in agreement with the fits. Finally, the Mgb and iron lines indicate that J-138717 is α-enhanced ([α/Fe] ∼ 0.15 dex), but the S/N is too low to put a stringent constrain on [α/Fe].

4.2. Probing the contribution of TP-AGB stars

The most sensitive spectral indices to TP-AGB stars are the molecular CO and CN bands (Maraston 2005; Riffel et al. 2007; Martins et al. 2013). In Fig. 8 we show these indices as observed for J-138717 and compared to the model predictions. We recall that M05 and M13 are semi-empirical models accounting for a heavy and mild contribution of TP-AGB stars, while EMILES, XSL, and CvD18 are based on fully empirical libraries and do not have a significant contribution.

Fig. 8. CN indices as a function of age. The marker and the lines are the same as in Fig. 7.

In the first two rows of Fig. 8 we show the CO1.56, CO1.58, CO1.60, CO1.64, and CO1.68 indices. In general, none of the libraries can simultaneously match all the observed CO indices. Models fully based on empirical stellar libraries (EMILES, XSL, and CvD18) generally provide more consistent measurements with observations. The inferred metallicities vary significantly from one index to the next and from one library to the next, however. For instance, the EMILES and XSL models predict a subsolar metallicity consistent with that estimated from full-spectral fitting for CO1.58 and CO1.60, while CO1.56 and CO1.68 suggest much lower metallicity, and finally, CO1.66 requires a supersolar metallicity to match the observations. The M13 models would be consistent with the observations only for the CO1.56 and CO1.58 index, assuming an age of ∼0.6 Gyr, as predicted by the fits, while the other predictions do not reproduce the observations. Finally, the M05 models do not match any of the observed CO indices. We verified that the M05 and M13 models could match all the observed CO indices only when assuming a very high metallicity (>0.35 dex), but this is in sharp contrast to the metallicity estimated from the fits, including those obtained with the M05 and M13 models (Fig. B.1).

In the third and fourth rows of Fig. 8 we show the CN1, CN2, CN3883, CN1.10, CN1.42, and CN1.46 indices. In this case, all libraries can match the optical CN indices (CN1, CN2, and CN3883). The EMILES models can also reproduce all the observed NIR indices of the CN with subsolar metallicities (although CN1.42 and CN1.46 require metallicities much lower than that estimated from the fits). The XSL models can reproduce the CN1.10 and CN1.42 (and only approximately CN1.46) indices, but the inferred metallicity again depends on the index considered. The CvD18 models tend to underpredict all the NIR CN indices and they would require supersolar metallicities to match the observations. The M13 models are consistent only with the CN1.10 in the NIR, but they would be consistent with all CN indices assuming an age of 0.6 Gyr. Finally, the M05 model predictions do not match any of the NIR CN indices. Similarly to the CO indices, the M05 and M13 models could match the CN indices assuming a very high metallicity.

For all CO and CN indices, we also checked the effect of C-enhancement, α-enhancement, and IMF variation and found contrasting results. For instance, a C-enhancement would decrease both the CO1.66 and CO1.68 predicted by the models, resulting in a worse match for the former and a better match for the latter when compared to observations. Similar results are found for the other indices and a similar effect is obtained when varying the IMF, while varying [α/Fe] has a negligible effect in all cases. Finally, we note that La Barbera et al. (2024) verified that the use of different isochrones does not significantly affect the values of the indices measured from the models.

Overall, none of the libraries can simultaneously match all the CO and CN observed indices. The models based on empirical stellar libraries (EMILES, XSL, and CvD18) generally predict CO and CN indices closer to observations than the theoretical ones (M05 and M13), but the inferred metallicity can vary significantly from index to index and from library to library. Instead, the theoretical models struggle to reproduce the observations, and in particular the NIR indices. The M13 models provide predictions consistent with the observations of all CN indices and two CO indices assuming an age of ∼0.6 Gyr, consistent with the estimates from full-spectral fitting, but systematically lower than the estimates of all the other libraries. Both the M05 and M13 models could match observations assuming a very high metallicity, but this is in contrast with the metallicity estimates from the fits. In conclusion, our analysis indicates that there is no heavy contribution of TP-AGB stars to the spectrum of J-138717.

4.3. Other NIR spectral indices

We describe the results based on other NIR spectral indices. As a result of editorial limitations, however, we have to show these results in Appendix C. In Fig. C.1 we show the comparison of the measured and predicted Paschen lines, which are age indicators of a galaxy, as they do not depend significantly on the metallicity. The Paβ line is consistent with all libraries. The Paγ is in agreement with the EMILES, XSL, and CvD18 models. Instead, the M05 and M13 models underpredict this index. Similarly to the CN indices, considering an age of 0.6 Gyr would provide consistent predictions for the M13 (and M05) models. Finally, none of the libraries can reproduce the Paδ line (barely, the EMILES models with solar metallicity). We note that the Paδ predictions for the M05 and M13 models are not shown because, due to the large wavelength sampling of the BaSeL library from which they are constructed, there are only two spectral pixels sampling the Paδ bandpass, resulting in an unrealistic flat shape of the index.

In Fig. C.2 we show other NIR spectral indices. Overall, none of the libraries can match all these indices. In general, the M13 models reproduce most of these indices, but the inferred metallicity changes significantly from index to index. The NaI1.14 and FeI1.16 indices are consistent with all libraries, except the M05 models. The AlI1.31 index is consistent with all libraries, but the EMILES (M05 and M13) models require very low (high) metallicities. The MgI1.50 index is consistent with the M13, XSL, and CvD18 libraries that have different metallicities. Finally, the AlI1.67 index is inconsistent with all libraries. When measured from the B24 spectrum, although it is consistent with the measurements from the DR3 spectrum (given the large uncertainty), it is much lower and consistent with all libraries, however. We conclude that the modeling of the NIR spectral indices at young ages is uncertain for all of the currrent models we considered.

5. Discussion

5.1. SFH and quenching of J-138717

The SFH shown in Fig. 6 indicates that J-138717 formed most of its mass (≥50%) recently ($t_{50}^{\mathrm{lbt}}=0.6$ Gyr) and then quenched recently and rapidly ($t_{90}^{\mathrm{lbt}}=0.4$ Gyr). The mechanisms of rapid quenching are generally attributed to gas-rich mergers and strong AGN feedback. Park et al. (2024) recently found that the majority of rapidly formed and recently quenched galaxies, such as J-138717, often exhibit strong AGN features (Davies et al. 2024; see also Taylor et al. 2024). The spectrum of J-138717 indeed shows high-ionization lines typical of AGNs, namely [NeIII], [OIII], and [NII]. The relatively low flux ratios log₁₀([OIII]/[OII]) ≃ −0.1 dex and log₁₀([NeIII]/[OII]) ≃ −0.5 dex suggest, however, that these emissions are sustained by some low level star formation. If this is the case, such residual star formation would be consistent with a scenario where stars formed very efficiently, quickly consuming most of the cold gas and then halt (Vazdekis et al. 1996, 1997; Weidner et al. 2013; Jeřábková et al. 2018; Salvador-Rusiñol et al. 2020), which could also explain the compact size of J-138717 (≃1.1 kpc).

On the other hand, the relatively high flux ratio log₁₀([NII]/Hα)≃0.5 dex suggests that J-138717 is a LINER galaxy (Kewley et al. 2006)⁵. This would be in agreement with other studies showing that many PSBs have LINER emissions (e.g., French et al. 2015), as also recently found in some other massive quiescent galaxies observed at comparable redshifts by Bugiani et al. (2024). The LINER emission we detect could be due to weak AGN activity, possibly preventing new star formation over the past 400 Myr, although post-AGB stars and shocks can also cause the LINER emissions (Belfiore et al. 2016; Yan 2018).

Quenching via strong AGN feedback in massive galaxies at high redshifts is often associated with massive multiphase outflows (Belli et al. 2024; Davies et al. 2024; Carniani et al. 2024a; Park et al. 2024; Taylor et al. 2024). These outflows are identified in the galaxy spectrum by a shift of the emission and absorption lines tracking the gas with respect to their systemic velocities or by asymmetric line profiles. For J-138717, the emission lines do not exhibit any asymmetry or shift in wavelength, as shown in Fig. 3. Then, in Fig. 9 we show the profile of the NaD absorption line, tracking the neutral gas. The observed lines are centered at the expected wavelengths, their shapes are Gaussian and symmetric, and their line strengths are comparable. We verified that the CaII (H&K) lines (shown in Fig. A.2) and MgII doublet at λ 2796,2804 Å, two other tracer of the neutral gas associated with outflows (Nestor et al. 2011; Belli et al. 2024; Taylor et al. 2024), show similar, well-behaved profiles.

Fig. 9. Observed profile of the NaD absorption line (black line). The dashed red lines indicate the expected wavelengths of the doublet at λ = 5890 Å and λ = 5896 Å.

Overall, our analysis suggests that there is no strong AGN feedback in J-138717. Because the AGN emission is not typically expected to be spherically symmetric, however, we might be missing this evidence due to the direction of the AGN emission and outflows being misaligned with respect to the observations. Other mechanisms associated with rapid quenching are wet mergers and gas-stripping. Both these cases would require a study of the environment of J-138717. Although we do not find any sources in its proximity at comparable redshift (the closest source among galaxies with 1.8<z_phot<2.3 is ∼14″ away, corresponding to ∼110 kpc), we do not attempt any detailed environmental study of J-138717, because no spectroscopic redshifts are available for close-by sources and the uncertainties on the photometric redshifts are too large (for example we measure a difference z_phot−z_spec∼0.3 for J-138717).

5.2. The low contribution of TP-AGB stars

According to the predictions of the fuel-consumption theorem (Renzini & Voli 1981), TP-AGB stars are expected to dominate the luminosity of an SSP at young ages (0.2−2 Gyr), accounting for up to 80% of the NIR and 40% of the total luminosity (Maraston 1998). The JWST observations of J-138717 offer a unique opportunity to probe the actual contribution of TP-AGB stars, as the resolution of the NIRSpec gratings allows for reliable measurements of the NIR spectral indices, and, given the young age and short SFH of this post-starburst galaxy, the contribution of TP-AGB stars should be at its peak. We highlight that such analysis has never been done before outside the local Universe. From our analysis, models based on empirical stellar libraries (EMILES, XSL, and CvD18), having a low contribution of TP-AGB stars, can reproduce a large number of CO and CN indices, tracing TP-AGB stars. The M13 models, having a mild contribution of TP-AGB stars, can reproduce many CO indices and CN indices assuming an age of ∼0.6 Gyr, which would be consistent with the fits (Fig. B.1). A younger age, combined with some residual star formation, could indeed dilute the signal of TP-AGB stars. On the other hand, we note that all the other models predict systematically older ages. Finally, the M05 models, having a heavy contribution of TP-AGB stars, do not reproduce any of the observed CO and NIR CN indices. These results suggest that the contribution of TP-AGB stars to the spectrum of J-138717 is not significant, in agreement with other PSBs studied by Kriek et al. (2010) at high redshifts and Conroy & Gunn (2010) in the local Universe.

We note that Eftekhari et al. (2022), found similar results for local massive ETGs, but we find a lower mismatch between observations and model predictions for the CO and CN indices. This might indicate that the TP-AGB contribution varies (or the modeling of the CO and CN features gets worse) at older ages, higher metallicities, and, possibly, steeper IMF slopes that characterize local massive quiescent ETGs. For example, Martínez-García et al. (2021) argued that the luminosity of TP-AGB stars may vary with metallicity as it controls their mass-loss rate. Overall, an improvement in the modeling of the NIR spectrum is still very much needed. Riffel et al. (2015) suggested that stars in advanced evolutionary phases, like red giant branch (RGB) and early-AGB stars, may help reduce the observed mismatch. We further discuss this issue in the next section.

Recently, Lu et al. (2025) reported the presence of strong absorption features, such as CN and CO, in the prism spectra of three young quiescent galaxies found at redshifts z = 1−2, similarly to J-138717, and attributed these features to a significant contribution of TP-AGB stars. Unfortunately, the resolution of the prism does not allow us to measure the NIR spectral indices. We have shown that, although J-138717 also shows strong CO and CN absorptions detectable in the prism data, the spectral indices analysis disfavors a strong contribution of TP-AGB stars. Therefore, observations with higher resolution and S/N are needed to confirm the results found by Lu et al. (2025). It is worth noting that the galaxy spectra shown by Lu et al. (2025) exhibit stronger CO, CN, and C₂ features than J-138717 as well as other features (such as TiO, VO, and ZrO) that are not detected in the spectrum of J-138717. Although this can be explained by the presence of other poorly modeled stars, like carbon-rich, early AGB, and RGB stars (Riffel et al. 2015), they are indeed three valuable candidates to further constrain the contribution of TP-AGB stars in galaxies. If confirmed, a strong contribution of TP-AGB stars in these three galaxies, at variance with J-138717, would imply that the significance of TP-AGB stars varies among galaxies. This would explain the conflicting results found in the literature and would call for considerable improvements in the modeling of evolved stars such as the TP-AGB.

5.3. General issues and caveats

This work offers the opportunity to discuss some issues related to the current data analysis and recent results in the literature. A major problem concerns the reliability of the data reduction pipeline of JWST. Although great endeavors have been made to provide reliable spectra, such as for the reduction of the JADES spectra (see D’Eugenio et al. 2025), some important problems still persist and can affect significantly the results. For example, we showed that there is a systematic offset between the prism and G395M spectrum (Fig. 1) that increases at longer wavelengths. This has been also reported in the DR3 of JADES. By comparing with B24, we showed that different data reductions provide significantly different spectra and thus different results. For example, the slope of the prism spectra of the DR3 is steeper than that of B24 (Fig. A.1), resulting in a younger age differing up to 0.5 Gyr (depending on the models). In both cases, we could not fit the wavelength region λ∼5000−6000 Å. Analogously, Lu et al. (2025) found that they were unable to fit the spectral region between λ∼5000−10 000 Å, and ascribed this issue to the poor modeling of TP-AGB stars in currrent models. This does not seem to be the case for J-138717, since our analysis based on the spectral indices still disfavors a strong presence of TP-AGB stars. Most importantly, the G140M spectrum, covering this wavelength range, does not suffer from this problem and models can easily fit the continuum. Given the consistency of the G140M and prism spectra, and the mismatch between the G395M and prism spectra (see Fig. 1), it is likely that the flux calibration of the prism is incorrect at longer wavelengths. Such systematics can have a huge impact, especially when working with young galaxies at very high redshifts such as JADES-GS-z14-0 and JADES-GS-z14-1 (Carniani et al. 2024b), as two extreme examples. So far, different teams have adopted their customized (private) pipelines, but an effort should be made by the community to provide a reliable public data reduction pipeline.

Another issue, concerns the reliability of the stellar population models. We derived stellar population parameters with many widely used models⁶. We found a generally good agreement in the estimated parameters when performing the full-spectral fitting, although the parameters have large uncertainties notwithstanding the relatively high S/N of the spectra. Instead, the analysis of the NIR spectral indices (see Sect. 4) reveals major inconsistencies, as none of the models we adopted can reproduce all the observations, and those indices that are well reproduced by models provide inconsistent metallicities. For example, models based on empirical libraries predict better CO, CN, and Pa indices (Figs. 8 and C.1), while M13 models better predict NaI, FeI, AlI, and MgI NIR absorption lines (Fig. C.2), but the metallicity derived varies significantly from index to index and from library to library. The inconsistency of the metallicity derived from the NIR indices was already pointed out by Riffel et al. (2015), who suggested that an improved treatment of TP-AGB stars, as well as red giants and early-AGB stars, might minimize the mismatch. We also note that the resolution of the M05 and M13 models (or other models in the literature using low-resolution NIR libraries) is potentially a source of major limitations in our data analysis (e.g., Dahmer-Hahn et al. 2018). For example, we have discussed how some indices, like Paδ, cannot even be measured because the FWHM is larger than the absorption line. An effort to construct models with higher resolution (in particular, those including a significant contribution of TP-AGB stars) is invaluable for future studies based on NIR spectral indices. Overall, current stellar population models should be used with caution when deriving stellar population properties (e.g., the metallicity) from the NIR, in particular at young ages.

Many recent works in the literature are discussing results based on the SFHs of high-redshift galaxies, having profound implications on our understanding of galaxy formation. The recovery of the SFH is very challenging, however, because of the intrinsic uncertainties that are difficult to overcome (e.g., Zibetti et al. 2024) and may depend on the methods adopted. For example, Glazebrook et al. (2024) found that ZF-UDS-7329, a massive (M_* = 10^11.4 M_⊙) quiescent galaxy observed at z≃3.2 that formed most of its mass already at z = 11, as also confirmed by Carnall et al. (2024). However, Turner et al. (2025) showed that the actual formation redshift can be lowered to z = 7 (i.e., doubling the age of the Universe), and the age can be by up to 0.5 Gyr younger, just by changing priors. We showed that the spectral range we fit can also provide differences of a few hundred million years in the SFH. For J-138717, this might be due to the presence of a non-negligible fraction (≥20%) of old (>1.2 Gyr) stars detected in the NIR spectrum, but, considering the analysis of the NIR spectral features, we cannot exclude that this is due to the difficulties current models have to reproduce the NIR spectrum. In the former case, this effect should be reduced at higher redshifts, when the stellar populations are homogeneously younger by construction. We note that Nanayakkara et al. (2025), except for some galaxies, found statistically consistent SFHs when fitting different stellar population models or spectral ranges. We point out, however, that the timescales they inferred from the different methods can have differences of about a few hundred million years, in agreement with our results. Considering all these uncertainties, more compelling evidence is required to shake the foundations of the ΛCDM.

6. Summary and conclusions

We studied the physical and stellar population properties of J-138717, a post-starburst galaxy observed at z_spec = 1.8845 with JWST, for which we estimated a stellar mass M_* = 3.5±0.2×10¹⁰ M_⊙ and velocity dispersion σ_* = 198±10 km s⁻¹. We summarize our results as follows.

• From the full-spectral fitting of the optical, NIR, and prism spectra, we estimated an age of 0.91±0.21, 0.95±0.25, and 0.89±0.11 Gyr, respectively, and a metallicity of −0.28±0.17, −0.38±0.15, and −0.35±0.08 dex with the EMILES models. The FIF method that fit the optical spectral indices and was fed with the EMILES models provided an age of $0.{81}_{-0.07}^{+0.05}$ Gyr and a metallicity of $-0.{19}_{-0.07}^{+0.11}$, which is consistent with the full-spectral fitting. We verified that other models generally provide consistent results.

• From the optical and prism spectra, we also estimated dust attenuation parameters A_V = 0.65±0.18 and 0.63±0.05 mag, and δ=−0.82±0.29 and −0.90±0.12 respectively. From EW(NaD), we also estimated A_V∼0.7, which is consistent with the fits. This agrees with other results in the literature that reported a significant dust component in high-redshift quiescent galaxies.

• From the Hα and [OII] emission lines, we estimated SFR = 0.16±0.03 M_⊙ yr⁻¹ and = 0.24±0.09 M_⊙ yr⁻¹, respectively, which implies an sSFR = 5±1×10⁻¹² yr⁻¹ and sSFR = 7±1×10⁻¹² yr⁻¹. This is lower by about one order of magnitude than a threshold of 1/5t_U(z) that defines quiescent galaxies.

• We reconstructed the SFH by fitting the optical spectrum and found that J-138717 started to form its stars at $t_{20}^{\mathrm{lbt}}=0.8\pm 0.1$ Gyr, then built half of its stellar mass by $t_{50}^{\mathrm{lbt}}=0.6\pm 0.1$ Gyr, and quenched at $t_{90}^{\mathrm{lbt}}=0.4\pm 0.1$ Gyr. This implies a formation timescale of about Δt = 0.4 Gyr, which corresponds to an average $\langle \text{SFR}\rangle =88$ M_⊙ yr⁻¹. The fits of the NIR spectrum provide an earlier $t_{20}^{\mathrm{lbt}}=1.2\pm 0.2$ Gyr, however, which might suggest an older population. This would imply a larger Δt and a lower $\langle \text{SFR}\rangle$.

• The low flux ratios log₁₀([OIII]/[OII])≃−0.1 dex and log₁₀([NeIII]/[OII])≃−0.5 dex suggest that these emissions are due to some residual low-level star formation. The log₁₀([NII]/Hα)≃0.5 dex indicates LINER emission, which might be due either to a weak AGN or to post-AGB stars. Furthermore, the symmetry of the emission and absorption lines observed at their systemic velocity argues against the presence of ongoing outflows.

• The analysis of the C₂4668 index indicated that J-138717 is enhanced in C, but the possible presence of a stellar population older than 1 Gyr is degenerate with it. From the FIF and the Mgb-Fe3 model grids, we derived [α/Fe] = 0.15±0.10 and 0.19±0.13 dex, respectively. This suggests that J-138717 is enhanced in α, but a higher S/N is needed to confirm this result.

• We measured the CO and CN absorption lines by tracing TP-AGB stars and compared them with different model predictions. We found that the EMILES, XSL, and CvD18 models, which are based on empirical stellar libraries and have a low contribution of TP-AGB stars, match a large number of these indices. The M13 models have a mild contribution of TP-AGB stars and match several indices assuming a younger age (0.6 Gyr), which is consistent with the fits. Finally, the M05 models, having a heavy contribution of TP-AGB stars, do not match any of the CO or CN indices. These results indicate that the contribution of TP-AGB stars to the spectrum of J-138717 is not significant.

• We measured several other NIR indices and compared them with the model predictions. We found that the M13 models can match almost all of them. EMILES, XSL, and CvD18 also reproduce many of these indices. In all cases, we inferred inconsistent metallicities that varied significantly from one index to the next and from one library to the next. Finally, the M05 models are only consistent with some NIR spectral indices. We highlight that a detailed analysis of the NIR spectral indices like this was performed outside the local Universe for the first time.

The SFH suggests that J-138717 formed most of its mass quickly and then quenched rapidly. The current data prevent us from clearly determining what caused its quenching. The analysis of the CO and CN spectral indices suggested that TP-AGB stars do not account for a significant fraction of stellar luminosity, which is in contrast to what is expected from theoretical predictions for young PSBs such as J-138717.

In general, it is reassuring to find consistent estimates of the stellar population properties from the full-spectral fitting when different wavelength ranges are fit and different model libraries are used. On the other hand, all models struggle to reproduce the NIR spectral indices consistently, which indicates the importance of improving the current stellar population models in the NIR, especially at young ages and low metallicities. This is highly relevant to studying high-redshift galaxies in the JWST era.

Data availability

All the data from the DR3 are publicly available and can be downloaded from the MAST archive at https://archive.stsci.edu/hlsp/jades. The B24 spectra are publicly available and can be retrieved from the Dawn JWST Archive.

Acknowledgments

We thank the referee for the constructive report. D.B. acknowledges support by the Archival Research Visitor Programme of ESA. D. B. is thankful to Claudia Maraston for kindly providing the M13 models and useful comments on the manuscript, and to Rashmi Gottumukkala for support with the B24 spectra. D.B., P.S., R.D.P., F.L.B., A.G., A.P., C.S., and S.Z. acknowledge support by the grant PRIN-INAF-2019 1.05.01.85. F.L.B. acknowledges support from the INAF grant 1.05.23.04.01. This work made use of and acknowledge the following software: NumPy (Harris et al. 2020); SciPy (Virtanen et al. 2020); Astropy (Astropy Collaboration 2022); matplotlib (Hunter 2007); EAZY (Brammer et al. 2008); pPXF (Cappellari & Emsellem 2004; Cappellari 2017, 2023); FAST (Kriek et al. 2009); emcee (Foreman-Mackey et al. 2013).

References

Appendix A: DR3 vs B24

A.1. Comparison of the observed spectra

In Fig. A.1 we show the comparison between the data reduction of the JADES DR3 and that of B24. In the upper panel we show the comparison of the prism spectra from the two reductions. In this case, there is a clear difference on the slope, varying from about a +20% to −20% between the two data reduction, with a peak of a 50% excess in the DR3 prism spectrum at rest-frame wavelengths below ∼3000 Å. As shown in the next section, this difference leads to different estimates of the stellar population parameters.

Fig. A.1. Comparison of the G140M (upper panel), G395M (middle panel), and prism (lower panel) spectra for the data reduction of the JADES DR3 (orange) and B24 (blue). The black lines are the flux ratios FDR3/FB24−1, the gray dotted lines and shaded regions are the median values and the quadratic errors on the flux, respectively, while the magenta line marks the zero.

In the middle panel, we show the G140M. The two spectra are very similar, with a median difference in the flux close to zero, and they are consistent within the errors all over the spectrum, except for the NUV (observed λ<1.2μm), where the noise dominates. All the measured spectral indices in this wavelength range are also consistent when measured from the DR3 and B24. As an example example, in the upper panel of Fig. A.2 we show the CaII-HK index.

Fig. A.2. CaII (H&K) index (left panel) and CO1.64 index (right panel) for the two data reductions, DR3 (blue) and B24 (orange). The measurements of the former (latter) index are consistent (inconsistent) when measured from the two data reductions. In both panels, the white central region is the wavelength region where the feature is measured, while blue and red regions are the blue and red continuum bandpasses, respectively.

Finally, in the lower panel, we show the G395M spectrum. In this case, the two spectra are consistent from ∼2.9−4μm (∼1−1.4μm rest-frame), with a median systematic offset of only ∼5%. At longer wavelengths, however, the two spectra are clearly different, as the flux of the B24 spectrum is significantly higher than that of the DR3, increasing with the wavelength, and up to ∼20% higher. We note that this difference is very similar to the one observed in the prism. Such differences do not affect only the continuum, but also the single absorption features. For example, in the lower panel of Fig. A.2 we show the CO1.64 index of the two reductions. Here, there is not only a systematic offset in the flux, but the shape of the feature, as well as of the side bands is clearly different, leading to inconsistent measurements of the same index.

A.2. Stellar population parameters

In Fig. A.3 we show the comparison between the stellar population parameters estimated for the DR3 and B24 prism spectra, using the same methods described in section 3.2 and different SSPs libraries (see also Appendix B). The different slopes of the DR3 and B24 spectra generally provide different estimated stellar population parameters. Generally, the B24 spectrum provides systematically older ages (up to 0.5 Gyr) and higher metallicities (up to 0.2 dex). This is consistent with the B24 spectrum being flatter than the DR3 spectrum. The dust attenuation parameters, instead, are generally more consistent.

Fig. A.3. Comparison of the different stellar population parameters estimated when fitting the DR3 and B24 prism spectra.

We verified that the G140M (fit in the wavelength range 3600−4400 Å) and G395M (fit in the wavelength range 1−1.4 μm) provide consistent results for all parameters. This is not surprising, given the consistency of the two spectra in the wavelength regions considered (Fig. A.1).

Appendix B: Stellar population parameters with different SSPs libraries

We fit the stellar population parameters using the models libraries described in section 2.2 (see also Table 1). Since the M05 and M13 models have a lower resolution than the G140M and G395M spectra, we lower the resolution of the latter to match that of the models before performing the fit. We note that different models adopt different isochrones, having different solar metallicities, so we rescaled all the metallicities to Z_⊙ = 0.0198 after the fit. In Fig. B.1 we show the results of the fits using different models libraries. In general, the estimates are consistent with those of the EMILES models. Considering all the different wavelength ranges and models libraries we fit, we estimate averages age 0.94 Gyr, metallicity −0.39 dex, A_V 0.65 mag, and δ -0.8 dex, consistent with those reported in 3.2. We note that the M13 models estimate ages ∼0.6 Gyr in all cases, that are the youngest values among the libraries. In section 4.2 and 4 we show the M13 can match better many indices assuming a younger age for J-138717. We also note that the XSL models provide typically older ages. Finally, we note that the CvD18 models provide solar metallicity when fitting the prism and NIR spectrum, although this could be due to the lack of younger models; this could also explain the high value of δ measured from the prism.

Fig. B.1. Comparison of the stellar population parameters estimated using different SSP libraries

Appendix C: Other NIR spectral indices

In this appendix we show the NIR spectral indices discussed in section 4.3. Fig. C.1 we show the Paschen lines and in Fig. C.2 we show some other NIR spectral indices.

Fig. C.1. Paschen series indices as a function of age. The marker and the lines are the same as in Fig. 7. Predictions from M05 and M13 for Paδ are not shown because this index cannot be measured safely due to the low spectral resolution (see section 4.3 for details).

Fig. C.2. Other NIR indices as a function of age. The marker and the lines are the same as in Fig. 7.

All Tables

All Figures

Fig. 1. Spectrum of the JWST of J-138717. The rest-frame wavelengths are calculated with respect to the spectroscopic redshift, z = 1.8845 (see sect. 3.1). The combined RGB imaging of the galaxy is shown in the upper right corner. The observed spectra are shown in the upper panels. The black lines show the medium-resolution G140M/F070LP and G395M/F290LP spectra, while the overplotted orange line shows the prism spectrum (we note that the prism spectrum also covers the rest-frame interval ∼6500−10 000 Å; see, e.g., Fig. A.1). The two lower panels show the comparison of the flux of the prism and grating spectra (black line). In both panels, the magenta line marks a difference of zero, and the gray lines and shadow regions indicate the median and the average errors, respectively. The vertical red line in the left panels marks λ = 1.27 μm, above which the spectrum is potentially affected by second-order contamination. In all panels, the vertical dashed lines correspond to different absorption and emission lines.

In the text

Fig. 2. Best-fit of the prism (upper panel), G140M (middle panel), and G395M (lower panel) of the EMILES models. The three spectra were fit separately. In each panel, the black line is the rest-frame observed spectrum. The red line is the best-fit to the stellar continuum. The orange line is the best-fit to the emission lines. The gray shaded regions are excluded from the fit. The green diamonds are the residuals and the blue horizontal line is their median value.

In the text

	Fig. 3. Zoom-in of the best fits for the emission lines. The lines and the symbols are the same as in Fig. 2. The magenta line in the lower right panel show the best-fit templates for the Hα and [NII] lines separately.
In the text

Fig. 4. Comparison of the observed indices (black dots show the spectral pixels) with the predictions of the FIF best-fit models (magenta line), with an age of 0.81 Gyr and a metallicity of −0.19 dex. The blue and red regions correspond to the wavelengths where the pseudo-continua are measured and the white regions are the wavelengths where the index is measured. In the bottom panels the residuals are plotted as green diamonds. The gray lines and shaded regions are the median of the residuals and the average observed errors on the indices, respectively, while the black dashed line marks the 0.

In the text

	Fig. 5. FAST best fit to the photometric data (red circles). Specifically, we used the convolved Kron magnitudes in the JWST filters F182M, F210M, and F444W, and the HST filters F606W, F775W, F814W, F850LP, as listed in the JADES photometric catalog.
In the text

	Fig. 6. Star formation history inferred from the regularized fit of the optical spectrum. Left panel: Map of the weights assigned by pPXF to each input template of given age and metallicity. Right panel: Cumulative sum of the weights. The green, blue, and red line correspond to $t_{20}^{\mathrm{lbt}},t_{50}^{\mathrm{lbt}},\text{and}{t}_{90}^{\mathrm{lbt}}$, respectively.
In the text

Fig. 7. C24668 spectral index as a function of age. The yellow star is the measured value from the G140M spectrum of J-138717. The green, red, gray, blue, and purple lines are the predictions of the M05, M13, XSL, CvD18, and EMILES models, respectively. Dark to light colors indicate higher to lower metallicities. Both observations and models have been smoothed to the instrumental resolution of the BaSeL stellar library. The brown (olive) line in the central (right) panel indicates by how much the tracks would shift towards higher index values due to a C-enhancement(α-enhancement) of 0.15 dex (0.4 dex), as predicted by the CvD18 (MILES) models.

In the text

	Fig. 8. CN indices as a function of age. The marker and the lines are the same as in Fig. 7.
In the text

	Fig. 9. Observed profile of the NaD absorption line (black line). The dashed red lines indicate the expected wavelengths of the doublet at λ = 5890 Å and λ = 5896 Å.
In the text

	Fig. A.1. Comparison of the G140M (upper panel), G395M (middle panel), and prism (lower panel) spectra for the data reduction of the JADES DR3 (orange) and B24 (blue). The black lines are the flux ratios FDR3/FB24−1, the gray dotted lines and shaded regions are the median values and the quadratic errors on the flux, respectively, while the magenta line marks the zero.
In the text

Fig. A.2. CaII (H&K) index (left panel) and CO1.64 index (right panel) for the two data reductions, DR3 (blue) and B24 (orange). The measurements of the former (latter) index are consistent (inconsistent) when measured from the two data reductions. In both panels, the white central region is the wavelength region where the feature is measured, while blue and red regions are the blue and red continuum bandpasses, respectively.

In the text

	Fig. A.3. Comparison of the different stellar population parameters estimated when fitting the DR3 and B24 prism spectra.
In the text

	Fig. B.1. Comparison of the stellar population parameters estimated using different SSP libraries
In the text

	Fig. C.1. Paschen series indices as a function of age. The marker and the lines are the same as in Fig. 7. Predictions from M05 and M13 for Paδ are not shown because this index cannot be measured safely due to the low spectral resolution (see section 4.3 for details).
In the text

	Fig. C.2. Other NIR indices as a function of age. The marker and the lines are the same as in Fig. 7.
In the text