The VLA/ALMA Nascent Disk and Multiplicity (VANDAM) Survey of Orion Protostars. IV. Unveiling the Embedded Intermediate-Mass Protostar and Disk within OMC2-FIR3/HOPS-370

The ALMA observatory is located on the Chajnantor plateau in northern Chile at an elevation of ∼5000 m. HOPS-370 was observed with ALMA at 1.3 mm (Band 6) on 2018 January 7 with 43 antennas operating and sampling baselines from 15 to 2500 m. The observations were executed within an 87-minute observation block, and HOPS-370 was observed along with 19 other Orion protostars. The total time spent on HOPS-370 was ∼2.42 minutes, and the precipitable water vapor was ∼2.3 mm. The complex gain, bandpass, and absolute flux calibrator was J0510 + 1800. The absolute flux calibration accuracy is expected to be better than 10%. The correlator was configured with the first baseband containing a 1.875 GHz continuum band centered at 232.5 GHz and observed in Time Division Mode (TDM) with 128 channels; the remaining three basebands were configured in Frequency Division Mode (FDM). The second baseband was split into two 58.6 MHz spectral windows with 1960 channels each (0.083 km s⁻¹ velocity resolution) and centered on ¹³CO (J = 2 → 1) and C¹⁸O (J = 2 → 1). The third baseband was split into four 58.6 MHz spectral windows with 980 channels each (0.168 km s⁻¹ velocity resolution) and centered on SO (J_N  = 6₅ → 5₄), H₂CO (J = 3_0,3 → 2_0,2), and H₂CO (J = 3_2,2 → 2_2,1); the final window was centered between H₂CO (J = 3_2,1 → 2_2,0) and CH₃OH (J = 4_2,2 → 3_1,2), enabling both lines to be observed. Finally, the fourth baseband was configured with two 234 MHz spectral windows (980 channels, 0.367 km s⁻¹ resolution), one centered on ¹²CO (J = 2 → 1) and the other centered between N₂D⁺ (J = 3 → 2) and ¹³CS (J = 5 → 4).

The data were reduced using the ALMA calibration pipeline within CASA version 4.7.2. In order to increase the signal-to-noise ratio (S/N) of the continuum and spectral lines, we performed self-calibration on the continuum. We performed three rounds of phase-only self-calibration: the first round used solution intervals that encompassed the length of an entire on-source scan, then the second round utilized 12.1 s solution intervals, and the third round used a 6.05 s solution interval, corresponding to a single integration. The phase solutions from the continuum self-calibration were also applied to the spectral line bands. The resultant rms noise in the 1.3 mm continuum was ∼0.22 mJy beam⁻¹ and ∼12 mJy beam⁻¹ in 0.33 km s⁻¹ channels for the spectral line observations. The continuum and spectral line data were imaged using the clean task within CASA version 4.7.2. The continuum image was deconvolved using Briggs weighting and a robust parameter of 0.5, while the spectral line observations were deconvolved using natural weighting. The continuum image only uses uv-points >25 kλ to mitigate striping resulting from bright large-scale emission that is not well sampled. The typical beam sizes of the continuum and molecular line images are 023 × 013 (90 au × 51 au) and 032 × 018 (125 au × 71 au), respectively. The observational setups are also detailed in Table 1, and the reduced data are available from the Harvard Dataverse (Tobin 2020). ¹⁸

The ALMA 0.87 mm observations were taken as three executions of the scheduling block, with two executions on 2016 September 6 and the third on 2017 July 19. The time on source during each execution was ∼0.3 minutes, for a total of ∼0.9 minutes on HOPS-370 at 0.87 mm. The correlator was set up with two basebands observed in TDM mode, each with 1.875 GHz of bandwidth and 128 channels, centered at 333 and 344 GHz. The other two basebands were observed in FDM mode, centered on ¹²CO (J = 3 → 2) at 345.79599 GHz and ¹³CO (J = 3 → 2) at 330.58797 GHz. The bandwidth and spectral resolution of the spectral windows were 937.5 MHz with 0.489 km s⁻¹ channels and 234.375 MHz with 0.128 km s⁻¹ channels. Additional details of the data reduction and imaging are provided in Tobin et al. (2019, 2020). The data are available from the Harvard Dataverse (Tobin 2019a).

In addition to the continuum data, in this paper we make use of integrated intensity maps from the ¹²CO (J = 3 → 2) data cubes observed as part of this work toward HOPS-370 and HOPS-66. The ¹²CO data cubes were generated with 1 km s⁻¹ channels using the CASA 4.7.2 clean task with robust = 2 weighting, uv-distances >50 kλ to avoid artifacts from large-scale emission, and tapering at 500 kλ to increase sensitivity to extended emission. Masks were created manually through interactive execution of the clean task. The integrated intensity maps were generated using the CASA task immoments selecting the channel ranges where ¹²CO emission was detected. The reduced data cubes are also available from the Harvard Dataverse (Tobin 2019b).

The VLA is located on the Plains of San Agustin in New Mexico, USA, at an elevation of 2100 m. The VLA observations of HOPS-370 were conducted on 2016 October 26, while the VLA was in A configuration with 26 antennas operating. The entire observation lasted 2.5 hr, with ∼1 hr on source. We used the Ka-band receivers with 3 bit samplers, providing a 4 GHz baseband centered at 36.9 GHz (8.1 mm) and the other baseband centered at 29 GHz (1.05 cm). Additional details of the data reduction and imaging are provided in Tobin et al. (2019, 2020); the reduced data are also available for download through the Harvard Dataverse (Tobin 2019c).

Near-infrared spectroscopy was obtained from the Astrophysics Research Corporation (ARC) 3.5 m telescope at Apache Point Observatory in New Mexico, USA. HOPS-370 was observed on 2017 October 13 using the TripleSpec spectrograph (Wilson et al. 2004). TripleSpec simultaneously records spectra from ∼0.9 to 2.5 μm with a resolution of R ∼ 3000 with a 11 × 45'' slit.

The slit was centered on the base of the near-infrared scattered light nebula associated with HOPS-370, oriented in the east−west direction. This slit orientation minimized contamination from a near-infrared point source (MIR 22) located ∼3'' to the south. HOPS-370 was observed in an ABBA pattern, and the integration time was 2.5 minutes at each nod position, with a total on-source time of 30 minutes. Nodding was done along the slit, and the separation of nod positions was ∼20''. The average air mass during the observation was 1.5. The telluric standard used was the A0 star HD 37887 with a magnitude of 7.74 in Ks band, and it was observed in an ABBA pattern with 30 s in each nod position with a total time on source of 4 minutes.

The data were reduced using the IDL package Tspectool, which is a modified version of Spextool (Vacca et al. 2003; Cushing et al. 2004). Wavelength calibration was performed using OH sky emission lines. Flat-fielding was done using exposures of quartz lamps on the telescope truss. The flat field was constructed by subtracting exposures with the lamps on from exposures with the lamps off. Flux calibration of the spectrum was performed using the telluric standard with its cataloged magnitude relative to Vega. The reduced spectrum is available from the Harvard Dataverse (Tobin 2020).

The ALMA and VLA continuum images of HOPS-370 at 0.87, 1.3, and 9 mm are shown in Figure 2. The 0.87 and 9 mm images were previously presented in Tobin et al. (2019), and they are reproduced here to emphasize the clear disk emission. The 0.87 and 1.3 mm images clearly show the presence of a disk-like structure in the dust continuum emission that we will simply refer to as a disk. This disk is orthogonal to the outflow and jet directions traced by both Herschel (González-García et al. 2016) and the continuum between 5 cm and 9 mm (Osorio et al. 2017). The 5 cm contours from Osorio et al. (2017) are overlaid on the 1.3 mm map in Figure 3, further demonstrating the orthogonal relationship between the disk and the jet.

Zoom In Zoom Out Reset image size

To measure the geometric parameters of the continuum emission and the integrated flux densities, we fitted Gaussians to the images using the imfit task of CASA 4.7.2. The single-component Gaussians are not perfect fits to the data, but the integrated flux densities agree reasonably well with comparison to the flux density measured within a polygon surrounding the source.

The results from the Gaussian fits to the continuum emission at 0.87 mm find a deconvolved FWHM of 034 × 011 with a position angle of 109° (Table 2). Adopting the 2σ value of the Gaussian fit as the disk radius (Tobin et al. 2020), ¹⁹ we find a disk radius of ∼113 au, and the inclination of HOPS-370 can be estimated to be ∼71° by assuming that it is a geometrically thin disk and then calculating the inverse cosine of the deconvolved minor axis divided by the deconvolved major axis.

