© The Authors 2026

1. Introduction

The formation of main-sequence stars from young stellar objects (YSOs) is accompanied by gas accretion, which often occurs in episodic outbursts (Hartmann & Kenyon 1996; Bae et al. 2014; Hartmann et al. 2016). During such events, the optical brightness of YSOs typically increases by 5−6 mag in the B/V bands, corresponding to an increase of about 2 orders of magnitude in bolometric luminosity, over a span of several months, with elevated brightness persisting for months, years, or even decades (Fischer et al. 2023). These episodic accretion sources are generally classified into two categories: FUors and EXors. FUors typically exhibit long-duration outbursts lasting from years to decades, whereas EXors display shorter, recurrent events on timescales of months to a few years. Regarding spectroscopic features, FUors exhibit absorption-dominated spectra, including prominent features from CO, water vapor, and various metal lines (Reipurth et al. 2012; Hillenbrand et al. 2018; Guo et al. 2024b). In contrast, EXors are characterized by strong emission lines, including hydrogen and helium recombination lines, metallic collisional excitation lines, and CO emission (Audard et al. 2014). A comprehensive observational overview of the spectroscopic characteristics of FUors is provided by Connelley & Reipurth (2018).

A long-standing issue in YSO studies is the discrepancy between the accretion rates inferred from luminosity measurements and those estimated from star-formation duration, commonly known as the luminosity problem (Evans et al. 2009). Episodic accretion has been proposed as a viable solution: during quiescent phases, YSOs exhibit lower accretion rates and, consequently, lower luminosities; however, during episodic outbursts, they can accrete substantial amounts of gas in short periods (Vorobyov & Basu 2006; Dunham et al. 2010). In FUor-type outbursts, luminosity is thought to be primarily powered by the accretion disk.

The radiation from YSOs mainly consists of active radiative components generated by stellar activity and accretion, and passive emission radiative components generated by the surrounding dust and gas after absorbing this energy. The active radiative components of YSOs arise from three main sources (Robitaille et al. 2006; Zhu et al. 2007; Robitaille 2017; Rodriguez & Hillenbrand 2022): (1) stellar contraction and nuclear processes such as deuterium burning, (2) viscous dissipation in the accretion disk, and (3) magnetospheric accretion columns. In contrast, passive radiative components originate from a dusty outer disk heated by stellar irradiation. Both active and passive radiation can be reprocessed and redistributed by the surrounding envelope (Ulrich 1976).

In the early days, discovery and follow-up observations of YSO outbursts were mainly carried out in the optical band (Evans et al. 2009; Findeisen et al. 2013). Nonetheless, significant dust extinction significantly influences the sample’s completeness, particularly impacting deeply embedded YSOs (Lucas et al. 2020; Wang et al. 2023). In recent years, with the development of infrared time-domain surveys, many projects have also begun to use infrared bands to search for and study these sources (Morales-Calderón et al. 2011; Contreras Peña et al. 2017; Park et al. 2021; Lucas et al. 2024; Li & Wang 2024). Li & Wang (2024) used the unTimely WISE catalog, a full-sky time-domain catalog based on unWISE (unblurred coadds from the Wide-field Infrared Survey Explorer, WISE), to search for variable sources in the sky area 295° < l < 350° and −1° < b < 1°; they identified numerous young stellar outburst sources by examining their color, the amplitude of their infrared light variations, and the pattern of this light fluctuation. In this study, we performed near-infrared spectroscopic follow-up observations on several FUor candidates identified by Li & Wang (2024), specifically those labeled V10, V191, V450, and V978, using the Immersion Grating Infrared Spectrometer (IGRINS) instrument on Gemini South. Section 2 describes the light curves and spectroscopic observations. In Sect. 3 we explain how we applied a disk model to near-infrared spectra and mid-infrared photometry to derive the essential physical parameters of these candidates. The fitting results, along with the properties of the sources identified via their infrared variability, are discussed in Sect. 4. Lastly, in Sect. 5 we present a summary of our main findings.

2. Light curves and spectra

2.1. Ancillary photometry

2.1.1. WISE

The NASA mission WISE performed an all-sky survey in the mid-infrared, repeatedly imaging the entire sky in four bands centered at 3.4, 4.6, 12, and 22 μm (W1–W4; Wright et al. 2010) during its first two years of operation, and subsequently in only the W1 and W2 bands following its reactivation in 2014. Its multi-epoch observing strategy makes WISE a valuable resource for time-domain photometry in the mid-infrared. In this study, we mainly utilized WISE observations to examine the infrared variability of our sources. Specifically, we used time-resolved photometric data from the Near-Earth Object Wide-Field Infrared Survey Explorer (NEOWISE) catalogs in the W1 and W2 bands, which offer a time baseline exceeding 10 years and a typical sampling interval of about six months (Mainzer et al. 2011, 2014).

2.1.2. VISTA

The Visible and Infrared Survey Telescope for Astronomy (VISTA) is a 4.1-m wide-field near-infrared survey telescope located at Paranal Observatory and operated by the European Southern Observatory (ESO). We used near-infrared photometric data obtained with VIRCAM, the VISTA InfraRed CAMera, which provides imaging in the Z, Y, J, H, and K_s bands.

The data employed in this work are drawn from the VISTA Variables in the Vía Láctea (VVV) survey, a public near-infrared survey targeting the Galactic bulge and adjacent disk regions (Saito et al. 2010; Minniti et al. 2010). For the variability analysis in this paper, we restricted the VVV data to the K_s band, which has the highest observational cadence. The remaining bands do not provide adequate temporal coverage for variability studies and are therefore not considered here; measurements in those bands are reported in Li & Wang (2024).

2.1.3. SPHEREx

