Effective temperatures of red supergiants estimated from line-depth ratios of iron lines in the YJ bands, 0.97--1.32 micron

Determining the effective temperatures (Teff) of red supergiants (RSGs) observationally is important in many fields of stellar physics and galactic astronomy, yet some significant difficulties remain due to model uncertainty originating majorly in the extended atmosphere of RSGs. Here we propose the line-depth ratio (LDR) method in which we use only Fe I lines. As opposed to the conventional LDR method with lines of multiple species involved, the LDR of this kind is insensitive to the surface gravity effects and expected to circumvent the uncertainty originating in the upper atmosphere of RSGs. Therefore, the LDR--Teff relations that we calibrated empirically with red giants may be directly applied to RSGs, though various differences, e.g., caused by the three-dimensional non-LTE effects, between the two groups of objects need to be kept in mind. Using the near-infrared YJ-band spectra of nine well-known solar-metal red giants observed with the WINERED high-resolution spectrograph, we selected 12 pairs of Fe I lines least contaminated with other lines. Applying their LDR--Teff relations to ten nearby RSGs, the resultant Teff with the internal precision of 30--70 K shows good agreement with previous observational results assuming one-dimensional LTE and with Geneva's stellar evolution model. We found no evidence of significant systematic bias caused by various differences, including those in the size of the non-LTE effects, between red giants and RSGs except for one line pair which we rejected because the non-LTE effects may be as large as ~250 K. Nevertheless, it is difficult to evaluate the systematic bias, and further study is required, e.g., with including the three-dimensional non-LTE calculations of all the lines involved.