The VLA 9 mm continuum image in Figure 2 shows a distinctly different morphology with respect to the ALMA images. Rather than a simple disk feature, the VLA 9 mm image shows a cross-like morphology. The extension in the northeast-to-southwest direction is longer and more prominent than the extension in the southeast-to-northwest direction. The longer axis is orthogonal to the major axis of the disk and in the presumed direction of the jet/outflow from HOPS-370. We also show the spectral index map determined from the VLA 9 mm data alone. The spectral index of the emission (S_ν  ∝ ν^α ) indicates that the presumed jet is emitting via the free–free emission process, with values between −0.1 and ∼1 (e.g., Anglada et al. 1998). The spectral index in the direction of the disk is close to 2 at the edges of the spectral index map and more consistent with dust emission than free–free emission. Toward the central region, the spectral index is ∼1, too low to be dust, and suggesting a significant free–free contribution. Therefore, at 9 mm we trace both the jet and disk emission. The direction of the jet emission at 9 mm is consistent with the resolved jet reported by Osorio et al. (2017) at longer wavelengths (see Figure 3) and the molecular outflow in ¹²CO reported by Tobin et al. (2019).

We fit the VLA data with four Gaussian components simultaneously to describe the cross-like structure observed at 9 mm. The first component is a point-like component to match the peak intensity, the second and third components model the surface brightness distribution of the extended jet, and the fourth component fits the emission along the position angle of the disk. We add parameter estimates to help imfit arrive at a solution that can fit the disk component. The two components are needed for the emission along the jet axis because it cannot be well described with a single Gaussian. The parameters of each component are given in Table 3; note that the optimization of the four components yields one component with a negative flux. When this negative component is combined with the brighter positive component, the surface brightness distribution is best reproduced.

The extent of the emission from the disk at 9 mm is estimated to be 027 × 009 and a position angle of 112°, and the flux density of the disk is ∼0.73 mJy. The position angle of the disk is consistent with the fit at 0.87 mm, and the angular size at 9 mm is modestly smaller than at 0.87 mm, which is typical when comparing such short- and long-wavelength data (e.g., Segura-Cox et al. 2016; Tobin et al. 2020).

We plot the radio spectrum of HOPS-370 in Figure 4, including the flux densities reported in this work and those from Osorio et al. (2017) at longer wavelengths. We fit power laws to the dust emission between 0.87 and 9 mm and remove this estimated contribution from the longer-wavelength data points. Then, the longer-wavelength data points are assumed to only trace the free–free emission from HOPS-370. The flux density of the dust emission scales ∝ ν^{2.84 ± 0.08}, and the flux density of the free–free emission scales ∝ ν^{0.19 ± 0.08}. The spectral index of the free–free emission indicates partially optically thick emission since the flux density of optically thin free–free emission is expected to scale ∝ ν^−0.1. This is consistent with the spectral index derived at the position of HOPS-370 by Osorio et al. (2017), while at larger distances from the protostar they found that the jet exhibits nonthermal spectral indices.

Zoom In Zoom Out Reset image size

The mass of the disk traced by dust continuum emission can be estimated under the assumption of isothermal and optically thin emission using the equation

$$\begin{eqnarray}&&{M}_{\mathrm{dust}}=\displaystyle \frac{{D}^{2}{F}_{\nu }}{{\kappa }_{\nu }{B}_{\nu }({T}_{\mathrm{dust}})}.\end{eqnarray} \tag{ 1 }$$

D is the distance (∼392 pc), F_ν is the observed flux density, B_ν (T_dust) is the Planck function, T_dust is the dust temperature, and κ_ν is the dust opacity at the observed wavelength. We adopt κ_{1.3 mm} = 0.899 cm² g⁻¹ and κ_{0.87 mm} = 1.81 cm² g⁻¹ (Ossenkopf & Henning 1994), appropriate for protostellar envelopes. We adopt κ_{9.1 mm}= 0.13 cm² g⁻¹ by using κ_{1.3 mm} as a reference point and extrapolating it to 9 mm assuming a dust opacity power-law index of 1.0. This ad hoc extrapolation from 1.3 to 9 mm is necessary because the Ossenkopf & Henning (1994) models predict an opacity that is too low at 9 mm to yield consistent dust masses with shorter-wavelength observations (Segura-Cox et al. 2016; Tobin et al. 2016b; Tychoniec et al. 2018b). T_dust is often assumed to be 20 or 30 K (Jørgensen et al. 2009; Tobin et al. 2015; Tychoniec et al. 2018a) for solar-luminosity protostars, consistent with temperature estimates on ∼100 au scales (Whitney et al. 2003b). However, as part of the VANDAM Orion survey (Tobin et al. 2020) we used radiative transfer models to estimate the average dust temperature dependence on radius and luminosity. Using those results, we estimate a dust average temperature for a 100 au embedded protostellar disk of ∼31 K around a 1 L_⊙ protostar. Assuming that T_dust scales as (L/L_⊙)^0.25, the expected average temperature for the HOPS-370 disk is ∼131 K, with L_bol = 314 L_⊙ (Furlan et al. 2016). We finally multiply the resulting value of M_dust by 100, assuming the typical dust-to-gas mass ratio of 1:100 (Bohlin et al. 1978) to arrive at an estimate of the total disk mass.

The continuum flux density at 1.3 mm is 0.207 Jy, corresponding to a disk mass of 0.084 M_⊙. At 0.87 mm, the continuum flux density is 0.533 Jy, corresponding to a disk mass of 0.048 M_⊙. Lastly, at 9 mm the continuum flux density from the disk is 0.732 mJy, corresponding to a disk mass of 0.098 M_⊙. The modest discrepancy between 0.87 and 1.3 mm could result from the dust continuum at 0.87 mm being more opaque and/or uncertainty in the relative dust opacity between 0.87 and 1.3 mm. Also, the larger beam at 1.3 mm may enable more emission to be recovered owing to more short baselines being included in those observations. The 9 mm measurement agrees well with the others considering the uncertainty in dust opacity and the separation of its emission from the jet using multicomponent Gaussian fitting.

The continuum imaging from both ALMA and the VLA strongly indicates the presence of a disk surrounding HOPS-370. However, the nature of the disk structure can be analyzed further using the kinematics of molecular line emission. To this end, a suite of molecular lines were observed toward HOPS-370 with ALMA at 0.87 and 1.3 mm. However, the longer integration of the 1.3 mm data and shorter baseline coverage leads to those data being more sensitive to molecular line emission than the few lines covered in the 0.87 mm observations. The C¹⁸O, ¹³CO, H₂CO, SO, NS, and CH₃OH molecular lines all trace a signature of rotation from the disk revealed by the dust continuum. NS and H¹³CN are the only lines highlighted here from the 0.87 mm observations, but other complex organic molecules are also detected toward the disk at lower S/N (Tobin et al. 2019). Of the other targeted molecular lines, the ¹³CS line was only weakly detected, N₂D⁺ was not detected, and ¹²CO traces the outflow. The nondetection of N₂D⁺ is expected because of the warm temperature of the disk and inner envelope resulting from the high luminosity of the protostar (e.g., Emprechtinger et al. 2009; Tobin et al. 2013).

We show the integrated intensity maps of the blue- and redshifted emission of each molecular line (except ¹²CO, N₂D⁺, and ¹³CS) in Figure 5. The integrated intensity maps were created using selected channel ranges where there is spectral line emission using the CASA task immoments. We selected the channels corresponding to blue- and redshifted emission using the system velocity of ∼11.0 km s⁻¹ (Tobin et al. 2019); the blueshifted maps are integrated between 4 and 11 km s⁻¹, and the redshifted maps are integrated between 11 and 18 km s⁻¹. These blue- and redshifted integrated intensity images are used to assess the kinematics traced by the particular molecular lines. We do not show the ¹²CO moment maps because they were presented in Tobin et al. (2019). We also do not show N₂D⁺ because it is not detected, and ¹³CS is not shown because it is only marginally detected.

Zoom In Zoom Out Reset image size

The consistency of the velocity gradient in all well-detected molecular lines indicates that the disk around HOPS-370 is clearly rotating and may be rotationally supported. We note that the peaks of the line emission can be northeast and/or southwest with respect to the center of the continuum position, and these peaks tend to avoid the midplane traced in dust continuum. Continuum opacity and/or line opacity are likely to cause this feature; molecular freeze-out is not likely because of the warm disk temperatures due to the high luminosity of the protostar.

The integrated intensity maps only show limited spectral information; therefore, we also show position–velocity (PV) diagrams for each molecular line in Figure 6. The PV diagrams were extracted from the molecular line data cubes using a custom Python code to obtain a complementary view of the kinematic structure. The PV diagrams are extracted along a 06-wide strip (15 pixels) along the major axis of the disk, at a position angle 100° east of north. The emission is summed within the 06 strip to produce a two-dimensional spectrum. The PV diagrams bear strong resemblance to other protostellar disks that have been detected and characterized (Tobin et al. 2012; Sakai et al. 2014; Oya et al. 2015; Aso et al. 2017). The PV diagrams bear the signature of a rotating Keplerian disk with a finite radius: a linear velocity gradient transitioning from blue- to redshifted on opposite sides of the protostar and higher-velocity emission from spatial scales closer to the central protostar following a Keplerian velocity profile.

Zoom In Zoom Out Reset image size

