A Bayesian Approach to Inferring Accretion Signatures in Young Stellar Objects: A Case Study with VIRUS^*

Understanding how quickly a young star accretes mass is important for constraining models of circumstellar disk evolution and planet formation. Over the first few million years of their life, low-mass stars (≲2M_⊙) possess a circumstellar disk that is dissipated mainly through photoevaporation, stellar winds, and accretion onto the star (R. Alexander et al. 2014). One commonly accepted paradigm is that young stellar objects (YSOs) in the class II phase accrete matter via magnetospheric accretion (Y. Uchida & K. Shibata 1985; A. Koenigl 1991; F. Shu et al. 1994; L. Hartmann et al. 2016). In this model, the stellar magnetic field couples to material from the inner circumstellar disk, and then guides gas along magnetic field lines onto the stellar surface. The gravitational energy released from this free-falling gas consequently heats up the gas and causes broad emission lines. The infalling material produces a “shock” when it encounters the stellar surface, heating up patches of the surface to about 10⁴ K and creating UV continuum emission (E. Gullbring et al. 2000). Class III YSOs can still exhibit line emission like the younger accreting class II stars, but often to a lesser degree, and the emission is attributed mainly to the chromosphere rather than to accretion processes (C. F. Manara et al. 2013b).

Accretion directly “uses up” some of the circumstellar disk material. Moreover, the energy released from the accretion drives dissipation of material through outflows and photoevaporation of the disk. Finally, accretion influences the temperature and chemical composition of the inner disk. Through these channels, accretion therefore augments what material in the protoplanetary disk is available for planet formation. Deriving the mass accretion rate (M_acc) of YSOs is therefore key to understanding both low-mass star formation and consequentially how planets may form around the star (C. F. Manara et al. 2019). In particular, the M_acc–M_⋆ empirical relation from class II YSOs has been heavily studied in order to probe changes in mass accretion over time (e.g., C. J. Clarke & J. E. Pringle 2006) and the related evolution in disks (e.g., B. Ercolano et al. 2014; C. F. Manara 2018, and references within).

L_acc is the total accretion luminosity: the luminosity from the YSO that can be attributed purely to the accretion process. L_acc can be treated as the release of gravitational potential energy during accretion, and therefore be converted into a mass accretion rate, M_acc (see Equation (6) in Section 5.4). Accreting YSO spectra present continuum emission in excess of a main-sequence stellar photosphere, most notably through a “jump” in emission in the Balmer continuum (a “Balmer jump”) for wavelengths shortward of ∼3600 Å. The strength of this excess continuum emission can be used to estimate L_acc. Although not the only valuable wavelength range (for example, see the study of near-IR lines in E. Fiorellino et al. 2021; E. Fiorellino et al. 2023), the optical and near-UV wavelength range of a YSO spectrum therefore contains particularly useful information for estimating the total accretion luminosity L_acc. The typically high extinction of the YSO can make this a challenging regime for targeted observations, however, compared to the near-IR.

The “direct” method for determining L_acc involves fitting a model of the excess continuum emission to the YSO spectrum over a wide wavelength range. This direct method has been applied to numerous YSOs (e.g., J. A. Valenti et al. 1993; E. Gullbring et al. 1998; G. J. Herczeg & L. A. Hillenbrand 2008; E. Rigliaco et al. 2012; C. F. Manara et al. 2016; J. M. Alcalá et al. 2014, 2017). However, the techniques from previous literature have not been able to thoroughly examine degeneracy among model parameters. For example, the procedure described in C. F. Manara et al. (2013a) uses a discrete grid of parameter values and then finds the best fit by minimizing a likelihood function similar to a χ² distribution (see Equation (4)). While this approach self-consistently determines L_acc, stellar spectral type, and extinction A_V, it does not allow one to examine in detail the covariance of the parameters. Moreover, the uncertainty in L_acc has generally only been roughly estimated. For example in J. M. Alcalá et al. (2017), they state that the approximate uncertainty in L_acc is about ∼0.25 dex for log(L_acc/L_⊙), and they qualitatively describe the contributors to this uncertainty. However, they do not include a precise quantitative accounting of the uncertainty in L_acc for each individual YSO. We aim to address these deficiencies using the same model as in C. F. Manara et al. (2013a), but in a Bayesian framework, to our collection of YSO spectra in Section 4. While we study a relatively small sample of YSOs, our work acts as a proof of concept for how the analysis of YSO spectra can be improved, and a demonstration of the new nuts-for-ysos tool.

The accretion luminosity L_acc can also be indirectly estimated by measuring the luminosities of a variety of emission lines (L_line) and then applying known empirical relationships between L_line and L_acc. For example, the strength of the He i lines, Ca ii lines, and both the visible and near-infrared hydrogen recombination lines has been shown to exhibit a correlation with the accretion luminosities of YSOs derived from the direct method (e.g., J. Muzerolle et al. 1998; G. J. Herczeg & L. A. Hillenbrand 2008; E. Rigliaco et al. 2012). Therefore, determining accretion luminosities from such emission lines is a reliable alternative, having yielded results that are consistent with the direct method, albeit with a higher scatter (G. J. Herczeg & L. A. Hillenbrand 2008; E. Rigliaco et al. 2012; J. M. Alcalá et al. 2014). Additional proxies besides emission-line fluxes include excess U-band emission (E. Gullbring et al. 1998), and the H-alpha 10% line width (A. Natta et al. 2004). These various “indirect” methods are useful for YSOs in which directly fitting a model of excess continuum emission near the Balmer jump is too difficult because the spectrum is low-signal-to-noise (due to high extinction in star-forming regions, or low throughput in transmission optics of the spectrograph, for example) or because the spectrum is only available over a narrow wavelength range.

Using either the direct method or indirect method, accretion rates have been acquired for YSOs in various star-forming regions including Taurus (E. Gullbring et al. 1998; J. Muzerolle et al. 2005), TWA (L. Venuti et al. 2019), ρ Oph (A. Natta et al. 2006), Chamaeleon I (J. Muzerolle et al. 2005; C. F. Manara et al. 2016, 2017b), η-Chamaeleon (M. Rugel et al. 2018), Lupus (J. M. Alcalá et al. 2014, 2017), Upper Scorpius (C. F. Manara et al. 2020), and NGC 1333 (E. Fiorellino et al. 2021), among others. M_* and M_acc for these stars have been shown to exhibit a loosely positive correlation, of roughly M_acc. Some stars of the same mass exhibit a discrepancy in M_acc of more than 3 dex (E. Gullbring et al. 1998; J. Muzerolle et al. 2005; A. Natta et al. 2006; K. Biazzo et al. 2012; S. Antoniucci et al. 2014; J. M. Alcalá et al. 2017, and references therein). There is currently no single explanation for why there is a relationship between stellar mass and accretion rate, or why there is a large scatter in this relationship, although several possibilities have been put forward (e.g., L. Hartmann et al. 1998; A. Natta et al. 2006; C. P. Dullemond et al. 2006; E. I. Vorobyov & S. Basu 2008; B. Ercolano et al. 2014). Gathering statistics on a larger population of YSOs in the future may help to parse out the reasons for the scatter in the M_*–M_acc relation (R. Alexander et al. 2023).

Despite the value in examining the Balmer jump of YSO spectra, there has nonetheless been a lack of large spectroscopic surveys dedicated toward studying YSOs in the UV range. This is partly due to instrumental limitations, and also because the extinction within star-forming regions often prohibits a high enough UV signal. The entire UV spectrum has only been accessible from space telescopes, such as through the ULYSSES survey with the Hubble Space Telescope (HST; J. Roman-Duval et al. 2020; C. C. Espaillat et al. 2022; C. V. Pittman et al. 2022). However, very large ground-based instruments have recently provided access to wavelengths as small as ∼3000 Å such as X-Shooter (J. Vernet et al. 2011). We utilize the Hobby-Eberly Telescope (HET), a 10 m aperture telescope located at the McDonald Observatory in the Davis Mountains in Texas (L. W. Ramsey et al. 1998; G. J. Hill et al. 2021). HET is one of the largest optical telescopes in the world, and it feeds the three spectrographs housed by the observatory: Visible Integral-field Replicable Unit Spectrograph (VIRUS; G. J. Hill et al. 2021), the second-generation Low Resolution Spectrograph (LRS2; T. S. Chonis et al. 2016), and the Habitable-zone Planet Finder (HPF; S. Mahadevan et al. 2012, 2014). Spectra of a 3500–5500 Å wavelength range were obtained using VIRUS, a low resolving power (R ∼ 800) integral-field spectrograph with a 18’ diameter field of view (FOV). In this paper, we present the results of our analysis of VIRUS observations that were taken in parallel with the other spectrograph observations at the HET. This parallel data consists of VIRUS exposures that take place whenever another instrument aboard the HET is exposing for more than 5 minutes.

The nuts-for-ysos code automates the model-fitting process such that in the future, it could be applied to numerous stars with little human input. Especially since VIRUS is a multiobject spectrograph, future targeted observations of star-forming regions by VIRUS could provide a large collection of YSO spectra to be studied using our new approach. Moreover, the wavelength region and spectral features inputted into nuts-for-ysos can be customized to observations from other spectrographs. A similar analysis framework can even be used for spectra in the near-infrared range, for example, with future large area near-IR surveys done by the upcoming Nancy Grace Roman Space Telescope (D. Spergel et al. 2015; R. Akeson et al. 2019). Our study therefore offers an exciting prototype of what can be achieved with both VIRUS observations or with other YSO spectroscopic surveys in the future.

