© ESO 2020

1 Introduction

The effective temperature, T_eff, is one of the basic quantities characterizing a star and, in particular, one of the fundamental parameters shaping its spectrum. It is crucial for the accurate determination of chemical abundances because a biased T_eff has a strong impact on the abundances derived, much larger than the surface gravity (log g), for instance. Temperatures are often derived from empirical relations employing photometry by tying colors to temperature scales (e.g., Bell & Gustafsson 1989; Alonso et al. 1999), which were previously defined using, for example, the infrared flux method (IRFM, Blackwell & Shallis 1977). In this technique, the combination of the integrated flux and the infrared flux in a given band enables one to derive the effective temperature and the angular radius of the star.

However, such methods are also sensitive to other parameters (e.g., log g) and can be strongly affected by unusual chemical compositions (Huang et al. 2015, and references therein). Moreover, they require accurate determinations of the interstellar extinction, which is always more difficult in the case of disk stars and especially for Cepheids since they are pulsating and located close to the plane of the Milky Way.

As far as Cepheids are concerned, many studies have tackled the T_eff issue. For instance, Pel (1978) did so using Walraven photometry and Kurucz (1975) models; Fernley et al. (1989) used the IRFM method and Kurucz (1979) models; and Kiss & Szatmáry (1998) used Strömgren photometry and Kurucz models. Welch (1994), Fouque & Gieren (1997), and Kervella et al. (2004a) calibrated optical and near-infrared surface-brightness (IRSB) relations anchored to interferometric measurements of the radius of Cepheids. Kervella et al. (2004b) reported a remarkable agreement between the previous two Cepheid studies as well as with their own calibration of the surface-brightness relation for dwarf stars (Kervella et al. 2004c). Period-luminosity relations can then be determined using distances derived from IRSB relations (Storm et al. 2011a,b). The Spectro-Photo-Interferometry of Pulsating Stars code (SPIPS, Mérand et al. 2015) now combines direct interferometric measurements of the Cepheid angular diameter, and its variations along the pulsation phase, with all other observables available in the same Cepheid model (Breitfelder et al. 2016). The method has improved the accuracy and precision of several model parameters.

Below 8000 K, the effective temperature dependence of the hydrogen lines makes them excellent T_eff indicators, especially given their independence on log g. In particular, the wings of the Balmer lines provide one of the most reliable T_eff indicators for FGK stars (Fuhrmann et al. 1993; Hanke et al. 2018; Amarsi et al. 2018). However, in Cepheids, their distorted profiles (e.g., Nardetto et al. 2008) rule them out so far, in the absence of appropriate atmosphere models for pulsating stars.

The excitation equilibrium of Fe I and Fe II lines is commonly used to derive the effective temperature of FGK stars, including in our own analyses of Cepheids (e.g., Lemasle et al. 2007, 2008, 2017), under the assumptions of hydrostatic and local thermodynamic equilibrium (LTE). The first hypothesis obviously does not hold for pulsating variables and biases the results as recently verified for Cepheids by Vasilyev et al. (2017, 2018). The LTE hypothesis does not hold either (e.g., Fuhrmann 1998; Thévenin & Idiart 1999) and this is also true for Cepheids (e.g., Vasilyev et al. 2019). Alternatively, line depth ratios (LDRs) of pairs of lines with different sensitivity to T_eff have been put forward as temperature indicators by, for example, Gray & Johanson (1991) and Gray (1994). The number of ratios specifically calibrated to derive T_eff of Cepheids has increased from 32 ratios in Kovtyukh & Gorlova (2000) to 131 in Kovtyukh (2007) and 257 in Proxauf et al. (2018). LDRs pose highly precise observables that are capable of reproducing the relative temperature variations throughout the pulsation phase. Yet, an accurate absolute calibration of the scale is still pending (e.g., Kovtyukh & Gorlova 2000; Mancino et al. 2020). Employing spectroscopic T_eff indicators bears the major advantage that they are free of influences from reddening, an important asset in the case of classical Cepheids which are located close to the Galactic plane.

In this paper, we take advantage of the large sample of high-quality spectra that are available for Cepheids to establish a data-driven approach using flux ratios (FRs) as means to predict T_eff, in a similar approach as developed for stable stars by Hanke et al. (2018, henceforth H18) for the new code ATHOS (A Tool for HOmogenizing Stellar parameters). Section 2 describes the data preparation, in particular how labels, the meaningful information (here, T_eff) attached to each spectrum, were derived (Sect. 2.1). The method for selecting FRs (from the bluest wavelengths to the calcium triplet (CaT) region) is outlined in Sect. 3. We report our results in Sect. 4 and discuss them in Sect. 5.

2 Data preparation

2.1 Stellar parameters labels from the IRSB technique

The near-infrared (J, H, K) surface-brightness (IRSB) method (Fouque & Gieren 1997) is a Baade-Wesselink-type method that uses the pulsational variations in radial velocity, color, and luminosity to determine precise distances and radii of pulsating stars. It is based on the geometrical relation between the distance, d, to the star, the radius, R, and angular diameter, θ of the star at different pulsation phases, ϕ, namely: $$\theta (\varphi )=2R(\varphi )/d=2(R_0+\Delta R(\varphi ))/d.$$(1)

With many data points at different phases this linear relation can be trivially solved for the stellar radius and distance.

The radius at a given phase is obtained by integrating the pulsational velocity curve for the star where the pulsational velocity is determined from the observed radial velocity by applying a suitable projection factor, p. The p-factor is largely a geometrical factor but depends also on limb darkening, the exact method used to derive the radial velocity (see Borgniet et al. 2019), and other technicalities. Storm et al. (2011a) calibrated the p-factor empirically and provide a detailed discussion and references on this issue.

The surface-brightness in the V -band, F_V, is defined as: $$F_V(\varphi )=4.2207-0.1V_0(\varphi )-0.5\log \theta (\varphi )$$(2)

where V₀ is the de-reddened visual magnitude. It was shown by Fouque & Gieren (1997) that the surface-brightness is well described by a linear relation in the near-infrared ${(V-K)}_0$ color. Storm et al. (2011a) adopted the calibration by Kervella et al. (2004b): $$F_V=-0.1336{(V-K)}_0+\mathrm{3.9530.}$$(3)

This calibration is based on interferometric angular diameters of Cepheids.

To determine T_eff (ϕ) from the Stefan-Boltzmann equation, in addition to the stellar radii just determined, it is necessary to determine the luminosity at each phase point given the absolute V -band magnitude M_V, and the bolometric correction. M_V is determined from the V -band light curve given the distance from above and the reddening. The bolometric corrections have been computed from the 2012 revised (ODFNEW¹) ATLAS9 fluxes of Bessell et al. (1998) based on Castelli & Kurucz (2003) models, given the metallicity, log g, and T_eff. As T_eff is needed toenter the table it was necessary to make a few (usually only one) iterations to make T_eff converge.

The surface gravity log g depends on radius, the adopted mass for the Cepheid, and at any given phase also the additional acceleration from the dynamical atmosphere. A mass of 6 M_⊙ was adopted in Storm et al. (2011a) for all the stars. The exact choice does not make a significant difference on the derived temperature. Had a mass of 8 M_⊙ been adopted then the value of T_eff would have changed by less than a Kelvin.

The adopted metallicity also affects the bolometric correction and thus the resulting T_eff but the effect is very weak. An error in the adopted metallicity of 0.5 dex would lead to an error in T_eff of only 15 K.

Similarly, the adopted p-factor for the conversion between the observed radial velocity and pulsation velocity has no effect on the derived temperature as the resulting change in radius is exactly compensated by the changed distance and thus in the luminosity.