The Spectro-Photometer for the History of the Universe, Epoch of Reionization and Ices Explorer (SPHEREx), a NASA space-based all-sky spectroscopic survey mission, provides low-resolution spectra over the wavelength range of approximately 0.75−5.0 μm with a pixel scale of 6.2 arcsec (Bock et al. 2026). We required a signal-to-noise ratio flux/flux_err > 10 and excluded measurements with high-order SPHEREx quality flags indicating spectrophotometric contamination (flags ≳ 4.3 × 10⁹), corresponding to cases where bad pixels were involved in the spectrophotometry calculation. An additional complication arises for the source V450. At an angular separation of 3.06 arcsec from V450, there is a nearby contaminating object (Gaia DR3 5885045277964196096). This object has magnitudes of H = 13.181 mag and K = 12.884 mag from 2MASS (Skrutskie et al. 2006), and I1 = 13.003 mag and I2 = 12.790 mag from the Spitzer Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE) survey (Werner et al. 2004). Given the SPHEREx pixel scale is 6.2 arcsec, the measured SPHEREx fluxes are substantially contaminated by this nearby source over several wavelength intervals. Consequently, we applied the following corrections: (1) we discarded the spectral regions from 0.70 μm and 1.30 μm, where the contaminant dominates the SPHEREx flux; (2) in the range 1.3−1.7 μm, we subtracted the contaminant contribution inferred from its H-band flux; (3) we removed the 2.3−3.2 μm interval because no reliable flux measurements for the contaminant are available there; (4) between 3.6 and 4.2 μm, we subtracted the contaminant contribution estimated from the Spitzer Infrared Array Camera (IRAC) I1 band; (5) between 4.2 and 5.0 μm, we subtracted the contaminant contribution estimated from the Spitzer IRAC I2 band. The SPHEREx photometric measurements used in this work have modified Julian dates (MJDs) primarily in the range 60 800−60 900, corresponding to approximately 500 days after the spectroscopic observations.

2.2. Light curve

Figure 1 illustrates the light curves for the four sources. Each of these sources experienced outbursts during the NEOWISE survey, with infrared brightness rising by more than 2 mag. This threshold, as indicated by Li & Wang (2024), is used to identify YSOs exhibiting large infrared outbursts. The large amplitudes and long duration of those outbursts align with the FUor classification.

Fig. 1. Near-infrared and mid-infrared light curves of the four FUor candidates. We show the light curves from Ks (green), W1 (light green), and W2 (orange) bands (Minniti et al. 2010; Mainzer et al. 2014; Schlafly et al. 2019; Meisner et al. 2023). The vertical red and blue lines mark the dates of spectroscopic observations presented in this paper and Guo et al. (2024a).

2.3. Spectra

To confirm this classification and characterize the properties of these sources, we conducted near-infrared spectroscopic observations using the Gemini South Telescope. The International Gemini Observatory comprises twin 8.1-meter optical/infrared telescopes located at two of the best observing sites on Earth (Mace et al. 2018). For our observations, we used IGRINS (Park et al. 2014), which simultaneously covers the full wavelength range from 1.45 to 2.45 μm in a single exposure, with a resolving power of R ∼ 45 000 and a slit size of 0.34″ × 0.5″. The spectra were reduced using the IGRINS Pipeline Package (Kaplan et al. 2024). They were then de-blazed and corrected for velocity as described in Sect. 3.2. We performed absolute flux calibration of the spectra by matching them to the fluxes measured by SPHEREx. For all targets except V978, the SPHEREx photometry provides a stable reference across the near-infrared bands. For V978, however, the SPHEREx photometric points in the H and K bands exhibit unusually large scatter, preventing a reliable flux calibration. In this case, we instead calibrated the spectrum using the telluric standard star observed during the same night. The resulting spectra are shown in Fig. 2. Among these objects, V10, V191, and V978 were observed with X-shooter on the Very Large Telescope, and their median-resolution near-infrared spectra were presented by Guo et al. (2024a); these sources correspond to L222_1, L222_10, and L222_18, respectively. Details of the observations are listed in Table 1.

Fig. 2. Near-infrared spectra of four FUor candidates. The spectra are vertically shifted for clarity and have been convolved to a resolving power of R = 1200. Prominent atomic absorption features are marked by vertical dashed lines. We show the Brackett series (red), the CO band heads Δν = 2 (orange) and 3 (green), and the broad H2O absorption band (horizontal gray bar).

All spectra exhibit the following common features, consistent with features typically observed in FUor-type objects (Connelley & Reipurth 2018; Siwak et al. 2025): (1) Strong CO Δv = 2 absorption bands. Prominent CO overtone (Δv = 2) absorption is detected in all spectra. This feature is widely regarded as one of the defining spectral characteristics of FUor–type sources, originating from the hot inner regions of the accretion disk. (2) Water vapor absorption bands. Strong water vapor bands are present at both ends of the H-band window, shaping the H-band continuum into a characteristic triangular or “peaky” profile. Classical FUors show such a triangular H-band continuum. This feature is also considered a hallmark of FUor-type sources. (3) Hydrogen absorption lines. In the H and K bands, the spectral coverage includes the hydrogen Brackett series. Clear Brackett absorption lines are observed, in particular Brγ, which is prominently detected in absorption. (4) Weak metal absorption lines. Weak absorption features from neutral metals, such as Na I (2.209 μm in vacuum) and Ca I (2.261 μm), are present. These features are characteristic of late-M–type spectra and are commonly observed in FUor objects. (5) Lack of emission lines. In contrast to EXor-type objects, which often display strong emission lines, FUor-type sources generally show little to no line emission. Consistent with this expectation, emission lines are largely absent in our spectra, although Brγ may exhibit a P Cygni–like profile in some cases.

3. Spectral fitting and results

3.1. Disk model

During episodic accretion, the near-infrared emission is dominated by the accretion disk, including both the viscous inner disk and the passively heated outer disk (Rodriguez & Hillenbrand 2022; Liu et al. 2022). We adopted the modeling approach described in Liu et al. (2022) to fit the observed spectra with theoretical models and to extract the relevant physical parameters. Below we provide a brief overview of the model.

3.1.1. Viscous accreting disk

We adopted the viscous accretion disk model of Shakura & Sunyaev (1973), which assumes that the disk is in steady state, geometrically thin and optically thick. As a consequence of the differential Keplerian rotation, the disk experiences viscous heating, which leads to a temperature gradient across the radius. This distribution of temperature is expressed by

$$\begin{array}{c}\hfill T^4(r)=\frac{3G{M}_{\ast }\dot{M}}{8\pi \sigma r^3}[1-{\left(\frac{R_{\mathrm{in}}}{r}\right)}^{1/2}],\end{array}$$(1)