We describe the VIRUS instrument and the observations in Section 2, and the selection criteria for class II and class III YSOs in Section 3. In Section 4, we then describe the procedure in nuts-for-ysos for directly fitting an accretion model to our data, and how our procedure allows us to simultaneously derive spectral type, luminosity, and A_V for each target. We then interpret results from the fitting procedure, calculating L_*, L_acc, M_*, and M_acc in Section 5. In Section 6 we measure various emission lines for these YSO spectra, and apply the empirical L_line–L_acc relationship from J. M. Alcalá et al. (2017) to the targets, so that we can separately determine L_acc,line and compare it to L_acc. In Section 7 we then discuss covariances found within the results, and compare to results from Lupus, Chamaeleon I, and NGC 1333. Finally, in Section 8 we summarize our main points and reiterate the usefulness of nuts-for-ysos for future possible applications. Appendix A contains details about the data-reduction process and a table of observations (Table 10), and Appendix B contains plots of the accretion model fit to each YSO in our current sample, and corresponding corner plots for the model parameters (Figures 21–50).

At its full capacity, VIRUS has 78 identical fiber-fed integral-field units (IFUs) of 448 fibers (each 51″  ×  51″) arrayed in a grid pattern over a field of 18’ diameter. Each IFU is connected to two spectrograph channels each with CCDs read out through two amplifiers. Each of the ∼35,000 total fibers in VIRUS has a 15 diameter, with 22 between each fiber (1/3 fill factor). The resulting spectra have R = 670 at 3900 Å, R = 850 at 4600 Å, and R = 990 at 5200 Å (A. Kelz et al. 2006; G. J. Hill et al. 2018; G. J. Hill et al. 2021). VIRUS is primarily built for making extragalactic observations for the HET Dark Energy Experiment (HETDEX; K. Gebhardt et al. 2021). However, VIRUS continues to take data when LRS2 or HPF is making targeted observations for more than 5 minutes. These “parallel” VIRUS pointings occasionally coincide with the galactic plane, in which case they may capture YSOs within the 18’ VIRUS FOV. Parallel VIRUS observations are not dithered to fill in the fiber pattern, so the field has a fill factor of only ≈33%. Both the location and exposure time of each parallel observation are solely determined by the observations of whichever other spectrograph at the HET is being used for the primary science program at that time. The spectra extracted from these parallel observations between 2019 January 1 and 2023 March 31 have recently been released within the ongoing HET VIRUS Parallel Survey (G. R. Zeimann et al. 2024). The details of the data-reduction process are presented in Appendix A.

VIRUS made 4269 parallel observations between 2019 January 1 and 2021 May 30, with each of these observations typically containing hundreds to thousands of extracted spectra. We isolated potential YSO spectra from this large collection using the catalog presented in G. Marton et al. (2019). This catalog adopts machine learning methods to assign a probability of a star being a YSO (as well as main-sequence star, evolved star, or extragalactic object) based on Gaia DR2 (Gaia Collaboration et al. 2016, 2018) and on AllWISE photometry (E. L. Wright et al. 2019; R. M. Cutri et al. 2021). G. Marton et al. (2019) only examines regions where the two-dimensional, 353 GHz R2.01 Planck dust opacity map (Planck Collaboration et al. 2016) yields a value of at least τ = 1.3 × 10⁻⁵, to purposely focus on dusty regions more likely to contain YSOs; G. Marton et al. (2019) found that 99% of YSOs known from the literature occupy such regions with τ > 1.3 × 10⁻⁵. They also do not consider any objects with multiple Gaia IDs or missing photometric bands within Gaia DR2 or AllWISE. In building our sample, we first crossmatch VIRUS-observed objects within 2″ of G. Marton et al. (2019) catalog entries having an assigned YSO probability LY  > 70%. We then retain only the VIRUS spectra with an average signal-to-noise ratio (SNR) >5 over the entire spectrum, and an average SNR of at least 1.5 between 3500 and 4000 Å. This is because the signal covering the Balmer jump is especially important for the continuum fitting we perform in Section 4. We then crossmatch this set of observations with the catalog presented in C. A. L. Bailer-Jones et al. (2021) within 2″, to determine the distance to each object. Within this catalog, we use the “geometric distances” derived from GAIA EDR3 parallaxes⁸ (Gaia Collaboration et al. 2021).

This resulting sample has VIRUS spectra with sufficiently high SNR, known distances, and photometry available in the AllWISE catalog. However, it is important to note that stars on the Galactic plane can often possess spurious WISE photometry. The WISE mission was mainly designed for studying near-Earth asteroids, infrared galaxies, and brown dwarfs. The process for identifying sources from WISE images (described in their explanatory supplement ⁹ and in K. A. Marsh & T. H. Jarrett 2012) is not well suited for regions containing dust or clouds bright in the mid-infrared. As a result, stars can often have W1 and W2 band photometry associated with a legitimate point source, but spurious W3 and W4 band photometry. In identifying YSOs by their infrared excess, WISE photometry must be used with great caution.

X. P. Koenig & D. T. Leisawitz (2014) investigates this issue with WISE photometry, and finds that real detections in AllWISE can be well-isolated by making cuts on four specific parameters within each of the W1–W4 bands:

1. 1.  
   wBsnr: The SNR in the band B  ∈  [1,2,3,4]
2. 2.  
   wBrchi2: the reduced chi2 of the profile fit
3. 3.  
   wBm: the number of exposures over which a profile-fit flux measurement for the source could be performed
4. 4.  
   wBnm: the number of profile-fit flux measurements for which the source was detected with wBsnr  > 3

G. Marton et al. (2019) addressed the problem of spurious AllWISE photometry by adopting the same general strategy of examining these AllWISE catalog parameters. However, instead of making hard cuts, they take a probabilistic approach. They make a training sample of 500 real and 500 spurious sources from visually inspecting WISE W3 and W4 images. They then use the AllWISE parameters from this training sample and apply the Random Forest method to assign to other sources a probability of being real, R.

We find that although this approach generally works, many objects have a probability R that hovers around 0.5, such that requiring R ≥ 0.5 still results in apparent false positives. We decide to impose additional requirements on AllWISE photometry, which we find ultimately isolates true YSOs more effectively. Following X. P. Koenig & D. T. Leisawitz (2014), we require that for WISE bands W1-W4, the photometric uncertainty wBsigmpro be non-null, and require that the signal-to-noise wBsnr and reduced chi-squared wBrchi2 meet the following conditions:

Though this set of equations is very effective in suppressing contamination from fake AllWISE detections, X. P. Koenig & D. T. Leisawitz (2014) noted that it also eliminates around two-thirds of real detections in bands W3 and W4. X. P. Koenig & D. T. Leisawitz (2014) thus changed their criteria to raise the retrieval rate in the W3 and W4 bands, at the cost of allowing more fake detections. We decide that instead of adopting this changed criteria, we select sources that either:

1. have AllWISE photometry with non-null photometric uncertainty, R ≥ 0.5 from G. Marton et al. (2019), and meet the criteria presented in Equation (1), or

2. have existing Spitzer photometry (IRAC 1-4 and MIPS1 bands with non-null photometric uncertainty) that can be used as an alternative.

We have thus only used AllWISE photometry when catalog parameters meet a number of strict requirements, and opted to search for Spitzer photometry in cases where the requirements are not met. The Spitzer photometry is acquired using the VizieR Photometry viewer tool¹⁰ with a Python interface,¹¹ using a search radius of 2.”

We ultimately select 16 YSOs for our current sample. Nine YSOs have usable AllWISE photometry determined with our criteria, and 15 of them have published Spitzer photometry. We use photometry from the c2d (Cores to Planet-forming Disks) Spitzer Legacy Program when available (N. J. I. Evans et al. 2003; C2D Team 2020), and use alternatives otherwise, as listed in Table 1. To classify the YSOs in our sample, we use the 4-class system introduced by T. P. Greene et al. (1994), which classifies via the spectral index α; essentially the slope of the SED between ∼2 and 20 µm:

To estimate α, we fit a line by least-squares fitting to the available AllWISE and/or Spitzer photometry. We then classify α with the cutoffs presented in T. P. Greene et al. (1994):

Information related to the VIRUS observations of each star is presented in Table 10 in Appendix A. Table 1 contains the YSO likelihood LY and the AllWISE validity likelihood R obtained from G. Marton et al. (2019), the classification and region for each object, the photometry used in determining each classification¹² , and the name of the object in the SIMBAD database. We find 12 of the stars are class II, and four are class III. Among them, Object 10 (EM* LkHA 351) is classified as class II but its spectral index α very nearly places it in the class III category. Object 3 (2MASS J20580138+4345201) is also considered class II, but is almost in the “Flat” category. Though it has a class II SED, the spectrum of Object 16 (ATO J052.3580+31.4444) could not be fit by a model of an accreting YSO, nor did it exhibit the Balmer emission lines expected of an accreting YSO. These two issues are discussed in Sections 4 and 6, respectively. Ultimately, this object was not included in the results (i.e., in Tables 5 and 6). By visual inspection, none of the objects had SEDs that resembled transition disk YSOs.

