© ESO 2022

1. Introduction

Stellar population models are fundamental in determining the basic properties of unresolved stellar systems. Those properties include the initial mass function (IMF), star formation rate, star formation history (SFH), total mass in stars, and stellar metallicity and abundance patterns (see the review by Conroy 2013). With next-generation wide-field spectroscopic facilities, such as the upcoming WEAVE for the William Herschel Telescope (Dalton et al. 2012, Jin et al., in prep.), MOONS for the Very Large Telescope (VLT; Cirasuolo et al. 2020), and 4MOST for the Visible and Infrared Survey Telescope for Astronomy (de Jong et al. 2019), spectroscopic information of different types of galaxies in various environments will increase in both quantity and quality. Furthermore, with recent advances in near-infrared (NIR) instrumentation on large telescopes, such as X-shooter (Vernet et al. 2011) and KMOS (Sharples et al. 2004, 2013) on ESO’s VLT or the forthcoming HARMONI on the ELT (Thatte et al. 2016), the domain of evolved cool stars in stellar populations will be increasingly accessible. Stellar spectral libraries and associated stellar population models need to keep up with these developments.

An increasing effort has been put into developing better stellar population models that are based on empirical stellar libraries. The goals are to build models with higher resolution and longer wavelength ranges based on stellar spectral libraries to cover the Hertzsprung–Russell diagram (HR diagram henceforth) more extensively than ever before. For example, the widely used UV–IR Bruzual & Charlot (2003) models are still based on theoretical spectra across large wavelength regions, while MILES stellar population models (Vazdekis et al. 2010, 2015), which are based on the fully empirical MILES library (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011), have been extended towards the NIR and UV over the years, resulting in the extended-MILES (E-MILES) models, which have a wavelength coverage of 1680–50 000 Å (Vazdekis et al. 2012, 2016; Röck et al. 2016). The recent SDSS MaStar stellar population models (Maraston et al. 2020), with a wavelength coverage of 3600–10 300 Å, are based on nearly 9000 stars, a tenfold increase on the previous generation of models, although they cover only the optical wavelength range. With these modern stellar population models, it is now common to find good agreement between the observed and the predicted NUV and optical properties of stellar clusters (e.g Peacock et al. 2011; Ricciardelli et al. 2012; Conroy et al. 2018). This consensus shows our good understanding of the main-sequence and early giant phases that constitute the near-ultraviolet (NUV) and optical light of stellar populations.

However, the existing optical-to-NIR stellar population models have problems. The NIR traces populations with a range of ages and suffers lower dust extinction than the optical. However, we are far from a complete understanding of some stellar evolutionary phases that strongly affect the spectral energy distributions (SEDs) of stellar populations in the NIR (e.g. Mouhcine & Lançon 2002; Vazdekis et al. 2016; Baldwin et al. 2018; Riffel et al. 2019). The asymptotic giant branch (AGB), especially the thermally pulsing AGB (TP-AGB), shapes the NIR spectra of 0.5–2 Gyr old stellar populations; the tip of red giant branch (RGB) defines the NIR spectra of older populations. Current stellar population models in the NIR are based on available empirical libraries, such as Pickles (1998), Lançon & Wood (2000), and (E-)IRTF (see below), or on theoretical stellar spectra, such as MARCS (Gustafsson et al. 2008), PHOENIX (Husser et al. 2013), or BaSEL (Lejeune et al. 1997, 1998; Westera et al. 2002) models. The IRTF Spectral Library of Rayner et al. (2009) and the extended-IRTF (E-IRTF) of Villaume et al. (2017) are empirical libraries of 0.8–5.0 μm and 0.7–2.5 μm (respectively) stellar spectra observed at a resolving power of R = 2000 with the SpeX spectrograph at the NASA Infrared Telescope Facility on Mauna Kea. The original IRTF library covers mainly solar-metallicity late-type stars (but also some oxygen-rich and carbon-rich AGB stars); the E-IRTF expands the metallicity coverage. The E-MILES, Conroy et al. (2018), and Meneses-Goytia et al. (2015) stellar population models take advantage of either the IRTF or E-IRTF library. The empirical library of Lançon & Wood (2000) has a spectral resolution R ∼ 1000 and is limited to cool giant and supergiant stars only. Mouhcine & Lançon (2002) and Maraston (2005) have included these spectra in their stellar population models.

None of these empirical libraries have extensive coverage of the important stellar evolutionary stages needed for stellar population modelling in the NIR. Furthermore, (O- and C-rich TP-AGB) and RGB stars are rarely segregated in stellar population modelling. This leads to a large variety of optical-to-NIR stellar populations, which in turn leads to discrepancies between the SFHs derived from optical and NIR spectral ranges, or from different models. An example is the Maraston (2005) set of simple stellar population (SSP) models, which have enhanced flux and strong molecular carbon and oxygen absorption features throughout the NIR spectra of intermediate-age populations compared to, for example, the Bruzual & Charlot (2003), E-MILES, and Conroy et al. (2018) models. Such strong molecular bands predicted by the Maraston (2005) models have been detected in some studies (e.g Lyubenova et al. 2012), but not in others (e.g Zibetti et al. 2013). Recent works in stellar evolution theory (Girardi et al. 2013; Pastorelli et al. 2020) explain this observational discrepancy with the ‘AGB boosting’ effect, which is linked to the physics of stellar interiors – stellar populations in a narrow 1.57 and 1.66 Gyr age range at Magellanic Cloud metallicities have a factor of ∼2 increase in the TP-AGB contribution to the integrated luminosity of the stellar population. Some of the Lyubenova et al. (2012) globular clusters of the Magellanic Clouds are in this very narrow age and metallicity range and show clear spectral features of TP-AGB stars, while the post-starburst galaxies of Zibetti et al. (2013) are probably not within this range. Besides the inclusion of the TP-AGB phase in the stellar population models, the overall quality and coverage of the stellar spectral library is important. Baldwin et al. (2018) found that the largest differences in derived SFHs are caused by the choice of stellar spectral library and suggested that the inclusion of high-quality NIR stellar spectral libraries in stellar population models should be a top priority for modellers.

Furthermore, theoretical stellar spectra cannot be used at present to make accurate predictions for the NIR spectra of stellar populations, as they have considerable problems in reproducing SEDs and molecular bands of observed cool stars (e.g. Martins & Coelho 2007; Kurucz 2011; Coelho 2014; Knowles et al. 2019; Martins et al. 2019; Coelho et al. 2020; Lançon et al. 2021). These stars are very difficult to model due to processes such as hot-bottom burning, stellar winds, long-period pulsations, the presence of circumstellar dust, and the third dredge-up (Pastorelli et al. 2019, 2020).

Another common limitation of existing stellar population models based on empirical stellar libraries is their low spectral resolution. Stellar population models that are based on empirical stellar libraries typically have resolutions of R ∼ 2000. For example, the commonly used spectral-line index system (LIS) of Vazdekis et al. (2010) suggests using LIS-5.0Å (R ∼ 1000 at the Mg triplet at 5170 Å, which corresponds to a velocity dispersion of 127 km s⁻¹) to study low-velocity-dispersion systems such as globular clusters or dwarf galaxies. Higher spectral resolution is required for more detailed modelling of emission and absorption lines; for example, higher-resolution spectroscopy can provide more accurate measurements of numerous absorption lines for many different chemical elements. Notable high-resolution empirical stellar populations model are the Pegase.HR stellar population models (Le Borgne et al. 2004, 2011), with R = 10 000 over the 4000–6800 Å wavelength range, which are based on the ELODIE stellar spectral library (Prugniel & Soubiran 2001, 2004; Prugniel et al. 2007). For example, Şen et al. (in prep.) defined a set of line indices with a resolution of σ = 25 km s⁻¹ using the Pegase.HR models and determined abundance ratios of 11 elements in small, unresolved galaxies outside the Local Group. They found that the majority of their dwarf galaxies have abundance ratios slightly less than solar.

Here we present SSP models¹ based on 639 stellar spectra from the X-shooter Spectral Library (XSL) data release 3 (Verro et al. 2022, DR3). This new library is designed for stellar population purposes, with unprecedented simultaneous wavelength coverage of 3500 Å–2.48 μm with a resolution of σ ∼ 13 km s⁻¹ (corresponding to R ∼ 10 000). XSL aims to cover the entire HR diagram as extensively as possible, with an emphasis on the advanced stellar evolutionary stages. We incorporated the spectra of 44 oxygen-rich (quasi-)static stars cooler than 4000 K, 39 oxygen-rich TP-AGB stars, and 26 spectra of carbon-rich TP-AGB stars into our new stellar population models. With this development, the XSL SSP models will help us to bridge the optical and the NIR studies of intermediate-age and old stellar populations and clarify the role of evolved cool stars in stellar population synthesis. The moderately high resolution of XSL population models creates new possibilities in the optical to NIR for absorption-line index studies.

This paper is structured as follows: We review the main ingredients for stellar population models in Sect. 2 and describe the model calculation in Sect. 3. We describe and analyse the general behaviour of the new XSL SSP models in Sect. 4. We compare colours and absorption-line indices with observed galaxies in Sects. 5 and 6, respectively. Furthermore, we provide the stellar mass-to-light ratios in Sect. 7 and further comment on the effects of cool giant stars on the population models in Sect. 8. Finally, in Sect. 9 we define the regions in age and metallicity where the XSL stellar population models are safe to use.

Throughout this paper, the UBVRI magnitudes are in the Johnson-Cousins system, and JHK magnitudes are in the homogenised Bessell system (Bessell & Brett 1988) (Vega system).

2. Main ingredients for stellar population models

The construction of SSP models is rather straightforward, as it consists of only three ingredients – stellar evolution theory (isochrones), an IMF, and a stellar spectral library. These ingredients are typically combined by Eq. (1):

$$\begin{array}{cc}\hfill f_{\mathrm{SSP}}(t,[\mathrm{Fe}/\mathrm{H}])& ={\displaystyle {\int }_{m_{\mathrm{low}}}^{m_{\mathrm{high}}(t)}{f}_{\ast }\left[T_{\mathrm{eff}}(M),logg(M)|t,[\mathrm{Fe}/\mathrm{H}]\right]}\hfill \\ \hfill & \times \mathrm{\Phi }(M)dM,\hfill \end{array}$$(1)

where M is the initial stellar mass, Φ(M) is the IMF, f_* is the spectrum of a star of mass M of effective temperature T_eff and surface gravity log g at metallicity [Fe/H], and f_SSP(t, [Fe/H]) is the resulting spectrum of a stellar population of a certain age (t) and metallicity [Fe/H]², and the integration runs from the lowest stellar mass in the IMF, m_low, to the highest stellar mass still living at time t, m_high(t). The nuances of the ingredients themselves are what make population modelling difficult in practice. We recommend the review by Conroy (2013) for an overview of this broad topic. We discuss the specific choices for the SSP models in this work below.

2.1. Isochrones

Due to the extension of these SSP models to the NIR, the advanced evolutionary stages of stars become extremely important. XSL contains a large number of evolved cool giants, which makes synthesising realistic stellar populations in the NIR possible, as long as we know how to integrate them into the models.

An isochrone determines the location of stars with the same age and metallicity on the HR diagram and is constructed from stellar evolution calculations. On one hand, our selection of isochrones is motivated by the thorough treatment of the advanced evolved stages. With that in mind, we used the PARSEC/COLIBRI isochrones. The PARSEC version 1.2S models (Bressan et al. 2012; Tang et al. 2014; Chen et al. 2014a, 2015) describe the evolution of stars from pre-main-sequence stars to the first thermal pulse in the helium shell, after forming an electron-degenerate carbon-oxygen core. Then the COLIBRI models (Marigo et al. 2013; Rosenfield et al. 2016; Pastorelli et al. 2019, 2020) add the TP-AGB evolution, from the first thermal pulse to the total loss of envelope. These isochrones have the most advanced handling of TP-AGB stars to date, based on high-quality observations of resolved stars in the Small Magellanic Cloud with detailed stellar population synthesis simulations computed with the TRILEGAL code (Girardi et al. 2005). On the other hand, we aim to calculate stellar population models from simpler and more widely used stellar evolution tracks as well. The Girardi et al. (2000, Padova00 henceforth) isochrones allow us to compare directly our models with the E-MILES stellar population models. These tracks include a simple but synthetic treatment until the tip of the AGB, but they do not include a third dredge-up. Therefore, there is no transition from O-rich to C-rich TP-AGB stars in the tracks, and so these isochrones are missing these stars.