where M_* is the mass of the central star, Ṁ is the disk accretion rate, and R_in is the inner boundary of the disk. We denote the outer radius of the viscous disk as R_out and treated it as a free parameter. The temperature profile described by the above equation formally peaks at R = 1.36 R_in. We therefore assumed that the temperature is capped at this maximum value, T_max, and remains constant for R < 1.36 R_in, following Zhu et al. (2007) and Rodriguez & Hillenbrand (2022). To obtain the spectrum of the gas component of the viscous disk, we employed the approach outlined by Liu et al. (2022). This involves using the BT–Settl grid of theoretical stellar spectra, which incorporates solar AGSS2009 abundances, to model the emission from annular regions at varying radii, each of which is associated with a particular temperature (Allard 2014). The atmosphere of the accretion disk should correspond to low surface gravity. In the temperature range of 1400 K to 2000 K, we utilized a surface gravity with log₁₀g = 3.5, which represents the lowest gravity templates available. Between 2000 K and 24 000 K, we employed templates characterized by log₁₀g = 1.5. For temperatures above 24 000 K or below 1400 K, we utilized a blackbody spectrum.

To account for the Doppler broadening caused by Keplerian rotation, it is essential to convolve the radiation emitted from each annulus of the accretion disk with a rotational velocity profile. The convolution kernel is specified as

$$\begin{array}{cc}& \varphi (\lambda \prime ,\lambda )={[1-{\left(\frac{\lambda \prime -\lambda }{{\lambda }_{\max }}\right)}^2]}^{-1/2}\hfill \end{array}$$(2)

$$\begin{array}{cc}& {\lambda }_{\max }=\lambda \frac{sini}{c}\sqrt{\frac{G{M}_{\ast }}{R}},\hfill \end{array}$$(3)

where i is the inclination angle of the disk. Assuming that the intrinsic stellar model spectrum at temperature T is given by F_{BT − Settl}(λ, T), the convolved flux at radius r becomes

$$\begin{array}{c}\hfill F_{\mathrm{conv}}(\lambda ,r)=\frac{\int \varphi (\lambda \prime ,\lambda )F_{\mathrm{BT}-\mathrm{Settl}}(\lambda \prime ,T(r))d\lambda \prime }{\int \varphi (\lambda \prime ,\lambda )d\lambda \prime }·\end{array}$$(4)

The total observed spectral luminosity from the viscous accretion disk is then calculated as

$$\begin{array}{c}\hfill L_{\mathrm{vd}}(\lambda )={\displaystyle {\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}}2\pi F_{\mathrm{conv}}(\lambda ,T(r))rdr.\end{array}$$(5)

3.1.2. Passively heated disk

For the outer disk, we adopted a passively heated, flared disk geometry, in which the vertical scale height follows

$$\begin{array}{c}\hfill \frac{H}{r}=0.2{\left(\frac{r}{R_{\mathrm{out}}}\right)}^{\beta },\end{array}$$(6)

where H is the disk scale height, r is the cylindrical radius, R_out is the outer radius of the viscous disk, and β is the flaring index (Carvalho et al. 2023, 2024a). In a passively heated disk, the temperature structure is set by radiative equilibrium between irradiation from the inner disk and re-emission by the disk surface. The incident flux at radius r can be written as

$$\begin{array}{c}\hfill F_{\mathrm{pd}}(r)\propto \frac{L_{\mathrm{irr}}}{4\pi r^2}sin\theta ,\end{array}$$(7)

where L_irr is the effective luminosity of the irradiation source and θ is the grazing angle of incidence.

For small angles, the grazing angle can be approximated as

$$\begin{array}{c}\hfill \theta \simeq \frac{\mathit{dH}}{\mathit{dr}}-\frac{H}{r}·\end{array}$$(8)

Substituting the assumed disk geometry, H ∝ r^1 + β, yields

$$\begin{array}{c}\hfill \theta \propto \beta r^{\beta }.\end{array}$$(9)

Balancing the absorbed irradiation with blackbody re-emission,

$$\begin{array}{c}\hfill \sigma T^4(r)\propto F_{\mathrm{pd}}(r),\end{array}$$(10)

the radial temperature profile of the passive disk follows

$$\begin{array}{c}\hfill T(r)\propto r^{-\alpha },\alpha =\frac{2-\beta }{4}·\end{array}$$(11)

In this framework, the temperature gradient α is not an independent parameter but is uniquely determined by the disk flaring index β. We treated α as a free parameter. We fixed the outer radius of the passive disk to 1000 R_⊙. Assuming the disk radiates locally as a blackbody, the flux density emitted by the passive disk at frequency λ can be computed by integrating the contribution from each annulus,

$$\begin{array}{c}\hfill L_{\mathrm{pd}}(\lambda )={\displaystyle {\int }_{R_{\mathrm{out}}}^{1000R_{\odot }}}2\pi [\pi B_{\lambda }(T(r))]r\mathrm{d}r,\end{array}$$(12)

where B_λ(T) is the Planck function, and R_out is the inner radius of the passive disk, which is set by the outer edge of the viscous disk.

3.1.3. Extinction from the envelope and interstellar medium

For early-stage YSOs, the central star and its accretion disk are embedded within a collapsing, dust-filled envelope, which causes substantial extinction. In addition, because most YSOs form inside molecular clouds, the radiation they produce experiences further extinction as it passes through the surrounding large-scale material.

We adopted the extinction curve of Fitzpatrick (1999) with R_V = 3.1 for wavelengths shorter than 3.3 μm. For wavelengths longer than 3.3 μm, we used the extinction prescription given in Gordon et al. (2021). The observed flux from the viscous accretion disk can therefore be written as

$$\begin{array}{c}\hfill F_{\mathrm{model}}(\lambda )=\frac{(L_{\mathrm{vd}}+L_{\mathrm{pd}})\cos i}{\pi d_{\ast }^2}\times {10}^{-A_V\times f(\lambda )},\end{array}$$(13)

where d_* denotes the distance from the YSO to Earth, A_V represents the extinction coefficient, i is the inclination angle of the disk, and f(λ) is the normalized extinction curve specified by the chosen model.

As discussed above, the model contains several unknown parameters, including the accretion rate (Ṁ), stellar mass (M_*), inner radius of the viscous disk (R_in), outer radius of the viscous disk (R_out), the radial temperature gradient of the passive disk parameterized by α, inclination angle (i), extinction (A_V), and distance (d_*). The distance can be derived by analyzing spectral lines to obtain radial velocities, which, when combined with the Galactic rotation curve, yield an estimate of the distance to the YSO (Wenger et al. 2018). However, this method is subject to significant uncertainties arising from both the uncertainty in the Galactic rotation curve and the peculiar velocity of the YSO relative to Galactic rotation. As shown by Eq. (13), the uncertainty in distance affects the overall flux uniformly across all wavelengths, acting as a global scaling factor similar to the effect of the inclination angle through the cos i term.