The four class III objects are 6, 7, 12, and 13. Two of these objects, 12 and 13 (2MASS J03283651+3119289 and 2MASS J03292815+3116285), have Spitzer IRAC and MIPS1 photometry obtained from separate catalogs since their MIPS1 photometry could only be found exclusively in Table 2 of H. Y. A. Meng et al. (2017). For both these objects, the MIPS1 photometry is labeled as “Not Individually Detected” by H. Y. A. Meng et al. (2017). Despite the dubious nature of their MIPS1 detections, we nonetheless include these two objects because their spectra confirm that they are young stars and because the MIPS1 photometry is not used for anything beyond categorizing the YSOs as class III.

In Section 3 we identified a sample set of YSO spectra, for which we now develop a procedure to fit a multicomponent YSO model. The goal is to acquire accretion luminosities L_acc of our class II sample using the “direct” method (fitting a model of excess continuum emission to the spectrum) with a Bayesian approach. We have written our procedure into a publicly available Python code called nuts-for-ysos. We first outline the methodology of nuts-for-ysos in Sections 4.1, 4.2, and 4.3. We then discuss in Section 4.4 the specific requirements of nuts-for-ysos for inputted YSO spectra, and which parts of the analysis can be customized by the user.

4.1. The Model

4.2. Criteria Used in Fitting

4.3. Initializing the Bayesian Fit

We use the Python package PyMC (J. Salvatier et al. 2016) to implement a Bayesian process for fitting the total model to each target spectrum. Within PyMC, we use the No U-Turn Sampler (NUTS) with 16 chains, each having a length of 2000. We set target_accept = 0.99, so that the sampler takes very small steps. We used NUTS because for a system with significant parameter correlations such as ours, a Hamiltonian-based Monte Carlo algorithm like NUTS is more efficient than the Metropolis–Hastings algorithm or Gibbs sampling (M. D. Hoffman & A. Gelman 2014).

To reduce the burn-in time, a starting point for the sampler is determined using a least-squares optimization. In finding the starting point, we attempt to match the spectral features listed in Table 3 by minimizing a likelihood function similar to a χ² function. This function was first introduced in C. F. Manara et al. (2013a) and is as follows:

where f_obs is the feature from the VIRUS spectrum, σ_obs is the respective uncertainty, and f is the same feature predicted the composite model. After the starting point for each parameter is found, it is used as the initial point for the chains of the NUTS sampler.

4.3.1. Interpolating between Photospheric Templates

NUTS requires each parameter in the model to have a continuous range, because it is a gradient-based sampler. NUTS takes the gradient of the likelihood with respect to each parameter, to reach convergence faster than other sampling methods (M. D. Hoffman & A. Gelman 2014). Six out of the seven total model parameters are inherently continuous, but T_eff is not, because we rely only on a set of 23 discrete photospheric templates of various spectral types. We thus made T_eff a continuous parameter by having the sampler linearly interpolate between photospheric templates. Each of the class III photospheric templates comes from real observations rather than simulation, and so each has an (unspecified) uncertainty on its stellar distance provided in C. F. Manara et al. (2013b) or C. F. Manara et al. (2017a). This introduces scatter in the perceived relative brightnesses between templates. For the sampler to smoothly interpolate between templates, we decided to scale the flux of each template so the templates are at the same distance, and then slightly rescale both the templates and their respective luminosities L_*,phot so that the brightness from each template would smoothly and monotonically increase with T_eff. We do so by forcing the median fluxes from 4500 to 5500 Å of each template to follow a fourth-degree polynomial fitted to the original medians. Figure 1 shows these original medians plotted against T_eff, with the best-fit polynomial overlaid as a dashed line. These adjustments ultimately give the sampler an easier set of templates to interpolate between without fundamentally changing the nature of the results. Without this cleaning of the templates, we would not have been able to make T_eff a continuous parameter in the model and instead would have needed to make an ad hoc selection of a template before the rest of the model-fitting process.

A similar methodology has been used by R. A. B. Claes et al. (2024) to interpolate between their grid of 57 class III photospheric templates. They first normalized their class III spectra to factor out the individual distances to each star. They then computed the median fluxes of the spectra within multiple wavelength ranges. Next, they performed nonparametric local polynomial fits to the medians as a function of SpT. The smooth curve resulting from each fit then served as a basis for an interpolable set, such that the grid of templates could be made as well-sampled as needed. Within the polynomial fitting procedure, they also account for uncertainties in the median fluxes, the extinction, and spectral type of each template. They do this by repeating the polynomial fit over 1000 Monte Carlo iterations with the error terms sampled from Gaussian distributions. The result of this process is that each interpolated template has an uncertainty in the spectral flux. This inclusion of uncertainty within the interpolation is an important addition, which was not in our work. In our work, we simply incorporated uncertainties of A_V and T_eff directly into the Bayesian inference of the parameters. However, propagating uncertainty during the template interpolation stage could be particularly relevant when using the late-M spectral type class III templates, which generally have the lowest SNRs in the grid. The interpolation procedure of R. A. B. Claes et al. (2024) is available through their FRAPPE tool on Github.¹³ While our own interpolation did not use this approach, nuts-for-ysos has recently been updated with a new version that gives the NUTS sampler the ability to include uncertainties in the spectral flux of an inputted template grid, such as those acquired from FRAPPE.

Zoom In Zoom Out Reset image size

4.3.2. Priors

After making the adjustment so the T_eff parameter is continuous rather than discrete, we are then able to create priors for every component of the YSO model. A schematic diagram is presented in Figure 2 showing how each prior contributes to the model. The priors associated with the slab model are bounded within a mostly typical range of parameter space: T_slab ranges from 5000–11,000 K, n_e ranges from 10¹⁰–10¹⁶ cm⁻³, and τ₀ ranges from 0.01–5.0. The T_slab and τ₀ priors are uniform distributions, and n_e is a uniform distribution in logarithmic space. The priors for K_slab and K_phot are set as Half-Flat distributions, constraining them to be above zero. The A_V prior is a uniform distribution bounded between 0 and 10.

Zoom In Zoom Out Reset image size

Table 4. Priors Used

Name of Prior	Distribution Type	Lower Bound	Upper Bound
Tslab	Uniform	5000K	11000K
log(ne)	Uniform	1010cm−3	1016cm−3
τ0	Uniform	0.01	5.0
Kslab	Half-Flat	0	∞
Kphot	Half-Flat	0	∞
Teff	Uniform	2615 K	5550 K
Teff Uncert	Bounded Normal (μ = 0K, σ = 100K)	−200K	200K
AV	Uniform	0	10
Distance dobs	GGDa	0 pc	∞ pc
Phot. template AV	Half-Normal (σ = 0.167)	0	∞
Phot. template log(L*,phot/L⊙) Uncert	Bounded Normal (μ = 0, σ = 0.2)	−0.5	0.5

Note. ^aMedian, 16th percentile, and 84th percentile of generalized gamma distribution defined for each target using the distance values in Table 5.

Download table as:  ASCIITypeset image

The benefit of a Bayesian model for the fitting is that systemic uncertainties in the photospheric templates can be seamlessly integrated into the model. For example, the templates have an uncertainty in their spectral type of roughly 1 subclass for K-type stars, and half a subclass for M-type stars (C. F. Manara et al. 2013b; C. F. Manara et al. 2017a), which we translate into an uncertainty in T_eff of 100 K. The prior for T_eff by itself is set to be a Uniform prior, bound between 2615 and 5550 K (roughly SpTs M8.5 to G8), but a Normal prior centered on 0K with σ = 100 K is then added to this, to reflect the uncertainty in T_eff of the photospheric templates. However, the T_eff uncertainty can only extend to the extrema of available templates (SpTs M9 and G5), beyond which we are unable to generate a model. We therefore bound the T_eff uncertainty prior to be within ±200K.

Additionally, the photospheric templates are uncertain in their extinction; although the templates were selected by C. F. Manara et al. (2013b) and C. F. Manara et al. (2017a) to have A_V close to zero, it is more accurate to include a possibility of nonzero extinction. According to C. F. Manara et al. (2017a), all templates compared to BT-Settl models (F. Allard et al. 2011) have an A_V <  0.5 mag, with the possibility of A_V = 0 mag being within 3σ at most. We add a Half-Normal prior centered on 0 with σ = 0.5/3 to the uniform A_V prior, to represent this uncertainty on A_V.

We also include other uncertainties not directly related to a fitted parameter. The photospheric templates have an uncertainty in their luminosity L_*,phot, given by C. F. Manara et al. (2013b) to be 0.2 dex. We represent this uncertainty with a Normal prior in log space, centered on 0 with σ = 0.2. For the distance to each target YSO, the catalog associated with C. A. L. Bailer-Jones et al. (2021) provides the median, 16th percentile, and 84th percentile of a posterior probability distribution. C. A. L. Bailer-Jones et al. (2021) used a generalized gamma distribution (GGD) as a prior when computing the “geometric distances” we used. We thus choose to represent the distance to each target YSO as a GGD, defined using the median, 16th percentile, and 84th percentile. Table 4 summarizes all of the priors.

