# Anatomy of the massive star-forming region S106

## The [O I] 63 µm line observed with GREAT/SOFIA as a versatile diagnostic tool for the evolution of massive stars*

N. Schneider<fnref n="1" />,<fnref n="2" />, M. Röllig<fnref n="1" />, R. Simon<fnref n="1" />, H. Wiesemeyer<fnref n="3" />, A. Gusdorf<fnref n="4" />, J. Stutzki<fnref n="1" />, R. Güsten<fnref n="3" />, S. Bontemps<fnref n="2" />, F. Comerón<fnref n="5" />, T. Csengeri<fnref n="3" />, J. D. Adams<fnref n="6" />, and H. Richter<fnref n="7" />

Received 20 December 2017 / Accepted 20 May 2018

## ABSTRACT

The central area (40'' × 40'') of the bipolar nebula S106 was mapped in the [O I] line at 63.2 µm (4.74 THz) with high angular (6'') and spectral (0.24 MHz) resolution, using the GREAT heterodyne receiver on board SOFIA. The spatial and spectral emission distribution of [O I] is compared to emission in the CO 16→15, [C II] 158 µm, and CO 11→10 lines, mm-molecular lines, and con- tinuum. The [O I] emission is composed of several velocity components in the range from –30 to 25 km s⁻<fnref n="1" />. The high-velocity blue- and red-shifted emission (v = –30 to –9 km s⁻<fnref n="1" /> and 8 to 25 km s⁻<fnref n="1" />) can be explained as arising from accelerated photodissociated gas associated with a dark lane close to the massive binary system S106 IR, and from shocks caused by the stellar wind and/or a disk-envelope interaction. At velocities from –9 to –4 km s⁻<fnref n="1" /> and from 0.5 to 8 km s⁻<fnref n="1" /> line wings are observed in most of the lines that we attribute to cooling in photodissociation regions (PDRs) created by the ionizing radiation impinging on the cavity walls. The velocity range from –4 to 0.5 km s⁻<fnref n="1" /> is dominated by emission from the clumpy molecular cloud, and the [O I], [C II], and high-J CO lines are excited in PDRs on clump surfaces that are illuminated by the central stars. Modelling the line emission in the different velocity ranges with the KOSMA-τ code constrains a radiation field χ of a few times 10⁴ and densities n of a few times 10⁴ cm⁻<fnref n="3" />. Considering self-absorption of the [O I] line results in higher densities (up to 10⁶ cm⁻<fnref n="3" />) only for the gas component seen at high blue- and red velocities. We thus confirm the scenario found in other studies that the emission of these lines can be explained by a two- phase PDR, but attribute the high-density gas to the high-velocity component only. The dark lane has a mass of ~275 M☉ and shows a velocity difference of ~1.4 km s⁻<fnref n="1" /> along its projected length of ~1 pc, determined from H¹³CO⁺ 1→0 mapping. Its nature depends on the geometry and can be interpreted as a massive accretion flow (infall rate of ~2.5 × 10⁻<fnref n="4" /> M☉ yr⁻<fnref n="1" />), or the remains of it, linked to S106 IR/FIR. The most likely explanation is that the binary system is at a stage of its evolution where gas accretion is counteracted by the stellar winds and radiation, leading to the very complex observed spatial and kinematic emission distribution of the various tracers.

**Key words.** accretion, accretion disks – ISM: atoms – ISM: molecules – ISM: clouds – HII regions – photon-dominated region (PDR)

## 1. Introduction

Massive stars form by accretion of mass in one of three ways: the core accretion model similar to low-mass stars via a prominent disk (McKee & Tan 2002), the competitive accretion scenario in a clustered environment with a very small accretion disk (Bonnell 2007; Bate 2012), or the fragmentation-induced starva- tion view (Peters et al. 2011) by fragmentation from gravitational instability in dense accretion flows. Only the second and third scenarios straightforwardly explain why massive stars often form in multiple systems and are found preferentially at the centre of stellar clusters (Bontemps et al. 2010). Regardless of the scenario, what is known so far about how the different stages can be observationally traced, in particular with regard to the far-infrared (FIR), can be summarized as follows (see also Tan et al. 2014):

**Pre-stellar phase:** massive dense clumps (masses from a few tens to a few thousands of M☉, size ~1 pc, density n ≥ 10⁶ cm⁻<fnref n="3" />, temperature T ~ 15 K) are the locations where high-mass stars form (e.g. Beuther 2002; Motte et al. 2007; Tan et al. 2013). They often show subfragmentation into several smaller cores, i.e. pre- and protostellar objects, of size scale ≪0.1 pc (e.g. Bontemps et al. 2010), but there are also examples for isolated massive cores (e.g. Duarte-Cabral et al. 2013; Csengeri et al. 2017). Observationally, these cold dense clumps and cores are best traced by dust emission and molecular lines (such as N₂H⁺) in the mm and sub-mm wavelength range. In contrast, molecular line emission in the FIR is not observed in the cold dense gas phase because the energy levels of rotational transitions are too high. Only light hydrides such as CH⁺, CH, NH, OH, etc., can be seen in more diffuse gas. Molecular lines with lower excitation than the background continuum temperature can be observed in absorption. A good example is NH₃ with its transitions from

<footnote n="1">¹ I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, 50937 Köln, Germany
e-mail: nschneid@phl.uni-koeln.de</footnote>

<footnote n="2">² OASU/LAB, Université de Bordeaux, 33615 Pessac, France</footnote>

<footnote n="3">³ Max-Planck Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote n="4">⁴ LERMA, Obs. de Paris, ENS, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ.Paris 06, Paris, France</footnote>

<footnote n="5">⁵ European Southern Observatory, Alonso de Córdova 3107, Vitacura, Santiago, Chile</footnote>

<footnote n="6">⁶ SOFIA, USRA, NASA/Armstrong Flight Research Center, 2825 East Avenue P, Palmdale, CA 93550, USA</footnote>

<footnote n="7">⁷ DLR, Rutherfordstraße 2, 12489 Berlin-Adlershof, Germany</footnote>

<footnote>* Movie associated to Fig. 16 is available at http://www.aanda.org</footnote>
non-metastable to metastable levels, which have low excitation temperatures. For example, Wyrowski et al. (2012) successfully observed the 1.8 THz $NH_3$ line in absorption in high-mass star-forming clumps.

**Protostellar/stellar phase:** massive protostars contract so quickly that there is virtually no pre-main sequence stage. Hydrogen burning starts almost instantaneously, and they are main sequence objects while still accreting. However, no well-defined Keplerian disk has been detected yet, only disk-like structures on size scales smaller than 2000 AU (e.g. Beltran et al. 2011; Sanchez-Monge et al. 2013). The envelope-disk interaction leads to an atomic jet of at least partly ionized gas orthogonal to the disk, which can then drive a molecular outflow (e.g. Kuiper et al. 2011). For low-mass protostars and massive protostars in the core accretion model, the bipolar radio jet and molecular outflow are well collimated with a time-dependent activity and opening angle (Bontemps et al. 1996). For massive protostars in the competitive accretion scenario, the jet/outflow is often poorly collimated. Surrounding lower-mass protostars can cause multiple overlapping, randomly aligned outflows so that the overall outflow pattern can be rather complex (Hunter et al. 2008). A gas distribution with a highly fragmented and filamentary structure is predicted in the fragmentation-induced starvation scenario where the dense accretion flows should be directly observable. Radiative heating leads to a higher Jeans mass so that fewer but more massive stars form (Peters et al. 2010a).

Generally, outflows from massive protostars are more massive and energetic (Beuther 2002; Duarte-Cabral et al. 2013; Maud et al. 2015) than those from low-mass protostars. Extreme UV photons (EUV, energy >13.6 eV) lead to ionization of the protostellar outflows and create an outflow-confined H II region. One class of such H II regions are those with a morphology of an hourglass shaped parsec-scale bipolar nebula such as S106. Another example is associated with the evolved star MWC349A (Gvaramadze & Menten 2012) which in addition shows a subparsec (sub-pc) nebula that is explained as being due to a biconical outflow of ionized gas pinched at the waist by a small disk seen edge-on (Cohen et al. 1985). Typical of many bipolar nebulae is also a belt of dense cold gas in the equatorial plane of the nebula (Gvaramadze et al. 2010). Though it is commonly accepted that a disk oriented perpendicular to the symmetry axis of the H II region plays a crucial role in shaping the small-scale nebula/H II region and the large-scale bipolar nebula, it is not clear how this process works in detail. In addition, the nature of the dense gas structure in the nebula waist, where the existing star or stars are embedded, is not clear. It might be the remains of an accretion flow, as was shown in numerical simulations (Peters et al. 2010a), or the result of a shock that leads to an expansion and ultimately the dispersal of the belt (Gvaramadze & Menten 2012).

The earliest stages of hyper-compact, ultra-compact, and compact H II regions are sometimes difficult to identify because the massive protostar is deeply embedded in dense gas, and even cm-emission can become optically thick. Ionized collimated jets have been seen in radio continuum for some massive protostars (e.g. Gibb et al. 2003; Guzman et al. 2014). At the stage when a larger cavity is blown out, entrained gas can also arise from the ablation of gas at the cavity borders to the surrounding molecular cloud, caused by the stellar wind. Generally, as soon as the massive star reaches the main sequence, all forms of radiative feedback (thermal heating, ionization, radiation pressure on dust) and mechanical feedback (stellar winds from the stellar surface, protostellar winds/outflows from the magneto-centrifugally driven flows powered by accretion)

become important. Observationally, it is challenging to disentangle the different processes and to establish an evolutionary sequence for the protostellar phase because it is short (a few $10^3$ yr, Duarte-Cabral et al. 2013) and because of the small number of objects to study. An additional complication is that most massive stars form as binary or multiple systems, so that the observed jet and outflow patterns are even more complex and depend on the properties of each system.

It is the objective of this study and follow-up studies on the star-forming region S106 to identify and potentially better characterize these detailed phases by combining the spatial, velocity, and intensity information of atomic and molecular line observations and continuum from the mm to the FIR. In this first paper we focus on interpreting the [O I] 63 µm emission observed in the immediate environment of the exciting source S106 IR (see below) in the framework of a photodissociation region (PDR). A possible shock origin of the [O I] emission will be discussed and modelled in a subsequent paper, followed by a study of the large-scale structure of the bipolar nebula using [C II] 158 µm SOFIA data.

The bipolar H II region S106 in Cygnus X is an enigmatic object that has received considerable attention owing to its eye-catching appearance. Its distance was estimated to be 1.7 kpc (Schneider et al. 2007) and then determined with maser parallax measurements to be 1.3 ± 0.1 kpc (Xu et al. 2013). In the recent Gaia DR2 catalogue, the distance is derived from the trigonometric parallax to be 1.671 kpc (+738, -392 pc). In this paper, we use a value of 1.3 kpc. S106 IR was always thought to be a single late O to early B star (see Hodapp & Schneider 2008 for a summary). However, recent observations (Comerón et al. 2018) show that it is a close (separation < 0.2 AU), massive binary system, most likely consisting of a late O and a late B star, being responsible for the bipolar emission nebula (Fig. 1). Hereafter, we refer to this binary system as S106 IR. A small (only slightly larger than the beam of ~30 mas), edge-on disk-like feature around the two stars was discovered by cm-interferometry (Hoare & Muxlow 1996; Gibb & Hoare 2007), but its exact evolutionary status is not clear (Adams et al. 2015). The system already shows signatures of main sequence stars, emitting copiously in the UV (where the O star dominates) and driving an ionized wind with a velocity of 100–200 km s⁻¹ (Simon & Fischer 1982; Bally et al. 1983) into the bipolar cavity. Associated with this system is an IR cluster with more than 160 members (Hodapp & Rayner 1991). The centroid of the cluster lies about 30" west and 15" north of S106 IR and is symmetrically distributed about this location.

Figure 2 shows a cartoon of S106 with its known features and those observed to date. The two lobes of the H II region are inclined with respect to a vertical north-south axis. In addition, the northern (southern) lobe is tilted away from (towards) the observer (Solf & Carsenty 1982). A large portion of the front part of the southern cavity facing the observer has been eroded so that the backside of the cavity is visible in the optical and IR. The northern lobe is obscured by foreground gas and dust. The strongest nebular emission in the optical and near-IR arises from the shell-like structures that represent the ionization front and from the limb-brightened edges of the cavity walls. High-velocity emission of the [C II] 158 µm fine structure line and high-J CO lines (Simon et al. 2012) was interpreted as arising from swept-up material at the front and back sides of the expanding wind-driven cavity.