As shown already by Fouque & Gieren (1997), the distance and radius estimates are very robust to errors in the adopted reddenings. On the other hand, M_V depends directly on the reddening so the estimated temperatures are indeed sensitive to the adopted reddening. The effect of artificially lowering the adopted reddening, E(B − V), by 0.05 mag is to reduce the estimated temperature by 150 to 200 K depending on the pulsation period of the Cepheid. This is the only significant systematic uncertainty apart from the adopted surface-brightness relation.

By construction (see Sect. 2.5), our spectroscopic dataset mostly consists of nearby Cepheids for which a wealth of external data is available. In particular, the values of the reddening, E(B − V), are as low as 0.003 mag for ζ Gem and below 0.2 mag for 46% of our sample. This ensures that the values of the reddening are either negligible or well-constrained in order to limit systematic errors on the labels. Our original sample also includes Cepheids in the Large Magellanic Cloud (LMC) and in the Small Magellanic Cloud (SMC) for which the E(B − V) values were determined using the reddening maps of Górski et al. (2020). E(B − V) values range from 0.081 to 0.175 mag for our LMC Cepheids and from 0.058 to 0.103 mag for our SMC Cepheids, with typical uncertainties of 0.005 mag.

The reddening value adopted is crucial for the determination of T_eff via the IRSB method. A simple test for the 7.6 day period Cepheid W Sgr showed that the estimated temperature is reduced by 190 K when E(B − V) is artificiallylowered by 0.05 mag with respect to the adopted value of 0.108 mag, whereas for the 41 day Cepheid RS Pup such a change of the reddening reduces the estimated temperature by 150 K. The reddening values adopted in our IRSB analysis, used in Storm et al. (2011a) and based on the Fouqué et al. (2007) compilation, match closely those tabulated by Laney & Caldwell (2007). They are on the photometric system of Laney & Caldwell (2007) as described inFouqué et al. (2007). Laney & Caldwell (2007) showed that externally determined reddenings (for instance, using main-sequence stars that either are companions of the Cepheid or belong to the same open cluster) fall on the same scale as those derived from multicolor photometry, hence providing a homogeneous reddening scale.

Fouqué et al. (2007) attempted to combine Cepheids’ reddenings obtained from the 1970’s on by different groups onto the system of Laney & Caldwell (2007) using linear transformations. From their Table 1, the scatter around the different transformations reaches 0.04 mag, while for the largest sample in that study (Fernie et al. 1995), it is only 0.03 mag. From the same table, it appears that the Dean Warren and Cousins system (Dean et al. 1978) is quite close to the Laney & Caldwell (2007) system, while most of the other studies² seem to cluster around similar values with a systematic offset of ~0.04 mag. Assuming two types of systems, one could claim systematic uncertainties of half of the offset or 0.02 mag. This is probably a too optimistic view and we adopt instead a more conservative value of 0.03 mag. It is worth mentioning that Turner (2016) using space reddenings reported that the values from Storm et al. (2011a) should be transformed linearly with an offset of +0.05 mag for a typical E(B − V) = 0.5 mag. Numerous studies have focused on the reddening scale in the Magellanic Clouds. We used the maps obtained by Górski et al. (2020), that have a systematic uncertainty of the average reddening value of 0.013 mag. A comparison with previous studies in the LMC, for instance, indicates a mean difference of 0.061 mag, with a standard deviation of 0.012 mag between the reddening values from Górski et al. (2020) and those from Haschke et al. (2011). The same comparison between the results of Górski et al. (2020) and those of Inno et al. (2016) leads to an average offset of 0.004 mag with a standard deviation of 0.087 mag.

Reddenings are by far the main source of systematic uncertainties in our study, and improving the accuracy of Cepheids’ reddenings will provide an exquisite accuracy for our flux ratio method.

2.2 Spectroscopic data

For training and testing the method, we collected a large number of high-resolution spectra at high signal-to-noise ratios (S/N). We limited ourselves to a small number of spectrographs, namely HARPS (Mayor et al. 2003) at the La Silla 3.6m telescope, HARPS-N (Cosentino et al. 2012) at the Telescopio Nazionale Galileo (TNG), UVES (Dekker et al. 2000) at the Very Large Telescope (VLT), FEROS (Kaufer et al. 1999) at the La Silla 2.2m telescope, and STELLA (Strassmeier et al. 2004; Weber et al. 2012) at the 1.2m telescope at Izana Observatory. We considered only spectra with S/N in excess of 50 pixel⁻¹ for Galactic Cepheids, but had to weaken this requirement to a S/N of 30 pixel⁻¹ for Magellanic Cepheids in order to be able to incorporate them in our sample³. Table 1 lists the properties of the spectrographs and the number of spectra gathered for each of them. Table 2 lists the Cepheids in the sample and the number of spectra gathered for each instrument.

Spectra that were taken with FEROS, HARPS, and UVES were retrieved from the ESO (European Southern Observatory) archive. The spectra taken with HARPS-N were taken from the TNG archive. HARPS and HARPS-N products were used in the 1d spectral data format where individual spectral orders are merged. The STELLA spectra were collected by J. Storm. Since coaddition of the spectral orders is not included in the STELLA pipeline, we performed this operation using a S/N-weighted sum in the spectral regions where orders overlap.

One of the strengths of our FR-based approach is that it does not rely on normalized spectra by restricting itself to small wavelength ranges for which the assumption of constant continua holds. For preparing the data, we thus only had to place the spectra in the stellar rest frame: we corrected the spectra for line-of-sight motions using a cross-correlation method with three different, synthetic template spectra that are representative of Cepheids and that were selected according to the temperature of the star at the time of the spectroscopic observation. The templates’ characteristics are listed in Table 3. Finally, we performed a simple cosmic ray rejection by interpolating over features showing strong variances on scales much smaller than the width of the instrumental profile.

2.3 Training at a resolving power of R = 5000: pros and cons

All the spectra were then rebinned to a common resolving power of R = 5000. Training at such a low resolution limits the number of available relations when compared to training at higher resolution (R ≫ 10 000). This is not a problem since the accuracy remains driven by the intrinsic accuracy of the labels, and the precision remains better than any other method available.

Training at R = 5000 comes instead with two indisputable advantages. Firstly, the already high S/N of our spectra reaches extremely high (>1000) values when down-sampled to R ~ 5000, which makes the training insensitive to photon shot noise. This also applies to unseen data, since high-resolution spectra (outside of the training sample) would need to be degraded to the resolving power R = 5000 in order to derive T_eff using the method developed in this study. Secondly, in the perspective of large spectroscopic surveys with, for instance, WEAVE⁴ (Dalton 2016), 4MOST⁵ (de Jong et al. 2019), or MOSAIC⁶ (Evans et al. 2015), it provides the same temperature scale for samples studied at low (on the order of R = 5000) or medium resolving power (on the order of R = 20 000).

However, training at lower resolving power further comes with the caveat of being in all generality more metallicity sensitive because of blending features. This is not problematic within the range of parameters covered by the training sample, as it is intrinsically accounted for, but limits the applicability of the method to this exact same range (which is anyway always the case when using machine-learning techniques). In case of varying [X/Fe] among targets, this very circumstance poses an additional potential bias. It should be noted that in Cepheids, despite star-to-star variations, Genovali et al. (2015) reported that [α/Fe] shows a flat distribution across the entire Galactic thin disk (as one would expect for a young population), and no trend with log P. As far as neutron-capture elements in Cepheids are concerned, da Silva et al. (2016) found a positive slope for [La/Fe], [Ce/Fe], [Nd/Fe], and [Eu/Fe] but no slope for [Y/Fe]. The same authors also reported no trend (except maybe for Ce) as a function of log P. Including a large number of Magellanic Cepheids in the sample is a strong requirement for improving the universality of the method.

