Modelling the eruptive young stellar object Re 50 N IRS 1 with ProDiMo

Andreas Postel

Abstract
Context: Episodic accretion plays an important role during the early phases of star-formation. The main processes responsible
arXiv:2312.16209v1 [astro-ph.SR] 23 Dec 2023

for the episodic accretion events remain, however, unclear.

Aims: Our main objective is to investigate the properties of FUors and EXors by analysing observational data, along with
numerical hydrodynamics simulations of protostellar disks, stellar evolution models of outbursting stars and thermo-chemical
models of star–disk systems in the outburst state. Our goal is to get a better understanding of the outburst processes and their
respective origin.
Methods: We used the radiation thermo-chemical code ProDiMo (PROtoplanetary DIsk MOdel) to match the dust emission
and gas emission lines originating from the environment surrounding the FUor star Re 50 N IRS 1. Our model focusses on
the observational data obtained by Herschel and Spitzer while we use archival photometry to complete the spectral energy
distribution.
Results: The modelling shows that the object is composed of a complex combination of several heating sources with different
spatial distribution. Our model uses a massive envelope with an mass infall rate of 1.35 ×10−5 M⊙ yr−1 to explain the continuum
emission in the (sub-)mm regime. At the same time we fit the CO and 13 CO ladders from 60 µm to 650 µm along with the two
[O i] lines centered at 63.18 and 145.53 µm. To explain the strong CO emission at shorter wavelengths and the oxygen lines, we
require a very warm disk due to a high disk accretion rate reaching 6 ×10−4 M⊙ yr−1 and an additional UV field of 3% of the
overall emission to heat the disk.
Keywords: stars: formation – stars: protostars – stars: pre-main sequence – protoplanetary disks – accretion, accretion disks –
infrared: stars

1 Introduction the transition between the two classes may be smooth (Audard
et al. 2014). The formation of massive stars also diplays signa-
Episodic accretion is recognized to play a central role in the tures of episodic accretion (e.g. De Wit et al. 2017; Meyer et al.
accretion history of young stars. One of the most outstanding 2017), suggesting that this phenomenon is somewhat universal
examples of episodic accretion are FUors, named after the among star formation.
archetype FU Orionis, which are low-mass young stars that A number of investigations at the infrared and the sub-
experience strong optical outbursts that can last several decades. millimeter regimes during the recent years finally allowed ac-
FU Orionis had its outburst in 1936, reaching an increase of 6 cessing of the dust emission and lines of FUors and EXors in
mag in the B-band at the peak, and has since then been slowly this wavelength range, which provided insights into the cir-
fading (Kenyon et al. 2000). A few dozens of objects with cumstellar disk and possible envelopes around these objects
similar behaviour have been identified since the outburst of (Lorenzetti 2005, Green et al. 2013b). It was found to be ex-
FU Ori (Hartmann & Kenyon 1996; Reipurth & Aspin 2010). tremely common that young stars show activity in the form
There have also been detections of some candidates (Quanz of significant IR variability. The inner regions of FUors and
et al. 2007; Semkov & Peneva 2010; Reipurth et al. 2012; EXors, which are of high interest to reveal the origin of an
Fischer et al. 2012; Fischer et al. 2013) that have not been outburst, are challenging for single-dish observations as these
observed in the pre-outburst state, which however show typical phenomena cannot be proven unambiguously to occur in the
characteristics of FUors. FUors appear to be predominantly inner disk alone. Recent studies with the Atacama Large Mil-
protostars which still accrete from their parental envelopes or limeter/submillimeter Array (ALMA) and Karl G. Jansky Very
are in the early T-Tauri star phase with only a remnant envelope. Large Array (JVLA) opened a new window for investigations
Characteristic of FUors are the very strong outbursts, with a of the outburst objects with very high spatial resolution and
luminosity change from 4 to 6 mag which occur on timescales allowed the observation of hot inner disks (Liu et al. 2017),
of days to weeks and last for several decades. Another type density waves in the disks (Pérez et al. 2016), outflows of the
of eruptive low-mass stars is characterised by weaker (3 to objects (Ruı́z-Rodrı́guez et al. 2017a,b), and shifts of the snow
5 mag) and shorter outbursts which however can occur more line during outburst (Cieza et al. 2016).
frequently (Herbig 1989). EX Lupi is the archetype of this The origin of accretion bursts is still unclear and could be
second class, which is accordingly called EXors. There are caused by viscous-thermal and/or magnetorotational instabil-
indications that these objects are more evolved than FUors, but ities that form in the inner disk, tidal effects caused by close

