APEX Observations of the CO Envelope around the Young FUor-type Star V883 Ori

TL;DR

This study uses APEX observations to analyze the CO gas envelope around the young FUor star V883 Ori, revealing an extended envelope and infall rates consistent with episodic accretion events.

Contribution

First detailed CO envelope mapping around V883 Ori, linking envelope structure to episodic accretion in FUor-type stars.

Findings

Detected extended CO emission out to >10,000 au

Estimated envelope mass of 0.2-0.4 solar masses

Inferred mass-infall rate of 1-2 x 10^{-6} solar masses per year

Abstract

The accretion-driven outbursts of young FU Orionis-type stars may be a common stage of pre-main sequence evolution and can have a significant impact on the circumstellar environment as it pertains to the growth of solids and eventually planets. This episodic accretion is thought to be sustained by additional gas in-falling from the circumstellar envelope and disk. We present APEX observations of the CO gas in the envelope around V883 Orionis, a young outbursting star. The observations mapped the $^{12}$CO(4-3), $^{12}$CO(3-2), and $^{13}$CO(3-2) lines with the FLASH$^{+}$ instrument and the $^{12}$CO(6-5) line with the SEPIA instrument. We detected high signal-to-noise emission extending out to radii $>10000$ au and calculated integrated fluxes of $1100~\rm Jy~km~s^{-1}$ for $^{12}$CO(6-5), $2400~\rm Jy~km~s^{-1}$ for $^{12}$CO(4-3), $1600~\rm Jy~km~s^{-1}$ for $^{12}$CO(3-2), and $450~\rm Jy~km~s^{-1}$ for $^{13}$CO(3-2). We used the thermo-chemical code P{\small RO}D{\small I}M{\small O} to test several models and find the data are best described by an envelope structure with $\rm M_{env}\approx 0.2-0.4\,M_{\odot}$ and a mass-infall rate of $\rm \dot{M}_{inf}=1-2\times10^{-6}\,M_{\odot}\,yr^{-1}$. We infer that the observed envelope and outflow structure around V883 Ori could be caused by multiple outbursts, consistent with episodic accretion.

Figures

Tables

Table 1. Table 1: Summary of our APEX observations. The line σ r ​ m ​ s subscript 𝜎 𝑟 𝑚 𝑠 \sigma_{rms} is from the line-free channels, re-binned to 1.0 ​ km ​ s − 1 1.0 km superscript s 1 1.0~{}\rm km~{}s^{-1} . The peak, FWHM, and envelope flux (Env. Flux) of a given line is taken from fitting a Gaussian to the lower-velocity peak in the integrated spectra in Fig.  1 . The uncertainties on the peak and FWHM are the 1 ​ σ 1 𝜎 1\sigma uncertainties from the co-variance matrix in the Gaussian fitting. The uncertainties on the envelope flux are σ r ​ m ​ s subscript 𝜎 𝑟 𝑚 𝑠 \sigma_{rms} , model fitting uncertainty, and absolute flux calibration uncertainty added in quadrature. The absolute flux calibration uncertainty is taken as 15 % percent 15 15\% for 12 CO(6-5) and 10 % percent 10 10\% for all other lines and is the dominate source of flux uncertainty. The beam efficiencies and Jy ​ K − 1 Jy superscript K 1 \rm Jy~{}K^{-1} conversion factors for each observational setup are given in Sec.  2 .

Species	Observation Date	Frequency	Beam	${\sigma }_{rms}$	Ch. Width	Peak	FWHM	Env. Flux
	YYYY MMM DD	[GHz]	[′′]	[Jy beam-1]	[$\mathrm{km}{s}^{-1}$]	[Jy beam-1]	[$\mathrm{km}{s}^{-1}$]	[Jy km s-1]
12CO(6–5)	2018 OCT 15/21	691.473	9.6	2.2	0.20	$1100\pm 45$	$0.38\pm 0.02$	$1100\pm 170$
12CO(4–3)	2016 SEP 06	461.041	14.4	1.4	0.050	$2800\pm 20$	$0.35\pm 0.01$	$2400\pm 240$
12CO(3–2)	2016 SEP 06	345.796	19.2	0.53	0.033	$1800\pm 10$	$0.36\pm 0.01$	$1600\pm 160$
13CO(3–2)	2016 SEP 06	330.588	20.1	0.81	0.033	$680\pm 6$	$0.26\pm 0.01$	$450\pm 45$