At the beginning of Sect. 2.3, we have exposed our reasons for training at R = 5000. Determining T_eff for an unseen Cepheid spectrum at any resolving power R > 5000 then requires that this spectrum is first degraded to the training resolving power (as mentioned above, no continuum normalization is required). It goes without saying that the method can therefore be applied only to spectra with R ≥ 5000. By construction, it does not loose generality with varying resolution or T_eff. This could instead be the case for LDRs (see Mancino et al. 2020) because some line could disappear or become blended at a given T_eff which may itself depend on the resolution of the spectrum.

We set a higher priority on our main goal (providing a homogeneous T_eff scale for low- and high-resolution Cepheid spectra in upcoming surveys) than on potentially optimizing accuracy and precision via a different choice for the training resolution. Moreover, in the perspective of a canonical chemical analysis, we are fully satisfied with the accuracy and precision reached, as is detailed in Sect. 4. A detailed analysis of the performances of flux ratios as a function of resolution is, therefore, not undertaken in this paper. We anticipate that training at higher resolution will increase the number of calibrating relations and therefore the precision.

2.4 The need for an accurate phasing of the spectra

We quickly realized during the training process that even small errors in the phasing (on the order of a few hundredths of the period) of the training spectra could lead to relatively large biases (up to a few 100 K in the labels, see Fig. 1), and that these biases would be the parameter limiting the accuracy of our temperatures. Here we want to stress that the phase of the observations is only used to attribute T_eff labels to the spectra employing values derived via the IRSB method, and is not used further in the process.

The phase of the spectroscopic observations was initially determined using the Modified Julian Date (MJD) at which the spectra were taken, together with the exact same values of reference epoch and pulsation period that were used in the IRSBanalysis. An extremely large (and constantly updated) collection of radial velocity measurements (spanning three decades) has been collected to apply the IRSB method (see, for instance, Storm et al. 2011a). To ensure that our high-resolution spectra (spanning 15 yr) are placed on the exact same phase scale as the IRSB analysis, we checked our own radial velocity measurements against the IRSB radial velocity curves and phase shifts were applied to ensure a consistent phasing between both datasets. It is no surprise that they concerned primarily:

• stars for which the high-resolution spectra were taken at very different epochs from the spectra used for the IRSB analysis;

• Cepheids with known period changes. These period changes can be related to secular evolution of the Cepheids across the instability strip (e.g., Berdnikov et al. 2019) or to cycle-to-cycle variations that affect, in particular, long-period Cepheids (Anderson et al. 2014, 2016; Anderson 2016).

An alternative approach would have been to take into account period drifts. However, data are often too scarce to derive them accurately, and it is usually difficult to reach a consensus on the value of even a simple, linear period drift (see Breitfelder et al. 2016, and references therein). To our knowledge, RS Pup is the only Cepheid for which the modulation of the period over several decades has been derived accurately, using a fifth-degree polynomial (Kervella et al. 2017). Given the current insufficient accuracy of period drifts, we adopted instead the method described above. We cannot exclude that a fraction of the remaining residuals discussed in Sect. 4.3 originates from a still inaccurate phasing of some spectra.

For a few stars, an accurate rephasing based on the radial velocity curve proved impossible. It is the case of the exceptional Cepheid X Sgr, where the signature of strong shock waves (emission, line doubling) has been reported by Mathias et al. (2006). In a few other stars (S Mus, EV Sct, U Aql, V Car, V Cen, and V350 Sgr), a strong orbital motion was detected, which limited the accuracy of the phasing, and these stars were therefore not further considered despite a good number of spectra available. Finally, a few other spectra were rejected because the radial velocity curves of the stars were of lower quality (possibly due to orbital motion as well) and/or with incomplete phase coverage (KQ Sco, RZ Vel: spectra between phases 0.9 and 0.1, S Nor, T Vel, UZ Sct, WZ Car, XX Sgr, Y Sgr).

Fig. 1 Impact of erroneous phasing of individual spectra on the Teff labels for the Cepheid RS Pup: top panel: MJD of the observation for individual spectra. Second panel: our individual radial velocity measurements (circles that are color-coded by MJD, cf., upper panel) versus radial velocity data used for the IRSB analysis (blue crosses, Storm et al. 2004). Spectra rephased to match the radial velocity curve are shown as pale brown squares. Using the same color coding, the third and fourth panels display the Teff variation, and the Teff offset between original labels and their corrected value after rephasing, respectively.

2.5 The dataset

A good dataset must have the following characteristics:

• the labels are accurate;

• they cover the entire parameter space and are homogeneously distributed, not only in label space (T_eff), but also in other parameters that may affect an FR (e.g., [Fe/H] and log g);

• a large number of spectra is available.

However, as with many if not all data-driven methods, it is close to impossible to generate such a dataset. Here we investigate whether the conditions listed above are satisfied.

We retained 24 Galactic Cepheids for our sample. They are listed in Table 2, together with the number of spectra used for each of them. We considered only stars for which a good number (>10) of spectra were available and for which accurate phases could be determined (see Sect. 2.4). We note in particular that a few stars are represented by a considerable number of spectra, namely l Car (191), ζ Gem (135), β Dor (125), δ Cep (122), and RT Aur (95). We end up with a total of 1324 spectra.

Although in the Magellanic Clouds only few spectra per star are available, 32 spectra of 13 Magellanic Cepheids were also added to the dataset. Including a large number of spectra of LMC and especially SMC stars in the sample would present the advantage that metallicity-dependent relations would naturally be excluded⁷. However, for these stars we only have few spectra, which in addition do not cover the entire wavelength range. This means that temperatures derived from flux ratios located in different spectral domains would be differently affected by possible metallicity dependences. This does not seem to be a satisfying outcome. We decided instead to train at high metallicity only and to check for T_eff consistency of individual relations at low metallicity (see Sect. 5.1).

The first condition we set for our sample is met: we selected bright, nearby Cepheids that were thoroughly investigated over the years. The radial velocity curves and the light curves are well sampled. The spectra were phased with great care. As a result, accurate labels could be derived following the method described in Sect. 2.1.

Our sample covers a broad range of periods (3.7 days ≤ P ≤ 78.8 days, see Table 2), which ensures a good coverage of the T_eff - log g space. Indeed, the coolest stars in our sample, l Car and RS Pup, have temperatures ranging from as low as 4578 to 5282 K (respectively from 4621 to 5560 K). The hottest Cepheid in our sample, RT Aur, has temperatures ranging from 5573 up to 6592 K. For these stars, log g varies between 0.41 dex and 1.30 dex (l Car), between 0.07 dex and 1.39 dex (RS Pup), and between 2.16 dex and 2.50 dex (RT Aur). The other Cepheids have intermediate values for T_eff and log g (Fig. 2).

