-
PDF
- Split View
-
Views
-
Cite
Cite
Takashi J Moriya, Ke-Jung Chen, Kimihiko Nakajima, Nozomu Tominaga, Sergei I Blinnikov, Observational properties of a general relativistic instability supernova from a primordial supermassive star, Monthly Notices of the Royal Astronomical Society, Volume 503, Issue 1, May 2021, Pages 1206–1213, https://doi.org/10.1093/mnras/stab622
- Share Icon Share
ABSTRACT
We present the expected observational properties of a general relativistic instability supernova (GRSN) from the 55 500 M⊙ primordial (Population III) star. Supermassive stars exceeding |$10^4\, \mathrm{M}_\odot$| may exist in the early Universe. They are generally considered to collapse through the general relativistic instability to be seed black holes to form supermassive (|$\sim 10^9\, \mathrm{M}_\odot$|) black holes observed as high-redshift quasars. Some of them, however, may explode as GRSNe if the explosive helium burning unbinds the supermassive stars following the collapse triggered by the general relativistic instability. We perform the radiation hydrodynamics simulation of the GRSN starting shortly before the shock breakout. We find that the GRSN is characterized by a long-lasting (550 d) luminous (|$1.5\times 10^{44}\, \mathrm{erg\, s^{-1}}$|) plateau phase with the photospheric temperature of around 5000 K in the rest frame. The plateau phase lasts for decades when it appears at high redshifts and it will likely be observed as a persistent source in the future deep near-infrared imaging surveys. Especially, the near-infrared images reaching 29 AB magnitude that can be obtained by Galaxy and Reionization EXplorer (G-REX) and James Webb Space Telescope (JWST) allow us to identify GRSNe up to z ≃ 15. Deeper images enable us to discover GRSNe at even higher redshifts. Having extremely red colour, they can be distinguished from other persistent sources such as high-redshift galaxies by using colour information. We conclude that the deep near-infrared images are able to constrain the existence of GRSNe from the primordial supermassive stars in the Universe even without the time domain information.
1 INTRODUCTION
Supermassive black holes (SMBHs) exceeding |$\sim 10^9\, \mathrm{M}_\odot$| are known to exist at z > 6 through high-redshift quasar observations (Fan et al. 2001, 2003; Willott et al. 2007; Mortlock et al. 2011; Morganson et al. 2012; Kashikawa et al. 2015; Wu et al. 2015; Bañados et al. 2016, 2018; Matsuoka et al. 2019). The age of the Universe at z > 6 is less than 1 Gyr. Forming SMBHs in 1 Gyr is very challenging and it is one of the frontiers in the modern astrophysics (see Haemmerlé et al. 2020; Inayoshi, Visbal & Haiman 2020 for recent reviews).
One possible path to form SMBHs in such a short time-scale is through Population III supermassive stars (SMSs) having |$10^4-10^6\, \mathrm{M}_\odot$|. Although the typical mass of Population III stars is predicted to be much less than |$10^4-10^6\, \mathrm{M}_\odot$| (e.g. McKee & Tan 2008; Hosokawa et al. 2011; Susa, Hasegawa & Tominaga 2014; Hirano et al. 2015; Sugimura et al. 2020), SMSs can be formed by preventing H2 cooling through, e.g. the intense Lyman–Werner ultraviolet background radiation photodissociating H2 (e.g. Omukai 2001; Oh & Haiman 2002; Bromm & Loeb 2003; Sugimura, Omukai & Inoue 2014; Chon et al. 2016; Regan et al. 2017; Chon, Hosokawa & Yoshida 2018). SMSs are predicted to collapse through general relativistic instability (e.g. Iben 1963; Chandrasekhar 1964; Fowler 1966; Osaki 1966; Shibata & Shapiro 2002; Shibata, Uchida & Sekiguchi 2016; Umeda et al. 2016; Uchida et al. 2017; Nagele et al. 2020), forming seed massive black holes (BHs). SMSs can be formed at z ≃ 15–20 (e.g. Agarwal et al. 2012; Dijkstra, Ferrara & Mesinger 2014; Habouzit et al. 2016) to leave the seed BHs and they can grow to SMBHs in 1 Gyr to explain the existence of SMBHs at z > 6.
The general relativistic instability of SMSs does not necessarily lead to their collapse to BHs (Fuller, Woosley & Weaver 1986; Montero, Janka & Müller 2012; Chen et al. 2014; Nagele et al. 2020). In particular, Chen et al. (2014) showed that a non-rotating Population III SMS with 55 500 M⊙ explodes as a supernova (SN) through the explosive helium burning following the collapse. Although SMSs causing such a general relativistic instability SNe (GRSNe) may be limited in a small mass range, they could be bright enough to be observed in the future near-infrared (NIR) transient surveys (Whalen et al. 2013). They also have distinctive chemical signatures that can be traced by the stellar archaeology of unusual extreme metal-poor stars (Johnson et al. 2013).
In this paper, we investigate the observational properties of the GRSN from the 55 500 M⊙ SMS presented by Chen et al. (2014). The observational properties of the GRSN were previously investigated by Whalen et al. (2013). They showed that future NIR transient surveys can discover the GRSN. In this work, we present our new radiation hydrodynamics simulation of the GRSN and suggest that it is not required to conduct transient surveys to discover the GRSN given its extremely long time-scale keeping its brightness. Having deep multiband images in NIR would be enough to identify them.
The rest of the paper is organized as follows. We first show our method of the GRSN investigation in Section 2. Then we discuss its properties in the rest frame in Section 3. Then we show its observational properties when it appears in the early Universe and discuss how to find it in Section 4. We conclude this paper in Section 5. We adopt the standard ΛCDM cosmology with |$H_0=70\, \mathrm{km s^{-1}}\, \mathrm{Mpc^{-1}}$|, ΩM = 0.3, and |$\Omega _\Lambda = 0.7$| throughout this paper.
2 METHOD
2.1 GRSN model
We take the GRSN explosion model of the 55 500 M⊙ Population III star presented by Chen et al. (2014). We adopt their 2D explosion model. The evolution of the 55 500 M⊙ star is followed by using the KEPLER 1D stellar evolution code (Weaver, Zimmerman & Woosley 1978). The star is found to explode with the 1D calculation with KEPLER (Chen et al. 2014). The radius of the star at the time of the explosion is 256 R⊙.
The 1D stellar structure at 600 s before the maximum compression at the centre in the 1D calculation, which is sufficiently long before the explosive burning leading to the explosion, is transferred to the multidimensional hydrodynamics code CASTRO (Almgren et al. 2010; Zhang et al. 2011). The subsequent hydrodynamic evolution with nucleosynthesis is followed until shortly before the shock breakout by 2D hydrodynamics. We adopt the result from the multidimensional simulation because of the strong mixing found in the GRSN which affect the abundance profile and, therefore, synthetic LCs (Chen et al. 2014). The explosion energy is |$9\times 10^{54}\, \mathrm{erg}$|. Because the core temperature (<109 K) and density (|$\lt 10^3\, \mathrm{g\, cm^{-3}}$|) is low at the explosive helium burning, little 56Ni is produced during the explosion and the total 56Ni mass in the ejecta is less than 0.1 M⊙.
We obtain the angle-averaged hydrodynamic and abundance profiles from the 2D simulation and use it as the initial condition for our 1D radiation hydrodynamics simulation described in the next section. For each physical quantity at a given radius, we take ten directions from polar angles of 0.1π, 0.2π, 0.3π, ..., π in the 2D result and then take the average value of the ten directions for the 1D calculation. In this way we approximately take the effect of multidimensional mixing found in the 2D calculation into account. Detailed chemical mixing patterns in the GRSN model are presented in fig. 4 of Chen et al. (2014). The angle-averaged density and abundance profiles are shown in Fig. 1. A hydrogen + helium layer exists from the surface to around 31 000 M⊙. The layers below are mostly composed of helium, oxygen, neon, and magnesium.