Another prominent feature related to S106 is a dark lane dividing the two lobes (indicated in dark blue in the cartoon) and lying apparently in front of S106 IR. It was initially
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Near-IR image (filters at J (1.25 µm), H (1.65 µm), and K' (2.15 µm)) of S106 taken with Subaru (Oasa et al. 2006), outlining the bipolar emission nebula, with contours of velocity integrated (–30 to 25 km s⁻<fnref n="1" />) [O I] emission (136 to 456 K km s⁻<fnref n="1" /> in steps of 64 K km s⁻<fnref n="1" />). The area mapped with SOFIA in [O I] and CO 16→15 emission is indicated with a white polygon. The star indicates the position of the S106 IR binary system and the triangle indicates the position of the young stellar object (YSO) S106 FIR.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Schematic view of the S106 region as seen on the sky. The nebula is tilted ≈25° to the east in the plane of the sky; the northern lobe, which is more obscured by foreground gas and dust, is inclined (≈15°) away from the observer. The southern lobe is inclined towards the observer and lacks the front side. The bipolar nebula and the H II region are embedded in an extended molecular cloud.</figcaption></figure>

interpreted as the shadow of a large-scale disk (Bally et al. 1983; Bieging 1984), but dust continuum observations (Vallée & Fiege 2005; Motte et al. 2007; Simon et al. 2012; Adams et al. 2015) revealed it to be an elongated high column-density feature associated with the extended molecular cloud (Schneider et al. 2007). Dust observations also show strong emission ≈15'' west of S106 IR, coinciding with two clusters of H₂O masers (Stutzki et al. 1982; Furuya et al. 1999). The dust peak was named S106 FIR by Richer et al. (1993), and interpreted as a Class 0 young stellar object (YSO). The molecular gas close to S106 IR is highly clumped (Schneider et al. 2002) and a number of authors interpreted their observations of ionic and fine-structure lines (van den Ancker et al. 2000; Schneider et al. 2003; Simon et al. 2012; Stock et al. 2015) as arising from photodissociation regions (PDRs) on the surfaces of these individual clumps around S106 IR. However, most of the studies of the FIR fine-structure cooling lines of S106 were hampered because only line integrated fluxes were observed and used for PDR modelling. Here we take a different approach in that we benefit from the unprecedented velocity resolution of the German REceiver

for Astronomy at Terahertz frequencies (GREAT) on board the Stratospheric Observatory for Far-Infrared Astronomy (SOFIA), and apply PDR modelling to line ratios in different velocity ranges so that we can determine more precisely the density and radiation field in individual portions of the gas.

In addition to getting a better understanding of the nature of the S106 star-forming complex itself, one of the objectives of this paper is putting our results into a larger context of how massive stars form and evolve. The paper is structured as follows. Section 2 gives observational details and Sect. 3 outlines how we realize the PDR modelling. In Sect. 4, we present maps of the observed lines and give results of PDR modelling. In Sect. 5, we develop a physical model of the S106 region and put the findings into a larger context of massive star formation. Section 6 summarizes the results of the paper.

## 2. Observations

### 2.1. SOFIA

The [O I] ³P₁ → ³P₂ atomic fine structure line at 4.74477749 THz (63.2 µm) and the CO 16→15 rotational line at 1.841345 THz (162.8 µm) were observed with the heterodyne receiver GREAT<fnref n="1" /> on board SOFIA during one flight on December 15, 2015, from Palmdale/California. The [O I] line was observed in the single-pixel H channel (now upgraded to a 7-pixel array), and the CO line in the single-pixel L2 channel. The average system temperatures over the total bandpass were 3641 and 5059 K for [O I] and CO, respectively. We employed a fast Fourier transform spectrometer (AFFTS). The centre intermediate frequency (IF) was 1455 MHz for the [O I] line and 1000 MHz for the CO line. The frequency resolution of the original [O I] and CO data

<footnote n="1">¹ German Receiver for Astronomy at Terahertz. GREAT is a development by the MPI für Radioastronomie and the KOSMA/Universität zu Köln, in cooperation with the MPI für Sonnensystemforschung and the DLR Institut für Planetenforschung.</footnote>
sets is 0.244 MHz, giving a velocity resolution of ~0.04 km s⁻¹ (16384 spectrometer channels).

Beam-switched on-the-fly maps with a scan length of 90", a step size of 3", and a dump time of 1 s were performed. The chopper throw was 360" with a chop angle of 45° anti-clockwise from +RA. The map centre position refers to S106 IR (RA, Dec)(J2000) = (20h27m26s.74, 37°22'47".9). The same area was covered three times while scanning in RA, and two times in Dec. Procedures to determine the instrument alignment and telescope efficiencies, antenna temperature and atmospheric calibration, as well as the spectrometers used, are described in Heyminck et al. (2012) and Guan et al. (2012). Tracking was done on the H-channel and the co-alignment between the H- and L2-channel is determined to be ~2". All line intensities are reported as main beam temperatures scaled with main-beam efficiencies of 0.69 and 0.68 for [O I] and CO, respectively, and a forward efficiency of 0.97. The main beam sizes are 15.3" for the L2 channel and 6.1" for the H channel.

The telluric [O I] line, originating from the mesosphere, contributes as a narrow feature at the velocity of the bulk emission of the cloud (approximately +3 km s⁻¹). We corrected for the absorption following the procedure described in Appendix A in Leurini et al. (2015), assuming that the profile can be characterized by a Gaussian. In summary, for the opacity correction the absorption strength was adjusted so as to achieve an adequate interpolation between adjacent unaffected spectral channels.

The calibrated [O I] and CO spectra were further reduced and analysed with the GILDAS software, developed and maintained by IRAM. From the spectra, a third-order baseline was removed and spectra were then averaged with 1/σ<fnref n="2" /> weighting (baseline noise). The mean r.m.s. noise temperatures per 0.6 km s⁻¹ velocity bin are 1.9 K for [O I] (6" beam) and 1.6 K for CO (15" beam). We estimate that the absolute calibration uncertainty is ~20%.

For some overlays and PDR modelling, we use [C II] 158 µm and $^{12}$CO 11→10 data observed with SOFIA from 2012 (Simon et al. 2012) and new observations from May 2015 to December 2016 (Simon et al., in prep.).

## 2.2. Complementary data

### 2.2.1. IRAM 30 m

We use molecular line data from the IRAM<fnref n="2" /> 30 m telescope, obtained in January 2016 (PI: N. Schneider). The observations were performed with two settings, one with the high-resolution spectrometer FTS50 attached to the EMIR E0 receiver (at 3 mm wavelength), and one with the lower resolution FTS200 in combination with EMIR E1 (at 1 mm wavelength). A large number of molecular lines (~20) was observed, covering a range from 85.3 to 106.1 GHz and from 89.8 to 110.5 GHz, respectively. We performed on-the-fly maps that consist of one horizontally and one vertically scanned map of size 300" at 3 mm and 168" at 1 mm (larger than is shown in Fig. 1). The centre position of all maps is the same as for the [O I] map. As emission-free reference position, we used a position with offset −700", 0" to the S106 map centre position. The quality of the data is very good, thus only baselines of first order were removed. For the H$^{13}$CO+ 1→0 data we employ here, the average system temperature during the observations was 88 K and the average baseline

rms for spectra smoothed to 15" angular resolution and with a velocity resolution of 0.17 km s⁻¹ is 0.16 K.

### 2.2.2. Other instruments

For overlays and PDR modelling, we use published or archived continuum data from Herschel (70, 160, 250, 350, 500 µm) and FORCAST (19.7, 25.3, 31.5, 37.1 µm) (Adams et al. 2015), Spitzer (3.6, 4.5, 5.8, 8.0 µm), MAMBO (1.2 mm) (Motte et al. 2007), the VLA (Bally et al. 1983), and Subaru (Oasa et al. 2006). The flux uncertainties for the continuum data are ~10% for PACS, SPIRE, and FORCAST, and ~20% for MAMBO (see references above for these values).

## 3. Photodissociation region modelling

We employ the KOSMA-τ PDR model (Störzer et al. 1996; Röllig et al. 2006) to derive local physical conditions from the observed line intensities and continuum fluxes. In the following, we explain the model and our approach in more detail because this methodology will also be used for subsequent studies.

### 3.1. Model description

The KOSMA-τ PDR model numerically computes the energetic and chemical balance in a spherical cloud that is externally irradiated. We here use a non-clumpy model approach, i.e. a single, spherical PDR with a density gradient similar to a Bonnor–Ebert sphere, and a cloud mass of $M = 100 M_\odot$ (the results depend only weakly on the model clump mass and we ran models with a parameter range between $10^{-3}$ and $10^3 M_\odot$). The full numerical computation scheme of KOSMA-τ involves three steps. First, the continuum radiative transfer code MCDRT (Szczerba et al. 1997) is used to compute the thermal balance of all dust components (see below) as well as the far-ultraviolet (FUV) radiative transfer within the model cloud, and the emergent continuum radiation (Röllig et al. 2013). We include a variety of different dust models, for example, the MRN model (Mathis et al. 1977) and the dust models by Weingartner & Draine (2001a). By assuming that the dust temperature is independent of the gas temperature, we then use the MCDRT output as input for the second step where the KOSMA-τ code computes the chemical and physical state of the gas. In Fig. 3, we show as an example a comparison between the observed mid-IR and FIR data at the position 3b (offsets 5", −25") and the best fitting dust continuum model. We determined in this way the dust continuum at all positions (nine in total) where we performed PDR modelling.

We assume a local chemical steady state, i.e. a balance between chemical formation and destruction processes of all involved species and compute their local densities and thermal equilibrium, i.e. a balance between all heating and cooling processes, resulting in local gas and dust temperatures. This is done for all radial grid positions in an iterative way to reach the numerical convergence criteria which is that the column density deviations between iterations are less than 1%. The result is the radial density and temperature profile of the model clump. In a last step the profile is used to perform the radiative transfer computations, giving the spectral line emission of the model cloud for comparison with observations (Gierens et al. 1992)<fnref n="3" />.

We use the most recent CO self-shielding functions by Visser et al. (2009) and assume a Doppler line width according to the

<footnote n="2">² IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain).</footnote>

<footnote n="3">³ The KOSMA-τ code was part of the PDR comparison benchmark study (Röllig et al. 2007).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** One example (position 3a, southern cavity) for KOSMA-τ PDR model continuum emission fitting the observed continuum observations (black points) of MAMBO (1.2 mm), Herschel (PACS: 70, 160, SPIRE: 250 µm), and FORCAST (19.7, 25.3, 31.5, 37.1 µm), assuming a unity beam filling factor. The error for the flux values is generally small (see Sect. 2.2). All data were smoothed to a common angular resolution of 15". The Herschel/SPIRE 250 µm data point has a resolution of 18".</figcaption></figure>

linewidth–size relation as given by Larson (1981). The photoelectric heating is computed according to Weingartner & Draine (2001b; model 4 from their Table 2). The formation of H₂ on grain surfaces follows Cazaux & Tielens (2002, 2004, 2010) and formation on polyaromatic hydrocarbon (PAH) surfaces is suppressed. For more details see Röllig et al. (2013). In addition, we use the UMIST Database for Astrochemistry chemical network (McElroy et al. 2013), including all isotopic reaction variants including ¹³C and/or ¹⁸O isotopes (Röllig & Ossenkopf 2013). In total 3766 reactions are considered and the 227 species that are included in the chemistry are listed in Table A.1. Numerically problematic reaction rate coefficients have been rescaled according to Röllig (2011). The formation of CH⁺ and SH⁺ is computed using state-to-state reaction rates (Agundez et al. 2010; Nagy et al. 2013).

## 3.2. Model fitting

Table 1 summarizes the model parameters. To derive local physical conditions, we compare observed line ratios<fnref n="4" /> with the corresponding model line ratios to find the best fitting model parameters. We use two line ratios, $[CII]_{158\mu m}/[OI]_{63\mu m}$ and CO(16–15)/CO(11–10), and compare them against $([CII]_{158\mu m} + [OI]_{63\mu m})/\Sigma_{FIR_{range}}$. We note that we do not use here the total FIR intensity $\Sigma_{FIR_{total}}$, which is commonly used in PDR models (Stock et al. 2015) because we work with velocity resolved data. The parameter $\Sigma_{FIR_{range}}$ is thus only the continuum intensity in one particular velocity range. We define five velocity intervals for the [C II] and [O I] emission as explained in Sect. 4. Observationally, however, we can only determine the total FIR intensity $\Sigma_{FIR_{total}}$ obtained over all velocity ranges. For that, we derive $\Sigma_{FIR_{total}}$ by numerically integrating the continuum data between 10 and 1000 µm, stemming from FORCAST, Herschel, and MAMBO (see Sect. 2.2) using a linear interpolation. Figure 3 shows an example of this procedure. For the determination of $\Sigma_{FIR_{range}}$ we then assume that gas and dust are well mixed and that the absolute [O I] 63 µm

<table><caption>**Table 1.** Overview of the most important model parameters.</caption><thead><tr><th colspan="2">Model input parameters</th></tr></thead><tbody><tr><td>He/H</td><td>0.0851</td><td>(1)</td></tr><tr><td>O/H</td><td>4.47 × 10<sup>-4</sup></td><td>(2)</td></tr><tr><td>C/H</td><td>2.34 × 10<sup>-4</sup></td><td>(2)</td></tr><tr><td><sup>13</sup>C/H</td><td>3.52 × 10<sup>6</sup></td><td>(3)<sup>a</sup></td></tr><tr><td><sup>18</sup>O/H</td><td>8.93 × 10<sup>-7</sup></td><td>(4)<sup>b</sup></td></tr><tr><td>N/H</td><td>8.32 × 10<sup>-5</sup></td><td>(2)</td></tr><tr><td>S/H</td><td>7.41 × 10<sup>-6</sup></td><td>(2)</td></tr><tr><td>F/H</td><td>6.68 × 10<sup>-9</sup></td><td>(2)</td></tr><tr><td>Z</td><td>1</td><td>solar metallicity</td></tr><tr><td>ζ<sub>CR</sub></td><td>2 × 10<sup>-16</sup> s<sup>-1</sup></td><td>CR ionization rate (5)</td></tr><tr><td>R<sub>V</sub></td><td>3.1</td><td>visual extinction/reddening (7, 8)</td></tr><tr><td>σ<sub>D</sub></td><td>1.75 × 10<sup>-21</sup> cm<sup>2</sup></td><td>UV dust cross section per H (8)</td></tr><tr><td>⟨A(λ)/A<sub>V</sub>⟩</td><td>3.339</td><td>mean FUV extinction</td></tr><tr><td>τ<sub>UV</sub></td><td>3.074 A<sub>V</sub></td><td>FUV dust attenuation</td></tr><tr><td>v<sub>b</sub></td><td>1 km s<sup>-1</sup></td><td>Doppler width</td></tr><tr><td>n<sub>0</sub></td><td>10<sup>2</sup>...7 cm<sup>-3</sup></td><td>total surface gas density</td></tr><tr><td>M</td><td>10<sup>2</sup> M<sub>☉</sub></td><td>cloud mass</td></tr><tr><td>χ</td><td>10<sup>0</sup>...6</td><td>FUV intensity w.r.t. (6)<sup>c</sup></td></tr><tr><td>σ</td><td>1.5</td><td>density power law index</td></tr><tr><td>R<sub>core</sub></td><td>0.2 R<sub>tot</sub></td><td>size of constant density core</td></tr><tr><td>N<sub>tot</sub>/A<sub>V</sub></td><td>1.62 × 10<sup>21</sup> cm<sup>-2</sup></td><td>(8)</td></tr></tbody></table>

**Notes.** All abundances are given with respect to the total H abundance. <sup>(a)</sup>Based on a ¹²C/¹³C ratio of 67, <sup>(b)</sup>based on a ¹⁶O/¹⁸O ratio of 500, <sup>(c)</sup>χ = 1.71 G₀ where G₀ is the mean ISRF from Draine (1978).