4.4. The nuts-for-ysos Workflow

The nuts-for-ysos package incorporates the Bayesian analysis described throughout Section 4. The code has been written to also be compatible with YSO spectra obtained by other instruments other than VIRUS. The inputted YSO spectrum can be of a higher resolution than that of VIRUS, but it must be of the same or lower resolution than the Class III templates, and must occupy a wavelength range covered by the Class III templates (for example, between 3300.0 Å and 10189.0 Å if using the UV and optical Class III templates from X-Shooter). The user can customize several aspects of the analysis. The grid of interpolable Class III templates can be changed from the default grid of templates listed in Table 2. For example, the new grid of templates presented in R. A. B. Claes et al. (2024) is better sampled in SpT and can be substituted within nuts-for-ysos. The only requirement for the templates themselves is that the wavelengths of the template spectra are no lower than 500.0 Å or higher than 25000.0 Å. As mentioned in Section 4.2, the features for which the model is evaluated can also be altered within nuts-for-ysos. The tool is by default capable of computing individual flux values, slopes, ratios, and photometric magnitudes from spectra. The number, types, and wavelength ranges of these features can be changed by the user as desired. As an example, we tested nuts-for-ysos with several low-resolution class II YSO spectra among the HST ULYSSES survey (J. Roman-Duval et al. 2020). ¹⁴ Specifically, we fit the model to data from the R ∼ 500 STIS G430L and G750L gratings. We did so between 3300Å and 6000Å, using the same Class III templates in Table 2 and using the same features as in Table 3, with three additional features included at ∼3310 Å, ∼5850 Å, and ∼6000 Å. We found the code runs successfully, and we show in Figure 3 an example fit for the Class II YSO Sz97. We performed compatibility tests of nuts-for-ysos with X-Shooter data of several Class II YSOs from J. M. Alcalá et al. (2014) and J. M. Alcalá et al. (2017) as well. Version 1.1 of nuts-for-ysos is available on GitHub¹⁵ and on Zenodo (L. H. Willett et al. 2024).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We check for convergence of each fit using the Gelman-Rubin statistic, by requiring ≤ 1.1 (A. Gelman & D. B. Rubin 1992) for every parameter. To mitigate autocorrelation, we only use every 20th value in each parameter posterior outputted by PyMC to compute further results. Since we ran 16 chains with length 2000 each, the original posterior has a length of 32,000, and the new thinned posterior has length of 1600.

The composite model can have a T_eff of at least 2615 K and at most 5550 K. We encountered two cases in which the fitted model to a target from our sample had an effective temperature on the edge of these boundaries. Object 13 (2MASS J03292815+3116285), a class III YSO, was found to have a T_eff posterior that peaks at ∼2600 K and is cut off for lower temperatures. Similarly, for Object 4 (EM* LkHA 188), a class II YSO, we found a T_eff posterior that peaks at ∼5500K and is cut off for higher temperatures. Because the sampler was unable to explore the entirety of the plausible model parameter space for these two objects, we take their T_eff to be upper and lower limits, respectively. The results for their other parameters are taken to only be approximate throughout further analysis.

Object 16 (ATO J052.3580+31.4444) was the only star for which we were unable to fit an accretion model, a finding that was briefly previewed at the end of the Sample Selection section (Section 3). This object has seemed to be a YSO in some ways, as it belongs to NGC1333 and has a class II SED. However, we found that our process was unable to plausibly fit an accretion model to the spectrum of Object 16. The spectrum also appears to have Balmer lines that are in absorption rather than emission, as discussed in Section 6. We were only able to fit a model of a reddened class III template having T_eff ≈  5500 K and A_V ≈  3.0. A plot of the data and model is shown in Figure 5.

Zoom In Zoom Out Reset image size

For the rest of the objects we were successful in fitting a model of an accreting YSO that converged and explored the full plausible range of T_eff. Plots of the median model fits for these 15 objects are shown in Figure 4. Figures 21–50 in Appendix B illustrate the full set of results for the 15 objects. This includes all median model fits, and underneath are also plots of probability distributions for each of the model parameters. For each object, a corner plot of the model parameters is also shown, as discussed in Section 5.2.

We find that in general, the least-constrained parameters tend to be those associated with the slab portion of the model. This is especially true of the electron density n_e, which has a posterior sometimes occupying a flat or nearly flat probability distribution in logarithmic space. The posterior probability distributions of slab parameters T_slab, n_e, and τ₀ also sometimes appear to peak at or near the bounds of their uniform priors; however, extending the bounds on these priors into more extreme territory was often found not to reveal any sudden dropoff in probability. Other model parameters K_phot, T_eff, and A_V are more constrained and tend to have posterior probability distributions that are more Gaussian in appearance. In Section 7 we discuss the posteriors and their correlations in more detail.

5.1. Discerning Upper Limits on Accretion

In objects that appear to have a minimal slab component in the model fit (e.g., Objects 8, 10, and 12, as seen in Figures 35, 39, and 43 of Appendix B), there is an additional technicality that needs to be considered. The K_slab posterior reaches very low values for these objects (for example, log(10⁶*K_slab) ∼ −5 for Object 12), but the lower bound for the uniform K_slab prior was set exactly to 0 in Section 4.3.2. The sampler has a small but nonnegligible step size, which keeps it from precisely reaching this K_slab = 0. We therefore must develop some way to check each object for whether the tapering of probability at low K_slab values is a legitimate property, or if it is caused by boundary behavior. We accomplish this by creating an altered K_slab prior that is uniform instead in logarithmic space. The lower bound of this prior is set to log(10⁶*K_slab) = −10, which is much lower than any K_slab resulting from the original model fits. If fitting the model with this new K_slab prior yields a posterior that goes all the way down to –10, then this means the original model fit was affected by the lower bound of the prior. In these cases, our original model fit should only be interpreted as an upper limit on accretion. The object could be accreting, but it also might not be accreting at all. In the opposite scenario, if the object is definitely accreting, the resulting K_slab will be biased toward lower values but will still taper off in probability at approximately the same value as the original result.

Using this criteria, we find that five YSOs have slab components that should only be interpreted as upper limits: the class III objects 6, 7, and 12, and class II objects 2 and 8. Figure 6 illustrates the difference between K_slab posteriors of a YSO that is definitely accreting, albeit relatively low (Object 10), and a YSO which has an upper limit only (Object 8). For Object 10, the posterior tapers at  ≈ −1 regardless of which prior is used. For Object 8, the posterior resulting from the altered prior extends all the way to the lower bound of log(10⁶*K_slab) = −10.

Zoom In Zoom Out Reset image size

5.2. Calculating Luminosity and Accretion Luminosity from Fit

The flux from accretion F_acc is derived simply by integrating the hydrogen spectrum portion of the model fit (the “slab” portion) over a very wide wavelength range; we integrate from 500 Å to 25000 Å. The integration is repeated for every possible slab created by the posterior. F_acc is then converted into an accretion luminosity using the relation for a distance d_obs from the target. The stellar luminosity of each target is derived using the class III template portion of the model fit (the photspheric portion). Every best-fit class III template has a precomputed luminosity L_*,phot as discussed in Section 4.3.1, which can then be combined with the determined template scaling, K_phot, to calculate the stellar luminosity L_* of the target. The formula is:

Since every parameter involved in these calculations occupies a probability distribution defined in advance or defined by posterior outputted by PyMC, we are essentially building probability distributions for L_* and L_acc.

Table 5 lists averages and standard deviations for the most physically relevant parameters of the sample: T_eff, A_V, log(L_*), and log(L_acc), as well as the distances to each of these targets (with their 16th and 84th percentiles included). Objects 4 and 13, for which T_eff is only a lower or upper limit, have all other model parameters listed as approximate. The objects with the lowest-SNR VIRUS spectra in Table 10 of Appendix A tend to have the largest uncertainties in T_eff and A_V, and the largest fractional uncertainty in K_phot. These uncertainties translate into wider error bars on log(L_*) and log(L_acc). For example, the class II Objects 1 and 3 (with average SNRs of 5.6 and 5.8, respectively) both have a standard deviation in log(L_acc) of σlog(L_acc) ∼ 0.5, larger than the σlog(L_acc) ∼ 0.15–0.30 found for other objects in the sample. In Table 5, we do not list numerical results for T_slab, n_e, and τ₀. As we briefly addressed in Section 4.1 (and discussed further in Section 7), the flux of the slab model is useful for determining L_acc, but the model itself is based on simplified physical assumptions, and so the numerical results for these three individual parameters should not be considered in great detail.

Table 5. Both Class II and Class III YSOs from the Sample, with Information on Their Distances Taken from C. A. L. Bailer-Jones et al. (2021), and Information on Temperature, Distance, Luminosity, and Accretion Luminosity Derived from the Model Fit to the Continuum