1
companions or close flybys of external stars, or accretion of
gaseous clumps in gravitationally unstable disks (see Audard
et al. 2014 for a review). We aim to investigate the properties
of FUors and EXors by analysing observational data (Postel
et al. 2019), along with numerical hydrodynamics simulations
of protostellar disks (Vorobyov & Basu 2015), stellar evolu-
tion models of outbursting stars (Elbakyan et al. 2019) and
thermo-chemical models of star–disk systems in the outburst
state (Rab et al. 2017). The goal of our project is to get a better
understanding of the outburst processes and their respective
origin.
In this paper, we focus on the FUor-like object Re 50 N
IRS 1 and present the results of a ProDiMo model (ProDiMo
is a radiation thermo-chemical code to model protoplanetary
disks) during a re-brightening phase. Interestingly, Chiang et al.
2015 proposed that the object is going through a dust clearing
phase in the optical and near-IR. Focussing on mid-/far-IR and
sub-mm data, we analyse what Herschel and Spitzer observa-
tions can tell us in the constraints of the model. ProDiMo has
been extended to handle envelope structures (Rab et al. 2017),
allowing it to not only model T-Tauri stars but also the earlier
evolutionary stages like FUors and EXors.
The main questions we aim to answer in this paper are how Figure 1: Image of Re 50 N IRS 1 in PACS. The colors refer to the three
the gas lines of the object can be modelled and where they different photometry channels of PACS, centered at 70 µm (blue), 100 µm
(green) and 160 µm. The target coordinates of the FUor are marked with a blue
come from, as well as whether the envelope plays a special
cross. The beam of the longest used wavelength is illustrated by the circle in
role in the scenario of Re 50 N IRS 1 and if we need episodic the lower left corner. We used a log scaling and removed values below zero.
accretion to explain the observational data. The spectra originate from the inner 39 arcsec.

solves the temperature, chemical abundances and radiation

1.1 Re 50 N IRS 1 field in a self-consistent way under constraints of axisymmetric
Re 50 N IRS 1 / HBC 494 is a known FUor-like object (Audard distribution of gas and dust in two dimensions. The protostar is
et al. 2014) in the Orion A cloud, which is 414±7 pc (Menten presented by a point source with corresponding luminosity in
et al. 2007) away from us, with a wide spread (Großschedl et al. the center of the structure, surrounded by a protoplanetary disk
2018) based on Galactic longitude. The object has a bolometric and spherical envelope with outflow cavities. The code is used
luminosity of 90 L⊙ , based on the SED of our recent work in to derive atomic and molecular abundances in different regions
Postel et al. (2019). It was observed in the past with several of the modelled object, and calculate the spectral lines (both in
instruments, covering the SED from the optical to the mm emission and absorption) and the spectral energy distribution
range with photometry and spectra, including Herschel (Sep (SED). This model was already used in Rab et al. 2017, where a
2013), Spitzer (7 Nov 2008) and ALMA (Aug 2015, see Ruı́z- detailed explanation, about how the disk and envelope structure
Rodrı́guez et al. 2017b) observations. Our Herschel data shows is built, is explained. In White et al. 2019 this model was used
many emission lines of Re 50 N IRS 1, mainly from CO from to fit APEX data with 3 lines but without Herschel data. More
low to high excitation temperatures, which makes the object a details about the disk accretion heating are shown in Appendix
promising candidate for detailed modelling. The object shows A.
strong silicate absorption in the bending and stretching bands
and CO2 absorption. Re 50 N IRS 1 is known for its very wide
outflow, detected by ALMA, and strong diffuse emission in two 3 Modelling
Herschel instruments (PACS/SPIRE: 60-671 µm) extending to
in the northeast and to the west, see Fig. 1. The procedure to find a good matching model can be separated
in three steps. The first phase was to check the literature for
already determined parameters. After that we matched the
2 Methods SED with ProDiMo. Once this was done, we matched the
gas emission lines. We investigated several configurations in
The radiation thermo-chemical code, ProDiMo (see Kamp et al. ProDiMo to probe a large parameter space. In Table 1, we
2017, Thi et al. 2011, Woitke et al. 2009, Woitke et al. 2016), provide the parameters of the model that we consider fits best

2
 Table 1: Main parameters of the best fitting ProDiMo model.
Quantity Symbol Value Reference
stellar mass M⋆ 0.5 M⊙ this work
central source eff. temp. T⋆ 6000 K this work
central source luminosity L⋆ 90 L⊙ this work, based on Postel et al. (2019)
amount of UV radiation fUV 3% this work
disk gas mass Mdisk 0.02 M⊙ Cieza et al. (2018)
disk accretion rate Ṁ 6 ×10−4 M⊙ yr−1 this work
disk inner radius Rin 1.0 au assumed
disk tapering-off radius Rtap 100 au Cieza et al. (2018)
inclination incl 15.3° this work
column dens. pow. ind. ϵ 1 assumed
reference scale height H(100au) 10 au assumed
flaring power index β 1.1 assumed
envelope mass Menv ≈ 2.57 M⊙ this work
mass infall rate Ṁif 1.35 ×10−5 M⊙ yr−1 fitted
outer radius Rout 15000 au calc., based on Postel et al. (2019)
cavity opening angle βcav 15° Ruı́z-Rodrı́guez et al. (2017b)
dust to gas mass ratio δ 0.01 Cieza et al. (2018)
min. dust particle radius amin 0.05 µm assumed
max. dust particle radius amax 3000 µm assumed
dust size dist. power ind. apow 3.5 assumed
dust composition Mg0.7 Fe0.3 SiO3 60% assumed
(volume fractions) amorph. carbon 20% assumed
vacuum 20% assumed
visual extinction AV 3.1 this work
distance d 414 pc Menten et al. (2007)
Notes: The model with no UV field has fUV=0 and the model with no disk accretion has Ṁ=0. The assumed parameters were fixed during the fitting process.

