Evidence for magnetic activity at starbirth: a powerful X-ray flare from the Class 0 protostar HOPS 383

Context. Class 0 protostars represent the earliest evolutionary stage of solar-type stars, during which the majority of the system mass resides in an infalling envelope of gas and dust and is not yet in the central, nascent star. Although X-rays are a key signature of magnetic activity in more evolved protostars and young stars, whether such magnetic activity is present at the Class 0 stage is still debated. Aims. We aim to detect a bona fide Class 0 protostar in X-rays. Methods. We observed HOPS 383 in 2017 December in X-rays with the Chandra X-ray Observatory ($\sim$84 ks) and in near-infrared imaging with the Southern Astrophysical Research telescope. Results. HOPS 383 was detected in X-rays during a powerful flare. This hard (E>2 keV) X-ray counterpart was spatially coincident with the northwest 4 cm component of HOPS 383, which would be the base of the radio thermal jet launched by HOPS 383. The flare duration was $\sim$3.3 h; at the peak, the X-ray luminosity reached $\sim$4 x 1E31 erg s --1 in the 2-8 keV energy band, a level at least an order of magnitude larger than that of the undetected quiescent emission from HOPS 383. The X-ray flare spectrum is highly absorbed (NH $\sim$ 7 x 1E23 cm --2), and it displays a 6.4 keV emission line with an equivalent width of $\sim$1.1 keV, arising from neutral or low-ionization iron. Conclusions. The detection of a powerful X-ray flare from HOPS 383 constitutes direct proof that magnetic activity can be present at the earliest formative stages of solar-type stars.