However, the condition is not met as far as [Fe/H] is concerned. The best-studied stars are naturally among the closest Cepheids (which also allowed to collect large amounts of spectra in a reasonable amount of observing time). Given the presence of abundance gradients in the Milky Way (e.g., Luck et al. 2011; Yong et al. 2012; Genovali et al. 2013), those Cepheids have metallicities close to Solar. Cepheids in the Milky Way span almost 1.0 dex in [Fe/H], from ≈+0.5 dex in the inner disk (e.g., Pedicelli et al. 2010; Andrievsky et al. 2016) to ≈–0.5 dex in the outer disk (e.g., Lemasle et al. 2008; Luck & Lambert 2011; Genovali et al. 2014). Due to their distances (>3 kpc) and their location close to the plane, less spectra are available for the stars located at the edges of the metallicity range, and the IRSB data is also scarce. To extend the metallicity coverage, a good alternative is to turn to the Magellanic Clouds, where a few tens of Cepheids have been studied spectroscopically at high resolution (Romaniello et al. 2008; Lemasle et al. 2017) and for which IRSB analyses are available (Storm et al. 2011b; Gieren et al. 2018). [Fe/H] roughly ranges from –0.2 to –0.6 dex in the LMC and from –0.6 to –0.9 dex in the SMC. Unfortunately, very few medium- or high-resolution spectra are available, and even less when requiring S∕N > 30 pixel⁻¹. Moreover, all stars are represented by a unique spectrum taken at random phase, which means that we do not have at our disposal spectra covering a large range of T_eff, log g. and periods for low-metallicity Cepheids.

The third condition is largely met: we have gathered a large number of spectra for each of the training Cepheids (see Table 2), namely a grand total of 1324 spectra, which provide a full phase coverage for each of the selected stars. The numberof nonpulsating benchmark stars for covering the entire Hertzsprung-Russell diagram is considerably lower (Jofré et al. 2014, Gaiabenchmark stars; Hanke et al. 2020, HD 20).

Fig. 2 Parameter space (Teff, log g, [Fe/H]) covered by the compiled sample of spectra. Highlighted by salmon and red circles are the Cepheids ζ Gem and l Car, for which we have the most spectra at hand (135 and 191, respectively). Gray squares indicate a small portion of the giants used in the ATHOS dataset for stable stars (H18).

3 Method

The approach pursued throughout this work closely follows the one described in H18. To this extent, in order to enable the computationally efficient handling of the data, all spectra were binned to a common log-linear dispersion scale, where the nth element is expressed as $${\log }_{10}{\lambda }_n={\log }_{10}{\lambda }_0+n\cdot \delta {\log }_{10}\lambda ,\text{for}n\in \{0,\mathrm{...},N\},$$(4)

with λ₀ being the starting wavelength of 3800 Å, δlog₁₀λ the dispersion spacing, 2.7 × 10⁻⁵, and N the total number of pixels, here taken to be 13 333. This way, our maximum considered wavelength is ~ 8705 Å. Employing a log-linear λ scale bears the main advantage that profile widths remain approximately constant in pixel space⁸ over the extensive range used here.

3.1 Feature selection

We recall here the definition from H18 of a spectral FR between two wavelength windows that are centered at λ_i and λ_j, respectively: $${\text{FR}}_{{\lambda }_i,{\lambda }_j}=\frac{\langle F_{{\lambda }_i}\rangle }{\langle F_{{\lambda }_j}\rangle },$$(5)

where $$\langle F_{{\lambda }_i}\rangle =\frac{1}{5}{\sum }_{n}F({\lambda }_n),n\in \{i-2,\mathrm{...},i+2\}$$(6)

is the mean flux of five pixels.

Each FR potentially qualifies as feature that encodes information about T_eff, which in principle leaves us with $\mathcal{O}({10}^{10})$ explanatory variables. This enormous number calls for an efficient and robust feature-selection process prior to the actual training and validation of feasible relations. Our first, and arguably most important, preselection criterion considered only those FRs for which the two flux windows are not separated by more than 30 pixels. While reducing the number of features by two orders of magnitude, this limit makes the FRs essentially scale-free, that is, normalization-independent, since it ensures that the continuum variation between the two closely-spaced regions of interest can be safely neglected (see Sect. 4.2). An additional applied criterion was that only those FRs were allowed as features that showed a maximum range of values in excess of 0.15. This was enforced to ensure that noise has a limited impact on the predicted labels from each FR (see also Sect. 4.1).

As a next step, we computed Pearson correlation coefficients, ${\rho }_{\text{FR},T_{\text{eff}}}$, between all remaining FRs and the T_eff labels associated to the underlying spectra. By requiring strong (anti-)correlations (i.e., $\left|{\rho }_{\text{FR},T_{\text{eff}}}\right|>0.93$) between the features and the labels, we further restricted the feature space to 33 291 candidates. We emphasize at this point that by proceeding like this we potentially rejected a large number of tight feature-label relations, in particular those that are strongly nonlinear and/or nonmonotonic. This is owed to using computationally cheap correlation coefficients that test for linearity and could only be circumvented by investing substantially larger computational resources.

Fig. 3 Temperature distribution of the spectra in our dataset, shown as a gray histogram. The temperature distributions of the training and testing samples are shown in blue and pale brown, respectively. We also indicated in red the temperature distribution of the spectra of stars with |Fe/H] < −0.3 dex that did notenter the training.

3.2 Training and testing

For fitting, the dataset was split into a training set (70% of the sample) and a holdout (or testing) dataset (30%), which was used for unbiasedestimates of the model fitness. In principle, it would be desirable to have at the same time both a homogeneous coverage in label space and the training and testing sets obeying identical probability distributions. In practice, however, since our T_eff distribution is skewed toward lower temperatures (see Fig. 3), we decided to foster T_eff -homogeneity for the training sample by randomly drawing with the inverse temperature density as weight. As a consequence, the holdout set follows a slightly different probability distribution, which is indicated in Fig. 3. The only alternative approach would be to disregard a substantial number of spectra in the intermediate- and low-temperature regime from the analysis to establish equal numbers throughout the whole T_eff range. Yet, we favor our choice for two reasons: first, some spectra have a discontinuous wavelength coverage due to chip gaps and others are limited to the visual range and do not extend to the near-infrared CaT region (8498–8662 Å). Hence, not all FRs can be measured in all spectra, which adds an additional layer of complication that is best mitigated by having a higher number of spectra from different spectrographs. The second reason for keeping all observed spectra is that it enables a better coverage in fundamental parameters other than T_eff, such as [Fe/H], log g, or period, which are key to spotting degeneracies with these higher-order parameters.

Even though we used correlation coefficients for feature selection and therefore technically prefiltered our features for close-to linear relations, we allowed our model to be more flexible by adding the possibility of curvature through second-order polynomials. It is noteworthy at this point that having such a simple model renders a validation set for tuning hyper-parameters superfluous. One such case is presented in Fig. 4, where we, on the one hand, exemplarily demonstratehow our most well-constrained spectral region (see next section for the selection criteria) around ~ 4405 Å behaves with temperature, and, on the other hand, we illustrate how the identified FR constrains T_eff both in the training and testing sets.

The actual training was performed on each of the 33 291 FRs that survived the feature selection. While the model was fit using thetraining set, its fitness in terms of accuracy (mean deviation) and precision (standard deviation, σ) was independently determined from the testing set. To lessen the impact of artifacts on individual FRs induced by, for example, erroneous phasing, wavelength calibration errors, and cosmics, we allowed a clipping of up to five of the most extremely deviating spectra both from the training and testing sets. The distribution of the respective findings for the individual features is depicted in Fig. 5.

We found that there is a slight accuracy bias that tends toward negative values with decreasing precision. This could be owed to the circumstance that for any set of two physically identical spectra in the training and testing sets, respectively, the testing spectrum shows an offset T_eff label. Such an effect could occur if these spectra are associated with different Cepheids, where inter-target systematics (e.g., reddening, projection factor, etc.) could hamper the labeling in the IRSB analysis. Nevertheless, given that the systematic bias remains wellbelow 25 K, we do not investigate this behavior further (but see Sects. 4.3 and 5.2).