the Re 50 N IRS 1 gas lines and SED (mainly the Herschel one, 3.2 SED matching
see below).
The second step in our analysis was to match the broad-band
SED, while focussing on the Herschel and Spitzer data. With
3.1 Literature the SED matching, we obtained better initial conditions for
Cieza et al. (2018) derived the dust mass from one ALMA changing the model parameters and fit the gas emission lines.
continuum image. Assuming a gas to dust mass ration of 100, This provides constrains for the disk and envelope dust struc-
they estimated the disk mass, which are both used in this work ture. We varied several parameters while we discuss here the
beside the tapering parameter. We use a reference scale height ones with the biggest impact on the SED. Other parameters
of the disk of 10 au which is slightly above the margin of that have been explored but are not discussed in detail are the
Cieza et al. (2018). Ruı́z-Rodrı́guez et al. (2017b) delivered a luminosity of the central source, effective temperature of the
measurement with ALMA of the opening angle of the cavity. central source and the cavity opening angle, as their impact
Other parameters like the effective temperature of the central on the SED was very limited. We focussed on Herschel and
source, the mass infall rate and the disk accretion rate as well Spitzer photometry and spectra between 10-700 µm. The Her-
as the amount of the UV radiation were derived in our analysis schel data were observed in 2013 (Postel et al. 2019), before
which is described later in more detail. Based on photometry the object was found to go through a brightening phase (Chiang
of Herschel (PACS) and Hubble, we estimate the outer radius et al. 2015). The data in the optical and near-IR, as well as
of the system of Re 50 N IRS 1 at up to 15,000 au. Within this the major part of the photometry in the sub-mm and mm are
distance, there is no significant overlap of the emission of the older (1988-2001). We take here the approach that such data,
object and the nearby reflection nebula. We assume that the line while informative, do not necessarily correspond to the SED
flux measured by Herschel from this region is not contaminated state during the Herschel observation, and therefore we do not
by another source for the modelling. optimize the parameters to match them. The Herschel data is
therefore used as our benchmark.
The first parameter we started with is the mass infall rate

3
onto the envelope, as this has the biggest impact on the emission centrated to the surface of the disk and did not lead to more
in the far-IR/sub-mm. In the model of the rotating envelope, this emission of high J CO lines at some point but photodissoci-
parameter has no heating effect on the structure. The change of ated the CO. We therefore introduced disk accretion to have a
the SED is shown in Fig. 2. For the mass infall rate for our best disk-internal heating process in addition to the UV field. Both
model (blue curve, Ṁ ≈ 10−5 M⊙ yr−1 ) the far-IR/sub-mm data processes, which affect only the gas, are disabled in the default
of Herschel ar matched very well, which was the main goal of settings but appear to be a necessity for this particular object.
this step, but there is an extreme variation in the wavelength
range of the Spitzer spectrum and for shorter wavelengths when
Mi f is varied. The next step was to match the mid-IR, for which 4 Results
the inclination of the object had a major impact, which left the
far-IR/sub-mm range mostly unchanged. We show in Fig. 3 In Fig. 7, we show three different models. Our best match result
the corresponding SED for a variation of the inclination. While of the SED and lines and two other models models (Mdot=0,
the model predicts higher flux between the Herschel PACS fUV=0) to compare the effects of disk accretion and UV field
spectrum and the Spitzer spectrum, a solution for the short of the central source: The other two models have either dis-
wavelength domain was covered within the parameter space abled disk accretion but an UV field or disk accretion but no
at an inclination between 10° and 16.25°, which we refine to UV field. Neither the UV field nor the disk accretion heating
15.3°. As the object is embedded and stellar parameters are mechanism are sufficient to raise the temperature of the disk
mostly unknown, we also ran a grid from 0.1-0.8 M⊙ where high enough to reproduce the observational measurements. We
we find minor differences in the SED. The range is reasonable believe therefore, that the model with both disk mass accretion
for such an object, and we eventually use a stellar mass of 0.5 and UV field is a better model to match the gas emission lines
M⊙ , however there is potential for large difference to the actual and keep the fit of the SED intact.
protostar. We further assume that there is extinction along the
line of sight that is not part of the model. This extinction is set
by the AV parameter. We find that we need some additional
extinction to match the short-wavelength regime. The result 10 8
Mdot=10 6
Mdot=10 5
of the entire SED matching process can be seen in Fig. 4, to- Mdot=10 4
10 9
gether with observational data and the spectrum of the protostar
without absorption caused by the envelope and disk. 10 10
F [ergcm 2 s 1]

3.3 Gas emission lines 10 11

The next step of our analysis was to match the gas emission 10 12

lines detected with Herschel while keeping the SED untouched.