Top: Initial density profile for the GRSN LC calculation by STELLA. Bottom: Abundance profile for the GRSN LC calculation. The top axis shows the velocity of the mass coordinate during the homologous expansion.
2.2 Light-curve calculation
The numerical LC calculation is performed by using the 1D radiation hydrodynamics code STELLA (Blinnikov et al. 1998, 2000, 2006). In addition to the hydrodynamics equations, STELLA implicitly treats time-dependent equations of the angular moments of intensity averaged over a frequency bin with the variable Eddington method. STELLA calculates the spectral energy distributions (SEDs) at each time-step and obtains multicolour LCs by convolving filter functions with the SEDs. We adopt the standard 100 frequency bins from 1 Å to 5 × 104 Å on a log scale. STELLA has been widely used for the LC modelling of hydrogen-rich SNe (e.g. Blinnikov et al. 2000; Baklanov, Blinnikov & Pavlyuk 2005; Moriya et al. 2011, 2016, 2018b; Tominaga et al. 2011; Goldberg, Bildsten & Paxton 2019, 2020), including hydrogen-rich Population III SNe (Tolstov et al. 2016; Moriya et al. 2019). STELLA is also shown to provide very similar SN II LCs to those obtained by the Monte Carlo radiation transfer approach (Tsang et al. 2020).
We put the initial condition that is shortly before the shock breakout provided by the CASTRO simulation and no energy is artificially injected. We do not set mass cut because no compact remnant remains in the GRSN explosion.
Because we have the rest-frame SED information at each time-step, we shift the SEDs at given redshifts and apply filter functions at the observer frame to estimate the observational properties of the GRSN at high redshifts. We adopt the following filter functions in NIR when we discuss the observational properties of the GRSN. The mean wavelength of each filter is also mentioned. The H band filter (|$1.77\, \mathrm{\mu m}$|) from Euclid, which is the reddest filter in Euclid is adopted (Maciaszek et al. 2016). We take the J129 (|$1.29\, \mathrm{\mu m}$|), H158 (|$1.58\, \mathrm{\mu m}$|), and F184 (|$1.84\, \mathrm{\mu m}$|) filters from the Nancy Grace Roman Space Telescope (RST, previously known as WFIRST).1 The K band filter (|$2.15\, \mathrm{\mu m}$|) is adopted from ULTIMATE-Subaru,2 which is planned to have a similar K-band filter to the KS band filter of Subaru/MOIRCS.3 We take three filters from Galaxy and Reionization EXplorer (G-REX), which is a proposed wide-field surveyor at |$2-5\, \mathrm{\mu m}$| (A.K. Inoue, private communication), i.e. the F232 (|$2.32\, \mathrm{\mu m}$|), F303 (|$3.03\, \mathrm{\mu m}$|), and F397 (|$3.97\, \mathrm{\mu m}$|). These filters were previously proposed for the Wide-field Imaging Surveyor for High-redshift (WISH) satellite4 and are currently considered for G-REX. Finally, we take three broad filters from James Webb Space Telescope (JWST)/NIRCam5: F277W (|$2.77\, \mathrm{\mu m}$|), F356W (|$3.56\, \mathrm{\mu m}$|), and F444W (|$4.44\, \mathrm{\mu m}$|).
3 GRSN PROPERTIES IN THE REST FRAME
Fig. 2 shows the velocity evolution at the outermost layers at the beginning of the LC calculation. The initial condition (t = 0) is shortly before the shock breakout and the shock front continues to travel towards the progenitor surface. The shock breakout occurs at |$t\simeq 500\, \mathrm{s}$|. After the shock breakout, the homologous expansion of the ejecta is set at |$t\simeq 5000\, \mathrm{s}$| and the velocity of each mass shell is fixed from the moment. The velocity of the mass coordinate during the homologous expansion is shown in Fig. 1.