Not surprisingly, the precision correlates with the correlation coefficient, although – due to allowing for curvature in the model – we deem the precision to be the more desirable quantifier for the goodness of the model fit.

Fig. 4 Properties of an individual relation using flux ratios around 4405 Å. Top left panel: evolution of the neighboring spectral domain with Teff, where darker colors indicate deeper features. Bottom left panel: coolest (dark brown) and hottest (orange) spectra at our disposal. For illustration purposes, the spectra were pseudo-normalized by the 99th percentile of the flux in the presented range. In both panels, the regions considered for computing the flux ratios are marked in clear blue. Top right panel: analytical relation (dashed black curve) together with the Teff labels for the training (blue) and testing (pale brown) samples. Using the same color coding, bottom right panel: difference between the temperatures computed using the individual analytical relation and the input labels. In both panels, the coolest and hottest spectra are again depicted in dark brown and orange, respectively.

Fig. 5 Precision and accuracy (in K) for all the relations retrieved by the analysis, color-coded by their correlationcoefficient. Only relations for which the precision is better than 100 K and less than 3 × the accuracy are further considered (red triangle). The dashed line shows a slight accuracy bias (see text) as a function of precision, which was computed using a boxcar median filter.

3.3 Final selection of flux ratios

After having fit all candidates for useful features, we imposed further quality cuts, in a sense that only such features were kept, for which the model assessment yielded precision∕accuracy > 3 and the precision remained below 100 K. These criteria are superimposed in Fig. 5 and resulted in 12 017 FRs for further analyses.

Next, we had to make sure that the final selection of features is linearly independent. Only then, the predicted labels from all models and unseen data can be averaged in a meaningful way, that is, without having to deal with inevitable (unknown) covariances⁹. We achieved this by ranking the features by their associated model precision and subsequently exclude those that have pixel overlaps with a feature of higher precision. This procedure further reduced the number of leftover features to 247, which indicates that the vast majority are clustered around a limited number of T_eff -sensitive line profiles and are only shifted by incremental λ steps with respect to each other.

Finally, the remaining 247 relations were visually checked for erratic behaviors in the residuals T_{eff, IRSB} − T_{eff, FR} to exclude (strong) influences on the FRs from parameters in excess of T_eff. Another considered factor was substructure that is indicative of an insufficient description by our quadratic model (e.g., discontinuities or nonflat residuals as a function of FR). For our best relation, the residuals are presented with respect to ϕ, log g, [Fe/H], and log P in Fig. 6. There, we also highlight the Cepheids ζ Gem, l Car, and δ Cep as representatives with good phase and log g coverage as well as very different pulsation periods. It is evident that there is a persistent, statistically significant bias of ~ 59 K between ζ Gem and l Car (⟨ΔT_eff⟩ = 56 ± 16 K and − 3 ± 19 K, respectively, both only computed from the testing set) irrespective of phase and log g. We would like to stress at this point that this shows that what we call precision of a relation (i.e., the scatter measured from test spectra of many different Cepheids) is largely governed by star-to-star systematics in the labels and therefore encompasses external errors. In fact, the internal precision of a single of our FR relations as deduced from one target alone (in the example above ζ Gem and l Car) resides below 20 K. The reader is referred to Sect. 4.3 where we discuss in detail the implications of labeling systematics on the full sample of relations. After our visual sanity check, we are left with 143 relations.

Fig. 6 Residuals Teff, IRSB − Teff, FR as a functionof phase (ϕ, top left), log g (top right), [Fe/H] (bottom left), and log P (bottom right), for the same relation as in Fig. 4. Colored open circles highlight spectra for the targets ζ Gem, l Car, and δ Cep.

4 Results

After carefully selecting and phasing the dataset, we trained our relations using IRSB T_eff labels and obtained 143 calibration relations that can be used to determine the effective temperature of Cepheids. Before describing the properties of the identified relations, we want to test the stability of individual relations with respect to the S/N and continuum variations, and briefly discuss the remaining systematics on the final temperature.

4.1 Signal-to-noise ratio stability

Given a high-S/N (e.g., >100 Å⁻¹) spectrum, even a single of our discovered features is capable of predicting T_eff at a precision that is more than sufficient for subsequent abundance analyses (i.e., ${\sigma }_{T_{\text{eff},\text{FR}}}\ll 100$ K). In practice, for larger samples this is only feasible for nearby and therefore bright Cepheids for which spectra can be obtained within reasonable exposure times and/or 8+ m-class telescopes. In light of the fact that upcoming spectroscopic campaigns such as WEAVE and 4MOST at 4 m-class facilities will only allow for a more restrictive spectrum quality (S∕N ≳ 30 Å⁻¹), our extended number of relations with a wide wavelength coverage is paramount. In order to optimize the averaging strategy for predicted labels from multiple features, knowledge about the S/N stability is desirable.

The parameter affecting the noise sensitivity most is ΔFR_{min, max}, which denotes the excursion an FR covers over the full range of T_eff. To test this, we have conducted a Monte Carlo (MC) analysis using a single-epoch spectrum of ζ Gem with attributed parameters (ϕ = 0.48, T_eff = 5206 K, log g = 1.29 dex, [Fe/H] = +0.10 dex) that are roughly consistent with the midpoint of the parameter coverage of our dataset (cf., Fig. 2). The spectrum was obtained with R ≈ 115 000 and S∕N = 276 pixel⁻¹ using HARPS. At 0.01 Å per pixel, this latter value corresponds to S∕N ≈ 2760 Å⁻¹, which leads us to the conclusion that for tests with artificial noise satisfying S∕N ≤ 1000 Å⁻¹ contributions from the real observed noise level can be safely neglected. To estimate both mean values and standard deviations at various S/N values between 25 and 1000 Å⁻¹, per S/N realization, we generated 500 spectra with injected random noise, degraded each to R = 5000, measured FRs, and predicted the T_eff labels.

The result can be seen in Fig. 7 where we indicate the course of ${\sigma }_{T_{\text{eff},\text{FR}}}$ as a function of S/N for three representative features with different realizations of ΔFR_{min, max}; the lowest allowed spacing of 0.15, an intermediate value (0.28), and a large excursion (0.50). For these three relations, the noise-induced standard deviation ranges from (366 K, 177 K, 77 K) at S∕N = 25 Å⁻¹ through (90 K, 47 K, 22 K) at S∕N = 100 Å⁻¹ to (9 K, 5 K, 2 K) at S∕N = 1000 Å⁻¹. Hence, even at S∕N = 25 Å⁻¹, and from a mathematical point of view – that is, neglecting artifacts that real spectra inevitably suffer from – averaging only 16 of our relations with the strongest (worst) S/N sensitivity already yields internal temperature precisions that are sufficient for meaningful abundance analyses ($366 \text{K}/\text{}\sqrt{16}=92$ K).

An important – though not unexpected – observation is that, considering the scatter, the mean T_{eff, FR} is perfectly stable with S/N. This highlights the superiority of our method over many already existing analysis tools that show systematically increasing over- or underestimations of the continuum level with varying S/N. As a consequence, further down the line such tools introduce potentially strong systematic scatter when comparing results (e.g., [Fe/H] or other chemical abundances) obtained from spectra of substantially varying quality (see for instance Fig. 2 in Reichert et al. 2020).