2.2. IMF

An IMF describes the initial distribution of masses for a population of stars formed at the same time. XSL stellar population models are calculated using a Salpeter (1955) or a Kroupa (2001) IMF. The Salpeter IMF is a single power law with an α = 2.35 slope, and is valid for 0.4 < m/M_⊙ < 10. The Kroupa IMF is a double power law, with α = 1.35 slope for m/M_⊙ < 0.5 stars, and α = 2.35 for higher-mass stars. In both cases, we use the relation in the mass range 0.09 < m/M_⊙ < 120, using extrapolation to lower masses. More XSL SSP models, calculated with various IMFs, will be presented and discussed in Verro et al. (in prep.).

2.3. The X-shooter Spectral Library (XSL)

XSL is a moderate-resolution (R ∼ 10 000) NUV–NIR stellar spectral library intended for stellar population modelling. We are using the XSL DR3 data to construct stellar population models. In Verro et al. (2022), we provided 830 spectra of 683 stars, which are corrected for galactic extinction and merged to the full wavelength range of X-shooter, 350–2480 nm. The data were homogeneously reduced and calibrated in Gonneau et al. (2020). XSL spectra are given in rest-frame, at a resolution σ = 13,  11,  16 km s⁻¹ in the UV-Blue (UVB), visible (VIS), and NIR arms, respectively (Gonneau et al. 2020). Arentsen et al. (2019) provided a uniform set of stellar atmospheric parameters – effective temperatures, surface gravities, and metallicities – for 754 spectra of 616 XSL stars. We adopted these stellar parameters for the DR3 spectra. Our sample has many stars with multiple observations. We regarded these observations as separate stars with slightly different stellar parameters, as determined in Arentsen et al. (2019).

Not all stars in DR3 are useful for stellar population modelling. With the exception of the red giants, we selected only non-variable non-peculiar stars with complete X-shooter spectra from the DR3 dataset. We excluded spectroscopic binary stars. We only included XSL spectra that have not been corrected for slit flux losses and galactic extinction when constructing the ‘static’, O-rich TP-AGB, and C-rich TP-AGB star sequences in Sect. 3.2, otherwise we use spectra that are corrected. We only included DR3 spectra for which stellar parameters have been estimated in Arentsen et al. (2019) or in Verro et al. (2022). Furthermore, as described in Sect. 3.2, we removed supergiants from the library because we aim to model stellar populations older than 50 Myr, in which supergiants do not occur. These selections resulted in 639 spectra of 534 stars (from the 830 stellar spectra of 683 stars of DR3), which were used to create the XSL stellar population models.

3. Model calculation

3.1. Spectral interpolator

Each point on the isochrone needs a representative stellar spectrum, when generating an SSP model. The limited coverage of the HR diagram by empirical libraries requires a method to assign the stars in the library to the isochrones. Commonly, an interpolator is used to do this. An interpolator creates a synthetic spectrum at a given set of parameters (e.g. T_eff, log g and [Fe/H]) from a library of empirical or theoretical spectra. A local interpolator (e.g. Vazdekis et al. 2003; Sharma et al. 2016; Dries 2018) interpolates spectra using its local neighbourhood: library stars in the vicinity of the point for which we want to create a spectrum are weighted and combined to create a representative spectrum for that point. A global interpolator (e.g. Prugniel & Soubiran 2001; Koleva et al. 2009; Prugniel et al. 2011; Wu et al. 2011) fits polynomials of T_eff, log g and [Fe/H] at each wavelength point to the whole or a large subset of the spectra in the library. Here we use a combination of the two in different areas of the HR diagram. Moreover, static giants, O-rich TP-AGB stars, and C-rich TP-AGB stars are treated separately.

3.1.1. Global interpolation

The global interpolator consists of polynomial expansions for each wavelength pixel in powers of the three stellar parameters. This type of interpolator was first introduced in Prugniel & Soubiran (2001) and used with the ELODIE stellar library in Prugniel & Soubiran (2001) and Wu et al. (2011) and the MILES stellar library in Prugniel et al. (2011) and Sharma et al. (2016). The polynomial coefficients β_i, i = 0, 1, …, 25, fitted for each wavelength pixel λ over the subset of spectra are as follows:

$$\begin{array}{cc}& Y(x,y,z,\lambda )={\beta }_0(\lambda )+{\beta }_1(\lambda )·x+{\beta }_2(\lambda )·z+{\beta }_3(\lambda )·y\hfill \\ \hfill & +{\beta }_4(\lambda )·x^2+{\beta }_5(\lambda )·x^3+{\beta }_6(\lambda )·x^4+{\beta }_7(\lambda )·x·z\hfill \\ \hfill & +{\beta }_8(\lambda )·x·y+{\beta }_9(\lambda )·x^2·y+{\beta }_{10}(\lambda )·x^2·z+{\beta }_{11}(\lambda )·y^2\hfill \\ \hfill & +{\beta }_{12}(\lambda )·z^2+{\beta }_{13}(\lambda )·x^5+{\beta }_{14}(\lambda )·x·y^2+{\beta }_{15}(\lambda )·y^3\hfill \\ \hfill & +{\beta }_{16}(\lambda )·z^3+{\beta }_{17}(\lambda )·x·z^3+{\beta }_{18}(\lambda )·y·z+{\beta }_{19}(\lambda )·y^2·z\hfill \\ \hfill & +{\beta }_{20}(\lambda )·y·z^2+{\beta }_{21}(\lambda )·x^3·y+{\beta }_{22}(\lambda )·x^4·y\hfill \\ \hfill & +{\beta }_{23}(\lambda )·x^3·z+{\beta }_{24}(\lambda )·x^2·y^2+{\beta }_{25}(\lambda )·x^2·y^3,\hfill \end{array}$$(2)

where x = log T_eff, y = log g and z = [Fe/H]. A weighted linear least squares solution is found for each wavelength pixel. The weight of each XSL spectrum in the sum of squared differences is made up of two factors: one being the signal-to-noise ratio (S/N) of the spectrum and the other reflecting how isolated this XSL spectrum is in parameter space. The latter is described by the inverse density of XSL spectra in a box of size (T_eff ± 2500 K, log g ± 1.5 dex, [Fe/H] ± 1.0 dex) surrounding the desired parameter. We assigned a uniform weight of S/N = 10 to spectra that have been corrected for slit flux losses in DR3 with a polynomial³. This interpolator assumes a smooth transition of spectra in the stellar atmosphere parameter space. Using information from a large subset of stars, an interpolated spectrum is little affected by issues of individual stars (e.g. poor flux calibration, dichroic contamination, inaccurate stellar parameter estimation, extinction correction issues, or XSL arm-merging inaccuracies). We use this type of interpolation in the ‘warm’ star regime (4000 K < T_eff < 8000 K; see Sect. 3.3 and Fig. 1 for more details)

Fig. 1. 2 Gyr, [Fe/H] = −0.6 PARSEC/COLIBRI isochrone and the locations on the HR diagram where spectra are generated by local interpolation, global interpolation, or potentially taken from the ‘static’ giant, O-rich TP-AGB, or C-rich TP-AGB sequences. We only switch to the respective sequences when we have reached the bluest average spectrum on that sequence (according to the colour–temperature relation). For example, only the coolest RGB stars (Teff ⪅ 4000 K) are represented by a spectrum originating from the static sequence, and the spectra of warmer RGB stars are created by the global interpolator.

3.1.2. Local interpolation

Unfortunately, global interpolation fails in poorly covered parameter-space regions of the library, including at the edges of the parameter space, due to the use of polynomials. In these regions, the local interpolation comes to assist. The local interpolator averages stellar spectra in a box of parameters around the desired point, so it works better in lower density regions of the HR diagram. We use this type of interpolation in the cool dwarf (T_eff < 4500 K, log g > 4.0 dex) and hot star (T_eff > 7000 K) regime (see Sect. 3.3 and Fig. 1 for more details).

The local interpolation combines weighted spectra in eight cubes in the stellar parameter space, all with one corner at θ₀ (≡5040/T_eff, 0), log g₀, [Fe/H]₀. The initial size of each three-dimensional cube of Δθ₀, Δlog g₀ and Δ[Fe/H]₀ is 3σ_θm × 3σ_{log gm} × 3σ_[Fe/H]m, where σ_paramm corresponds to the minimum uncertainty in the determination of the respective stellar parameters. Following stellar parameter estimations from Arentsen et al. (2019), we adopted the following values as uncertainties [σ_θm, σ_{log gm}, σ_[Fe/H]m]: hot stars – [0.018, 0.20, 0.1]; cool dwarfs – [0.012, 0.14, 0.1]. If no stars are found, the box is enlarged along each of its axes in steps of 1σ, until at least one star is found. This interpolation scheme is described in detail in Vazdekis et al. (2003, Appendix B) and Dries (2018).

3.1.3. Interpolation quality

To test the local and global interpolator, we created an interpolated spectrum for each star in XSL DR3 dataset (excluding cool giants, which are incorporated into the models separately), in such a way that the original XSL star is excluded from the dataset that we use to build the interpolator. We calculated the median residual R_S between the original S_or and the interpolated spectrum S_int:

$$\begin{array}{c}\hfill R_{\mathrm{S}}=\mathrm{median}\left(\frac{\mathrm{abs}(S_{\mathrm{or}}-S_{\mathrm{int}})}{S_{\mathrm{or}}}\right).\end{array}$$(3)

The main telluric bands, as well as the XSL dichroic areas (shown in Table 1) are masked out when calculating R_S. The median residuals R_S for the XSL stars used in stellar population modelling though the usage of the global and local interpolators are shown in Fig. 2 as a function of positions on the HR diagram for the full wavelength range, and for the X-shooter UVB, VIS, and NIR arms separately. Figure A.1 shows a similar plot for the median residuals around four spectral line indices (CaHK, Hβ, CaT, and CO1.6). We note, however, that we only discuss the incorporation of the very cool giant stars to the models in Sect. 3.2; therefore the very cool giants are missing from Fig. 2.

Fig. 2. Weighted median residual, RS, between the original spectrum and the interpolated spectrum for the full wavelength range, and for the X-shooter UVB, VIS, and NIR arms separately, as a function of position on the HR diagram. The colour bar is logarithmic. Histograms show the distributions of RS calculated within these spectral ranges at the full XSL resolution. For ease of visualisation, spectra with RS > 0.1 are placed into the RS = 0.1 bin in the histograms.

The median residuals should be taken as a rough quality measure, considering the variety of spectral types and long wavelength range of XSL – cool stars have relatively less signal (lower S/N) in the UVB than hot stars and hot stars have relatively less signal in the NIR than cool stars. A large mismatch between an interpolated spectrum and the corresponding original spectrum may be the result of very low S/N. Large median residual values can also indicate poorer reproduction of the star by the interpolation due to uncertain stellar parameters, extinction correction, DR3 merging factors, peculiarity of the spectrum, residual telluric lines, or due to the interpolation scheme itself. We expect poorer matches at the edges of the parameters space and areas with low density, due to the nature of the interpolators. As seen from Fig. 2, the residuals in all wavelength ranges are mostly less than 5%. For the VIS and NIR spectra, the median residuals are of the order of a few percent. For the regions around Hβ, CaT, and CO1.6 lines in the VIS and NIR arms, the median residuals are also of the order of a few percent, while region around the CaHK line in the UVB arm has higher residuals. The UVB-arm spectrum is the most difficult to interpolate for most spectral types, due to the multitude of spectral features compared to the VIS and NIR arms and the continuum shape changes rapidly with stellar parameters and cool stars have near-zero fluxes in the UVB arm. The CaHK line (3800–4000 Å) is in the region of the UVB spectrum where these difficulties are most evident. We gave some examples of XSL DR3 and its interpolated counterpart in Verro et al. (2022), where we used the same interpolation scheme.