Table 2. Table 2: Comparison of our models to the continuum model and selected photometric fluxes of Furlan et al. ( 2016 ) . The parameters are described in the Sec.  3 and B and in Furlan et al. ( 2016 ) .

Model	$L$	${\dot{M}}_{\inf }$	350 $\mu m$ Flux	$870\mu m$ Flux
	$\left[L_{\odot }\right]$	$\left[M_{\odot }{\mathrm{yr}}^{-1}\right]$	[Jy]	[Jy]
Furlan et al. (2016)	194.67	$\approx 4\times {10}^{-6}$	$11.3\pm 4.5$	$1.2\pm 0.2$
Model 1	200	$2\times {10}^{-6}$	$3.4$	$0.2$
Model 1 + disk	200	$2\times {10}^{-6}$	$14.3$	$1.6$
Model 2	200	$1\times {10}^{-6}$	$1.9$	$0.1$
Model 2 + disk	200	$1\times {10}^{-6}$	12.5	1.46

Full text

circumstellar matter, stars: individual (V883 Orionis), stars: pre-main sequence, submillimeter: stars

††journal: ApJ††facilities: APEX(SEPIA, FLASH+)††software: GILDAS, PRODIMO (Woitke et al., 2009; Woitke et al., 2016; Rab et al., 2017)

1 Introduction

During a star’s pre-main sequence evolution, material will be accreted from the disk onto the protostar, fueling its growth. A common stage of this evolution may be periods of significantly enhanced mass infall, known as episodic accretion. This is evidenced by outbursts that increase a system’s apparent brightness by several magnitudes and total luminosity by as much as $10-100\times$. Understanding the underlying mechanisms that cause the observed outbursts of young stellar objects (YSOs) is important not only for building a complete picture of star formation, but also for the potential implications on the planet formation process. Young outbursting systems can be categorized as either FUors or EXors. FUors, named after FU Orionis (Herbig, 1966), are embedded systems that have very strong outbursts that can last $>100$ yr. These stars will have optical spectra similar to rapidly rotating earlier-type giant stars and infrared (IR) spectra more similar to M-giant stars (Hartmann & Kenyon, 1996). EXors, named after EX Lupi (Herbig, 1989), have outbursts that are weaker and shorter lived as compared to FUors. Indeed, FUors and EXors could represent different evolutionary stages of the same types of objects as EXors typically have lower disk masses and are no longer embedded. The likely cause of both these types of outbursts is episodic accretion of material from the circumstellar disks onto the protostars. The transport mechanism(s) that delivers additional material onto protostars, however, is still unknown (see Audard et al., 2014, for a review of outbursting YSOs).

Winds and molecular outflows are common components of star formation and disk evolution (Williams & Cieza, 2011). These components can carry away angular momentum, allowing material in the circumstellar disk to accrete onto the star. As material is transported inward, the disk is thought to be replenished by material from the large envelopes, which can extend to 1000s of au (e.g., Cassen & Moosman, 1981; Terebey et al., 1984). This “refueling” of material can support the episodic nature of the accretion events. Therefore, studying the envelopes of FUors is crucial to understanding the nature of accretion events in young stars (e.g., Quillen et al., 2005). The envelopes of FUors can be inferred through absorption seen in the $10~{}\mu$m silicate feature (Quanz et al., 2007) or through IR excess (e.g., Green et al., 2006). Direct observations of the envelope gas, however, are less common but crucial to constraining the mass and morphological structure (Kóspál et al., 2017). An accurate characterization of the envelope around an FUor can determine the timescale and strength of the mass infall rate, will yield insight into how many future outburst events can be sustained, and can help distinguish the difference between the circumstellar environments of FUors and that of other YSOs.

V883 Ori is an embedded $1.3\pm 0.1~{}\rm M_{\odot}$ (Cieza et al., 2016) pre-main sequence star located in the Lynds 1641 cloud in the Orion star forming region at a distance of $414\pm 7$ pc (Menten et al., 2007). Strom & Strom (1993) first suggested that V883 Ori may be an FUor as it shares similar optical spectral features to FU Orionis and has a flat spectral energy distribution (SED). High resolution near-IR spectra further indicate that V883 Ori is best classified as an FUor (Greene et al., 2008). While V883 Ori has not been observed in a pre-outburst state, analysis of photographic plates indicate the outburst likely began before 1888 (Pickering, 1890). Coupled with a potential ALMA-resolved water snowline (Cieza et al., 2016) and the presence of complex organic molecules (van’t Hoff et al., 2018; Lee et al., 2019), V883 Ori is a key testbed system for understanding outbursts and the early stages of star formation and disk evolution.