Fig. 7 Stability against S/N per Å for three relations of different maximum FR excursion (cf., color coding and legend in the lower panel). The corresponding S/N per resolution element for four different resolving powers (see legend) are given in the top panel. Shown are the measured mean values (colored dashed lines) and 1 σ standard deviations (solid lines enclosing filled areas) from an MC simulation with artificially injected noise (see main text for details). The zoom-in view emphasizes the (marginal) relation-to-relation inaccuracies (see middle panel of Fig. 10). Vertical lines indicate the minimum S/N requirements for different 4MOST consortium surveys: S∕N = 50 Å−1 in the magnesium or calcium triplet regions (4MIDABLE-LR10); S∕N = 30 Å−1 in the magnesium or calcium triplet regions (1001MC11, LR); S∕N = 100 Å−1 at 6200 Å (1001MC, HR; 4MIDABLE-HR12).

4.2 Stability against continuum variations

Another investigated source of error is the stability of T_eff predictions from our models against extreme variations of the continuum level between two involved wavelength ranges due to, for instance, a mistreated removal of the instrument’s blaze function. As opposed to noise-induced errors – which can be reduced by averaging a number of relations – potential errors from this are of systematic origin and ultimately affect our accuracy. In Fig. 8 we illustrate how the effect differs for two FRs at the extreme ends of the distribution of wavelength spacings, Δλ. We employed the same spectrum of ζ Gem introduced in the previous section and multiplied linear continua with slopes as drastic as ± 0.1 % Å⁻¹. Not even thesteepest slopes exceed a ± 25 K effect for the closely-spaced case, whereas a maximum excursion of ± 50 K was found for the largest Δλ. One might desire to drop or downweight the ratios with largest Δλ when analyzing high-dispersion spectra with narrow orders, since they are potentially affected by somewhat stronger continuum variations than other devices.

Fig. 8 Impact of continuum variations on FRs from close- (red) and wide-spaced (blue) wavelength components. Upper panel: 15 times the same ζ Gem spectrum as was used in Sect. 4.1 with varying continuum slopes (see color bar in lower panel). The missing part in the presented spectrum corresponds to HARPS’ chip gap. Middle panel: zoomed-in view of the region around ~ 5295 Å where the λ regions are highlighted by vertical spans. Lower panel: deduced Teff, FR values as a function of the measured FRs at different continuum slopes. The models are represented by red and blue solid curves while dashed red and blue lines and filled circles of the same colors indicate the finding from the original spectrum without any modification of the continuum.

Fig. 9 Representative examples of the three distinct residual patterns as discussed in the text: δ Cep (IRSB issuesseemingly constant with MJD, left panel); l Car (IRSB issueswith MJD-dependence, middle panel); ζ Gem (constant offset possibly due to erroneous reddening, right panel). Shown are the median residuals of the best 100 trained relations that were measured in individual spectra with respect to the input labels. The abscissa indicates the pulsation phase, ϕ. As in Fig. 1, the color coding resembles the time of observation of the individual spectra (cf., upper panels). Typical values for the 1 σ scatter are depicted by gray error bars. The inlay in the middle panel highlights five spectra of l Car that were taken with HARPS back to back within less than 30 min. The rms scatter of 1.6 K underlines the internal precision of our method.

4.3 Size and origin of remaining systematics

In the two previous sections, we investigated individual relations. From now on, we are dealing with the average T_eff derived from a largenumber of relations, and not anymore from a temperature derived from a single relation.

Thanks to our efforts to minimize systematics, we are left with small-scale residuals (a few tens K) that show different patterns, displayed in Fig. 9. The first pattern is a phase-dependent offset that seems independent from the MJD of the observation, at least for the relatively small range of MJDs considered here (case of, e.g., δ Cep). The residuals arelikely related to the IRSB method. We checked, however, that there is no phase shift between the optical and near-infrared photometry, nor with the radial velocity data (see Storm et al. 2011a,b, and references therein). Since the peak of the deviation happens at maximum and minimum radii, we also investigated whether the radius amplitude could be influenced by the value selected for the p-factor. Using p = 1.3 instead of the nominal value (p = 1.41), we found that the effect on T_{eff, IRSB} is less than 2 K over the entire phase range,thus excluding the p-factor from the possible source for the residuals. The second pattern is a similar phase-dependent offset, this time MJD-dependent (case of, e.g., l Car). Although the residuals are again likely related to the IRSB analysis, their variation probably originates from an inaccurate rephasing or from cycle-to-cycle variation of the temperature curve. The third pattern resembles a constant offset between the retrieved temperature, T_{eff, FR}, and the labels (case of, e.g., ζ Gem), possibly dueto an erroneous value for a global, phase-independent parameter (reddening?) in the IRSB analysis (see Sect. 2.1).

4.4 Properties of the identified relations

It is not a surprise that the majority of the flux ratios involves features in the blue region of the spectra (Fig. 10, top panel), given the wealth of potentially temperature-sensitive metallic lines in this spectral domain. However there is still a good number of relations in the [5000–6000] Å range, and more specifically around the magnesium triplet (MgT) and the Hα line. A handful of ratios will certainly prove very useful in the CaT region.

There are only three relations overlapping the spectral range of the Gaia Radial Velocity Spectrometer (RVS, Sartoretti et al. 2018). As shown in the second and third panels of Fig. 10, their precision and accuracy are among the lowest within our set of relations, probably because only a small fraction of the spectra in our sample cover this wavelength range. Such a paucity would also prevent a proper outlier rejection in case of undesirable artifacts in the analyzed spectra. At this stage, we would therefore caution against using these three relations solely, for instance in order to analyze Gaia RVS spectra. However, training at the resolution of the RVS (R = 11 500) might provide more calibration relations with better accuracy and precision. This will be investigated in a forthcoming paper.

In Fig. 11, we display a few interesting spectral regions harboring flux ratios:

• (a)

  The first feature we show is located in the blue spectral region, in the vicinity of Hβ. FRs involving the wings of hydrogen lines are typical T_eff indicators as already found by H18. A T_eff-sensitive feature (the blend of a Ni I and a Fe I line) is normalized by another spectral domain in the blue wing of Hβ, which varies much less with T_eff.

• (b)

  The second example is located in the MgT region (5167–5184 Å). One of the spectral domains is located directly in the bluemost of the Mg I lines of the Mgb triplet and varies notably with T_eff while the other one is only slightly temperature-sensitive.

• (c)

  Another example involves metallic lines, where the blend of a Fe I line and a Ti I line shows a strong dependence on T_eff, which seems to be compensated by the wing of a Fe I line whose sensitivity on T_eff is much less pronounced and might also correct for the effect of another parameter like log g or [Fe/H].

• (d)

  Finally, we show a flux ratio with two components located close to a H I line in the near-infrared. It does not have a typical behavior, as both components have an opposite sensitivity to T_eff : the first one is located in the wing of the H I line, which becomes stronger when T_eff increases, while the other one lies in the wing of a Fe I line, which becomes more prominent when T_eff decreases.

Before examining the second panel of Fig. 10, we wish to recall the argument already made at the end of Sect. 3.3: what we call precision is the scatter measured from test spectra of many different Cepheids. It combines two effects, of which star-to-star systematics is the dominant one, while the internal precision of an individual relation as deduced from a single target falls below 20 K. Keeping this in mind, we see that the precision on individual relations, capped at 100 K by construction, is actually much better (Fig. 10, second panel). Indeed roughly half of the relations have an individual precision better than 80 K, and the best relations have precisions ranging from 40 to 60 K; they cover, however, a more restricted wavelength domain.

When combining all relations, we are able to trace T_eff variations of a few Kelvins. For instance, the middle bottom panel of Fig. 9 indicates the outcome of our analysis when investigating consecutive short exposure spectra of l Car (P ~ 35.5 d) taken with HARPS within less than 30 min. We achieved a rms scatter of 1.6 K (using only the 100 best relations), which underlines the internal precision of our method.