In addition, the SPHEREx photometric observations were obtained approximately 500 days after the spectroscopic observations. Given the intrinsic variability of FUor-type objects, the use of non-contemporaneous SPHEREx photometry for flux calibration may introduce additional systematic offsets. We therefore absorbed the combined effects of distance uncertainty, inclination, and potential temporal variability between photometric and spectroscopic observations into a single multiplicative scale parameter in the fitting procedure. We therefore express the observed flux as

$$\begin{array}{c}\hfill F_{\mathrm{model}}(\lambda )=\frac{(L_{\mathrm{vd}}+L_{\mathrm{pd}})}{\pi d_{\ast }^2}\times {10}^{-A_V\times f(\lambda )}\times scale.\end{array}$$(14)

The parameters M_* and Ṁ have a degenerate effect on the temperature distribution and total flux. Hence, we replace them with a single combined parameter M_∗Ṁ to minimize this degeneracy. In summary, there are seven free parameters in the model: M_∗Ṁ, $sini\sqrt{M_{\ast }}$, R_in, R_out, α, A_V, and scale.

The ranges of the free parameters adopted in the fitting procedure are guided by physical considerations and previous studies of FUor-type systems. The combined parameter M_∗Ṁ, which is units of M_⊙² yr⁻¹, is allowed to vary such that log₁₀(M_∗Ṁ) ranges between −7 and −3, covering the expected range of accretion rates during FUor outbursts for low-mass YSOs. The inner radius of the viscous disk, R_in, is constrained to 1.5 − 6 R_⊙, consistent with disk truncation near the stellar surface, while the outer radius R_out is limited to 60 − 300 R_⊙. The parameter $sini\sqrt{M_{\ast }}$ is restricted to 0 − 1.5 M_⊙^1/2, corresponding to typical stellar masses of low-mass YSOs. The visual extinction A_V is allowed to vary between 0 and 25 mag, encompassing both lightly reddened and deeply embedded sources. The temperature gradient parameter α, describing the radial temperature profile of the passive disk (T ∝ r^−α), is permitted to vary between 0 and 3, enabling the model to explore a wide range of disk flaring and irradiation geometries. Finally, the multiplicative scale parameter is allowed to vary between 0.4 and 1.2 to account for uncertainties in distance, inclination effects, and potential systematic offsets introduced by non-contemporaneous photometric and spectroscopic observations.

3.2. Radial velocity and distance

To determine the line-of-sight velocity of the spectra, we first visually inspected a set of viscous disk model spectra and selected those whose CO absorption profiles exhibit similar line spacing to the observed spectra. Both the model and observed spectra were then continuum-normalized, and a cross-correlation function (CCF) analysis was performed. The CCF was computed over the wavelength range 2.2850−2.3500 μm, which corresponds to the CO Δv = 2 (2−0) absorption band. The procedure is illustrated in Fig. 3.

Fig. 3. Example of wavelength alignment of observed and template spectra for radial velocity determination. The observed spectrum (gray) and the viscous model spectrum (orange) with $M_{\ast }\dot{M}={10}^{-4.5}{M}_{\odot }^2$/yr, $sini\sqrt{M_{\ast }}=0.02\sqrt{M_{\odot }}$, Rin = 3.5 R⊙, and Rout = 270 R⊙ are continuum-normalized and shown over the wavelength range 2.2850−2.3200 μm, covering the CO Δv = 2 (2−0) absorption band. The observed spectrum after applying the radial velocity correction derived from the cross-correlation analysis is shown in blue. The inset displays the CCF as a function of radial velocity, with the dashed red line indicating the velocity corresponding to the maximum correlation.

The measured radial velocities were first corrected to the barycentric reference frame based on the time of observation, and subsequently transformed to the local standard of rest (LSR) frame. With the source’s Galactic coordinates and the V_LSR derived from the observed spectra, we estimated the kinematic distance using a Galactic rotation model following Persic et al. (1996), with updated Galactic parameters from Reid et al. (2014, see Wenger et al. 2018 for further details). The derived kinematic distances for each source are listed in Table 2. For several sources, both near and far kinematic distance solutions are feasible; however, V10 admits only the far distance solution.

3.3. Derived physical parameters

The model incorporates seven free parameters: M_∗Ṁ, $sini\sqrt{M_{\ast }}$, R_in, A_V, scale, R_out, and α. To derive the physical properties for each source, we fit the observed near-infrared spectra and the SPHEREx measurements with the corresponding model spectra. For a consistent comparison, we re-binned the observational and model spectra onto a common logarithmic wavelength grid, where the ratio between successive wavelength pixels is fixed at 1.00002. For the SPHEREx data, after applying the quality cuts described in Sect. 2, we evaluated the model prediction for each SPHEREx spectral element by averaging the model flux across the wavelength interval associated with that SPHEREx bandpass. This approach ensures that the comparison between the model and the SPHEREx observations correctly incorporates the finite spectral resolution of SPHEREx.

Parameter estimation is performed using a Markov chain Monte Carlo (MCMC) approach. The details of the likelihood function, priors, and sampling strategy are described in Appendix A. For all sources except V10, two possible kinematic distances (near and far) are available; therefore, we perform the fitting procedure independently for the two distance solutions. In the subsequent fitting, we excluded the far distance for V450 and the near distance for V978.

To validate our fitting methodology, we apply the same procedure to the well-studied FU Orionis. We use a high-resolution spectrum observed with IGRINS at the 4.3-m Discovery Channel Telescope (DCT) at Lowell Observatory on 19 January 2019 (Sawczynec et al. 2025). The spectrum is processed following the same steps described above. Since all SPHEREx data points for FU Orionis are flagged with values exceeding 4.29 × 10⁹, indicating unreliable spectrophotometry, we do not use SPHEREx measurements for this source. Instead, the flux calibration is performed using the telluric standard star observed during the same night, following the same procedure adopted for V978. In addition, FU Orionis falls within the masked region of the kinematic distance uncertainty model (160° < l < 200°), where kinematic distances are known to be highly uncertain. Consequently, we do not attempt a kinematic distance determination. Instead, we adopt the distance inferred from the parallax measurement reported by Gaia Collaboration (2021), as compiled in the SIMBAD database (Wenger et al. 2000), yielding a value of 408 pc for FU Ori. This estimate is consistent with recent literature determinations, such as those reported by Kounkel et al. (2018) and Roychowdhury & Hillenbrand (2024).