ID	Class	dobs (pc)	dobs,L (pc)	dobs,U (pc)	Mean Teff	σTeff	Mean AV	σAV	Mean log(L*)	σlog(L*)	Mean log(Lacc)	σlog(Lacc)
1	II	694.89	660.88	746.28	2971	220	2.7158	0.9156	−0.2056	0.3459	−1.2957	0.5477
2	II	734.37	720.62	746.51	4291	107	2.8836	0.1629	0.5163	0.1657	≤−1.5617
3	II	785.71	754.77	814.22	3968	387	3.4723	0.786	0.241	0.268	−0.1417	0.5203
4	II	746.32	737.9	753.86	≥5166		≈ 1.4186		≈ 0.5773		≈ 0.3249	0.2741
5	II	411.97	398.19	430.07	3761	116	3.6198	0.3877	−0.3555	0.2052	−1.3807	0.3005
6	III	430.93	424.05	438.52	4341	114	2.1424	0.1462	−0.3679	0.1564	≤−1.5028
7	III	347.26	340.54	352.53	5012	163	0.1269	0.0994	−0.0044	0.166	≤−1.6016
8	II	335.56	332.25	338.38	3582	98	0.1224	0.0904	−0.6987	0.1687	≤−3.1574
9	II	282	277.84	286.37	3060	93	0.1863	0.1282	−0.5797	0.1889	−2.7998	0.1467
10	II	285.52	283.18	287.5	3350	94	0.5525	0.1467	−0.2532	0.1578	−2.813	0.2962
11	II	303.16	297.88	307.58	3147	111	2.0873	0.399	−0.0704	0.2079	−1.2693	0.2532
12	III	281.19	275.61	286.59	3274	92	0.0824	0.0707	−1.0181	0.1653	≤−6.2905
13	III	297.2	285.6	308.59	≤2938		≈ 0.2679		≈ −1.1377		≈ −3.413
14	II	308.01	296.58	321.96	2828	114	2.1155	0.2844	−0.7414	0.2095	−2.2991	0.177
15	II	286.53	282.49	289.76	3676	106	1.0444	0.2114	−0.8576	0.1712	−2.542	0.2611

Download table as:  ASCIITypeset image

Appendix B (containing Figures 21–50) additionally includes corner plots of the model fit results for each object: both the parameters of the model itself, as well as log(L_*) and log(L_acc). We find several correlations between parameters within the corner plots, which are discussed in Section 7. Though all corner plots are in Appendix B, we also show an example for Object 11 in Figure 7. Figure 8 shows log(L_acc) plotted against log(L_*) for the 15 objects, where the error bars signify the standard deviations in log(L_acc) and in log(L_*). Downward arrows indicate an upper limit on log(L_acc) and open circles indicate only an approximate result for log(L_acc) and log(L_*).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.3. Estimating the Mass Posterior

We estimate masses for the class II and class III sample by interpolating pre-main-sequence (PMS) evolutionary models. Each object has a posterior outputted by PyMC for T_eff and L_*, which can be fed directly into an interpolation of evolutionary models to generate a probability distribution for the stellar mass M_⋆. We interpolate two sets of PMS models, from I. Baraffe et al. (2015) and L. Siess et al. (2000). The I. Baraffe et al. (2015) models are available for masses ≤1.4M_⊙ and ages ≥ 0.5 Myr. The L. Siess et al. (2000) models extends to younger ages and to a higher mass of 7M_⊙, but the lowest mass available is 0.1M_⊙. An HR diagram for all of the targets is plotted in Figures 9 and 10 with the I. Baraffe et al. (2015) and L. Siess et al. (2000) tracks, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the I. Baraffe et al. (2015) tracks, there are a number of objects (all class II) that lie above the youngest isochrone. We use a nearest interpolation for these objects, as opposed to linear interpolation for the rest. When determining which mass track a data point is nearest to on the HR diagram, we create a metric where 1 dex in log(L_*) scales to 1000K in T_eff. This scaling matches the ratio of the full ranges in log(L_*) and T_eff covered by the I. Baraffe et al. (2015) tracks.

Object 4 (EM* LkHA 188), a class II object with an SpT of roughly G8, was briefly discussed in Section 5 and considered to have only a lower limit on T_eff. It is too massive to lie on any track from I. Baraffe et al. (2015). For the tracks from L. Siess et al. (2000), it has an estimated mass of ≈2M_⊙. Object 13 (2MASS J03292815+3116285), a class III object with an SpT of roughly M8.5, was also briefly discussed in Section 5 and considered to have only an upper limit on T_eff. This T_eff is slightly too low-mass to lie on the 0.1M_⊙ track from L. Siess et al. (2000), whereas for the tracks from I. Baraffe et al. (2015), it has an estimated mass of ∼0.1M_⊙ or less.

Four of the objects lying above the youngest (0.5 Myr) I. Baraffe et al. (2015) isochrone (Objects 10, 11, 13, and 14) are from NGC 1333, which, in its entirety, is estimated to have an age of only <1 Myr (E. Fiorellino et al. 2021). NGC 1333 also includes Object 13, which also demonstrates a young age of (∼1Myr) according to I. Baraffe et al. (2015) models despite its class III SED. The other class III object from NGC 1333, Object 12, has an older estimated age of 3 Myr from I. Baraffe et al. (2015) and 6 Myr from L. Siess et al. (2000). On the other hand, Object 15 appears to be the oldest object from NGC 1333 in the sample (∼10 Myr), despite being a class II object, raising the possibility that it is subluminous. The rest of the class II sample has ages ≲10 Myr according to both the L. Siess et al. (2000) and I. Baraffe et al. (2015) models.

5.4. Estimating the Mass Accretion Rate Posterior

With the estimated accretion luminosities and radii from Table 5 and masses estimated with the I. Baraffe et al. (2015) and L. Siess et al. (2000) tracks, we calculate mass accretion rate M_acc with the formula

where R_⋆ is the stellar radius (calculated with ), and R_in is the inner radius of the circumstellar disk (E. Gullbring et al. 1998). R_in is assumed to be 5R_⋆ as has been done in a number of other studies (e.g., E. Gullbring et al. 1998; L. Hartmann et al. 1998; J. M. Alcalá et al. 2014; E. Fiorellino et al. 2021).

Masses and mass accretion rates (or their upper limits) are presented in Table 6 for both the I. Baraffe et al. (2015) and L. Siess et al. (2000) models. We find that for objects that have a mass that can be determined from both models, the masses and accretion rates from each model tend to lie within each others’ errors. Figures 11 and 12 show log(M_acc) plotted against log(M_*) for the 15 objects, derived using the I. Baraffe et al. (2015) and L. Siess et al. (2000) models, respectively.

Table 6. Mass and Mass Accretion Rate Information Derived from both the L. Siess et al. (2000) and I. Baraffe et al. (2015) Stellar Models

ID	Class	log(M⋆) (S00)	σ log(M⋆)	log(Macc)	σ log(Macc)	log(M⋆) (B15)	σ log(M⋆)	log(Macc)	σ log(Macc)
		(M⊙ yr−1)	(M⊙ yr−1)	(M⊙ yr−1)	(M⊙ yr−1)	(M⊙ yr−1)	(M⊙ yr−1)	(M⊙ yr−1)	(M⊙ yr−1)
1	II	−0.72	0.12	−7.68	0.62	−0.64	0.11	−7.57	0.65
2	II	0.02	0.07	≤−8.45		−0.02	0.06	≤−8.42
3	II	−0.16	0.17	−6.99	0.57	−0.18	0.15	−6.91	0.66
4	II	≈ 0.23		≈ −6.93
5	II	−0.30	0.07	−8.29	0.35	−0.30	0.08	−8.30	0.35
6	III	−0.02	0.04	≤−8.82		−0.03	0.04	≤−8.81
7	III	0.06	0.06	≤−8.94		0.07	0.06	≤−8.94
8	II	−0.41	0.06	≤−10.06		−0.39	0.08	≤−10.08
9	II	−0.73	0.08	−9.2	0.19	−0.75	0.06	−9.19	0.17
10	II	−0.53	0.05	−9.33	0.33	−0.54	0.06	−9.32	0.32
11	II	−0.64	0.06	−7.54	0.32	−0.56	0.08	−7.62	0.29
12	III	−0.65	0.08	≤−13.06		−0.62	0.09	≤−13.08
13	III					≈ −0.93		≈ −10.04
14	II	−0.93	0.09	−8.52	0.24	−0.87	0.1	−8.58	0.21
15	II	−0.36	0.07	−9.63	0.28	−0.3	0.08	−9.69	0.3

Download table as:  ASCIITypeset image

Having determined accretion rates for the YSOs using a model fit to the continuum, we now check that they agree with results independently derived from the emission lines. J. M. Alcalá et al. (2017) studied X-shooter spectra of 81 class II or transition disk YSOs in Lupus, using the “direct” method on each YSO by fitting the continuum model to each spectrum. They then updated empirical linear relationships between L_line and L_acc to calibrate the “indirect” method. The 3500–5500 Å wavelength range of VIRUS is suited for detecting numerous emission lines characteristic to YSOs, most notably the hydrogen Balmer lines beginning with the Hβ line (4860 Å) and higher-level lines in the series. For each class II YSO, we measured fluxes (or estimated upper limits) of emission lines within the VIRUS wavelength range for which in J. M. Alcalá et al. (2017) there are empirical relations given between line luminosity L_line and L_acc. This includes the hydrogen Balmer lines, from the Hβ line up to H15, the Ca ii K line, and various helium lines.