In this paper, we present APEX FLASH*+* and SEPIA B9 observations of the CO in V883 Ori’s envelope. In Sec. 2, we describe the data reduction and observational results. In Sec. 3, we detail the envelope model. In Sec. 4, we discuss the implications of the model. In Sec. 5, we summarize the results.

2 Observations

We observed V883 Ori with the Atacama Pathfinder EXperiment (APEX), a 12 m diameter millimeter-wave telescope located on the Llano de Chajnantor in the Chilean Atacama desert (Güsten et al., 2006). We used the FLASH+ receiver to measure the 12CO(3–2), 13CO(3–2), and 12CO(4–3) lines on 2016 September 6 (Project 098.F-9505(A), PI: Th. Henning), and the SEPIA B9 receiver to measure the 12CO(6–5) line on 2018 October 15 and October 21 (Project 0102.F-9802, PI: J.A. White). The SEPIA observations were part of the Science Verification of the B9 receiver. For each line, we used an On-The-Fly (OTF) mapping procedure centered at the V883 Ori J2000 coordinates of R.A. = 05h38m18s.097 and DEC. = $-07^{\circ}02^{\prime}25\farcs 96$. The main parameters of our observations are summarized in Table 1.

FLASH+ is a dual-frequency heterodyne receiver, operating simultaneously in the 345 GHz and 460 GHz atmospheric windows, providing 4 GHz bandwidth in each sideband (Klein et al., 2014). The lower frequency channel was tuned to 344.2 GHz in the upper sideband (USB) to simultaneously cover the 13CO(3–2) and 12CO(3–2) lines at 330.588 GHz and 345.796 GHz, respectively. The higher frequency channel was tuned to the 12CO(4–3) line at 461.041 GHz in USB. We used the XFFTS backends, providing a 38 kHz ($0.033~{}\rm km~{}s^{-1}$) spectral resolution for both the 3–2 lines and 76 kHz ($0.050~{}\rm km~{}s^{-1}$) for the 4-–3 line. We obtained 120${}^{\prime\prime}{\times}$120*′′* OTF maps using a relative reference off position 2500*′′* away in R.A. Before the observation, we checked in total power mode that the off position was clean of CO emission. The precipitable water vapor (PWV) was on average 0.5 mm during our observations.

SEPIA is the Swedish-ESO PI receiver for APEX (Belitsky et al., 2018). The SEPIA Band 9 receiver is a dual polarization 2SB receiver (Hesper et al., 2018; Baryshev et al., 2015) built to the specifications of the ALMA Band 9 receiver (600–720 GHz). Currently, the XFFTS backend records 4 GHz bandwidth in each sideband. The USB was centered on the 12CO(6–5) line at 691.473 GHz with 450 kHz ($0.20~{}\rm km~{}s^{-1}$) resolution. We obtained 150${}^{\prime\prime}{\times}$150*′′* OTF maps, while the relative reference off position was 1000*′′* away in R.A. The PWV was 0.5–0.7 mm during our two observations.

The data were reduced using the GILDAS CLASS package111http://www.iram.fr/IRAMFR/GILDAS/ and a first-order baseline was removed from the spectra. The main beam efficiencies determined by observing Mars are 0.73, 0.60, and 0.49 (APEX Team, private communication) at 352, 464, and 691 GHz, respectively. The data were converted to Jy using 41 Jy K*-1*, 48 Jy K*-1*, and 79 Jy K*-1* at 352, 464, and 691 GHz, respectively. The telescope’s beam size is 20$\farcs$1, 19$\farcs$2, 14$\farcs$4, and 9$\farcs$6 for the 13CO(3-–2), 12CO(3–2), 12CO(4–3), and 12CO(6–5) lines, respectively. The $\sigma_{rms}$ for each line is calculated by re-binning each spectrum to a $1~{}\rm km~{}s^{-1}$ channel width, then calculating the root-mean-square of all the line-free channels. The resulting values are summarized in Table 1. The absolute flux uncertainty is still unconstrained for these instruments but is taken as $10\%$ at 352 GHz and 464 GHz and as $15\%$ at 691 GHz.

2.1 Observational Results