For illustration, we present the posterior distributions (corner plot) for source V10 in Fig. 4. V10 has only a far-distance solution, while the corner plots for the remaining sources and distance solutions are provided in Figs. B.1–B.4. The best-fitting model spectra for all sources are shown in Fig. 5, and zoomed-in views of selected spectral regions for the four sources and FU Orionis are shown in Figs. 6 and B.5–B.8. We note that for parameters with asymmetric or extended tails, the median and associated scatter should be interpreted with caution.

Fig. 4. Parameter distributions and corner plot of V10 assuming the far distance. The marginalized one-dimensional distributions are shown along the diagonal, where the dashed blue lines indicate the median and the 16th and 84th percentiles ±1σ. The two-dimensional panels show the joint posterior distributions, with contours enclosing 39%, 68%, 95%, and 99% credible regions.

Fig. 5. Comparison between the observed spectra and the corresponding best-fit model spectra for the sample objects. The gray lines display the observed near-infrared spectra in the H and K bands. The colored lines show the model spectra derived under different distance assumptions: the distance estimated from the parallax measurement compiled in the SIMBAD database (red), the near-distance solution (orange), and the far-distance solution (blue). Black circles with error bars mark the SPHEREx spectral data points included in the fit, whereas blue circles indicate SPHEREx measurements excluded from the fitting procedure because they were used for flux calibration. All spectra are plotted on a logarithmic flux scale.

Overall, the models reproduce the observed spectra – both the SPHEREx photometry and our high-resolution data – quite well, with good consistency in the absorption features and the general spectral shape for most objects. A notable exception is source V10, where the modeled mid-infrared flux is systematically lower than the SPHEREx measurements. This mismatch likely arises from several causes. First, the mid-infrared emission is largely produced by the passively heated outer disk, for which we assume a simplified geometric prescription, H/r = 0.2(r/R_out)^β, which may not accurately represent the true disk structure or temperature profile. Second, V10 may be embedded in a surrounding envelope whose warm dust adds mid-infrared emission that is not accounted for in our disk-only model. Third, additional processes such as small-scale surface inhomogeneities in the disk or localized heating events may also contribute. Given non-contemporaneous observations and the simplified nature of our passive disk treatment, we refrain from a detailed analysis of the origin of this mid-infrared flux discrepancy in the present work. We also conclude that the far-distance solution for V450 fails to adequately reproduce the observed spectrum. As further illustrated by the associated corner plot in Fig. B.2, this solution implies log₁₀(M_∗Ṁ) ≳ −3, which is markedly higher than the values typically derived for FUor-type objects, and is therefore rejected. In the case of V978, the near-distance solution produces model absorption features that are consistently weaker than those observed, and the corresponding accretion rate (log₁₀(M_∗Ṁ) ∼ −6.4 in Fig. B.3) is far below what is expected for FUor outbursts; this solution is likewise discarded. As a result, V191 is the only source for which both the near- and far-distance solutions remain compatible with the observational constraints. In summary, after removing the nonphysical or disfavored distance solutions, the final fitted physical parameters for all sources are reported in Table 2.

4. Discussion

4.1. Validation with FU Orionis

FU Orionis is one of the best-studied FUor systems and therefore serves as an ideal benchmark source for validating our spectral fitting framework. By applying the same fitting procedure to its near-infrared spectrum, we were able to assess the reliability and limitations of our method in comparison with well-established results from the literature (Hartmann et al. 2011; Carvalho et al. 2024b).

For the mass accretion rate, previous studies by Zhu et al. (2020, see their Fig. A.1) adopt a stellar mass of about 0.6 M_⊙ and infer an accretion rate of about 3.8 × 10⁻⁵ M_⊙ yr⁻¹, corresponding to log₁₀(M_∗Ṁ) 𢑘 −4.64. Our best-fit value, log₁₀(M_∗Ṁ) = −4.87, is lower by only about 0.1 dex. This small offset can be naturally explained by differences in flux calibration. In this work, the spectrum is flux-calibrated using the telluric standard star observed on the same night, so the calibrated flux reflects the instantaneous brightness of FU Orionis at the time of the spectroscopic observation. Converting our calibrated spectrum to equivalent 2MASS magnitudes yields H = 5.92 mag and K = 5.37 mag, which are fainter than the original 2MASS measurements (H = 5.70 mag, K = 5.16 mag). This behavior is consistent with the long-term photometric fading observed in the All Sky Automated Survey for SuperNovae (ASAS-SN) survey, indicating that the difference in log₁₀(M_∗Ṁ) primarily reflects intrinsic variability rather than a systematic bias in our modeling (Jayasinghe et al. 2018).

In contrast, parameters that are more sensitive to short-wavelength constraints show larger discrepancies. Our fitted inner disk radius (R_in ≈ 1.5 R_⊙) is smaller than the ∼3.5 R_⊙ reported by Zhu et al. (2020), and the derived extinction (A_V ≈ 2.87) is higher than the commonly adopted range of A_V ∼ 1.5 − 2.0 (Zhu et al. 2007; Carvalho et al. 2024b). These differences are most likely caused by the lack of optical photometric and spectroscopic data in our fitting. Extinction is only weakly constrained in the near-infrared, and R_in, which traces the hottest inner disk regions, is strongly coupled to the optical continuum level and line diagnostics. Without optical constraints, both parameters may therefore suffer from systematic uncertainties.

Adopting an inclination range of 37° −55° (Pérez et al. 2020; Zhu et al. 2007) and a stellar mass of 0.6 M_⊙, the expected range of $sini\sqrt{M_{\ast }}$ is 0.47−0.63, in excellent agreement with our fitted value of 0.55. In addition, the distance to FU Orionis is derived from parallax measurement compiled in the SIMBAD database rather than derived from the Galactic rotation curve, resulting in a relatively small systematic uncertainty in the scale parameter. The expected scale range of 0.57−0.80 is consistent with our best-fit value of 0.79.

FU Orionis lacks usable SPHEREx measurements, and we also examined contemporaneous NEOWISE photometry at the epoch of the spectroscopic observation. However, the NEOWISE measurements are severely inconsistent with the unTimely WISE data, most likely due to saturation effects. We therefore did not include mid-infrared constraints for this source, which may lead to unreliable estimates of the outer disk radius and viscosity parameter.

Overall, our benchmark test demonstrates that our spectral fitting approach robustly constrains key parameters such as log₁₀(M_∗Ṁ) and $sini\sqrt{M_{\ast }}$ using near-infrared spectra alone. Parameters such as R_in and A_V remain more susceptible to systematic uncertainties in the absence of optical data. In the case of FU Orionis specifically, R_out and α could not be reliably determined due to the lack of mid-infrared coverage, whereas for the other sources in our main sample, SPHEREx data provide sufficient constraints on these parameters.