Velocity evolution at the surface of the progenitor shortly after the beginning of the STELLA calculation.
The bolometric LC of the GRSN is shown in Fig. 3. Fig. 4 presents the velocity, temperature, mass coordinate, and radius of the photosphere which is where the Rosseland mean optical depth becomes 2/3. We first see the shock breakout signal peaking at 14 s with |$10^{47}\, \mathrm{erg\, s^{-1}}$|. The shock breakout signal lasts for several hundred seconds which roughly correspond to the light-crossing time of the progenitor (600 s).

Bolometric LC of the GRSN. The horizontal axes in the top panel are in the log scale and those in the bottom panel are in the linear scale. The top panel shows the analytically estimated bolometric LC decay rate after the shock breakout (∝ t−0.35) and the bottom panel shows the nuclear decay rate of 56Co.

The bolometric LC evolution after the shock breakout roughly follows ∝ t−0.35 as expected from the analytical models (e.g. Chevalier & Fransson 2008; Rabinak & Waxman 2011; Kozyreva et al. 2020). The bolometric LC keeps declining until 17 d when the outermost layers get close to the hydrogen recombination temperature at |$\simeq 6000\!-\!5000\, \mathrm{K}$|. The subsequent bolometric LC evolution can be interpreted simply with L ∝ R2T4, where L is the bolometric luminosity, R is the photospheric radius, and T is the photospheric temperature. The bolometric luminosity increases from 17 to 250 d when the recession velocity of the photosphere is slower than the hydrodynamic expansion velocity of the ejecta. In other words, the photospheric temperature is set at the recombination temperature but the photospheric radius increases (Fig. 4), making the bolometric luminosity to increase.
The bolometric luminosity increases until 250 d when the recession velocity of the photosphere and the hydrodynamic expansion velocity matches. From this time, the photosphere is kept at the same radius and it also keeps the same recombination temperature. Thus, the bolometric luminosity becomes constant and the plateau phase appears. This physical condition is the same as that found in Type IIP SNe during their plateau phase (Grassberg, Imshennik & Nadyozhin 1971; Falk & Arnett 1977). The plateau luminosity (|$\simeq 1.5\times 10^{44}\, \mathrm{erg\, s^{-1}}$|) is much larger than those observed in Type IIP SNe (Bersten & Hamuy 2009) because of the extremely larger explosion energy. Exceeding |$10^{44}\, \mathrm{erg\, s^{-1}}$|, the GRSN can be regarded as a member of superluminous SNe (Moriya, Sorokina & Chevalier 2018a; Gal-Yam 2019, for recent reviews). The plateau continues up to around 800 d when the photosphere reaches the bottom of the layers containing hydrogen (see the mass coordinate of the photosphere in Fig. 4). The hydrogen recombination ends at this point and the photospheric temperature drops suddenly. The plateau duration is about 550 d in the rest frame which is also much longer than those found in Type IIP SNe (about 100 d, Bersten & Hamuy 2009; Anderson et al. 2014). The total radiated energy during the plateau phase amounts to |$\simeq 7\times 10^{51}\, \mathrm{erg}$|.
The plateau phase is followed by the sudden drop in luminosity and then the ‘tail’ phase starts. The bolometric luminosity keeps decreasing during the tail phase. The tail phase luminosity is not powered by the radioactive decay of 56Co unlike the case of Type IIP SNe because of the small amount of 56Ni (less than 0.1 M⊙) in the ejecta. The LC decline rate is different from that of 56Co as shown in Fig. 3. The luminosity source in this tail phase is still thermal energy from the initial explosion as in the previous hydrogen recombination phase. About a half of magnesium, which is the most abundant element in the hydrogen-free core (Fig. 1), remain ionized even after the plateau phase because of its low ionization energy. Therefore, the opacity at the hydrogen-free core remains high (of the order of |$0.01\, \mathrm{cm^2\, g^{-1}}$|) even after the hydrogen recombination phase and the photosphere is kept at the surface of the hydrogen-free layers during the epochs presented in Fig. 3. This is why the photospheric radius increase linearly after the hydrogen recombination phase, i.e. the photosphere is kept at the surface of the hydrogen-free core. The bolometric luminosity decreases as the temperature decreases due to the adiabatic cooling. In the photospheric temperature evolution, we find a sharp increase at around 1150 d. The exact reason for the jump is not clear and it could be a numerical artefact.
Our bolometric LC behaviour is different from that presented in the previous study in Whalen et al. (2013) in which the same initial explosion model is used but the numerical LC calculation is performed in a different code. The bolometric LC from Whalen et al. (2013) shows a slow luminosity decline starting from the shock breakout peak. Their LC model does not show the plateau phase caused by the hydrogen recombination as we found in our model. We believe that the existence of the plateau phase is expected from the presence of the massive hydrogen-rich envelope as in the case of Type IIP SNe (Grassberg et al. 1971; Falk & Arnett 1977). We have compared our results with predictions of the analytic model by Nagy & Vinkó (2016) (see also Szalai et al. 2019 for the update of that model). We find a good agreement of synthetic light curves produced by STELLA with this simple model on the plateau stage. The earlier bolometric LC behaviour in our model is also consistent with the analytic models as discussed earlier in this section. The reasons for the difference between our LC model and the model in Whalen et al. (2013) remain unclear.
4 DISCOVERING GRSNE IN THE FUTURE NIR SURVEYS
We discuss the observational properties of GRSNe when they appear in the early Universe based on the GRSN model presented in the previous section.
4.1 GRSNe as persistent sources
Fig. 5 shows the GRSN LCs at z = 15 with the G-REX and JWST filters. The observed wavelengths correspond to |$0.1-0.3\, \mathrm{\mu m}$| in the rest frame. The SEDs at |$0.1-0.3\, \mathrm{\mu m}$| are strongly affected by metal absorption and have steep gradients. This causes large differences in the observed magnitudes even if the differences in the filter central wavelengths are small. This effect of absorption becomes larger with higher redshifts (cf. Fig. 6).

