<figure><img src="image_1.png" /></figure>

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

J. A. White¹, Á. Kóspál¹,², C. Rab³, P. Ábrahám¹, F. Cruz-Sáenz de Miera¹, T. Csengeri⁴, O. Fehér¹,⁵, R. Güsten⁴, T. Henning², E. Vorobyov⁶,⁷, M. Audard⁸, and A. Postel⁸

¹Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences, Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary; jacob.white@csfk.mta.hu

²Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

³Kapteyn Astronomical Institute, University of Groningen, Postbus 800, 9700 AV Groningen, The Netherlands

⁴Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

⁵IRAM, 300 Rue de la Piscine, F-38046 Saint Martin d'Hères, France

⁶University of Vienna, Department of Astrophysics, Türkenschanzstr. 17, 1180 Wien, Austria

⁷Research Institute of Physics, Southern Federal University, Stachki 194, Rostov-on-Don, 344090, Russia

⁸Department of Astronomy, University of Geneva, Ch. d'Ecogia 16, 1290 Versoix, Switzerland

Received 2019 January 25; revised 2019 April 10; accepted 2019 April 12; published 2019 May 20

## 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 infalling 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 ¹²CO(4–3), ¹²CO(3–2), and ¹³CO(3–2) lines with the FLASH⁺ instrument and the ¹²CO(6–5) line with the SEPIA instrument. We detected high signal-to-noise emission extending out to radii >10,000 au and calculated integrated fluxes of 1100 Jy km s⁻¹ for ¹²CO(6–5), 2400 Jy km s⁻¹ for ¹²CO(4–3), 1600 Jy km s⁻¹ for ¹²CO(3–2), and 450 Jy km s⁻¹ for ¹³CO(3–2). We used the thermochemical code PRODIMO to test several models and find the data are best described by an envelope structure with Mₐvil ≈ 0.2–0.4 M☉ and a mass-infall rate of M'ₐinf = 1–2 × 10⁻⁶ M☉ yr⁻¹. We infer that the observed envelope and outflow structure around V883 Ori could be caused by multiple outbursts, consistent with episodic accretion.

**Key words:** circumstellar matter – stars: individual (V883 Orionis) – stars: pre-main sequence – submillimeter: stars