10 13
Fig. 5 shows the emission line fluxes from Herschel together
with the ProDiMo modelled fluxes based on the SED-matching 10 14
model, with a variation of the stellar luminosity. The emission 10 1 100 101 102 103
wavelength [ m]
lines at long wavelengths (⪆ 300 µm) are matched quite well
by the model, but the shorter wavelength lines, being orders
Figure 2: Change of the SED with different mass infall rates (10−6 − 10−4
of magnitudes too low, are not properly reproduced. Since the M⊙ yr−1 ) for the envelope. The goal of this grid was to match the far-IR/sub-
short wavelength CO lines, originating from high J transitions, mm data. The Herschel and Spitzer spectra which worked as a benchmark for
and the oxygen emission are caused by material which is much the SED matching are shown in colors (Spitzer: brown; Herschel PACS: green;
hotter than the material which emits the long wavelength lines, Herschel SPIRE SSW: purple; Herschel SPIRE SLW: red)
we increased the effective temperature of the central source
of the model to heat some material further. The result can be
seen in Fig. 6. The impact was negligible. Looking at the gas
temperature, we found that the CO and oxygen of the disk was
simply not getting warm enough, and even raising the effective 4.1 Spectral properties and surrounding emis-
temperature of the central source to 12000 K did not help. sion
We therefore investigated another parameter and introduced The matching of the SED (Fig. 4) with ProDiMo, mainly the
a UV field around the central source, assuming that accretion Herschel and Spitzer data, was precise enough to recreate the
shocks near the protostar would create a significant UV radi- results from the observations reasonably well. At long wave-
ation field, indicated by the high level of [O i] emission, that lengths there is a difference to some ground based observations.
would in turn heat the disk and envelope temperature suffi- Also, there is a gap to observational data in the range of 800
ciently high. The temperature increase was, however, con- nm to 3 µm. These observations were performed before the

4
 close to the cavity, and emissions at ≥325 µm have their origin
only in the envelope. We therefore deduce that Re 50 N IRS 1
10 8 incl=10°
incl=22.5° must consist of a relatively complex set of temperatures, caused
incl=35°
10 9 by multiple heating processes which are spread over a larger
spatial area. In the Fig. 8, the level of the continuum emission
F [ergcm 2 s 1]

10 10 at the respective wavelength (visible in the top left corner of the

panels with CO emission, in micrometer) is shown, in addition
10 11
to the area where most of the emission in vertical and radial
10 12 direction is created. An important result for observations at
these specific wavelengths. The CO density is strongly reduced
10 13 in the envelope in the line of sight with the disk and the central
source, up to a distance of ≈ 4300 au. This effect is caused
10 14 1 by the freeze-out of the CO at low temperature with a weaker
10 100 101 102 103
wavelength [ m] effect at low densities (i.e. larger distance). Closer to the
outflow cavity, photo-desorption occurs as the UV radiation
Figure 3: Change of the SED with inclination angles, varying from 10-35°. can penetrate into then envelope, which is however blocked by
The goal of this grid was to match the mid-IR data. The colour code is the
the disk. The behaviour of the CO abundance structure is as
same as in Fig. 2.
expected, see Rab et al. (2017) for more details.

10 8
4.2 Rotational diagrams
Fig. 10 shows the rotational diagram for CO emission for our
10 9 best model and the Herschel observation of the object. This
figure is an alternative way to show the line emission from
F [ergcm 2 s 1]

10 10 Fig. 7, but in a normalized way. The diagram can be used to

obtain the excitation temperatures by measuring the slope of the
10 11 curve, and column density of the gas by measuring the offset for
observational dat (see Goldsmith & Langer 1999, Green et al.
10 12 2013a, Dionatos et al. 2013). The figure shows that the slope,
i.e. the temperature distribution of the curve is matched by the
10 13 model. However, we can see that the column density of the
10 1 100 101 102 103 model is slightly higher than what we observed. The approach
wavelength [ m]
to split the emission lines into three domains with assumed LTE
Figure 4: SED of the model (blue) and the observational data of photometry from Postel et al. (2019) turns out to be a solid approximation,
and spectroscopy (grey). The spectrum of the central source is shown in black. albeit used due to the different Herschel instruments (PACS,
SPIRE SSW, SPIRE SLW).

re-brightening of the object. The difference to the ground-based

observations in the (sub-)mm range are explained by our model 4.3 UV field and accretion
by a change of the mass infall rate of the envelope by about one
order of magnitude (see Fig. 2). As we modelled the condition of Re 50 N IRS 1 during the
Fig. 8 shows that the number density of CO is particularly phase of higher activity, it is possible that the UV field and the
high in the envelope close to the outflow cavity (panels with accretion rate are lower now. As ProDiMo does not have the
spatial extent up to 3000 au) and in the disk (panelss with capabilities to model directly the effects of shocks leading to
spatial extent up to 50 au). A very high temperature of the sufficiently high temperatures in certain areas, especially along
gas causes a dissociation and depopulation of the CO (see Fig. the outflow cavity, we bypass this with the introduction of the
9) which consequently leads to a reduced emission from that UV field originating from the central source. For reaching
region. Too low temperature, however, does not result in the sufficiently high temperatures for the CO transitions with short
strong emission lines we see from the Herschel observation, as wavelengths and the [O i], the respective parameters appear to
the short lines would be then too faint. ProDiMo allows us to be inevitable in the constraints of ProDiMo. It turns out that the
pinpoint down the origin of the different CO lines, which is also UV field and the disk accretion are both a necessity to create
shown in Fig. 8. While the short-wavelength lines like the 72.8 the short emission lines in the observed combination, and the
µm emission still originates predominantly from the disk, the accretion of the envelope is the main parameter in our model,
emission at 186 µm already comes mainly from the envelope responsible for the overall shape of the SED for this object.