Introduction
Low-mass objects that have evolved beyond the Class 0 stage are conspicuous X-ray emitters (Dunham et al. 2014), that is, Class I protostars with remnant envelopes and massive accretion disks, and Class II and III pre-main sequence stars with and without circumstellar disks (T Tauri stars). Their high luminosities (∼10 28−31 erg s −1 ) compared to the solar maximum (∼10 27 erg s −1 ) and their intense flaring activity in Xrays make them appear as extremely magnetically active young suns (Feigelson & Montmerle 1999;Güdel & Nazé 2009). In Class 0 protostars (André et al. 1993(André et al. , 2000, the hydrostatic core is deeply embedded within its envelope and the molecular cloud, making its detection difficult at most wavelengths (Giardino et al. 2007). The most deeply-embedded X-ray sources reported in star-forming regions (Hamaguchi et al. 2005;Getman et al. 2007 and references therein; Kamezaki et al. 2014) have evolved beyond the Class 0 stage or their bolometric luminosity is not accurate enough for a robust conclusion (Appendix A). Moreover, any non-thermal radio emission from putative magnetic activity at Class 0 protostars is very likely absorbed by the bases of their ionized outflows (Güdel 2002;Dzib et al. 2013).
The object of this study in the Orion Molecular Cloud 3 (OMC-3) was identified as a protostar by a Spitzer survey (Megeath et al. 2012) and included in the Herschel Orion protostar survey (HOPS; Furlan et al. 2016) as source 383. The submillimeter to bolometric luminosity ratio of HOPS 383, 1.4%, (Safron et al. 2015) well surpasses the (0.5%) threshold for bona fide Class 0 designation (André et al. 1993(André et al. , 2000. Notably, HOPS 383 is the first Class 0 protostar known to have undergone a mass-accretion-driven eruption (Safron et al. 2015), which peaked by 2008 and ended by 2017 September (Appendix B).

Observations
We observed HOPS 383 three times with the Chandra X-ray Observatory (Weisskopf et al. 2002) from 2017 December 13 to 14 with simultaneous near-infrared imaging on 2017 December 14 using the 4.1 m Southern Astrophysical Research (SOAR) telescope (Krabbendam et al. 2004) in Chile.
We used Chandra's Advanced CCD Imaging Spectrometer (ACIS-I) in the very faint mode with a frame time of 3.141 s (Appendix C.1). Chandra ACIS-I has on-axis a full-width halfmaximum (FWHM) angular resolution of 0 . . 492 pixels display logarithmically the 0.5-2 keV and 2-7 keV Xrays, respectively, which were detected from 2017 December 13 to 14 during three ACIS-I exposures, totaling 83 877 s. The logarithmic grayscale K-band image was obtained on 2017 December 14 with the Spartan infrared camera on the SOAR telescope. North is up and east is left. The pink pluses and diamond, the green pluses, and cyan crosses are near-and mid-infrared sources (Megeath et al. 2012), H 2 -shocked emission (Stanke et al. 2002), and X-ray sources, respectively. The white dashed box shows the field-of-view of Fig. 1b. Panel b: Spartan H 2 -filter image obtained on 2017 December 14 in linear grayscale. The white contours show the 12 CO blueshifted outflow of HOPS 383 (Feddersen et al. 2020). The blue pluses are radio sources (Galván-Madrid et al. 2015;Rodríguez et al. 2017). The blue and green arrows show the direction of the thermal-jet candidate and the proper motion of the H 2 emission knot SMZ 1-2B, respectively. The white dashed box shows the field-of-view of Fig. 1c.
spectral resolution of 261 eV at 6.4 keV. Data reduction is described in Appendix C.2. The SOAR Spartan infrared camera (Loh et al. 2012) is composed of 4 CCDs of 2, 048 × 2, 048 pixels. We selected the widefield mode with a single-detector field-of-view of 2 .25 × 2 .25 (i.e., 5 .04 × 5 .04 edge-to-edge) and a pixel scale of 0 . 0661 with the K-band and H 2 narrow-band filters. We used a five-point dithering pattern where, in the first exposure, HOPS 383 was put near the center of the northeastern detector (det3), which was then moved by 30 from its initial position sequentially toward the south, north, west, and east in the four other exposures. This sequence was repeated until the requested exposure was achieved (Appendix D.1). Data reduction is described in Appendix D.2.

Near-infrared results
The near-infrared nebulosity that was prominent during the outburst (Safron et al. 2015) is not detected in our SOAR K-band image (Fig. 1a), indeed, it vanished by 2015 December 30 (Fischer & Hillenbrand 2017). The H 2 -narrow filter image shows an emission knot of shocked molecular hydrogen at ∼15 southeast from the mid-infrared location of HOPS 383. We identify this H 2 source with the H 2 emission knot SMZ 1-2B (Stanke et al. 2002) that was observed on 1997 September 13, yielding a proper motion of ∼1 . 8±1 southeast in 20.25 yr (Fig. 1b). The corresponding velocity is ∼(95±53)×(sin i/sin 69.5 • ) −1 ×(d/420 pc) km s −1 , where d is the distance to HOPS 383 and i is the inclination an-gle (90 • corresponds to an edge-on view; Table 1 of Furlan et al. 2016). This is typical of the velocity values observed in the outflows from Class 0 protostars (Bontemps et al. 1996;Reipurth & Bally 2001). Assuming a constant velocity, the angular distance of ∼16 between HOPS 383 and this emission knot corresponds to a kinematical timescale of ∼180 ± 100 yr. The proper motion direction is consistent with the orientation of the outburst nebula, the 12 CO blueshifted outflow (Feddersen et al. 2020), and the position angle of the binary radio counterpart of HOPS 383 at 4 cm (Galván-Madrid et al. 2015). We note that the unresolved 3 cm counterpart (Rodríguez et al. 2017) corresponds to the faint northwest component at 4 cm, JVLA-NW, which suggests that JVLA-NW would be the base of the radio thermal jet launched by HOPS 383, whereas the brighter southeast component at 4 cm, JVLA-SE, would be an emission knot along this thermal jet (Fig. 1c).   dra exposure that we obtained in 2017 (Fig. 2a). The first photon is detected ∼2.9 h after the start of the observation and the X-ray burst duration is ∼3.3 h (Fig. 2b). We estimated the shape of the light curve by smoothing the event times of arrival (Appendix C.4). The count rate peaks to ∼5.6 counts ks −1 in ∼0.9 h after the first photon detection, and it decays gradually in ∼2.5 h until the last photon detection (red curve in Fig. 2c). This time variation, with rapid rise and slow decay, is similar to the typical light curves of magnetic flares from young stellar objects (YSOs; Getman et al. 2008). X-rays may also be emitted by accretion shocks at the protostellar surface and/or shocks in outflows (Güdel & Nazé 2009), but the former produces plasmas that are not bright and hot enough to emit the observed hard Xrays and the latter is constant on week timescales (Kastner et al. 2005), which is much longer than the hour variability that we observed. Therefore, this burst of X-rays indicates magnetic activity from HOPS 383 that is associated with JVLA-NW.

X-ray flare spectrum
We modeled the X-ray spectrum with isothermal coronal plasma emission suffering from intervening absorption (Appendix C.5). The best-fit model has an excess and a deficit relative to the observed spectrum below and above 6.4 keV, respectively, suggesting an additional component at 6.4 keV. This deviation is likely associated with the emission line arising from neutral or lowionization iron. Therefore, we included in the model a Gaussian emission line centered at 6.4 keV. We performed a Markov chain Monte Carlo analysis (MCMC) to compute the probability density functions of the physical parameters, from which we determined the median parameter values and 90% confidence intervals. We find that CXOU J053529 is highly embedded with a hydrogen column density (N H,X ) of 7.0 +3.5 −3.1 × 10 23 cm −2 . This result for N H is consistent with the large median energy of the 0.5-8 keV counts (MedE=5.402 keV) 1 , given the log N H versus MedE relations found for YSOs in the Orion nebula cluster (ONC; Fig. 8 of Feigelson et al. 2005) and M17 (Fig. 10a of Getman et al. 2010).
We detected the 6.4 keV emission line with a high probability of 98.9%. Our spectral analysis suggests that we have detected five photons from the Fe 6.4 keV emission line; these photons are likely the closest in energy to 6.4 keV. Indeed, we observed five photons with energies that are densely packed around the predicted emission feature at 6.4 keV (diamonds in Fig. 2b inset). Moreover, the spatial distribution of these Fe line photons is point-like (Fig. 2a), suggesting all of the observed Fe fluorescence emission arises from within the Chandra point spread function, that is to say, a radius from HOPS 383 that is lower than 105 au at the assumed distance to the protostar. The iron line equivalent width, EW = 1.1 +2.1 −0.9 keV, is large (even at the lower limit) relative to what was expected from the possible emission processes. The neutral iron can be ionized by photons or electrons that have energies larger than the Fe K-shell ionization potential of 7.112 keV. The contributions of these two ionization channels to the observed X-ray line in accreting young stars and Class I protostars (Güdel & Nazé 2009;Hamaguchi et al. 2005;Czesla & Schmitt 2010;Hamaguchi et al. 2010;Pillitteri et al. 2019) are not firmly settled, but photoionization is less energetically challenging (Czesla & Schmitt 2010). The plasma temperature is not well constrained mainly due to the low count rate, but a hot plasma temperature is favored (kT = 4.1 +13.8 −2.8 keV), which is consistent with a magnetic flare and the photoionization of iron. . Probing the gas and dust toward HOPS 383 with the X-ray spectrum. Panels a-c: 2D density distribution from small to large spatial scales and viewing angle from the model of HOPS 383 (Furlan et al. 2016). Dots indicate regions with dust. Panel d: cumulative hydrogen column density along the line-of-sight from the X-ray source. The dashed lines show the contribution of the accretion disk, the bipolar cavity, the infalling envelope, and the cloud. The data point with an error bar (90% C.I.) is the hydrogen column density that we obtained from the X-ray spectrum modeling.

X-ray flare energetics
The mean value of the absorption-corrected luminosity during the flare time interval is 1.7 +5.5 −1.1 × 10 31 erg s −1 in the 2-8 keV energy band. Multiplying this value by the flare peak to flare mean count rate ratio (2.39) provides the flare peak value of ∼4.2 +13.0 −2.6 × 10 31 erg s −1 . For comparison purposes, the maximum peak value of the 84 X-ray flares with typical shapes observed from YSOs in the ONC (Getman et al. 2008) is 8.3 × 10 32 erg s −1 in the 0.5-8 keV energy band; only ∼34% of these flares have peak values that are higher than the HOPS 383 one. During our observations, the non-detection of HOPS 383 outside the flare time interval implies that its quiescent absorption-corrected luminosity is lower than 2.0 × 10 30 erg s −1 in the 2-8 keV energy band (Appendix C.6). Therefore, the flare that we detected peaked at least 21 +65 −13 times above the quiescent level. The lower limit on the total energy released in X-rays during this flare, as computed from the mean luminosity times the flare duration, was 2.0 +6.6 −1.3 × 10 35 ergs.

Discussion
Starting from a model of HOPS 383 consisting of an infalling envelope with bipolar cavities carved by outflows and a small (5 au radius) accretion disk (Fig. 3 a-c), constrained by the postoutburst spectral energy distribution (SED) using a Monte Carlo radiative transfer code (Furlan et al. 2016), we computed the predicted N H toward the central protostar (Appendix E). Due to the high inclination angle in this model, the line-of-sight intercepts the upper layers of the accretion disk, which produce the bulk of N H ; moreover, this value is extremely sensitive to the viewing angle (Fig. 3d). The observed N H,X favors a slightly lower inclination. However, the available grid of protostellar models (Furlan et al. 2016) does not include a lower-mass accretion disk or a model with no accretion disk at all, which may reproduce both the SED and N H,X . These model trials illustrate the importance of X-ray constraints on N H to the overall consistency of protostar modeling.
In the optically thin case at 6.4 keV, that is to say, as long as the H column density by a cold absorber with a solar elemental abundance (N H ) is lower than ∼10 24 cm −2 , the iron line EW produced by photoionization is EW = 2.5 (Ω/4π)(N H /10 22 cm −2 ) eV, where Ω is the subtended angle of the irradiated material from the X-ray irradiating source (Tsujimoto et al. 2005). Since Ω/4π < 1, the large observed value of EW in HOPS 383 rules out the optically thin condition. In the optically thick case of the photoionization of photospheric iron by an X-ray flare, the iron line EW that is computed by Monte-Carlo radiative transfer simulations is limited to ∼130 eV for solar iron abundance and can only be increased by factors < 2 by disk flaring (Drake et al. 2008). However, large iron line EW can be obtained in both optical depth conditions when the photoionizing X-ray source is partly eclipsed (Drake et al. 2008). In the model, the X-ray source is located on the obscured region of the star and near the base of the accretion funnels (Fig. 3a). This source irradiates the accretion funnels, the upper layers of the accretion disk (Fig. 3b), and the circumstellar envelope (Fig. 3c).

Conclusions
Our detection of an X-ray flare from HOPS 383 provides conclusive evidence that strong magnetic activity is present at this bona fide Class 0 protostar. The resulting X-ray irradiation contributes to the ionization of the base of the outflow (Shang et al. 2004), and the magnetic reconnection that triggers powerful flares, similar to the one we detected, most likely drives energetic particle ejections (Feigelson et al. 2002). Such protostellar cosmic rays have been proposed to collide with refractory dust grains located at the inner edge of the accretion disk, inducing spallation reactions that could yield short-lived radionuclei, such as are observed in the refractory inclusions of chondritic meteorites (Gounelle et al. 2013;Sossi et al. 2017). Mass ejections have similarly been invoked in the case of the bona fide Class 0 protostar OMC-2 FIR 4 to explain the production of free electrons via collisions in the envelope, thereby enhancing the abundances of molecular ions (Ceccarelli et al. 2014). The Class 0 stage (∼10 000 yr) has an ∼10 times shorter duration than the Class I stage; however, during the earliest stage, half of the envelope mass is accreted, building the central star and the accretion disk. Therefore, the determination of the internal irradiation level in Class 0 protostars is paramount for the understanding of protostellar chemistry.
Longer observations are required to determine the flaring activity level and to collect more photons from such deeplyembedded, nascent stars. Athena X-IFU (Barret et al. 2018), which is scheduled to be launched at the beginning of the 2030s, will improve the X-ray throughput and spectral resolution at 6.4 keV, enabling a potential leap in our knowledge of the onset of protostellar magnetic activity.
Acknowledgements. We thank the anonymous referee for the careful reading of the manuscript. We are grateful to Joshua Wing (Chandra X-ray Center; CXC) and Steve Heathcote (Director of the Cerro Tololo Inter-American Observatory) for the scheduling of the simultaneous Chandra and SOAR observations. K.H., D.P., and N.G. are thankful to Jonathan (Jay) Elias (SOAR Director) and Patricio Ugarte (Observer Support) for technical support during the Spartan observation. N.G. thanks Patrick Broos and Leisa Townsley for Acis Extract support. This work was supported by "the Programme National de Physique Stellaire" (PNPS) of CNRS/INSU co-funded by CEA and CNES. Support for this work was provided by the National Aeronautics and Space Administration (NASA) through Chandra Award Numbers GO7-18012A and B issued by the CXC, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the NASA under contract NAS8-03060. The scientific results reported in this article are based on observations made by the Chandra X-ray Observatory (ObsIDs 18927, 20882 and 20883). This research has made use of software provided by the CXC in the application packages CIAO. SOAR is a joint project of the Ministério da Ciência, Tecnologia, e Inovação (MCTI) da República Federativa do Brasil, the U.S. National Optical Astronomy Observatory (NOAO), the University of North Carolina at Chapel Hill (UNC), and Michigan State University (MSU). This publication makes use of data products from WISE, which is a joint project of the University of California (Los Angeles) and the Jet Propulsion Laboratory/California Institute of Technology (JPL/Caltech), and NEOWISE which is a project of the JPL/Caltech. WISE and NEOWISE are funded by the NASA. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/web/gaia, processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of Python and astropy, corner, matplotlib, numpy, pyregion, scipy and statsmodels modules.     (Saraceno et al. 1996;André et al. 2000;André et al. 2008). In this scenario, the Class 0 and I border zone is defined by M env = 0.1 L bol (constant accretion, see André & Montmerle 1994) and M env = 0.15 L bol 0.6 (exponentially decaying accretion, see Bontemps et al. 1996). Propagating the model and flux uncertainties of HOPS 383 (Safron et al. 2015), we gather that: M env = 0.70 ± 0.14 M , L post-outburst bol = 7.5 +6.5 −1.5 L , and L pre-outburst bol = 0.21 +0.19 −0.05 L (obtained from the preoutburst-tooutburst flux ratio at 24 microns of 35 ± 3). The post-outburst position of HOPS 383 is within this Class 0 and I border zone due to its high accretion luminosity, but the pre-outburst position of HOPS 383 is well above it, where M env is larger than M , confirming its bona fide Class 0 stage. For comparison purposes, the infrared source IRS 7B in the R Coronae Australis star-forming core, which has a deeply embedded (N H ∼ 3 × 10 23 cm −2 ) X-ray counterpart that was proposed as a Class 0 protostar candidate (Hamaguchi et al. 2005), is on the side of the Class I stage in this plot, based on the bolometric luminosity reported in Groppi et al. (2007). The bolometric luminosity of the Class 0 candidate CMMS4 (Kamezaki et al. 2014) is not accurate enough for a robust conclusion.
We computed the protostellar evolutionary tracks assuming the following (Bontemps et al. 1996;André et al. 2000;André et al. 2008): an initial envelope mass, M env (0), decays exponentially with a characteristic timescale, τ = 10 5 yr; a mass accretion rate,Ṁ acc (t) = M env (t)/τ, where the local star-formation efficiency is = 0.5 (André et al. 2008); and L bol (t) = L acc (t) + L (t), where L acc (t) = GṀ acc (t)M (t)/R (t) is the accretion luminosity, where G is the gravitational constant, and M (t), R (t), and L (t) are the protostellar mass, radius, and interior stellar luminosity, respectively, on the birthline of Fig. 1 in Palla & Stahler (1999) where a constant accretion rate of 10 −5 M yr −1 is assumed to compute the radius-mass relation (Stahler 1988;Palla & Stahler 1992). We estimate with this toy model that the final stellar mass of HOPS 383 will be ∼0.4M and that less than ∼10% of it has already been accreted onto the central protostar.  Notes. The profile-fit photometry was computed from the singleexposure (L1b) source tables using the mean and standard deviation of the N profile-fit fluxes in the W2 band (w2mpro) obtained during each observing epoch. We note that ∆MJD is the time interval between the first and the last frame for each epoch.    Table B.1).

Appendix B: WISE and NEOWISE
For each observing epoch, we selected the W2-band frames of a good quality (qual_ f rame > 0) and stacked them with the WISE/NEOWISE coadder 2 using simple area weighting and a pixel scale of 1 . 375 for the final image. HOPS 383 is still visible after 2017 February in these stacked images as a faint source. To obtain aperture photometry, we used a standard 8 . 25 radius aperture centered on the Spitzer position of HOPS 383 (Megeath et al. 2012) and a custom 11 radius backgroundaperture located 22 west of it, on a region free of point sources with lower extended emission compared to the 50 − 70 radius annulus used in the single exposure source tables. The aperture correction (w2mcor) is taken from the single exposure (L1b) source table. The uncertainty of the magnitude zero-point (MAGZP) is fixed to the value of the WISE all-sky single exposures (σ MAGZP = 0.037 mag). The aperture flux decays by 2.40±0.06 mag in 7.5 years (white dots in Fig. B.1 and Table B.2), which is consistent with the behavior of the profile-fit flux, despite contaminating sources in the source aperture. The aperture flux increases by 0.25±0.06 mag between 2017 September and 2018 February, before coming back in 2018 September to the lower level that was observed previously. We conclude that the decay of the 2008 accretion outburst ended by 2017 Septem-2018-06-11) 6 . The AE source apertures were centered on source positions with a default PSF fraction of 90% at 1.49 keV, which was determined from point source simulations computed by raytracing with MARX . The AE background regions were built to get enough background counts in order to have the source photometry errors within 3% of the values that would be obtained with no uncertainty in the background. New positions were determined for isolated sources with an off-axis position lower than 5 using the mean position of the extracted events. These positions were used to determine the astrometric position offsets between the ObsIDs and a reference catalog. We used the Gaia data release 2 (DR2) 7 positions and position errors (Gaia Collaboration et al. 2016, 2018) that we propagated with topcat 8 (Taylor 2019) from the J2015.5 reference epoch to the J2017.948 Chandra epoch (using the same zero radial velocity and radial velocity error as in the Celestia 2000, Hipparcos & Tycho catalog) and added 2MASS sources without the Gaia DR2 best neighbor (Marrese et al. 2019); we excluded a 2MASS source that was resolved by Gaia and Chandra. The new positions were also determined for isolated sources with an off-axis position larger than 5 using the PSF correlation position. For crowded sources (PSF fraction < 87%), the peak in a maximumlikelihood image reconstruction was used. The correlation and reconstruction were performed using a multi-ObsID event image and a composite PSF image. These steps were repeated until converging to no offsets between the ObsIDs.
(Eq. (6.9) in Feigelson & Jogesh Babu 2012, and cyan curve in Fig 2c); and the optimal bandwidth, h cv , maximizing the crossvalidation (Eq. (6.10) in Feigelson & Jogesh Babu 2012, and green curve in Fig 2). The latter is ∼3.6 larger than the former, producing a smoother light curve in particular at the end of the X-ray burst.
Then, we used an adaptive (bandwidth) kernel estimator (Eqs. (6.14)-(6.15) in Feigelson & Jogesh Babu 2012) with the recommend value of 0.5 for the sensitivity parameter α. This was computed from the pilot densities previously obtained with h rot (blue curve in Fig 2c) and h cv (red curve in Fig 2c).

Appendix C.5: Spectral analysis
The source-plus-background spectrum was limited to the flare time interval ranging from the first to the last detected event and the 0.5-9.9 keV energy range. An energy-dependent aperture correction was applied by AE to the ancillary response file, which calibrates the extraction. The background spectrum was extracted from the full exposure of the first observation in the 0.1-10 keV energy range.
The simultaneous spectral fitting of the background and source-plus-background spectra was made with XSPEC 9 (version 12.10.1) using the C-statistic applied to unbinned data. The background cumulative spectrum was modeled with the AE cplinear model, a continuous piecewise-linear function with ten vertices ). The source spectrum was modeled with an interstellar medium absorption (Wilms et al. 2000) and an isothermal collisional-radiative plasma (Smith et al. 2001) with the typical coronal abundances of pre-main sequence stars (Güdel et al. 2007) plus an emission line with zero width at 6.4 keV (TBabs*(vapec+gaussian)) to which the background model, re-scaled to the source aperture, was added to fit the source-plus-background cumulative spectrum (data flow diagram in Fig. 10 of Broos et al. 2010). The model parameters, namely the hydrogen column density, the plasma temperature and normalization, the line normalization, and the ten parameters of cplinear, were first estimated by C-statistic minimization (C-stat=677.0 using 1309 degrees of freedom; Fig. C.1a). The emission line was added to improve the modeling below ∼6.4 keV (C-stat=671.8 using 1308 degrees of freedom; Fig. C.1b).
We fixed the cplinear parameters since the background contribution to the source-plus-background spectrum (0.1 count) was negligible. We performed a Markov chain Monte Carlo (MCMC) (Hogg & Foreman-Mackey 2018) to compute the probability density functions of the four physical parameters of this model using the Jeremy Sanders' xspec_emcee program 10 that used a Python implementation of the MCMC ensemble sampler with affine invariance proposed by Goodman & Weare (2010), emcee 11 (version 2.2.1; Foreman-Mackey et al. 2013). We customized xspec_emcee to output subsidiary Markov chains for the acceptance fraction (Foreman-Mackey et al. 2013) and the following physical values: soft and hard fluxes; soft and hard absorption-corrected fluxes; and line equivalent width and counts. We used uniform priors for the physical parameters, N w = 46 walkers that were started clustered around the previous parameter estimates, and N iter = 10 6 iterations.
To test for convergence of the Markov chains, we followed Hogg & Foreman-Mackey (2018) by using the integrated au- tocorrelation time that we computed for each walker and then averaged these estimates 12 for each parameter, and retained the largest value, τ (top panel of Fig. C.2). The MCMC burn-in phase was supposed to last less than N burn-in = 10τ = 10 008 iterations (Foreman-Mackey et al. 2013), as this was illustrated by the mean of the walkers' acceptance fractions converging to ∼0.3, which is well within the recommended range (Foreman-Mackey et al. 2013) of 0.2-0.5 (bottom panel of Fig. C.2b).
This MCMC burn-in phase was cut-off before producing the corner plots (Fig. C.3a) from which the parameter median values and 90% confidence interval 13 were determined (Table C.2). For comparison purposes, Fig. C.4 shows the sourceplus-background extracted cumulative counts spectrum (gold), the model for the parameter median values (green), and a random subset of 200 MCMC models (gray).
In this MCMC, the number of effectively independent data points is N = (N iter − N burn-in ) × N w ≈ 4.6 × 10 7 and the relative 12 See Dan Foreman-Mackey's blog on autocorrelation time estimation (2017 October 16): https://dfm.io/posts/autocorr/ . 13 We followed Hogg & Foreman-Mackey (2018), who are in favor of reporting median and quantiles, which are based on integrals, and against mode, which is not. accuracy is 100 √ τ/N = 0.5%. The probability of the emission line was computed as the probability to have at least 0.5 count from this line (Fig. C.3b).

Appendix C.6: Quiescent luminosity
We used the CIAO tool aprates to compute the source count rate in the 2-8 keV energy band during our observations and excluding the flare time interval, which led to a live time of ∼72,000 s. The source-plus-background aperture centered on the source position had an area of ∼2.4 2 and corresponding to a PSF fraction of 90% at 1.49 keV. The background aperture was a ∼2 . 5-∼30 annulus centered on the source position where the neighbor source region was excluded with an area of ∼2700 2 , corresponding to a PSF fraction of 0.9% at 1.49 keV. The total counts in source-plus-background and background apertures were 1 (at 3.0 keV) and 304, respectively. Inside the source-plusbackground aperture, we estimated a quiescent source count rate of 0.011 count ks −1 (0-0.057 count ks −1 , 90% C.I.).
From a tbabs*vapec model with N H,X = 7.0 × 10 23 cm −2 and a typical plasma quiescent temperature (kT = 1 keV), we determined, using XSPEC, that the aperture corrected, observed flux and absorption-corrected luminosity were lower than 2.0 × 10 −15 erg cm −2 s −1 and 2.0 × 10 30 erg s −1 in the 2-8 keV energy band, respectively.  (Safron et al. 2015).   with SExtractor (Bertin & Arnouts 1996) and SWarp (Bertin et al. 2002). We wrote a custom IDL program to automatically detect and correct, in the dark and flat-field corrected frames, any residual bulk images that were produced by bright stars, close to saturation in the raw frames. In each raw frame, this program identified any saturated regions centered on bright stars and looked at these detector positions in the subsequent calibrated frames for any excess above the sky level by fitting a constant level or a parabola (based on an F-test). These excesses were then subtracted or masked when too close to a source. These calibrated and corrected images were registered using the 2MASS catalog and stacked using median values and a pixel scale of 0 . 5.
In the H 2 narrow-band filter image, the shocked molecular hydrogen emission has a higher S/N compared to the K-band image. We produce for validation a pure H 2 emission image by subtracting the K-band image scaled by a factor that we obtain by minimizing the residuals on stars in the final image (Navarete et al. 2015).

Appendix E: Density model
We started from the protostellar model of HOPS 383, which is composed of an accretion disk in an infalling envelope with bipolar cavities (Furlan et al. 2016). The fixed parameters in this model grid (Table 3 of Furlan et al. 2016) that was computed with the Monte Carlo radiative transfer code, ttsre 15 (2008 version; Whitney et al. 2003) were: the stellar mass (M = 0.5 M ), the stellar effective temperature (T = 4000 K), the disk mass (M disk = 0.05 M ), the magnetospheric truncation radius of the gas disk (R trunc = 3 R ), the dust-disk inner radius (set to the dust sublimation radius, R sub /R = (T sub /T ) −2.1 = 6.8 where T sub = 1600 K is the dust sublimation temperature; Whitney et al. 2003), the scale height of the disk at the dust sublimation radius (set to the hydrostatic equilibrium solution; Eq. (35) of D'Alessio et al. 1998), the mean molecular mass per hydrogen nucleus (µ = 2.3), the radial (A = 2.25) and vertical (B = 1.25) exponent in disk density law, the fractional area of the hot spots on the star ( f spot = 0.01), the envelope outer radius (R env = 10 000 au), the centrifugal radius of the envelope (R c = R disk ), the exponent of the cavity polynomial shape (b cav = 1.5), and the vertical offset of the cavity wall (z cav = 0 au).
The envelope density was given by Eqs. (1) and (2) of Whitney et al. (2003). We note that the bipolar cavities in Furlan et al. (2016) are free of dust (and gas); for consistency, we filled them in with gas, using a constant density (n H 2 = 2 × 10 4 cm −3 ; Whitney et al. 2003).
The disk density was given by Eq.
(3) of Whitney et al. (2003) and normalized to M disk from R trunc to R disk using a spherical outer boundary (Fig. 3b this work and Fig. 2b of Whitney et al. 2003). We added accretion funnels at R trunc , which were dust-free as located inside the dust sublimation radius, assuming for simplicity a dipole geometry for the magnetic field (Hartmann et al. 1994), and using f spot to compute a latitude band between the smallest and largest magnetic loops at the stellar surface. The accretion-funnel density was given by Eq. (9) of Hartmann et al. (1994).
The values of N H along the line-of-sight were computed by numerically integrating these density equations. In this model grid, the gas-to-dust ratio is R = 100 and the protostellar dust opacity in the visible is κ ext,V = 5.3 × 10 4 cm 2 g −1 (Fig. 6 of Ormel et al. 2011, (ic-sil,gra) at 0.3 Myr), leading to: N H /A V = ln 10 × (1 + R)/(2.5 µ m H κ ext,V ) = 4.6 × 10 20 cm −2 mag −1 where m H is the hydrogen mass. Therefore, the observed N H,X corresponds to A V ∼ 1500 mag.