3.2. Incorporating cool evolved stars

In the XSL stellar population models, we give extra attention to the cool (T_eff < 4000 K) evolved giants. We divided these stars into (O-rich) static giants, O-rich TP-AGB stars, and C-rich TP-AGB (‘carbon’) stars. Individual observed spectra of such stars vary strongly in their (NIR) SEDs as a function of time, and from star to star, so they cannot be used directly in the synthesis of galaxy spectra. An ideal stellar library should include spectra of individual variable stars observed over their pulsation cycles. In reality, the light curves and phases are generally not accurately known. Furthermore, the stellar parameters are not accurately known. The stellar parameter estimation should be done based on spectral type and temperature-sensitive spectral features. Full-spectrum fitting with theoretical (Lançon et al. 2019, 2021) or interpolated empirical spectra (Arentsen et al. 2019) for these stars is unreliable. Re-evaluating the stellar parameters for these stars and conducting additional observations of stars in different pulsation stages are well beyond the scope of the current paper. Instead, we used the approach described in Lançon & Mouhcine (2002) – using average spectra of static giants, O-rich TP-AGB stars, and C-rich TP-AGB stars, binned by broadband colour, and relying on empirical relations to dictate where an average spectrum of a star of a certain colour should occur – to incorporate these stars into our stellar population models. This method allows us to also use XSL giant stars with T_eff < 4000 K, for which the parameter estimation by Arentsen et al. (2019) is inadequate.

3.2.1. Differentiating between O-rich static giants, supergiants, and variable stars

Differentiating between O-rich static or quasi-static giants, supergiants and variable stars is difficult. The spectra of long period variables with small visual amplitudes are very similar to those of static stars. Supergiants and giants can have similar optical features, although supergiants have redder SEDs (Lançon & Wood 2000; Alvarez et al. 2000). All of these type of stars can have the same broadband colours, so binning by a certain temperature sensitive broadband colour without separating by spectral type first could result in a very red old stellar population model with strong supergiant or TP-AGB features. We separated the (quasi-)static from the high amplitude variables using the (I − K) colour and the H-band H⁻/H₂O feature.

The H⁻/H₂O feature is a combination of the 1.6 μm ‘bump’ in the minimum opacity of H⁻ ion and H₂O vapour absorption bands around 1.4 μm and 1.9 μm. These H₂O bands create the curved shape of the spectrum illustrated in Fig. 3, which is a characteristic feature of long-period variable M stars (Bessell et al. 1989; Matsuura et al. 1999; Alvarez et al. 2000). Although these water bands are contaminated by telluric lines, the overall feature is still distinctive. We created an H⁻/H₂O index to describe the feature, defined in Table 2. We measured the index at the native XSL resolution, in magnitudes, following the index equation of Worthey et al. (1994):

$$\begin{array}{c}\hfill I=-2.5log\left[\left(\frac{1}{{\lambda }_1-{\lambda }_2}\right){\displaystyle {\int }_{{\lambda }_1}^{{\lambda }_2}}\frac{F_{\lambda }}{F_{\mathrm{cont}}}d\lambda \right],\end{array}$$(4)

Fig. 3. Spectral H−/H2O feature. It is a combination of the 1.6 μm ‘bump’ in the minimum opacity of H−ion and H2O vapour absorption bands around 1.4 μm and 1.9 μm. The index bands are marked in red. The telluric absorption is marked in grey. This spectrum is an average of O-rich TP-AGB stars with (I − K) = 3.19.

where F_c is the pseudo-continuum flux defined by drawing a straight line from the midpoint of the blue continuum level to the midpoint of the red continuum level, and F_λ is the flux of the index.

Figure 4 shows the XSL stars on this colour–index plane, colour-coded by log T_eff from Arentsen et al. (2019). The stronger the H-band feature, the more negative the index. While the rest of the stars in the XSL follow a linear relation in this plane, stars redder than (I − K) = 2, corresponding to stars cooler than 4000 K, show a wide variety of H⁻/H₂O index strengths. C-rich TP-AGB stars are plotted in this figure but have been removed when dividing red giants into static and variable stars. C-rich stars are very recognizable due to their distinctive SEDs, with bands of carbon compounds and an absence of oxide bands. XSL C-rich TP-AGB stars have been studied in detail by Gonneau et al. (2016, 2017). The NIR bands of oxygen-rich H₂O and carbon-rich CN and C₂ overlap in wavelength. Carbon stars have strong H⁻/H₂O index strengths due to CN and C₂ absorption, not because of H₂O.

Fig. 4. H−/H2O index strengths of XSL stars as a function of their (I − K) colours. Points are colour-coded by their effective temperatures from Arentsen et al. (2019), saturated at 10 000 K. Carbon stars are marked with circles. Supergiants are marked with diamonds. Both C-rich TP-AGB stars and supergiants are excluded from the division into static and variable stars. M dwarfs have weak H−/H2O index strengths but can have red (I − K) colours, and they have been marked with triangles.

To remove supergiants from this dataset, we used the CN1.10 index defined by Röck (2015). Supergiants display prominent NIR CN absorption (in particular at 1.10 μm), while other O-rich giants do not. However, the 1.15 μm H₂O band, which is also heavily blended with TiO and VO bands in the coolest long-period variables, can be confused with the CN band (Lançon & Wood 2000). Those long-period variables should have strong H⁻/H₂O features, while supergiants should not, and the two can be separated. Furthermore, supergiants may not have a strong CN1.10 feature. We also removed three spectra of stars that are in the Massey (2002) catalogue of supergiants.

3.2.2. Average spectra of static red giants

We selected 44 oxygen-rich stars from Fig. 4, which are (quasi-)static. We created a sequence of average spectra, with each average spectrum consisting of stars with similar (I − K) colour. (I − K) colour is known to correlate with the effective temperature in M stars (e.g Bessell et al. 1989, 1998; Lançon & Mouhcine 2002; Lançon et al. 2019). We combined stars into one bin using weighted averaging (with S/N in the I band as the weight). This ‘static sequence’ is shown in Fig. 5. The selected 44 stars are listed in Table B.1 and shown individually in Fig. B.1.

Fig. 5. Sequence of static giant spectra, sorted according to (I − K) (on the left of each spectrum). The number of spectra the average is composed of is given in the brackets. The spectra have been smoothed to lower resolution for clarity.

3.2.3. Average spectra of O-rich TP-AGB stars

The (I − K) colour is known to correlate with the effective temperature also for TP-AGB stars; (V − K) and (R − K) could be used for this purpose as well (Ridgway et al. 1980; Lançon & Mouhcine 2002). However, TP-AGB stars do not follow a simple colour–temperature relation. Stars with different pulsation properties can be found at the same TP-AGB temperature. When their temperatures decrease, TP-AGB luminosities rise, their radii increase, their masses decrease due to mass loss, and their stellar pulsation properties change (e.g. the DARWIN models for M-type AGB stars Bladh et al. 2019). Therefore, using a simple colour–temperature relation means we assume that the spectrum of an individual variable star, averaged over its cycle, is similar to an average spectrum of many stellar spectra of various masses, amplitudes and phases, but with a common colour-inferred temperature.

We selected 39 XSL spectra as O-rich TP-AGB stars. These stars are listed in Table C.1 and shown individually in Fig. C.1. This selection produces a sequence of O-rich spectra with continually evolving properties, as seen in Fig. 6. We note the deepening of the NIR H-band H⁻/H₂O feature with increasing colour. We call this the oxygen-rich ‘variable sequence’. We note that the reddest average spectrum on this sequence, marked with grey in Fig. 6, consists of only one star, X0145 (OGLEII DIA BUL-SC41 3443), and has a very extreme colour of (I − K) = 5.76. We do not use this star due to it being the only star with such an extreme colour. Moreover, according to the colour–temperature relation of Worthey & Lee (2011), we do not need it (see Sects. 3.3 and 8.2). The average spectrum (I − K) = 4.88 is also a single star, X0020 (ISO-MCMS J005714.4-730121). There are no other stars with such red colours. Because we need a spectrum with such extreme colours, we do use X0020 in the stellar population models.

Fig. 6. Sequence of O-rich TP-AGB spectra, sorted according to (I − K) (on the left of each spectrum). The number of spectra the average is composed of is given in the brackets. The spectra have been smoothed to lower resolution for clarity. The reddest average spectrum on this sequence, displayed in grey, is not used in the stellar population models due to its extreme SED.

3.2.4. Average spectra of C-rich TP-AGB stars

Some O-rich TP-AGB stars will become C-rich through convective dredge-up of newly synthesised carbon from their cores. This third dredge-up is induced by thermal pulses (e.g. Iben & Renzini 1983) and depends on the initial mass and metallicity of these stars. Their spectra differ radically from those of other of cool giants. The spectrum of a C-rich TP-AGB star is characterised by bands of carbon compounds, such as CN and C₂ bands, and by the absence of oxide bands such as TiO and H₂O. As with O-rich AGB stars, C-rich TP-AGB stars are variable in nature and so difficult to include in a stellar population model. However, they are essential contributors to the NIR light of 1–3 Gyr old stellar populations, especially at sub-solar metallicities (e.g. Ferraro et al. 1995; Lançon et al. 1999; Mouhcine & Lançon 2003; Pastorelli et al. 2019, 2020).

Similar to O-rich TP-AGB stars, Lançon & Mouhcine (2002) suggested using a NIR colour as a classification parameter for the C-rich TP-AGB star spectra in stellar population models but warned that this disregards information such as metallicity, carbon-to-oxygen ratio or pulsation properties. Although there are 51 spectra of C-rich TP-AGB stars in XSL, we selected 26 of them. The chosen spectra have the full spectrum available and are corrected for flux losses in Gonneau et al. (2020). These stars are listed in Table D.1. Loidl et al. (2001) showed that (R − J) and (R − H) are among the best effective temperature indicators for these stars, but also (V − K), (J − K), and (H − K) have been shown to correlate well with temperature (Bergeat et al. 2001). We constructed six average spectra of XSL C-rich TP-AGB stars based on these 26 XSL carbon-rich stars, sequenced and averaged based on their (R − H) colours. We prefer (R − H) broadband colour, as this results in the cleanest sequence of XSL C-rich TP-AGB stars.

We show the sequence of the average C-rich TP-AGB star spectra in Fig. 7 and the spectra inside individual bins in Fig. D.1. The number of stars in each bin varies, as the 26 spectra do not cover the (R − H) colour sequence uniformly and we aim to combine together the closest spectra in this broadband colour.

Fig. 7. Sequence of C-rich TP-AGB spectra, sorted by their (R − H) colours (on the left of each spectrum). The number of spectra the average is composed of is given in the brackets. The spectra have been smoothed to lower resolution for clarity.

3.3. Combining the interpolation methods and the average spectra of evolved giants

Figure 1 shows an example of how the global, local, and the three sequences of evolved giant star spectra are used to generate the representative spectra in different regions in an HR diagram. The cool dwarf stars are generated by the local interpolator below 4000 K. Between 4000 K and 4500 K, the resulting spectrum is the linear combination of the spectra produced by local interpolation and global interpolation, weighted by

$$\begin{array}{c}\hfill q=\frac{log(T_{\mathrm{eff}})-log(T_{\mathrm{lower}})}{log(T_{\mathrm{higher}})-log(T_{\mathrm{lower}})},\end{array}$$(5)

where T_lower = 4000 K and T_higher = 4500 K. Spectra of stars with effective temperatures between 4500 K and 7000 K are generated by the global interpolator, and star hotter than 8000 K by the local interpolator. The transition from the warm (global) to the hot (local) regime is from 7000 K to 8000 K using the weights in Eq. (5) of T_lower = 7000 K and T_higher = 8000 K.

We used isochrone keywords to determine where the isochrone track enters the relevant evolutionary stage where the static, O-rich TP-AGB, or C-rich TP-AGB spectra are used. On Padova00 isochrones, we modelled the bottom of the RGB (‘RGBb’) until the first thermal pulse (‘1TP’) with the static sequence; and from the first thermal pulse and beyond with the O-rich TP-AGB sequence. On the PARSEC/COLIBRI isochrones, we modelled stages 3 (RGB) to 7 (early-AGB, including) using the static sequence, and stage 8 (TP-AGB) with the O-rich TP-AGB sequence until the given carbon over oxygen ratio becomes one. Stars with C/O ≥ 1 were modelled using the C-rich TP-AGB sequence.