4.2. Properties and spectral characteristics of our FUor candidates

For the accretion rates, we could only constrain the combination M_∗Ṁ from our spectral fitting. The four sources in our sample have log₁₀(M_∗Ṁ) values ranging from −5.1 to −3.7. Assuming a stellar mass of 0.5 M_⊙, the inferred accretion rates fall within the range typical for classical FUor sources, indicating that our discoveries are consistent with expected FUor accretion activity.

Regarding extinction, our sources exhibit significantly higher A_V values compared to FUor sources identified via optical surveys, with even the lowest-extinction object having A_V = 10 mag. This highlights the importance of infrared variability searches, which enable the identification of heavily obscured FUor systems that would otherwise be missed in optical surveys. Consistent with this picture, only one source in our sample (V10) has a Gaia counterpart (G = 20.72 mag), while the remaining sources lack Gaia detections, underscoring the severe optical obscuration affecting these objects (Gaia Collaboration 2023).

The H- and K-band spectra of our sources display characteristic FUor features, including CO Δv = 2 absorption, strong H₂O absorption, Brackett line absorption, and weak metallic lines such as Na I (2.209 μm in vacuum) and Ca I (2.261 μm). Comparison with the multi-temperature photosphere models shows good agreement across the H–K spectral range, capturing both line strengths and subtle continuum variations. While V10 and FU Orionis exhibit weak P Cygni profiles in Brγ (see Figs. 6 and B.8), which may indicate disk winds, confirmation requires additional J-band or optical spectroscopic data (see Carvalho et al. 2024a; Hillenbrand et al. 2023, 2025).

Fig. 6. Near-infrared H-band spectra (top panel) and K-band spectra (second panel) of V10 compared with best-fitting disk models. The gray curve shows the observed spectrum after telluric and flux calibration, while colored curves represent model spectra computed using different distance assumptions: red for the distance derived from parallax measurement compiled in the SIMBAD database, orange for the near kinematic distance, and blue for the far kinematic distance (when applicable). Vertical dashed lines mark prominent atomic transitions (gray), Brackett lines (red), and CO overtone band heads (green), as labeled. The shaded region in the H band indicates strong H2O absorption. Lower panels: Zoomed-in views of selected spectral lines in velocity space, centered on the rest wavelength of each transition. The vertical dashed line at zero velocity marks the line center, and nearby atomic and molecular transitions are indicated by vertical markers. Atomic lines are shown in red and orange to visually separate closely spaced transitions and their labels, while the molecular transitions (blue) shown in the first and eighth lower panels correspond to different rotational transitions within the same vibrational band. Observed spectra are shown in gray, and model spectra follow the same color coding as in the upper panels.

We additionally assessed the FUor classification of our targets using the equivalent width (EW) diagnostic introduced by Connelley & Reipurth (2018), which compares the summed EW of the Na+Ca lines to that of the CO v = 2 − 0 overtone band head. This EW diagram has the advantage of being largely independent of inclination and distance, as the selected lines arise from similar radii and are subject to similar veiling. Adopting the procedure of Liu et al. (2022), we first degraded our spectra to a resolution of 1200 prior to measuring the EWs. As illustrated in Fig. 7, all objects fall on the FUor side of the empirical boundary EW_CO = 3 × EW_Na + Ca + 5 Å, with the exception of V978, which lies close to the dividing line. We caution that our approach to determining EW_CO tends to yield slightly lower values than those reported in earlier studies (see Liu et al. 2022); therefore, we regard all sources as consistent with the FUor EW criterion.

Fig. 7. EW diagram of Na I + Ca I versus CO for the FUor template and the observed candidates. The dashed line represents the empirical criterion defined by Connelley & Reipurth (2018) to separate FUor-type objects from other YSOs. Objects below the dashed line show enhanced CO absorption relative to Na I and Ca I and are therefore consistent with FUor-like near-infrared spectra.

4.3. Confirmation rate of FUor candidates selected from infrared light curves

Li & Wang (2024) identified a total of 18 FUor candidates, together with one EXor candidate, based on distinctive infrared light-curve morphologies. To date, near-infrared spectra are available for 8 of these 18 FUor candidates. Among them, one source is classified as an EXor (V1298). Another object, V148, was previously characterized as an EXor (see Guo et al. 2021) but has recently been found to show a transition from CO emission to CO absorption, a hallmark of FUor-type spectra (Li et al., in prep.). In addition, V1707 is a rare object exhibiting strong aluminum monoxide (AlO) absorption while still retaining otherwise typical FUor spectral characteristics. The remaining sources – V10, V191, V450, V978, and V1349 – exhibit conventional FUor spectral features.

In our analysis, we confirm FUor-type spectra for six sources including four sources in this paper. Although V148 does exhibit CO absorption in our observations, our preliminary reduction suggests that the absorption lines are significantly broader than those observed in the FUor sources analyzed in this work. If interpreted within the FUor framework, this would imply a relatively large value of $sini\sqrt{M_{\ast }}$. Moreover, the infrared variability amplitude of V148 appears weaker compared to the confirmed FUor sources in our sample. Given these ambiguities, and in light of newly approved K-band spectroscopic observations, we refrain from drawing a definitive classification for V148 at this stage. A comprehensive analysis incorporating these new data will be carried out in a future study. We therefore consider the number of FUor-type spectra identified among the eight observed candidates to be k = 6–7, where k = 7 applies if V148 is ultimately confirmed as a FUor.

Assuming a parent sample of N = 18 FUor candidates, among which K objects truly exhibit FUor-type spectra, the probability of identifying k FUors in a randomly selected subsample of n = 8 objects follows a hypergeometric distribution,

$$\begin{array}{c}\hfill P_{A|B}(n,k,N,K)=\frac{C_K^k{C}_{N-K}^{n-k}}{C_N^n},\end{array}$$(15)

where event B denotes the existence of K FUor-type sources in the full sample and event A represents the identification of k FUor spectra in the observed subsample. Assuming a uniform prior probability for K, the posterior distribution is given by

$$\begin{array}{c}\hfill p_{B|A}(n,k,N,K)=\frac{p_{A|B}(n,k,N,K)p_B(K,N)}{p_A(n,k,N)},\end{array}$$(16)

with