5
 line flux [Wm 2] line flux [Wm 2] line flux [Wm 2]

10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10

18
17
16
15
20
19
18
17
16
15
20
19
18
17
16
13CO 13CO 13CO 15
544.1 544.1 544.1
13CO 6 13CO 6 13CO 6
453.5 453.5 453.5
13CO 0 13CO 0 13CO 0
388.7 388.7 388.7
4 4 4

best
Obs.
Obs.
Obs.
13CO 13CO 13CO

fUV=0
340.1 340.1 340.1

Mdot=0
8 8 8

T=8000K
T=6000K
L=95 Lsun
L=90 Lsun
L=85 Lsun
13CO 13CO 13CO

(fUV = 3%) in combination.

T=10000K
L=100 Lsun
302.4 302.4 302.4
13CO 1 13CO 1 13CO 1
272.2 272.2 272.2
0 0 0
CO 6 CO 6 CO 6
50.25 50.25 50.25
CO 5 CO 5 CO 5
20.23 20.23 20.23
CO 4 CO 4 CO 4
33.56 33.56 33.56
CO 3 CO 3 CO 3
71.65 71.65 71.65
CO 3 CO 3 CO 3
25.23 25.23 25.23
CO 2 CO 2 CO 2
89.12 89.12 89.12
CO 2 CO 2 CO 2
60.24 60.24 60.24
CO 2 CO 2 CO 2
36.61 36.61 36.61
CO 2 CO 2 CO 2
16.93 16.93 16.93
CO 2 CO 2 CO 2
00.27 00.27 00.27
CO 1 CO 1 CO 1
86.00 86.00 86.00
CO 1 CO 1 CO 1
73.63 73.63 73.63
CO 1 CO 1 CO 1
62.81 62.81 62.81

6
CO 1 CO 1 CO 1
53.27 53.27 53.27
CO 1 CO 1 CO 1
44.78 44.78 44.78
CO 1 CO 1 CO 1
37.20 37.20 37.20
CO 1 CO 1 CO 1
30.37 30.37 30.37
CO 1 CO 1 CO 1
24.19 24.19 24.19
CO 1 CO 1 CO 1
18.58 18.58 18.58
CO 1 CO 1 CO 1
13.46 13.46 13.46
CO 1 CO 1 CO 1
08.76 08.76 08.76
CO 1 CO 1 CO 1
04.44 04.44 04.44
3 margin around the observational data. The change of the luminosity did not have a strong impact on the lines.

CO 9 CO 9 CO 9
3.35 3.35 3.35
CO 9 CO 9 CO 9
0.16 0.16 0.16
CO 8 CO 8 CO 8
7.19 7.19 7.19
CO 8 CO 8 CO 8
4.41 4.41 4.41
CO 8 CO 8 CO 8
1.81 1.81 1.81
CO 7 CO 7 CO 7
9.36 9.36 9.36
CO 7 CO 7 CO 7
7.06 7.06 7.06
CO 7
4.89 Figure 6: Same as Fig. 5, but for different central source effective temperatures, ranging from 6000 to 10000 K. CO 7
4.89
CO 7
4.89
CO 7 CO 7 CO 7
2.84 2.84 2.84
CO 7 CO 7 CO 7
0.91 0.91 0.91
CO 6 CO 6 CO 6
9.07 9.07 9.07
OI 1 4 OI 1 4 OI 1 4
5.53 5.53 5.53
OI 6 3 OI 6 3 OI 6 3
.18 .18 .18