The spectra of each molecular species is presented in Fig. 1, spatially integrated over a region of 10000 au. The spectra from V883 Ori’s envelope are expected to have a single peak following an approximately Gaussian profile given the resolution of the observations. However, there are clearly 2 peaks in the emission seen in 12CO(6–5) and 13CO(3–2) and 3 peaks seen in 12CO(4–3) and 12CO(3–2). The best fit positions of the lower-velocity peak in each spectrum are all centered at $\sim 4.1~{}\rm km~{}s^{-1}$ (the approximate systemic velocity). The other emission features are likely a combination of outflow material (which flows along the SE/NW axis; Ruíz-Rodríguez et al., 2017a), a non-spherical envelope structure, or diffuse emission from the surrounding molecular cloud. This is best seen in the velocity channel maps of each line (see Appendix A). For example, Fig. 4 shows the 12CO(6–5) line channel map. In this figure, the emission from the envelope (lower-velocity peak in the spectra) is largely contained within radii $<10000$ au at velocities $\rm<5.0~{}km~{}s^{-1}$ while the emission at velocities $\rm>5.0~{}km~{}s^{-1}$ is more consistent with diffuse emission at larger radial separations. The emission of the higher-velocity components fills the field of view in several velocity channels and does not show a concentration or a peak towards the position of V883 Ori. Those structures might be velocity coherent structures within the molecular cloud complex L 1641, similar to the “fibers” in the Orion Integral Filament mapped in detail by Hacar et al. (2018).

The compact emission in the channels near to the systemic velocity (e.g., $3.4~{}\mathrm{km~{}s^{-1}}$ and $4.6~{}\mathrm{km~{}s^{-1}}$ in Fig. 1) suggests in-falling material, as higher velocities are expected closer to the center (i.e., free-fall velocity $\propto 1/\sqrt{r}$). To calculate the peak, FWHM, and total envelope flux for each species, we first integrated each spectrum over a 10000 au radius aperture centered on V883 Ori. A multi-component Gaussian was fit to each spectrum in order to properly differentiate between emission from the envelope and the other peaks in the spectra. The best fit parameters from the lower-velocity Gaussian as seen in Fig. 1 are summarized in Table 1. The integrated intensity maps ($0^{th}$ order moment) for each species are shown in Fig. 2. Each of these maps only shows the integrated intensity from the lower-velocity peak of each spectra.

3 Envelope Model

We present models for the envelope structure applying the radiation thermo-chemical code PRODIMO (Woitke et al., 2009; Woitke et al., 2016). This code self-consistently solves for the local radiation field, temperature, and chemical abundances for a given fixed axis-symmetric two-dimensional dust and gas density structure, and produces synthetic SEDs (Thi et al., 2011) and spectral lines (Woitke et al., 2011). To model the density structure, we use the rotating in-falling envelope model of Ulrich (1976) which has the mass of the protostar, $M_{\mathrm{*}}$, and the mass infall rate of the envelope, $\dot{M}_{\mathrm{inf}}$, as input parameters (see Rab et al., 2017, for further details).

For V883 Ori, Cieza et al. (2016) derived $\rm M_{*}=1.3\,M_{\odot}$ from ALMA $\rm C^{18}O(2-1)$ observations that resolve the inner Keplerian disk. These observations also constrain the disk’s inclination to $38\arcdeg$ and the systemic velocity to $4.3\,\mathrm{km\,s^{-1}}$. We do not model the disk structure as the APEX observations are not sensitive to the disk emission because of beam dilution. Our model also includes an empty outflow cavity with an opening angle of $75\arcdeg$ and has an extension of $r=18000\,\mathrm{au}$.

The outburst of V883 Ori is represented by a point-like source with $\rm L=200\,L_{\odot}$ in the center of the structure (see also, Rab et al., 2017). The radial temperature and density profiles through the midplane of our models are shown in Fig. 3. These profiles are very similar to the Class I outburst models of Rab et al. (2017) and the 1D outburst models of Visser et al. (2015) for embedded targets. To compare our model to the data, we produced synthetic line cubes, convolved them with a Gaussian beam (size according to the data), and re-binned them to the velocity resolution of the observations. Further details on the modeling procedure and a comparison to the dust modelling of Furlan et al. (2016) are presented in Appendix B.