## 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×. 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 thousands 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 µ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 ± 0.1 M☉ (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 ± 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
<table><caption>**Table 1**
Summary of Our APEX Observations</caption><thead><tr><th>Species</th><th>Observation Date<br>YYYY MMM DD</th><th>Frequency<br>(GHz)</th><th>Beam<br>(")</th><th>σ<sub>rms</sub><br>(Jy beam<sup>-1</sup>)</th><th>Ch. Width<br>(km s<sup>-1</sup>)</th><th>Peak<br>(Jy beam<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>Env. Flux<br>(Jy km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td><sup>12</sup>CO(6–5)</td><td>2018 Oct 15/21</td><td>691.473</td><td>9.6</td><td>2.2</td><td>0.20</td><td>1100 ± 45</td><td>0.38 ± 0.02</td><td>1100 ± 170</td></tr><tr><td><sup>12</sup>CO(4–3)</td><td>2016 Sep 6</td><td>461.041</td><td>14.4</td><td>1.4</td><td>0.050</td><td>2800 ± 20</td><td>0.35 ± 0.01</td><td>2400 ± 240</td></tr><tr><td><sup>12</sup>CO(3–2)</td><td>2016 Sep 6</td><td>345.796</td><td>19.2</td><td>0.53</td><td>0.033</td><td>1800 ± 10</td><td>0.36 ± 0.01</td><td>1600 ± 160</td></tr><tr><td><sup>13</sup>CO(3–2)</td><td>2016 Sep 6</td><td>330.588</td><td>20.1</td><td>0.81</td><td>0.033</td><td>680 ± 6</td><td>0.26 ± 0.01</td><td>450 ± 45</td></tr></tbody></table>

**Note.** The line $\sigma_{rms}$ is from the line-free channels, rebinned to 1.0 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 Figure 1. The uncertainties on the peak and FWHM are the $1\sigma$ uncertainties from the covariance matrix in the Gaussian fitting. The uncertainties on the envelope flux are $\sigma_{rms}$, model fitting uncertainty, and absolute flux calibration uncertainty added in quadrature. The absolute flux calibration uncertainty is taken as 15% for $^{12}\text{CO}(6-5)$ and 10% for all other lines and is the dominate source of flux uncertainty. The beam efficiencies and Jy K$^{-1}$ conversion factors for each observational setup are given in Section 2.

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 indicates 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 Section 2, we describe the data reduction and observational results. In Section 3, we detail the envelope model. In Section 4, we discuss the implications of the model. In Section 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 $^{12}\text{CO}(3-2)$, $^{13}\text{CO}(3-2)$, and $^{12}\text{CO}(4-3)$ lines on 2016 September 6 (Project 098.F-9505(A), PI: Th. Henning), and the SEPIA B9 receiver to measure the $^{12}\text{CO}(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. = $0.5^{h}38^{m}18^{s}09$7 and decl. = $-0.7^{\circ}02'25''/96$. The main parameters of our observations are summarized in Table 1.

FLASH$^{+}$ is a dual-frequency heterodyne receiver, operating simultaneously in the 345 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 $^{13}\text{CO}(3-2)$ and $^{12}\text{CO}(3-2)$ lines at 330.588 GHz and 345.796 GHz, respectively. The higher-frequency channel was tuned to the $^{12}\text{CO}(4-3)$ line at 461.041 GHz in USB. We used the XFFTS backends, providing a 38 kHz (0.033 km s$^{-1}$) spectral resolution for both the 3–2 lines and 76 kHz (0.050 km s$^{-1}$) for the 4–3 line. We obtained 120" × 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 built (Baryshev et al. 2015) to the specifications of the ALMA Band 9 receiver (600–722 GHz). Currently, the XFFTS backend records 4 GHz bandwidth in each sideband. The USB was centered on the $^{12}\text{CO}(6-5)$ line at 691.473 GHz with 450 kHz (0.20 km s$^{-1}$) resolution. We obtained 150" × 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 package<fnref n="9" /> 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 2019, private communication) at 352, 464, and 691 GHz, respectively. The data were converted to Jy using 41, 48, and 79 Jy K$^{-1}$ at 352, 464, and 691 GHz, respectively. The telescope’s beam size is $20''/1, 19''/2, 14''/4$, and $9''/6$ for the $^{13}\text{CO}(3-2)$, $^{12}\text{CO}(3-2)$, $^{12}\text{CO}(4-3)$, and $^{12}\text{CO}(6-5)$ lines, respectively. The $\sigma_{rms}$ for each line is calculated by rebinning each spectrum to a 1 km s$^{-1}$ channel width, then calculating the rms 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 and 464 GHz and as 15% at 691 GHz.

### 2.1. Observational Results

The spectra of each molecular species are presented in Figure 1, spatially integrated over a region of 10,000 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 two peaks in the emission seen in $^{12}\text{CO}(6-5)$ and $^{13}\text{CO}(3-2)$ and three peaks seen in $^{12}\text{CO}(4-3)$ and $^{12}\text{CO}(3-2)$. The best-fit positions of the lower-velocity peak in each spectrum are all centered at ~4.1 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 nonspherical 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, Figure 4 shows the $^{12}\text{CO}(6-5)$ line channel map. In this figure, the emission from the envelope (lower-velocity peak in the spectra) is largely contained within radii <10,000 au at velocities <5.0 km s$^{-1}$ while the emission at velocities >5.0 km s$^{-1}$ is more consistent with diffuse

<footnote n="9">⁹ http://www.iram.fr/IRAMFR/GILDAS/</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Integrated spectra of all the CO species over a 10,000 au radius aperture. The black lines are the observed spectra. The red line is Model 1, with a mass-infall rate of $\dot{M}_{\text{inf}} = 2 \times 10^{-6} M_{\odot} \text{yr}^{-1}$. The blue line is Model 2, with $\dot{M}_{\text{inf}} = 1 \times 10^{-6} M_{\odot} \text{yr}^{-1}$. The lower-velocity peak for each species is from the envelope. The additional peaks seen for all species at ~5.5 km s⁻¹ and for ¹²CO(4–3) and ¹²CO(3–2) at ~7.5 km s⁻¹ are likely due to either outflow material or contamination from the surrounding molecular cloud. This is further illustrated in the velocity channel maps in the Appendix (Figures 4–7).</figcaption></figure>

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 toward 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 the systemic velocity (e.g., 3.4 and 4.6 km s⁻¹ in Figure 1) suggests infalling material, as higher velocities are expected closer to the center (i.e., freefall velocity $\propto 1/\sqrt{r}$). To calculate the peak, FWHM, and total envelope flux for each species, we first integrated each spectrum over a 10,000 au radius aperture centered on V883 Ori. A multicomponent 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 Figure 1 are summarized in Table 1. The integrated intensity maps (0th-order moment) for each species are shown in Figure 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 thermochemical code PRODIMO (Woitke et al. 2009, 2016). This code self-consistently solves for the local radiation field, temperature, and chemical abundances for a given fixed axisymmetric 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 infalling envelope model of Ulrich (1976), which has the mass of the protostar, $M_*$, and the mass-infall rate of the envelope, $\dot{M}_{\text{inf}}$, as input parameters (see Rab et al. 2017 for further details).

For V883 Ori, Cieza et al. (2016) derived $M_* = 1.3 M_{\odot}$ from ALMA C¹⁸O(2 − 1) observations that resolve the inner Keplerian disk. These observations also constrain the disk’s inclination to 38° and the systemic velocity to 4.3 km s⁻¹. 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° and has an extension of $r = 18,000$ au.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Integrated intensity maps for all of the observed CO lines. Each spectrum was integrated from 2.4 to 4.8 km s⁻¹, which corresponds to the lower-velocity peak as seen in Figure 1. The position of the star is denoted by the black plus and the beam size by the black circle. The 3σ contours are indicated by the black lines. We note that the contours are not visible for ¹²CO(4–3) and ¹²CO(3–2) as the envelope emission is detected at a very high signal-to-noise ratio.</figcaption></figure>

The outburst of V883 Ori is represented by a point-like source with $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 Figure 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 rebinned them to the velocity resolution of the observations. Further details on the modeling procedure and a comparison to the dust modeling 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 CO(6–5) and ¹³CO(3–2), which show approximately symmetric emission in the lower-velocity component (see Figure 2). CO(4–3) and CO(3–2) show also asymmetric structures especially toward the northeast (Figure 2), which we cannot account for in our symmetric envelope model. Therefore, we only compare our model results within a radius of 10,000 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 Figure 1, we show the results for two models with $\dot{M}_{\text{inf}} = 2 \times 10^{-6} M_\odot \text{yr}^{-1}$ (Model 1) and $\dot{M}_{\text{inf}} = 1 \times 10^{-6} M_\odot \text{yr}^{-1}$ (Model 2). Model 1 is a good match for the lower-velocity peak of the ¹²CO(4–3) and ¹²CO(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 infalling 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 $v_{\text{sys}} \approx 5 \text{ km s}^{-1}$ is required, which is inconsistent with the value derived from ALMA disk observations. Furthermore, the FWHM of this single component ($\sim 2 \text{ km s}^{-1}$) is too broad for a single infalling 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 Figure 1, Model 1 underestimates the ¹³CO(3–2) by a factor of $\sim 5$ and overestimates the ¹²CO(6–5) by a factor of $\sim 2$.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Radial number density (top panel) and temperature (bottom panel) profiles in the midplane of the two-dimensional envelope models. Model 1 is indicated by the red curve, and Model 2 is indicated by the blue curve.</figcaption></figure>

Model 2 on the other hand is a better match for ¹²CO(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 10,000 au. For these emitting regions we can estimate averaged quantities such as the particle number density *n*, temperature *T*, and CO molecular abundances *ϵ*. For Model 1, we find for the ¹²CO(3–2) line: *n* ≈ 8.7 × 10³ cm⁻³, *ϵ* ≈ 2.5 × 10⁻⁵, *T* ≈ 26.0 K; for ¹³CO(3–2): *n* ≈ 2.7 × 10⁴ cm⁻³, *ϵ* ≈ 1.1 × 10⁻⁴, *T* ≈ 33.3 K; for ¹²CO(4–3): *n* ≈ 1.3 × 10⁴ cm⁻³, *ϵ* ≈ 4.6 × 10⁻⁵, *T* ≈ 28.5 K; and for ¹²CO(6–5): *n* ≈ 2.7 × 10⁴ cm⁻³, *ϵ* ≈ 1.1 × 10⁻⁴, *T* ≈ 33.7 K. This shows that in the model ¹³CO(3–2) and ¹²CO(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 underpredict the ¹³CO(3–2) significantly, which might indicate that ¹³CO(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 Section 4.

The envelope masses in our models are approximately 0.4 and 0.2 $M_\odot$ for Models 1 and 2, respectively. The mass is determined by $\dot{M}_{\text{inf}}$ and the outer radius of the structure. The outer radius, however, is not well constrained. An $r = 12,000$ au model matches the data as well, but the mass of the envelope would then only be ≈0.14 $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_{\text{bol}} = 492$ 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 multiwavelength data to precisely constrain the evolutionary stage of V883 Ori.

## 4. Discussion

Our models are consistent with a large ($r \sim 10,000$ au) envelope-like structure traced by the lower-velocity feature centered at the systemic velocity of the CO lines in the APEX data (Figure 1). We note that the ¹²CO APEX data show strong emission at radial distances $\gtrsim 10,000$ au, which may be related to the central envelope or outflow material. A larger envelope model, however, would be inconsistent with the more compact ¹²CO(6–5) emission. A model that matches ¹³CO would require twice as massive of an envelope, but is inconsistent with the ¹²CO emission. This more massive envelope produces a ¹³CO 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 ¹³CO(3–2) line is underpredicted because our model does not include outflows. The ALMA data trace the inner $<10^{\prime\prime}$ region (~5000 au) and show a prominent outflow structure. However, due to limited baseline coverage, the ALMA observations are prone to spatial filtering. For example, the redshifted outflow visible in ¹³CO ALMA observations is not seen in ¹²CO. The APEX observations do not resolve this inner region but clearly show large-scale envelope structure ($r \gtrsim 5000$ au), which is likely spatially filtered out in the ALMA data. The ¹³CO(3–2) line is optically thin in our model with $\tau_{\text{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 ¹³CO(2–1) but not in ¹²CO(2–1). If the lower-velocity ¹³CO(3–2) component in the envelope is indeed optically thin, then the APEX observation of ¹³CO(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 ¹³CO emission as it predominantly traces outflows not included in our modeling. 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 km s⁻¹ (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 ($200 L_\odot$) is indeed dominated by accretion, then the current mass accretion rate onto the star is $\dot{M}_{\text{accr}} \sim 1-2 \times 10^{-5} M_\odot \text{ yr}^{-1}$ (assuming the protostar is $\sim 2.3 R_\odot$ and $\sim 1.3 M_\odot$ as determined by Cieza et al. 2016). ALMA observations suggest an outflow rate of
$\dot{M}_{\text{outf}} \sim 2 \times 10^{-6} M_{\odot} \text{yr}^{-1}$ (Ruíz-Rodríguez et al. 2017b) and the mass-infall rate onto the envelope from our model is $\dot{M}_{\text{inf}} \approx 1-2 \times 10^{-6} M_{\odot} \text{yr}^{-1}$. This is consistent with the episodic accretion scenario, where $\dot{M}_{\text{accr}} < \dot{M}_{\text{inf}}$ in the quiescent stage allowing material to accumulate in the disk and trigger the outburst. Consequently, $\dot{M}_{\text{accr}}$ abruptly increases yielding $\dot{M}_{\text{accr}} > \dot{M}_{\text{inf}}$. The observed $\dot{M}_{\text{outf}}/\dot{M}_{\text{accr}}$ ratio is consistent with the canonical value of $\sim 0.1$ (Pelletier & Pudritz 1992). However, the high value of $\dot{M}_{\text{outf}}$ suggests the outflow is a result of an accretion outburst, as otherwise $\dot{M}_{\text{outf}}/\dot{M}_{\text{accr}} \gtrsim 1$, assuming that the quiescent accretion rate is $\lesssim \dot{M}_{\text{inf}}$. Ruíz-Rodríguez et al. (2017b) constrained the dynamical timescale of the outflow to $\sim 10^4$ 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$ 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 Figure 2, the $^{12}$CO(6–5) line traces a more compact region and does not show as much large-scale emission as the other $^{12}$CO lines. Assuming that $^{12}$CO(6–5) is a better tracer of the envelope implies that the other $^{12}$CO 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 overpredicts 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, $\gtrsim 30''$ 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 $^{12}$CO and $^{13}$CO 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 $^{12}$CO(6–5), $^{12}$CO(4–3),

$^{12}$CO(3–2), and $^{13}$CO(3–2) lines of the envelope. Each line was resolved and detected at a high signal-to-noise ratio. We find the envelope is well represented by a spherical model with a $\gtrsim 10,000$ au radius, a mass of $M_{\text{env}} \approx 0.2–0.4 M_{\odot}$, and a mass-infall rate of $\dot{M}_{\text{inf}} = 1–2 \times 10^{-6} M_{\odot} \text{yr}^{-1}$. The model does not, however, reproduce $^{13}$CO(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 $^{12}$CO lines, which 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 Radioastromie, 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. C.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.C. 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.

**Facilities:** APEX (SEPIA, FLASH+).

**Software:** GILDAS, PRODIMO (Woitke et al. 2009, 2016; Rab et al. 2017).

## Appendix A
### Velocity Channel Maps

The velocity channel maps for each species are shown in Figures 4–7. As the $^{12}$CO(6–5) line has the lowest velocity resolution of the line observed with $0.20 \text{ km s}^{-1}$, the other lines are rebinned so that every species has the same velocity resolution. In each figure, $3\sigma$ contours are shown where $\sigma$ is the $\sigma_{\text{rms}}$ of the line-free channels. For each species, we find the emission around $4.1 \text{ 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 Section 3. The higher-velocity emission peaks seen in Figure 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.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** $^{12}$CO(6–5) velocity channel map. The selected channels highlight how the emission from V883 Ori’s envelope is predominately located at velocities <5 km s$^{-1}$ (the lower-velocity peak in all spectra in Figure 1), while the emission at velocities >5 km s$^{-1}$ is likely due to outflow material or diffuse emission from the molecular cloud. The 3σ contours are given by the black lines with σ = 30 Jy beam$^{-1}$ (0.38 K beam$^{-1}$). The position of the star is denoted by the black plus, and the beam size is denoted by the black circle.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** $^{12}$CO(4–3) velocity channel map with channels rebinned to 0.2 km s$^{-1}$. The selected channels highlight how the emission from V883 Ori's envelope is predominately located at velocities <5 km s$^{-1}$ (the lower-velocity peak in all spectra in Figure 1), while the emission at velocities >5 km s$^{-1}$ is likely due to outflow material or diffuse emission from the molecular cloud. The 3$\sigma$ contours are given by the black lines with $\sigma = 42$ Jy beam$^{-1}$ (0.88 K beam$^{-1}$). The position of the star is denoted by the black plus, and the beam size is denoted by the black circle.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** $^{12}$CO(3–2) velocity channel map with channels rebinned to 0.2 km s$^{-1}$. The selected channels highlight how the emission from V883 Ori's envelope is predominately located at velocities <5 km s$^{-1}$ (the lower-velocity peak in all spectra in Figure 1), while the emission at velocities >5 km s$^{-1}$ is likely due to outflow material or diffuse emission from the molecular cloud. The 3$\sigma$ contours are given by the black lines with $\sigma = 22$ Jy beam$^{-1}$ (0.54 K beam$^{-1}$). The position of the star is denoted by the black plus, and the beam size is denoted by the black circle.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** $^{13}$CO(3–2) velocity channel map with channels rebinned to 0.2 km s$^{-1}$. The selected channels highlight how the emission from V883 Ori’s envelope is predominately located at velocities <5 km s$^{-1}$ (the lower-velocity peak in all spectra in Figure 1), while the emission at velocities >5 km s$^{-1}$ is likely due to outflow material or diffuse emission from the molecular cloud. The 3$\sigma$ contours are given by the black lines with $\sigma = 18$ Jy beam$^{-1}$ (0.44 K beam$^{-1}$). The position of the star is denoted by the black plus, and the beam size is denoted by the black circle.</figcaption></figure>

# Appendix B

## Modeling Procedure and Comparison to Furlan et al. (2016)

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

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

For the two chosen models, we explored including a disk component with $M_{\text{disk}} = 0.2 M_{\odot}$ and a radius of 200 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$ µ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 modeling 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 modeling 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 modeling and assumed that the data are not sensitive to the disk. In Table 2 we compare $L$, $\dot{M}_{\text{inf}}$, and the continuum fluxes at 350 and 870 µm.
<table><caption>**Table 2**
Comparison of Our Models to the Continuum Model and Selected Photometric Fluxes of Furlan et al. (2016)</caption><thead><tr><th>Model</th><th><i>L</i><br>(L<sub>☉</sub>)</th><th><i>M</i><sub>inf</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>350 µm Flux<br>(Jy)</th><th>870 µm Flux<br>(Jy)</th></tr></thead><tbody><tr><td>Furlan et al.<br>(2016)</td><td>194.67</td><td>≈4 × 10<sup>-6</sup></td><td>11.3 ± 4.5</td><td>1.2 ± 0.2</td></tr><tr><td>Model 1</td><td>200</td><td>2 × 10<sup>-6</sup></td><td>3.4</td><td>0.2</td></tr><tr><td>Model 1 + disk</td><td>200</td><td>2 × 10<sup>-6</sup></td><td>14.3</td><td>1.6</td></tr><tr><td>Model 2</td><td>200</td><td>1 × 10<sup>-6</sup></td><td>1.9</td><td>0.1</td></tr><tr><td>Model 2 + disk</td><td>200</td><td>1 × 10<sup>-6</sup></td><td>12.5</td><td>1.46</td></tr></tbody></table>

**Note.** The parameters are described in Section 3 and Appendix B and in Furlan et al. (2016).

To derive $\dot{M}_{inf}$ from Furlan et al. (2016), we used their Equation (2) with $M_{*} = 1.3 M_{☉}$. Our results are in good agreement; however, their $\dot{M}_{inf}$ is 2–4 times that of our models. In Furlan et al. (2016), $\dot{M}_{inf}$ is not a direct input parameter but depends on the stellar mass, which was fixed to $M_{*} = 0.5 M_{☉}$, and results in $\dot{M}_{inf} \approx 2.5 \times 10^{-6} M_{☉} \text{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_{bol} = 217.88 L_{☉}$ and $T_{bol} = 492$ K from the observational data considered for their modeling. We note that a bolometric luminosity of 400 $L_{☉}$ (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_{bol}$ of V883 Ori dropped over the past decade. Regardless, our model is consistent with the values reported by Furlan et al. (2016).

ORCID iDs

J. A. White https://orcid.org/0000-0001-8445-0444

Á. Kóspál https://orcid.org/0000-0001-7157-6275

C. Rab https://orcid.org/0000-0003-1817-6576

P. Ábrahám https://orcid.org/0000-0001-6015-646X

F. Cruz-Sáenz de Miera https://orcid.org/0000-0002-4283-2185

T. Csengeri https://orcid.org/0000-0002-6018-1371

O. Fehér https://orcid.org/0000-0003-3453-4775

T. Henning https://orcid.org/0000-0002-1493-300X

E. Vorobyov https://orcid.org/0000-0002-6045-0359

References

Arce, H. G., & Sargent, A. I. 2006, ApJ, 646, 1070

Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, in Protostars & Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 387

Baryshev, A. M., Hesper, R., Mena, F. P., et al. 2015, A&A, 577, A129

Belitsky, V., Lapkin, I., Fredrixon, M., et al. 2018, A&A, 612, A23

Cassen, P., & Moosman, A. 1981, Icar, 48, 353

Cieza, L. A., Casassus, S., Tobin, J., et al. 2016, Natur, 535, 258

Cieza, L. A., Ruíz-Rodríguez, D., Perez, S., et al. 2018, MNRAS, 474, 4347

Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5

Green, J. D., Hartmann, L., Calvet, N., et al. 2006, ApJ, 648, 1099

Greene, T. P., Aspin, C., & Reipurth, B. 2008, AJ, 135, 1421

Güsten, R., Nyman, L. Å., Schilke, P., et al. 2006, A&A, 454, L13

Hacar, A., Tafalla, M., Forbrich, J., et al. 2018, A&A, 610, A77

Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207

Herbig, G. H. 1966, VA, 8, 109

Herbig, G. H. 1989, in ESO Workshop on Low Mass Star Formation and Pre- Main Sequence Objects, ed. B. Reipurth (Garching: ESO), 233

Jørgensen, J. K., Visser, R., Williams, J. P., et al. 2015, A&A, 579, A23

Klein, T., Ciechanowicz, M., Leinz, C., et al. 2014, ITTST, 4, 588

Kóspál, Á., Ábrahám, P., Csengeri, T., et al. 2017, ApJ, 836, 226

Kounkel, M., Covey, K., Suarez, G., et al. 2018, AJ, 156, 84

Lee, J. E., Lee, S., Baek, G., et al. 2019, NatAs, 3, 314

Menten, K. M., Reid, M. J., Forbrich, J., et al. 2007, A&A, 474, 515

Millan-Gabet, R., Monnier, J. D., Akeson, R. L., et al. 2006, ApJ, 641, 547

Pelletier, G., & Pudritz, R. E. 1992, ApJ, 394, 117

Pickering, E. C. 1890, Obs, 13, 80

Quanz, S. P., Henning, T., Bouwman, J., et al. 2007, ApJ, 668, 359

Quillen, A. C., Thorndike, S. L., Cunningham, A., et al. 2005, ApJ, 632, 941

Rab, C., Elbakyan, V., Vorobyov, E., et al. 2017, A&A, 604, A15

Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2017a, MNRAS, 466, 3519

Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2017b, MNRAS, 468, 3266

Strom, K. M., & Strom, S. E. 1993, ApJL, 412, L63

Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529

Thi, W.-F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711

Ulrich, R. K. 1976, ApJ, 210, 377

van't Hoff, M. L., Tobin, J. J., Trapman, L., et al. 2018, ApJL, 864, L23

Visser, R., Bergin, E. A., & Jørgensen, J. K. 2015, A&A, 577, A102

Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115

White, J. A., Audard, M., Ábrahám, P., et al. 2018, in ASP Conf. Ser. 517, Science with a Next Generation Very Large Array, ed. E. Murphy (San Francisco, CA: ASP), 177

Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67

Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103

Woitke, P., Riaz, B., Duchêne, G., et al. 2011, A&A, 534, A44