We developed a routine in Python separate from nuts-for-ysos, to measure these emission-line fluxes from the VIRUS spectra. For each spectrum, we first subtract an approximate continuum determined via a least-squares polynomial fitting function. For each individual emission line, we then estimate a baseline to match the flux of the local continuum with greater precision. In order to determine whether each emission line is detected in a spectrum, we use a “threshold finding” function from the package Specutils, which looks for deviations from the continuum above a given noise factor (Astropy-Specutils Development Team 2019). We deem a line “detected” if the function locates an emission feature above spectrum uncertainty by a factor of at least 2, and within 3 Å of the expected centroid. Finally, we measure the emission-line flux using an integrating function combined with a Monte Carlo procedure. We integrate the emission line 50 times with random noise added that is consistent with spectrum variance at these locations. The average of these 50 integrations is taken to be the final flux measurement, and the standard deviation becomes the associated error. If an emission line is not detected, we estimate 3σ upper limits on the flux using the same general approach as J. M. Alcalá et al. (2014). We take the upper limit to be 3 × F_noise × Δλ, where F_noise is the rms flux uncertainty over a 20 Å region centered on the expected line centroid, and the line width Δλ of the undetected line is assumed to be 5 Å at most. While this routine is not a part of the main nuts-for-ysos code, we do include it in a separate folder of the GitHub repository for the interested reader. There also exist plenty of other general line flux measurement tools in the community (eg. Specutils), and tools specific to YSOs such as the STAR-MELT package (J. Campbell-White et al. 2021).

The full list of emission lines used, and the corresponding extinction-corrected flux measurements for our current sample, can be found in Tables 7, 8, and 9. J. M. Alcalá et al. (2017) included L_line–L_acc relationships for the H and Ca ii H lines, but we did not attempt to measure fluxes for these emission lines. These two lines are partially blended in X-Shooter spectra, but fully blended with one another in our relatively low-resolution VIRUS spectra. Therefore, they could not be de-blended, as was done in J. M. Alcalá et al. (2017). Several line fluxes were omitted within Tables 7, 8, and 9 due to data-reduction-related issues over narrow wavelength regions. For example, Object 11 (2MASS J03285101+3118184) has no Hβ emission-line flux listed in Table 7. We also measured line fluxes for the class III YSOs, though we found that Objects 6 and Object 7 demonstrated very few detected emission lines. Object 16 (ATO J052.3580+31.4444), for which we were unable to fit a continuum accretion model to its spectrum, displays almost all of the applicable lines as either unmeasurable or in absorption rather than emission. The only exception is the Ca ii K line, which is in emission. We therefore did not measure line fluxes for Object 16.

We correct the emission-line fluxes for extinction using the A_V derived from the continuum-fitting process in Section 4. We use the same A_V so that the direct and indirect methods can be consistently compared. The extinction correction takes into account the uncertainty in A_V. This is accomplished by repeating the extinction correction over the entire posterior distribution of A_V values. We then take the mean to be the resulting extinction-corrected flux reported in Tables 7, 8, and 9, with the standard deviation taken to be the corresponding uncertainty.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The extinction-corrected emission-line fluxes are converted into line luminosities using the distance to each YSO. We propagate uncertainty in distance by using the entire posterior distribution in distance of each target (described in Section 4.3.2) for each calculation. In Figures 13 and 14, the resulting line luminosities and uncertainties are plotted against the continuum-derived L_acc for our sample. The empirical relationship from J. M. Alcalá et al. (2017) for each line is then plotted as a dotted red line. It is important to note that each line is not fitted to our data; it is completely independently derived by J. M. Alcalá et al. (2017) but shows remarkable agreement with our data. If we apply the relationships from J. M. Alcalá et al. (2017) to each L_line of our sample, we then acquire another estimate for the accretion luminosity. We then take the average of the estimates from each emission line, deriving a new quantity we call L_acc,line. We calculated L_acc,line for every object except for the class III Objects 6 and 7, since almost all of their line fluxes are only upper limits. Ultimately, we find that applying the updated relationships from J. M. Alcalá et al. (2017) to the line luminosities of our class II sample yields L_acc,line in good agreement with each L_acc derived in Section 5.2. Figure 15 shows L_acc,line plotted against the continuum-derived L_acc. This plot visually demonstrates an approximate 1:1 correlation, suggesting trustworthiness of our new Bayesian approach. The only exception is the class III Object 12, which has L_acc,line = 0.36 and an upper limit on L_acc≤ −6.28 that lies well outside of the plot limits. This very likely nonaccretor has significant emission lines (unlike Objects 6 and 7), but the lines could possibly be attributed mainly to chromospheric activity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Emission lines are produced by chromospheric activity in addition to accretion processes. This adds a bias to L_acc,line, which needs to be assessed. We check the L_acc,line of each star relative to L_* to evaluate whether any of the YSOs in our current sample can be considered “weak accretors” having emission-line fluxes dominated by chromospheric activity. Figure 16 shows the ratio log(L_acc,line/L_*) plotted against T_eff for each object. The red dashed line shows the level at which chromospheric noise is expected to become important, as derived by C. F. Manara et al. (2013b) using a set of class III YSOs. The orange dashed line shows the update for 4000K < T_eff < 5800K objects, made by C. F. Manara et al. (2017a). Finally, the blue dashed line shows the chromospheric noise level derived by R. A. B. Claes et al. (2024), with the T_eff scale converted from G. J. Herczeg & L. A. Hillenbrand (2014) to the one used in this work (K. L. Luhman et al. 2003 for M-type stars, and S. J. Kenyon & L. Hartmann 1995 for K-type stars).

Zoom In Zoom Out Reset image size

Among the two class III YSOs for which L_acc,line was calculated, Object 12 is the closest to the chromospheric threshold of C. F. Manara et al. (2013b), and is below that of R. A. B. Claes et al. (2024). This implies that the log(L_acc,line) value for this object is overestimated due to chromospheric line emission. Indeed, the upper limit of log(L_acc) found by fitting to the continuum is lower than the log(L_acc,line) value by approximately 2.5 dex. The two class II YSOs with upper limits on the L_acc derived in the continuum fit (Objects 2 and 8) also both demonstrate emission-line fluxes at or near to the chromospheric level. In Figure 16, Object 8 has a log(L_acc,line/L_*) that lies under the chromospheric threshold of R. A. B. Claes et al. (2024). Object 2 has a higher log(L_acc,line/L_*) lying slightly above the threshold. These objects are both plausible “weak accretors” by separate consideration of both their continuum emission and their emission lines.

7.1. Parameter Probabilities and Correlations

We have demonstrated an improved method to fitting an accretion model to YSO spectra by using a Bayesian framework. One main motivation in using a Bayesian approach is to understand hidden correlations among the parameters. From our small sample of YSOs, we have noticed several trends in behavior.

The quantity log(L_acc) often exhibits correlation with A_V. Uncertainties in A_V can therefore have a particularly strong influence on the determination of log(L_acc). The left-hand side of Figure 17 shows the overlaid posteriors of log(L_acc) versus A_V for every object in the class II sample for which log(L_acc) is not an upper limit. Within individual posteriors, a strong correlation is often evident between A_V and log(L_acc). Moreover, there is a scattered trend that the objects with a higher median log(L_acc) value tend to have the higher median A_V values. This degeneracy between the intrinsic luminosity of the accretion slab and the extinction could be partially improved by analyzing YSO spectra with a wider wavelength coverage. C. F. Manara et al. (2013a) pointed out that if one fits the model to a broader wavelength range (such as that of X-Shooter) or even just includes more photometric data from other bands in addition to the spectroscopic data, the reddening effect of extinction and the blue enhancement by the accretion slab may be better disentangled, although never completely. We also find that K_phot and A_V typically appear correlated. This makes sense, since K_phot and A_V both pertain directly to the scaling of the class III photospheric template and can be adjusted somewhat interchangeably, resulting in degeneracy. This can introduce a degeneracy of log(L_*) with A_V as well, although we find that the correlation between these two is less pronounced, as shown in the right-hand side of Figure 17.

Zoom In Zoom Out Reset image size

Examining the corner plots in Appendix B, there is often a notable lack of correlation between log(L_acc) and the parameters that make up the slab portion of the model: T_slab, n_e, and τ₀. This is true both within individual posteriors and also over the entire sample; for example, objects with a high median T_slab value do not demonstrate higher or lower median L_acc values than other objects. As previously noted in E. Rigliaco et al. (2012), there are degeneracies in the slab model, in which trade-offs of different slab parameters are able to produce similar accretion luminosities. E. Rigliaco et al. (2012) found that reasonable slab model fits to a YSO spectrum yielded L_acc estimates within only  ∼10% of each other, even over a wide variety of slab parameters. C. F. Manara (2014) included a full discussion of the interplay of these three parameters, and how they each affect the shape of the slab model. The Bayesian method shows us that these are nuisance parameters that can be marginalized over. Instead, only the overall scaling of the slab model K_slab tends to demonstrate some correlation with log(L_acc) in the corner plots, alongside the extinction A_V.