We note that our model does not cover to the entire area mapped by the observations, but does cover the region traced by $\mathrm{CO(6-5)}$ and $\mathrm{{}^{13}CO(3-2)}$ which show approximately symmetric emission in the lower-velocity component (see Fig. 2). $\mathrm{CO(4-3)}$ and $\mathrm{CO(3-2)}$ show also asymmetric structures especially towards the north-east (Fig. 2), which we cannot account for in our symmetric envelope model. Therefore, we only compare our model results within a radius of $10000\,\mathrm{au}$. However, the observations might still pick up emission along the line of sight from structures that may or may not be part of V883 Ori and are not included in our models. This could lead to our model overestimating the flux from the inner region.

Using our model framework, we can test if a simple envelope structure with the above constraints is consistent with the APEX data. In Fig. 1, we show the results for two models with $\dot{M}_{\mathrm{inf}}=2\times 10^{-6}\,\rm M_{\odot}\,yr^{-1}$ (Model 1) and $\dot{M}_{\mathrm{inf}}=1\times 10^{-6}\,\rm M_{\odot}\,yr^{-1}$ (Model 2). Model 1 is a good match for the lower-velocity peak of the 12CO(4–3) and 12CO(3–2) lines. The modeled velocity FWHM and peak are in good agreement with the data. This indicates that the lower-velocity peak of the profile is indeed consistent with an in-falling envelope positioned on top of the disk of V883 Ori. Our model also shows that a single envelope structure is not sufficient to explain all the peaks in the spectra. Assuming that self-absorption produces the double-peaked profile $\rm v_{sys}\approx 5\,km\,s^{-1}$ is required, which is inconsistent with the value derived from ALMA disk observations. Furthermore, the FWHM of this single component ($\rm\sim 2\,km\,s^{-1}$) is too broad for a single in-falling envelope structure and the two components are not symmetric (e.g., the FWHM of the right component is a factor of 1.5 too wide).

As seen in Fig. 1, Model 1 underestimates the 13CO(3–2) by a factor of $\sim 5$ and overestimates the 12CO(6–5) by a factor of $\sim 2$. Model 2 on the other hand is a better match for 12CO(6–5), but underestimates all other lines. From the models, we can determine the main line emitting volume (i.e., the origin of 50% of the emission) for each individual line. For all lines, this region is within the inner $10000\,\mathrm{au}$. For these emitting regions we can estimate averaged quantities such as the particle number density $n$, temperature $T$ and CO molecular abundances $\epsilon$. For Model 1, we find for the 12CO(3–2) line: $n\!\approx\!8.7\!\times\!10^{3}\,\mathrm{cm^{-3}}$, $\epsilon\!\approx\!2.5\!\times\!10^{-5}$, $T\!\approx\!26.0\,\mathrm{K}$; for 13CO(3–2): $n\!\approx 2.7\!\times\!10^{4}\,\mathrm{cm^{-3}}$, $\epsilon\!\approx\!1.1\!\times\!10^{-4}$, $T\!\approx\!33.3\,\mathrm{K}$; for 12CO(4–3): $n\!\approx\!1.3\!\times\!10^{4}\,\mathrm{cm^{-3}}$, $\epsilon\!\approx\!4.6\!\times\!10^{-5}$, $T\!\approx\!28.5\,\mathrm{K}$ and for 12CO(6–5): $n\!\approx\!2.7\!\times\!10^{4}\,\mathrm{cm^{-3}}$, $\epsilon\!\approx\!1.1\!\times\!10^{-4}$, $T\!\approx\!33.7\,\mathrm{K}$. This shows that in the model 13CO(3–2) and 12CO(6–5) trace similar regions that are warmer and denser compared to the other two lines. We explored varying other parameters in the model (e.g., central luminosity, dust opacities) but did not find a significantly better model that reproduces the data. This may be because our envelope density structure is too simple or our model does not correctly predict the temperature in the inner regions. All of the models tested under-predict the 13CO(3–2) significantly, which might indicate that 13CO(3–2) does not only trace the envelope but also the SE/NW outflows reported by Ruíz-Rodríguez et al. (2017a). We discuss possible interpretations of this in Sect. 4.

The envelope masses in our models are approximately 0.4 and $0.2\,\rm M_{\odot}$ for Model 1 and 2, respectively. The mass is determined by $\dot{M}_{\mathrm{inf}}$ and the outer radius of the structure. The outer radius, however, is not well constrained. An $r=12000\,\mathrm{au}$ model matches the data as well, but the mass of the envelope would then only be $\approx 0.14\,\rm M_{\odot}$. Our mass estimates are similar to the derived envelope masses for Class I objects reported by Kóspál et al. (2017) in their sample of outbursting sources. In addition, Furlan et al. (2016) derived $T_{\mathrm{bol}}=492\,\rm K$, which further indicates that V883 Ori is a Class I object. Our simple model is consistent with the Class I evolutionary stage (see also Appendix B), but current mass estimates of the envelope and the disk are uncertain and require more detailed modeling of multi-wavelength data to precisely constrain the evolutionary stage of V883 Ori.