Although our criterion was slightly looser (±40 K), the accuracy of our relations mostly falls within ±10 K (Fig. 10, third panel) and rarely exceeds 20 K, thus ensuring that we reproduce almost exactly the T_eff scale of the input labels. From a practical point of view, we would recommend to make use of only a selected number of the best relations (since their individual accuracy and precision are known) rather than sigma-clipping temperatures derived from individual relations, if the goal of the operation is to improve the precision of the final T_eff. This will probably come at the cost of the wavelength coverage, as shown in the second and third panel of Fig. 10. An alternative, possibly better, approach would be to perform a weighted median using the (known) inverse variance (i.e., squared precision) of individual relations as weight.

The last two panels of Fig. 10 display the excursion ΔFR_{min, max} covered by flux ratios, and the wavelength spacing Δλ of the two components of the flux ratios, respectively. They both can be interpreted by the larger occurrence of temperature-sensitive features in the blue regions of Cepheids’ spectra.

Fig. 10 Distribution of the properties of FRs as a function of wavelength. Blue, green, and red vertical spans depict the three spectral windows of 4MOST (filled) and WEAVE (hatched) high-resolution spectrographs. Upper panel: number of independent relations to derive Teff from flux ratios per wavelength range. Second panel: precision of individual relations to derive Teff from flux ratios as a function of wavelength. Third panel: accuracy of individual relations to derive Teff from flux ratios as a function of wavelength. Fourth panel: flux ratio excursion, ΔFRmin, max, as a functionof wavelength. Bottom panel: wavelength spacing Δλ of the two components of a ratio as a function of wavelength.

Fig. 11 Illustration of a few spectral domains where Teff-sensitive flux ratios have been recovered: a) blue wing of Hβ, b) MgT region, c) metallic lines around 5510 Å, d) H I line in the near-infrared. See Sect. 4.4 for details. Spectra are color-coded by temperature, where darker colors indicate cooler temperatures. For illustration purposes, the spectra were pseudo-normalized by the 99th percentile of the flux in the presented range.

Fig. 12 Same as Fig. 6, but where Teff derived for low-metallicity SMC Cepheids are shown by red crosses.

5 Discussion

5.1 Applicability at low metallicity

We used thesmall number of spectra of SMC Cepheids at our disposal to test the applicability of our method at low metallicity. Since the spectra do not cover the entire wavelength domain, not all the relations could be investigated. For those that could be tested, we see, however, a clear pattern exemplified in Fig. 12: the relations hold down to [Fe/H] ~−0.6 dex and the residuals quickly increase below this value, reaching ≈800 K at [Fe/H] ~−0.8 dex.

This is not a surprise, as the training sample does not contain any low-metallicity SMC Cepheids (see Sect. 2.5). In its current stage, the method is then applicable to Galactic Cepheids in the outer disk and to a large fraction of the LMC Cepheids, but not to the bulk of the SMC ones. Here we want to stress that the current values of [Fe/H] have been obtained using a T_eff scale based on line depth ratios, and could therefore be modified when adopting temperatures derived from our flux ratios.

Finally, we note that the lack of SMC Cepheid spectra not only deprives us from a training sample at low metallicity, but also decreases the quality of the spectral phasing: with only a few consecutive spectra at hands, we cannot investigate possible period changes or cycle-to-cycle variations, which are to be expected since the current sample is biased toward brighter, longer-period Cepheids.

Extending our T_eff scale to low-metallicity Cepheids would require a large number of spectra of SMC Cepheids. Including such spectra in the training will certainly lead to a different (presumably smaller) set of (more) universal relations. A larger number of relations may however be recovered by using higher-order/more-dimensional models.

5.2 A potentially new (spectroscopic) method for determining the reddening of Cepheids

Different parameters can limit the efficiency of the various flavors of IRSB methods in deriving accurate distances for Cepheids, namely the limb darkening correction (see Sabbey et al. 1995; Neilson & Lester 2013), the possible dependence of the projection factor on period (see Groenewegen 2013 and Kervella et al. 2017 for instance) or metallicity (discarded by e.g., Nardetto et al. 2011), the presence of a circumstellar envelope (e.g., Kervella et al. 2006; Hocdé et al. 2020), or the reddening.

If one assumes that reddening is the dominant uncertainty in the determination of T_eff as a by-product of IRSB methods, those Cepheids with minimal reddening (which we used in our training and testing sample) provide the absolutezero-point of the temperature scale. The other Cepheids should then align on the same scale. As a consequence, any systematic offset between the temperature derived (with unprecedented precision) using flux ratios and the temperature determined using the IRSB method could be attributed to uncertainties on the reddening. Conversely, it is possible to compute such a temperature offset and, in turn, to find out how the reddening would be modified for both T_eff scales to match. For the sake of the argument, we conducted such an exercise for the Cepheid δ Cep and managed to compensate for a systematic shift of 23 K between T_{eff, IRSB} and T_{eff, FR} by lowering E(B − V) by −0.006 mag (0.069 instead of the initial 0.075 mag). However, modifying E(B − V) does not result in a simple, global shift of T_{eff, IRSB}: Fig. 13 shows the outcome of our δ Cep test case. Itindicates that the difference in the IRSB temperatures computed either with E(B − V) = 0.075 mag or with E(B − V) = 0.069 mag is phase-dependent. Incidentally, since this phase dependence is different from the one seen in Fig. 9, and since its magnitude is different by a factor of ≈20, it cannot smooth out the current residuals. It is beyond the scope of this paper to investigate this behavior further. These findings may lead to a better understanding of the current scatter of IRSB-based PL relations and/or the p-factor controversy.

Fig. 13 Phase-dependent shift (modified − adopted) of the IRSB temperature of δ Cep when computed either with the modified E(B − V) = 0.069 mag or with E(B − V) = 0.075 mag as adopted throughout this work.

6 Conclusions

Using an innovative, data-driven method, we inferred a T_eff scale for classical Cepheids tying flux ratios to temperatures based on IRSB labels. For the first time in such studies, special emphasis was placed on accurately phasing the spectra when matching them to labels.

The next step to consider is to extend the T_eff scale toward lower, SMC-like, metallicities, which could not be done here simply because of the lack of spectra to build a training/testing sample. We also want to recover other atmospheric parameters like log g or [Fe/H] using a similar technique. On a more technical side, one could envision a more refined weighting scheme, involving all four quantities: precision, accuracy, ΔFR, Δλ.

Since it is purely spectroscopic, the method is independent from reddening to the same extent as the IRSB labels are, and could actually offer a new path to estimate reddening values for Cepheids. The derived temperatures are precise down to a few K, and we have full control of the accuracy, We estimate that the accuracy remains better than 150 K at all phases, a value which will be greatly reduced once the uncertainties on the reddening scale, hence the T_eff /E(B − V) degeneracy will be lifted. The uncertainties on E(B − V) obtained via the SPIPS algorithm are on the order of 0.025 mag (Breitfelder et al. 2016). Their use may then provide a different set of labels with reduced systematic uncertainties. On the other hand, since the number of Cepheids currently within reach of interferometric measurements is very limited, this would come at the expense of number statistics, and in particular period coverage. We desire not to train an inhomogeneous sample with different sources for the labels. Rescaling relations trained on a given set of labels would be the safer way to pursue, in case there is a constant or linear shift between IRSB and interferometric labels.