These features are produced because the outer radius of the disk causes the emission from the largest radii of the disk to have a similar velocity profile to a rotating ring; in this case the ring is the outer disk. The finite radius of the outer disk means that Keplerian rotation will not continue to spatial scales that extend beyond the disk. Then, radii closer to the protostar rotate more quickly, producing the higher-velocity emission that is typically associated with Keplerian rotation toward embedded protostars. The linear transition from blue to red is most easily seen in the PV diagrams for NS and SO molecules. However, the transition is not always easily detectable toward protostellar disks (e.g., Tobin et al. 2012; Murillo & Lai 2013; Ohashi et al. 2014) owing to spatial filtering and blending with the infalling envelope and the molecular cloud near the system velocity. But certain molecular transitions that specifically trace the disk and not the surrounding cloud/envelope show the low-velocity emission of the disk with higher fidelity. This is due to a variety of possible reasons, but mainly critical density, chemistry, and excitation temperature. The lines SO, CH₃OH, NS, and H₂CO (J = 3_2,2 → 2_2,1), (J = 3_2,1 → 2_2,0) with T_ex ∼ 60 K best trace the disk of HOPS-370.

Previous studies have been able to use point-like line emission in velocity channels away from the system velocity toward disks around lower-mass Class 0 and Class I protostars to fit the emission centroids from PV diagrams, data cubes, and/or the uv-data themselves to map the rotation curves in one dimension. In these cases, the emission centroids systematically changed with each velocity channel, high velocities near the continuum source and lower velocities centered farther from the continuum source (Tobin et al. 2012; Yen et al. 2013; Ohashi et al. 2014). However, the clear non-Gaussian emission in the images prevents such analyses from being viable for HOPS-370. Thus, the rotation curve must be examined using modeling of the molecular line emission, which will be discussed in Section 5.

The near-infrared spectrum of HOPS-370 is shown in Figure 7 from 2 to 2.4 μm. HOPS-370 is also detected shortward of 2 μm but with lower S/N and is not shown. The main spectral features are the prominent H₂ emission lines detected throughout the band. We also detect weaker emission from the Brackett γ (Brγ) atomic hydrogen recombination line and CO band head emission. Continuum emission is detected, but there are no obvious photospheric absorption features detected in this medium-resolution spectrum. The H₂ emission lines are likely associated with shocks in the outflow from HOPS-370 given that most of the near-infrared emission detected in this spectrum is from the outflow cavity owing to the central protostar being too highly extincted. The properties of the Brγ emission are of the greatest interest to extract owing to their frequent association with accretion processes in young stars (e.g., Muzerolle et al. 1998; Connelley et al. 2009). The equivalent width of the Brγ line is −1.42 Å, with an integrated flux of 2.9 × 10⁻¹⁵ erg s⁻¹ cm⁻². This line flux translates to a Brγ luminosity of 5.6 × 10²⁸ erg s⁻¹ or 1.5 × 10⁻⁵ L_⊙. The equivalent width and line luminosity of the Brγ line are within the ranges typically observed toward Class I protostars by Connelley & Greene (2010). However, these values are on the low end for a protostar that is expected to be accreting rapidly. The interpretation of the line emission will be further discussed in Section 6.3.

Zoom In Zoom Out Reset image size

The molecular line images and PV diagrams show strong rotation signatures (Figures 5 and 6). To quantitatively determine whether the rotation is tracing a Keplerian disk and, if so, to measure the protostar mass, we must make use of radiative transfer modeling to fully utilize the constraints offered by the channel maps for multiple molecular lines. The pdspy package has distinct modes for fitting molecular line kinematic data and continuum data. The basis for both models is an analytic physical model for a protostellar system with a surrounding disk embedded within an infalling envelope.

5.1.1. Physical Model