Effective temperatures of red supergiants
The effective temperature (T eff ) of the red supergiant (RSG) is one of the essential parameters in many fields of stellar physics, e.g., the theory of stellar evolution (e.g. Massey & Olsen 2003;Ekström et al. 2012;Choi et al. 2016) and theory of the maximum initial mass of the type-II supernova pro-genitor (e.g. Fraser et al. 2011;Smartt 2015). Moreover, the accuracy of T eff directly affects that of chemical abundances of RSGs, which could trace the abundance distribution of young stars in the Milky Way and nearby galaxies, leading to the galactic evolution theory (Patrick et al. 2017;Alonso-Santiago et al. 2019;Origlia et al. 2019). Hence importance of accurately determining T eff of RSGs with observations cannot be overstated.
However, the accuracy of the reported T eff of RSGs is still under debate. Interferometry and lunar occultation are the two simplest ways to measure T eff of nearby stars, but these methods are subject to uncertainties in some or all of the following three issues: (1) limb darkening (Dyck & Nordgren 2002;Chiavassa et al. 2009), (2) the complex extended envelope called MOLsphere (Tsuji 2006;Montargès et al. 2014) and (3) interstellar reddening Walmswell & Eldridge 2012). A more recent and improved way is to measure the stellar radius at continuum wavelengths to circumvent the effects of the MOLsphere with the spectro-interferometric technique (e.g. Ohnaka et al. 2013;Arroyo-Torres et al. 2013); however, the result is still affected by limb darkening and possibly uncertainty in interstellar reddening unless it is negligible.
Another type of method to estimate T eff of the RSG is to model-fit molecular bands in relatively low-resolution spectra in the optical (e.g. Scargle & Strecker 1979;Oestreicher & Schmidt-Kaler 1998). The reliability of this approach has been improved thanks to the advancements of the model of the stellar atmosphere such as, most notably, the MARCS model with improved handling of molecular blanketing (Gustafsson et al. 2008). The most comprehensive work of this approach was presented by Levesque et al. (2005) for RSGs in the Milky Way, in which they made use of the optical TiO bands (see also Levesque et al. 2006;Massey et al. 2009;Massey & Evans 2016, where RSGs in other galaxies were studied with the same method); their method is hereafter referred to as the TiO method.
A similar approach is to compare photometric colours of RSGs with the prediction of the MARCS model (e.g. Levesque et al. 2005;Drout et al. 2012;Neugent et al. 2012), in which molecular lines are taken into account to calculate the spectral energy distribution (SED) of the RSG. This approach tends to give a lower T eff in literature than the TiO method by up to ∼ 200 K (Levesque et al. , 2006, and hence there remain systematic uncertainties in fitting molecular absorption in optical spectra (Davies et al. 2013).
Alternatively some authors estimated T eff and spectral types of RSGs from molecular CO, H2O, OH and/or CN lines in near-infrared spectra (e.g. Carr et al. 2000;Cunha et al. 2007;Origlia et al. 2019). However, Lançon et al. (2007) and Lançon & Hauschildt (2010) found that the strengths of the CN molecular bands and the ratios of the strengths of the CO bandheads could not be reproduced with modern models of the static photosphere for the CNO abundances as predicted by stellar evolution models and a microturbulent velocity (vmicro) of 2 km s −1 .
The temperature estimated on the basis of molecular bands are affected more or less by a few factors, including chemical abundances and discrepancy between the atmosphere of the real star and that estimated with a simplified static model assuming one-dimensional local thermodynamic equilibrium (LTE) without MOLsphere. For example, Chiavassa et al. (2011b) and Davies et al. (2013) found that the TiO-band strength depends not only on T eff and the luminosity but also on the temperature structure, which is significantly affected by granulation.
A promising way to circumvent the problems raised so far is to use signals like weak atomic lines from the photosphere defined by the optical continuum or its outside close vicinity, given that they are expected to be less affected by some of the above-mentioned factors (Davies et al. 2013;Tabernero et al. 2018). A few works have presented T eff of RSGs estimated on the basis of photospheric signals, e.g., the SED of wavelengths less contaminated by molecular lines (the SED method; Davies et al. 2013), equivalent widths of some strong atomic lines in the J band (the Jband technique; Davies et al. 2015) and some atomic lines around the Calcium triplet (Tabernero et al. 2018). In these three works T eff of common RSGs in the Magellanic Clouds were derived. Nevertheless some systematic offsets are apparent between their results, possibly due to chemical abundances, contamination from some molecular lines to the signals and/or some effect of strong lines which is in part formed in upper layers of the atmosphere. Moreover, the departure from the 1D LTE condition may be important in RSGs (e.g. Bergemann et al. 2012b;Kravchenko et al. 2019). Most of the above works assumed the 1D LTE, and this shortcoming could contribute to the offsets.
Considering these discrepancies and the aforementioned factors that possibly add some systematic error to the temperatures estimated in most of the past methods, we here take a different method that satisfies the following three conditions. First, the method should use only relatively weak atomic lines that are not severely contaminated by molecular lines. Such lines do not originate in atmospheric layers far above the photosphere, where the temperature structure is not well constrained and moreover is expected to show variability. Second, the method should be independent of stellar parameters, abundances and reddening of both interstellar and circumstellar origins as much as possible. Third, ideally, the method should be independent of uncertainties in theoretical models, in particular with regard to the factors originating from the parameters in the upper atmosphere layers. Therefore, we use the 'line-depth ratio' (LDR) method, which relies on the empirical calibration of the relations between the LDR and T eff and satisfies the above-mentioned three requirements. In the next subsection, we explain what the LDR and LDR-method are and review relevant past studies.

Line-depth ratio as an indicator of the effective temperature
In cool stars (with a temperature of the solar temperature or lower), the depths of low-excitation lines of neutral atoms tend to be sensitive to T eff whereas those of high-excitation lines are relatively insensitive (Gray 2008a). Therefore, the ratios of the depths between the low-and high-excitation lines (that is, the LDRs) are good temperature indicators (Gray & Johanson 1991;Teixeira et al. 2016;López-Valdivia et al. 2019, and references therein). The basic procedure is as follows: (1) search high-resolution spectra for the line pairs whose depth ratios are sensitive to T eff , (2) establish the empirical relations between LDRs and T eff for stars with well-determined T eff and (3) apply the relations to target objects to determine their T eff . We focus on the near-infrared Y and J bands. These bands host a relatively small number of molecular lines in the RSG spectra (Coelho et al. 2005;Davies et al. 2010), and we can find sufficiently many atomic lines that are not contaminated by molecular lines (see Section 3.1 for detail and real examples). Taniguchi et al. (2018, hereafter Paper I) recently investigated LDRs in the Y and J bands of ten red giants with well-determined effective temperatures (3700 < T eff < 5400 K). They calibrated 81 LDR-T eff relations with the neutral atomic lines of various elements, and found that a precision of ±10 K is achievable in the best cases, i.e., high-resolution spectra of early M-type red giants with a good signal-to-noise ratio (S/N). This precision rivals those achieved in the previous results in which the LDRs in optical high-resolution spectra were employed (e.g. Kovtyukh et al. 2006).
Although Paper I gave well-defined LDR-T eff relations for solar-metal red giants, some significant systematic uncertainty remains in their work, where a simple relation between the LDR and T eff was assumed to hold universally. In reality, LDRs depend on other stellar parameters than T eff , even though T eff is the dominant parameter that determines the LDR of a neutral atom. Jian et al. (2019), for example, demonstrated that the LDR-T eff relations in the H band depend on the metallicity ([Fe/H]) by 100-800 K/dex. Also, Paper I suggested some effect of [Fe/H] and abundance ratios on their LDR-T eff relations. In this work, we limit the sample to objects with [Fe/H] around solar (Section 2.1) to circumvent the problem of the [Fe/H] dependence.
Surface gravity effect on LDRs has been also detected. Jian et al. (2020, hereafter Paper II) compared the LDRs of the line pairs, reported by Paper I, between 20 dwarfs, 25 giants and 18 supergiants and reported the LDRs' dependency on the surface gravity. This dependency is explained with difference in the ionization stages among the elements involved in the line pairs. All the past attempts to derive a number of relations utilized any combination of usable species among selected elements; accordingly, their estimates might be significantly affected by the surface gravity. The LDRs of two lines with common species are, in contrast, insensitive to the surface gravity. This kind of LDRs is also insensitive to the chemical abundance ratios. Therefore, the LDR-T eff relations derived on the basis of line pairs of the same species which are calibrated with the spectra of red giants are, when applied to RSGs, expected to yield T eff without introducing large systematic errors. Specifically, we use Fe i lines in this work because they are the only species that gives a sufficient number of the lines, hence a sufficient number of the line pairs, with which well-constrained LDR-T eff relations can be derived. Nevertheless, the difference in the size of the 3D non-LTE effect between red giants and RSGs could introduce systematic errors to T eff of RSGs obtained with our method, and should be examined.
In this work, using the high-resolution spectra in the near-infrared Y and J bands of solar-metal red giants with well-determined stellar parameters observed with the WINERED spectrograph (Section 2), we construct a set of many empirical LDR-T eff relations with only Fe i lines for the first time (Section 3). Then, we apply these relations to the WINERED spectra of nearby RSGs and deter-mine their T eff in an unprecedented accuracy as non-spectrointerferometric measurements (Section 4).

Sample
We consider two groups of targets for this work. The first group consists of nine well-known red giants and are used for calibration of the LDR-T eff relations in this work. Five of these stars (ε Leo, Pollux, µ Leo, Aldebaran and α Cet) were selected from Gaia benchmark stars (Blanco-Cuaresma et al. 2014), whose T eff were well constrained with interferometry by Heiter et al. (2015). The other four stars are well-studied nearby red giants selected from the MILES sample (Sánchez-Blázquez et al. 2006). The temperatures of these stars were determined using the ULySS program (Koleva et al. 2009) by Prugniel et al. (2011) with the MILES empirical spectral library as a reference. The temperature scale of Prugniel et al. (2011) relies on the compilation of literature stellar parameters, for which most spectroscopic works assumed the 1D LTE condition, and we have checked the consistency between the temperature scale of Prugniel et al. (2011) and that of Heiter et al. (2015) as follows. There are 13 stars 1 with [Fe/H] higher than −1.0 dex that were investigated in both papers. The temperatures from the two studies differ by 28 K on average with the unbiased standard deviation of the differences being 68 K, which is comparable with the combined measurement errors. Thus, the temperature scale of Prugniel et al. (2011) is the same as that of Heiter et al. (2015) within the errors, and we make use of the temperatures by Prugniel et al. (2011). The stellar parameters of these nine red giants are summarized in Table 1. We also adopted elemental abundances from Jofré et al. (2015) if available. We note that these red giants were also used in Paper I to establish the LDR-T eff relations from pairs of lines with all possible combinations of the selected elements. The parameters T eff of the red giants range from 3700 K to 5400 K and are used in the calibration of the LDR-T eff relations (Section 3.2). The other stellar parameters are used only to check blending of absorption lines and to select the useful lines (Section 3.1).

Observations
All the objects were observed using the near-infrared highresolution spectrograph WINERED installed on the Nasmyth platform of the 1.3 m Araki Telescope at Koyama Astronomical Observatory of Kyoto Sangyo University in Japan (Ikeda et al. 2016). We used the WINERED WIDE mode to collect spectra covering a wavelength range from 0.90 to 1.35 μm (z , Y and J bands) with a spectral resolution of R ∼ 28, 000. We selected the nodding pattern of A-B-B-A or O-S-O. All our targets are bright (−3.0 ≤ J ≤ 3.0 mag), and the total integration time for each target within the slit ranged between 3-240 sec, with which a S/N per pixel of 100 or higher was achieved. Telluric standard stars (slow-rotating A0V stars in most cases; see Sameshima et al. 2018) were also observed, the spectra of which were used to subtract the telluric absorption. Table 2 summarizes the observation log.
As in Paper I, we utilized the echelle orders 57th-52nd (Y band; 0.97-1.09 µm) and 48th-43rd (J band; 1.15-1.32 µm) only among the orders 61st-42nd because stellar atomic lines in the unselected orders are severely contaminated by other lines. The orders 61st-59th, 50th-49th and 42nd are dominated by telluric lines in our spectra collected in Kyoto (see Figure 4 in Sameshima et al. 2018). The orders 58th and 51st are dominated by stellar CN molecular lines of the (1-0) bandhead and those of the (0-0) band of the CN A 2 Π-X 2 Σ + red system, respectively (e.g. Kurucz 2011;Brooke et al. 2014). These two orders are heavily contaminated also by telluric lines in the spectra taken in summer. The wavelength ranges of the individual orders of WINERED are given in Table 4.

Data reduction
The initial steps of the spectral reduction were performed with the WINERED data-reduction pipeline ( (Table 6). We have not identified the origins of the three stellar lines labelled with '?', among which the line at 10338.5Å is catalogued in the list of unidentified lines by Matsunaga et al. (2020). Table 3. Stellar parameters used to simulate spectra of red giants and RSGs; we use these spectra for various purposes (see text). We assumed the solar chemical abundance pattern (Asplund et al. 2009 al. in preparation). Some details of the pipeline are given in Section 3.1 in Paper I. Then, the telluric absorption lines were removed using the observed spectra of A0V stars after their intrinsic lines had been removed with the method described in Sameshima et al. (2018), except for the 55th-53rd orders (1.01 to 1.07 μm) of the objects taken in winter, in which almost no significant telluric lines were present. The radial velocities were measured and corrected by comparing the observed and synthesized spectra. Finally the continuum was re-normalized and the realistic S/N was estimated as we describe in the following subsections. An example of the reduced spectra of RSGs is presented in Figure 1.

Continuum Normalization
The nominal continuum normalization was automatically performed in the data-reduction pipeline, but the automatic normalization often yields unsatisfactory results, especially for line-rich stars like RSGs. Therefore, we further optimized the normalizations for our spectra in the following procedure. First, we prepared the target spectra with the telluric absorption lines subtracted and with the wavelength-scale adjusted to the one in the air at rest. Second, we selected a set of continuum regions for each group of the red giants and RSGs, where the spectra are not significantly affected by the stellar lines. For the nine red giants, the continuum regions were selected with visual inspection. In contrast, we searched for the continuum regions of the RSGs by choosing the wavelength ranges with the depths from unity smaller than 0.3% in the model spectrum with the stellar parameters of 'RSG3' in Table 3 (see Section 3 about the spectral synthesis). These steps left a few tens of evenly distributing continuum segments in each order. Finally, we normalized the continuum of each order of each target star with the interactive mode of the IRAF continuum task, in which we mainly used the selected continuum regions but sometimes further combining some continuum regions visually selected for each order of each star.

Signal-to-Noise Ratios of Telluric-Subtracted Spectra
We need to know the errors in line depth for calibrating the LDR-T eff relations and for estimating T eff . The depth errors are considered to originate in three sources: (1) statistical noises in the telluric-subtracted spectra, which can be estimated with a combination of S/N of the target and telluric standard spectra (200 per pixel or higher in most of the cases), (2) residuals in the telluric subtraction and (3) systematic errors caused by slightly wrong continuum placement. Though the errors introduced by the first source can be estimated using a certain method, e.g., the one by Fukue et al. (2015, see their Section 3.2), those by the subsequent sources cannot be estimated in a straightforward manner.
In order to estimate realistic errors in line depths including all the above-mentioned three sources, we considered the 'effective' S/N per pixel of the telluric-subtracted spectra as follows. For this, we used the ready-to-use target spectra with normalized continuum. First, we measured the deviations, δi, from unity of each pixel, i, within the continuum regions selected in Section 2.3.1. The pixels {i} that satisfy |δi| > 0.05 might not be in real continuum regions, and thus were excluded in the subsequent analysis. Of course, such a threshold cannot be used for low-S/N spectra, but our spectra are sufficiently high 'effective' S/N, ∼ 50 at least. We note that some of such pixels might be due to the absorption lines in the observed RSG spectra that were not predicted in the synthesized spectra, and others might be induced by the bad continuum normalization. The 'effective' S/N were then estimated according to, after three iterations of three-sigma clipping, where N pixel is the number of the used pixels. We note that the sigma-clipping may underestimate the error, but by only 2% at most, providing that δi follows the Gaussian distribution. The resultant S/N of the RSGs are summarized in Table 4, whereas those for the red giants are found in Table 2 in Paper I. The reciprocal of the S/N measured here are henceforth regarded as the error in line depth. We note that the 'effective' S/N of RSGs could be under-or overestimated to some extent because, for example, weak stellar absorption lines that are not recognized in the chosen continuum regions may exist in reality. Moreover, the wavelength ranges in which many absorption features exist may be affected by the normalization error more than indicated by the 'effective' S/N because it is harder to trace the continuum at these regions than the selected regions without absorption lines.

CALIBRATION OF THE TEMPERATURE SCALE USING RED GIANTS
In this section, we establish the key relations of this work, i.e., a set of the reliable empirical relations between the LDR and T eff , in the following strategy. First, we choose Fe i lines that are relatively free from blending by other lines. Then we find the best pairs whose LDRs show a well-defined correlation with T eff . For the spectral synthesis employed in this section, we developed a wrapper software in Python3 of the spectral synthesis code MOOG (Sneden 1973;Sneden et al. 2012) 2 . MOOG synthesizes spectra, assuming the 1D LTE with the plane-parallel geometry. We remark that this assumption Table 4. 'Effective' S/N per pixel of the reduced spectra of RSGs in each echelle order (57th to 52nd in the Y band and 48th to 43rd in J ) and the wavelength coverage, λ min < λ air < λmax. does not affect our final determination of T eff , which relies only on the empirical LDR-T eff relations, though the selection of the candidate lines could be affected to some extent.
The grid of the model atmospheres of RSGs and red giants was taken from the MARCS (Gustafsson et al. 2008) with the spherical geometry (M = 2M for most cases) assumed, and the model atmosphere at a set of finer-scale stellar parameters was obtained with linear interpolation. We note that although the geometry of the radiative transfer code (plane parallel) is different from that of the model atmospheres assumed here (spherical), the difference in their geometries is known to induce only a small effect in general on the synthesized spectra (Heiter & Eriksson 2006). As for the line list, we used the third release of the Vienna Atomic Line Database (VALD3; Ryabchikova et al. 2015) as the main source. The molecular lines in the Y and J bands contained in the VALD3 are limited to CH, OH, SiH, C2, CN and CO, among which only CN gives a significant absorption in the spectra of our target RSGs and red giants. In addition, we adopted the line list of FeH by B. Plez (private communication) 3 . As a result of the extensive examination, we found that lines of FeH appeared in the spectra of the RSGs with a depth up to ∼ 10% and that those lines were well reproduced by synthesized spectra with a dissociation energy of 1.59 eV (Schultz & Armentrout 1991).

Line selection
We chose the Fe i lines that are sufficiently strong and relatively free from blending out of the list of the Fe i lines in the VALD3. Since our aim is to estimate T eff of RSGs using the LDR-T eff relations calibrated with red giants in spite of the difference in the surface gravity between the two groups, it is necessary to choose the Fe i lines that are least sensitive to the surface gravity. The LDR of a pair of lines from two different elements, especially those with different ionization energies, is known to show dependence on the surface gravity (Paper II). Also, depths of molecular lines are significantly different between RSGs and red giants because of their strong dependence on the surface gravity (Lançon et al. 2007), non-solar CNO abundances (Lambert et al. 1984;Ekström et al. 2012) and existence of MOLspheres In order to evaluate the effect of the potential blending, we used synthetic spectra of the two model RSGs (RSG1 and RSG2 in Table 3) and those of the five coolest red giants (µ Leo, Alphard, Aldebaran, α Cet and δ Oph). The stellar parameters of RSG1 and RSG2 roughly correspond to those of the warmest and coolest RSGs in our sample, respectively, and the line blendings in RSG2 are expected be severer than those in all the sample observed RSGs. With each set of stellar parameters we measured four types of depths (dAll, dOneOut, dSameElIonOut and dAtomOut) for each Fe i line in VALD3 (let us use λ0 to denote the centre wavelength listed in VALD3) as follows. First, we synthesized the respective four types of spectra around the target line with different groups of the lines included: (1) All-all the atomic and molecular lines, (2) OneOut-all the lines except for the target Fe i line, (3) SameElIonOut-all the lines except any Fe i lines and (4) AtomOut-only the molecular lines (see an example in Figure 2). Next, we identified the absorption around λ0 in the synthesized spectrum All and determined the peak wavelength, λc, the value of which is different from λ0 if the target line is blended, or identical to it if not. Finally, the respective depths from unity at λc in all the four types of the synthesized spectra were measured. The reason why we measured the depths at λc rather than λ0 in the synthesized spectra is that we measured the depth of each line in the observed spectra at the wavelength of the peak position around the line rather than the wavelength in the line list λ0 (Section 3.2).
With these depths, we considered two types of criteria for selecting lines. First, we imposed dAll > 0.02 on the synthetic spectra of RSG1. When two or more lines were assigned to the same wavelength of minimum λc in RSG1, we would consider only the line that contributes most to the peak, i.e. the line with the smallest dOneOut, as a candidate.  Figure 2. Examples of the synthesized spectra of RSG3 around Fe i λ10742.550Å (λ 0 ) with different lines included. Cyan, orange, red and pink lines show the spectra of All, OneOut, SameElIonOut and AtomOut, respectively, defined in the main text. The peak wavelength in All (λc) was indicated with the vertical orange dashed line. We note that this line in this figure shows one of the most complex blended lines among the final set of the selected lines.
Then, we imposed the following three criteria, which should be satisfied in all the seven synthetic spectra of the RSGs and the red giants: (2) Applying the combination of all the above-mentioned criteria to the current line list left 76 Fe i lines in total (41 in the Y band and 35 in J ). We note that the index dOneOut/dAll had been used in evaluating line blending in some of our recent papers (Kondo et al. 2019;Matsunaga et al. 2020;Paper II); the difference in this work was that two slightly different indices, dSameElIonOut and dAtomOut, were added to the condition to give tighter constraints. Furthermore, we considered the following six points, with which some candidate lines would be filtered out. First, some of the observed lines were found to be contaminated by unexpected stellar lines. For example, though Fe i λ12393.067Å is expected to be blended only with Fe i λ12393.281Å in the synthetic spectrum and to meet the criteria described above, the peak wavelengths λc in the observed spectra of all the red giants were found to be systematically shorter by ∼ 0.25Å than the wavelength (λ0) listed in the VALD3 line list; this fact implies that this line is blended with an unknown line(s) with its peak wavelength at around 12392.7Å. Second, some of the observed lines were contaminated by residuals from telluric subtraction or suffered from imperfect continuum normalization due to many stellar lines, especially CN molecular lines, concentrated in the close vicinity of the line. Third, some of the observed lines were found to be much weaker than those in the synthetic ones possibly due to the inaccurate oscillator strength values in the VALD3 line list. For example, the depths of Fe i λ10462.155Å were 0.016-0.028 in the syn-thesized spectra of the four coolest red giants and this line meets the criteria, but the observed depths are too small, 0.010-0.019. Fourth, we excluded two lines, Fe i λ10780.694 and λ12340.481Å, because they were affected with diffuse interstellar bands λ10780 and λ12337, respectively, reported by Hamano et al. (2015). Fifth, hydrogen Paschen series and helium 10830Å lines are known to show large variations between stars, especially in supergiant stars (e.g. Obrien & Lambert 1986;Huang et al. 2012), and hence may affect the blending significantly. However, we found that none of the selected lines were affected by them in our case. Sixth, Matsunaga et al. (2020) detected dozens of lines that are not found in the available line lists in the YJ -band spectra of supergiants (4000 < T eff < 7200 K), and these lines may affect our result significantly. However, all the lines selected here were separated by more than 30 km s −1 from these uncatalogued lines (Table 3 of Matsunaga et al. 2020), and thus we safely ignored this potential influence.
Applying all these points to the candidate lines, we obtained in the end 52 Fe i lines in total (28 in the Y band and 24 in J ) and summarize the result in Table 5. We note that Paper I considered another criterion requiring that the depth of each line be smaller than 0.5 in the observed spectra of Arcturus; this condition is satisfied in all the lines that we selected.
Then, we measured the depths of the selected 52 Fe i lines in each observed spectrum by fitting a quadratic function to three (or four) pixels centred at the peak of each line (Strassmeier & Schordan 2000). Here we define d i , which includes both the systematic errors (e.g., in the telluric subtraction and the continuum normalization) and the statistical errors (Poisson, read-out and background noises).

Line-pair selection: method
We searched for the best set of line pairs of which the LDR-T eff relations yield the most precise estimates of T eff . Here, we followed the procedures in Paper I and Paper II except that we combined the lists of Fe i lines in different echelle orders in each of the Y and J bands. The basic idea of the process is to select the set of line pairs that meet the following conditions best: (1) high precision in reproducing T eff of our sample stars and (2) small difference in wavelength between the two lines of each line pair. The latter condition is desirable in the sense that different instruments from WINERED in future observations would have a better chance to be capable of detecting both lines of the line pair simultaneously.
First, we evaluated the relation between LDR and T eff of all the possible line pairs individually. In this case, 211 and 128 pairs in the Y and J bands, respectively, were found to have two lines that were both detected in more than four stars and have excitation potentials (EPs) separated by more than 1 eV. For each pair, denoted with the subscript j, of all the 211 + 128 line pairs, we calculated the LDR denoted as r (n) j and its error. We plotted T eff against the common logarithms of the LDRs (log rj) for each pair and determined the regression line T eff = aj log rj + bj, using the Weighted Total Least Squares method (see Markovsky & Van Huffel 2007, for a review). For this regression, the weight of each point was given by where ∆T for each star, respectively. The weighted dispersion with regard to each regression line was then calculated as (4) We set the threshold of σj < 150 K for the pair to be valid. As a result, 41 and 6 line pairs in the Y and J bands, respectively, were left in the subsequent analyses. We also set another threshold aj > 0, but no line pairs were rejected with this condition. Second, we searched for the optimal set of the line pairs for each of the Y and J bands as follows. Each set of line pairs can be considered as an undirected graph whose nodes correspond to absorption lines and whose edges connecting the nodes correspond to line pairs. We require that each line be used only once in each set, i.e. no duplication of a line in separate pairs are allowed. This ensures that the errors in r (n) j values in different LDR-T eff relation be independent of each other and thus makes it straightforward to calculate the statistical errors of the combined T eff (Equation 8 and Equation 9). We determined the optimal matching, M , that meets our requirements as follows. We considered the maximum matching, where the number of the edges of the matching is as large as possible under the condition where no nodes are connected with more than one edge. In the ideal case, the size of the maximum matching, |M |, corresponds to a half of the number of all the selected lines. However, it did not in our case because many edges were rejected due to, for example, large scatters around the corresponding LDR-T eff relations. For a given maximum matching, M k , of this undirected graph, we calculated the weighted mean temperature T (n) M k on the basis of the LDR-T eff relations for each star. We also considered the difference in wavelength between the two lines of each line pair j k,m , denoted as ∆λ m,k , and calculated the evaluation function E(M k ; e) defined as where ET (M k ) represents the size of the error in redetermining T eff of the nine stars for a given matching M k , and E λ (M k ) represents the wavelength difference of the line pairs and works as a penalty term. There are different allowed combinations of the line pairs which form different maximum matchings, and for a given e value we selected the one that gives the least E(M k ; e) as the optimal matching, M k(e) . We would determine the coefficient e by considering how ET (M k(e) ; e) and E λ (M k(e) ; e) depend on e, as explained in the following subsection.

Line-pair selection: results
We applied the procedure described in the previous subsection to our 28 and 24 nodes (i.e. lines) and 41 and 6 edges (i.e. line pairs) in the Y and J bands, respectively. We searched for the optimized matching, M k(e) , that gives the smallest E at each e value between 0 and 2 K/Å to see how e would affect the solutions. Figure 3 plots the values of ET and E λ for M k(e) with varying e. Larger the parameter e is, larger the weight on E λ is relative to ET in the total evaluation function E, and then the optimal matching and the values of ET and E λ change. We decided to employ the same value, e = 0.5 K/Å, as Paper I because ET at e = 0.5 K/Å is similar to that at e = 0 K/Å, which optimizes the precision in the redetermined temperature, whereas E λ is sufficiently small. Moreover, the same matching would have been selected if another value of e among a wide range of e values (0.39 < e < 1.81 K/Å) was employed.
Consequently, we obtained 10 and 2 line pairs in the Y and J bands, respectively, i.e., 12 line pairs in total (Table 6 and Figure 4). This set of line pairs made use of 24 unique lines in the Y and J bands combined (nb., no duplications of the lines are allowed; see Section 3.2). This condition has only a weak effect on the statistical error, which could have been better by only ∼ 10% even if we had used all the 41 + 6 pre-selected pairs. The number of the selected lines (52 in total, or 28 in Y and 24 in J ; see Section 3.1) is about onefourth of that used by Paper I (224 in total, or 125 in Y and 99 in J ), where atomic lines of various elements were used. Accordingly, the total number of the final line-pairs, 12, is smaller than 81 in Paper I. Nevertheless, the number 12 is larger than the number, 8, of the Fe i line pairs in Paper I.
The evaluation functions with e = 0.5 K/Å with all the line pairs in the Y and J bands combined are E all T = 43.9 K and E all λ = 289.3Å. The former E all T in this work is only 1.5 times larger than the counterpart in Paper I despite the smaller number of the line pairs and smaller depths, both of which could have increased the statistical errors in T (n) M k . The relatively small difference in E all T indicates that the errors in literature T eff (19-65 K) mainly contribute to E all T in Paper I. In contrast, the latter E all λ in this work is 4.3 times larger than that by Paper I and is similar to that by Fukue et al. (2015) for the H band. This is expected because Fukue et al. (2015) and we treated lines in each band together, whereas Paper I treated those in each echelle order independently. These differences in ET and E λ are well demonstrated in Figure 3.

Re-determination of the effective temperatures of red giants
With each line pair j, the effective temperature and its error of a target star (RSG in our case) would be estimated to be ∆Tj = (aj∆ log rj) 2 + Saa(log rj) 2 + 2S ab log rj + S bb , e>.Å@ where Saa S ab S ab S bb is the variance-covariance matrix of the coefficients (aj, bj) of the regression line. The value referred to as the LDR effective temperature TLDR of the target star is defined as the weighted mean of the temperatures estimated from the available line pairs among the 12 line pairs (Tj ± ∆Tj; for j = 1, · · · , Npair). We define two forms of the error for TLDR, given by the following two formulae, and calculated both of these for each target, Table 6. List of low-and high-excitation Fe i lines and the LDR-T eff relations; a and b denote the coefficients in T eff = a log r + b, the three S values (Saa, S ab and S bb ) represent the variance-covariance matrix of a and b, N is the number of stars used in the fitting, and σ is the weighted dispersion given by Equation 4. The line pair (6)  log r j T eff >.@  Table 6). The line-pair ID together with the wavelengths (Å) of low-and high-excitation lines are indicated at the top of each panel. Blue data points show the relation between the observed log(LDR) and effective temperatures in literature for red giants. Red lines indicate the best-fitting relations that we obtained, T eff = a log(LDR) + b, with light-red shaded areas for the 1σ confidence intervals.
The former (∆T Eq8 LDR ) is the weighted standard error of TLDR used by Paper I, and the latter (∆T Eq9 LDR ) is the error prop-agated from the errors in the LDRs and coefficients of the relations.
In order to examine the dependence of the relations on T eff and [Fe/H], we determined TLDR of the group of the nine red giants (Table 2). Table 7 summarizes the re-determined effective temperatures and their standard errors in this work, together with those estimated on the basis of the relations Table 7. Effective temperatures of nine red giants in kelvin. T eff is the literature effective temperatures adopted in this work. ∆T Eq9 LDR in group (T18) is calculated in this work, whereas T LDR and ∆T Eq8 LDR in group (T18) are adopted from Paper I, in which they calculated the values in the same way as in this work.

Paper I (T18)
This work (TW) in Paper I. In the table and hereafter, the former and latter cases are abbreviated as 'TW' (as of This Work) and 'T18', respectively, and the parameters based on either of them are distinguished with the corresponding abbreviated word in parentheses as the suffix; e.g., ∆T Eq8 LDR (TW) denotes the error of the LDR effective temperature on the basis of Equation 8 calculated with the relations in this work (TW). To estimate the standard errors according to Equation 8 is, however, problematic in the case of the relations in this work for two reasons; (1) the error of ∆T Eq8 LDR (TW) may be large because the number of the line pairs is small, and (2) the dispersion of Tj, and hence ∆T Eq8 LDR (TW), may be underestimated because the nine red giants themselves were used to calibrate the LDR-T eff relations. The errors according to Equation 9 are more likely to be closer to the true values as long as the errors of the LDRs are accurately estimated. In fact, ∆T Eq8 LDR (TW) tends to be significantly smaller than ∆T Eq9 LDR (TW) probably because of the aforementioned reasons. Therefore, we concluded that the latter, ∆T Eq9 LDR (TW), is more robust and reliable than ∆T Eq8 LDR (TW). In contrast, the error ∆T Eq8 LDR (T18) was found to be similar to ∆T Eq9 LDR (T18) probably because of the large number of line pairs, and thus we simply adopted ∆T Eq8 LDR (T18) as in Paper I. Figure 5 demonstrates the differences between the redetermined effective temperatures TLDR(TW) and those in literature T eff . It shows no apparent correlation between the difference and either of T eff and [Fe/H]. Our TLDR was consistent with that in literature within ∼ 50 K over the entire range of T eff and [Fe/H] of our sample red giants for the calibration. Besides, the error bars on the basis of ∆T Eq9 LDR (TW) and ∆T eff explain the scatter around zero well, indicating that the error estimates were reasonable.

LDR temperatures of red supergiants
Our new LDR-T eff relations for Fe i-Fe i line pairs are supposed to be less affected by the surface gravity and line broadening than the relations by Paper I as a result of our careful selection of the isolated Fe i lines. We here apply the relations directly to ten nearby RSGs in the same way as to Table 8. Effective temperatures of ten RSGs using different sets of LDR-T eff relations in kelvin. T TiO (L05) shows the effective temperatures estimated with the optical TiO bands by Levesque et al. (2005) for comparison.

Paper I (T18) This work (TW)
Object  red giants in Section 3.4. Table 8 lists the resultant TLDR of our targets together with the temperatures derived using the relations of Paper I.
The accuracy of our TLDR of RSGs strongly depends on the assumption that empirical LDR-T eff relations of Fe i lines are well consistent between giants and supergiants. Paper II concluded that, at T eff = 4500 and 5000 K, the LDRs of Fe i lines of dwarfs, giants and supergiants agree with each other. In order to test whether the agreement is also found between red giants and RSGs with 3500 T eff 4000 K, we first compare TLDR(TW) and TLDR(T18) to examine the dependence of the LDRs of multiple species on the surface gravity. The effect of the surface gravity on the LDRs is, at least to some extent, governed by the difference in the ionization energies, ∆χ, of the elements forming low-and high-excitation potential lines (Paper II). In fact, a weak anti-correlation of ∼ −80 K/eV between ∆χ and T eff from individual LDR relations (Tj) is seen for the ten RSGs (Figure 6). This trend is similar to the one observed by Paper II for giants and supergiants at higher temperatures. It indicates that the surface gravity effect on the LDRs remains similar down to the low-temperature ranges of the RSGs and, therefore, TLDR(T18) of the RSGs have systematic er-rors. Moreover, the residuals from the ∆χ-Tj correlation for most of the line pairs can be explained with the measurement error. Therefore, the combination of the ∆χ effect and the measurement error may have a larger impact on Tj estimation for most line pairs with multiple species than other effects (e.g. 3D non-LTE correction).
Concerning the error in TLDR, ∆T Eq9 LDR (T18) is found to be systematically (1.5-2 times) smaller than ∆T Eq8 LDR (T18); this fact suggests that Tj for T18, hence TLDR(T18), has systematic errors introduced by the effect of the surface gravity. In contrast, ∆T Eq9 LDR (TW) is similar to or larger than ∆T Eq8 LDR (TW). This fact implies that the scatter of Tj of Fe i line pairs with each star can be explained mostly with the expected statistical errors except for a few cases; three stars (TV Gem, BU Gem and ψ 1 Aur) have larger line broadening than the other RSGs, which could slightly bias their TLDR. However, the effect could be negligible because the lines that we used are well isolated ones. We also examine the error budget of TLDR(TW) using Tj − TLDR, which is a measure of the error in Tj (blue dots in Figure 7). Though Tj − TLDR with each Fe i line pair has a distribution with a non-zero offset, the offsets are within ∼ ±150 K for the pairs except the pair ID (2), which has an unpredictable large offset (∼ 500 K). These offsets are smaller than the measurement error, which indicates, again, that the error in TLDR(TW) is dominated by the expected statistical error. Moreover, the small systematic offsets may be partly explained by other lines contaminating the Fe i lines that we use.
Consequently, our new LDR-T eff relations with only iron lines (∆χ = 0 eV) are not affected by the surface gravity. Since ∆T Eq9 LDR may be slightly overestimated (see the discussion on the 'effective' S/N in Section 2.3.2), we adopt ∆T Eq8 LDR as the more accurate error of TLDR for the RSGs than ∆T Eq9 LDR . Note that since ∆T Eq9 LDR is calculated using the scatter of Tj, this error includes a part of the systematic error on Tj. The inclusion of the systematic error and the larger measurement error may explain the larger ∆T Eq9 LDR (TW) of RSGs than that of cool red giants.
In order to further examine the consistency between the LDR-T eff relations of red giants and RSGs, we compare the Fe i LDRs of the two groups by making use of their model spectra synthesized with MOOG (see Section 3 for spectral synthesis). We calculated the differences ∆ log rj in the logarithm of the LDRs log rj between RSG3 and RGB1, which have the same T eff of 3850 K (see Table 3 for their stellar parameters). Though ∆ log rj calculated with the imperfect model spectra cannot reproduce observed one in a quantitative sense, aj∆ log rj can be treated as a measure of the systematic error on Tj due to the difference in the LDR-T eff relations between red giants and RSGs. Figure 7 shows aj∆ log rj and Tj − TLDR for each line pair. Except the pair ID (2), aj∆ log rj calculated with three types of synthetic spectra and observed Tj − TLDR are similar to each other in a qualitative sense. Moreover, since aj∆ log rj calculated with MOOG has a nearly symmetric distribution around the zero, the systematic errors for individual line pairs are expected to be cancelled out by taking an average of Tj. Though the reason of the inconsistency of the pair ID (2) is unclear, it might be safe to reject this pair because this pair is suggested to involve a large ( 200 K) systematic error both observationally and theoretically. We, thus, regard the weighted mean of Tj after excluding the pair ID (2) as the final estimate of the LDR temperature of RSGs (Table 9).
Since some lines of red giants and RSGs suffer from the non-LTE effects in the extended atmospheres of these objects (e.g. Bergemann et al. 2012b;Lind et al. 2012), we finally examine how much the non-LTE effects could affect the LDR temperatures. For this purpose, we used the online spectral synthesis tool developed by M. Bergemann's group (Kovalev et al. 2019) 4 to synthesize model spectra with and without the non-LTE effect. This tool can account for the non-LTE effect for many Fe i lines calculated by Bergemann et al. (2012a). However, the line list of the tool contains only half of the Fe i lines that we have selected, and we could calculate the LDRs of only four of the final line pairs with the tool (IDs 2, 3, 4 and 7). The difference in aj∆ log rj between LTE and non-LTE model spectra gives an estimate of the systematic bias of the LDR temperature caused by the non-LTE effect: ∼ 200 K for ID (2), ∼ 250 K for ID (3), ∼ 0 K for ID (4) and ∼ 100 K for ID (7) (see Figure 7). While the observational results indicated by blue 4 Last accessed on 2020 June 1. /'5,' 6\VWHPDWLFELDVRQTj>.@ 022*/7( %HUJHPDQQ/7( %HUJHPDQQQRQ/7( 2EVHUYHGTj TLDR Figure 7. Blue dots show the empirical deviation between the temperatures from a given line pair, T j , and the LDR temperatures, T LDR (TW), for each Fe i line pair for each star. Orange squares, green triangles and red inverted triangles show theoreticallyestimated systematic bias (a j ∆ log r j ) of the LDR-based temperatures T j caused by the difference between red giants and RSGs calculated with MOOG, Bergemann's tool with LTE and that with non-LTE, respectively. dots in Figure 7 show consistent trends for the pairs ID (4) and (7), we found conflicting results for the pairs ID (2) and (3). The former two line pairs with the small non-LTE effect indicate that the systematic bias caused by the non-LTE effect is ∼ 100 K or less, though we cannot draw a strong conclusion based on the two pairs only. Moreover, it is hard to predict the systematic bias caused by the non-LTE effects without the non-LTE calculations for all the lines fully performed. We note that the 3D effect could also have an impact on line depths of late-type stars (e.g. Collet et al. 2007;Kučinskas et al. 2013), but the lack of a comprehensive grid of 3D model atmospheres of RSGs prevents us from estimating its impact.

Comparison with previous T eff determinations
Betelgeuse is one of the most studied RSGs, whose T eff is suggested to be free from strong time variability (White & Wing 1978;Bester et al. 1996;Gray 2008b;Levesque & Massey 2020), and thus it must be a good standard RSG for comparing T eff by different methods. Many previous works have estimated T eff of Betelgeuse using a variety of methods, e.g., broad-band interferometry (e.g. Dyck et al. 1998;Haubois et al. 2009), SED fitting (e.g. Tsuji 1976Scargle & Strecker 1979), the TiO method (e.g. Levesque et al. 2005;Levesque & Massey 2020) and the excitation balance of CO molecular lines (e.g. Lambert et al. 1984;Carr et al. 2000), as summarized in Dolan et al. (2016). These methods may, however, be affected by the systematic uncertainties described in Section 1.1. Davies et al. (2010) used the Jband technique, which they claim is less affected than some other methods by these systematic uncertainties, and obtained 3520 ± 160 and 3660 ± 170 K for spectra with spectral resolutions of R = 2, 000 and 6, 000, respectively. Our result of 3611 ± 38 K for Betelgeuse is consistent with theirs. Their Table 9. Final estimates of effective temperatures and luminosities of ten RSGs derived with 11 LDR-T eff relations except ID (2)  estimate for V424 Lac, 3580 ± 230 K, based on its spectrum with R = 2, 000 is also consistent with ours, 3749 ± 49 K, within the errors, but the variability of V424 Lac is poorly known. We note that T eff of no other RSGs in our sample were measured by Davies et al. (2010) and other recent works using the J -band technique.
Another promising approach is the spectrointerferometric observation, which is less affected by the outer envelopes of RSGs. Ohnaka et al. (2011) determined Betelgeuse's T eff to be 3690 ± 54 K with this approach. In addition Arroyo-Torres et al. (2013) estimated it to be 3620 ± 137 K, using the angular diameter measured by Ohnaka et al. (2011). Our result of 3611 ± 38 K for Betelgeuse is consistent with theirs, and hence supports the validity of the T eff scale that we obtained. Unfortunately, no other RSGs than Betelgeuse in our sample have the spectro-interferometric T eff in literature to compare with. Figure 8 compares the effective temperatures of all the RSGs in this work (TLDR ±∆T Eq8 LDR ) and those determined by Levesque et al. (2005) (TTiO(L05)), where the TiO method was employed. Whereas comparing T eff of each star is not so practical, given that the time variations of T eff for RSGs are as large as several hundreds kelvins in the worst cases (e.g. Levesque et al. 2007;Clark et al. 2010;Wasatonic et al. 2015), the comparison of these as a whole tells something significant; the two sets of T eff were found to be well consistent, considering the errors of the former (∆T Eq8 LDR ) and latter (∼ 50 K), but with a slope of 0.70 ± 0.14, which is slightly different from one. This consistency indicates that the TiO method yields not strongly biased T eff of RSGs that have the solar abundance pattern. However, it is unclear whether the TiO method can yield unbiased T eff of RSGs that have non-solar abundances, considering the uncertainties discussed in Section 1.1. In fact, Levesque et al. (2006) and Davies et al. (2013) determined T eff of common ∼ 20 RSGs in the Magellanic Clouds, using the TiO method, but T eff by Levesque et al. (2006) are systematically ∼ 100 K higher than those by Davies et al. (2013). Moreover, Davies et al. (2013) compared T eff derived with the TiO method to that with the SED method, and found that the latter method yielded ∼ 400 K higher T eff than the former.

Red supergiants on the HR diagram
In many previous works, more than a decade ago, observationally-determined T eff of the brightest RSGs (M bol < −7 mag) were often lower than expected from theoretical models at a given luminosity (see, e.g., Figure 8 in Massey 2003). Since then, several observational studies have compared their T eff with theoretical models, mainly with Geneva's stellar evolution model (Ekström et al. 2012;Georgy et al. 2013), on the Hertzsprung-Russel (HR) diagram. Some of them found good consistency between them (e.g. Levesque et al. 2005;Davies et al. 2013;Wittkowski et al. 2017), whereas some others found offsets between the measurements and theoretical models (e.g. Levesque et al. 2006;Tabernero et al. 2018).
In order to plot our observed RSGs on the HR diagram, we derived the bolometric luminosity, using the parallax and bolometric magnitude of each source as follows. As for the parallaxes, we used the values in Gaia EDR3 (Gaia Collaboration et al. 2016Collaboration et al. , 2021 for all our RSG samples except Betelgeuse, for which we used the parallax in the Hipparcos catalogue by van Leeuwen (2007) because Gaia EDR3 has no entry for it. Since there is a systematic bias in the Gaia parallax (Lindegren et al. 2021a), we followed the recipe by Lindegren et al. (2021b) to reduce the bias. In brief, we inputted G-band magnitude, effective wavenumber and ecliptic latitude tabulated in Gaia EDR3 into the code developed by them 5 to calculate the bias estimate. Estimated biases, summarized in Table 9 and ranging from −43 to −22 µas, only change the final log L by up to 0.07, which is smaller than the statistical error in it. Nevertheless, it should be kept in mind that (i) all the sample RSGs have G < 6 mag, which is out of the range of the bias calibration by Lindegren et al. (2021b), and (ii) the Gaia parallaxes of very bright stars could contain large systematic biases; for example, several stars with G 4 mag and with parallaxes larger than 5 mas have systematic biases of 1 mas in Gaia DR2 (Drimmel et al. 2019) though the bias would be smaller in EDR3 (Torra et al. 2021), and Gaia DR3 parallaxes of all our sample RSGs are smaller than 4 mas. We note that the statistical errors of the parallaxes of our targets tend to be larger than those expected from their magnitudes in both the Hipparcos (van Leeuwen 2007) and Gaia EDR3 (Lindegren et al. 2021a) catalogues. It is most likely, at least partly, due to strong granulation in RSGs, which fluctuates the positions of their photometric centroids by ∼ 0.1 AU (Chiavassa et al. 2011a;Pasquato et al. 2011).
As for the bolometric magnitude, we estimated it for each RSG sample from the 2MASS K s-band magnitude (Skrutskie et al. 2006), considering the extinction and the bolometric correction. The sampled epochs of the catalogued K s-band magnitude that we used are different from those of their WINERED observations. However, the difference in the epochs would not be a significant problem, given that the magnitude variation of the RSG in the infrared bands is known to be negligible (e.g. Yang et al. 2018;Ren et al. 2019), in contrast to the optical variation, which is as large as a few magnitudes (Kiss et al. 2006;Soraisam et al. 2018, and references therein). Concerning the interstellar and circumstellar extinction, we converted the V -band extinction A(V ) of the RSGs measured by Levesque et al. (2005), with the typical precision of 0.15 mag, to the K s-band extinction A(Ks), using the reddening law A(Ks)/A(V ) = 0.116 by Cardelli et al. (1989), assuming the total-to-selective extinction ratio RV = 3.1 as in Levesque et al. (2005). We estimated the bolometric correction for the K s band, BCK s , interpolating the BCK value with regard to the obtained T eff (Section 4.1) in the tabulated relation between T eff and BCK for RSGs in Table 6 of Levesque et al. (2005). Since the extinction to our sample RSGs are small, A(V ) = 2.17 mag at most, the choice of the reddening law and the RV value and the difference between K s and K do not significantly affect the resultant luminosities. Finally, the bolometric luminosity (L) of each target RSG and its error were calculated using the Monte Carlo method (Table 9).
In the HR diagram in Figure 9, we compared the dis-5 https://gitlab.com/icc-ub/public/gaiadr3_zeropoint T eff >.@ tribution of the data points of our sample RSGs in our estimates with those expected from the latest Geneva's stellar evolution model with the solar metallicity (Ekström et al. 2012). Although the sample size is limited, T eff obtained in this work are consistent with the range of T eff in which RSGs are expected to stay relatively long time. A larger sample of RSGs, especially those with a higher luminosity, would enable us to investigate the properties of the Galactic RSGs, e.g., comparing the T eff distribution more closely with various evolutionary models, determining the relation between the spectral type and T eff and so on.

SUMMARY AND FUTURE PROSPECTS
In this paper, we calibrated the empirical relations between twelve LDRs of Fe i lines and T eff , using nine solar-metal red giants observed with WINERED. These relations enabled us to determine T eff of red giants to a precision of ∼ 30 K in the best cases, i.e., early-M type giants with good S/N. We applied these relations to ten nearby RSGs and obtained T eff to a precision of ∼ 40 K. The estimated T eff are in good agreement with the values estimated with different methods in relevant literature: those with the TiO method by Levesque et al. (2005) and that of Betelgeuse on the basis of spectrointerferometry by Ohnaka et al. (2011) and Arroyo-Torres et al. (2013), as well as those expected in Geneva's stellar evolution model (Ekström et al. 2012). Our method uses only Fe i lines in the near-infrared (the Y and J bands) high-resolution spectra. Because of it, our method is expected to give more unbiased T eff of RSGs than the other published spectroscopic methods, which rely on lines of molecules and/or multiple atoms and inevitably suffer from some significant systematic uncertainties. First, our method is independent of chemical abundance ratios because our method relies on only one species, Fe i. Second, our method is expected to be less affected by the granulation than the methods using molecular bands because the granulation tends to vary the temperature structure of the upper atmospheric layers in particular and, hence, the strengths of molecular bands. Finally, our method is less affected by the surface gravity effect from which conventional LDR methods with multiple species suffer.
Though there are possible systematic biases of the temperatures derived with individual line pairs due to the surface gravity, microturbulent velocity and/or line broadening effects, it is expected that they cancel each other out when the temperatures derived with several pairs are averaged. In contrast, it is unclear whether the systematic bias due to the non-LTE effect, as large as ∼ 250 K for each pair, can be reduced by taking the average of individual estimates. Nevertheless, the final LDR temperatures using all the 11 pairs are well consistent within ∼ 100 K with the temperatures derived with the LDR ID (4), which is insensitive to the non-LTE effect. This fact indicates that the actual systematic bias due to the non-LTE effect is as small as ∼ 100 K. To further examine the possible systematic errors, in-depth consideration with reliable 3D non-LTE model spectra is desired.
Our sample is limited to solar-metal objects, red giants and RSGs, and the relations we obtained in this work can be applied to solar-metal objects only. In order to get rid of this limitation, the metallicity dependence of the LDRs in the Y and J bands needs to be examined. Observations with recent and/or future near-infrared high-resolution spectrographs with high sensitivities (e.g. WINERED; Ikeda et al. 2016Ikeda et al. , 2018 will provide sufficient-statistics data of RSGs at large distances, e.g., those in the inner/outer-Galaxy and some local-group galaxies like the Magellanic Clouds and more. Then T eff of these sources can be determined in this method, whereas interferometric measurements of such distant sources would be practically impossible. Obtaining T eff of RSGs in various environments is desirable to examine the potential environmental dependency of T eff on [Fe/H] and to test stellar evolution models further.