However, we only switched to the sequences when we reached the bluest average spectrum on the sequence. Hence, only the coolest (T_eff ⪅ 4000 K) giants are represented by a spectrum originating from the static, O-rich TP-AGB, or TP-AGB star sequences. Warmer stars were created with a global interpolator. There is no transition region when switching from global interpolation to the static sequence, or from the static to the O-rich TP-AGB variable star sequence, or from the O-rich TP-AGB sequence to the C-rich TP-AGB sequence. We linearly interpolated between the spectra on each sequence to infer a representative spectrum for a point on an isochrone with a given colour.

The choice of the colour–temperature relation is important in NIR stellar population modelling, and can change the NIR colours of SSPs considerably (see Sect. 8 for a discussion). For the O-rich TP-AGB sequence, we used the empirical surface-gravity-dependant (I − K) colour–temperature relation of Worthey & Lee (2011). We used the colour–temperature relation of Bergeat et al. (2001) to assign a (J − K) colour for the C-rich TP-AGB sequence stellar parameters. We note that the broadband colour we use to construct the C-rich TP-AGB sequence differs from the broadband colour we use here, because Bergeat et al. (2001) does not provide a (R − H)–temperature relation. We prefer their relation, because it is based on the measurements of angular diameters of 52 stars available from lunar occultations and interferometry, the largest set to date.

Old solar-metallicity and metal-rich populations need a template spectrum at the tip of the RGB, which is redder than the reddest spectrum on the static sequence. However, as the tip of the RGB dominates the NIR light of these populations, we cannot switch to the redder TP-AGB spectra, as this would introduce strong TP-AGB features into the population models. This issue is discussed further in Sect. 8.

3.4. Bolometric corrections

We employed the V-band bolometric corrections (BC_V henceforth) given in Worthey & Lee (2011) for all stars except the C-rich TP-AGB stars. Worthey & Lee (2011) reviews literature bolometric corrections and uses a combination of sources: VandenBerg & Clem (2003) for the middle of the temperature range, supplemented by the Vacca et al. (1996) formula for 4.40 < log(T_eff K⁻¹) < 4.75 for the hottest dwarfs and supergiants; Bessell et al. (1998) for giants; and Leggett et al. (2001) for cool dwarfs. For giants with T_eff < 4000 K, we switched from the V to I band using the (V − I) colour provided by Worthey & Lee (2011), because these stars can have little to no flux in the V band. We used the Kerschbaum et al. (2010)K-band bolometric correction for carbon-rich giants. Within the bolometric corrections, we adopted BC_V, ⊙ = −0.09, BC_I, ⊙ = 0.61, BC_K, ⊙ = 1.42 and a bolometric magnitude of 4.72 for the Sun (Torres 2010).

4. General behaviour of the models

In this section we focus on predictions of colours and absorption-line indices of our SSP models. We compare them with the E-MILES models (Vazdekis et al. 2010, 2015; Röck et al. 2016), the Maraston et al. (2009, M09 hereafter) models and the Conroy et al. (2018, C18 hereafter) models. Example spectra of these models are shown in Fig. 8. Further examples of XSL SSP models are shown in Appendix G. Here, we use the XSL SSP models calculated using the Salpeter IMF.

Fig. 8. XSL (‘PC’: PARSEC/COLIBRI; ‘P00’: Padova00), E-MILES P00, Maraston et al. (2009, M09), and Conroy et al. (2018, C18) SSP model spectra of 1 Gyr (panel a) and 10 Gyr (panel b) solar-metallicity stellar populations. The spectra are smoothed to R = 500. The M09 spectra are displayed at their original resolution (R ≈ 500). All spectra are normalised to a common I-band flux. The residual of the XSL P00 model from the E-MILES P00 model is shown in grey in each panel.

The comparison with the E-MILES is relevant because these models are widely used in the study of intermediate-age and old stellar populations (Neumann et al. 2021; Rodriguez Beltran et al. 2021; Barbosa et al. 2021; Lonoce et al. 2021, to name a few recent works). Here, we used the E-MILES models calculated using the Salpeter IMF and both the Padova00 and BaSTI (Pietrinferni et al. 2004, 2006; Cordier et al. 2007; Percival et al. 2009) isochrones. The E-MILES models cover an extensive 1680–50 000 Å wavelength range. Unlike the XSL models, they do not consist of stars observed simultaneously at all wavelengths. Instead, E-MILES models are a combination of separately generated UV, optical and NIR population models, merged at overlapping wavelengths. They make use of the Indo-US (Valdes et al. 2004), MILES (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011), Calcium II Triplet (CaT) library (Cenarro et al. 2001a,b, 2002), NGSL (Gregg et al. 2006; Koleva & Vazdekis 2012), and IRTF (Rayner et al. 2009) stellar libraries at different wavelengths. Because of this mixture, the resolution of E-MILES models varies from a constant FWHM = 2.5 Å in the NUV and optical to a constant σ = 60 km s⁻¹ in the NIR wavelengths.

The importance of TP-AGB stars in SSP models was emphasised by Maraston (2005), Maraston et al. (2006, 2009). This is why we included M09 models in some comparisons here. These solar-metallicity models extend from the UV to NIR (1150–25 000 Å) and have low resolution (R ≈ 500). The M09 models make use of the Pickles (1998) library of empirical stellar spectra. The M09 models were calculated for the isochrone sets of Cassisi & Salaris (1997), Cassisi et al. (1997, 2000). M09 used the ‘fuel consumption theorem’ with the average spectra of TP-AGB stars of Lançon & Mouhcine (2002) to include the TP-AGB stars into the SSP models. They calibrated the flux contribution of this phase against optical and NIR photometry of globular clusters in the Magellanic Clouds. Due to this particular treatment of TP-AGB stars, the NIR flux of stellar populations of ages between 0.5 and 1.5 Gyr is enhanced. This can be seen in Fig. 8a – M09 models have clearly stronger carbon star features than other SSP models.

The C18 models are similar to the E-MILES models, as they also use the MILES and the E-IRTF spectral libraries to synthesise the optical and NIR part of the SSP. However, the C18 models are based on isochrones of the MIST project (Dotter 2016; Choi et al. 2016) and use different interpolation methods than the E-MILES models. C18 models have a constant σ = 100 km s⁻¹ resolution.

The hottest turnoff star determines the shape of the optical population model. The stars on the tip of the RGB dominate the NIR light of the 10 Gyr population, but the TP-AGB stars dominate the NIR light of the 1 Gyr population. This leads to larger differences between different SSP models for 1 Gyr populations than for the 10 Gyr populations that are seen in Fig. 8 – TP-AGB stars are more difficult to incorporate into SSP models than RGB stars. This will be further discussed in Sect. 8.

4.1. Colours measured from our models

Figure 9 shows the behaviour of the optical/NIR colours measured from our Padova00- and PARSEC/COLIBRI-based SSP models and from the other models discussed above. We show the colour behaviour as a function of age (left panels) and metallicity (right panels). Ages span from 1 Gyr to 16 Gyr, and metallicities span from [Fe/H] = −2.2 dex (XSL PARSEC/COLIBRI) or [Fe/H] = −1.7 dex (other) to [Fe/H] = +0.2 dex. We note that E-MILES models are not safe to use in the NIR below [Fe/H] = −0.4 dex, but are included for illustrative purposes.

Fig. 9. Comparison of the behaviour of colour as a function of age and metallicity for the E-MILES, M09, C18, and XSL SSP models (see legend). Left panels: colour as a function of age. Right panels: colour as a function of metallicity. Top row: (B − V). Top-middle row: (V − I). Bottom-middle row: (I − J). Bottom row: (J − K). Ages span from 1 Gyr to 16 Gyr, and metallicities span from [Fe/H] = −0.4 dex to [Fe/H] = +0.2 dex (left panels); [Fe/H] = −2.2 dex (XSL PARSEC/COLIBRI) or [Fe/H] = −1.7 dex (other) to [Fe/H] = +0.2 dex (right panels). We note that E-MILES models are not safe to use in the NIR below [Fe/H] = −0.4 dex but are included for illustrative purposes. Solar-metallicity E-MILES models are shown in heavier line strengths than the sub- and super-solar models of the left panels. XSL SSP sub- and super-solar models are represented by shaded areas, centred on the solar metallicity. We note the different colour-scale values between the same colour panels.

Age–colour relations: All models follow the same trends as our SSP models and become redder in (B − V), (V − I) and (I − J) with increasing age. The NIR (J − K)–age relation is flat for old ages. There are some notable differences between models: at super-solar metallicities, E-MILES BaSTI models have (B − V) and (V − I) colours similar to XSL, C18 and M09 solar models. NIR colours of E-MILES super-solar models are redder (Δ(I − J)≈0.1) than XSL models. Even the (I − J) colours of E-MILES solar-metallicity models are ∼0.05 redder than other models. Furthermore, the M09 and C18 solar-metallicity models are somewhat similar to the [Fe/H] = −0.4 models of XSL and E-MILES in (J − K). (I − J) and (J − K) have model-dependent behaviour in the TP-AGB regime (ages < 3 Gyr).

Metallicity–colour relations: All models follow the same trends as our SSP models and become redder in all colours with increasing metallicity. While the (B − V)–metallicity relation is almost identical for XSL and E-MILES models, differences arise towards the NIR. XSL Padova00 models have bluer (V − I), (I − J) and (J − K) colours than XSL PARSEC/COLIBRI models. Considering the range of metallicities where E-MILES models are safe to use ([Fe/H] ∈ [ − 0.4, 0.0, −0.2]), (I − J) colour stands out having a steeper metallicity–colour relation than XSL models.

It is hard to pinpoint a single reason for these model discrepancies, specially in the NIR. Differences in used empirical libraries is one of them, but E-MILES and C18 models do not agree as well. Issues arising from E-MILES or C18 SSP model merging or XSL DR3 merging of stellar spectra are another possible source of disagreements between models. Moreover, we include cool giants into SSP models differently than other groups, with the use of the static and variable sequences. The NIR colour differences between XSL Padova00 and PARSEC/COLIBRI reflect the usage of isochrones with different levels of sophistication for the description of the TP-AGB phase. Sub-solar XSL PARSEC/COLIBRI SSP models, which have more thorough description of the TP-AGB phase than the XSL Padova00 models, show bluer NIR colours (Δ(I − J)≈0.06 and Δ(I − J)≈0.07 at [Fe/H] = −1.0); differences are small for solar-metallicity models, but noticeable in Fig. 8.

We concentrate on the comparison with E-MILES, as those models are widely applied and their behaviour studied. Röck et al. (2016) has presented a thorough analysis of E-MILES optical and NIR colours.

4.2. Optical absorption-line indices measured from our models

We compared the widely used optical absorption-line indices measured from the XSL and E-MILES SSPs, using diagnostic plots such as Hβ versus Mgb, Ca4455, Fe5015, NaD (Trager et al. 1998), CaHK (Serven et al. 2005), and [MgFe] indices in Fig. 10. [MgFe] is defined by Thomas et al. (2003) as

$$\begin{array}{c}\hfill [\mathrm{MgFe}]\equiv \sqrt{\mathrm{Mg}b\times (0.72\times \mathrm{Fe}5270+0.28\times \mathrm{Fe}5335)}.\end{array}$$(6)

Fig. 10. Comparison of the behaviour of model Mgb, CaHK, Ca4455, [MgFe], Fe5015, and NaD absorption-line index strengths as a function of the model Hβ index strength. The shaded area represents XSL Padova00 models with varying spectral resolution (σ) from σ = 13 km s−1 (the native XSL resolution) to σ = 60 km s−1 (the minimum E-MILES resolution). Black lines represent E-MILES Padova00 model predictions, with dotted, dashed, and solid lines representing [Fe/H] = +0.2, 0.0, and −0.4 dex, respectively, measured at the original E-MILES resolution. The solid blue line represents predictions from the C18 solar-metallicity models.