4 Discussion

Our models are consistent with a large ($r\sim 10000$ au) envelope-like structure traced by the lower-velocity feature centered at the systemic-velocity of the CO lines in the APEX data (Fig. 1). We note that the 12CO APEX data show strong emission at radial distances $\gtrsim 10000\,\mathrm{au}$, which may be related to the central envelope or outflow material. A larger envelope model, however, would be inconsistent with the more compact 12CO(6–5) emission. A model that matches 13CO would require twice as massive of an envelope, but is inconsistent with the 12CO emission. This more massive envelope produces a 13CO line profile that is broader than observed and the line becomes optically thick.

Considering the high-resolution ALMA data of Ruíz-Rodríguez et al. (2017b) and the APEX observations in context of our model implies that the 13CO(3–2) line is under-predicted because our model does not include outflows. The ALMA data trace the inner $<10{{}^{\prime\prime}}$ region ($\sim 5000\,\mathrm{au}$) and show a prominent outflow structure. However, due to limited baseline coverage, the ALMA observations are prone to spatial filtering. For example, the red-shifted outflow visible in 13CO ALMA observations is not seen in 12CO. The APEX observations do not resolve this inner region but clearly show large scale envelope structure ($r\gtrsim 5000\,\mathrm{au}$), which is likely spatially filtered out in the ALMA data. The 13CO(3–2) line is optically thin in our model with $\tau_{\mathrm{line}}\lesssim 1$, but the modeled emission is nearly a factor of five too weak compared to the data. The main isotopologue lines have optical depths $\tau>10$ in our model, consistent with the observations of Ruíz-Rodríguez et al. (2017b) where the inner disk structure becomes visible in 13CO(2–1) but not in 12CO(2–1). If the lower-velocity 13CO(3–2) component in the envelope is indeed optically thin, then the APEX observation of 13CO(3–2) might be more dominated by outflow emission than the other optically thick lines. This indicates that our model is consistent with the ALMA data but fails to reproduce the 13CO emission as it traces predominantly outflows not included in our modelling. Regardless, the observed CO emission of V883 Ori is likely a superposition of outflow and envelope emission as the characteristic outflow velocity is only $0.65\,\mathrm{km\,s^{-1}}$ (Ruíz-Rodríguez et al., 2017b), hence most of the outflow emission is likely in the lower-velocity component that we model as a pure envelope structure. However, as the ALMA observations miss the large scale emission and the APEX data do not resolve the outflow region, it is not possible with the currently available data to clearly disentangle and quantify those two contributions.

If V883 Ori’s total luminosity ($\rm 200~{}L_{\odot}$) is indeed dominated by accretion, then the current mass accretion rate onto the star is $\rm\dot{M}_{accr}\sim 1-2\times 10^{-5}~{}M_{\odot}~{}{yr^{-1}}$ (assuming the protostar is $\sim 2.3\,\rm R_{\odot}$ and $\sim 1.3\,\rm M_{\odot}$ as determined by Cieza et al., 2016). ALMA observations suggest an outflow rate of $\rm\dot{M}_{outf}\sim 2\times 10^{-6}~{}M_{\odot}~{}{yr^{-1}}$ (Ruíz-Rodríguez et al., 2017b) and the mass-infall rate onto the envelope from our model is $\rm\dot{M}_{inf}\approx 1-2\times 10^{-6}~{}M_{\odot}~{}{yr^{-1}}$. This is consistent with the episodic accretion scenario, where $\rm\dot{M}_{accr}<\dot{M}_{inf}$ in the quiescent stage allowing material to accumulate in the disk and trigger the outburst. Consequently, $\rm\dot{M}_{accr}$ abruptly increases yielding $\rm\dot{M}_{accr}>\dot{M}_{inf}$. The observed $\rm\dot{M}_{outf}/\dot{M}_{accr}$ ratio is consistent with the canonical value of $\sim 0.1$ (Pelletier & Pudritz, 1992). However, the high value of $\rm\dot{M}_{outf}$ suggests the outflow is a result of an accretion outburst, as otherwise $\rm\dot{M}_{outf}/\dot{M}_{accr}\gtrsim 1$, assuming that the quiescent accretion rate is $\rm\lesssim\dot{M}_{inf}$. Ruíz-Rodríguez et al. (2017b) constrained the dynamical timescale of the outflow to $\sim 10^{4}\,\mathrm{yr}$. This timescale implies the outflow structure may not be related to the ongoing outburst. We infer that the outflow may instead be a consequence of a previous outburst from $\gtrsim 10^{4}\,\mathrm{yr}$ ago with properties similar to the current outburst. Such a scenario is consistent with episodic accretion models (e.g., Vorobyov & Basu, 2015) and derived burst frequencies (e.g., Jørgensen et al., 2015). This implies the observed combination of envelope, outflow, and diffuse material may be predominately shaped by multiple outbursts.