$$\begin{array}{c}\hfill p_A(n,k,N)={\displaystyle \sum _{K=k}^N}{p}_{A|B}(n,k,N,K)p_B(K,N).\end{array}$$(17)

For N = 18, n = 8, and k = 6 − 7, we infer an expected value of K/N ≃ 0.72 − 0.83. The posterior probability that more than 12 objects in the parent sample exhibit FUor-type spectra reaches 63%−90%, supporting the conclusion that infrared variability is an efficient way to identify FUor systems.

5. Conclusions

In this work, we have presented near-infrared spectroscopic follow-up observations of four FUor candidates identified via infrared variability by Li & Wang (2024), using the IGRINS instrument on Gemini South. Our main conclusions are as follows:

1. Spectroscopic confirmation of FUor characteristics: All four sources (V10, V191, V450, and V978) exhibit classical FUor spectral signatures, including strong CO Δv = 2 absorption, pronounced H₂O absorption shaping the H-band continuum, Brackett line absorption, and weak metallic absorption features (Na I and Ca I), with little to no emission lines. Equivalent width diagnostics further support their FUor classification.

2. Physical parameters from spectral modeling: Using a combined viscous and passively heated disk model fitted to the near-infrared spectra and SPHEREx photometry, we derived key properties of the sources. The inferred accretion rates (Ṁ ∼ 10⁻⁵−10⁻⁴ M_⊙ yr⁻¹, assuming M_* ∼ 0.5 M_⊙) and high extinctions (A_V ≳ 10) are consistent with heavily embedded FUor systems, emphasizing the power of infrared surveys to uncover obscured outbursts missed by optical searches.

3. Validation of the fitting methodology: Application of the same modeling procedure to FU Orionis demonstrates that our approach reliably constrains fundamental parameters such as log₁₀(M_∗Ṁ) and $sini\sqrt{M_{\ast }}$ from near-infrared spectra alone, while parameters sensitive to optical or mid-infrared coverage (e.g., R_in, A_V, and R_out, α) can retain systematic uncertainties.

4. Efficiency of infrared variability selection: Of the eight FUor candidates with available near-infrared spectra, six to seven exhibit bona fide FUor-type spectra, implying a high success rate (72 − 83%) for selecting FUor outbursts based on infrared light curves. This highlights the effectiveness of mid-infrared variability surveys in identifying episodic accretion events, particularly for sources inaccessible to optical monitoring.

Overall, our study confirms that infrared time-domain selection, combined with high-resolution spectroscopic follow-up and detailed disk modeling, provides a robust framework for identifying and characterizing FUor outbursts in YSOs. Future systematic studies of infrared-selected FUor populations, with broader wavelength coverage, will particularly benefit from optical spectra to better constrain extinction and inner disk properties, as well as mid-infrared spectra to model the passive disk.

Acknowledgments

This research is supported by NSFC funding (NSFC-12393814). Based on observations obtained at the international Gemini Observatory, a program of NSF NOIRLab (Program ID: GS-2024A-Q-113), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation on behalf of the Gemini Observatory partnership: the U.S. National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This work used The Immersion Grating Infrared Spectrometer (IGRINS) was developed under a collaboration between the University of Texas at Austin and the Korea Astronomy and Space Science Institute (KASI) with the financial support of the US National Science Foundation under grants AST-1229522, AST-1702267 and AST-1908892, McDonald Observatory of the University of Texas at Austin, the Korean GMT Project of KASI, the Mt. Cuba Astronomical Foundation and Gemini Observatory. The RRISA is maintained by the IGRINS Team with support from McDonald Observatory of the University of Texas at Austin and the US National Science Foundation under grant AST-1908892. This publication also makes use of data products from the Spectro-Photometer for the History of the Universe, Epoch of Reionization and Ices Explorer (SPHEREx), which is a joint project of the Jet Propulsion Laboratory and the California Institute of Technology, and is funded by the National Aeronautics and Space Administration, and from NEOWISE, a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the Planetary Science Division of the National Aeronautics and Space Administration. In addition, this work makes use of data products from the Two Micron All Sky Survey (2MASS), a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. We also acknowledge the use of the SIMBAD database, operated at CDS, Strasbourg, France, and near-infrared photometry data provided by the VISTA VVV survey.

References

Appendix A: Fitting procedure

The fitting procedure is based on Bayesian inference. The observational data used in the fitting consist of two components: (i) the ground-based near-infrared spectra and (ii) the SPHEREx low-resolution spectroscopic measurements. The two datasets were jointly used to constrain the model parameters.

For a model with parameter θ_j, the likelihood of the combined dataset is defined as

$$\begin{array}{c}\hfill p(\mathrm{data}|{\mathit{\theta }}_j)=p({\mathrm{data}}_{\mathrm{spec}}|{\mathit{\theta }}_j)p({\mathrm{data}}_{\mathrm{SPHEREx}}|{\mathit{\theta }}_j),\end{array}$$(A.1)

where data_spec denotes the ground-based near-infrared spectra and data_SPHEREx denotes the SPHEREx spectroscopic measurements. Assuming independent Gaussian uncertainties, each likelihood term can be written as

$$\begin{array}{c}\hfill p({\mathrm{data}}_k|{\mathit{\theta }}_j)=\frac{1}{\sqrt{{(2\pi )}^{n_k}}{\displaystyle {\prod }_{i=1}^{n_k}}{\sigma }_{k,i}}exp[-\frac{1}{2}{\displaystyle \sum _{i=1}^{n_k}}{\left(\frac{F_{\mathrm{model},i}^k({\mathit{\theta }}_j)-F_i^k}{{\sigma }_{k,i}}\right)}^2],\end{array}$$(A.2)

where the subscript k refers to either the ground-based spectra or the SPHEREx data, $F_{\mathrm{model},i}^k$ is the model flux averaged over the corresponding wavelength bin, $F_i^k$ is the observed flux, and σ_k, i is the associated uncertainty. The total likelihood can therefore be expressed as

$$\begin{array}{c}\hfill p(\mathrm{data}|{\mathit{\theta }}_j)\propto exp(-\frac{1}{2}{\chi }_{\mathrm{tot}}^2),\end{array}$$(A.3)

with

$$\begin{array}{c}\hfill {\chi }_{\mathrm{tot}}^2={\chi }_{\mathrm{spec}}^2+{\chi }_{\mathrm{SPHEREx}}^2.\end{array}$$(A.4)

The corresponding log-likelihood function is then written as