**References.** (1) Asplund et al. (2005), (2) Simon-Diaz & Stasinska (2011), (3) Langer & Penzias (1990), (4) Polehampton et al. (2005), (5) Hollenbach et al. (2012), (6) Draine (1978), (7) Röllig et al. (2013), (8) Weingartner & Draine (2001a).

line integrated intensity in all velocity<fnref n="5" /> ranges is higher than that of [C II] and CO (which is the case in each velocity range for each observed/analysed position, see Table 2). With these two assumptions, the fraction of [O I] line integrated intensity in one velocity range to that in the total velocity range is the same as the fraction of $\Sigma_{FIR_{range}}$ to $\Sigma_{FIR_{total}}$.

All maps are on the same angular resolution of 15" (also the mid- and FIR continuum data), except the CO 11→10 map, which has a resolution of ~20", and the Herschel/SPIRE 250 µm data, which have a resolution of 18". In the following, we thus focus on the values for density and UV field obtained from the $[CII]_{158\mu m}/[OI]_{63\mu m}$ ratio since beam filling factors are eliminated to the first order and because only the [O I] and [C II] lines show emission in all velocity ranges, in contrast to the CO lines that show no emission (or only marginal amounts) at the highest velocities. The model continuum intensities are provided by MCDRT; the model line intensities are computed by KOSMA-τ.

# 4. Results

## 4.1. [O I] Spectral line map and average spectrum

Figure 4 shows the [O I] spectral line map of the region outlined in Fig. 1 on a 6" grid, representing the original resolution of the

<footnote n="4">⁴ All ratios are computed from intensities in erg s⁻¹ cm⁻² sr⁻¹. See Table 2 for conversion factors.</footnote>

<footnote n="5">⁵ To convert the FORCAST data from MJy pixel⁻¹ to MJy sr⁻¹ we use the conversion factor 7.213 × 10¹⁰ pixel sr⁻¹.</footnote>
<table><caption>**Table 2.** Velocity integrated intensities in a 15'' beam (20'' for CO 11→10) at the locations indicated in Fig. 4.</caption><thead><tr><th>Pos.</th><th>Offset<br>('', '')</th><th>[O I]<sup>a</sup></th><th>[O I]<sup>a</sup></th><th>%</th><th>[C II]<sup>b</sup></th><th>[C II]<sup>b</sup></th><th>%</th><th><sup>12</sup>CO<sup>b</sup><br>16→15</th><th><sup>12</sup>CO<sup>b</sup><br>16→15<br>(total)</th><th>%</th><th><sup>12</sup>CO<sup>b</sup><br>11→10</th><th><sup>12</sup>CO<sup>b</sup><br>11→10<br>(total)</th><th>%</th><th>Total<br>FIR<sup>c</sup></th></tr><tr><td></td><td></td><td colspan="3">(total<sup>d</sup>)</td><td colspan="3">(total)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></thead><tbody><tr><td colspan="18" style="text-align: left; font-weight: bold;">High-velocity blue v = -30 to -9 km s<sup>-1</sup></td></tr><tr><td>1</td><td>(-2, -2)</td><td>13.4</td><td>28.9</td><td>46</td><td>45.0</td><td>308</td><td>15</td><td>3.6</td><td>45.7</td><td>8</td><td>4.7</td><td>36.3</td><td>13</td><td>24.1</td></tr><tr><td>2</td><td>(7, 12)</td><td>8.0</td><td>24.2</td><td>33</td><td>30.0</td><td>314</td><td>10</td><td>1.0</td><td>25.8</td><td>4</td><td>3.9</td><td>30.0</td><td>13</td><td>23.4</td></tr><tr><td colspan="18" style="text-align: left; font-weight: bold;">Outflow blue v = -9 to -4 km s<sup>-1</sup></td></tr><tr><td>3a</td><td>(2, -20)</td><td>4.3</td><td>21.5</td><td>20</td><td>82.8</td><td>292</td><td>28</td><td>4.6</td><td>25.4</td><td>18</td><td>6.2</td><td>24.8</td><td>25</td><td>15.3</td></tr><tr><td>3b</td><td>(-15, -5)</td><td>4.9</td><td>22.4</td><td>22</td><td>65.0</td><td>245</td><td>27</td><td>7.1</td><td>24.7</td><td>29</td><td>8.5</td><td>24.3</td><td>35</td><td>17.5</td></tr><tr><td>3c</td><td>(8, 8)</td><td>6.7</td><td>24.7</td><td>27</td><td>69.9</td><td>331</td><td>21</td><td>22.3</td><td>34.1</td><td>65</td><td>15.4</td><td>34.0</td><td>45</td><td>24.8</td></tr><tr><td>4</td><td>(0, 15)</td><td>7.2</td><td>23.2</td><td>31</td><td>50.5</td><td>273</td><td>18</td><td>13.6</td><td>16.1</td><td>84</td><td>12.0</td><td>25.5</td><td>47</td><td>19.7</td></tr><tr><td colspan="18" style="text-align: left; font-weight: bold;">Cloud bulk emission v = -4 to 0.5 km s<sup>-1</sup></td></tr><tr><td>3b</td><td>(-15, -5)</td><td>9.6</td><td>22.4</td><td>43</td><td>96.4</td><td>245</td><td>39</td><td>4.6</td><td>24.6</td><td>19</td><td>10.9</td><td>24.3</td><td>45</td><td>17.5</td></tr><tr><td>5a</td><td>(10, -5)</td><td>4.8</td><td>22.1</td><td>21</td><td>98.9</td><td>337</td><td>29</td><td>4.6</td><td>43.1</td><td>11</td><td>14.8</td><td>37.7</td><td>39</td><td>20.4</td></tr><tr><td>5b</td><td>(-10, 15)</td><td>7.4</td><td>19.0</td><td>39</td><td>57.4</td><td>205</td><td>28</td><td>-</td><td>-</td><td>-</td><td>8.1</td><td>20.3</td><td>40</td><td>12.4</td></tr><tr><td colspan="18" style="text-align: left; font-weight: bold;">Outflow red v = 0.5 to 8 km s<sup>-1</sup></td></tr><tr><td>3a</td><td>(2, -20)</td><td>4.7</td><td>21.5</td><td>22</td><td>38.2</td><td>292</td><td>13</td><td>6.9</td><td>25.7</td><td>27</td><td>4.0</td><td>24.5</td><td>16</td><td>15.3</td></tr><tr><td colspan="18" style="text-align: left; font-weight: bold;">High-velocity red v = 8 to 25 km s<sup>-1</sup></td></tr><tr><td>6</td><td>(10, 25)</td><td>5.4</td><td>14.7</td><td>37</td><td>16.8</td><td>170</td><td>10</td><td>2.1</td><td>19.5</td><td>11</td><td>0.6</td><td>1.50</td><td>4</td><td>9.8</td></tr><tr><td>3a</td><td>(2, -20)</td><td>3.2</td><td>21.5</td><td>15</td><td>46.1</td><td>292</td><td>16</td><td>0.5</td><td>25.7</td><td>2</td><td>0.3</td><td>24.5</td><td>1</td><td>15.3</td></tr></tbody></table>

**Notes.** Each position reflects a peak of emission in a certain velocity range and one or more tracers and/or a source position: Pos. 1 is the southern HV-blue [O I] peak, closest to the central binary system S106 IR; Pos. 2 is the northern HV-blue [O I] peak; Pos. 3a is the southern cavity blue outflow [O I] peak and the southern HV-red [O I] peak; Pos. 3b corresponds to the YSO S106 FIR, associated with H$_{2}$O maser emission; Pos. 3c is the CO 16→15 emission peak; Pos. 4 is the northern cavity blue outflow [O I] peak; Pos. 5a is the bulk CO 2→1 peak located in the dark lane; Pos. 5b is the NW bulk [O I] peak; Pos. 6 is the northern HV-red [O I] peak. The first column for each line gives the intensity in the respective velocity range, the second column the total line integrated (−30 to +25 km s⁻¹) intensity, and the third column the percentage of how much of the individual velocity integrated line contributes to the total line flux. The last column gives the total FIR flux, determined as explained in Sect. 3.2. The values for the red outflow velocity range (v = 0.5 to 8 km s⁻¹) should be treated with care because the [O I] and [C II] lines suffer from absorption and the CO lines are weak or show no emission. The overall uncertainty of the values given in this table is ~20%. (a)In units of (10⁻³ erg s⁻¹ sr⁻¹ cm⁻²). Conversion I(OI) (erg s⁻¹ sr⁻¹ cm⁻²) = 1.0943 10⁻⁴ I(OI) (K km s⁻¹). (b)In units of (10⁻⁵ erg s⁻¹ sr⁻¹ cm⁻²). Conversion I(CII) (erg s⁻¹ sr⁻¹ cm⁻²) = 7.0354 10⁻<fnref n="6" /> I(CII) (K km s⁻¹). Conversion I(CO16-15) (erg s⁻¹ sr⁻¹ cm⁻²) = 6.3953 10⁻<fnref n="6" /> I(CO16-15) (K km s⁻¹). Conversion I(CO11-10) (erg s⁻¹ sr⁻¹ cm⁻²) = 2.0835 10⁻<fnref n="6" /> I(CO11-10) (K km s⁻¹). (c)In units of (erg s⁻¹ sr⁻² cm⁻²). (d)The total velocity range is from −30 to +25 km s⁻¹.

data. The line profiles are very complex, showing several velocity components and extended wing emission, particularly in the blue range. Some of the spectral features can be explained by self-absorption (see below). Figure 5 displays the average spectrum (across the 40'' × 40'' area indicated in Fig. 1) of [O I] together with other tracers such as [C II] and molecular lines. Based upon the [O I] line profile, the emission from other tracers, and earlier studies of S106 (Schneider et al. 2002, 2003) we define and label five distinct velocity ranges to characterize the velocity distribution in S106. We note that the terms “high-velocity” and “outflow” used here are descriptive, comprising all gas dynamics (molecular and atomic) provoked by stellar feedback (radiation, wind, shocks). In the following sections, we single out the possible origin of each of the different observed features.

- *High-velocity blue emission (HV-blue)*. v = −30 to −9 km s⁻¹, significant emission only observed in [O I] and [C II], no emission in optically thin lines.

- *Blue outflow emission*. v = −9 to −4 km s⁻¹, prominent emission in all optically thick lines in the form of extended wings.

- *Bulk emission*. v = −4 to 0.5 km s⁻¹, commonly defined velocity range for the bulk emission of the associated molecular cloud, peak emission for all observed lines, substructure in line profiles for optically thick lines, strong self-absorption features in the [O I] line.

- *Red outflow emission*. v = 0.5–8 km s⁻¹, broad line wings for CO lines, individual components in [C II] and [O I]. At +3 km s⁻¹, there is a prominent dip in the [O I] line and in all other line tracers that is due to absorption<fnref n="6" />.

- *High-velocity red emission (HV-red)*. v = 8–25 km s⁻¹, significant emission only observed in [C II] and [O I].

<footnote n="6">⁶ As already discussed and shown in more detail in Schneider et al. (2003), the dip seen in [O I], [C II], and $^{12}$CO 2→1 spectra gets filled in with emission in the corresponding $^{13}$CO 2→1 spectra. The ratio $^{12}$CO/$^{13}$CO drops to values lower than 1 for this velocity range over a large area on the western side of S106 IR (see Fig. 8 in Schneider et al. 2003). The lack of $^{12}$CO (and [O I] and [C II]) emission around +3 km s⁻¹ is thus due to absorption of these lines becoming optically thick in cold foreground material. This cold gas originates from a molecular clump, associated with the extended molecular cloud.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Spectral map of [O I] emission with a velocity resolution of 1 km s⁻¹ on a 6'' beam-sampled grid in the velocity range from −40 to 40 km s⁻¹ and main beam brightness temperature range from −4 to 40 K.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Spatially averaged spectrum of molecular and atomic lines (the spectral resolution is typically 0.2–0.6 km s⁻¹). All spectra within the central area of 40'' × 40'' around S106 IR were included. Based upon the different spectral features, we define five major velocity ranges. The local standard of rest (lsr) velocity of the S106 molecular cloud is approximately −3 km s⁻¹.</figcaption></figure>

## 4.2. Overview of [O I], [C II], and CO emission in S106

Because the line profiles of [O I] and other atomic and molecular line tracers are very complex and show a strong positional variation, we show the spectra for the positions where we performed PDR modelling in Fig. 7, an [O I] channel map (Fig. 8), and overlays of different tracers in the main velocity ranges (Appendix B). The spectra are extracted at the peak of emission in different line tracers in the respective velocity range and reflect different physical environments. They are indicated in Fig. 6 where we plot the overall line integrated [O I] emission as contours on a dust map at 19 µm (Adams et al. 2015). This overlay indicates that the total [O I] emission correlates very well with the warm dust, while the most prominent peaks of [O I] emission (positions 2, 3a, 3b) are slightly shifted with regard to the peaks of 19 µm emission. The line integrated intensities and the total

FIR continuum (used for PDR modelling) are listed in Table 2. We note that positions 3a and 3b appear more than once in the table because there are peaks of emission in [O I] at different velocities.

**High-velocity blue emission.** The most prominent feature of the [O I] spectra is the high-velocity blue (HV-blue) emission that appears as a single component with an extended blue wing at temperatures up to 20 K at positions 1, 2, and 3c in Fig. 7. The dip in emission at −9 km s⁻¹ is most likely not caused by self-absorption because the [C II] line has a very similar emission profile (though it cannot be excluded that both lines are affected by self-absorption). The channel map (Fig. 8) clearly shows how the HV-blue emission starts south-west of S106 IR at $v = -20.6$ km s⁻¹ and then gradually develops a north-eastern component. Both peaks together are best visible at $v = -12.6$ km s⁻¹.

**Blue outflow emission.** The channel map shows how the northern and southern cavity lobes are outlined in [O I] emission in the blue outflow velocity range (panels −8.6 to −4.6 km s⁻¹). Position Pos. 3c represents the peak of blue outflow emission in the CO 16→15 line, and Pos. 4 is the peak position of the blue outflow in [O I] emission in the northern cavity.