Within our current sample, we find only one instance in which the model fit shows clear preference for a certain T_slab value between the upper and lower bounds: Object 13 has a T_slab posterior that peaks at ∼7000K with a symmetric distribution. In other cases, we mainly find broad probability distributions for T_slab that are either relatively flat or skewed to low temperatures. The least-constrained parameter tends to be n_e. We find that n_e often does not display a strong peak in probability at any one specific value between the imposed lower and upper bounds of 10¹⁰ cm⁻³ and 10¹⁶ cm⁻³. In most cases, the distribution for log(n_e) gradually rises to peak probability at  ∼10 or ∼16, or remains nearly flat. For τ₀, we find only two instances where the model settles on a specific value between the lower and upper bounds of 0.01 and 5.0. These are Objects 11 and 14, for which the most probable τ₀ is  ∼3.5 and  ∼4.0, respectively. We more often find τ₀ skewed to the lowest or highest values possible, and otherwise a flat distribution or somewhat bimodal distribution between both extremes.

Interestingly, we notice that there is sometimes a correlation between T_slab and log(n_e) for T_slab values below  ∼10,000K. This can be seen most clearly in the corner plots for Objects 1, 5, 9, 10, 11, and 15. One possible reason for the correlation is the way that T_slab and n_e jointly affect the height of the Balmer jump (i.e., the ratio of flux at 3600 Å and 4000 Å). C. F. Manara (2014) showed that a specific value for the Balmer jump can be achieved for higher and higher T_slab, as long as n_e is also increased. However, the ratio gradually becomes independent of n_e at high T_slab values. A similar phenomenon occurred when C. F. Manara (2014) examined the ratio of Balmer continuum emission to Paschen continuum emission.

Three of the lowest-mass class II stars (Objects 1, 11, and 14) demonstrate both high T_slab and τ₀ values, and display an apparent anticorrelation between these two parameters for T_slab > 10,000 K. In their work, C. F. Manara (2014) showed how either increasing τ₀ or increasing T_slab can suppress the strength of features such as the height of the Balmer jump and the slope of the Balmer continuum. This trade-off may be one possible reason for the anticorrelation. Such a phenomenon may be particularly noticeable for these three objects because the Balmer jump is a fairly dominant component in their spectra, and because their T_slab and τ₀ posterior distributions happen to both prefer higher values.

C. F. Manara (2014) discussed the trustworthiness of the slab model for the sole purpose of determining L_acc. As we alluded to in Section 4.1, the slab model is based off an older paradigm for accreting YSOs, where particles cross a boundary layer between the inner circumstellar disk and the star (e.g., D. Lynden-Bell & J. E. Pringle 1974). This is different from the presently more widely accepted magnetospheric accretion model. Therefore, the three parameters T_slab, n_e, and τ₀ of the slab model should not be interpreted as physical quantities due to their lack of basis in reality. According to C. F. Manara (2014), however, there are a number of reasons why the slab model is still useful other than its simplicity. One reason is that more complex “shock models” based in magnetospheric accretion (e.g., that of N. Calvet & E. Gullbring 1998) have not been able to fully match the veiling of spectra for accreting class II stars at long wavelengths (L. Ingleby et al. 2013). However, the slab model has been able to match observations (e.g., J. A. Valenti et al. 1993; G. J. Herczeg & L. A. Hillenbrand 2008). Moreover, for smaller wavelengths of ≲3000 Å, using either model has yielded similar bolometric correction factors for the accretion luminosity L_acc (e.g., N. Calvet & E. Gullbring 1998; G. J. Herczeg & L. A. Hillenbrand 2008; C. F. Manara 2014). Therefore, the slab model, which is simple to implement with only a few parameters, is often still used. Nonetheless, it is sensible to treat the model with some caution given its dubious physical origins. Continued tests on YSO spectra at a variety of photometric wave bands can further verify the usefulness of the slab model for calculating L_acc.

7.2. Comparison to Specific SFRs

We have studied only a small sample of YSOs that do not all belong to one star-forming region. We therefore refrain from making physical interpretations of our results for M_acc. Instead, our study offers a prototype of what can be accomplished with optical/UVB spectroscopy of YSOs, combined with our Bayesian framework for fitting the accretion model, which can be applied to larger surveys in the future. With our small sample, we can at least note that our results for the L_*–L_acc plot (Figure 8) and M_*–M_acc plots (Figures 11 and 12) occupy a similar range as in previous studies of various star-forming regions. We briefly compare to previous results of class II stars from Lupus (J. M. Alcalá et al. 2017), Chamaeleon I (C. F. Manara et al. 2016), and NGC 1333 (E. Fiorellino et al. 2021).

Figure 18 shows that the L_* and L_acc of our class II sample follow a loose linear relationship in agreement with these three star-forming regions. A few objects in our current sample (Objects, 8, 9, and 10) have fairly low L_acc/L_* compared to most objects, but they are still certainly plausible, as several stars from Lupus and Chamaeleon I demonstrate similarly low L_acc/L_* values or upper limits.

Zoom In Zoom Out Reset image size

Figure 20 shows how the M_*–M_acc relationship from our current sample of 15 YSOs compares to Lupus, Chamaeleon I, and NGC1333. Past studies have demonstrated a loose correlation between M_* and M_acc. Given our small sample size, we do not attempt to fit a line to log(M_*)–log(M_acc) plot. However, we note that our sample, although scattered, occupies a space similar to these three star-forming regions. The highest accretion rates in our current sample are attained by objects 3 and 4 (respectively, log(M_acc) = −6.99 and log(M_acc)  ≈ −6.93 in M_⊙ yr⁻¹; via the L. Siess et al. 2000 evolutionary model). Although high, these accretion rates are similar to several objects from Chamaeleon I. We also observe a relatively high M_acc/M_* ratio for objects 1, 11, and 14. Although they lie along the high-M_acc border of the M_*–M_acc relationships for Lupus, Chamaeleon I, and NGC 1333, they are not implausibly high. Objects 1 and 11 are both from NGC 1333 and are within expectations when compared to similar-mass class II YSOs studied by E. Fiorellino et al. (2021).

There are two class II YSOs in the sample for which we derive only upper limits on mass accretion rate. Object 8 (CVSO 1897) demonstrates the lowest mass accretion rate among all of the class II stars, with an upper limit of log(M_acc) ≤–10.0 in M_⊙ yr⁻¹. This upper limit lies along the low-M_acc border of the M_*–M_acc plots found for Lupus, Chamaeleon I, and NGC 1333. As can be seen in Figure 19, Object 8 has a lower accretion rate than Object 10 (EM* LkHA 351), even though Object 10 has a lower spectral index α (being nearly in the class III category). Though Object 8 has a class II SED, its spectrum appears similar to a class III star. As shown in Figure 16 and previously discussed in Section 6, Object 8 has emission-line fluxes close to the expected chromospheric level for its T_eff, according to the thresholds defined by C. F. Manara et al. (2013b) and R. A. B. Claes et al. (2024). One other measurement for the accretion rate of Object 8 has been made by E. Manzo-Martìnez et al. (2020), and they found a log(M_acc) of −9.12 in M_⊙ yr⁻¹. They derived this accretion rate using the equivalent width of the Hα emission line. While their estimated M_acc is above ours, their estimates for A_V, L_* and M_* are similar. They found an A_V = 0.14 mag, log(L_*/L_⊙) = −0.80, and log(M_*/M_⊙) = −0.48, while we derive an A_V = 0.12  ±  0.09 mag, log(L_*/L_⊙) = −0.70  ±  0.17, and log(M_*/M_⊙) = −0.41  ±  0.06 using the L. Siess et al. (2000) models. The other upper limit in M_acc within our sample is Object 2, one of the highest-mass stars in the sample with M_* ≈ M_⊙. Object 2 demonstrates an upper limit of log(M_acc) ≤ −8.4. Similarly to Object 8, Object 2 also has emission-line fluxes close to the expected chromospheric level in Figure 16.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There is one other star in our sample, Object 14 (2MASS J03284782+3116552), for which we found a mass accretion rate measurement in previous literature. This object is present in the E. Fiorellino et al. (2021) study of NGC1333. E. Fiorellino et al. (2021) found that the star has log(L_*/L_⊙) = −0.72, in agreement our result of log(L_*/L_⊙) = −0.74  ±  0.2. Using the fluxes of the Paβ and Brγ lines and applying the empirical line luminosity relations from J. M. Alcalá et al. (2017), they found log(L_acc/L_⊙) = −1.95  ±  0.36. This is larger than the log(L_acc/L_⊙) = −2.30  ±  0.18 found in our study. Using the L. Siess et al. (2000) models, we estimate M_* = 0.12M_⊙, in close agreement to the M_* = 0.11M_⊙ mass from E. Fiorellino et al. (2021). However, we find a lower accretion rate of log(M_acc) = −8.81  ±  0.25 in M_⊙ yr⁻¹, as opposed to their log(M_acc) = −8.13. The disagreement could potentially be ascribed to the difference in A_V estimated. E. Fiorellino et al. (2021) estimated A_V = 4.0 mag by fitting the J and H bands of the class II spectrum with a class III template by hand. We estimate A_V = 2.1 though, using a similar process but in the near-UV and optical range. The disagreement may also be due to the underlying chromospheric activity, or differences in the choice of photospheric template between our work and E. Fiorellino et al. (2021). Besides Objects 8 and 14, there were no other literature values to compare with our sample, as far as we are aware.