GRSN LCs at z = 15 with the G-REX and JWST filters in the observer frame.

Magnitudes of the GRSN at the plateau phase as a function of redshift. The SED at 500 d after shock breakout presented in Section 3 is used to estimate the plateau magnitudes. The horizontal dashed line shows a typical limiting magnitude of planned imaging surveys (29 mag).
The most prominent feature is the long plateau phase lasting for three decades. It originate from the luminous plateau phase lasting for 550 d in the rest frame (Section 3). Because of the time dilution, the long plateau phase becomes even longer and lasts for decades when GRSNe appear at high redshifts.
The long plateau lasting for more than a decade would be observed as a persistent source rather than a transient during the NIR survey, because the major NIR survey missions are conducted by satellite telescopes and their operation period is usually set to 5 yr. Fig. 6 shows the predicted magnitudes of the GRSN plateau as a function of redshift. We find that GRSNe up to z ≃ 15 can be observed with G-REX and JWST if a deep imaging survey reaching at least |$\simeq 29\, \mathrm{AB\, mag}$| is conducted. The deeper images enable us to reach higher redshifts (Fig. 6).
To identify non-transient GRSNe, we need to distinguish them from other persistent sources such as high-redshift galaxies. Fig. 7 shows the colour–colour diagrams in which GRSNe at plateau, high-redshift galaxies, high-redshift quasars, and dwarf stars in our Galaxy. The high-redshift galaxy templates are produced with the stellar radiation provided by BPASS (v2.1; Eldridge et al. 2017) in addition to nebular emission (both lines and continua) as calculated in a self-consistent manner with Cloudy (see Nakajima et al. 2018 for more details). Here we adopt ‘young’ and ‘evolved’ galaxies, with a continuous star formation history at the ages of 10 and 500 Myr, the metallicity of 0.1 and 0.5 solar metallicity, and the dust reddening of E(B − V) = 0.0 and 0.2 assuming the Calzetti et al. (2000)’s extinction law, respectively. The young population shows intense emission lines with a rest-frame equivalent width of [O iii]λ5007 + H β as large as ∼1000 Å, which corresponds to extreme emission line galaxies such as green pea galaxies (e.g. Cardamone et al. 2009). The evolved population represents continuum-selected galaxies like Lyman-break galaxies (e.g. Steidel et al. 2014). Both galaxy populations are then placed at redshifts of z = 6–12 with an IGM attenuation prescribed by Inoue et al. (2014). At these redshifts, strong optical emission lines such as H α and [O iii] λ5007 fall in the NIR filters of G-REX and JWST/NIRCam. The nebular component in our galaxy-SEDs is thus essential to obtain realistic predictions of NIR colours of high-redshift galaxies, which are inferred to present intense H α and [O iii]+H β emission lines (e.g. Smit et al. 2014; Roberts-Borsani et al. 2016). The quasar models are based on a composite of spectra of quasars at z = 1–2 observed in the rest-frame ultraviolet to near-infrared wavelength (Selsing et al. 2016). The composite is redshifted and IGM-attenuated in the same way as the galaxy models to mock up the high-redshift quasars. The dwarf stars are plotted by using the NIR spectra obtained by SpeX (Rayner et al. 2003) and presented in Cushing, Rayner & Vacca (2005), Rayner, Cushing & Vacca (2009).6 We took all the dwarf spectra with the spectral coverage up to |$4\, \mathrm{\mu m}$| for G-REX and |$5\, \mathrm{\mu m}$| for JWST and convolved the filter functions. When there is a gap in the dwarf spectra caused by the telluric absorption, we simply interpolated the gap.