**Bulk emission.** Generally, the bulk emission velocity range shows a very clumpy distribution around S106 IR (panels −2.6 and −0.6 km s⁻¹). The “hole” in the emission around S106 IR is partly caused by self-absorption, but also reflects a real lack of gas since all molecular line tracers show no emission (see Sect. 5.3). The [O I] spectra (Fig. 7) are very complex for this velocity range, showing mostly several components, which can be due to several PDR layers on clump surfaces along the line of sight and/or intrinsic self-absorption effects. Future observations of the more optically thin [O I] 145 µm line may help to decide between the different possibilities. The [O I] 145/63 µm ratio determined from velocity unresolved data
<figure><img src="image_6.png" /></figure>

**Fig. 6.** [O I] emission of the total velocity range from −30 to 25 km s⁻¹ in contours overlaid on a map of 19 µm emission (FORCAST; Adams et al. 2015). Positions for which we show individual spectra of various lines and report the line integrated intensities in Table 2 (in a 15'' beam) are indicated with a cross and labelled 1–6. The angular resolution of the 19 µm emission and the [O I] beam are indicated in the panel. For comparison, the 15'' resolution we used for PDR modelling is also given. Positions 1–6 do not always correspond to peaks in [O I] emission, but also to peak emission in high- and low-J CO lines (see text for details).

(Schneider et al. 2003; Stock et al. 2015) yield a value of 0.15–0.17, indicating that one or both of the lines is optically thick. We note that the molecular cloud velocity is −1 km s⁻¹, indicated by the ¹³CO 2→1 line in (Fig. 7), while the [O I] line has a centre velocity of typically −4 km s⁻¹. A very prominent peak of emission is found at position Pos. 3b which arises from a single clump SW of S106 IR (position S106 FIR, location of the H₂O maser), clearly visible in panel −2.6 km s⁻¹ in the channel map. Position 5a represents the “dark lane”, which is strong in ¹³CO emission and weaker in the atomic and high-J CO lines, and 5b is the NW peak of [O I] emission.

**Red outflow.** The velocity range from 0.5 to 8 km s⁻¹ is strongly affected by absorption, basically all [O I] spectra show a dip between 1 and 3 km s⁻¹ (see footnote on page 6). Significant emission is only found around $v = 8$ km s⁻¹ in the southern cavity at Pos. 3a.

**High-velocity red emission.** The [O I] emission at velocities higher than $v = 8$ km s⁻¹ is very different from the HV-blue emission; it is more diffuse and extends further away from S106 IR, mostly outlining the eastern cavity walls. Position 6 represents the [O I] peak in the northern cavity.

### 4.3. Shocks in S106

We did not detect the SiO (2→1) or (5→4) lines, classical tracers of stationary C-type (e.g. Schilke et al. 1997, Gusdorf et al. 2008a), J-type (e.g. Guillet et al. 2009), or non-stationary CJ-type (e.g. Gusdorf et al. 2008b) shocks. However, the presence of shocks that can give rise to the [O I] emission observed at high blue and red velocities and at outflow velocities cannot be excluded. These non-detections could indeed mean that the shocks are either self-irradiated (see the study of shocks with radiative precursors by Hollenbach & McKee 1989) or externally irradiated (e.g. Lesaffre et al. 2013). In the latter case, the irradiation could come from the FUV field of the central protostars. In shocks where irradiation plays a significant role, SiO is expected to be converted to Si and Si⁺. These species can be observed through three transitions: the ³P₁⁻³P₀ at 129.68 µm, the ³P₂⁻³P₁ at 68.47 µm for the Si lines, and the ²P₃/₂⁻²P₁/₂ at 34.82 µm. The ground-state Si transition at 129.68 µm has been marginally detected by Herschel-PACS, contrary to its excited counterpart at 68.47 µm. None of the Si lines has been detected by ISO in contrast to the Si⁺ transition, which was observed. The presence of atomic and ionized silicon in the gas phase, combined with the detection of optical and near-infrared lines reported by Riera et al. (1989), hints at the presence of irradiated shocks in the S106 region. Though we briefly discuss possible shocks in Sect. 5.1.1, a more sophisticated attempt to characterize their type and the role they play in the excitation of [O I] and [C II] lies out of the scope of the present study, and will be investigated in a separate publication.

## 5. Analysis and discussion

### 5.1. Emission from different velocity ranges

In this section, we qualitatively and quantitatively discuss the emission of the observed [C II], [O I], and CO lines in the different velocity intervals. Each velocity range reflects different properties of the gas in terms of radiation field, temperature, and density, and the emission can have a different origin. In particular the [O I] and high-J CO lines can originate from PDRs, i.e. UV-heated clump surfaces and cavity walls, and/or from shocks due to disk-envelope interactions (accretion shock) and local shocks when the stellar winds hit the cavity walls. One way to determine the origin of the emission is by studying its velocity structure (see Sects. 5.1.1 and 5.1.2). The excitation energies for [C II] and [O I] emission are 91 and 228 K, respectively. These transitions are easily excited by collisions with H and H₂ in a PDR with different temperature layers. The critical densities for [O I] and [C II] are then $5 \times 10^5$ and $3 \times 10^3$ cm⁻³, respectively. The high-J CO lines (16→15 and 11→10) have higher excitation energies, ~750 and ~365 K, and when the emission is purely thermal these lines probe hot gas associated with the PDR/molecular cloud layer at critical densities of $2 \times 10^6$ and $4 \times 10^5$ cm⁻³, respectively. Though we discuss a possible origin from shocks in the following subsections, we use the observed FIR line ratios and the FIR flux only for PDR modelling.

The observed line intensities in the different velocity ranges and the total line integrated intensity and the fraction of intensity are listed in Table 2. All results from PDR modelling for the individual velocity ranges are summarized in Fig. 9 and Fig. 10 where we show the KOSMA-$\tau$ results for the [CII]₁₅₈µm/[OI]₆₃µm and CO(16–15)/CO(11–10) line ratios as iso-contours. Dashed lines show contours of constant FUV
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Individual spectra with 15" angular resolution at the positions indicated in Fig. 5.</figcaption></figure>

fields, while black contours represent lines of constant density. In Tables 3 and 4, we summarize the observed ratios and give the model parameters corresponding to their position in the diagnostic Figs. 9 and 10.

Using the observed line ratios as diagnostics for PDR modelling involves several caveats which need to be kept in mind. Firstly, we use a non-clumpy PDR model because a more sophisticated set-up, as was used when studying the Orion Bar

(Andree-Labsch et al. 2017), is beyond the scope of this paper. Secondly, the [O I] 63 µm line and to a lesser extent the [C II] 158 µm line are affected by self-absorption. This becomes obvious in the line profiles (Fig. 7), and is indicated by the low [O I] 145/63 µm ratio of < 0.2 (Schneider et al. 2003; Stock et al. 2015). If the [O I] 63 µm line is only moderately optically thick ($\tau \sim 1-2$), Liseau et al. (2006) showed that self-absorption is responsible for the low ratio. We thus additionally constrained
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Channel map of [O I] emission in 2 km s⁻¹ velocity bins. The contour levels are 5, 10, 20, 45, and 60 K km s⁻¹. The logarithmic colour scale ranges from 1.5 to 60 K km s⁻¹. The star indicates the binary system S106 IR and the triangle S106 FIR.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 9.** PDR diagnostic diagrams for observed fine-structure lines based on the KOSMA-τ PDR model for a clump of M = 100 $M_{\odot}$. [CII]158 µm/[OI]63 µm against ([CII]158 µm + [OI]63 µm)/ΣFIR, where ΣFIR is the continuum intensity between 10 and 1000 µm in the respective velocity range determined as described in Sect. 3.2. The various markers indicate the different positions corresponding to different velocity ranges. The density n and the Draine field χ are given for each position.</figcaption></figure>

the PDR models using an [O I] intensity for all velocity ranges that is a factor 2, 4, and 10 higher. This is only a first approximation, but it allows us to estimate how the density and the radiation field change with varying [O I] line intensity. Using a much higher [O I] intensity (4 and 10 times the observed value) results in higher densities (typically $10^{5–6}$ cm⁻³) for all positions but in a much lower radiation field (typically χ ~ $10^{2–3}$). Such a low radiation field is not supported by observations; for example, we determine from Herschel flux measurements (see Sect. 5.2) a

value of χ ~ (2–4) × $10^4$ as a lower limit. We thus only use the results of a PDR model with the nominal [O I] intensity, and one that is a factor of 2 higher (given in Cols. 7 and 8 in Tables 3 and 4 in parentheses.)

### 5.1.1. High-velocity blue and red range

The HV-blue ($v < -9$ km s⁻¹) [O I] emission is spatially concentrated in two regions northeast (Pos. 1) and southwest (Pos. 2)
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** PDR diagnostic diagrams for observed fine-structure lines based on the KOSMA-$\tau$ PDR model for a clump of $M = 100~M_{\odot}$. CO(16–15)/CO(11–10) against [(CII)158$\mu$m + [OI]63$\mu$m)/$\Sigma_{\text{FIR}}$, where $\Sigma_{\text{FIR}}$ is the continuum intensity between 10 and 1000 $\mu$m in the respective velocity range determined as described in Sect. 3.2. The various markers indicate the different positions corresponding to different velocity ranges. The density $n$ and the Draine field $\chi$ are given for each position.</figcaption></figure>

of S106 IR (Fig. 11). Position 1 is correlated with a 19 $\mu$m dust peak, and Pos. 2 is located very close to S106 IR. It is remarkable that at Pos. 1, the [O I] emission at these high blue velocities contributes ~50% of the whole line integrated emission at this position (see Table 2). The [O I] HV-blue emission corresponds well to the [C II] emission distribution (Fig. B.1), though the latter is more extended and outlines the cavity better. The $^{12}$CO 2$\to$1 emission in this velocity range (Fig. B.1) traces the molecular outflow and the cold gas of the dark lane, but is probably more affected by beam dilution (11'' resolution with respect to 6'' for [O I]). The $^{12}$CO 16$\to$15 and 11$\to$10 lines show only very weak emission at some positions (Fig. 7).

The HV-red ($v > 8$ km s$^{-1}$) [O I] emission is more dispersed and has two peaks along the eastern limb-brightened cavity wall, and a northern peak that is not correlated with the dust. While the overall [C II] emission in this velocity range corresponds well with [O I], the $^{12}$CO 2$\to$1 line is only visible at the position of the northern clump (Fig. B.5).

There are three possible explanations for the origin of the high-velocity blue and red [O I] emission: (1) disk wind in the context of an accretion shock, (2) shocks caused by the stellar wind hitting the inside of the dark lane, or (3) classical PDR emission from the inner surface of the dark lane. Further details of each possibility are given below:

1. An accretion shock can occur because S106 IR is surrounded by a disk-like structure. However, considering the distance of less than 0.2 AU between the two stars, the existence of a circumbinary disk is very unlikely; it is probably only one disk or its remains that were seen in cm-interferometry. Adopting what is known for low-mass stars (see Hartmann et al. 2016 for a review), accretion of material from the disk takes place in small-scale (sub-pc) flows along magnetic field lines (see e.g. Fig. 1 in Hartmann et al. 2016). The material has high free-fall velocities (on the order of several hundred km s$^{-1}$), that then causes shocks at the stellar surface. However, the resulting shock heats the gas to very high temperatures (~$10^6$ K), which is far above the excitation of the [O I] 63 $\mu$m line. Moreover, the HV [O I] emission could originate

from the disk that produces a bipolar flow or jet driven by accretion energy. The observed HV-blue [O I] emission is indeed very collimated and would fit into this scenario of an atomic jet, emanating from the disk-envelope-star system. The projected velocity (towards the observer) is $v_{\text{lsr}} \sim 25$ km s$^{-1}$, so that the real velocity is $v = v_{\text{lsr}}/\sin(\alpha)$, where $\alpha$ is the inclination angle of the ionized lobes. Hippelein & Münch (1981) and Solf & Carsenty (1982) determined the angle to be around 15° so that the velocity is approximately 100 km s$^{-1}$ and thus consistent with high-velocity dissociative J-type shocks (see Hollenbach 1985; Hollenbach & McKee 1989 for theory, and Nisini et al. 2015; Leurini et al. 2015 for observations). However, there are two arguments against this shock scenario on a sub-pc scale. Firstly, the HV-red [O I] emission is not collimated in the same way as the HV-blue emission (Fig. 11) which would be expected in the case of a symmetric atomic jet. Moreover, the emission is much more extended (up to several parsecs) and correlates well with the cavity walls. Secondly, the CO 16$\to$15 line as a typical J-type shock tracer has a very low line intensity (Table 2) in this velocity range and is only visible as an extended wing (Fig. 7).

2. S106 IR consists of an OB star system that has a strong stellar wind of around 200 km s$^{-1}$ (Hippelein & Münch 1981). The wind is a source of high kinetic energy and develops a shock when it expands into the ambient medium, in this case hitting the inner part of the dark lane and the molecular cloud fragments close to S106 IR. The shocked gas can then cool via the [O I] 63 $\mu$m line.