In Fig. 10, we show measurements from the XSL Padova00 models and E-MILES P00 models older than 1 Gyr and with metallicities [Fe/H] ∈ [ − 0.4, 0.0, +0.2] dex. Similar grids for XSL PARSEC/COLIBRI and E-MILES BaSTI models are shown in Fig. E.1, but with an extension towards the lowest metallicities of XSL SSP models. Furthermore, we added absorption-line indices measured from solar C18 models with ages between 1 and 13 Gyr, measured at its original σ = 100 km s⁻¹ resolution.

The optical absorption-line index trends of different models are similar. The comparison of XSL and E-MILES Padova00 models shows that some differences arise from the different stellar spectra and interpolation methods. However, the differences between C18 and E-MILES models illustrate how well models using the same stellar library (MILES/IRTF) but different stellar population modelling techniques compare.

There are a few notable differences between the grids. On one hand, the Ca II absorption-line index CaHK shows different behaviours at older ages. There is a saturation seen for the oldest population models, but this saturation happens at different metallicities for the varying models. XSL spectra have stronger index values than E-MILES models, while the C18 models are closer to XSL in this index. This is a prominent spectral feature in SSPs, coming from F, G, and K stars. On the other hand, there is roughly less than a 0.15 Å disagreement between the models for Ca I line index Ca4455. There is also a crossing of Mg (Mgb and [Mg/Fe]) index values of young SSP models with different metallicities. Furthermore, XSL models show a larger spread in NaD index values at a given Hβ strength than the other models; however, this line lies in the dichroic contamination region in the XSL spectra, and should be used with (extreme) caution.

5. Colours of Coma cluster galaxies

To show the potential of our new XSL models, we compare model colour predictions with photometry of galaxies in the Coma cluster on the colour–colour planes in Fig. 11. We use these galaxies since many galaxies in a rich cluster show colours consistent with old SSP models (e.g. Bower et al. 1992). The photometry is taken from Eisenhardt et al. (2007, Table 9). We only show the photometry of galaxies with redshift within 3σ of the average redshift of the Coma cluster galaxies from Upadhyay et al. (2021), z = 0.0224 ± 0.0033. This results in 180 galaxies, the majority of which are early-type galaxies (ETGs). We redshifted the models to z = 0.0224 and have used the response functions provided by Eisenhardt et al. (2007) for the spectrophotometry.

Fig. 11. Colour–colour diagrams of galaxies in the Coma cluster with model predictions overlaid. Grey points show the Eisenhardt et al. (2007) data for galaxies in the Coma cluster. Panel a: XSL and E-MILES Padova00 models. (b) XSL PARSEC/COLIBRI and E-MILES BaSTI models. Only [Fe/H] = +0.2, 0.0, and −0.4 dex E-MILES models are shown. C18 models with [Fe/H] = +0.2, 0.0, and −0.5 dex models are denoted by solid blue lines. XSL models extend to even lower metallicities.

The XSL SSP models reproduce optical colours of ETGs well in general. However, there are some cases where models do not match the data. For example, all models shown in Fig. 11 have redder (V − I) (∼0.1 mag) colours than the galaxies at fixed (B − R). In the solar and metal-rich regime, the XSL SSP models are most cases bluer than E-MILES and C18 in the NIR. This is very apparent in the (I − J) or (V − K) colours at fixed (B − R), where XSL SSP models are roughly ∼0.1 mag bluer. The colours containing the I- and J-filter are particularly interesting since they encompass the joining region the UV-optical MIUSCAT SSP models (based on MILES, Indo-US and CaT libraries) and IRTF-based NIR SSP models in the E-MILES models and the merging of the individual spectral arms in the XSL models.

The colour offsets, especially in (I − J), can be due to merging of the XSL DR3 VIS–NIR spectra. It can also be due to inclusion of cool giant stars using separate giant sequences and the colour–temperature relation. The NIR colours are constrained by the reddest static and variable giant templates.

Model offsets in optical colours have been discussed in detail by Ricciardelli et al. (2012), who tested the MIUSCAT models, the optical part of the E-MILES models, on nearby ETGs. None of the MIUSCAT SSP models are able to match some of the observed optical colour distribution (namely (u − g) or (r − i) colours at fixed (g − r)) of nearby ETGs, while the colours of Milky Way globular clusters are reproduced remarkably well. They suggest that the ETGs of their sample are not necessarily simple old stellar populations, and need small contributions from either young or/and metal-poor stellar populations. Furthermore, the impact of α-enhancement and the choice of IMF on galaxy colours cannot be neglected.

6. Optical/NIR absorption-line indices

To date, stellar population studies of unresolved galaxies have used mainly the optical absorption-line indices, but the NIR spectral features can provide insights into the stellar populations dominated by cool stars (Lançon et al. 1999, 2008; Mouhcine et al. 2002; Riffel et al. 2007, 2008, 2015, 2019; Mármol-Queraltó et al. 2009; Kotilainen et al. 2012; Lyubenova et al. 2012). Riffel et al. (2019) presented 47 correlations among the different absorption features in the optical and NIR for 16 star-forming galaxies (SFGs) and for 19 ETGs. They found that the models consistently agree with the observations for the optical absorption features, but not so much for the NIR indices.

Motivated by this discrepancy, we looked at some of the suggested indices from Riffel et al. (2019) and compare them with the XSL PARSEC/COLIBRI SSP model predictions and those from the E-MILES BaSTI models. Although Riffel et al. (2019) found correlations among the different absorption features in the optical and NIR, seemingly suggest an evolution from an SFG to an ETG, multiple stellar populations are likely to be an issue when attempting to compare optical and NIR indices of SFGs.

We selected six NIR indices and plot them against [MgFe]. These index–index diagrams are shown in Fig. 12. We limited the comparisons to XSL PARSEC/COLIBRI Salpeter and E-MILES BaSTI Salpeter SSP models only, as they have a more up-to-date handling of cool giant evolutionary phases. We note that Riffel et al. (2019) used E-MILES Padova00 models in their comparisons.

Fig. 12. Selected index–index comparisons from Riffel et al. (2019). Shaded areas represent XSL PARSEC/COLIBRI SSP model predictions, with red, yellow, and teal indicating [Fe/H] = +0.2, 0.0, and −0.40 dex, respectively. Shaded areas represent models with a spectral resolution of σ = 16 km s−1 (the native XSL resolution) to σ = 228 km s−1 (the resolution of the Riffel et al. 2019 spectra) centred on σ = 60 km s−1 (the E-MILES NIR resolution). Black lines represent E-MILES BaSTI model predictions, with dotted, dashed, and solid lines representing [Fe/H] = +0.26, 0.06, and −0.35 dex, respectively (roughly the same metallicities as the XSL models). We note that the age range differs for the E-MILES models as the NIR spectra of E-MILES are only reliable above 1 Gyr. Panel b: the CN1.10* (Röck 2015) index definition is used instead of the Riffel et al. (2019) definition, which is affected by residuals from telluric absorption correction. We omitted the SFG and ETG measurements of CN1.1 due to differences in index definitions. Panels b–d: XSL static sequences, O-rich TP-AGB sequences, C-rich TP-AGB sequences, and XSL supergiants (which are not included in the XSL SSP models) are shown in grey at arbitrary optical index values (as these stars lack optical features) and their median values with larger black symbols. Indices of SFGs are marked in blue and indices of ETGs in red. These values are taken from Tables 6–7 and B1–B3 of Riffel et al. (2019), respectively.

The ZrO–[MgFe] comparison in Fig. 12a shows that this line is affected by the CN features of C-rich TP-AGB stars in XSL SSP models. Metal-poor stellar population models with ages less than 3 Gyr show a steep increase in the strength of this index due to dominance of TP-AGB stars, especially the C-rich stars. Carbon stars have very high ZrO index values, around 70 Å, but we omitted the giants from Fig. 12a for clarity. The E-MILES models show a different behaviour of ZrO – SSP models with ages less than 2 Gyr showing a steep decrease in the strength of this index.

We also included the CN1.10–[MgFe] comparison to illustrate the behaviour of this important NIR index. However, in Fig. 12b we use the CN1.10 index definition of Röck (2015). We omitted the SFG and ETG measurements of CN1.10 due to differences in index definitions between Röck (2015) and Riffel et al. (2019). Riffel et al. (2019) defined the CN1.10 red continuum band at 11310–11345 Å, coinciding with the region of severe telluric absorption. We see some residuals in the telluric absorption region of XSL spectra, which also affects the NaI1.14 index. The CN1.10 (Röck 2015) index definition, the same definition we used to remove supergiants, has the red continuum placed at 11100–11170 Å, away from the telluric contamination. The CN1.10–[MgFe] comparison in Fig. 12b shows a systematic offset between the E-MILES and XSL models in CN1.10 index values. The smaller XSL predictions are a direct consequence of separating C-rich TP-AGB stars and removing supergiant stars, as described in Sect. 3.2. The supergiants are not included in the XSL SSP models, but are shown in Fig. 12b-c for illustrative purposes. The XSL SSP models have mostly shallower CN1.10 features. However, stellar population models with ages less than 3 Gyr show a steep increase in the strength of this index due to C-rich stars.

As seen in Fig. 12c–d, SFGs show similar, if not stronger, CO1.5a and CO2.2 index features compared to ETGs. However, none of the SSP models reproduce these strong CO features. Both carbon and oxygen are abundant elements in cool giants, and these molecules are formed and can be observed in both M and C stars. But CO2.2 lines and CO1.5 lines originate from different regions within the extended atmospheres of cool stars (Nowotny 2005). Furthermore, stellar population H-band CO lines are blends. The CO1.5a line is a blend of CO and Mg I. The CO2.2 line is an almost pure CO feature (Riffel et al. 2019). The static sequence spectra have strong CO2.2 indices, influencing older models to have strong CO2.2 values. TP-AGB stars show a variety of CO2.2 line strengths – the TP-AGB phase does not substantially influence the CO index. This was also concluded by Röck (2015), and the same can be seen for the CO1.5a index. CO is expected to be enhanced in younger (< 50 Myr) stellar populations (e.g. Lançon et al. 2008; Riffel et al. 2007, 2015) due to the presence of supergiant stars, which are not included in the XSL SSP models.

If an SFG hosts even a small population of supergiants, the NIR CO and CN indices will be affected, but the optical indices might not be affected by this younger population component. This is clearly seen from Fig. 12b–d, where the CO and CN index strengths of the XSL supergiants are much stronger than of the other cool giants.

It is now possible to perform in-depth studies of spectral features in the NIR. The XSL SSP models are useful tools due to the moderate-to-high-resolution spectra of the XSL. Furthermore, the models include a large number of spectra of cool giant stars. On the one hand, XSL SSP models improve the model range of some lines, such as the MgI1.7 line in Fig. 12e. On the other hand, XSL models expand the range of predicted values of the NaI2.2 index (Fig. 12f), but towards lower index values, contrary to the strong index values of SFGs and ETGs. Individual elemental abundance variations, velocity dispersion broadening, wavelength shifts, residuals from telluric absorption correction, S/N, flux calibration, IMF, inclusion of cool giant stars, and the presence of multiple stellar populations can all influence NIR spectral line indices. Indeed, Röck et al. (2017) and La Barbera et al. (2017) suggested that for ETGs the large values obtained for the NaI2.2 index are due to a combination of a bottom-heavy IMF and enhanced sodium abundances. Further research is needed for the majority of the NIR spectral features, using purposefully defined NIR indices, such as those of Eftekhari et al. (2021). A full analysis of the colours and indices of the galaxies of Riffel et al. (2019) over the X-shooter range of wavelengths requires models with non-trivial SFHs and configurations (for instance a prescription for the spatial distribution of dust relative to young and old stars), and lies outside the scope of this paper.

7. Stellar mass-to-light ratios

The stellar mass-to-light (M_*/L hereafter) ratio is an important characteristic of a stellar population. Many of the population properties (e.g. morphology or SFH) are correlated with the stellar mass (e.g. Bernardi et al. 2020; Telford et al. 2020; Ge et al. 2021; de Graaff et al. 2021; D’Eugenio et al. 2021, to name some recent works). The stellar mass of a population is not a directly observable quantity but its luminosity is. One way of estimating population mass is through the synthetic M_*/L ratio: in such a case, the population light is converted into a mass using a stellar M_*/L ratio derived from stellar population models (see the reviews by Conroy 2013; Courteau et al. 2014).