Colour–colour diagrams of the GRSN at the plateau phase (500 d after the shock breakout). We also put high-redshift galaxies and dwarf stars in our Galaxy for comparison.
The colour–colour diagram (Fig. 7) shows that GRSNe are much redder than high-redshift galaxies, high-redshift quasars, and dwarf stars and they can be easily distinguished from GRSNe at the plateau phase. This is because the temperature of the plateau is only |$\simeq 5000-6000\, \mathrm{K}$| in the rest frame and they become extremely red if they are at high redshifts. This temperature at the plateau phase is determined by the hydrogen recombination temperature. Thus, the colour during the plateau phase is a solid prediction that does not depend much on the actual mass and energy of the GRSN explosions.
4.2 Discovering GRSNe with NIR transient surveys
As discussed previously by Whalen et al. (2013), it is also possible to find GRSNe with NIR transient surveys. As presented in Fig. 3, GRSNe are brightest shortly after the shock breakout and their luminosity keeps declining until 17 d. GRSNe at this phase can be observed at z ≳ 15 thanks to their large luminosity and high temperature. Fig. 8 shows examples of the GRSN LCs during the declining phase. The LCs rise in about a month and declines. Fig. 9 shows the peak luminosity of the first peak and the rise time to the peak in the NIR filters as a function of redshift. NIR transient surveys reaching 28.5 AB mag are required to find high-redshift GRSNe. The time-scale is about a month. Thus, deep NIR transient surveys reaching down to 28.5 AB mag with a cadence of around 10 d are required to identify high-redshift GRSNe shortly after the explosion.