3. Given the proximity of the exciting S106 IR binary system, the HV [O I] emission can be explained by PDRs on the surface of surrounding clumps and from the backside of the dark lane. From the PDR model (Sect. 3.2), we derive a density of $5-6 \times 10^4$ cm$^{-3}$ and an UV field of $\chi \sim 7 \times 10^4$ from the ratio [C II]/[O I] for S106 IR and the northern HV-blue and HV-red peak emission positions. The southern HV-red [O I] peak has a lower density of $10^4$ cm$^{-3}$ and an UV field of $\chi = 1.6 \times 10^4$ from [C II]/[O I]. These density values increase to $2 \times 10^6$ cm$^{-3}$ (while the radiation field becomes slightly
<table><caption>**Table 3.** Data points and model results from Fig. 9 for $[CII]_{158\mu m}/[OI]_{63\mu m}$ against $([CII]_{158\mu m} + [OI]_{63\mu m})/\Sigma_{FIRrange}$.</caption><thead><tr><th>Position</th><th>Velocity range<br>(km s<sup>-1</sup>)</th><th>Offset<br>(", ")</th><th>[CII]/[OI]<br>(×10<sup>-2</sup>)</th><th>[CII]+[OI]<br>Σ<sub>FIRrange</sub><br>(×10<sup>-3</sup>)</th><th>log n<br>(cm<sup>-3</sup>)</th><th>log χ<sub>Draine</sub></th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">High-velocity blue</td></tr><tr><td>1</td><td>-30/-9</td><td>S106 IR<br>(-2, -2)</td><td>3.4</td><td>1.24</td><td>4.8 (6.4)</td><td>4.8 (4.4)</td></tr><tr><td>2</td><td>-30/-9</td><td>northern HV-blue [O I]<br>(7, 12)</td><td>3.8</td><td>1.08</td><td>4.7 (5.0)</td><td>4.9 (4.6)</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Outflow blue</td></tr><tr><td>3a</td><td>-9/-4</td><td>southern cavity blue [O I]<br>(2, -20)</td><td>19.4</td><td>1.68</td><td>4.1 (4.5)</td><td>4.0 (4.0)</td></tr><tr><td>3b</td><td>-9/-4</td><td>S106 FIR<br>(-15, -5)</td><td>13.3</td><td>1.45</td><td>4.3 (4.6)</td><td>4.3 (4.2)</td></tr><tr><td>3c</td><td>-9/-4</td><td>CO 16→15 peak<br>(8, 8)</td><td>10.5</td><td>1.10</td><td>4.3 (4.6)</td><td>4.5 (4.4)</td></tr><tr><td>4</td><td>-9/-4</td><td>northern cavity blue [O I]<br>(0, 15)</td><td>7.0</td><td>1.26</td><td>4.5 (4.8)</td><td>4.6 (4.4)</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Cloud bulk emission</td></tr><tr><td>3b</td><td>-4/0.5</td><td>S106 FIR<br>(-15, -5)</td><td>10.1</td><td>1.41</td><td>4.4 (4.7)</td><td>4.4 (4.3)</td></tr><tr><td>5a</td><td>-4/0.5</td><td>dark lane<br>(10, -5)</td><td>20.8</td><td>1.31</td><td>4.1 (4.4)</td><td>4.1 (4.1)</td></tr><tr><td>5b</td><td>-4/0.5</td><td>NW bulk [O I]<br>(-10, 15)</td><td>7.8</td><td>1.67</td><td>4.5 (4.8)</td><td>4.4 (4.2)</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Outflow red</td></tr><tr><td>3a</td><td>0.5/8</td><td>southern cavity red [O I]<br>(2, -20)</td><td>8.1</td><td>1.52</td><td>4.4 (4.8)</td><td>4.5 (4.3)</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">High-velocity red</td></tr><tr><td>6</td><td>8/25</td><td>northern HV-red [O I]<br>(10, 25)</td><td>3.1</td><td>1.54</td><td>4.8 (6.3)</td><td>4.7 (4.4)</td></tr><tr><td>3a</td><td>8/25</td><td>southern HV-red [O I]<br>(2, -20)</td><td>14.3</td><td>1.61</td><td>4.2 (4.6)</td><td>4.2 (4.1)</td></tr></tbody></table>

**Notes.** The abbreviation “range” signifies that we scale the total FIR flux with the fractional [O I] intensity in this velocity range (see text for further details). The values for density *n* and radiation field *χ* (Colss. 7 and 8) are given for the nominal observed [O I] intensity and for an [O I] intensity increased by a factor of two (in parentheses).

<table><caption>**Table 4.** Data points and model results from Fig. 10 for CO(16–15)/CO(11–10) against $([CII]_{158\mu m} + [OI]_{63\mu m})/\Sigma_{FIRrange}$.</caption><thead><tr><th>Position</th><th>Velocity range<br>(km s<sup>-1</sup>)</th><th>Offset<br>(", ")</th><th>CO(16-15)/CO(11-10)</th><th>[CII]+[OI]<br>Σ<sub>FIRrange</sub><br>(×10<sup>-3</sup>)</th><th>log n<br>(cm<sup>-3</sup>)</th><th>log χ<sub>Draine</sub></th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">High-velocity blue</td></tr><tr><td>1</td><td>-30/-9</td><td>S106 IR<br>(-2, -2)</td><td>0.76</td><td>1.24</td><td>4.4 (4.6)</td><td>4.6 (4.3)</td></tr><tr><td>2</td><td>-30/-9</td><td>northern HV-blue [O I]<br>(7, 12)</td><td>0.27</td><td>1.08</td><td>4.2 (4.3)</td><td>4.4 (4.1)</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Outflow blue</td></tr><tr><td>3a</td><td>-9/-4</td><td>southern cavity blue [O I]<br>(2, -20)</td><td>0.74</td><td>1.68</td><td>4.5 (4.6)</td><td>4.4 (4.1)</td></tr><tr><td>3b</td><td>-9/-4</td><td>S106 FIR<br>(-15, -5)</td><td>0.84</td><td>1.45</td><td>4.5 (4.7)</td><td>4.5 (4.2)</td></tr><tr><td>3c</td><td>-9/-4</td><td>CO 16→15 peak<br>(8, 8)</td><td>1.44</td><td>1.10</td><td>4.7 (4.9)</td><td>4.9 (4.6)</td></tr><tr><td>4</td><td>-9/-4</td><td>northern cavity blue [O I]<br>(0, 15)</td><td>1.13</td><td>1.26</td><td>4.6 (4.8)</td><td>4.7 (4.4)</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Cloud bulk emission</td></tr><tr><td>3b</td><td>-4/0.5</td><td>S106 FIR<br>(-15, -5)</td><td>0.42</td><td>1.41</td><td>4.3 (4.4)</td><td>4.4 (4.1)</td></tr><tr><td>5a</td><td>-4/0.5</td><td>dark lane<br>(10, -5)</td><td>0.31</td><td>1.31</td><td>4.2 (4.3)</td><td>4.3 (4.0)</td></tr><tr><td>5b</td><td>-4/0.5</td><td>NW bulk [O I]<br>(-10, 15)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Outflow red</td></tr><tr><td>3a</td><td>0.5/8</td><td>southern cavity red [O I]<br>(2, -20)</td><td>1.74</td><td>1.52</td><td>4.9 (7.0)</td><td>4.9 (3.2)</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">High-velocity red</td></tr><tr><td>6</td><td>8/25</td><td>northern HV-red [O I]<br>(10, 25)</td><td>3.52</td><td>1.55</td><td>6.3 (6.3)</td><td>6.0 (6.0)</td></tr><tr><td>3a</td><td>8/25</td><td>southern HV-red [O I]<br>(2, -20)</td><td>2.11</td><td>1.61</td><td>5.3 (5.4)</td><td>4.8 (4.5)</td></tr></tbody></table>

**Notes.** The values for density *n* and radiation field *χ* (Colss. 7 and 8) are given for the nominal observed [O I] intensity and for an [O I] intensity increased by a factor of two (in parentheses).
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** [O1] line integrated emission overlayed as contours on 19 µm emission from dust (Adams et al. 2015). The left panel shows the high-velocity blue emission and the right panel the high-velocity red emission. The dust emission ranges between 0 and 2.6 Jy pix⁻¹, the [O1] contour lines go from 30 to 100% of maximum intensity (397 K km s⁻¹ for HV-blue and 121 K km s⁻¹ for HV-red).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** [O1] line integrated emission overlayed as contours on cm emission from the VLA (Bally et al. 1983). The left panel shows the blue outflow range and the right panel the red emission. The [O1] contour lines go from 36 to 156 in steps of 10 K km s⁻¹ for the blue outflow and 40 to 180 in steps of 10 K km s⁻¹ for the red outflow. The star indicates the binary system S106 IR and the triangle S106 FIR.</figcaption></figure>

lower, typically $\chi$ around a few times $10^4$) when the [O1] intensity is increased by a factor of two. The density and radiation field obtained from the CO(16→15)/CO(11→10) line ratio for the HV-blue emission are lower ($n$ = 1.5 and $2.8 \times 10^4$ cm⁻³ and $\chi$ = 2–4 × $10^4$) and stay below $10^5$ for the density even if the [O1] intensity increases. Modelling the HV-red CO emission is tricky because of the low line intensities (see Table 2). For Pos. 6, we derive a density of $n \sim 2 \times 10^6$ cm⁻³ for both original and increased [O1] intensity, but the radiation field with a value of $\chi = 10^6$ is uncharacteristically high.

On first sight, it is puzzling that in the high-velocity blue range we observe strong [O1] emission but nearly no CO 16→15 emission, even though the two lines have similar critical densities.

Because CO 16→15 has a higher excitation temperature (750 K) than [O1] (228 K), and emits much less prominently in this velocity range, we assume that this high-velocity gas is not strongly heated and consists only of gas with high densities ($>10^6$ cm⁻³) at T ~ 300 K. The CO 16→15 line is then mostly subthermally excited.

### 5.1.2. Outflow blue and red [O1] emission

The [O1] emission in the blue and red outflow velocity range (~9 to ~4 km s⁻¹ and 0.5 to 8 km s⁻¹, respectively) outlines precisely the cavity of the bipolar nebula, illustrated in Fig. 12, which shows [O1] contours overlaid on a cm VLA image (Bally et al. 1983), and in Fig. B.2. The cm emission delineates the
<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Optical image from the Gran Telescopio Canaries (colour) with an overlay of CO 16→15 line integrated emission (contour levels 18.6 to 83.8 in steps of 6.5 K km s⁻¹) in the blue outflow velocity range (−9 to −4 km s⁻¹). The dark lane is indicated by one black contour (3 Jy) of SHARC 350 µm emission (see Fig. 13).</figcaption></figure>

H II region but avoids the dark lane. Bally et al. (1983) argue that the density in the lane is so high that all ionizing radiation is absorbed (see Sect. 5.3). Similar to the [C II] emission distribution (Simon et al. 2012), the red [O I] outflow emission is absent in the northern lobe. Both [O I] and [C II] red outflow emission thus stem mostly from swept-up gas of the cavity walls from the backside of the southern lobe (since the front side does not exist). In addition, we observe limb-brightened emission at the cavity walls and, closer to S106 IR, the [O I] emission in the blue outflow range is also correlated with the clump containing S106 FIR. What is causing the gas dynamics is not clear. The stellar wind of S106 IR can provoke shocks when it hits (mechanically) the cavity walls. However, the radiation of S106 IR ionizes the interface layer between the H II region and the molecular cloud and creates a classical PDR. A mixed scenario involving “irradiated shocks” (Lesaffre et al. 2013) and/or “moving PDRs” (Störzer & Hollenbach 1998) is also possible. From our modelling, we obtain average densities at the three cavity positions (3a,b, and 4) of 1.3–5.6 × 10⁴ cm⁻³ at a rather constant radiation field of χ ~ a few 10⁴, depending on position and ratio ([C II] and [O I] or CO). Interestingly, the values for density and radiation field do not change much if we increase the [O I] intensity. Though the [C II] line and to a lesser extent the [O I] line can also be excited by collisions with electrons, we assume that only a very small fraction of [C II] emission arises from the H II region because we detected only very weak [N II] 205 µm emission. Generally, this line arises only from the H II region, so that its detection only together with a sufficient [N II] to [C II] line ratio leads us to conclude that a significant amount of [C II] also emerges from the ionized gas phase (Heiles 1994). This is not the case for S106.

In contrast to [O I], the emission distribution of the CO 16→15 line (Fig. 13) in this velocity range shows a strong spatial correlation with the “dark lane”, peaking at position 3c, which shows the highest values of CO 16→15 emission (Table 2). Figure B.2 shows that the peak of emission in the low-J CO lines (¹²CO and ¹³CO 2→1) is slightly shifted towards S106 IR. Because the lane appears prominently in the optical

<figure><img src="image_14.png" /><figcaption>**Fig. 14.** SHARC 350 µm emission (colour wedge in Jy) with overlays of CO 16→15 line integrated emission (levels 14 to 47 by 3 K km s⁻¹, grey contours) in the bulk velocity range (−4 to 0.5 km s⁻¹).</figcaption></figure>

image as a dark feature and consists of cold gas (see SHARC image, Fig. 14, and overlays with the low-J CO lines, Fig. B.2), the CO 16→15 emission can only arise from a PDR layer of warm and dense gas on the backside of the lane towards S106 IR. The rather low line velocities and small line width argue against a high-velocity shock origin of the high-J CO emission, but low-velocity shocks cannot be excluded. From PDR modelling, using the CO 16→15/11→10 ratio at position 3c, we derive a density of 5.6 × 10⁴ cm⁻³ and a radiation field of χ = 7.4 × 10⁴. For the red outflow velocity range (0.5–8 km s⁻¹) all lines with high excitation temperature and high critical density ([O I], [C II], CO 16→15) show a similar emission distribution (Fig. B.4) with an extended clump south of S106 IR. This velocity range is affected by absorption so that the derived line intensities and ratios are not reliable and we do not use those for PDR modelling.

Summarizing, we interpret the emission in the outflow velocity ranges for all lines as arising from PDR surfaces of heated and UV illuminated clumps around S106 IR, of ablated gas from the cavity walls, and from a PDR surface on the backside of the dark lane.

### 5.1.3. Bulk emission

The velocity range from −4 to 0.5 km s⁻¹ is dominated by widespread but clumpy [O I] emission (see channel map, Fig. 8) in which the area around S106 IR is mostly devoid of emission (best visible in panel *v*=-0.6 km s⁻¹). The most prominent clump is the western one associated with S106 FIR. Two other clumps northwest and south of S106 IR (Fig. B.3) are also clearly defined, and the overall emission distribution corresponds very well with that of [C II] in this velocity range (Fig. B.3). The peaks of emission in the low-J CO lines and CO 16→15 emission are close to the dark lane and at the position of S106 FIR (Figs. 14 and B.3). The CO 16→15 peak corresponds to Pos. 5a, the NW peak of [O I] emission to Pos. 5b, and the western [O I] peak to Pos. 3b, i.e. the continuum source S106 FIR (Richer et al. 1993; Little et al. 1995), possibly a low- or high-mass Class 0 YSO. The absolute line intensities of [O I] and [C II] emission are high and the line fractional emission is 20–40%, indicating that these lines are the most important cooling lines for this portion of the gas. The CO 16→15 line is less strong at these positions (in contrast to the CO 11→10 line) and excited only in the outermost surface PDR layer of the dark lane.
facing S106 IR. Because the density and radiation field are not
very high (see below), it is a fair assumption that this line is
subthermally excited.