$$\begin{array}{c}\hfill ln\mathcal{L}(\mathit{\theta })=-\frac{1}{2}{\chi }_{\mathrm{tot}}^2+\mathrm{const},\end{array}$$(A.5)

where the additive constant is irrelevant for parameter inference and is therefore omitted in the MCMC sampling with 32 chains (walkers). The posterior distributions of the model parameters are explored using a MCMC sampler. For each source, the chains are evolved for 20,000 steps, and the first 5,000 steps are discarded as burn-in. The convergence of the chains is assessed using the integrated autocorrelation time (τ), which is found to be in the range of about 10 − 20 for all parameters across all sources. For visualization, we used the corner.corner function from the Python corner package to generate the posterior distribution plots. A mild smoothing is applied with the parameters smooth = 1 and smooth1d = 1, corresponding to Gaussian smoothing of the two-dimensional and one-dimensional distributions, respectively. This smoothing is applied solely to suppress high-frequency sampling noise due to the finite length of the MCMC chains. The relatively small smoothing scale ensures that key features of the posterior distributions, such as peak locations and widths, are preserved.

Appendix B: Distribution, corner diagrams, and spectral fits

We present the parameter distributions, corner diagrams, and near-infrared spectral fits for all sources. The layout and notation follow those described in Figs. 4 and 6.

Fig. B.1. Distribution of parameters and corner diagram for V191. Top panel: Near solution. Bottom panel: Far solution. The symbols have the same meanings as in Fig. 4.

Fig. B.2. Same as Fig. B.1 but for V450.

Fig. B.3. Same as Fig. B.1 but for V978.

Fig. B.4. Same as Fig. B.1 but for FU Orionis.

Fig. B.5. Same as Fig. 6 but for V191.

Fig. B.6. Same as Fig. 6 but for V450.

Fig. B.7. Same as Fig. 6 but for V978.

Fig. B.8. Same as Fig. 6 but for FU Orionis.

All Tables

All Figures

	Fig. 1. Near-infrared and mid-infrared light curves of the four FUor candidates. We show the light curves from Ks (green), W1 (light green), and W2 (orange) bands (Minniti et al. 2010; Mainzer et al. 2014; Schlafly et al. 2019; Meisner et al. 2023). The vertical red and blue lines mark the dates of spectroscopic observations presented in this paper and Guo et al. (2024a).
In the text

	Fig. 2. Near-infrared spectra of four FUor candidates. The spectra are vertically shifted for clarity and have been convolved to a resolving power of R = 1200. Prominent atomic absorption features are marked by vertical dashed lines. We show the Brackett series (red), the CO band heads Δν = 2 (orange) and 3 (green), and the broad H2O absorption band (horizontal gray bar).
In the text

Fig. 3. Example of wavelength alignment of observed and template spectra for radial velocity determination. The observed spectrum (gray) and the viscous model spectrum (orange) with $M_{\ast }\dot{M}={10}^{-4.5}{M}_{\odot }^2$/yr, $sini\sqrt{M_{\ast }}=0.02\sqrt{M_{\odot }}$, Rin = 3.5 R⊙, and Rout = 270 R⊙ are continuum-normalized and shown over the wavelength range 2.2850−2.3200 μm, covering the CO Δv = 2 (2−0) absorption band. The observed spectrum after applying the radial velocity correction derived from the cross-correlation analysis is shown in blue. The inset displays the CCF as a function of radial velocity, with the dashed red line indicating the velocity corresponding to the maximum correlation.

In the text

	Fig. 4. Parameter distributions and corner plot of V10 assuming the far distance. The marginalized one-dimensional distributions are shown along the diagonal, where the dashed blue lines indicate the median and the 16th and 84th percentiles ±1σ. The two-dimensional panels show the joint posterior distributions, with contours enclosing 39%, 68%, 95%, and 99% credible regions.
In the text

Fig. 5. Comparison between the observed spectra and the corresponding best-fit model spectra for the sample objects. The gray lines display the observed near-infrared spectra in the H and K bands. The colored lines show the model spectra derived under different distance assumptions: the distance estimated from the parallax measurement compiled in the SIMBAD database (red), the near-distance solution (orange), and the far-distance solution (blue). Black circles with error bars mark the SPHEREx spectral data points included in the fit, whereas blue circles indicate SPHEREx measurements excluded from the fitting procedure because they were used for flux calibration. All spectra are plotted on a logarithmic flux scale.

In the text

Fig. 6. Near-infrared H-band spectra (top panel) and K-band spectra (second panel) of V10 compared with best-fitting disk models. The gray curve shows the observed spectrum after telluric and flux calibration, while colored curves represent model spectra computed using different distance assumptions: red for the distance derived from parallax measurement compiled in the SIMBAD database, orange for the near kinematic distance, and blue for the far kinematic distance (when applicable). Vertical dashed lines mark prominent atomic transitions (gray), Brackett lines (red), and CO overtone band heads (green), as labeled. The shaded region in the H band indicates strong H2O absorption. Lower panels: Zoomed-in views of selected spectral lines in velocity space, centered on the rest wavelength of each transition. The vertical dashed line at zero velocity marks the line center, and nearby atomic and molecular transitions are indicated by vertical markers. Atomic lines are shown in red and orange to visually separate closely spaced transitions and their labels, while the molecular transitions (blue) shown in the first and eighth lower panels correspond to different rotational transitions within the same vibrational band. Observed spectra are shown in gray, and model spectra follow the same color coding as in the upper panels.

In the text

	Fig. 7. EW diagram of Na I + Ca I versus CO for the FUor template and the observed candidates. The dashed line represents the empirical criterion defined by Connelley & Reipurth (2018) to separate FUor-type objects from other YSOs. Objects below the dashed line show enhanced CO absorption relative to Na I and Ca I and are therefore consistent with FUor-like near-infrared spectra.
In the text

	Fig. B.1. Distribution of parameters and corner diagram for V191. Top panel: Near solution. Bottom panel: Far solution. The symbols have the same meanings as in Fig. 4.
In the text

	Fig. B.2. Same as Fig. B.1 but for V450.
In the text

	Fig. B.3. Same as Fig. B.1 but for V978.
In the text

	Fig. B.4. Same as Fig. B.1 but for FU Orionis.
In the text

	Fig. B.5. Same as Fig. 6 but for V191.
In the text

	Fig. B.6. Same as Fig. 6 but for V450.
In the text

	Fig. B.7. Same as Fig. 6 but for V978.
In the text

	Fig. B.8. Same as Fig. 6 but for FU Orionis.
In the text