Conversely, if the outflow/inflow rates are taken at face value, and assuming the quiescent accretion rate is well represented by the mass infall rate in the envelope, we find the ratio to be of order unity, much higher than the canonical value of $\sim 0.1$ (Pelletier & Pudritz, 1992). This implies that either the mass outflow rate has increased significantly or the mass infall rate has decreased significantly. The discrepancy could, however, arise from a deviation from the assumed morphology or contamination from the surrounding molecular cloud. As can be seen in Fig. 2, the 12CO(6–5) line traces a more compact region and does not show as much large scale emission as the other 12CO lines. Assuming that 12CO(6–5) is a better tracer of the envelope implies that the other 12CO lines may be contaminated by large scale emission from the molecular cloud. This is consistent with Model 1, which better represents the large scale structure, but over-predicts the flux in the central region.

High spatial and velocity resolution maps are necessary to further constrain the envelope parameters and obtain a single model that reproduces all the molecular species. Single pointings with ALMA (e.g., Cieza et al., 2016; Ruíz-Rodríguez et al., 2017b) will resolve out most of the large scale structure in the envelope, therefore $\gtrapprox 30^{\prime\prime}$ maps with uniform sensitivity are necessary to ensure large scale structure is preserved. Furthermore, the spatial resolution of the APEX data is not sufficient to see a clear difference in the shapes of the spectral profiles that might indicate the different origin of the 12CO and 13CO lines. Having a combination of high spatial and velocity resolution over a large area with ALMA and future facilities such as the ngVLA (e.g., White et al., 2018) will allow for the detection of any asymmetries in the envelope as well as accurately disentangle the contribution from inflows, outflows, and disk emission at all velocities.

5 Summary

We presented APEX FLASH*+* and SEPIA B9 science verification observations of the young outbursting star V883 Ori. We mapped and modeled the 12CO(6–5), 12CO(4–3), 12CO(3–2), and 13CO(3–2) lines of the envelope. Each line was resolved and detected at a high high signal-to-noise. We find the envelope is well represented by a spherical model with a $\gtrsim 10000\,\mathrm{au}$ radius, a mass of $\rm M_{env}\approx 0.2-0.4\,M_{\odot}$, and a mass-infall rate of $\rm\dot{M}_{inf}=1-2\times 10^{-6}\,M_{\odot}\,{yr^{-1}}$. The model does not, however, reproduce 13CO(3–2) and we conclude that this line is likely dominated by emission from the outflow structure apparent in high-resolution ALMA data (Ruíz-Rodríguez et al., 2017a). The large-scale structure in the 12CO lines, that we interpret as an envelope, is likely spatially filtered out in the ALMA observations. The combination of the inferred accretion rate and the observed envelope and outflow structures may be predominately shaped by multiple outbursts consistent with an episodic accretion scenario.

We thank the anonymous referee for feedback that improved the manuscript. We thank A. Hacar for discussions on the interpretation of the data. This publication is based on data acquired under programs 0102.F-9802(A) and 098.F-9505(A) with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the Max-Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. This project received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program, grant agreement No 716155 (SACCRED). CH.R. acknowledges funding from the Netherlands Organisation for Scientific Research (NWO) TOP-1 grant project number 614.001.552. CH.R. and E.V. acknowledge funding by the Austrian Science Fund (FWF): project number I2549-N27. M.A. and A.P. acknowledge funding by the Swiss National Science Foundation (SNSF): project number: 200021L_163172. T.Cs. acknowledges support from the Deutsche Forschungsgemeinschaft, DFG via the SPP (priority program) 1573 “Physics of the ISM”. We thank the Center for Information Technology of the University of Groningen for their support and access to the Peregrine high performance computing cluster.