GRSN LCs at z = 15 and 20 during the cooling phase after the shock breakout with the selected NIR filters.

Peak magnitudes and rise times of GRSNe during the cooling phase shortly after the shock breakout as a function of redshift.
5 CONCLUSIONS
We have presented the observational properties of the GRSN from the 55 500 M⊙ Population III star. The progenitor has the radius of 256 R⊙ and the explosion energy is |$9\times 10^{54}\, \mathrm{erg}$|. We take the results of the GRSN explosion simulation up to shortly before the shock breakout conducted previously by Chen et al. (2014). We put it as the initial condition to the radiation hydrodynamics code STELLA to investigate its observational properties. The overall observational properties of the GRSN is similar to those of Type IIP SNe, but the GRSN has the much longer (550 d in the rest frame) and more luminous (|$1.5\times 10^{44}\, \mathrm{erg\, s^{-1}}$|) plateau phase. The photospheric temperature during the plateau is characterized by the hydrogen recombination temperature (|$\simeq 5000-6000\, \mathrm{K}$|) as is the case for Type IIP SNe. The plateau phase was not predicted to exist in the previous study by Whalen et al. (2013) in which the same initial explosion model is used to investigate the GRSN properties, but the existence of the plateau phase is expected from the presence of the massive hydrogen-rich envelope as in the case of Type IIP SNe. The plateau phase lasts for many decades when the GRSN appears at high redshifts and the GRSN is not recognized as a transient in the future NIR imaging surveys. However, it can be distinguished from other persistent sources such as high-redshift galaxies and local dwarf stars by using colour information. The deep NIR images reaching 29 AB mag obtained by G-REX and JWST will allow us to search for the GRSN up to z ≃ 15. The deeper images allow us to reach higher redshifts. The NIR transient surveys reaching |$28-28.5\, \mathrm{AB mag}$| with a cadence of 10 d may identify the GRSN at z ≳ 15 during the adiabatic cooling phase after the shock breakout as a transient.
Because the exact mass range of the GRSN explosions is unclear, it is difficult to predict the expected number of the GRSN detection in the deep NIR images obtained by the future NIR surveys. Still, the extremely red colour of high-redshift GRSNe makes them easy to identify. The existence or lack of the GRSNe can be constrained in the future deep NIR images without conducting transient surveys. Such a constraint provides valuable information on the fate of SMSs possibly leading to the formation of SMBHs in the early Universe.
ACKNOWLEDGEMENTS
We thank the anonymous referee for constructive comments that improved this paper. TJM thanks the organisers of the First Star VI conference at Universidad de Concepción, Chile, where this work was initiated. TJM is supported by the Grants-in-Aid for Scientific Research of the Japan Society for the Promotion of Science (JP18K13585, JP20H00174). KC acknowledges support from the Ministry of Science and Technology (Taiwan, R.O.C.) grant number MOST 107-2112-M-001-044-MY3. SB is supported by grant RSF 19-12-00229 for the development of STELLA code. This research has been supported in part by the RFBR (19-52-50014)-JSPS bilateral program. Numerical computations were in part carried out on PC cluster at Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.
Footnotes
We obtained the spectra at http://irtfweb.ifa.hawaii.edu/ spex/IRTF_Spectral_Library/.