Existing stars in a stellar population contribute to the mass and luminosity of that population. But stars progressively die and turn into stellar remnants (white dwarfs, neutron stars, and black holes) as the stellar population ages. Those remnants contribute to the mass but not to the luminosity. The total mass of an SSP with a certain age and metallicity is the sum of stellar and remnant masses, weighted by the IMF. The weight of the IMF is determined by the initial mass of the star, but the mass of the star/remnant at that time is what contributes to the mass budget.

In Fig. 13 we present the synthetic M_*/L ratios derived from the XSL PARSEC/COLIBRI models with Salpeter IMFs in the V, I, and K photometric bands. The luminosity is given in units of solar luminosity in the respective photometric band. The solar magnitudes used are: (V, I, K) = (4.81,4.11,3.30) mag, measured from the Solar spectrum of Colina et al. (1996). We used the relation in the mass range 0.09 < m/M_⊙ < 120. The PARSEC/COLIBRI models describe the mass loss of stars and provide both initial and actual stellar masses for existing stars. We used the metallicity-dependent initial–remnant mass relation descriptions provided in Fryer et al. (2012) for massive stars (9–120 M_⊙) for non-solar metallicities and the Sukhbold et al. (2016) relation for solar metallicity. For low- and intermediate-mass stars (0.87 < M_*, init < 8.2 M_⊙), we used the PARSEC-based white dwarf initial–final mass relation of Cummings et al. (2018), extrapolating the relation to 8.2–9 M_⊙. We assume that the mass lost in the form of ejected gas is blown out of the stellar population and does not contribute to the mass budget.

Fig. 13. Evolution of synthetic (log) stellar mass-to-light ratios for the V band (left), I band (middle), and K band (right), measured from the XSL PARSEC/COLIBRI Salpeter models.

The dominant driver of SSP luminosity is its age, as the most-massive stars have short lifetimes but are orders of magnitude more luminous than the less massive stars. The luminosity of an SSP changes rapidly with time. The mass of an SSP is dominated (for the Salpeter IMF considered here) by the least-massive stars. These stars live a long time, and thus the mass of an SSP changes little after the first few gigayears. As seen from Fig. 13, the M_*/L ratio changes rapidly until about 2 Gyr, with the most massive and luminous stars dying off. The effect of metallicity on the M_*/L ratio is weaker. For stellar populations older than a few gigayears, the higher the stellar population’s metallicity, the higher the M_*/L ratios for optical passbands, but the (slightly) lower the NIR M_*/L ratio.

The differences between the M_*/L ratios in the V, I, and K photometric bands are expected, as the hottest turnoff star determines the V-band luminosity; the stars at the tip of the RGB determine the K-band luminosity of old stellar populations, and the TP-AGB stars determine the K-band luminosity of 50 Myr to 2 Gyr populations. Furthermore, the influence of TP-AGB stars on the M_*/L ratio peaks for populations with ages between 0.4 to 1.58 Gyr and metallicities between [Fe/H] = −0.6 and 0. These stars emit mostly in the NIR, increasing the NIR luminosity and lowering the NIR M_*/L ratio. This can be clearly seen in Fig. 13c. The M_*/L ratio is dependent on the stellar evolutionary phases accounted for in the modelling. Without the TP-AGB stars, the M_*/L would increase monotonically with age.

The M_*/L ratio is strongly dependent on the IMF. We provide discussion of the M_*/L ratios from XSL PARSEC/COLIBRI models calculated with other IMFs in an upcoming paper (Verro et al. in prep.). Furthermore, there are differences between M_*/L ratios determined from different models. We discuss this briefly in Appendix F.

8. On the separation of static and variable giants

Separating static cool giant (from RGB to early-AGB) stars from the variable TP-AGB stars with the use of the static and variable sequences in XSL SSP models is an important step towards understanding the source of NIR flux in stellar populations. These stars lie very close to each other on the HR diagram, but their spectral shapes can be very different, as discussed in Sect. 3.2. There is an ongoing debate as to their impact on the integrated spectra of even SSPs. Clear C-rich TP-AGB signatures have been detected in some of the J- and H-band spectra of globular clusters in the LMC (Lyubenova et al. 2012). These globular clusters are of intermediate age (1–2 Gyr) and have metallicities around [Fe/H] = −0.4 dex. However, Zibetti et al. (2013) explored a set of post-starburst galaxies, with luminosity-weighted ages between 0.8 and 1.6 Gyr and metallicities between [Fe/H] = −0.68 and +0.3 dex and found no strong spectral signatures of these stars. This discrepancy has been explained by Girardi et al. (2013) by AGB boosting effect, which is linked to the physics of stellar interiors – stellar populations in a narrow 1.57 and 1.66 Gyr age range at MC metallicities have TP-AGB contribution to the integrated luminosity of the stellar population increase by a factor of ∼2. This was recently confirmed by Pastorelli et al. (2020); their modelling showed a 80% peak in K-band flux coming from (mainly C-rich) TP-AGB stars.

8.1. RGB and TP-AGB light fractions in the NIR

We show the contribution of RGB stars and TP-AGB stars to the total K-band luminosity of the XSL models in Fig. 14. In the XSL PARSEC/COLIBRI SSP models, the RGB contribution changes from low in young populations to high in old populations, with a strong transition around 2 Gyr. A high contribution of TP-AGB stars to the K-band flux extends roughly from 0.5 to 1.6 Gyr, contributing 40% or more of the flux in the K band at these ages, peaking around 0.8 Gyr and [Fe/H] = −0.2 dex, contributing 55%–60% of the K-band flux in this population. Mainly C-rich TP-AGB stars contribute to this peak. C-star H-band signatures can be recognised in Fig. 8a for the 1 Gyr solar-metallicity XSL PARSEC/COLIBRI SSP models. However, we do not see the AGB boosting peak in our 1.58 Gyr and [Fe/H] = −0.4 dex SSP models. For younger and older ages, the predicted TP-AGB contribution is almost entirely due to O-rich stars and increases with metallicity, together with the O-rich TP-AGB lifetimes. There is another peak in the TP-AGB contribution in very young and very metal-poor SSP models, where the flux contribution from other stars is lower.

Fig. 14. Contribution of RGB stars and TP-AGB stars to the total K-band luminosity of the XSL models. Panels a and b: Contribution of RGB and TP-AGB stars, respectively, in the PARSEC/COLIBRI models. Panels c and d: contribution of RGB and TP-AGB stars, respectively, in the Padova00 models.

These behaviours are expected. The TP-AGB phase for low- and intermediate-mass stars (M = 2–7 M_⊙) culminates in stellar populations of ages between 0.5 and 2 Gyr. These stars emit mainly in the NIR spectral range, given their low temperatures. The SSP models calculated in Pastorelli et al. (2020) predict a TP-AGB contribution peak at around 1 Gyr (roughly between 0.3 and 2 Gyr) that does not exceed 55% in the K-band luminosity. In comparison, the peak is as high as 80% in the K band at [Fe/H] = −0.3 dex for the M05 models.

On the other hand, the XSL Padova00 models have a completely different TP-AGB fraction behaviour with SSP parameters: the younger and more metal-poor models have higher TP-AGB fraction. For example, a 0.2 Gyr, [Fe/H] = −1.7 dex model has 80% of its K-band flux coming from TP-AGB stars. This is why we discourage the usage of XSL Padova00 models outside of the narrow safe zone suggested in Sect. 9.

8.2. Colour–temperature relations

We used colour–temperature relations of Worthey & Lee (2011) and Bergeat et al. (2001), shown in Fig. 15, to assign an average spectrum of a O-rich static, variable or C-rich star to a point on an isochrone when generating an SSP model. This allows us to bypass stellar parameter estimation for the complex XSL stars that make up these average spectra. This assignment comes with some caveats.

Fig. 15. Colour–temperature relations for C-rich and O-rich stars. The grey horizontal lines (upper panel) represent the colours of the bluest and reddest average spectra on the C-rich TP-AGB sequence. The horizontal red and dark red lines represent the bluest and reddest spectra from the O-rich static and O-rich TP-AGB sequences, respectively.

We made a version of the models using the colour–temperature relation of Lejeune et al. (1997) to compare with the Worthey & Lee (2011) relation used in our default models. Figure 16 shows that the Lejeune et al. (1997) relation undesirably enhances the NIR fluxes and introduces stronger O-rich TP-AGB features (such as the H-band H⁻/H₂O feature) to the models.

Fig. 16. 1 Gyr old, super-solar-metallicity XSL SSP generated using the Lejeune et al. (1997) (red) and Worthey & Lee (2011) (black) colour–temperature relations.

We also note here that the static sequence does not have spectra that appropriately represent RGB or early-AGB stars with colour-inferred temperatures less than 3100 K. This affects the older metal-rich models the most, where the tip of the RGB dominates the NIR light. This is the reason why the super-solar-metallicity XSL SSP models cannot reach NIR colours as red as E-MILES models in the colour comparisons in Sects. 4 and 5. The E-MILES models might achieve these colours by effectively mixing the redder spectra of TP-AGB with the spectra of RGB stars (within the local interpolator scheme) at the tip of the RGB in these populations, as those models do not distinguish between RGB and post-RGB (early-AGB and TP-AGB) stars. Our empirical separation into static and O-rich TP-AGB stars based in Fig. 4, might exclude from our static sequence a few spectra of very cool stars that would in fact be acceptable representations of the coolest RGB stars. However, it remains difficult to match these cool spectra with synthetic ones (see Lançon et al. 2019, 2021) and hence to separate effects of temperature, metallicity, circumstellar extinction and variability. A star-by-star study that would incorporate variability and mid-infrared information where available, may help improve this separation in future versions.

Furthermore, the XSL does not have spectra that appropriately represent O-rich TP-AGB stars with colour-inferred temperatures less than 2700 K and C-rich TP-AGB stars with colour-inferred temperatures less than 1700 K. The Worthey & Lee (2011) colour-temperature relation for O-rich giants does not go to lower temperatures. Also, none of the colour-temperature relations discussed in Worthey & Lee (2011) go to lower temperatures. By contrast, the PARSEC/COLIBRI models include extremely cool TP-AGB stars (e.g the 1 Gyr model in Fig. 17). The reddest average spectrum of the O-rich or C-rich TP-AGB sequence will represent these stars in our models.

Fig. 17. Selected PARSEC/COLIBRI solar-metallicity isochrones illustrating the coverage of isochrones by the XSL DR3. The locations of the tip of the RGB for the 1 Gyr and 10 Gyr isochrones are indicated with horizontal black lines. We have employed the colour-temperature relations used in this study to add the giants to this figure to illustrate how well the cool RGB and AGB stars are represented (at arbitrary log g values). Because luminous cool stars are particularly important contributors to the red and NIR light of galaxies, XSL was designed to contain a large number of such objects. We are limited by the lack of hot stars and as such cannot create XSL SSP models younger than 50 Myr.

8.3. Metallicity effects

The metallicities of the majority of the stars from which the static, O-rich, and C-rich TP-AGB sequences were constructed are unknown or not accurately known. These stars come from a variety of environments – the solar neighbourhood, star clusters, the Galactic bulge, and the Magellanic Clouds. Furthermore, the range of ages in those stars is likely to correspond to a range of metallicities, as determined by the chemical evolution of these environments. Hence, we combined together spectra with various metallicities. We note that the SEDs of C-rich TP-AGBs are not as sensitive to metallicity as to the effective temperature, or the C/O ratio (Lançon & Mouhcine 2002), so we do not consider these stars in this discussion.

The combination of O-rich TP-AGB stars with different metallicities will have three separate consequences on the resulting SSP modelling. Firstly, the relation between effective temperature and spectrophotometric properties changes due to molecular opacities. Lower metallicity cool giants have bluer spectra with weaker molecular bands (Hauschildt et al. 1999; Lançon & Mouhcine 2002). This would mean one set of static and O-rich TP-AGB averages would not be enough and we would need to create static and O-rich TP-AGB sequences for different metallicity bins separately. Considering that each average on the static and O-rich sequence consists of a handful of spectra, further division would be impossible. Furthermore, neither the Bergeat et al. (2001) or Worthey & Lee (2011) relations we use in the modelling are metallicity dependent.