From PDR modelling, we derive a density of 1.1 (1.6) 10⁴ cm⁻³ (from the [CII]/[OI] and CO ratios, respectively), and a radiation field of χ = 1.1 (2.1) 10⁴ for Pos. 5a, i.e. the CO 16→15 peak at the inner working surface of the dark lane. The density and radiation field (n = 3.0 × 10⁴ cm⁻³ and χ = 2.7 × 10⁴) are similar at the position of the NW [OI] peak. S106 FIR shows peak emission in hot dust, but the gas density and radiation field is not elevated, i.e. n ~ 2 × 10⁴ cm⁻³ and χ ~ 2.5 × 10⁴. As for the blue and red outflow velocity range, the density and radiation field do not change significantly when the [OI] intensity is increased in the PDR model.

In summary, the [O I] emission in this velocity range seems
to arise from the inner cavity walls in combination with PDR
surfaces of individual clumps located within the cavity, close to
S106 IR.

## 5.2. Comparison between SOFIA results and earlier studies

This is the first time that a spectrally resolved map of the [O I]
63 µm FIR cooling line is available for S106. A [C II] 158 µm
map obtained with GREAT on SOFIA was presented by Simon
et al. (2012). Other observations (not spectrally resolved) of
the [C II] and [O I] lines with the Infrared Space Observatory
(ISO), the Kuiper Airborne Observatory (KAO), and Herschel
were reported by van den Ancker et al. (2000), Schneider et al.
(2003), Stock et al. (2015), respectively. These studies, how-
ever, have a much lower angular resolution (40" to 80") and
thus contain the different physical regimes in S106 (H II region,
molecular cloud bulk emission, cavity, etc.) in one beam.
Nevertheless, all studies agree that there are different gas phases
with at least two temperature regimes. Schneider et al. (2003)
proposed a two-phase gas model with small (<0.2 pc) high-
density (n = 3 × 10⁵ to 5 × 10⁶ cm⁻³) clumps with high surface
temperature (T ~ 200–500 K) embedded in lower density
(n ~ 10⁴ cm⁻³) gas. The radiation field in these two regimes is
χ ~ 10⁵ and χ ~ 10³, respectively. Stock et al. (2015), based on
their PACS and SPIRE Herschel study, support this scenario.
They find a two-phase regime with a hot (T ~ 400 K) region
characterized by a high radiation field (χ > 10⁵) and high density
(n > 10⁵ cm⁻³), based on CO rotational temperatures, and a less
dense region (n ~ 10⁴ cm⁻³) with moderate UV field (χ ~ 10⁴)
and a temperature of around 300 K, dominating the emission of
the atomic cooling lines. They argue that a third PDR gas phase
may be present with very high temperatures/radiation field and
densities to explain their excess emission in the CO rotational
ladder for the highest-J CO lines (J > 20). From fitting the high-
J CO lines within a PDR model, they obtain a radiation field of
χ ~ 10⁵ and densities of 10⁸ cm⁻³. The assumption of such high
densities is somewhat problematic as they are not observed in
any other way, thus excitation by shocks can also be the reason
for the strong CO emission. However, the contribution of shocks
is not clear. van den Ancker et al. (2000) and Stock et al. (2015)
estimate that their contribution is small, while Noel et al. (2005)

attribute their H₂ observations to shocks. In a subsequent paper
on S106, we will apply irradiated shock models and investigate
whether they reproduce our observed line intensities and ratios
more closely.

In this study, we explain the observed [O I], [C II], and CO line ratios with a PDR model of a single gas phase with densities of a few 10⁴ cm⁻³, exposed to a radiation field χ of a few 10⁴. This would be consistent with the lower density PDR component that was found by Schneider et al. (2003) and Stock et al. (2015). If we consider self-absorption of the [O I] line and approximate the missing emission by doubling the [O I] intensity, we find that the HV-blue and HV-red emission has its origin in much denser gas of n ~ 10⁶ cm⁻³, but with a similar radiation field<fnref n="8" /> of χ = 2–3 × 10⁴. This can correspond to the high-density gas component found by van den Ancker et al. (2000), Schneider et al. (2003), and Stock et al. (2015). Because they are velocity unresolved observations, it was not possible to attribute one gas component to a certain velocity range. Our data now point toward a possible scenario in which the FIR line emission in the outflow and bulk velocity ranges arises from PDRs at lower density (a few 10⁴ cm⁻³) exposed to a radiation field of a few χ = 10⁴, while the high-density gas component is only found at the high blue and red velocities. This finding emphasizes the importance of velocity resolved observations of [O I], [C II], and high-J CO lines.

Our findings also show that we do not need to invoke a third,
very high-density gas phase (Stock et al.) nor a high-temperature
two-gas phase regime (at T ~ 300 K and T ~ 700 K) as was
needed to model outflow emission from low-mass protostars
(Kristensen et al. 2017).

## 5.3. The dark lane: an accretion flow or an evaporating filamentary structure?

It was shown by dust observations in the mm-wavelength range
(Vallée & Fiege 2005; Motte et al. 2007), in the mid-IR with
FORCAST (Adams et al. 2015), and in the FIR with SHARC
(Simon et al. 2012) at 350 µm, that the dark lane east of S106 IR
is a high column-density feature (see Fig. 14 for the SHARC
350 µm map). From a column density map derived from a SED
fit to the 160–500 µm flux maps of Herschel (Schneider et al.
2016), we obtain an average column density<fnref n="9" /> of 4.5 × 10²² cm⁻²,
a total mass of 275 M☉, and an average density of 3 × 10⁴ cm⁻³
for the lane.

Our IRAM 30 m observations of H¹³CO⁺ 1→0 reveal the dynamics of the dark lane more precisely than earlier ¹²CO and ¹³CO data of the region (Schneider et al. 2002). Figure 15 shows a channel map of H¹³CO⁺ 1→0 emission in the velocity range between -1.3 and -3.2 km s⁻¹ overlaid on the Subaru IR image (Fig 1). The emission distribution has a close link to S106 IR. It starts as a more widespread emission around -1.3 km s⁻¹ with two lobes of emission and develops into a single feature between -2.3 and -3.2 km s⁻¹. The peak of emission moves gradually from a position further away northeast of S106 IR at -1.8 km s⁻¹ to a peak close to S106 IR at -3.2 km s⁻¹. This velocity distribution is thus consistent with two possible geometries that we illustrate in Fig. 16. Either the far end of the lane is slightly tilted away from the observer and the gas flows towards

<footnote n="7">⁷ As a consistency check, we determined the total [O I] 63 µm and [C II] 158 µm line intensity in a 40'' beam around the position of S106 IR and compared it to the value obtained with Herschel (Stock et al. 2015). Our line intensities are slightly lower but agree well, considering the calibration uncertainties for both instruments (15–20%). Our values of [O I] and [C II] emission are 2.2 × 10⁻² and 3.6 × 10⁻³ erg cm⁻² sr⁻¹ s⁻¹, respectively, compared to 2.9 × 10⁻² and 3.9 × 10⁻³ erg cm⁻² sr⁻¹ s⁻¹ from Stock et al. (2015).</footnote>

<footnote n="8">⁸ An independent measure of the radiation field is to use the Herschel fluxes at 70 and 160 µm (see e.g. Schneider et al. 2016). With this method, we obtained a field of typically χ = 2–4 × 10⁴ or S106.</footnote>

<footnote n="9">⁹ The visual extinction in the lane has been estimated between A_ν = 12 and 21 mag (Eiroa et al. 1979; Hodapp & Rayner 1991; van den Ancker et al. 2000).</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. 15.** Channel map of H$^{13}$CO$^{+}$ 1→0 emission (contours in K km s$^{-1}$ from 0.25 to 1.75 in steps of 0.25) obtained with the IRAM 30 m telescope. The background image is the Subaru IR image from Fig. 1. The blue contour is the 3 Jy level of the SHARC 350 µm emission. The star marks the position of the binary system S106 IR; the dashed red line indicates the run of the dark lane, i.e. the possible accretion flow.</figcaption></figure>

S106 IR (the accretion flow scenario, in green in Fig. 16) or the far end of lane is oriented towards the observer and the gas is streaming off the lane (the dispersal scenario, in blue in Fig. 16)$^{10}$. In both scenarios, the lane is located in front of the HII region because it is clearly visible as a dark feature in front of a bright HII region (e.g. Fig. 15). While this is straightforward to understand in the dispersal scenario because the lane is inclined towards the observer, the flow scenario requires that the lane should be only slightly tilted away from the observer so that it can wrap around the equatorial plane of the hourglass-shaped HII region, but still be located in front of the ionized gas phase.

The dispersal scenario was originally proposed by Hodapp & Rayner (1991). They explain the bipolar HII region as arising from an anisotropic circumstellar wind, absorbing the ionizing radiation in the equatorial plane, and postulate that the single massive star S106 IR must be in a mass-loss evolutionary state and not in an accretion phase. However, such a scenario requires an energy transfer of the stellar wind onto the dense and cold gas deep in the dark lane. We observe in [O I] and [C II] high-velocity gas around –30 km s$^{-1}$ that stems from PDRs and possibly shocks, caused by the stellar wind and radiation hitting the inner working surface of the dark lane. However, the cold gas within the dark lane occurs around the bulk emission of the cloud, from –1 to –3 km s$^{-1}$, and is thus decoupled from the high-velocity gas phase.

In contrast, an accretion flow scenario$^{11}$ was proposed by Bally et al. (1983) based on their VLA 5 GHz data (see Fig. 12). The area just around S106 IR shows only very weak 5 GHz radio emission, indicating the absence of ionized gas (because the 5 GHz emission is not attenuated by dust). The gas at very high densities is then able to absorb the emitted ionizing radiation. Our data show that the lane can be a flow that is illuminated by S106 IR.

Both scenarios were also proposed in Balsara et al. (2001), based upon interferometric $^{13}$CO observations. They emphasize the importance of the magnetic field that can channel material in filamentary structures onto the highest gravitational potential well. The magnetic field was found to run parallel to the dark lane (Vallée & Fiege 2005), thus consistent with an inflow of gas. A number of recent observational studies support this dynamic scenario (e.g. Schneider et al. 2010; Kirk et al. 2013; Peretto et al. 2013; Rayner et al. 2017). Despite the higher angular resolution of their $^{13}$CO data compared to our H$^{13}$CO$^{+}$ data at 30'' resolution, it is not possible to discriminate between the two views, mostly because a low-J $^{13}$CO line is not the best tracer for dense gas. It is also not possible to tell whether accretion is still occurring. Interferometric observations are required with a high-velocity resolution of high-density tracers such as H$^{13}$CO$^{+}$, H$^{13}$CN, and in particular N$_2$H$^{+}$ (a probe for cold dense gas). This

<footnote n="10">¹⁰ See also the animated version of this cartoon available online.</footnote>

<footnote n="11">¹¹ We talk here of a parsec-scale flow that connects the molecular cloud with the disk(s) or the remains thereof, and not of possible sub-pc flows that channel material from the disk(s) onto the star (Sect. 5.1.1).</footnote>
<figure><img src="image_16.png" /></figure>