The disk structure for the molecular line models uses an exponentially tapered disk density profile (Lynden-Bell & Pringle 1974). The exponentially tapered density profile is described by

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{0}{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{-\gamma }\,\exp \left[-{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{(2-\gamma )}\right],\end{eqnarray} \tag{ 2 }$$

where r is the disk radius defined in cylindrical coordinates. The power-law index of the surface density profile is defined by γ. The normalization constant, Σ₀, is given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}=\displaystyle \frac{(2-\gamma ){M}_{d}}{2\pi {r}_{c}^{2}},\end{eqnarray} \tag{ 3 }$$

where M_d is the disk mass and the other parameters are as defined above. The vertical structure of the disk is set by hydrostatic equilibrium, assuming that the disk is vertically isothermal such that

$$\begin{eqnarray}&&h(r)={\left(\displaystyle \frac{{k}_{{\rm{b}}}{r}^{3}{T}_{g}(r)}{G{M}_{* }{\mu }_{{\rm{m}}}{m}_{{\rm{H}}}}\right)}^{1/2},\end{eqnarray} \tag{ 4 }$$

where k_b is the Boltzmann constant, G is the gravitation constant, M_* is the central stellar mass, μ_m is the mean molecular weight of 2.37 (Lodders 2003), m_H is the mass of a hydrogen atom, and T_g (r) is the temperature profile defined as

$$\begin{eqnarray}&&{T}_{g}(r)={T}_{0}{\left(\frac{r}{1\mathrm{au}}\right)}^{-q},\end{eqnarray} \tag{ 5 }$$

where both T₀ and q can be free parameters.

The rotating, infalling envelope is described by the density profile of a rotating collapse model (Ulrich 1976; Cassen & Moosman 1981; Terebey et al. 1984). The envelope density profile is defined as

$$\begin{eqnarray}&&\rho =\frac{\dot{{M}_{\mathrm{env}}}}{4\pi }{\left({{GM}}_{* }{r}^{3}\right)}^{-\tfrac{1}{2}}{\left(1+\frac{\mu }{{\mu }_{0}}\right)}^{-\tfrac{1}{2}}{\left(\frac{\mu }{{\mu }_{0}}+2{\mu }_{0}^{2}\frac{{R}_{c}}{r}\right)}^{-1},\end{eqnarray} \tag{ 6 }$$

where ${\dot{M}}_{\mathrm{env}}$ is the mass infall rate of the envelope onto the disk, R_c is the centrifugal radius where the infalling material has enough angular momentum to orbit the star, μ = cos θ, and μ₀ is the cosine polar angle of a streamline out to r → ∞. The density profile inside of R_c will be ρ_env ∝ r^−1/2, and outside R_c it will be ρ_env ∝ r^−3/2. For the sake of our modeling, we consider R_c to be equivalent to the outer disk radius for a truncated disk model and the critical radius for an exponentially tapered disk. Thus, R_c  = r_c from Equations (2)–(4).

The envelope model includes outflow cavities with reduced envelope density. The widths of the outflow cavities are parameterized as

$$\begin{eqnarray}&&z\gt 1\,\mathrm{au}+{r}^{\xi },\end{eqnarray} \tag{ 7 }$$

and the envelope density is reduced by a factor of f_cav. The outflow cavity opening angle will be less than 45° for ξ < 1 and greater than 45° for ξ > 1. The outflow cavity opening angle, ψ, can be directly calculated from

$$\begin{eqnarray}&&\psi =2{\tan }^{-1}(\xi ).\end{eqnarray} \tag{ 8 }$$

The velocity profile of the infalling envelope is also adopted from the rotating collapse model (Ulrich 1976), where

$$\begin{eqnarray}&&{v}_{r}(r,\theta )=-{\left(\displaystyle \frac{{{GM}}_{* }}{r}\right)}^{1/2}{\left(1+\displaystyle \frac{\cos \theta }{\cos {\theta }_{0}}\right)}^{1/2},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{v}_{\theta }(r,\theta )=-{\left(\displaystyle \frac{{{GM}}_{* }}{r}\right)}^{1/2}\,(\cos {\theta }_{0}-\cos \theta ){\left(\displaystyle \frac{\cos {\theta }_{0}+\cos \theta }{\cos {\theta }_{0}{\sin }^{2}\theta }\right)}^{1/2},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{v}_{\phi }(r,\theta )=-{\left(\displaystyle \frac{{{GM}}_{* }}{r}\right)}^{1/2}\left(\displaystyle \frac{\sin {\theta }_{0}}{\sin \theta }\right){\left(1-\displaystyle \frac{\cos \theta }{\cos {\theta }_{0}}\right)}^{1/2}.\end{eqnarray} \tag{ 11 }$$

The velocity profile of this equation results in rotational velocity in the equatorial plane of the envelope that is equivalent to the Keplerian orbital velocity at a radius of R_c . As such, material located within the disk at radii smaller than R_c has velocities described by Keplerian rotation with

$$\begin{eqnarray}&&{v}_{\phi }(r)={\left(\displaystyle \frac{{{GM}}_{* }}{r}\right)}^{1/2}.\end{eqnarray} \tag{ 12 }$$

The v_r and v_θ components are expected to cancel out upon incorporation into the disk.

5.1.2. Parameters

5.1.3. Model Fitting

5.1.4. Modeling Results

We present the results from the model fitting without including the envelope in Figure 8. The disk-only models find similar protostar mass and disk radii to the models that include an envelope; the model results that include an envelope are presented in the Appendix. Figure 8 shows the results from fitting three H₂CO lines and the SO line; we also show the results from CH₃OH and NS when simultaneously fit with H₂CO. The fit parameters for each line or combination of lines are given in Table 5. The best-fitting values are determined from the median of the posterior distribution for each parameter, and the uncertainties reflect the standard deviation of the posterior distribution. In calculating the best-fitting values and their uncertainties, we filtered outliers from the posterior distributions by rejecting walkers that did not conform to a ξ² distribution. These outliers were typically walkers that never converged to the best-fitting values and typically represent less than 5% of the 200 walkers used. To avoid filtering too aggressively, we relaxed the filtering criteria such that we included at least 95% of the walkers in the final statistics.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. HOPS-370 Molecular Line Models—Exponentially Tapered Disks without Envelope

Transition(s)	Mass	Disk Radius	Disk Mass	Vlsr	Pos. Angle	Δx	Δy	γ	T(1 au)	q	log10(Mol./H2) b
	(M☉)	(au)	(M☉)	(km s−1)	(deg)	(arcsec)	(arcsec)		(K)
H2CO (J = 30,3 → 20,2)	2.72 ± 0.03	82.28 ± 0.7	0.31 ± 0.004	11.07 ± 0.02	353.86 ± 0.34	0.004 ± 0.002	0.022 ± 0.002	1.08 ± 0.01	995.92 ± 5.59	0.35	⋯
H2CO (two lines)	2.33 ± 0.03	82.79 ± 0.8	0.43 ± 0.007	11.08 ± 0.01	352.69 ± 0.39	0.011 ± 0.002	0.012 ± 0.002	1.16 ± 0.02	998.06 ± 1.72	0.35	⋯
H2CO (three lines)	2.52 ± 0.02	82.99 ± 1.5	0.54 ± 0.012	11.09 ± 0.01	352.64 ± 0.25	0.006 ± 0.001	0.014 ± 0.001	1.10 ± 0.02	998.82 ± 1.63	0.35	⋯
CH3OH (J = 42,2 → 31,2)	2.37 ± 0.03	109.29 ± 1.3	0.39 ± 0.007	10.99 ± 0.02	358.02 ± 0.58	0.010 ± 0.002	0.004 ± 0.002	0.89 ± 0.02	998.70 ± 1.89	0.35	⋯
CH3OH a , H2CO (three lines)	2.50 ± 0.02	90.36 ± 1.0	0.49 ± 0.007	11.06 ± 0.01	353.45 ± 0.24	0.007 ± 0.001	0.013 ± 0.001	1.05 ± 0.01	998.38 ± 1.09	0.35	−0.88 ± 0.013
SO (JN  = 65 → 54)	2.51 ± 0.05	104.78 ± 2.1	0.19 ± 0.005	11.42 ± 0.03	345.10 ± 0.74	0.023 ± 0.002	0.041 ± 0.002	0.80 ± 0.03	999.21 ± 1.11	0.35	⋯
SO a and CH3OH	2.37 ± 0.04	120.53 ± 0.9	0.29 ± 0.005	11.03 ± 0.03	353.00 ± 0.58	0.011 ± 0.002	0.029 ± 0.002	0.67 ± 0.02	998.44 ± 1.85	0.35	−0.88 ± 0.005
SO a , H2CO (three lines)	2.55 ± 0.02	99.29 ± 0.3	0.42 ± 0.003	11.12 ± 0.01	351.46 ± 0.24	0.009 ± 0.001	0.018 ± 0.001	0.93 ± 0.01	999.63 ± 0.47	0.35	−2.84 ± 0.203
NS	3.63 ± 0.16	79.77 ± 20.8	0.01 ± 0.006	11.05 ± 0.09	350.56 ± 1.69	−0.007 ± 0.004	0.014 ± 0.004	0.31 ± 0.24	776.35 ± 246.00	0.35	⋯
NS, SO a	1.83 ± 0.03	82.87 ± 1.3	0.02 ± 0.001	10.98 ± 0.03	352.76 ± 0.54	0.016 ± 0.002	0.062 ± 0.001	1.06 ± 0.02	998.44 ± 1.23	0.35	−0.88 ± 0.011
NS a , CH3OH	2.44 ± 0.03	109.43 ± 1.7	0.38 ± 0.008	11.00 ± 0.02	356.56 ± 0.56	0.008 ± 0.002	0.006 ± 0.002	0.89 ± 0.02	996.65 ± 3.85	0.35	−0.88 ± 0.018
NS a , H2CO (three lines)	2.56 ± 0.02	88.86 ± 1.1	0.50 ± 0.008	11.09 ± 0.01	352.46 ± 0.27	0.007 ± 0.001	0.014 ± 0.001	1.06 ± 0.01	998.87 ± 1.23	0.35	−0.88 ± 0.012
H2CO (three lines, q fit)	2.37 ± 0.02	69.69 ± 0.6	1.00 ± 0.003	11.13 ± 0.01	353.76 ± 0.31	0.001 ± 0.001	0.009 ± 0.001	1.39 ± 0.01	997.92 ± 2.28	0.001 ± 0.001	⋯

Notes. The results from each model list the molecular line(s) fit with the pdspy models. For H₂CO, "three lines" refers to all three (J = 3_n,n  → 2_n,n ) transitions, while "two lines" refers to the higher-excitation (J = 3_2,2 → 2_2,1) and (J = 3_2,1 → 2_2,0) transitions (see Table 1). The parameter q, the power-law index of the disk temperature profile, is fixed for all models except when "q fit" is listed in the description. The "Mass" column refers to the protostar mass, Δx and Δy refer to the offset of the model center of mass with respect to the image phase center, γ is the power-law index of the disk surface density profile, and SO Abund. refers to the SO abundance adopted or fit by the models.

^a Abundance of denoted molecule was allowed to vary as part of the fitting process to enable better convergence. ^b Abundance provided is for the molecules in the first column that have a superscript a.

Download table as:  ASCIITypeset image

As can be seen in Table 5, there is variation in the best-fitting protostar mass and disk radius, depending on the molecule(s) being fit. We note that the uncertainties listed in Table 5 are smaller than the differences between best-fit values for the different molecules fit; 1σ uncertainties are shown, but differences are even in excess of 3σ uncertainties. The variation of best-fit parameters can reflect both the adopted model not being fully representative of the prototstellar disk and envelope and the molecular line emission from the various molecules tracing different spatial extents. Differences in the spatial distribution of molecular line emission have been found toward several protostars on comparable spatial scales (Sakai et al. 2014; Yen et al. 2014). The models assume a constant molecular abundance with radius, and if this is not true, systematic differences between modeling of different molecules can arise.

Averaging the 12 independent fits, we find an average protostar mass of 2.5 ± 0.15 M_⊙. The full range of best-fitting protostar masses is between 1.8 and 3.6 M_⊙, with the two highest and lowest masses being quite different from the other 10 fits, which are much closer to the average value. The best-fitting disk radii are between 70 and 121 au. Averaging the 12 independent fits, we find an average radius of ∼94 ± 13 au. The uncertainties of the average are calculated using the median absolute deviation (MAD) from the collection of 12 model fits, scaled such that the MAD would correspond to one standard deviation of a normal distribution. Other average fitting parameters are V_lsr = 11.1 ± 0.04 km s⁻¹, PA = 3527 ± 14, and γ = 0.92 ± 0.14, and the position offsets are also very small, typically ∼001. We discuss the other fitted parameters and the possible ramifications of various assumptions for the modeling in the Appendix.

We compare our best-fitting protostar mass of ∼2.5 M_⊙ to the PV diagrams in Figure 6. We overlaid a Keplerian rotation curve for a mass of 2.5 M_⊙ at an inclination of 722 (see Section 5.2). The Keplerian rotation curves encompass the high-velocity emission, indicating that the average fitted mass of 2.5 M_⊙ is quite consistent with the observed data. Thus, the measured protostar mass clearly confirms that HOPS-370 is an intermediate-mass protostar.

To determine the physical structure of the disk around HOPS-370, the dust continuum emission must also be modeled using radiative transfer, but using a more realistic treatment of the disk and envelope temperature structure than the molecular line modeling employed. We model the disk and envelope around HOPS-370, also with pdspy and RADMC-3D, but taking advantage of the radiative equilibrium mode where the photons are propagated from a central luminosity source and the temperature structure of the disk and envelope are calculated self-consistently. This mode is much more computationally intensive than using prescribed temperature profiles, and as such we could not make use of this mode for the molecular line modeling. We follow the methodology employed by Sheehan & Eisner (2017), which utilizes the pdspy package, employing the same underlying MCMC sampling of the parameter space as the molecular line fitting, using 200 walkers.

The dust continuum and SED modeling within pdspy adopts the same envelope physical model as the molecular line modeling. The principal differences, however, are in the temperature and density structure of the disk. The volume density structure of the disk is defined as

$$\begin{eqnarray}&&\rho (r)={\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{-\alpha }\displaystyle \frac{{\rm{\Sigma }}(r)}{\sqrt{2\pi }\ h(r)}\,\exp \left(-\displaystyle \frac{1}{2}{\left[\displaystyle \frac{z}{h(r)}\right]}^{2}\right),\end{eqnarray} \tag{ 13 }$$

where Σ(r) is defined in Equation (3), α is the power-law index of the disk volume density profile, z is the height above the disk midplane in cylindrical coordinates, and h(r) is the vertical scale height of the disk as a function of radius. The vertical density structure for the disk in the continuum model is not defined by hydrostatic equilibrium, but parameterized as

$$\begin{eqnarray}&&h(r)={h}_{0}{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{\beta }.\end{eqnarray} \tag{ 14 }$$

Here β refers to the power-law index that defines how the vertical scale height of the disk varies with radius in the disk, not the power-law index of the dust opacity curve. The power-law index of the disk surface density profile is equivalent to γ = α − β, resulting from the multiplication of the volume density profile and the vertical density profile. The temperature structure for the continuum model is not prescribed as it is for the molecular line model. The temperature of the dust is set by the radiative equilibrium calculation performed by the RADMC-3D code. The disk dust properties are taken from Woitke et al. (2016), where the maximum size of the dust grains is parameterized as a_max and the power-law index of the dust grain size distribution as p (n(a) ∝ n^−p ), both of which are free parameters; a_min is fixed to be 0.05 μm. The envelope dust properties are taken from Pollack et al. (1994) with a_min = 0.005 μm, a_max = 1 μm, and p = 3.5.

We simultaneously fit the 0.87 mm continuum emission in the uv-plane and the SED from the near-infrared to the millimeter. The parameter space explored for the dust continuum is significantly larger than that of the molecular line kinematic modeling because now we fit the emission of the envelope, disk, and the overall SED of the system. Our goodness-of-fit metric is a weighted χ² where

$$\begin{eqnarray}&&{\chi }^{2}={\omega }_{0.87\mathrm{mm},\mathrm{vis}}{\chi }_{0.87\,\mathrm{mm},\mathrm{vis}}^{2}+{\omega }_{\mathrm{SED}}{\chi }_{\mathrm{vis}}^{2}.\end{eqnarray} \tag{ 15 }$$

The terms ω_{0.87 mm,vis} and ω_SED are determined empirically to be 0.2 and 1.0. The χ_vis are calculated by directly comparing the real and imaginary visibility components between the data and the model; the comparison is done with the two-dimensional visibility data and not using azimuthally averaged one-dimensional profiles; such profiles are only used for visual comparison of the models and data.

The best-fitting models compared to the data are shown in Figures 9 and 10. The circularly averaged visibility amplitude profile at 0.87 mm demonstrates that the model fits the 0.87 mm visibility data quite well at uv-distances greater than ∼50 kλ. Also plotted in Figure 9 is the estimated contribution of just the disk alone to the visibility amplitude profile; this shows that the envelope contributes significantly to the emission at uv-distances less than 300 kλ. The fit to the SED, also shown in Figure 9, is not perfect, but it is a close approximation to the shape and flux densities of the SED. Some flux density points are over- or underpredicted, but the relatively low angular resolution of Spitzer and Herschel at wavelengths longer than 10 μm makes it difficult to construct an SED that only includes emission from HOPS-370. However, the fact that it is the most luminous protostar within an arcminute means that contributions from other sources will not have an extremely negative impact on the SED. However, extended emission at wavelengths longer than 100 μm can cause the luminosity to be overestimated at those wavelengths.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We show an image-plane comparison of the data, model, and residuals in Figure 10. The residual images are generated from the residual visibility amplitudes and not an image-plane subtraction. The main disk feature in the continuum is well modeled and removed from the residual image, but the 0.87 mm residual image does show structure around the protostar that is not captured in the model. There are clear oversubtractions in the outer disk and center, while there are undersubtractions within the disk as well. This residual emission may represent complexities in the disk and inner envelope density structure that are not reflected in our physical model. We note that in the 0.87 mm residual map there is compact emission southeast of the protostar (−04, −03), but its nature is unclear and could stem from heating along the outflow cavity wall.

We list the parameters and their best-fitting values for the disk of HOPS-370 in Table 6, but we discuss the most relevant parameters here, which are the disk radius, disk mass, surface density power-law index, and luminosity. The best-fitting disk mass is 0.035 M_⊙, which is somewhat lower than the mass calculated under the assumption of optically thin emission and an average temperature of 131 K. However, the maximum dust grain size fit by the modeling is 440 μm, so the dust of the model will emit much more efficiently at 0.87 mm owing to its opacity than the Ossenkopf & Henning (1994) dust adopted for the simple mass estimate. We do note that the temperature of the disk in the model fit at a radius of 100 au is ∼122 K (Figure 11), which implies that the disk average temperature is comparable to the value estimated from the bolometric luminosity. The luminosity of the protostar is fit to be 276 L_☉, which is a bit less than L_bol of 314 L_☉. The more extended emission at wavelengths longer than 100 μm, which is not fit well by the model SED, could result from warm dust surrounding the protostar and may cause the bolometric luminosity to be overestimated.

Zoom In Zoom Out Reset image size

Table 6. HOPS-370 Continuum Modeling Results

Parameter Description	Parameter	Fit Value
System luminosity	L* (L⊙)	${279.9}_{-13.1}^{+14.1}$
Disk mass	Mdisk (M⊙)	${0.035}_{-0.003}^{+0.005}$
Disk outer radius	Rdisk (au)	${61.9}_{-1.2}^{+0.7}$
Disk inner radius	Rin,disk (au)	${0.54}_{-0.09}^{+0.27}$
Envelope mass	Menv (M⊙)	${0.12}_{-0.01}^{+0.01}$
Envelope radius	Renv (au)	${1881}_{-137}^{+78}$
Envelope inner radius	Rin,env (au)	=Rin,disk
Envelope centrifugal radius	Rc (au)	=Rdisk
Surface density power-law index	γ	$-{0.47}_{-0.01}^{+0.02}$
Flaring power-law index	β	${0.66}_{-0.03}^{+0.02}$
Scale height at 1 au	h0 (au)	${0.132}_{-0.008}^{+0.009}$
Inclination	i	${74.0}_{-0.2}^{+0.2}$
Density scaling in outflow cavity	fcav	${0.33}_{-0.02}^{+0.01}$
Outflow cavity shape power-law index	ξ	${1.12}_{-0.05}^{+0.06}$
Maximum dust grain size	amax (μm)	${432}_{-47}^{+31}$
Dust grain size distribution power-law index	p	${2.63}_{-0.06}^{+0.04}$
Position angle east of north	PA (deg)	${109.7}_{-0.2}^{+0.1}$
Distance	d (pc)	400.0

Note. The disk and envelope parameters derived from continuum modeling are reported here. More explanation of the parameters is provided in Sections 4.1 and 4.2. However, we explain a few that are less intuitive. The parameter f_cav is the reduction of envelope density within the outflow cavity, and the outflow cavity shape and opening angle are defined by the parameter ξ; translating to opening angle can be done by calculating 2tan⁻¹(ξ), corresponding to a full opening angle of ∼98°. The maximum dust size and grain size distribution are defined by a_max and p, where the dust grain size distribution follows a power law n(a)^−p . Absolute flux calibration uncertainty was not taken into account in the modeling.

Download table as:  ASCIITypeset image

The best-fitting envelope mass is 0.12 M_⊙ with a radius of ∼1900 au. However, the poor sampling of short uv-distances will limit the robustness of the envelope model fit to the ALMA visibility data, despite the additional constraints from the SED. This is because the SED shortward of 100 μm is most sensitive to the inner envelope density and not the overall mass or radius. The overall mass and radius can affect the longer-wavelength data more, but there is a degeneracy between dust temperature and mass. Thus, it is possible that the total mass and radius of the envelope are not accounted for in our model fit. The full outflow cavity opening angle is fit to be ∼98°, as computed from the value of ξ in Table 6. While this may seem large at first glance, it appears comparable to the width of the outflow cavities near the protostar viewed in low-velocity ¹²CO emission and shown in Tobin et al. (2019). However, the shape of the full outflow cavity extent may be more parabolic, meaning that the apparent opening angle at larger radii will appear smaller.

The disk radius from the continuum fit, 62.1 au, is smaller than our estimate of the radius from Gaussian fitting; however, Figure 11 does show that the fitted disk surface density is ∼1.3 g cm⁻² at 100 au because the exponentially tapered disk extends beyond the critical radius r_c . The continuum disk radius is also comparable to the range of exponentially tapered disk radii fit with the molecular line modeling. However, it is known that the gas disk tends to be larger than the dust disk from Class II disks (Ansdell et al. 2018). This is thought to be caused by radial drift of dust particles due to gas drag experienced because the gas orbits the star at slightly sub-Keplerian velocities (Weidenschilling 1977). Thus, the dust in protostellar disks may also experience radial drift (Birnstiel et al. 2010), which would cause a disagreement between the dust and gas disk radii.

While the disk outer radius is reasonable, there are some peculiarities with the disk structure. The surface density profile increases with radius as Σ ∝ r^0.47. However, this may result from the high opacity of the disk and its high inclination, leading to a suboptimal model fit. Moreover, given that the steepness of the exponential cutoff depends on γ, the smaller γ leads to the disk surface density falling off more quickly. Thus, the best-fitting γ may tell us more about the sharpness of the disk's "edge" than the surface density profile. The negative residuals from oversubtraction in Figure 9 could indicate that the disk needs a sharper cutoff than the exponential taper provides. In addition to Σ increasing with radius, the inner disk radius is fit to be 0.51 au. This radius is slightly smaller than the radius at which dust is expected to be destroyed, where T ∼ 1400 K, which occurs at ∼1 au in our model.

The disk vertical height is also not highly flared; the best fitting in disk scale height with radius is h(r) ∝ r^0.66, normalized to 0.131 au at 1 au, which indicates that the disk height increases slowly with radius. This fit likely reflects the large grain population of the disk and may not accurately reflect the disk properties of the smaller grains and/or gas disk.

Stellar mass is the most fundamental property of a stellar system, given that the entire evolution of a star is determined by its mass. Thus, the protostar mass measurement of ∼2.5 M_⊙ for HOPS-370 solidifies its status as an intermediate-mass protostar still in the early stages of formation. While there is some uncertainty in the mass when fitting the molecular line emission for kinematics of different molecules, the differences are most likely systematic because SO, CH₃OH, NS, and H₂CO do not trace the same material in the disk and inner envelope, meaning that their abundance profiles in the radial and vertical directions are not equivalent.

We can see this particularly in the SO and NS line emission, which is more compact than the H₂CO (J = 3_0,3 → 2_0,2) line emission (Figures 5 and 6), with H₂CO (J = 3_0,3 → 2_0,2) emission extending to slightly higher velocities on the blueshifted side. The masses from fitting are between 1.8 and 3.6 M_⊙; however, considering the full range of masses fit, the mean is ∼2.5 M_⊙ with a fractional uncertainty of ∼7%.

To determine the mass accretion rate, we use the equation for accretion luminosity

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}=\displaystyle \frac{{{GM}}_{\mathrm{ps}}\dot{M}}{{R}_{\mathrm{ps}}},\end{eqnarray} \tag{ 16 }$$

where G is the gravitational constant, M_ps is the protostar mass, $\dot{M}$ is the mass accretion rate from the disk to the protostar, and R_ps is the protostar radius. Thus, M_ps is now the most highly constrained of the parameters needed to determine $\dot{M}$. Palla & Stahler (1993) calculated pre-main-sequence evolution of intermediate-mass stars, indicating that the luminosity of the protostellar object itself (L_*) will be less than ∼10 L_⊙, meaning that the vast majority of the L_bol = 314 L_⊙ will be dominated by luminosity from mass accretion. We assume simplistically that

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}\approx {L}_{\mathrm{tot}}-{L}_{* }.\end{eqnarray} \tag{ 17 }$$

Palla & Stahler (1993) also calculated that the likely stellar radius for a protostar like HOPS-370 is ∼5 R_⊙. We note that these protostellar stellar birthline models require assumptions that are still subject to debate, the effect of accretion in particular, so there is uncertainty in the most appropriate stellar radius to assume for a given mass (Baraffe et al. 2009; Hosokawa et al. 2010).

If we assume L_bol ≈ L_tot, then the mass accretion rate from the disk to the protostar is likely 2.25 × 10⁻⁵ M_⊙ yr⁻¹. This value compares well to the outflow rate of ∼2.3 × 10⁻⁶ M_⊙ yr⁻¹ calculated from the [O i] jet emission by González-García et al. (2016) if one assumes that the outflow rate is ∼10% of the accretion rate (Shu et al. 1994; Frank et al. 2014). Due to beaming of the luminosity and foreground extinction, the actual L_tot may differ from L_bol (Whitney et al. 2003b; Offner et al. 2012). Our model fit to the continuum visibilities and SED provides a luminosity of 276 L_⊙, which is slightly lower than L_bol, implying an accretion rate of ∼1.7 × 10⁻⁵ M_⊙ yr⁻¹. The SED fit by Furlan et al. (2016) provides a much higher estimate of the luminosity at 511 L_⊙, implying an accretion rate of ∼3.2 × 10⁻⁵ M_⊙ yr⁻¹. Another SED modeling effort by Adams et al. (2012) inferred a luminosity of 300 L_⊙, which would imply that the accretion luminosity is comparable to the bolometric luminosity. In summary, the mass accretion rate appears to be well constrained to be between 1.7 × 10⁻⁵ M_⊙ yr⁻¹ and 3.2 × 10⁻⁵ M_⊙ yr⁻¹; for comparison, these accretion rates are ∼1000× the typical accretion rates found in low-mass T Tauri stars (e.g., Ingleby et al. 2013; Alcalá et al. 2017).

The inferred infall rate from the envelope to the disk is 3.2 × 10⁻⁵ M_⊙ yr⁻¹ from our best-fitting model. This is very comparable to previous estimates derived from SED fitting of 2.2 × 10⁻⁵ M_⊙ yr⁻¹ (Furlan et al. 2016) and 4.4 ×10⁻⁵ M_⊙ yr⁻¹ (Adams et al. 2012), both of which are scaled to reflect the protostar mass of 2.5 M_⊙. The inferred range of accretion rates from the disk to the protostar are comparable to the envelope to disk infall rates derived from SED fitting and our best-fitting model. This indicates that the envelope is supplying the disk with mass at a comparable or higher rate as compared to how rapidly the disk material drains onto the protostar.

All of the aforementioned analytic models, however, assume simplified geometries for the envelope, disk, and outflow cavity. The differences in the parameters from different models suggest that the underlying physical models may not accurately describe the true structure of the protostellar system. However, the analytic models ignore the potential effects of turbulence and magnetic pressure support, which may also play a role in regulating the infall from the envelope to disk (e.g., Seifried et al. 2013; Li et al. 2014). For this reason, the mismatches in the infall and accretion rates from different model fits could arise from physical model inadequacies and degeneracies due to only fitting the SED in some cases. Finally, while we are quoting the infall rates derived from these models, the precise numbers should be regarded with caution for the reasons outlined in this paragraph.

Based on our knowledge of the disk mass and protostar mass from modeling, we can estimate how important self-gravity is to the HOPS-370 disk using Toomre's Q and its relationship to the disk and protostar mass

$$\begin{eqnarray}&&Q\approx 2\displaystyle \frac{{M}_{\mathrm{ps}}}{{M}_{d}}\displaystyle \frac{H}{r}.\end{eqnarray} \tag{ 18 }$$

Note that this is not the conventional representation for Q, but is rewritten for a rotationally supported disk around a central gravitating body, in our case a protostar (Kratter & Lodato 2016; Tobin et al. 2016a). The disk scale height is H (equivalent to c_s /Ω, where c_s is the disk sound speed and Ω is the Keplerian angular velocity), M_ps is the protostar mass, and M_d is the disk mass. We calculate Q at a radius of 50 au and find that Q ∼ 6, assuming the inferred disk gas mass from modeling (0.035 M_⊙), protostar mass (∼2.5 M_⊙), and a typical disk temperature of 131 K (c_s  = 0.56 km s⁻¹) at a radius of 50 au (Ω ∼ 8.9 × 10⁻¹⁰ for M_* = 2.5 M_⊙). Thus, the disk around HOPS-370 is not expected to be highly self-gravitating. We note that Q could be lower if the disk mass was higher, but the disk mass would have to be significantly underestimated for Q to approach 1.

We also calculated Q using the temperature and density structure derived from the best-fitting radiative transfer model to the continuum data. We extracted the disk surface density and midplane temperature profile from the best-fitting radiative transfer model (plotted in Figure 11). Then, using the surface density and temperature, we calculated Q as a function of radius using

$$\begin{eqnarray}&&Q=\frac{{\rm{\Omega }}(r){c}_{s}(r)}{\pi G{\rm{\Sigma }}(r)},\end{eqnarray} \tag{ 19 }$$

where Σ is the radial surface density profile and G is the gravitation constant. Figure 11 shows Q as a function of radius. The radial distribution of Q shows that the disk does not approach instability at any radius and the minimum value of Q is ∼14, even larger than the approximate calculation from Equation (18). This further demonstrates that self-gravity is not likely important in the disk of HOPS-370 at the present time.

The near-infrared spectral features were presented in Section 4.4, where the principal features of importance with respect to the accretion rate and luminosity of the protostar are the Brγ line emission and CO band head emission. Connelley & Greene (2010) found a strong correlation between Brγ emission, veiling, and CO band head emission. The authors inferred that when Brγ and CO band heads are in emission and veiling is high, mass accretion is also high. Modeling of the CO band head spectra for different stellar types and accretion rates by Calvet et al. (1991) predicted when CO band head emission should appear in absorption or emission for a given mass accretion rate and stellar effective temperature. For example, very high accretion rate systems, like FU Ori systems, have CO band head absorption due to the hot accreting disk midplane and cooler disk atmosphere.

For the accretion rate needed to produce the observed bolometric luminosity with the expected stellar temperature of ∼4500 K (Palla & Stahler 1993) for HOPS-370, the models of Calvet et al. (1991) indicate that CO band head absorption should be expected rather than the observed emission. On the other hand, Najita et al. (1996) suggested that high accretion rates could still lead to CO band head emission because accretion through the disk could cause a temperature inversion in the inner disk, making the surface hotter than the midplane, leading to emission. In addition, disk accretion is not the only possible mechanism for producing CO band head emission; other studies suggest that winds could also produce the emission (e.g., Chandler et al. 1995), but higher spectral resolution and S/N is required to differentiate between a wind and disk origin.

While the exact origin of the CO band head emission in HOPS-370 is uncertain, the CO band head and Brγ emission demonstrate that despite the high inferred accretion rate from the luminosity, the spectrum of HOPS-370 is decidedly not FU Ori-like since FU Ori-type spectra have no detectable Brγ emission and the CO band heads are in absorption (Connelley & Greene 2010; Fischer et al. 2012). Brγ line luminosity has been used to infer the accretion luminosity of young stars, including protostars with relationships defined by Muzerolle et al. (1998, L_acc = 1.26 log₁₀(L_Brγ ) + 4.43). Applying this relationship, our inferred L_acc for HOPS-370 from Brγ emission is 2.2 × 10⁻² L_⊙, which is at odds with the accretion luminosity inferred in Section 6.1 by a factor of 10,000. The Brγ line luminosity has not been corrected for attenuation by dust extinction, and correction will only lead to higher accretion luminosities. Furthermore, the Brγ emission we detect is from scattered light, and radiative transfer models from Whitney et al. (2003a) indicate that the amount of emergent flux at 2 μm in the outflow cavities can be between 2 and 4 orders of magnitude lower than the input stellar spectrum. Therefore, it is plausible that the Brγ emission we observe is originating from accretion, but it is difficult to accurately infer an accretion rate from the line emission owing to the combined effects of extinction and observing the spectrum in scattered light.

The most current compilation of protostar masses is found in Yen et al. (2017), which contains protostar masses measured by ALMA, the Plateau de Bure Interferometer, and the Submillimeter Array. HOPS-370 is one of the most massive protostars to have a kinematic mass measurement. Comparable protostars are HH111 MMS (1.8 M_⊙), Elias 29 (2.5 M_⊙), R CrA IRS7B (2.3 M_⊙), Oph IRS 43 (1.9 M_⊙), and L1489 IRS (1.6 M_⊙). However, these are all Class I protostars, except for R CrA IRS7B, which is a borderline Class 0/I protostar, and HOPS-370 is the only one with L_bol > 100 L_⊙. Of these protostars with comparable masses, the highest-luminosity system is HH111 MMS at 17.4 L_⊙. Thus, HOPS-370 is unique and requires a very high accretion luminosity to explain its L_bol. The accretion rate for HOPS-370 inferred from the luminosity is greater than an order of magnitude larger than the other protostars with a similar mass, and it has a higher inferred accretion rate than all other protostars listed in Yen et al. (2017).

The high rate of accretion in HOPS-370 begs the question, is this an outbursting source or a higher-mass star being formed through sustained infall from the envelope to disk? Outbursts and variability seem to be common among protostars (e.g., Hartmann et al. 1996; Dunham et al. 2010; Audard et al. 2014; Safron et al. 2015; Fischer et al. 2017), which are thought to reflect higher accretion rates than average for short intervals of time (hundreds to thousands of years). However, the study of intermediate-mass star formation has been complicated by multiplicity for systems other than HOPS-370. We will discuss the merits of the two possible accretion scenarios for HOPS-370.

6.5.1. An Outbursting Protostar in a High Accretion State?

6.5.2. Sustained Accretion?

The alternative is that we are not witnessing an outburst, but sustained high accretion rates that could be typical for the formation of an intermediate-mass protostar. Models of intermediate- to high-mass star formation (e.g., McKee & Tan 2003; Wuchterl & Tscharnuter 2003) predict the luminosities over time during the formation of such systems. While these analytic models do not assume accretion through a disk, but rather direct infall from the envelope onto the protostar, they demonstrate a scenario in which a protostar system that will ultimately form a high-mass star will have a significantly higher overall luminosity during its formation, as compared to the stellar luminosity of the central protostar due to a sustained high accretion accretion rate. Moreover, competitive accretion models for protostars forming within clustered environments, not unlike the OMC2 region, also predict that the protostars that ultimately have higher final masses will accrete at higher rates (Bonnell et al. 2001; Hsu et al. 2010; Bate 2012).

The estimates of the accretion rate from the disk to the protostar, based on the bolometric luminosity and protostar mass, range from 1.7 × 10⁻⁵ M_⊙ yr⁻¹ to 3.2 × 10⁻⁵ M_⊙ yr⁻¹. Then, the estimates of the infall rate from the envelope to the disk range between 2.2 × 10⁻⁵ M_⊙ yr⁻¹ and 4.4 × 10⁻⁵ M_⊙ yr⁻¹. Thus, the disk is being fed with material rapidly enough that accretion could be sustained at its high rate, regardless of the spread in the estimates. However, the infall rates are model dependent, and their accuracy is more questionable than the disk-to-star accretion rates.

These high accretion rates for HOPS-370 indicate that it could gain another solar mass of material within ∼3 × 10⁴ yr to ∼9 × 10⁴ yr. However, the envelope surrounding HOPS-370 may not have sufficient mass to sustain high accretion indefinitely. The envelope in our model fit only had ∼0.12 M_⊙ of material out to ∼1900 au. However, our modeling does not account for the total mass available in the envelope; other observational measurements find ∼2.5 M_⊙ from 3 mm continuum data (Kainulainen et al. 2017). Even if the immediate surroundings of HOPS-370 have limited mass, it is embedded within the dense molecular environment of the northern integral-shaped filament. This means that the total reservoir that HOPS-370 could accrete from is larger than the envelope mass fit from modeling.

Even though we expect the disk to be gravitationally stable, material must still accrete through the disk with a sustained high accretion rate that keeps the luminosity around the present value of ∼300 L_⊙. If the disk self-gravity is negligible, then gravitational instability cannot efficiently transport angular momentum anywhere in the disk (e.g., Rice et al. 2010). Even if disk self-gravity is not large enough to drive accretion, there are alternative ways to promote the accretion of material through the disk toward the protostar. One such way is the excitation of spiral density waves in the disk due to infalling material (Lesur et al. 2015). In this scenario, the infalling material has lower specific angular momentum than the disk at the point where material arrives at the disk. This creates an unstable accretion shock that promotes the formation of spiral arms and outward angular momentum transport, enabling efficient accretion from the outer disk to the inner disk.

While the MRI or some other viscous process may be responsible for accretion between the inner and outer disk, the exact mechanism remains debated. However, the challenge for MRI is to have enough ionization such that the magnetic field is coupled strongly enough to the gas to transport angular momentum. Work by Offner et al. (2019) has shown that cosmic rays produced by accretion could increase the ionization enough to enable MRI in the inner ∼10 au of the disk. Assuming that accretion proceeds owing to turbulent viscosity, we can use the observationally inferred accretion rate to constrain the necessary inner disk properties using the equation

$$\begin{eqnarray}&&\dot{M}\sim 3\pi \alpha {c}_{s}H{\rm{\Sigma }},\end{eqnarray} \tag{ 20 }$$

following Hartmann (2008). The term α refers to the α-viscosity of the disk (Shakura & Syunyaev 1973), Σ is the surface density of the disk at the radius in question, and the other parameters are as defined previously. Rearranging these terms to solve for Σ at 1 au with ${\dot{M}}_{\mathrm{acc}}=1\times {10}^{-5}$ M_⊙ yr⁻¹, we adopt α = 0.1, which is typical for a high accretion rate (Zhu et al. 2009), c_s  = 1.86 × 10⁵ cm s⁻¹ (from T = 2000 K), and H = 8.98 × 10¹¹ cm (from the continuum model, Table 6). With these values, we find that Σ at 1 au should be ∼3700 g cm⁻² for the disk to accrete at 1 × 10⁻⁵ M_⊙ yr⁻¹. Our best-fitting pdspy continuum model only has a surface density at 1 au of ∼4.6 cm² g⁻¹, which is inconsistent with the need for a high surface density in the inner disk to facilitate accretion. Even if we extrapolate the maximum disk surface density of 19 g cm⁻² at 40 au to 1 au assuming γ = 1, we would find a surface density of ∼760 g cm⁻², several times lower than the value needed to sustain accretion at ∼1 × 10⁻⁵ M_⊙ yr⁻¹. Thus, this can be taken as further evidence that the surface density profile of our best-fit disk model is inconsistent with other characteristics of the system and is likely not well constrained from the current model. Furthermore, the disk surface density depends on the dust grain size distribution, and if the dust opacities are too high (the best-fitting a_max ∼ 440.4 μm), then the surface density will be too low. With only a single wavelength, a_max can have a degeneracy with disk mass. Despite the uncertainties in the underlying density structure of the disk model and dust opacity, the surface density would need to be an order of magnitude larger than currently observed for self-gravity to be important in the disk of HOPS-370. It is unlikely that the true surface density could be so large and fit the observed data.

Disk winds have been promoted as an alternative mechanism to promote accretion in the inner disks of young stars (e.g., Pudritz et al. 2007; Bai & Stone 2013). However, it is unclear whether disk winds can promote accretion at rates as high as ∼10⁻⁵ M_⊙ yr⁻¹. This is because disk winds can only be active in a thin, ionized layer in the disk, and the overall surface density of this layer is significantly lower than the overall disk surface density.

The peculiarities of the inner disk fit aside, it seems plausible that HOPS-370 could be forming an intermediate-mass star from steady accretion supplied from the envelope to the disk and the disk to the star, without an outburst being necessary.

The current protostar mass of HOPS-370 is comparable to that of disk-hosting Herbig Ae/Be stars, which are intermediate-mass pre-main-sequence stars, with masses between ∼1 and ∼4 M_⊙. Herbig Ae/Be systems typically have disk radii that are comparable to or larger than that of the HOPS-370 disk. Thus, HOPS-370 could be a progenitor of these Herbig Ae/Be systems if the protostar does not grow significantly larger in mass and a reasonably massive disk remains for roughly a few megayears once infall from the envelope has stopped.

Kama et al. (2020) examined the disk masses of 15 Herbig Ae/Be systems, which have stellar masses ranging between ∼1.5 and 3.0 M_⊙ using the gas mass upper limits provided by HD (J = 1 → 0) observations from the Herschel Space Observatory. The upper limits from HD were then interpreted in the context of dust disk masses measured from continuum emission and gas masses measured using CO molecular line mission. Many of these disks show ring or asymmetric structures in their dust continuum emission that have been interpreted as signatures of giant planet formation (e.g., Isella & Turner 2018; Keppler et al. 2018; Pinte et al. 2018; Teague et al. 2018). However, the upper limits of the gas masses of these disks indicate that they are most likely gravitationally stable at present. The projected planet masses require conversion of disk mass to planet masses at efficiencies >10%, and these efficiencies are high enough that the presence of these putative giant protoplanets is easier to explain if they formed earlier during the embedded protostar phase, possibly enabled by gravitational instability (Tychoniec et al. 2020). But HOPS-370 in its current state with a high bolometric luminosity is gravitationally stable and is expected to still be stable even if its luminosity was 10× lower.

Between the two possible mechanisms for giant planet formation, core accretion and gravitational instability, core accretion seems to be the most favorable in HOPS-370. Gravitational instability in the disk of HOPS-370 is unlikely because the disk would have to remain massive, while the system luminosity decreased significantly, presumably due to lower accretion. Giant planet formation via core accretion is also plausible; the cores could form during the protostar phase, and they would be able to begin gas accretion while the disk is still being fed with infalling material. Constant replenishment of the disk with gas and dust would enable protoplanets to grow to Jupiter masses. The efficiency of disk mass to planet mass is not required to be high because there will be a large mass flux onto and through the disk.

Furthermore, the current disk mass of ∼0.035 M_⊙ and the current accretion rate from the disk to the star of ∼10⁻⁵ M_⊙ yr⁻¹ indicate that the entire gas disk could be accreted in less than 10⁴ yr. However, the properties of the current disk are still relevant for the further evolution of solid material because the gas disk may be accreted, but the solid material within the disk may not be completely accreted with the gas. This is because pressure bumps in the disk can trap large dust particles, promoting further dust growth (e.g., Pinilla et al. 2012). It has been demonstrated that infalling material to the disk can trigger Rossby wave instabilities and provide a mechanism for particle trapping (Bae et al. 2015). Thus, while the disk will continue to evolve through the protostellar phase, its current properties are not completely removed from the initial conditions of planet formation.

On the other hand, HOPS-370 may continue to gain mass and become a higher-mass star, which could render it significantly different from the Herbig Ae/Be systems discussed here. Moreover, HOPS-370 is currently within a relatively dense cluster of protostars, while the Herbig Ae/Be systems studied by Kama et al. (2020) are much more isolated. Thus, it is unclear whether a long-lived disk around HOPS-370 would be able to survive a few megayears in a clustered environment.

We also fitted the molecular line emission including an envelope component whose mass and radius were free parameters in the fitting. The best-fitting parameters from these model fits are presented in Table 7, and we show representative model fits with an envelope in Figure 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The range of protostar masses fit is comparable to the disk-only fits. The average mass from the 12 independent fits is 2.6 ± 0.2 M_⊙, only larger than the disk-only fit by 0.1 M_⊙. The best-fitting protostar mass values range between 1.85 and 3.83 M_⊙. Like the disk-only fits, the largest outliers are driven by the NS-only fit (3.83 M_⊙) and the NS and SO fit (1.85 M_⊙).

The best-fitting disk radii are between 14 and 114 au. The average disk radius is ∼68 ± 21 au, smaller than the average from the disk-only fits by 26 au. The smallest disk radius comes from the NS and SO fit, which also has the smallest protostar mass. Other average fitting parameters are V_lsr = 11.2 ± 0.07 km s⁻¹, PA = 3526 ± 13, and γ = 1.08 ± 0.12, and the position offsets are also very small, typically ∼001.

The protostar masses derived from fitting with an envelope component appear robust and agree well with the disk-only component. Indeed, the averages are consistent within their associated uncertainties. However, the range of disk radii fit points to degeneracy when attempting to fit both an envelope component and a disk simultaneously. This is likely due to the fact that the spatial scales of both structures are not well sampled in the observations and the molecular abundance is a constant for both the envelope and disk. For these reasons, we regard the range of disk radii derived from the disk-only fits with greater confidence. It is important to point out that some of the combinations of molecular lines fitted with both a disk and envelope component do fit disk radii that are consistent with the disk-only fits. The line fits with the three H₂CO transitions and another line in particular performed well.

Several other parameters that were fit during modeling most likely do not reflect the actual physical parameters that they correspond to owing to model limitations and limitations of the data themselves. This applies to both the disk-only fits and the disk-plus-envelope fits. The disk masses are unreliable in an absolute sense because they depend on the assumed abundance of the molecule, which we assume is constant throughout the disk. The envelope masses and radii are considered unreliable because of the limited sensitivity to spatial scales beyond ∼4'' and the fixed molecular abundances. The temperature normalization of the disk (T(1 au)) is also likely unreliable because the fitted molecules have little temperature sensitivity, leading to a degeneracy between mass and temperature. Finally, the molecular abundance fits for the cases where we fit more than one species are also unreliable since these were only varied to improve fitting with a constant disk mass.

While we fixed the power-law index of the temperature profile q for all the models, we did examine the impact of this assumption. We ran another fit to all H₂CO transitions because their different excitation energies can enable their line ratios to provide constraints on temperature. Allowing q to vary resulted in a best-fitting protostar mass that was 0.29 M_⊙ lower than when q was fixed at 0.35. However, this fit is likely unrealistic because q was fit to be consistent with 0, indicating no temperature change with radius, and this was also the lower limit for the parameter. Thus, the overall impact of fixing q results in ∼10% higher masses, but the best-fitting value for q was not realistic, indicating that our choice of a fixed q is more reasonable; this particular model is not included in the average protostar mass and radius.

With regard to the protostar mass and disk radius, the model fits with the most discrepant fitting parameters were NS alone, NS, and SO. The fit to NS alone resulted in the highest protostar mass of all the fits for both the disk-only and disk-plus-envelope models. Then, the lowest protostar masses fit were from the NS and SO fits. Furthermore, the NS and SO fits with an envelope component resulted in a model fit with a very small disk radius. The reason for these discrepant fits could be that both NS and SO had the smallest spatial extent of line emission in our data. This potentially explains the small disk radii. The more concentrated emission from NS likely resulted in a discrepant protostar mass determination because there was less of a rotation curve to model. Thus, we regard the fits to multiple molecular lines as the most robust, particularly when the majority of lines sample the full range of radii in the disk.

To perform the fit for CH₃OH, we had to first subtract its extended, off-center emission from the visibilities. This is because the brightest CH₃OH emission is not centered on HOPS-370, but comes from a spot on the edge of the outflow cavity wall several arcseconds south of the disk. This is outside the field of view in Figure 5 and not shown because it is not relevant for this analysis. This emission was preventing model convergence because it dominated the imaginary component of the visibility amplitude data. To remove this emission, we interactively cleaned the CH₃OH emission, only cleaning the emission that was not associated with the disk. This created a model for the CH₃OH emission not associated with the disk. We used clean to save the model to the measurement set, and we subtracted it from the data using the CASA task uvsub. This made it possible to reliably model the CH₃OH emission only coming from the vicinity of the disk.