Secondly, the AGB evolutionary tracks shift to lower effective temperatures. Again, signatures of cooler giant stars will be less pronounced at lower metallicities. We do take this effect into account, as the lower metallicity isochrones we use here shift into cooler temperatures and we select giants with lower colour-temperature from the static and O-rich and TP-AGB sequences.

The third effect is the metallicity dependence of the mass loss. It affects the AGB lifetimes and the efficiency of the production of C-rich TP-AGB stars. This is addressed in the stellar evolution calculations, which we use as an input (e.g. Pastorelli et al. 2020, and other papers from the PARSEC/COLIBRI group).

9. Applicability of the XSL stellar population models

The coverage of the spectral library, the capabilities of interpolator(s), and the selected isochrones determine the stellar population models that we can create. Figure 18 shows the parameter coverage of the XSL stars for dwarfs and giants (separated at log g = 3). At solar metallicity, all types of stars are well represented, and the coverage is good even for the lower metallicities. Reliable models can be computed down to [Fe/H]= − 2.2 dex, for very old populations, due to lack of hot metal-poor (6500 < T_eff < 10 000 K) stars, and down to [Fe/H]= − 1.6 dex for intermediate-age stellar population models. The coverage of metal-rich dwarf and giant stars allows us to safely compute stellar population models up to [Fe/H]= + 0.20 dex.

Fig. 18. Fundamental parameter coverage of the XSL stars, for dwarfs (left panel) and giants (right panel). Dwarfs and giants are separated at log g = 3. The stars that are used to construct static, variable (O-rich TP-AGB), and carbon (C-rich TP-AGB) star sequences are not included, as we ignore the information about their metallicity in the models (Sect. 3.2).

An HR diagram, such as those shown in Fig. 17 with 50, 80, 100, 150, 200 Myr, 1 Gyr, and 10 Gyr solar-metallicity PARSEC/COLIBRI isochrones, gives another perspective. We are limited by the lack of hot stars in XSL for very young ages even at solar metallicity. Furthermore, we have removed supergiants. As the hottest turnoff star determines the shape of the optical population model and the supergiants dominate the NIR light of young populations, we can create stellar population models of 50 Myr and older. However, the lack of blue loop stars limits the models to ages older than roughly 80 Myr at solar metallicity or older than roughly 100 Myr at sub-solar metallicities.

Because luminous cool stars are particularly important contributors to the red and NIR light of galaxies, XSL was designed to contain a large number of such objects. The PARSEC/COLIBRI models in Fig. 17 include TP-AGB stars cooler than any XSL TP-AGB star, given the colour-temperature relations of Worthey & Lee (2011). We believe the XSL sample has allowed for significant progress in inclusion of TP-AGB stars in SSP models. That is why we do not base our judgement of the age and metallicity limits of the models on the XSL coverage of this extreme region of the HR diagram. However, the limits of the XSL Padova00 models are justified based on the discussion in Sect. 8.1 regarding the handling of TP-AGB stars. Figure 19 shows the age and metallicity limits of the XSL SSP models, based on the age and metallicity coverage of the isochrones and the temperatures and metallicity of the XSL stars included in the model.

Fig. 19. Safe zones of XSL SSP models. ‘P/C’ indicates PARSEC/COLIBRI-based models, with white and grey indicating safe and unsafe zones, respectively. ‘P00’ indicates Padova00-based models, with the grey outline delimiting the safe and unsafe zones.

Stellar population modelling smears out the individual issues of the constituent stellar spectra to some extent. Nevertheless, we warn that certain areas in the SSP model spectrum may be influenced by the dichroic contamination or residuals from telluric absorption correction present in many of the XSL DR3 spectra (Verro et al. 2022). Table 1 defines the problematic spectral regions.

XSL SSP models are based on an empirical stellar library, which, like any empirical stellar library, has limited coverage of the HR diagram. Here we give rough limits as to what extent we believe the XSL SSP models to be reliable to use. We note that Coelho et al. (2020) described extensively the effect of the coverage of the stellar library on the SSP model predictions. They found that predicted colours are more affected by the coverage effect than the choice of a synthetic versus empirical library. Derived galaxy ages can be underestimated when stellar population synthesis models with limited parameter coverage are used. On the other hand, metallicities are robust against limited HR diagram coverage but are underestimated when using synthetic libraries.

10. Conclusions

We present the XSL SSP models, which are based on 639 stellar spectra from XSL DR3. These SSP models have various improvements compared to other available models. XSL SSP models cover a wide wavelength range, from the NUV (350 nm) to the NIR (2480 nm); they have moderate-to-high resolution throughout the wavelength range, with original σ = 13/11/16 km s⁻¹ in the UVB/VIS/NIR arms of X-shooter; and they are constructed from stars for which the spectra have been observed simultaneously at all wavelengths and extend over a metallicity range of −2.2 < [Fe/H] < +0.2 dex and an age range of 0.05 < t_SSP < 16 Gyr.

To construct these models, we have used recent PARSEC/COLIBRI stellar evolutionary tracks, which include the TP-AGB phases that control the evolution of NIR colours. Particular care was taken to include the RGB, early-AGB, and TP-AGB stars in the stellar population models, with the use of average spectra of static giants, variable O-rich TP-AGB stars, and C-rich TP-AGB stars. Instead of relying on the stellar parameter estimation of these stars, we used the established colour–temperature(–metallicity) relations of Worthey & Lee (2011) (for O-rich cool giants) and Bergeat et al. (2001) (for C-rich cool giants). We also provide XSL SSP models constructed with older Padova00 stellar evolutionary tracks, but we discourage the use of these models outside the narrow safe zone defined in Sect. 9.

We have gone through an extensive characterisation of the stellar population models. We have compared colours and absorption-line indices with existing stellar population models (E-MILES, C18, and M09). We have also compared our model predictions with the colours of Eisenhardt et al. (2007) ETGs from the Coma cluster and spectral features of SFGs and ETGs from Riffel et al. (2019) and find encouraging agreement with the observations. The XSL SSP models can reproduce the optical colours of ETGs in the Coma cluster, which is comparable to the success of the E-MILES and C18 models. The differences between the models are largest at the NIR super-solar metallicities of old populations, which may be due to the inclusion of cool giant stars using separate giant sequences and the colour–temperature relation. The NIR colours are constrained by the reddest static and variable giant template. Offsets in (I − J) might also come from the inaccurate merging of the VIS and NIR arms of the spectra of cool giants. While the behaviour of optical absorption-line indices is similar in E-MILES and C18, there are discrepancies between models for NIR indices. The XSL models improve the range of predicted values for many NIR indices, such as MgI1.7. Careful separation of XSL RGB, early-AGB, and TP-AGB stars and including them in the XSL SSP models will allow us to analyse NIR indices more systematically in the future.

The extended wavelength coverage and high resolution of the new XSL-based stellar population models will help us to bridge optical and NIR studies of intermediate-age and old stellar populations and clarify the role of evolved cool stars in stellar population synthesis.

Acknowledgments

A.A. acknowledges funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No. 834148). R.F.P. acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 721463 to the SUNDIAL ITN network. A.V. and J.F.-B. acknowledge support through the RAVET project by the grant PID2019-107427GB-C32 from the Spanish Ministry of Science, Innovation and Universities (MCIU), and through the IAC project TRACES which is partially supported through the state budget and the regional budget of the Consejería de Economía, Industria, Comercio y Conocimiento of the Canary Islands Autonomous Community. P.C. acknowledges support from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under grant 310041/2018-0 and from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) process number 2018/05392-8. P.S.B. acknowledges the financial support from the Spanish National Plan for Scientific and Technical Research and Innovation, through the grant PID2019-107427GB-C31. L.M. thanks FAPESP (grant 2018/26381-4) and CNPQ (grant 306359/2018-9) for partial funding of this research.

References

Appendix A: The weighted median residuals around four spectral line indices

Fig. A.1. Weighted median residual, RS, between the original spectrum and the interpolated spectrum for four spectral line indices (CaHK, Hbeta, CaT, and CO1.6) as a function of position on the HR diagram. We normalised the spectra over the wavelengths of interest before calculating the line-level median residuals. The colour bar is logarithmic. Histograms show the distributions of RS calculated within these spectral ranges at the full XSL resolution. For ease of visualisation, spectra with RS > 0.1 are placed into the RS = 0.1 bin in the histograms.

Appendix B: Static giants

Fig. B.1. XSL spectra of O-rich, cool static giant stars from which the static sequence is constructed.

Appendix C: O-rich TP-AGB star average bins

Fig. C.1. XSL spectra of O-rich TP-AGB stars from which the variable O-rich TP-AGB sequence is constructed.

Appendix D: Carbon star average bins

Fig. D.1. Spectra of XSL C-rich TP-AGB stars from which the C-rich TP-AGB sequence was constructed.

Appendix E: XSL PARSEC/COLIBRI and E-MILES BaSTI absorption-line index grids

Fig. E.1. Comparison of the behaviour of the Mgb, CaHK, Ca4455, [MgFe], Fe5015, and NaD absorption-line indices as a function of the Hβ index. The shaded areas represent XSL PARSEC/COLIBRI models, with spectral resolution varying from σ = 13 km s−1 (the XSL native resolution) to σ = 60 km s−1 (the minimum E-MILES resolution in the optical). Black lines represent E-MILES BaSTI model predictions, with dotted, dashed, and solid lines representing [Fe/H] = +0.2, 0.0, and −0.4 dex, respectively, measured at the original E-MILES resolution. The solid blue line represents the predictions of the C18 solar-metallicity models.

Appendix F: M_*/L ratio model comparisons

Figure F.1 shows log(M_*/L) in the V, I, and K bands derived from the XSL PARSEC/COLIBRI Salpeter models, XSL Padova00 Salpeter models and E-MILES BaSTI Salpeter models. The M_*/L ratios for E-MILES models are taken from the model prediction tables provided by the MILES collaboration on their web page⁴ (Vazdekis et al. 2012; Ricciardelli et al. 2012; Vazdekis et al. 2016).

Fig. F.1. M*/L ratio derived from XSL PARSEC/COLIBRI (P/C) Salpeter models (orange), XSL Padova00 (P00) Salpeter models (green), and E-MILES BaSTI Salpeter models (black). The shaded areas represent the spread in M*/L arising from the metallicity spread from [Fe/H] = −0.4 dex (lower log M*/L) to [Fe/H] = +0.2 dex (higher log M*/L), centred on the solar-metallicity values.

The predicted M_*/L ratios differ among SSP models (assuming they use the same IMF). This discrepancy reflects, in part, a different accounting of stellar remnants and possible stellar mass span of the IMF. However, the biggest differences originate from the modelling of the advanced evolutionary stages of stars and found in the NIR. The XSL PARSEC/COLIBRI models include more stellar evolutionary phases, namely the TP-AGB phase, than the Padova00 models. These NIR-bright stars give their contribution to the light budget of a population, increasing its NIR luminosity and lowering its log(M_*/L_NIR).

This is clear from Fig. F.1, where the XSL PARSEC/COLIBRI models show a dip at log(M_*/L_K) 0.5-1 Gyr (coinciding with the SSP model ages where TP-AGB start becoming dominant; see Fig. 14), but XSL Padova00 and E-MILES BaSTI models do not. XSL Padova00 models monotonically decrease log(M_*/L) towards younger ages with the increase in the K-band flux fraction of TP-AGB stars, as seen in Fig. 14.

Over all ages, the XSL PARSEC/COLIBRI models have higher log(M_*/L) compared to XSL Padova00 models, due to having more mass in existing stars. The TP-AGB stars after the first thermal pulse are non-existent in Padova00 isochrones (replaced by remnants), but PARSEC/COLIBRI models still have them. The K-band luminosities of older XSL PARSEC/COLIBRI models are also lower, which is responsible for the increased M_*/L_K compared to XSL Padova00 models. Hence, the M_*/L_NIR ratio depends on the isochrones used as well as on the coverage of the spectral library and the modelling methods.