The scatter in the log(M_*)−log(M_acc) relationship for previous works suggests that M_acc depends on more than just M_*. The aformentioned studies of Lupus, Chamaeleon I, and NGC 1333 each used samples of YSOs from single star-forming regions and, therefore, were all borne from the same local conditions at roughly the same time (although possibly with some spread in age; see the discussion in C. F. Manara et al. 2023 regarding age gradients within regions). As a result, scatter in the log(M_*)–log(M_acc) found by these studies can be more confidently attributed to individual differences in disk mass or the angular momentum of prenatal cores, or perhaps to time variability in accretion. On the other hand, the YSOs of our current sample have the additional complication of being from several SFRs rather than just one. Within singular star-forming regions, we do still note this definite diversity in M_acc. A prime example is the difference in M_acc between Objects 9 and 11. Both of these class II YSOs belong to NGC 1333 and have nearly the same mass (0.23 M_⊙ and 0.19 M_⊙, respectively; via L. Siess et al. 2000 models). However, they display M_acc that differ by  ≈1.5 dex. Object 11, the stronger accretor, has a spectrum that displays a visibly much more dramatic Balmer jump compared to Object 9.

We have developed a Python-based package called nuts-for-ysos and used it to study the spectra of 15 photometrically identified YSOs from VIRUS parallel data. With a wavelength range of 3500–5500 Å, VIRUS captures several important accretion indicators in YSOs, including the excess Balmer continuum and a variety of optical emission lines. Within nuts-for-ysos, we applied the No U-Turn Sampler (NUTS) to fit a model consisting of a theoretical accretion spectrum added to a class III YSO template, in order to replicate the continuum of our 15 VIRUS spectra. In doing so, we were able to simultaneously derive values for several stellar parameters and accretion parameters while examining uncertainties and covariances to an extent not achieved in the past. From these parameters, we were then able to determine mass accretion rates for each YSO. One main takeaway is the strong correlation between A_V and L_acc derived from the model fits. This highlights the importance of precisely determining A_V in order to constrain L_acc. This can better be achieved by fitting the model over a wide wavelength range, to differentiate the effects of extinction and accretion on the spectrum as much as possible. We then compared our results to those acquired using emission lines of the spectra, and found strong agreement. In comparing our study to previous studies of Lupus, Chamaeleon I, and NGC 1333, we also found our results to occupy a typical range in L_acc and M_acc.

Our results demonstrate the promise of applying a Bayesian framework like nuts-for-ysos for analyzing YSO spectra. The nuts-for-ysos code is flexible with regards to the number and types of spectral features it fits for, as well as the wavelengths of these features. It can therefore be applied to both future VIRUS observations or to observations from other spectrographs. Regardless, a future cohesive analysis of a larger collection of YSOs will be especially important for better characterizing the relationship between M_* and M_acc. This, in turn, could help elucidate the physics behind star formation and by association the formation of planets.

We thank the anonymous referee for providing insightful comments that informed edits to our manuscript. Observations were obtained with the Hobby-Eberly Telescope (HET), which is a joint project of the University of Texas at Austin, the Pennsylvania State University, Ludwig-Maximilians-Universität München, and Georg-August-Universität Göttingen. The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly.

VIRUS is a joint project of the University of Texas at Austin, Leibniz-Institut für Astrophysik Potsdam (AIP), Texas A&M University (TAMU), Max-Planck-Institut für Extraterrestrische Physik (MPE), Ludwig-Maximilians-Universität Muenchen, Pennsylvania State University, Institut für Astrophysik Göttingen, University of Oxford, and the Max-Planck-Institut für Astrophysik (MPA). In addition to Institutional support, VIRUS was partially funded by the National Science Foundation, the State of Texas, and generous support from private individuals and foundations.

The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing high-performance computing, visualization, and storage resources that have contributed to the research results reported within this paper.

Computations for this research were performed on the Pennsylvania State University’s Institute for Computational and Data Sciences’ Roar supercomputer. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University and the Eberly College of Science.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

This work is based in part on observations made with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.

This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in F. Ochsenbein et al. (2000). This research has also made use of the SIMBAD database, operated at CDS, Strasbourg, France.

We also thank the Pennsylvania State University Erickson Discovery Grant for providing support to this project during Summer 2021. Funding from this grant paved the path to complete a portion of the analysis.

Software: PyMC (J. Salvatier et al. 2016), PyTensor (PyMC Team 2023; R. Al-Rfou et al. 2016), Matplotlib (J. D. Hunter 2007), NumPy (C. R. Harris et al. 2020), SciPy (P. Virtanen et al. 2020), Astropy (Astropy Collaboration et al. 2013, 2018, 2022), Specutils (Astropy-Specutils Development Team 2019), PySpecKit (A. Ginsburg & J. Mirocha 2011).

The parallel observations were reduced using Remedy, a data processing pipeline for VIRUS.¹⁶ Remedy is responsible for calibrating the spectra from each individual fiber, and then combining these spectra into final “extracted spectra” of bright continuum sources. A full description of parallel observation reductions is presented in G. R. Zeimann et al. (2024), but we present relevant details to our study below. For each observation, Remedy first performs the gain multiplication, bias and dark frame subtractions, and masking of hot pixels identified in the master dark frame. The locations of the fibers on the detector are then determined using a master twilight frame compiled over several days. Wavelength calibration for each fiber is performed using Cd and Hg arc lamps. In order to extract complete spectra of continuum sources from the fiber spectra, a fiber normalization must be performed, which accounts for variations in CCD quantum efficiency and both fiber and spectrograph throughput variations. A master twilight frame is used to evaluate the relative throughput between each fiber and their average, and the fibers are normalized appropriately. A fiber is omitted from the extracted spectrum if its normalization factor is below 10%. For each exposure, a single sky model is constructed using identified “blank” fibers in the field. This first sky model is individually subtracted from every normalized fiber spectrum. Then a more localized residual sky model is subtracted, employing a Gaussian kernel in the fiber and wavelength direction. This kernel has a standard deviation of seven fibers in the fiber direction and 14 Å in the wavelength direction. Fibers with continuum emission >2σ are masked in this residual sky process. This entire sky subtraction process tends to fail in crowded fields with a lack of fibers pointed on blank sky. For example, several VIRUS parallel observations contained YSOs in nebulous regions such as the Orion Nebula. These observations were omitted from our study because the sky subtraction was unreliable. After sky subtraction, the VIRUS observations are then astrometrically calibrated using Pan-STARRS (Chambers & Pan-STARRS Team 2017) Data Release 2. Remedy matches point sources from each exposure with Pan-STARRS, and then shifts and rotates the VIRUS astrometry accordingly; shifts are usually <5” and rotations ≲01.

Remedy then extracts continuum sources from the normalized, sky-subtracted, and astrometrically calibrated fiber spectra. The seeing at the HET is usually around the same as the 15 diameter of the VIRUS fibers, meaning that the weight for each fiber in receiving light from a source can easily change with the presumed position of the source. Another complication is the presence of differential atmospheric refraction, which causes variation in the throughput with wavelength for a given fiber. Both of these phenomena are especially problematic for parallel VIRUS observations, which are not dithered, and therefore sparsely sample the sky compared to the usual 3-point dither pattern of VIRUS observations. Remedy uses the other stars in the VIRUS FOV (usually more than 20) to model the spatial point-spread function (PSF) and the differential atmospheric refraction, so that it can address these issues and accurately extract the spectrum of the source.

Each VIRUS parallel observation can ultimately contain hundreds to thousands of extracted spectra. Flux calibration is first performed on the extracted spectra using the throughput curve for the HETDEX collaboration (K. Gebhardt et al. 2012). Then, these preliminary calibrated spectra are convolved with the Pan-STARRS g filter. For each entire parallel observation, the overall biweight-estimated offset between the Pan-STARRS1 DR2 g magnitudes and VIRUS g magnitudes is then used to normalize this preliminary calibration. The finished flux calibration has a standard deviation of usually 0.1–0.15 mag. In this study, however, we chose to individually re-normalize extracted spectra using their Pan-STARRS1 DR2 “Mean PSF” g magnitudes (Chambers & Pan-STARRS Team 2017).¹⁷ This was because stars in star-forming regions may be young and variable in brightness, and these variations can affect a normalization done over an entire shot at once. Normalizing using photometry one at a time does not remove the possibility of each target’s variability affecting flux calibration, but it does allow for a more consistent approach not affected by the number and variability of neighboring stars. Moreover, we use the Pan-STARRS1 DR2 Mean PSF r and i magnitudes when fitting the accretion model to the VIRUS spectra in Section 4, so it is best that the g magnitudes of the VIRUS spectra are kept completely consistent with Pan-STARRS photometry. We are therefore assuming that the YSO spectrum has not drastically changed between the Pan-STARRS photometric observation and VIRUS observation. While this might be an unsuitable assumption for sources significantly varying in luminosity, like erupting FUors or EXors stars, none of our targets have previously been identified with these categories. The uncertainty in the Pan-STARRS photometry is taken into account when calculating the final uncertainty for each VIRUS spectrum. Table 10 lists the details of the VIRUS parallel observation corresponding to each object, as well as the name of the object in the SIMBAD database, the average SNR of the extracted spectrum, and the Pan-STARRS ID for each object.

Included in Figures 21–50 are plots of the median model fit for each object, and the probability distributions for each model parameter. Moreover, for each object there are corner plots of the model parameters and of log(L_*) and log(L_acc).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size