Our method to derive the temperature of classical Cepheids is very flexible as it works even at low S/N, and without requiring a normalization stage. Since it is always possible to down-sample high-resolution spectra to R = 5000 (at which the method is trained), we ensure homogeneous temperatures for both low- and high-resolution spectra, a very interesting feature with respect to large spectroscopic surveys. For instance, the stellar, circumstellar, interstellar component (SCIP) survey with WEAVE, the Milky Way Disc and Bulge Low-Resolution (4MIDABLE-LR, Chiappini et al. 2019), the Milky Way Disc and Bulge High-Resolution (4MIDABLE-HR, Bensby et al. 2019) and the One Thousand and One Magellanic Fields (1001MC, Cioni et al. 2019) 4MOST surveys will all observe classical Cepheids at R = 6500 and/or R = 20 000.

Our novel method paves the way for highly accurate and precise metallicity estimates, which will allow us to investigate the possible metallicity dependence of period-luminosity relations (e.g., Romaniello et al. 2008; Bono et al. 2010; Storm et al. 2011b) and ultimately solidify our measurement of the Hubble constant H₀.

Acknowledgements

We thank the referee, Dr. Simon Borgniet, for his useful comments which helped improving the quality of this paper. This work is partially based on observations collected at the European Southern Observatory. This work is partially based on data obtained with the STELLA robotic telescopes in Tenerife, an AIP facility jointly operated by AIP and IAC. This research hasmade use of the services of the ESO Science Archive Facility. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subprojects A05,A08).

References

All Tables

All Figures

Fig. 1 Impact of erroneous phasing of individual spectra on the Teff labels for the Cepheid RS Pup: top panel: MJD of the observation for individual spectra. Second panel: our individual radial velocity measurements (circles that are color-coded by MJD, cf., upper panel) versus radial velocity data used for the IRSB analysis (blue crosses, Storm et al. 2004). Spectra rephased to match the radial velocity curve are shown as pale brown squares. Using the same color coding, the third and fourth panels display the Teff variation, and the Teff offset between original labels and their corrected value after rephasing, respectively.

In the text

	Fig. 2 Parameter space (Teff, log g, [Fe/H]) covered by the compiled sample of spectra. Highlighted by salmon and red circles are the Cepheids ζ Gem and l Car, for which we have the most spectra at hand (135 and 191, respectively). Gray squares indicate a small portion of the giants used in the ATHOS dataset for stable stars (H18).
In the text

	Fig. 3 Temperature distribution of the spectra in our dataset, shown as a gray histogram. The temperature distributions of the training and testing samples are shown in blue and pale brown, respectively. We also indicated in red the temperature distribution of the spectra of stars with |Fe/H] < −0.3 dex that did notenter the training.
In the text

Fig. 4 Properties of an individual relation using flux ratios around 4405 Å. Top left panel: evolution of the neighboring spectral domain with Teff, where darker colors indicate deeper features. Bottom left panel: coolest (dark brown) and hottest (orange) spectra at our disposal. For illustration purposes, the spectra were pseudo-normalized by the 99th percentile of the flux in the presented range. In both panels, the regions considered for computing the flux ratios are marked in clear blue. Top right panel: analytical relation (dashed black curve) together with the Teff labels for the training (blue) and testing (pale brown) samples. Using the same color coding, bottom right panel: difference between the temperatures computed using the individual analytical relation and the input labels. In both panels, the coolest and hottest spectra are again depicted in dark brown and orange, respectively.

In the text

	Fig. 5 Precision and accuracy (in K) for all the relations retrieved by the analysis, color-coded by their correlationcoefficient. Only relations for which the precision is better than 100 K and less than 3 × the accuracy are further considered (red triangle). The dashed line shows a slight accuracy bias (see text) as a function of precision, which was computed using a boxcar median filter.
In the text

	Fig. 6 Residuals Teff, IRSB − Teff, FR as a functionof phase (ϕ, top left), log g (top right), [Fe/H] (bottom left), and log P (bottom right), for the same relation as in Fig. 4. Colored open circles highlight spectra for the targets ζ Gem, l Car, and δ Cep.
In the text

Fig. 7 Stability against S/N per Å for three relations of different maximum FR excursion (cf., color coding and legend in the lower panel). The corresponding S/N per resolution element for four different resolving powers (see legend) are given in the top panel. Shown are the measured mean values (colored dashed lines) and 1 σ standard deviations (solid lines enclosing filled areas) from an MC simulation with artificially injected noise (see main text for details). The zoom-in view emphasizes the (marginal) relation-to-relation inaccuracies (see middle panel of Fig. 10). Vertical lines indicate the minimum S/N requirements for different 4MOST consortium surveys: S∕N = 50 Å−1 in the magnesium or calcium triplet regions (4MIDABLE-LR10); S∕N = 30 Å−1 in the magnesium or calcium triplet regions (1001MC11, LR); S∕N = 100 Å−1 at 6200 Å (1001MC, HR; 4MIDABLE-HR12).

In the text

Fig. 8 Impact of continuum variations on FRs from close- (red) and wide-spaced (blue) wavelength components. Upper panel: 15 times the same ζ Gem spectrum as was used in Sect. 4.1 with varying continuum slopes (see color bar in lower panel). The missing part in the presented spectrum corresponds to HARPS’ chip gap. Middle panel: zoomed-in view of the region around ~ 5295 Å where the λ regions are highlighted by vertical spans. Lower panel: deduced Teff, FR values as a function of the measured FRs at different continuum slopes. The models are represented by red and blue solid curves while dashed red and blue lines and filled circles of the same colors indicate the finding from the original spectrum without any modification of the continuum.

In the text

Fig. 9 Representative examples of the three distinct residual patterns as discussed in the text: δ Cep (IRSB issuesseemingly constant with MJD, left panel); l Car (IRSB issueswith MJD-dependence, middle panel); ζ Gem (constant offset possibly due to erroneous reddening, right panel). Shown are the median residuals of the best 100 trained relations that were measured in individual spectra with respect to the input labels. The abscissa indicates the pulsation phase, ϕ. As in Fig. 1, the color coding resembles the time of observation of the individual spectra (cf., upper panels). Typical values for the 1 σ scatter are depicted by gray error bars. The inlay in the middle panel highlights five spectra of l Car that were taken with HARPS back to back within less than 30 min. The rms scatter of 1.6 K underlines the internal precision of our method.

In the text

Fig. 10 Distribution of the properties of FRs as a function of wavelength. Blue, green, and red vertical spans depict the three spectral windows of 4MOST (filled) and WEAVE (hatched) high-resolution spectrographs. Upper panel: number of independent relations to derive Teff from flux ratios per wavelength range. Second panel: precision of individual relations to derive Teff from flux ratios as a function of wavelength. Third panel: accuracy of individual relations to derive Teff from flux ratios as a function of wavelength. Fourth panel: flux ratio excursion, ΔFRmin, max, as a functionof wavelength. Bottom panel: wavelength spacing Δλ of the two components of a ratio as a function of wavelength.

In the text

Fig. 11 Illustration of a few spectral domains where Teff-sensitive flux ratios have been recovered: a) blue wing of Hβ, b) MgT region, c) metallic lines around 5510 Å, d) H I line in the near-infrared. See Sect. 4.4 for details. Spectra are color-coded by temperature, where darker colors indicate cooler temperatures. For illustration purposes, the spectra were pseudo-normalized by the 99th percentile of the flux in the presented range.

In the text

	Fig. 12 Same as Fig. 6, but where Teff derived for low-metallicity SMC Cepheids are shown by red crosses.
In the text

	Fig. 13 Phase-dependent shift (modified − adopted) of the IRSB temperature of δ Cep when computed either with the modified E(B − V) = 0.069 mag or with E(B − V) = 0.075 mag as adopted throughout this work.
In the text