Appendix G: Examples of XSL PARSEC/COLIBRI SSP models

Fig. G.1. Examples of XSL PARSEC/COLIBRI Salpeter SSP models for four metallicities: [Fe/H] = +0.2, 0.0, −1.0, and −1.6 dex. We note that for [Fe/H] = −1.6 dex, models younger than 1 Gyr are not shown because they are outside the safe zones of XSL SSP models. Spectra are show in original resolution. The telluric absorption is marked in grey.

All Tables

All Figures

Fig. 1. 2 Gyr, [Fe/H] = −0.6 PARSEC/COLIBRI isochrone and the locations on the HR diagram where spectra are generated by local interpolation, global interpolation, or potentially taken from the ‘static’ giant, O-rich TP-AGB, or C-rich TP-AGB sequences. We only switch to the respective sequences when we have reached the bluest average spectrum on that sequence (according to the colour–temperature relation). For example, only the coolest RGB stars (Teff ⪅ 4000 K) are represented by a spectrum originating from the static sequence, and the spectra of warmer RGB stars are created by the global interpolator.

In the text

Fig. 2. Weighted median residual, RS, between the original spectrum and the interpolated spectrum for the full wavelength range, and for the X-shooter UVB, VIS, and NIR arms separately, as a function of position on the HR diagram. The colour bar is logarithmic. Histograms show the distributions of RS calculated within these spectral ranges at the full XSL resolution. For ease of visualisation, spectra with RS > 0.1 are placed into the RS = 0.1 bin in the histograms.

In the text

	Fig. 3. Spectral H−/H2O feature. It is a combination of the 1.6 μm ‘bump’ in the minimum opacity of H−ion and H2O vapour absorption bands around 1.4 μm and 1.9 μm. The index bands are marked in red. The telluric absorption is marked in grey. This spectrum is an average of O-rich TP-AGB stars with (I − K) = 3.19.
In the text

Fig. 4. H−/H2O index strengths of XSL stars as a function of their (I − K) colours. Points are colour-coded by their effective temperatures from Arentsen et al. (2019), saturated at 10 000 K. Carbon stars are marked with circles. Supergiants are marked with diamonds. Both C-rich TP-AGB stars and supergiants are excluded from the division into static and variable stars. M dwarfs have weak H−/H2O index strengths but can have red (I − K) colours, and they have been marked with triangles.

In the text

	Fig. 5. Sequence of static giant spectra, sorted according to (I − K) (on the left of each spectrum). The number of spectra the average is composed of is given in the brackets. The spectra have been smoothed to lower resolution for clarity.
In the text

	Fig. 6. Sequence of O-rich TP-AGB spectra, sorted according to (I − K) (on the left of each spectrum). The number of spectra the average is composed of is given in the brackets. The spectra have been smoothed to lower resolution for clarity. The reddest average spectrum on this sequence, displayed in grey, is not used in the stellar population models due to its extreme SED.
In the text

	Fig. 7. Sequence of C-rich TP-AGB spectra, sorted by their (R − H) colours (on the left of each spectrum). The number of spectra the average is composed of is given in the brackets. The spectra have been smoothed to lower resolution for clarity.
In the text

Fig. 8. XSL (‘PC’: PARSEC/COLIBRI; ‘P00’: Padova00), E-MILES P00, Maraston et al. (2009, M09), and Conroy et al. (2018, C18) SSP model spectra of 1 Gyr (panel a) and 10 Gyr (panel b) solar-metallicity stellar populations. The spectra are smoothed to R = 500. The M09 spectra are displayed at their original resolution (R ≈ 500). All spectra are normalised to a common I-band flux. The residual of the XSL P00 model from the E-MILES P00 model is shown in grey in each panel.

In the text

Fig. 9. Comparison of the behaviour of colour as a function of age and metallicity for the E-MILES, M09, C18, and XSL SSP models (see legend). Left panels: colour as a function of age. Right panels: colour as a function of metallicity. Top row: (B − V). Top-middle row: (V − I). Bottom-middle row: (I − J). Bottom row: (J − K). Ages span from 1 Gyr to 16 Gyr, and metallicities span from [Fe/H] = −0.4 dex to [Fe/H] = +0.2 dex (left panels); [Fe/H] = −2.2 dex (XSL PARSEC/COLIBRI) or [Fe/H] = −1.7 dex (other) to [Fe/H] = +0.2 dex (right panels). We note that E-MILES models are not safe to use in the NIR below [Fe/H] = −0.4 dex but are included for illustrative purposes. Solar-metallicity E-MILES models are shown in heavier line strengths than the sub- and super-solar models of the left panels. XSL SSP sub- and super-solar models are represented by shaded areas, centred on the solar metallicity. We note the different colour-scale values between the same colour panels.

In the text

Fig. 10. Comparison of the behaviour of model Mgb, CaHK, Ca4455, [MgFe], Fe5015, and NaD absorption-line index strengths as a function of the model Hβ index strength. The shaded area represents XSL Padova00 models with varying spectral resolution (σ) from σ = 13 km s−1 (the native XSL resolution) to σ = 60 km s−1 (the minimum E-MILES resolution). Black lines represent E-MILES Padova00 model predictions, with dotted, dashed, and solid lines representing [Fe/H] = +0.2, 0.0, and −0.4 dex, respectively, measured at the original E-MILES resolution. The solid blue line represents predictions from the C18 solar-metallicity models.

In the text

Fig. 11. Colour–colour diagrams of galaxies in the Coma cluster with model predictions overlaid. Grey points show the Eisenhardt et al. (2007) data for galaxies in the Coma cluster. Panel a: XSL and E-MILES Padova00 models. (b) XSL PARSEC/COLIBRI and E-MILES BaSTI models. Only [Fe/H] = +0.2, 0.0, and −0.4 dex E-MILES models are shown. C18 models with [Fe/H] = +0.2, 0.0, and −0.5 dex models are denoted by solid blue lines. XSL models extend to even lower metallicities.

In the text

Fig. 12. Selected index–index comparisons from Riffel et al. (2019). Shaded areas represent XSL PARSEC/COLIBRI SSP model predictions, with red, yellow, and teal indicating [Fe/H] = +0.2, 0.0, and −0.40 dex, respectively. Shaded areas represent models with a spectral resolution of σ = 16 km s−1 (the native XSL resolution) to σ = 228 km s−1 (the resolution of the Riffel et al. 2019 spectra) centred on σ = 60 km s−1 (the E-MILES NIR resolution). Black lines represent E-MILES BaSTI model predictions, with dotted, dashed, and solid lines representing [Fe/H] = +0.26, 0.06, and −0.35 dex, respectively (roughly the same metallicities as the XSL models). We note that the age range differs for the E-MILES models as the NIR spectra of E-MILES are only reliable above 1 Gyr. Panel b: the CN1.10* (Röck 2015) index definition is used instead of the Riffel et al. (2019) definition, which is affected by residuals from telluric absorption correction. We omitted the SFG and ETG measurements of CN1.1 due to differences in index definitions. Panels b–d: XSL static sequences, O-rich TP-AGB sequences, C-rich TP-AGB sequences, and XSL supergiants (which are not included in the XSL SSP models) are shown in grey at arbitrary optical index values (as these stars lack optical features) and their median values with larger black symbols. Indices of SFGs are marked in blue and indices of ETGs in red. These values are taken from Tables 6–7 and B1–B3 of Riffel et al. (2019), respectively.

In the text

	Fig. 13. Evolution of synthetic (log) stellar mass-to-light ratios for the V band (left), I band (middle), and K band (right), measured from the XSL PARSEC/COLIBRI Salpeter models.
In the text

	Fig. 14. Contribution of RGB stars and TP-AGB stars to the total K-band luminosity of the XSL models. Panels a and b: Contribution of RGB and TP-AGB stars, respectively, in the PARSEC/COLIBRI models. Panels c and d: contribution of RGB and TP-AGB stars, respectively, in the Padova00 models.
In the text

	Fig. 15. Colour–temperature relations for C-rich and O-rich stars. The grey horizontal lines (upper panel) represent the colours of the bluest and reddest average spectra on the C-rich TP-AGB sequence. The horizontal red and dark red lines represent the bluest and reddest spectra from the O-rich static and O-rich TP-AGB sequences, respectively.
In the text

	Fig. 16. 1 Gyr old, super-solar-metallicity XSL SSP generated using the Lejeune et al. (1997) (red) and Worthey & Lee (2011) (black) colour–temperature relations.
In the text

Fig. 17. Selected PARSEC/COLIBRI solar-metallicity isochrones illustrating the coverage of isochrones by the XSL DR3. The locations of the tip of the RGB for the 1 Gyr and 10 Gyr isochrones are indicated with horizontal black lines. We have employed the colour-temperature relations used in this study to add the giants to this figure to illustrate how well the cool RGB and AGB stars are represented (at arbitrary log g values). Because luminous cool stars are particularly important contributors to the red and NIR light of galaxies, XSL was designed to contain a large number of such objects. We are limited by the lack of hot stars and as such cannot create XSL SSP models younger than 50 Myr.

In the text

	Fig. 18. Fundamental parameter coverage of the XSL stars, for dwarfs (left panel) and giants (right panel). Dwarfs and giants are separated at log g = 3. The stars that are used to construct static, variable (O-rich TP-AGB), and carbon (C-rich TP-AGB) star sequences are not included, as we ignore the information about their metallicity in the models (Sect. 3.2).
In the text

	Fig. 19. Safe zones of XSL SSP models. ‘P/C’ indicates PARSEC/COLIBRI-based models, with white and grey indicating safe and unsafe zones, respectively. ‘P00’ indicates Padova00-based models, with the grey outline delimiting the safe and unsafe zones.
In the text

Fig. A.1. Weighted median residual, RS, between the original spectrum and the interpolated spectrum for four spectral line indices (CaHK, Hbeta, CaT, and CO1.6) as a function of position on the HR diagram. We normalised the spectra over the wavelengths of interest before calculating the line-level median residuals. The colour bar is logarithmic. Histograms show the distributions of RS calculated within these spectral ranges at the full XSL resolution. For ease of visualisation, spectra with RS > 0.1 are placed into the RS = 0.1 bin in the histograms.

In the text

	Fig. B.1. XSL spectra of O-rich, cool static giant stars from which the static sequence is constructed.
In the text

	Fig. C.1. XSL spectra of O-rich TP-AGB stars from which the variable O-rich TP-AGB sequence is constructed.
In the text

	Fig. D.1. Spectra of XSL C-rich TP-AGB stars from which the C-rich TP-AGB sequence was constructed.
In the text

Fig. E.1. Comparison of the behaviour of the Mgb, CaHK, Ca4455, [MgFe], Fe5015, and NaD absorption-line indices as a function of the Hβ index. The shaded areas represent XSL PARSEC/COLIBRI models, with spectral resolution varying from σ = 13 km s−1 (the XSL native resolution) to σ = 60 km s−1 (the minimum E-MILES resolution in the optical). Black lines represent E-MILES BaSTI model predictions, with dotted, dashed, and solid lines representing [Fe/H] = +0.2, 0.0, and −0.4 dex, respectively, measured at the original E-MILES resolution. The solid blue line represents the predictions of the C18 solar-metallicity models.

In the text

	Fig. F.1. M*/L ratio derived from XSL PARSEC/COLIBRI (P/C) Salpeter models (orange), XSL Padova00 (P00) Salpeter models (green), and E-MILES BaSTI Salpeter models (black). The shaded areas represent the spread in M*/L arising from the metallicity spread from [Fe/H] = −0.4 dex (lower log M*/L) to [Fe/H] = +0.2 dex (higher log M*/L), centred on the solar-metallicity values.
In the text

	Fig. G.1. Examples of XSL PARSEC/COLIBRI Salpeter SSP models for four metallicities: [Fe/H] = +0.2, 0.0, −1.0, and −1.6 dex. We note that for [Fe/H] = −1.6 dex, models younger than 1 Gyr are not shown because they are outside the safe zones of XSL SSP models. Spectra are show in original resolution. The telluric absorption is marked in grey.
In the text