Appendix A Velocity channel maps

The velocity channels maps for each species are shown in Figs. 4, 5, 6, and 7 below. As the 12CO(6–5) line has the lowest velocity resolution of the line observed with $\rm 0.20~{}km~{}s^{-1}$, the other lines are re-binned so that every species has the same velocity resolution. In each figure, $3\sigma$ contours are shown where $\sigma$ is the $\sigma_{RMS}$ of the line-free channels. For each species, we find the emission around $\rm 4.1~{}km~{}s^{-1}$ to be more spatially concentrated around the location of V883 Ori, which is indicated by the black plus in the center of each image. This is consistent with emission from an envelope structure as detailed in Sec. 3. The higher-velocity emission peaks seen in Fig. 1 are significantly less spatially concentrated for all species and may be consistent with a combination of diffuse emission from the molecular cloud and outflow material.

Appendix B Modelling procedure and Comparison to Furlan et al. (2016)

In this section we provide further details on our modelling procedure for the APEX data and a comparison to models of the continuum emission of V883 Ori. For the description of the modelling code see Sect 3.

For our modeling, we at first computed a small grid of models varying the mass infall rate $\dot{M}_{\mathrm{inf}}$ in the range of $5\times 10^{-7}-5\times 10^{-5}\,\rm M_{\odot}~{}yr^{-1}$ and the central luminosity $L$ in the range of $100-400\,\rm L_{\odot}$. We next sampled parameters around the best matching model of the initial grid varying only the $\dot{M}_{\mathrm{inf}}$. Our models do not include a disk component as the APEX data is not sensitive to the disk emission (see Sec. 3) and to minimize the required computational time. From this exploration of parameter space, we chose the Model 1 and Model 2, as described in Sec. 3.

For the two chosen models, we explored including a disk component with $M_{\mathrm{disk}}=0.2\,\rm M_{\odot}$ and a radius of $200\,\mathrm{au}$. The impact of the disk component on the synthetic line observables we compare to the APEX data is not significant and does not affect any of our conclusions. However, the disk component is the dominant source of continuum emission at wavelengths $\gtrsim 100\,\mathrm{\mu m}$. We note that it is outside the scope of this paper to fit all the available data (including the continuum), and instead aim to provide a reasonable model for the gaseous envelope structure of V883 Ori. Nevertheless, it is useful to compare our modelling results also to pure continuum models.

Furlan et al. (2016) modeled the SED of a large number of protostars in Orion also including V883 Ori (HOPS-376). Their modelling approach also assumes a fixed rotating envelope structure and a disk component. However, they note that they did not vary the disk parameters in their modelling and assumed that the data is not sensitive to the disk. In Table 2 we compare $L$, $\dot{M}_{\mathrm{inf}}$, and the continuum fluxes at $350\,\mathrm{\mu m}$ and $870\,\mathrm{\mu m}$. To derive $\dot{M}_{\mathrm{inf}}$ from Furlan et al. (2016), we used their Eq. 2 with $M_{\mathrm{*}}=1.3\,\rm M_{\odot}$. Our results are in good agreement, however their $\dot{M}_{\mathrm{inf}}$ is $2-4$ times that of our models. In Furlan et al. (2016), $\dot{M}_{\mathrm{inf}}$ is not a direct input parameter but depends on the stellar mass, which was fixed to $M_{*}=0.5\,\rm M_{\odot}$, and results in $\dot{M}_{\rm inf}\approx 2.5\times 10^{-6}\,\rm M_{\odot}\,{yr^{-1}}$. From Table 2, it is also apparent that in our model the continuum fluxes are dominated by the disk. Furlan et al. (2016) do not discuss if this is also the case for their V883 Ori model.

Furlan et al. (2016) classified V883 Ori as a flat spectrum source and derived $L_{\mathrm{bol}}=217.88\,\mathrm{L_{\odot}}$ and $T_{\mathrm{bol}}=492\,\mathrm{K}$ from the observational data considered for their modelling. We note that a bolometric luminosity of $400\,\mathrm{L_{\mathrm{\odot}}}$ (Strom & Strom, 1993) is commonly adopted for V883 Ori, however it is not apparent how Strom & Strom (1993) derived this value. It might be the case that $L_{\mathrm{bol}}$ of V883 Ori dropped over the last decade. Regardless, our model is consistent with the values reported by Furlan et al. (2016).