**Fig. 16.** Schematic view of two possible scenarios (shown at the same time) explaining the nature of the dark lane. The front view (**left panel**), as S106 is seen on the sky, does not allow us to distinguish between infall scenario (dark lane in green) or expansion scenario (dark lane in blue) of the gas in the dark lane. The middle view shows best the “accretion flow” scenario, i.e. how the (green) dark lane wraps around the equatorial waist of the hourglass-nebula. The right view shows how the (blue) dark lane, which is more detached from the nebula, is tilted towards the observer and dispersed. An animated version of this cartoon is available [online](#).

sort of observation would enable us to study the velocity structure of the dark lane and to explore the immediate environment around S106 IR.

Assuming that we indeed observe a flow towards S106 IR, the projected length *l* of the flow is roughly 2.2', corresponding to ~1 pc at a distance of 1.3 kpc. The velocity difference $\Delta v$ along the flow, determined from the H$^{13}$CO$^{+}$ map, is ~1.4 km s$^{-1}$. We assume a random distribution of orientation angles and thus an average angle of the flow to the line of sight of 57.3° so that the lifetime *t* of the flow calculates as $t = l/(\Delta v \tan(57.3))$. The lifetime is then approximately $1.1 \times 10^6$ yr, leading to a mass input rate of $2.5 \times 10^{-4} M_{\odot} \text{yr}^{-1}$ (using the *Herschel* determined mass of the lane). This rate is lower than that found for the massive subfilaments of the DR21 ridge (a few $10^{-3} M_{\odot} \text{yr}^{-1}$, Schneider et al. 2010), but similar to that observed for Mon R2 (a few $10^{-4} M_{\odot} \text{yr}^{-1}$; Rayner et al. 2017). We note, however, that this is a very crude approximation, and does not consider the true inclination and bending of the flow.

The flow scenario is consistent with simulations of Peters et al. (2010a,b) who modelled the collapse of rotating, massive cloud cores including radiative heating by both non-ionizing and ionizing radiation. The simulations show fragmentation from gravitational instability in the dense accretion flows in which either a single massive star is formed or several massive stars with many low-mass stars within the regions of the accretion flows. Because there is competition of mass, but not in the classical meaning of “competitive accretion” (Bonnell 2007; Bate 2012), they named this process “fragmentation-induced starvation”. The simulations (Peters et al. 2010b) also show the non-uniform expansion of H II regions with the formation of bipolar lobes, very similar to what is observed for S106. The bipolar outflow structure of the H II region visible at early time steps of their runs A (single sink) and B (multiple sinks) reflects qualitatively what is seen for S106. From what is known so far, S106 IR is a binary system (Comerón et al. 2018) and is associated with a large cluster of low-mass stars, detected in the IR (Hodapp & Rayner 1991). In addition, several dense pre-stellar cores have formed in the dark lane (S20 to S24 in Motte et al. 2007), indicating ongoing fragmentation.

In the scenario of the dark lane as an ionized accretion flow in the fragmentation-induced starvation model, the observed emission distribution of [O I], [C II], and other tracers fit very well. In particular the high-velocity blue emission is consistent with

gas that moves downward perpendicularly to the accretion flow, down the steepest density gradient. This causes the very focused emission seen in various tracers (Figs. 11 and B.1). Interestingly, the very complex morphology and velocity distribution of jets and outflowing gas – as we observe it for S106 – is well reproduced in the models of Kuruwita et al. (2017). They produced simulations of the outflow pattern of close (separation < 10 AU) and wide (>10 AU) binaries and showed that the geometry and velocity of the jets and outflows are strongly modified with respect to a single star.

## 6. Summary: a scenario for S106 IR

Summarizing our observational results and what is already known about S106, we develop the following scenario:

S106 IR is a binary system with two stars that are sources of a strong UV field and stellar wind. The very small disk-like structure that was detected in cm-interferometry can be an intact accretion disk that is connected to a large-scale accretion flow, known as the dark lane or the remains of a disk without a link to the dark lane. The lane is illuminated by the more massive star of the system, presumably an O9 star with $20 M_{\odot}$ (the companion has a preliminary classification as a B8 star with ~3 $M_{\odot}$; Comerón et al. 2018), and forms a dense hot PDR that cools mostly via [O I] 63 µm and CO 16→15 emission. The PDR gas is highly dynamic; it flows fast and follows the steepest density gradients of the dark lane, and “escapes” the lane close to S106 IR, giving rise to the very collimated emission distribution in the [O I] lines at velocities from -30 to -9 km s$^{-1}$ and 8 to 25 km s$^{-1}$. We obtain a radiation field of $\chi = 2-7 \times 10^4$ at a density of $1.5-6 \times 10^4$ cm$^{-3}$ from PDR modelling. Generally, modelling the [O I] and [C II] emission is difficult because both lines show self-absorption features. By increasing the [O I] emission by a factor of 2, we obtain that the densities in the HV emission ranges increase up to $10^6$ cm$^{-3}$, which is more consistent with earlier findings (Schneider et al. 2003; Stock et al. 2015). If the emission in this velocity range could also arise from a shock (disk-envelope interaction and/or radiation and stellar wind hitting locally the dark lane) cannot yet be answered and will be addressed in an upcoming study.

The gas in the blue and red outflow velocity ranges (from -9 to -4 km s$^{-1}$ and 8 to 25 km s$^{-1}$) has a very similar emission distribution compared to the optical visible lobes. The [O I]
(and [CII]) emission distributions at the blue and red velocities indicate that the emission mostly arises from the back side of the southern lobe and the front side of the northern lobe. The outflow gas is entrained in the wind of S106 IR and ablated by radiation (and possibly shocks) from the cavity walls. PDR modelling gives densities typically of a few $10^4$ cm$^{-3}$ at a radiation field of $\chi$ a few $10^4$. The low value for the density indicates that the CO $16 \rightarrow 15$ line is most likely subthermally excited.

Molecular cloud clumps and possibly fragments of what once was the larger-scale circumbinary disk around S106 IR are seen in the close environment of the star. These clumps form PDRs on their surfaces (the most prominent one is the western clump) that emit in all lines at the cloud bulk velocity around $-3 \text{ km s}^{-1}$. The clumps are exposed to a radiation field $\chi$ of a $\sim 2 \times 10^4$ and the density is $\sim 2 \times 10^4 \text{ cm}^{-3}$.

S106 IR and its bipolar HII region are embedded in a larger molecular cloud that provides the gas reservoir for a possible accretion flow onto S106 IR. Mapping of the high-density tracer $H^{13}CO^+$ $1 \rightarrow 0$ revealed a velocity gradient across the dark lane that is consistent with either a flow onto S106 IR or gas streaming away, depending on the geometry of the region. Only interferometric observations can elucidate the nature of the dark lane. The flow scenario is more consistent with the fragmentation-induced starvation scenario of Peters et al. (2010a,b) than with the monolithic collapse model of McKee & Tan (2002). As yet unclear is whether shocks driven by the ionizing stellar wind that hits the accretion flow and the cavity walls can also cause the observed emission of the FIR lines. It is also not clear to what extent the [OI] $63 \text{ µm}$ line is self-absorbed. Higher line intensities will change the [CII]/[OI] line ratios and thus modify the outcome of the PDR models. Observations of the [OI] $145 \text{ µm}$ line, if this line is optically thin, may help to tackle this problem.

Summarizing, this study shows that the new detection of high-velocity emission in the [OI] line, and the identification of various velocity components in other FIR lines ([CII], high-J CO), which are only now possible with the (up)GREAT receiver on SOFIA, help to diagnose more precisely the physical properties of different gas phases in complex star-forming regions.

**Acknowledgements.** This work was supported by the Agence National de Recherche (ANR/France) and the Deutsche Forschungsgemeinschaft (DFG/Germany) through the project 'GENESIS' (ANR-16-CE92-0035-01/DFG1591/2-1). N.S. acknowledges support from the BMBF, Projekt Number 500R1714 (MOBS - MODellierung von Beobachtungsdaten SOFIA). This work is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. This work is based on observations carried out under project number 140-15 with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). N.S. acknowledges support from the Deutsche Forschungsgemeinschaft, DFG, through project number Os 177/2-1 and 177/2-2, and central funds of the DFG-priority program 1573 (ISM-SPP). This work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number SFB 956. We thank B. Rumph for carefully reading the manuscript.

## References

Adams, J. D., Herter, T. L., Hora, J. L., et al. 2015, ApJ, 814, 54

Agundez, M., Goicoechea, J. R., Cernicharo, J., et al. 2010, ApJ, 713, 662

Andree-Labsch, S., Ossenkopf-Okada, V. Röllig, M. 2017, A&A, 598, A2

Asplund, M., Grevesse, N., & Sauval, A. J. 2005, ASPC, 336, 25

Bally, J., Snell, R. L., & Predmore, R. 1983, ApJ, 272, 154

Balsara, D., Ward-Thompson, D., & Crutcher, R. M. 2001, MNRAS, 327, 715

Bate, M. R. 2012, MNRAS, 419, 3115

Beltran, M., Cesaroni, R., Neri, R., & Codella, C. 2011, A&A, 525, A151

Beuther, H., Schilke, P., & Menten K. M. 2002, ApJ, 566, 945

Bieging, J. H. 1984, ApJ, 286, 591

Bonnell, I. 2007, KITP Conference: Star Formation Then and Now, Santa Barbara, 25

Bontemps, S., André, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858

Bontemps, S., Motte, F., Csengeri, T., & Schneider, N. 2010, A&A, 524, A18

Cazaux, S., & Tielens, A. G. G. M. 2002, ApJ, 575, L29

Cazaux, S., & Tielens, A. G. G. M. 2004, ApJ, 604, 222

Cazaux, S., & Tielens, A. G. G. M. 2010, ApJ, 715, 698

Cohen, M., Bieging, J. H., & Welch, W. J. 1985, ApJ, 292, 249

Comerón, F., Schneider, N., Djupvik, A. A., & Schnugg, C. 2018, A&A, 615, A2

Csengeri, T., Bontemps, S., Wyrowski, F., et al. 2017, A&A, 600, L10

Draine, B. T. 1978, ApJS, 36, 595

Duarte-Cabral, A., Bontemps, S., Motte, F., et al. 2013, A&A, 550, A8

Eiroa, C., Elsässer, H., & Lahulla, J. F. 1979, A&A, 74, 89

Felli, M., Massi, M., Staude, H. J., et al. 1984, A&A, 135, 261

Furuya, R. S., Kitamura, Y., Saito, M., et al. 1999, ApJ, 525, 821

Gibb, A. G., & Hoare, M. 2007, MNRAS, 380, 246

Gibb, A. G., Hoare, M., Little L., & Wright M. 2003, MNRAS, 339, 1011

Gierens, K. M., Stutzki, J., & Winnewisser, G. 1992, A&A, 259, 271

Guan, X., Stutzki, J., Graf, U., et al. 2012, A&A, 542, L4

Guillet, V., Jones, A. P., & Pineau Des Forêts, G. 2009, A&A, 497, 145

Gusdorf, A., Cabrit, S., Flower, D. R., & Pineau Des Forêts, G. 2008a, A&A, 482, 809

Gusdorf, A., Pineau Des Forêts, G., Cabrit, S., & Flower, D. R. 2008b, A&A, 490, 695

Guzman, A. E., Garay, G., Rodriguez, L. F., et al. 2014, ApJ, 796, 117

Gvaramadze, V. V., & Menten, K. M. 2012, A&A, 541, A7

Gvaramadze, V. V., Kniazev, A. Y., & Fabrika, S. 2010, MNRAS, 405, 1047

Hartmann, L., Herczeg, G., & Calvet, N. 2016, Ann. Rev. Astron. Astrophys., 54, 135

Heiles, C. 1994, ApJ, 436, 720

Heyminck, S., Graf, U. U., Güsten, R., et al. 2012, A&A, 542, L1

Hipplelein, H., & Münch, G. 1981, A&A, 99, 248

Hoare, M., & Muxlow, T. 1996, Radio Emission from the Stars and the Sun, eds. A. R. Taylor, & J. M. Paredes (ASP), 93, 47

Hodapp, K., & Rayner, J. 1991, AJ, 102, 1108

Hodapp, K., & Schneider, N. 2008, Handbook of Star Forming Regions, Vol. I, ed. B. Reipurth (The Northern Sky ASP Monograph Publications, Vol. 4), 90

Hollenbach, D. 1985, Icarus, 61, 36

Hollenbach, D., & McKee, C. F. 1989, ApJ, 342, 306

Hollenbach, D., Kaufman, M. J., Neufeld, D., et al. 2012, ApJ, 754, 105

Hunter, T. R., Brogan, C. L., Indebetouw, R., et al. 2008, ApJ, 680, 1271

Kirk, H., Myers, P. C., Bourke, T. L., et al. 2013, ApJ, 766, 115

Kristensen, L. E., van Dishoeck, E. F., Mottram, J. C., et al. 2017, A&A, 605, A93

Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20

Kuruwita, R., Federrath, C., & Ireland, M. 2017, MNRAS, 470, 1626

Langer, W. D., & Penzias, A. A. 1990, ApJ, 357, 477

Larson, R. B. 1981, MNRAS, 194, 809

Leurini, S., Wyrowski, F., Wiesemayer, H., et al. 2015, A&A, 584, A70

Liseau, R., Justtanot, K., & Tielens, A. G. G. M. 2006, A&A, 446, 561

Lesaffre, P., Pineau des Forêts, G., Godard, B., et al. 2013, A&A, 550, A106

Little, L. T., Kelly, M. L., Habing, R., et al. 1995, MNRAS, 277, 307

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Maud, L. T., Moore, T. J. T., Lumsden, S. L., et al. 2015, MNRAS, 435, 645

McElroy, D., Walsh, C., & Markwick, A. J. 2013, A&A, 550, A36

McKee, C. F., & Tan, J. C. 2002, Nature, 416, 59

Minier, V., Tremblin, P., Hill, T., et al. 2013, A&A, 550, A50

Motte, F., Bontemps, S., Schilke, P., et al. 2007, A&A, 476, 1243

Nagy, Z., Van der Tak, F. F. S., & Ossenkopf, V. 2013, A&A, 550, A96

Nisini, B., Santangelo, G., Giannini, T., et al. 2015, ApJ, 801, 121

Noel, B., Joblin, C., Maillard, J. P., & Paumard, T. 2005, A&A, 436, 569

Oasa, Y., Tamura, M., Nakajima, Y., et al. 2006, AJ, 131, 1608

Peretto, N., Fuller, G. A., Duarte-Cabral, A., et al. 2013, A&A, 555, A112

Peters, T., Klessen, R. S., Mac Low, M.-M., et al. 2010a, ApJ, 725, 134

Peters, T., Banerjee, R., Klessen, R. S., et al. 2010b, ApJ, 711, 1017

Peters, T., Banerjee, R., Klessen, R. S., et al. 2011, ApJ, 729, 72

Polehampton, E. T., Baluteau, J.-P., & Swinyard, B. M. 2005, A&A, 437, 957

Rayner, T., Griffin, M., Schneider, N., et al. 2017, A&A, 607, A22

Richer, J., Padman, R., Ward-Thompson, D., et al. 1993, MNRAS, 262, 839

Riera, A., Mampaso, A., Vilchez, J. M., & Phillips, J. P. 1989, A&A, 210, 351

Röllig, M. 2011, A&A, 530, A9

Röllig, M., & Ossenkopf, V. 2013, A&A, 550, A56

Röllig, M., Ossenkopf, V., Jeyakumar, S., et al. 2006, A&A, 451, 917
Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187

Röllig, M., Szczerba, R., Ossenkopf, V., & Glück, C. 2013, A&A, 549, A85

Sanchez-Monge, A., Cesaroni, R., Beltran, M. T., et al. 2013, A&A, 552, L10

Schilke, P., Walmsley, C. M., Pineau des Forets, G., & Flower, D. R. 1997, A&A, 321, 293

Schneider, N., Simon, R., Kramer, C., Stutzki, J., & Bontemps, S. 2002, A&A, 384, 225

Schneider, N., Simon, R., Kramer, C., et al. 2003, A&A, 406, 915

Schneider, N., Simon, R., Bontemps, S., et al. 2007, A&A, 474, 873

Schneider, N., Csengeri, T., Bontemps, S., et al. 2010, A&A, 520, A49

Schneider, N., Bontemps, S., Motte, F., et al. 2016, A&A, 591, A40

Solf, J., & Carsenty, U. 1982, A&A, 116, 54

Simon, M., & Fischer, J. 1982, BAAS, 14, 925

Simon, R. Schneider, N., Stutzki, J., et al. 2012, A&A, 542, L12

Simon-Diaz, S., & Stasinska, G. 2011, A&A, 526, A48

Stock, D., Wolfire, M., Peeters, E., et al. 2015, A&A, 579, A67

Störzer, H., & Hollenbach, D. 1998, ApJ, 495, 853

Störzer, H., Stutzki, J., & Sternberg, A. 1996, A&A, 310, 592

Stutzki, J., Ungerechts, H., & Winnewisser, G. 1982, A&A, 111, 201

Szczerba, R., Omont, A., Volk, K., et al. 1997, A&A, 317, 859

Tan, J. C., Kong, S., Butler, M., et al. 2013, ApJ, 779, 96