disabled UV field (Ṁ = 4.2 × 10−4 M⊙ yr−1 , fUV = 0, orange diamonds). The best model uses both, disk accretion (Ṁ = 6 × 10−4 M⊙ yr−1 and the UV field
Figure 7: Line fluxes of the best model (green diamonds) compared to a model with disabled disk accretion (Ṁ = 0, fUV = 1%, blue diamonds) and with
Herschel line fluxes. The numbers next to the molecules/atoms refer to the wavelength of the respective lines in micrometer. The grey boxes represent a factor of
Figure 5: Emission lines after the SED matching was performed. Here we show the effect of different luminosities of the central source in comparison to the
 25 25
CO 72.8428

20 8.0 20 12.0

6.0 10.0

logn H [cm 3]
lognCO [cm 3]
15 15
z [au]

z [au]
4.0 8.0
10 10
2.0 6.0

5 0.0 5 4.0

0 0
10 20 30 40 50 10 20 30 40 50
r [au] r [au]

3000 3000
CO 650
CO 325.0
2500 CO 200.2725 3.0 2500
CO 186.0
9.0
2000 1.5 2000 8.0

logn H [cm 3]
lognCO [cm 3]
0.0
z [au]

1500 z [au] 1500 7.0

-1.5
1000 1000 6.0
-3.0
500 500 5.0
-4.5
0 0
0 500 1000 1500 2000 2500 3000 0 500 1000 1500 2000 2500 3000
r [au] r [au]

10000 10000
CO 650
CO 325.0
CO 200.2725 3.0
8000 8000 9.0
1.5
8.0

logn H [cm 3]
lognCO [cm 3]

6000 6000
0.0
z [au]

z [au]

7.0
4000 -1.5 4000
6.0
-3.0
2000 2000 5.0
-4.5
0 0
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
r [au] r [au]

Figure 8: Spatial distribution of the number density of CO and Hydrogen. The first line shows the smallest scale where the accretion disk is dominating. For
the CO emission, there is a gap between the disk and the envelope which extends to a radius of about 37 au at a height of 11 au. This gap is caused by numerical
effects at the transition region. The second and third line show plots up to 3000 and 10000 au distance from the central object, where the disk does not play a
significant role any more. The dashed lines show the level of the continuum emission at the respective wavelength. The dotted line shows the spatial region
where 15-85 % of the emission in vertical direction originates from. The box shows the spatial region where 15-85 % of the emission in radial direction are
originating from.

5 Discussion The matching of the emission lines shows that a simple enve-
lope + disk model of Re 50 N IRS 1 is not sufficient, but the
The major difficulty is the distribution of the temperature around system requires a higher complexity of additional heating mech-
the source, which cannot only be provided by the central source. anisms. We therefore use the UV field to model the effects of

7
 25 25
4.2
20 20 2.8
3.6

15 15 2.4
3.0

logTg [K]

logTd [K]
z [au]

z [au]
2.4 2.0
10 10

1.8 1.6
5 5
1.2 1.2
0 0
10 20 30 40 50 10 20 30 40 50
r [au] r [au]

10000 10000
4.2
8000 8000 2.8
3.6

6000 6000 2.4

3.0
logTg [K]

logTd [K]
z [au]

z [au]
2.4 2.0
4000 4000

1.8 1.6
2000 2000
1.2 1.2
0 0
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
r [au] r [au]

Figure 9: Spatial variation of the gas component (left) and dust (right) temperature. The highest temperatures of the gas are reached on the outflow cavity on
small distance to the central source. The hot area between the disk and envelope is caused by numerical effects due to the reduced number density of CO in this
region (see Fig. 8).

shocks in jets/outflows to recreate the high-J CO lines, since et al. (2013) which targeted the behaviour of CO.
ProDiMo does not have the capabilities to model them directly. For some of the spectral features in the SED, e.g. the narrow
Evidence for the presence of a protostellar jet is provided by absorption of CO2 ice at 15 µm, no matching happened since
the strong [O i] lines in the Herschel spectra and was shown in ProDiMo currently does not have the capabilities to model
other works to be a source for a radiation component (see Dion- these features.
atos & Güdel 2017 and Dionatos et al. 2020). The inclination The high accretion rate of the circumstellar disk of 6 ×10−4
and the opening angle of the outflow cavity of the object play M⊙ yr−1 would cause the accretion disk to dissipate, if it would
a main role in the visibility of the gas emission lines, as they stay at this high level for a long time, since the mass infall
affect the visibility of the accretion disk, where a large fraction rate of the envelope with 1.35 ×10−5 M⊙ yr−1 which feeds
of the emission originates from. The envelope of Re 50 N IRS the disk over a long time period is significantly lower. This
1 is likely still dense enough to absorb most of the emission means that the object cannot keep the current properties but
lines in the line of sight of the disk. Most of the CO emission the disk accretion is currently in a burst phase which wil last
originates in the inner disk, at some tens of au, where high J for a limited amount of time before a lowering of the disk
transition lines in our model originate from. For the intermedi- accretion rate must start. The Spitzer data have been observed
ate J transition lines, the envelope close to the outflow cavity some time before the Herschel data (2010 vs. 2013), so during
turns out to be the major source, where another high population the brightening of the object the SED in the near- and mid-IR
of CO is located. The cavity is also where the highest tempera- could have changed. A check with recent photometry or even a
tures of the gas occur, which is also where the radiation field spectrum would be interesting to confirm this. For a possible
of shocks would be created at. Both, the CO abundance and observation, a major part of the SED in the near- and mid-IR
temperature distribution is in agreement with former works on should be traced while the spectral resolution plays a minor
outbursts, like the hydrodynamical simulations from Vorobyov role.

8
 19 with the temporal dust clearing effect at shorter wavelengths.
Herschel This could also be a reason for the higher optical extinction and
18 ProDiMo
variation they mention, ranging from AV = 50 − 26 mag. Cieza
et al. (2018) derives an inclination of 70◦ , a result that is in
17
disagreement with our model. Based on the Spitzer data and as
Log10(Nu/gu)

16
we show in Fig. 3, the inclination in our model is not expected
to be above 22.5◦ while 15.3◦ reproduce the observational data
15 in the best way.
Assuming that the envelope temperature is in the order
14 of 10-20 K, most models of FUors infer a lower envelope
infall rate (a few 10−6 M⊙ /yr) if the core was initially close
13 to a Bonnor-Ebert sphere. The significantly higher infall rate
therefore means that the temperature of the core is much higher,
12
0 1.000 2.000 3.000 4.000 possibly due to a UV source, or the initial pre-stellar core was
Eu(K) off the equilibrium. The latter could happen if the core was
created by a shock wave.
Figure 10: CO rotational diagram for the best model compared to the Postel et al. (2019) include several objects with a similar
observational data obtained with Herschel. The gap around 2000 K is caused shape of the SED as Re 50 N IRS 1, e.g. PP 13 S, V346 Nor,
by a non-detection of the lines in the Herschel data. V883 Ori and HH 381 IRS. However, the emission lines among
these objects appear very different, raising the question of the
When calculating the stellar radius based on the values of reason for this. From our results with ProDiMo for Re 50 N
the modelled object via R star = 0.5 ṀGM star /Laccr (assuming IRS 1, it appears that these objects, if they show similar line
that half of the energy goes into the star and half escapes), properties, namely PP 13 S and V346 Nor, must host a similar
we would end up with a stellar radius of 52 R⊙ during the distribution of the temperature in the disk and outflow cavity.
outburst, a too large value compared to models. This equation Thus, objects like V883 Ori and HH 381 IRS, which show very
also assumes that the luminosity is dominated by the accretion faint emission lines, could consist of a colder disk, which would
luminosity and is therefore a lower limit. Given the small imply that the activity of these objects is way lower than the
impact on the SED, we think that a reduction of the stellar mass one of Re 50 N IRS 1. Also, inclination could play a role for
to 0.3 M⊙ is feasible, leading to a reduced radius of 31 R⊙ these objects, with an envelope shadowing emission lines from
which appears still high. A larger distance than the used 414 pc the disk, as proposed by Green et al. (2006). There could be
would lead to higher luminosities, and thus, a lower radius. A also a reduced outflow cavity, obscuring the disk. This would
10% increase (≈ 450 pc) would lead to 110 L⊙ , and a radius of mean that they are in an earlier evolutionary stage if the cavity
less than 26 R⊙ , which is more of a reasonable value according did not become large yet, but the envelope shows a similar mass
to recent numerical simulations of stellar characteristics during infall rate as later objects. Objects with a similar SED as Re
accretion bursts (Elbakyan et al. 2019). Großschedl et al. (2018) 50 N IRS 1 would be of interest for further modelling with
showed that the Orion A cloud is inclined from the plane of sky ProDiMo, to test these assumptions, and also multi-dish (sub-
and reaches up to 470 pc, depending on the Galactic longitude. )mm observations with high angular resolution would provide
Based on their Figure 3, at l = 211.6◦ , YSOs generally are valuable information about the inner regions of the bjects, in
have distances of 425-450 pc (the range is consistent with their particular the outflow cavity and the disk, which have not yet
derived median distance for the Tail of Orion A, L1641, ≈ 430 been observed, especially if they can provide information about
pc, although they derive an average distance of 409.36 ± 31.9 the temperature.
for the above longitude, using a 1-degree longitude bin, albeit While the observed features of Re 50 N IRS 1 could be well
with a large uncertainty), suggesting that a higher distance for reproduced in this work, some line features in other objects
Re 50 N IRS 1 is possible. could be an interesting challenge. HH 354 IRS shows a different
For some parameters, there seem to be strong differences shape of the SED, which indicates another evolutionary stage.
among other works. Gramajo et al. (2014) listed for the object a Informative here would be to check the connection of the SED
distance of 460 pc and a luminosity of 50 L⊙ while Chiang et al. with the presence of the strong oxygen emission at 63 µm and
(2015) came to a luminosity of 250 L⊙ . Also other parameters the relatively weak emission of CO and oxygen at 145.5 µm.
in their model, like the accretion rate, differ significantly from
our results. As we mentioned earlier, Chiang et al. (2015) 6 Summary
assume that a dust clearing effect occured, based on near-IR
data. However, we see from mid-IR to sub-mm data that there We provide a ProDiMo model of the young stellar object Re
must still be a very strong accretion, possibly in combination 50 N IRS 1 / HBC 494 which matches observational data of

9
Herschel and Spitzer during a phase of re-brightening of the Gramajo, L. V., Rodón, J. A., & Gómez, M. 2014, AJ, 147, 140
object, covering both the SED and emission lines of CO, 13 CO
and two [O I] lines. We find a way to recreate the gas lines Green, J. D., Evans, II, N. J., Jørgensen, J. K., et al. 2013a, ApJ,
of the target by an accreting envelope and a circumstellar disk 770, 123
with high accretion rate and an additional UV field to heat up Green, J. D., Evans, II, N. J., Kóspál, Á., et al. 2013b, ApJ, 772,
the disk and outflow cavity to form the necessary environmental 117
conditions. The gas lines appear to have their origin in the disk
and the transition region of the envelope and the outflow cavity, Green, J. D., Hartmann, L., Calvet, N., et al. 2006, ApJ, 648,
created by shocks in jets/outflows. The envelope turns out to be 1099
dense enough to absorb most emission in the line of sight to the
disk. The system requires episodic accretion to be consistent Großschedl, J. E., Alves, J., Meingast, S., et al. 2018, A&A,
over long time with a high disk accretion rate and a significantly 619, A106
smaller mass infall rate of the envelope. Hartmann, L. & Kenyon, S. J. 1996, ARA&A, 34, 207
Herbig, G. H. 1989, in European Southern Observatory Confer-
Acknowledgements ence and Workshop Proceedings, Vol. 33, European South-
ern Observatory Conference and Workshop Proceedings, ed.
We thank A. Kospal for providing the ALMA data
B. Reipurth, 233–246
for Re 50.
Kamp, I., Thi, W.-F., Woitke, P., et al. 2017, A&A, 607, A41
References Kenyon, S. J., Kolotilov, E. A., Ibragimov, M. A., & Mattei,
J. A. 2000, ApJ, 531, 1028
Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, Proto-
stars and Planets VI, 387–410 Liu, H. B., Vorobyov, E. I., Dong, R., et al. 2017, A&A, 602,
A19
Chiang, H.-F., Reipurth, B., Walawender, J., et al. 2015, ApJ,
805, 54 Lorenzetti, D. 2005, Space Sci. Rev., 119, 181

Cieza, L. A., Casassus, S., Tobin, J., et al. 2016, Nature, 535, Menten, K. M., Reid, M. J., Forbrich, J., & Brunthaler, A. 2007,
258 A&A, 474, 515

Cieza, L. A., Ruı́z-Rodrı́guez, D., Perez, S., et al. 2018, MN- Meyer, D. M.-A., Vorobyov, E. I., Kuiper, R., & Kley, W. 2017,
RAS, 474, 4347 arXiv e-prints [arXiv:1710.02320]

D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016,
500, 411 Science, 353, 1519

De Wit, W.-J., Caratti, A., & Kraus, S. 2017, in Star Formation Postel, A., Audard, M., Vorobyov, E., et al. 2019, A&A, 631,
from Cores to Clusters, 15 A30

Dionatos, O. & Güdel, M. 2017, A&A, 597, A64 Quanz, S. P., Henning, T., Bouwman, J., Linz, H., & Lahuis, F.
2007, ApJ, 658, 487
Dionatos, O., Jørgensen, J. K., Green, J. D., et al. 2013, A&A,
558, A88 Rab, C., Elbakyan, V., Vorobyov, E., et al. 2017, A&A, 604,
A15
Dionatos, O., Kristensen, L. E., Tafalla, M., Güdel, M., &
Persson, M. 2020, arXiv e-prints, arXiv:2006.01087 Reipurth, B. & Aspin, C. 2010, in Evolution of Cosmic Objects
through their Physical Activity, ed. H. A. Harutyunian, A. M.
Elbakyan, V. G., Vorobyov, E. I., Rab, C., et al. 2019, MNRAS, Mickaelian, & Y. Terzian, 19–38
484, 146
Reipurth, B., Aspin, C., & Herbig, G. H. 2012, ApJ, 748, L5
Fischer, W. J., Megeath, S. T., Stutz, A. M., et al. 2013, As-
tronomische Nachrichten, 334, 53 Ruı́z-Rodrı́guez, D., Cieza, L. A., Williams, J. P., et al. 2017a,
MNRAS, 468, 3266
Fischer, W. J., Megeath, S. T., Tobin, J. J., et al. 2012, ApJ, 756,
99 Ruı́z-Rodrı́guez, D., Cieza, L. A., Williams, J. P., et al. 2017b,
MNRAS, 466, 3519
Goldsmith, P. F. & Langer, W. D. 1999, ApJ, 517, 209

10
Semkov, E. & Peneva, S. 2010, The Astronomer’s Telegram,
2801
Thi, W.-F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711
Vorobyov, E. I., Baraffe, I., Harries, T., & Chabrier, G. 2013,
A&A, 557, A35
Vorobyov, E. I. & Basu, S. 2015, ApJ, 805, 115
White, J. A., Kóspál, Á., Rab, C., et al. 2019, ApJ, 877, 21
Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103

Woitke, P., Thi, W.-F., Kamp, I., & Hogerheijde, M. R. 2009,

A&A, 501, L5

A Disk accretion heating

For the viscous heating we use the approach of (D’Alessio et al.
1998) and adapted it to account for the tapered outer edge of
the disk.
 r   !(2−ϵ) 
3GM∗ Ṁ  R∗  r
Fvis (r) =
 
· 1 −  ·exp −  [erg cm−2 s−1 ],
8πr 3 r Rtap
(1)
where r is the distance to the star, G is the gravitational constant,
M∗ and M∗ the stellar mass and radius, Rtap the disk tapering-off
radius and ϵ the disk column density power index. To distribute
the heating per unit surface over the disk heat we use

ρP (r, z)
Γvis (r, z) = Fvis (r) R [erg cm−3 s−1 ], (2)
ρP (r, z′ )dz′

where ρ is the gas density and P = 1.25 to avoid unstoppable

heating in the upper disk layers where ρ drops towards zero.

11