Tan, J. C., Beltran, M. T., Caselli, P., et al. 2014, Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, T. Henning (Tucson: University of Arizona Press), 914, 149

Vallée, J. P., & Fiege, J. D. 2005, ApJ, 627, 263

van den Ancker, M., Tielens, A. G. G. M., & Wesselius, P. 2000, A&A, 111, 201

Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323

Weingartner, J. C., & Draine, B. T. 2001a, ApJ, 548, 296

Weingartner, J. C., & Draine, B. T. 2001b, ApJS, 134, 263

Wyrowski, F., Güsten, R., Menten, K. M., et al. 2012, A&A, 542, L15

Xu, Y., Li, J. J., Reid, M.J., et al. 2013, ApJ, 769, 15
## Appendix A: Chemical network

<table><caption>Table A.1. Chemical species included in the model.</caption><thead><tr><th></th><th>e<sup>-</sup></th><th>H</th><th>H<sub>2</sub></th><th>H<sup>+</sup></th><th>H<sub>2</sub><sup>+</sup></th><th>H<sub>3</sub><sup>+</sup></th><th>He<sup>+</sup></th><th>He</th><th>O<sup>+</sup></th></tr></thead><tbody><tr><td></td><td>O</td><td>C<sup>+</sup></td><td>C</td><td>13C<sup>+</sup></td><td>13C<sup>+</sup></td><td>OH<sup>+</sup></td><td>OH</td><td>O<sub>2</sub></td><td>CO<sup>+</sup></td></tr><tr><td></td><td>CO</td><td>CH<sup>+</sup></td><td>CH</td><td>13CO<sup>+</sup></td><td>13CO</td><td>13CH<sup>+</sup></td><td>13CH</td><td>HCO<sup>+</sup></td><td>H<sub>2</sub>O<sup>+</sup></td></tr><tr><td></td><td>H<sub>2</sub>O</td><td>H<sup>13</sup>CO<sup>+</sup></td><td>CH<sub>2</sub><sup>+</sup></td><td>13CH<sub>2</sub><sup>+</sup></td><td>H<sub>3</sub>O<sup>+</sup></td><td>C<sub>2</sub><sup>+</sup></td><td>C<sub>2</sub></td><td>C<sup>13</sup>C<sup>+</sup></td><td>C<sup>13</sup>C</td></tr><tr><td></td><td>C<sub>3</sub><sup>+</sup></td><td>C<sub>3</sub></td><td>C<sup>13</sup>CF<sup>+</sup></td><td>C<sup>13</sup>CH<sub>2</sub><sup>+</sup></td><td>O<sub>2</sub><sup>+</sup></td><td>HF<sup>+</sup></td><td>HF</td><td>H<sub>2</sub>F<sup>+</sup></td><td>F<sup>+</sup></td></tr><tr><td></td><td>F</td><td>CF<sup>+</sup></td><td>C<sup>13</sup>CF<sup>+</sup></td><td>CH<sub>2</sub></td><td>C<sub>2</sub>H<sup>+</sup></td><td>C<sub>2</sub>H</td><td>C<sup>13</sup>CH<sup>+</sup></td><td>C<sup>13</sup>CH</td><td>13CH<sub>2</sub></td></tr><tr><td></td><td>CH<sub>3</sub><sup>+</sup></td><td>CH<sub>3</sub></td><td>C<sub>3</sub>H<sup>+</sup></td><td>C<sub>3</sub>H</td><td>C<sub>2</sub>H<sub>2</sub><sup>+</sup></td><td>C<sub>2</sub>H<sub>2</sub></td><td>C<sup>13</sup>CH<sup>+</sup></td><td>C<sup>13</sup>CH</td><td>C<sup>13</sup>CH<sub>2</sub><sup>+</sup></td></tr><tr><td></td><td>C<sup>13</sup>CH<sub>2</sub></td><td>13CH<sub>3</sub><sup>+</sup></td><td>13CH<sub>3</sub></td><td>CH<sub>4</sub><sup>+</sup></td><td>CH<sub>4</sub></td><td>13CH<sub>4</sub><sup>+</sup></td><td>13CH<sub>4</sub></td><td>CH<sub>5</sub><sup>+</sup></td><td>CH<sub>2</sub>CCH</td></tr><tr><td></td><td>CH<sub>2</sub>C<sup>13</sup>CH</td><td>CH<sup>13</sup>CCH</td><td>C<sub>3</sub>H<sub>3</sub><sup>+</sup></td><td>C<sub>2</sub>H<sub>4</sub><sup>+</sup></td><td>C<sub>2</sub>H<sub>4</sub></td><td>C<sup>13</sup>CH<sub>3</sub><sup>+</sup></td><td>C<sup>13</sup>CH<sub>4</sub><sup>+</sup></td><td>C<sup>13</sup>CH<sub>4</sub></td><td>13CH<sub>5</sub><sup>+</sup></td></tr><tr><td></td><td>13CH<sub>2</sub>CCH</td><td>C<sub>2</sub>H<sub>5</sub><sup>+</sup></td><td>C<sub>2</sub>H<sub>5</sub></td><td>C<sup>15</sup>CH<sub>5</sub><sup>+</sup></td><td>C<sup>15</sup>CH<sub>5</sub></td><td>C<sup>13</sup>CH<sub>5</sub><sup>+</sup></td><td>HOC<sup>+</sup></td><td>HO<sup>13</sup>C<sup>+</sup></td><td>HCO<sub>2</sub><sup>+</sup></td></tr><tr><td></td><td>HCO</td><td>H<sub>3</sub>CO<sup>+</sup></td><td>H<sup>13</sup>CO<sup>+</sup></td><td>H<sub>2</sub>CO<sup>+</sup></td><td>H<sub>2</sub>CO</td><td>H<sup>13</sup>CO<sup>+</sup></td><td>H<sup>13</sup>CO</td><td>H<sup>13</sup>CO</td><td>CO<sub>2</sub><sup>+</sup></td></tr><tr><td></td><td>CO<sub>2</sub></td><td>CH<sub>3</sub>OH<sub>2</sub><sup>+</sup></td><td>CH<sub>3</sub>OH<sup>+</sup></td><td>CH<sub>3</sub>OH</td><td>13CO<sub>2</sub><sup>+</sup></td><td>13CO<sub>2</sub></td><td>13CH<sub>3</sub>OH<sup>+</sup></td><td>13CH<sub>3</sub>OH</td><td>HCOOCH<sub>3</sub></td></tr><tr><td></td><td>HCOO<sup>13</sup>CH<sub>3</sub></td><td>HC<sub>2</sub>O<sup>+</sup></td><td>HC<sup>13</sup>CO<sup>+</sup></td><td>H<sub>5</sub>C<sub>2</sub>O<sub>2</sub><sup>+</sup></td><td>H<sub>5</sub>C<sup>13</sup>CO<sub>2</sub><sup>+</sup></td><td>H<sup>13</sup>COOCH<sub>3</sub></td><td>COOCH<sub>4</sub><sup>+</sup></td><td>CH<sub>3</sub>CHO<sup>+</sup></td><td>CH<sub>3</sub>CHO</td></tr><tr><td></td><td>CH<sup>13</sup>CHO<sup>+</sup></td><td>CH<sup>13</sup>CHO</td><td>13CH<sub>3</sub>CHO<sup>+</sup></td><td>13CH<sub>3</sub>CHO</td><td>OCN<sup>+</sup></td><td>OCN</td><td>O<sup>13</sup>CN<sup>+</sup></td><td>O<sup>13</sup>CN</td><td>NO<sub>2</sub><sup>+</sup></td></tr><tr><td></td><td>NO<sub>2</sub></td><td>NO<sup>+</sup></td><td>NO</td><td>NH<sub>4</sub><sup>+</sup></td><td>NH<sub>3</sub><sup>+</sup></td><td>NH<sub>3</sub></td><td>NH<sub>2</sub><sup>+</sup></td><td>NH<sub>2</sub></td><td>NH<sup>+</sup></td></tr><tr><td></td><td>NH</td><td>N<sup>+</sup></td><td>N</td><td>HNO<sup>+</sup></td><td>HNO</td><td>HNCO<sup>+</sup></td><td>HNC</td><td>HN<sup>13</sup>CO<sup>+</sup></td><td>HN<sup>13</sup>C</td></tr><tr><td></td><td>HCN<sup>+</sup></td><td>HCN</td><td>H<sub>2</sub>NC<sup>+</sup></td><td>H<sub>2</sub>N<sup>13</sup>C<sup>+</sup></td><td>H<sub>2</sub>CN</td><td>H<sup>13</sup>CN</td><td>H<sup>13</sup>CN<sup>+</sup></td><td>H<sup>13</sup>CN</td><td>CN<sup>+</sup></td></tr><tr><td></td><td>CN</td><td>13CN<sup>+</sup></td><td>13CN</td><td>N<sub>2</sub>O<sup>+</sup></td><td>N<sub>2</sub>O</td><td>N<sub>2</sub>H<sup>+</sup></td><td>N<sub>2</sub><sup>+</sup></td><td>N<sub>2</sub></td><td>SO<sub>2</sub><sup>+</sup></td></tr><tr><td></td><td>SO<sub>2</sub></td><td>SO<sup>+</sup></td><td>SO</td><td>S<sup>+</sup></td><td>S</td><td>OCS<sup>+</sup></td><td>OCS</td><td>O<sup>13</sup>CS<sup>+</sup></td><td>O<sup>13</sup>CS</td></tr><tr><td></td><td>NS<sup>+</sup></td><td>NS</td><td>HS<sup>+</sup></td><td>HS</td><td>HCS<sup>+</sup></td><td>HCS</td><td>H<sub>2</sub>S<sup>+</sup></td><td>H<sub>2</sub>S</td><td>H<sub>2</sub>CS<sup>+</sup></td></tr><tr><td></td><td>H<sub>2</sub>CS</td><td>H<sup>13</sup>CS<sup>+</sup></td><td>H<sup>13</sup>CS</td><td>H<sup>13</sup>CS<sup>+</sup></td><td>H<sup>13</sup>CS</td><td>CS<sup>+</sup></td><td>CS</td><td>C<sup>13</sup>CS<sup>+</sup></td><td>C<sup>13</sup>CS</td></tr><tr><td></td><td>13CS<sup>+</sup></td><td>13CS</td><td>S<sub>2</sub><sup>+</sup></td><td>S<sub>2</sub></td><td>Si</td><td>Si<sup>+</sup></td><td>SiOH<sup>+</sup></td><td>SiH<sub>5</sub><sup>+</sup></td><td>SiH<sub>4</sub><sup>+</sup></td></tr><tr><td></td><td>SiH<sub>4</sub></td><td>SiH<sub>3</sub><sup>+</sup></td><td>SiH<sub>3</sub></td><td>SiH<sub>2</sub><sup>+</sup></td><td>SiH<sub>2</sub></td><td>SiH<sup>+</sup></td><td>SiH</td><td>SiC<sup>+</sup></td><td>SiC</td></tr><tr><td></td><td>Si<sup>13</sup>C<sup>+</sup></td><td>Si<sup>13</sup>C</td><td>HCSi<sup>+</sup></td><td>HCSi</td><td>H<sup>13</sup>CSi<sup>+</sup></td><td>H<sup>13</sup>CSi</td><td><sup>18</sup>O<sup>+</sup></td><td><sup>18</sup>O</td><td>O<sup>18</sup>O</td></tr><tr><td></td><td>H<sup>18</sup>O<sup>+</sup></td><td>H<sup>18</sup>O<sup>+</sup></td><td>H<sup>18</sup>O<sup>+</sup></td><td>H<sup>18</sup>O<sup>+</sup></td><td><sup>18</sup>OH<sup>+</sup></td><td>HC<sup>18</sup>O<sup>+</sup></td><td>H<sup>13</sup>C<sup>18</sup>O<sup>+</sup></td><td>C<sup>18</sup>O<sup>+</sup></td><td>C<sup>18</sup>O</td></tr><tr><td></td><td>13C<sup>18</sup>O<sup>+</sup></td><td>13C<sup>18</sup>O</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

## Appendix B: Overlays of different tracers in individual velocity ranges

<figure><img src="image_17.png" /><figcaption>Fig. B.1. High-velocity blue emission (–30 to –9 km s⁻¹): contours of [C II] (30–100 K km s⁻¹), and ¹²CO 2→1 (3.1–44.4 K km s⁻¹) emission overlaid on [O I] emission in colour scale (wedge in K km s⁻¹).</figcaption></figure>

Figures B.1–B.5 show overlays between [O I] emission in colour scale and other tracers ([C II], CO 16→15, ¹²CO 2→1, and ¹³CO 2→1) in the five main velocity ranges.
<figure><img src="image_18.png" /><figcaption>**Fig. B.2.** Outflow blue emission (–9 to –4 km s⁻¹): contours of [C II] (59–192 K km s⁻¹), ¹²CO 2→1 (3.1–44.4 K km s⁻¹), ¹²CO 16→15 (13.1–191.1 K km s⁻¹), and ¹³CO 2→1 (2.0–28.9 K km s⁻¹) emission overlaid on [O I] emission in colour scale (wedge in K km s⁻¹).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. B.3. Molecular cloud bulk emission (–4 to 0.5 km s⁻¹): contours of [C II] (72.5–235.5 K km s⁻¹), ¹²CO 2→1 (24.7–357.8 K km s⁻¹), ¹²CO 16→15 (9.5–66.9 K km s⁻¹), and ¹³CO 2→1 (5.8 to 84.5 K km s⁻¹) emission overlaid on [O I] emission in colour scale (wedge in K km s⁻¹).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. B.4.** Outflow red emission (0.5–8 km s⁻¹): contours of [C II] (58.5–190.2 K km s⁻¹), ¹²CO 2→1 (3.2–46.6 K km s⁻¹), and ¹²CO 16→15 (15.0–48.9 K km s⁻¹) emission overlaid on [O I] emission in colour scale (wedge in K km s⁻¹).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. B.5.** High-velocity red emission (8–25 km s⁻¹): contours of [C II] (34.4–111.7 K km s⁻¹) and ¹²CO 2→1 (3.1–44.4 K km s⁻¹) emission overlaid on [O I] emission in colour scale (wedge in K km s⁻¹).</figcaption></figure>