© The Authors 2024

1 Introduction

Bright-rimmed clouds (BRCs) are located near the edges of evolved H II regions and are exposed to the UV radiation from the exciting stars, creating a photon-dominated region (PDR) at their surface (the bright rim). Their evolution and morphology, which sometimes includes young stellar objects (YSOs), are often discussed in the context of triggered star formation from a radiation-driven implosion (Bertoldi 1989; Lefloch & Lazareff 1994; Miao et al. 2009). Early analytical studies (Oort & Spitzer 1955) and more sophisticated 3D radiation-magnetohydrodynamic simulations (Henney et al. 2009) predict a photoevaporation flow from the illuminated surface and a rocket effect to the cloud as its back reaction. Sugitani et al. (1991) classified BRCs into three types based on the morphology (rim curvature), which is explained by the combination of the initial self-gravitational state and the evolutionary sequence (Miao et al. 2009). The term “BRC” has been used in the literature to refer to a surface illuminated by exciting sources, while the term “pillars” (a column-like shape) and “globules” (isolated and having a head-tail structure) are used to specify their morphology (Schneider et al. 2016).

IC 1396 is an H II region with a diameter of a few degrees on the sky (Fig. 1; 1 degree corresponds to 16 pc at IC 1396’s distance of 925 pc; Pelayo-Baldárrago et al. 2023), located at the edge of the Cepheus Bubble. The whole region is predominantly excited by the O6.5V star in the binary system HD 206267 (Abt 1986; Sota et al. 2011), which is a member of the cluster Tr 37. An optical image of this region shows a number of cometary silhouettes, some of them with bright rims facing the exciting source. Most of them are bright in the mid-infrared (MIR) and far-infrared (FIR; Fig. 1). This rich selection of BRCs in an apparently simple geometry with a single dominant exciting star makes IC1396 a textbook example of an evolved H II region and ideal for studying the evolution of BRCs (e.g., Sugitani et al. 1991; Patel et al. 1995; Weikard et al. 1996; Sicilia-Aguilar et al. 2019).

The [C II] 158 µm and [O I] 63 µm FIR fine-structure lines are the dominant cooling lines in PDRs and thus are excellent tracers of the physical conditions of BRCs and, if observations are velocity-resolved, their kinematics. Recent studies demonstrate that large, velocity-resolved [C II] 158 µm maps are a powerful tool for investigating the radiative and mechanical feedback as well as the physical properties of the illuminated surface, together with other PDR and molecular gas tracers such as [O I] and CO in nearby H II regions (Schneider et al. 2012; Mookerjea et al. 2019; Pabst et al. 2019; Luisi et al. 2021; Tiwari et al. 2021; Schneider et al. 2021; Kabanovic et al. 2022).

So far, the [C II] emission from IC 1396A has been analyzed using only an under-sampled map in Okada et al. (2012). In the current study, we present 4–10′ scale fully sampled [C II] maps around four BRCs in IC 1396 and a deep [O I] 63 µm map of IC 1396A, which allowed us to investigate the dynamics and physical properties of the regions. We focused on the detailed structure and dynamics of individual BRCs. The excitation and column densities of the carbon-bearing species and the modeling using PDR models will be presented in a follow-up paper.

The organization of this paper is as follows. In Sect. 2 we give an overview of IC 1396 based on previous studies, and in Sect. 3 we present the data used in this study. Section 4 presents the results of our analysis. After providing an overview of IC 1396 as a whole (Sect. 4.1) and explaining the methods for calculating integrated line intensities (Sect. 4.2), we reveal the details of the spatial distributions and dynamics of the four observed sub-regions (Sects. 4.3.1 to 4.3.4). The [¹³C II] analysis of IC 1396A is presented in Sect. 4.4, followed by a discussion of the heating efficiency and the [C II]–mid-infrared correlation in Sects. 4.5 and 4.6. Signs of photoevaporation flows, 3D kinematic structures, and a comparison of different BRCs are discussed in Sect. 5. Our main results are summarized in Sect. 6.

Fig. 1 Overview of IC 1396. In this three-color composite image, blue is DSS-гed, green is the 9 µm map from the MIR all-sky survey by the Infrared Camera (IRC; Onaka et al. 2007; Ishihara et al. 2010) on board AKARI (Murakami et al. 2007), and red is the WIDE-L (140 µm) map from the FIR all-sky survey by the Far-Infrared Surveyor (FIS; Kawada et al. 2007; Doi et al. 2015) on board AKARI. The maps are not convolved to the common spatial resolution, and the coarsest spatial resolution is 88″ (Takita et al. 2015) for FIS (red), which is about a third of the size of the box around at IC 1396A. The fields observed with upGREAT in the [C II] line emission are overlaid as white boxes. The regions are labeled following the notation of Weikard et al. (1996). The position of the exciting source, HD 206267, is marked.

2 IC1396

The Cepheus Bubble is thought to have been created by a first generation of stars formed 13–18 Myr ago, followed by the formation of the Cep OB2 association in the compressed expanding shell (Patel et al. 1998). HD 206267 is estimated from its stellar properties to have been born about 4 Myr ago (Patel et al. 1995, and references therein), creating an H II region. The BRCs are found at the interface between the H II region and the associated molecular cloud. Patel et al. (1995) proposed an expanding ring model based on the different CO velocities in the surrounding BRCs. They show that the momentum of the BRCs cannot be explained by the stellar wind, and model the dynamical evolution with the “rocket effect,” suggesting a dynamical age of IC 1396 of 2–3 Myr. Using Gaia data, Pelayo-Baldárrago et al. (2023) found that the stellar component in IC 1396 is composed of several substructures that are clearly distinct in proper motion and partly in age, suggesting a complex star formation history with multi-episode formation processes.

Figure 1 shows a three-color composite image of optical, MIR, and FIR emission. As already shown in previous studies (Weikard et al. 1996; Barentsen et al. 2011), BRCs often appear as a silhouette with a bright rim in optical images. Weikard et al. (1996) made an extensive ¹²CO/¹³CO(1−0) survey covering an area of more than 6 deg² in IC 1396, together with focused maps around a number of BRCs in the J = 2−1 and 3−2 transitions. They show the presence of warmer molecular gas in the BRCs and discuss the detailed kinematics of individual clouds. Whether or not the detected clouds are associated with IC 1396 has been investigated based on several criteria such as excitation conditions and kinematics. In general, clouds with central velocities of −8 km s⁻¹ to +7 km s⁻¹ associate with IC 1396. There are a number of clouds that are not associated with IC 1396 despite having similar central velocities. For example, there is a prominent cometary optical dark cloud between HD 206267 and IC 1396E (Fig. 1) that was identified in CO with a velocity of −0.7 km s⁻¹ (Weikard et al. 1996, B161, part of the broad dark lane Khav 161). Despite its appearance, they found no indication that it is associated with IC 1396. This area is quite dark in the MIR image and at shorter wavelength in the FIR image (Fig. 1). Some structure, including B161, starts to become visible between 140 µm and 160 µm.

Carlsten & Hartigan (2018) explained the observed positional offset of H₂ and Brγ in some of the BRCs of IC 1396 with photoevaporative flow models. Sicilia-Aguilar et al. (2019) also suggested that the velocity of CN is indicative of photoevaporative flows. Soam et al. (2018) studied the magnetic field and the shape of BRCs in IC 1396 and confirmed the theoretical prediction (Henney et al. 2009) that BRCs with small angles between the direction of the magnetic field and the ionizing radiation have broad heads, while BRCs with large angles of these have a more flattened and elongated shape.

Of the four BRCs that we study in this paper, IC 1396A and IC 1396E have been particularly well studied, and we summarize them in the following two subsections. IC 1396B is connected to IC 1396A, and is located further away from the exciting source (Fig. 1). IC 1396D has the largest projected distance from the exciting source and is weaker in the infrared among the four BRC samples. It contains a cometary-shaped BRC with an embedded Infrared Astronomical Satellite (IRAS) source identified by Sugitani et al. (1991), SFO34.

2.1 IC 1396A

IC 1396A is the prominent BRC closest in projection to the exciting source (Fig. 1). A cavity within the BRC is visible in ¹²CO and ¹³CO maps (Wootten et al. 1983; Nakano et al. 1989, see also Fig. 5 in Sect. 4.3.1), and is associated with the YSOs LkHα 349a and c. The cavity is also clearly visible in the MIR (Reach et al. 2004, 2009) and FIR (Sicilia-Aguilar et al. 2014) images. Reach et al. (2009) showed that it is plausible that the stellar wind from LkHα 349a, more massive than LkHα 349c, created a bubble inside the BRC. LkHα 349a is likely a very young Herbig Ae/Be star and will be a late B- to early A-type star on the main sequence (Reach et al. 2009). It shows P Cygni profiles in Hα, indicating ejection of material at velocities above 300 km s⁻¹ (Hernández et al. 2004). The rim of the cavity is bright in shock-tracing optical [S II] emission (Sicilia-Aguilar et al. 2013), reflecting interaction between the YSOs formed inside the cavity and the surrounding cloud. In addition to this most prominent cavity, the Spitɀer image shows that IC 1396A hosts a number of cavities and that all of the major cavities contain YSOs, suggesting that the BRC is being reshaped from the inside via outflows from protostars (Reach et al. 2009). IC 1396A must be located in front of the ionizing region based on its appearance as a dark cloud with a bright rim in optical wavelengths (Barentsen et al. 2011), a sharp edge of Hα (Barentsen et al. 2011) and [S II] (Sicilia-Aguilar et al. 2013) emission, and its bluest velocity among the clouds in IC 1396 (Weikard et al. 1996). Sicilia-Aguilar et al. (2014) found a Class 0 object (IC 1396A-PACS-1) with Herschel/Photodetector Array Camera and Spectrometer (PACS) observations. It is located inside the BRC, ~0.1 pc from the bright rim defined by a sharp boundary in the PACS map, surrounded by a cold region, and is consistent with its formation being triggered via radiation-driven implosion (Sicilia-Aguilar et al. 2014). Combining the Institut de radioastronomie millimétrique (IRAM) 30m observations and Gaia data, Sicilia-Aguilar et al. (2019) confirm that IC 1396A has undergone several episodes of star formation. From the Gaia proper motion of the stars and the radial velocity of the gas, they derived a total (3D) velocity of IC 1396A of ~8 km s⁻¹ away from Tr 37. This relative velocity is too low compared to the expected velocity due to rocket acceleration if IC 1396A had been much closer to HD206267 in the past. Moriarty-Schieven et al. (1996) found an H I tail behind IC 1396A (opposite to the exciting stars) extending across IC 1396B and further west. The blueshifted H I velocity relative to CO at the head of the BRC and the relatively constant velocity of the tail is consistent with a scenario where H I traces ablated material from the BRCs, drifting into the shadow of the BRC so that it is not further accelerated. Okada et al. (2012) presented the first velocity-resolved [C II] map in IC 1396A, revealing the prominent [C II] wings at the northern and southern edges of the rim. They performed a simple KOSMA-τ PDR model fit to the [C II], [C I], and CO(4−3) data and suggest that the UV radiation is shielded at the rim.

2.2 IC 1396E

Another well-studied BRC is IC 1396E (also called IC 1396N), which hosts prominent protostellar outflows, located in the northern part of IC 1396 (Fig. 1). Sugitani et al. (1989) mapped the central outflow in CO(2−1) and show that it is associated with IRAS 21391+5802. The density profile derived from CS observations is steep and centered at this IRAS source (Serabyn et al. 1993), well away from the ionization front, suggesting a preexisting density enhancement rather than compression by external shocks. IRAS 21391+5802 was resolved into three sources in the millimeter continuum (BIMA1, 2, and 3; Beltrán et al. 2002, see also Fig. 15 in Sect. 4.3.4), and the brightest, BIMA2, which drives the outflow, was further resolved into three sources in millimeter interferometric observations (Neri et al. 2007), forming an intermediate-mass protocluster. Fuente et al. (2009) used molecular line observations and suggested that the brightest identified source in BIMA2, IRAM 2A, is a low-mass or Herbig Ae star, or a Class 0/I transition object that has already formed a small PDR. Codella et al. (2001) found another outflow in the northern part of IC 1396E, and Beltrán et al. (2012) showed that it collides with yet another outflow. The outflows have a clumpy morphology, suggesting that they are the result of episodic mass-loss events (Codella et al. 2001; Beltrán et al. 2012). H₂ 2.12 µm knots are observed in collimated stellar jets, excited by the shock produced by the outflows (Nisini et al. 2001; Beltrán et al. 2009). The interaction of the outflows with the ambient gas is also seen in the distorted shape of the western lobe of the central outflow. Codella et al. (2001) suggested a significant contribution of swept-up material to the central outflow mass based on the low velocity component (±3 km s⁻¹ with respect to the ambient velocity). CO emission up to J = 25 was detected with Infrared Space Observatory (ISO)/Long Wavelength Spectrometer (LWS) (Saraceno et al. 1996), likely tracing shock heated gas (Molinari et al. 1998; Nisini et al. 1998). Patel et al. (1995) revealed the dynamics of the CO gas using¹³ CO(1−0); there is a pair of symmetrical extended features from the cometary head to the northeast and northwest, likely moving away from the cometary head. Their velocities differ by 1.7 km s⁻¹ (also visible in the ¹³CO(3−2) channel maps in Fig. E.5). They proposed that the velocity difference in the two tails is due to either the different inclinations or the different densities, the former resulting in different line-of-sight velocities and the latter causing different shock velocities.

3 Observations and data reduction

3.1 [O I] 63µm, 145µm, and [C II] 158µm observations with up GREAT on board SOFIA in IC 1396A

We observed the [O I] lines at 4744.77749 GHz (63 µm) and 2060.06886 GHz (145 µm) and the [C II] line at 1900.53690 GHz (158 µm) with the German REceiver for Astronomy at Terahertz Frequencies¹ (upGREAT; Heyminck et al. 2012; Risacher et al. 2018) on board the Stratospheric Observatory for Infrared Astronomy (SOFIA; Young et al. 2012) in IC 1396A, as a guaranteed time project. The High Frequency Array (HFA) with 7 pixels in a hexagonal configuration was tuned to [O I] 63 µm and all observations were made in single phase chopped on-the-fly (OTF) mode with 3″ step size (see Table 1 for the chop parameters). We observed the [C II] emission in parallel, either with the 7 pixels of the Low Frequency Array (LFA) or, earlier on, with the single pixel L2 channel. On some flights, we split the two polarizations (H and V) in the LFA and observed [O I] 145 µm in parallel (see Table 1 for receiver combinations). We used several generations of the Fast Fourier Transform Spectrometer (FTS) as backends, which were continuously upgraded over several years of observations in this project.

The data were calibrated along the following steps. First, we performed a frequency correction on the HFA data. This is a standard procedure for the HFA, because the absolute frequency stability of the quantum cascade laser local oscillator is limited by the temperature and current stability of the laser (Risacher et al. 2018). We corrected the frequency scale scan by scan using the known narrow telluric [O I] line as a reference. We then corrected for a gain drift by scaling the raw counts of individual observations to those immediately after each load measurement for individual pixels. The obtained data were calibrated using the standard GREAT pipeline (Guan et al. 2012), which converts the observed counts to the main beam temperature scale (T_mb) and corrects for atmospheric attenuation using an atmospheric radiative transfer model. We then resampled the spectra to a resolution of 0.5 km s⁻¹. To improve the baseline quality, we extracted principal components from the OFF spectra and subtracted these components from the ON spectra of [C II] and [O I] 63 µm (see our Appendix A as well as Tiwari et al. 2021 and Buchbender et al. in prep. for details of the methods and evaluation of the results). We then convolved the [C II] and [O I] 63 µm maps to a spatial resolution of 16″ or 25″. We use 16″ resolution to discuss the velocity profiles and spatial structure in channel maps, integrated intensity maps, velocity profiles in selected positions and position-velocity (p–v) diagrams. For the analysis of the line widths and the [¹³C II] line, and the derivation of the heating efficiencies, we used 25″ maps to have abetter signal-to-noise ratio (S/N). The median rms noise of the spectra at a given pixel on the 16″ resolution map is 0.4 K and 0.18 K for [C II] and [O I], respectively.

The [O I] 145 µm is located on a shoulder of a significant atmospheric ozone feature in the signal sideband (Fig. D.1). The line is detected at only a few positions on the 16″ map. Because the atmospheric transmission is smooth in the velocity range of interest, we consider the velocity profile of [O I] 145 µm as reliable as that of the other lines (Appendix D).

Since we obtained a fully sampled map for [O I] (beam size of 6.3″), the parallel [C II] map (beam size of 14.1″) has a S/N ten times higher than the data shown in Okada et al. (2012), observed with a single pixel. The [C II] map obtained in this study thus reveals much more detailed structures.

3.2 [C II] 158µm observations with (up)GREAT on board SOFIA in IC 1396E, D, and B

[C II] mapping observations were also made in three other regions: IC 1396E, D, and B, as part of the guaranteed time observations. IC 1396E was observed with the single pixel L2 channel in May 2014 and January 2015, and with the 7 pixel upGREAT LFA in May 2015 during the upGREAT commissioning flights. We also observed IC 1396B and D in May 2015. IC 1396D observations were extended in December 2015. All observations were made in total-power OTF mode with a 5–6″ step size. The OFF positions are (RA, Dec.)(J2000) = (21:34:13.0, 57:23:18.0), (21:33:39.4, 57:58:20.0), and (21:40:39.0, 58:13:28.0) for IC 1396B, IC 1396D, and IC 1396E, respectively.

The data were calibrated using the standard GREAT pipeline as described above. We subtracted linear baselines and spectrally resampled to a resolution of 0.5 km s⁻¹. We did not apply the principal component analysis (PCA) to these data because the analysis with the chopped data in IC 1396A (Sect. 3.1) shows that the baseline improvement by using PCA is not significant in the LFA data and the total power observations do not provide a sufficient number of OFF measurements to identify the principal components. We convolved the [C II] map to 16″ or 25″ resolution as for IC 1396A (except that we use 25″ maps for p–v diagrams in IC 1396B and IC 1396D for better S/N). The median rms noise on the 16” resolution map is 1.2 K, 0.74 K, and 0.77 K for IC1396 B, IC1396 D, and IC1396 E, respectively.

3.3 Complementary observations: CO, ¹³CO, and C¹⁸O from JCMT archive data

We obtained calibrated CO spectral maps from the Heterodyne Array Receiver Program (HARP; Buckle et al. 2009) on the James Clerk Maxwell Telescope (JCMT²) from their data archive³. Table 2 summarizes the available isotopes and transitions. The original beam size is 14″ for CO(3−2) and 19.7″ for CO(2−1). In this paper, we focus on CO(3−2) and its isotopologues.

We subtracted a zeroth-order baseline and spectrally resam-pled to a resolution of 0.5 km s⁻¹. We then convolved the maps to the same spatial resolution and grid as the upGREAT [C II] data (16″ or 25″ resolution). To convert the intensity to the main beam scale, we utilized main beam and forward efficiencies of 0.63 and 0.88, respectively, for CO(3−2) and its isotopologues (following JCMT’s heterodyne data reduction cookbook ver1.0).

3.4 Complementary observations: Far-infrared continuum

To derive the global distribution of the dust temperature and column density (Sect. 4.1), we obtained FIR maps from the FIR all-sky survey by the Far-Infrared Surveyor (FIS; Kawada et al. 2007; Doi et al. 2015) on board AKARI. They have four bands, N60, WIDE-S, WIDE-L, and N160, centered at 65 µm, 90 µm, 140 µm, and 160 µm. While the beam sizes along the in-scan and cross-scan directions are not the same (Takita et al. 2015), we use a simplified circular Gaussian with the mean full width at half maximum to convolve the N60 and WIDE-S maps (63″ and 78″, respectively) to the resolution of WIDE-L and N160 (88″). We also applied the color correction tabulated in Shirahata et al. (2009) to take into account the discrepancy between the integrated flux in an imaging band and the intensity at the central wavelength of the band. We used the correction factor corresponding to a gray-body with β = 2 and T = 30 K.

In the analysis of the heating efficiencies in IC 1396A (Sect. 4.5), we used higher spatial resolution data from Herschel⁴ and the JCMT⁵. We downloaded the following data from their data archives; PACS 70 µm and 160 µm data presented in Sicilia-Aguilar et al. (2014)⁶, and SCUBA-2 850 µm data⁷ (Holland et al. 2013; Chapin et al. 2013; Dempsey et al. 2013). For the SCUBA-2 data, we selected maps observed on four different nights and combined them with equal weight. We convolve the PACS and SCUBA-2 maps to 25″ resolution with the same spatial grid as the line emission data, with which we performed the heating efficiency analysis. In the SCUBA-2 850 µm map, the cavity around LkHα 349 appears to have negative intensities. The most negative intensity in the convolved map is −26 MJy sr⁻¹, which was added as an offset to the whole map and is considered as a systematic uncertainty in the absolute flux.

4 Results

4.1 Global picture from dust emission

We fit the spectral energy distribution (SED) from the four FIS bands with a single temperature modified gray-body with β = 2.

This corresponds to a simplified dust model in which the wavelength dependence of the absorption coefficient follows λ^−β in the infrared (for the dust emission) and is constant in the optical (for the stellar light absorption), β = 2 falls in the range of the dust models of Jones et al. (2013), Draine & Li (2007) or Hensley & Draine (2023, see Fig. 16 therein). Dense molecular clouds show a shallower spectral index β. We fit the dust temperature (T) and the optical depth at a reference wavelength (100 µm, τ₁₀₀) in the equation τ₁₀₀(λ/100 µm)^−βB_v(T). Due to an insufficient number of data points to differentiate between different dust models, here we focus on the relative spatial distribution of the dust temperature and optical depth. Figure 2 shows the resulting temperature and dust optical depth at 100 µm for the entire IC 1396 region. The area with T_dust ≥ 20 K has a more or less round distribution centered at HD 206267. It does not have a clear gradient toward the center (cross in Fig. 2) and there is no significant structures except for a colder dust lane running from the southeast to the north. This cold dust lane includes Khav 161, and probably traces foreground gas that is not heated by HD 206267. We calculated the expected dust temperature as a function of distance from HD 206267, assuming a luminosity of 10^5.2 L_⊙ for an O6.5V star (Martins et al. 2005) and the dust model described above. The derived dust temperature 16 pc away from the source (1 deg at the distance of 925 pc) is 25 K, which is shown as a white circle in Fig. 2. This is a few kelvins higher than the temperature obtained from the SED fit, corresponding to a factor of 1.9 in energy (∝T⁶), but considering that the fit model is very simplified and the theoretical estimate ignores energy loss within the H II regions and considers only the projected distance, we conclude that they are in good agreement and confirm that HD 206267 is the dominant heating source. In the dust optical depth map, the BRCs are clearly peaking out, but otherwise the dynamic range of the dust column density is not high. The small variation of the dust temperature and the dust optical depth, except for the BRCs, and the fact that the region is well visible in Hα (Barentsen et al. 2011) give the picture of a spherical H II region breaking out only in the direction toward the observer.

Fig. 2 Dust temperature (left) and the continuum optical depth at 100 µm (right) fitted with a single temperature model to the AKARI/FIS data (see the main text). The white cross marks the position of HD 206267, and the white circle in the left panel indicates a 1 deg radius, where the dust temperature is 25 K from a purely geometric estimate (see the main text).

4.2 Integrated line intensities

In IC 1396A and IC 1396E, the integrated line intensities were calculated by integration over a fixed velocity range in each region (−13 to −4 km s⁻¹ and −5 to 8 km s⁻¹, respectively). On the other hand in IC 1396B and IC 1396D, the velocity structure is more complex and different velocity components coexist in the observed regions, which requires a broader velocity range to integrate (−15 to 10 km s⁻¹ and −12 to 10 km s⁻¹, respectively) and leads to significantly higher noise (see the lowest panel of Fig. B.1). For these two subregions, we derived the integrated intensities using dynamic integration windows based on a 3D dendrogram analysis. Dendrograms are a hierarchical structure representation (Rosolowsky et al. 2008). There are three control parameters: minimum intensity, minimum number of pixels, and minimum significance. We limited our analysis to the area with a baseline noise of 2 times the median noise over the map (σ_med) for each subregion, so that the dendrogram routine runs over a cube with a relatively homogeneous noise distribution. We use the minimum intensity limit of 2σ_med at each velocity channel and at each position. We set the minimum significance for a leaf to be considered an independent entity to 4σ_med. The number of minimum pixels for an independent leaf is set to 3³ = 27. Three spatial pixels correspond to the spatial resolution of the map. Then we created a 3D mask using all detected leaves, branches, and trunks (we call them clusters in the following). Each spatial and spectral pixel belonging to one of the detected clusters has a value of 1 in this 3D mask. We then derived the integrated intensity map by integrating the 3D spectral cube over this mask. This is equivalent to defining a matching velocity range for integration pixel by pixel, thus avoiding integration over a noisy baseline range without signal. However, this systematically underestimates the integrated intensities by ignoring wing contributions (spatially and spectrally) that fall under the minimum intensity (see Appendix B for more discussion and comparisons of the dendrogram-based methods and the masked moment methods with different threshold parameters).

4.3 Spatial distributions and dynamics

4.3.1 IC 1396A

Figures 3 and 4 show the channel maps of [C II] and [O I], Fig. 5 shows the integrated intensity maps for each observed species, Fig. 6 plots the spectra of selected positions marked in Fig. 5, and Fig. E.1 shows the p–v diagrams for two cuts through IC 1396A. [O I] 145 µm is detected at only a few positions, and thus shown only in Fig. 6. Position 1 is at the [C II] peak, and position 2 is at the [O I] peak as well as the ¹³CO(3−2) local peak (Fig. 5). Position 3 is the position of the Class 0 object PACS-1 (Sicilia-Aguilar et al. 2014), position 4 is the Spitɀer InfraRed Array Camera (IRAC) 8 µm peak at the tip of the BRC. Position 5 is inside the BRC deeper than the ¹³CO peak and [C I] peak (Fig. E.2), and position 6 is where [C II] shows a strong blue wing in the spectrum (Fig. 6).

At the bottom of Fig. 5, we sketch the structures identified in IC 1396A. In the following, we explain these structures in detail.

(i) At −11 and −10.5 km s⁻¹ in the channel maps of [C II] (Fig. 3), blueshifted components are visible in the northern and southern regions, appearing as a clear blueshifted tail in the spectra (position 6 in Figs. 5 and 6).

(ii) At −10 to −8 km s⁻¹, a shell around the cavity with LkHα 349a/c is visible with a velocity gradient: at −10 km s⁻¹, the southwest of LkHα 349 shows a clearer void, while at −8 km s⁻¹ the void is more northeast of LkHα 349. Sicilia-Aguilar et al. (2013) detected [S II] emission at the rim of this void, with particularly strong emission from the northwest region, suggesting that low-mass stars contribute to shaping of the surrounding cloud. The velocity gradient in [C II] can be explained by expansion of the gas, likely also caused by the low-mass stars located inside the cavity. At these velocities, more void-like structures are observed in IC 1396A, with thin bridges and shells surrounding them.

(iii) The outer rim that is prominent in the IRAC 8 µm and PACS 70 µm maps (Fig. 5) is the interface to the ionized region illuminated by the UV radiation from HD 206267. The [C II] peak velocity at the tip of the rim (the eastern side) is at −8 to −9 km s⁻¹.

(iv) As shown in Okada et al. (2012) with lower S/N [C II] data, the integrated intensity map of [C II] correlates well with the Spitɀer IRAC 8 µm map, with a peak in the northern part of the BRC (position 1). Inside the northern part of the BRC (position 2), the velocity profile of CO(3−2) shows a clear flattop shape in contrast to that of ¹³CO(3−2), probably due to its high optical depth. In the northern part of the BRC, [C II] has its peaks closer to the rim, and [O I] 63 µm is stronger inside the BRC, similar to CO(3−2). This indicates a clumpy structure in the northern part to allow the UV radiation to reach the middle of the BRC and excite the [O I].

(v) Position 3 is the position of the Class 0 object PACS-1 with cold dense gas on its western side (Sicilia-Aguilar et al. 2019). We also see that C¹⁸O(3−2) is concentrated west of the illuminated rim, while ¹³CO(3−2) and CO(3−2) are more widely distributed in the BRC. The [O I] 63 µm spectra at positions 3 and 4 show a flat top shape with a slope (the apparent peak of the [O I] 63 µm is blueshifted compared to the [O I] 145 µm peak), and the [O I] 145 µm has the same peak velocity as those of ¹³CO(3−2) and C¹⁸O(3−2), clearly indicating the self-absorption of [O I] 63 µm. The peak velocity of [C I] at lower spatial resolution (Okada et al. 2012) also matches those of ¹³CO(3−2) and C¹⁸O(3−2). We conclude that there is a high column of hot and dense gas at these positions. This is consistent with the strong 70 µm emission there and the temperature distribution derived by Sicilia-Aguilar et al. (2019).

(vi) At position 5, the flat-top spectra of CO(3−2) and even ¹³CO(3−2) indicate a high CO column density. While the [C II] spectrum also appears to have a flat-top, the [¹³C II] analysis shows no evidence of significant optical thickness (Sect. 4.4). Inside the southern part of the BRC, including position 5, [C II] and [O I] 63 µm are generally weak. The [O I] 145 µm is not detected at position 5 (as mentioned above, the baseline around −15 km s⁻¹ is affected by the atmospheric absorption line, but there is no sign of a strong emission line at −8 km s⁻¹), and the [O I] 63 µm is also weak, although self-absorption is not fully excluded due to the limited S/N of [O I] 145 µm. Position 5 is near the peak of [C I] 492 GHz (Fig. E.2; Okada et al. 2012), which is located deeper in the BRC than the C¹⁸O(3−2) peak, close to the CO(3−2) peak (roughly along the lines of positions 3–5). All these are consistent with the scenario that the UV radiation in the southern part is strongly shielded by the dense molecular gas west of the rim, and the PDR emissions are weak in the southern part of the BRC. The cold temperature derived by Sicilia-Aguilar et al. (2019) behind the rim in the southern part also supports our interpretation.

(vii) At a redshifted velocity of −6 km s⁻¹ there is a ring-like structure tracing the northern edge of the BRC and continuing inside the BRC toward the south (west of the C¹⁸O integrated intensity peak). This redshifted ring exists only in [C II], not in the [O I] channel maps (Fig. 4).

The spectra in Fig. 6 show that the peak velocity of [C II] is often blueshifted compared to ¹³CO(3−2) and C¹⁸O(3−2) (position 1–4 and 7). The displacement of the [C II] and ¹³CO(3−2) peak velocity is also clearly visible in the p–v diagrams (Fig. E.1). The difference is small, while it exists, near the edge of the rim (cut a), but it is significant on the backside of the BRC (cut b). Cut (b) also clearly traces the blueshifted tails on the northern and southern edges of the BRC that are visible in the [C II] channel maps (Fig. 3), and ¹³CO(3−2) also follows this curve to a lesser extent. We also found that the velocity profile of [C I] 492 GHz (Okada et al. 2012) is overall very close to that of C¹⁸O(3−2). As mentioned above, [O I] 145 µm at positions 3 and 4 has the same peak velocity as that of ¹³CO(3−2) and C¹⁸O(3−2). This confirms that there is no peak velocity shift between [O I] and CO at these positions. Positions 1 and 2 were not covered by the [O I] 145 µm observations, but a steep decrease in the [O I] 63 µm at the peak velocity of ¹³CO(3−2) at position 2 is also consistent with the [O I] 63 µm being self-absorbed. We discuss the interpretation of the blueshifted [C II] spectra compared to other emission lines in Sect. 5.1.

We quantitatively analyzed the velocity profile in the same manner as in Sicilia-Aguilar et al. (2019); we fit the line with multiple Gaussians and determined the peak velocity and the width of the line where the intensity is above 10% of the peak intensity (Fig. 7). As discussed in Sicilia-Aguilar et al. (2019), we used the fitted spectra to derive these quantities so that they are not affected by the noisy data, and the derived quantities are not affected by the non-uniqueness of the fit because we do not interpret the individual Gaussian components but only use the sum of the fitted Gaussians. The peak velocity here represents the absolute peak, independent of the structure of the line wings (i.e., the skewness of the profile) in contrast to the standard moment 1. The width above 10% intensity better represents weak wings than the width above 50%. We used the map with a resolution of 25″ and limited the analysis to spectra with S/N greater than 8 in order to avoid an over-interpretation of low S/N spectra. Figure 7 shows the derived line width and peak velocity for ¹³CO(3−2), C¹⁸0(3−2), and [С II]. The peak velocity maps clearly show that [CII] is overall blueshifted compared to ¹³CO(3−2) and C¹⁸0(3-2). The [C II] peak velocity map shows that the northwestern and southwestern part of the BRC have bluer peak velocities (< −9.0 km s⁻¹), which is also visible in ¹³CO(3−2). As mentioned above, this is not due to the blue wing, but to the velocity shift of the peak emission, and it is also presented as the curved shapes of the p–v diagram of cut (b) in Fig. E.1. The jump in the peak velocity at the northwest edge of [C II] is due to a stronger peak in the velocity component at −6.5 km s⁻¹ (Fig. 3). The maps of ¹³CO(3−2) and [С II] width have a similar spatial distribution, in the sense that the line width increases toward the southwest and northwest, as shown with C¹⁸O(2−1) by Sicilia-Aguilar et al. (2019), but the [С II] shows an overall broader line width (note that the color scales of the [C II] width are different from those of the other lines). Very wide [C II] widths at the northern and southern edges of the map are due to a blue wing (see Fig. E.l). The region where the width of ¹³CO(3−2) and C¹⁸O(3−2) increases west of the densest part of the BRC match with the peak of the [C I] 492 GHz integrated intensity (Fig. E.2).

Fig. 3 Channel maps of [C II] in IC 1396A (16″ resolution). The color scale is Tmb averaged over a velocity bin centered at the value noted in each panel with a width of 0.5 km s−1. Red crosses mark LkHα 349a (southeast) and c (northwest). HD 206267 is east of the BRC at (RA, Dec)(J 2000) = (21:38:57.619, 57:29:20.54).

Fig. 4 Same as Fig. 3 but for the [O I] 63 µm emission (colors) overlaid with the [C II] 158 µm channel map contours (white with a spacing of 5 K) of IC 1396A.

Fig. 5 Integrated intensity maps (colors, 16″ resolution) of the different tracers (denoted on top of each panel) and overlaid with contours of the [C II] integrated intensity (solid white lines, 16″ resolution) in IС 1396A. Integrated intensities were calculated via a direct integration over the velocity range (−13, −4) km s−1. White areas indicate that a region was not observed, and black areas indicate that no line was detected (positions where the integrated intensities are below 3σ of the baseline are masked). Red crosses mark LkHα 349a (southeast) and с (northwest). Black crosses with numbers mark the positions of the spectra shown in Fig. 6. The bottom-right panel shows the schematic decomposition. Different colors indicate a difference in velocity (red represents the most redshifted component and blue is the most blueshifted), and the Roman numerals correspond to the description in the text. The black box border indicates the size of the images in the other panels.

Fig. 6 Spectra at selected positions in 1С 1396A, as marked in Fig. 5 (16″ resolution). The [C II] and [O I] spectral baselines are offset for a better visibility of the individual lines, and the vertical dotted lines are for guidance when comparing the peaks between different lines. [O I] 145 µm spectra are affected by the atmospheric ozone feature at velocities of < −10 km s−1 (see the main text).

Fig. 7 Peak velocity (upper panels) and line width at 10% of the peak intensity (lower panels) in units of km s−1 for 13CO(3−2), C18O(3−2), and [C II] (from left to right) in IС 1396A at 25″ resolution. The contours are the integrated intensity of [C II] at 25″ resolution. Note the different scale of the color wedge used for the [C II] line width, as the spatial distribution of the [C II] width would not be visible at the color scale used for the other two lines.

Fig. 8 Channel maps of [С II] in 1С 1396В (16″ resolution). The bottom-right panel shows the Hα image from the archive of the INT (2.5 m Isaac Newton Telescope)/WFC (Wide Field Camera) Photometric Hα Survey of the Northern Galactic Plane (IPHAS; Barentsen et al. 2014) for comparison.

4.3.2 1С 1396B

Figure 8 shows the channel maps of [C II], Fig. 9 shows the integrated intensity maps, Fig. 10 plots the spectra of selected positions marked in Fig. 9, and Fig. E.3 shows the p-v diagrams for two cuts through 1С 1396B.

The spatial distribution of the [C II] integrated intensity not only follows AKARI/Infrared Camera (IRC) 9 µm, but also generally agrees with the CO(3−2) map (Fig. 9), which is not the case for the other BRCs we study in this paper. The [C II] channel maps of 1С 1396B show a rich velocity structure in this region (Fig. 8). In the right panel of Fig. 9, we sketch the structures identified in 1С 1396B.

(i) The southern region appears prominent at a velocity of −9 to −8 km s⁻¹, similar to the rim in 1С 1396A. In the MIR and optical images (Fig. 1), a bridging structure between 1С 1396A and В is visible, and a similar velocity between the southeastern part of 1С 1396B and 1С 1396A indicates that they are indeed connected.

(ii) At the velocity of −8 to −7km s⁻¹, there is a strong linear structure, where the southeast is more blueshifted, and it becomes more redshifted toward the northwest, connecting to the strongest cloud components at the velocity of < −7 km s⁻¹. Figure E.3 shows the p–v diagrams along this structure (cut (b)). This linear structure is visible in both [C π] and CO(3−2) (Fig. 9 and positions 5 in Fig. 10). There are no known YSOs along this structure (Barentsen et al. 2011; Nakano et al. 2012). Similar straight features have been observed in other regions, for example the “Stick” in Orion (Kong et al. 2021) or a linear feature in NGC 7538 (Beuther et al. 2022). Unlike the Stick, which is modeled by collision-induced magnetic reconnection (Kong et al. 2021), the linear feature here in IС 1396B has no obvious ring or double line structure. It is also unlikely that the direction of the magnetic field is parallel to this straight feature, because three out of the four BRCs measured by Soam et al. (2018) have the magnetic field parallel to the Galactic plane, which is rather perpendicular to this line. This may just be the UV-illuminated edge of a part of the structure connecting IС 1396B and A, also visible as a dark lane in the optical (see the Hα panel of our Fig. 8 as well as Fig. 2 of Barentsen et al. 2011), which happens to look straight. The [C II] spectrum along this straight feature (position 5 in Fig. 10) shows a slightly blueshifted profile and a blue wing compared to CO(3−2).

(iii) At a velocity of −6 to −3.5 km s⁻¹, the channel maps show two bright [C II] clouds (Fig. 8). The Hα map (the bottom right panel of Fig. 8, Barentsen et al. 2011) outlines these clouds as a silhouette with bright rims, strongly suggesting that these clouds are in front of the ionized gas along the line of sight.

(iv) There is a clearly separated velocity component in both [C II] and CO(3−2) at +1.5 to +3 km s⁻¹ (e.g., position 1 in Fig. 10). The p–v diagrams in Fig. E.3 do not indicate any connection between the positive and negative velocity components. The positive velocity components form a large rim structure, facing east-northeast, although the direction of the exciting star is more east-southeast. The p–v diagram of cut (b) in Fig. E.3 does not show a clear velocity gradient of the positive velocity component along the rim. There is no corresponding structure visible in the Hα map. This may be a separate cloud behind the ionized gas, illuminated from the front by the exciting star. The northern part of the velocity component >0 km s⁻¹ is dominated by the [C II] emission and CO(3−2) is very faint (west side of the cut (b) in Fig. E.3 and position 7 in Fig. 10).

Except for a blue wing in the [C II] emission around position 5 and 6 (Fig. 10), the [C II] profile shows no significant broadening compared to CO(3−2).

Fig. 9 Integrated intensity maps (colors, 16″ resolution) of the different tracers (denoted on the top of each panel) overlaid with contours of [C II] integrated intensity (solid white lines) in IС 1396B. White and black areas have the same meaning as in Fig. 5. Integrated intensities were calculated over the components detected by the dendrogram analysis (see the main text). Black crosses with numbers mark positions whose spectra are shown in Fig. 10. The right panel shows the schematic decomposition as in Fig. 5.

Fig. 10 Spectra at the selected positions in IС 1396B marked in Fig. 9 (16″ resolution).

4.3.3 IС 1396D

IС 1396D shows a complex structure with multiple velocity components as shown in the channel maps (Fig. 11). The CO and [C II] integrated maps (Fig. 12) look very different, but individual velocity components are often common in the spectra (Fig. 13). At the bottom of Fig. 12, we sketch the structures identified in IС 1396D.

(i) A cometary-shaped cloud in the north is a BRC with an IRAS source identified by Sugitani et al. (1991) (SF034). [С II] at the rim is faint but detected (Fig. 11 around −5km s⁻¹). The rim is also visible in the IRAC 8 µm image of Rebull et al. (2013) and AKARI9 µm (Fig. 12). The Hα image (Fig. 11 bottom right panel) also shows a bright rim just outside the dark cometary cloud, indicating that they are located in front of the ionized gas along the line of sight.

(ii) CO filaments extend to the south, and the same structure is visible in [C II] (−8 to −5 km s⁻¹). This filamentary structure appears as dark clouds in Hα (Fig. 11 bottom right panel), with a sharp boundary on the eastern side toward the exciting source, indicating that they are also located in front of the ionized gas.

(iii) The integrated intensity peak of [C II] is located west of the CO filament, where the 9 µm emission is also bright, and further away from the star in projection than the CO filament. The velocity structure is complex (Fig. E.4) and some components have a very broad line profile (positions 4–7 in Fig. 13). This indicates that the region with strong [C II] receives the UV radiation from the exciting source without being shielded by the BRC and the CO filament in the north-south direction, but they are displaced along the line of sight. This is also supported by the gradual increase in Hα emission from the CO filaments toward the west, and the fact that the boundary to the next cloud (at the [C II] peak) does not appear as sharp as the blueshifted BRC and the CO filament. The presence of multiple velocity components in [C II] also suggests a 3D distribution of C⁺. An ionized gas origin of the [C II] emission is unlikely, or at least not significant, because of the associated weak CO emission (−5 to 0 km s⁻¹ seen at position 6) and a good correlation with the 9 µm emission, which traces the polycyclic aromatic hydrocarbon (PAH) emission.

In a simple picture of an expanding shell structure created by the evolving H II region, the blueshifted clouds are in front of the H II region along the line of sight, moving toward us, and the redshifted clouds are behind the H II region. This picture is consistent with the structure observed in IC 1396D: the blueshifted BRC and the associated north-south CO filamentary structure are located closer to the observer, while the western cloud with strong [C II] is located behind the ionized gas.

When CO and [C II] velocity components are associated with each other (e.g., position 1 in Fig. 13), in most cases no significant broadening of the [C II] emission is visible. At position 1, this means that in the BRC with a cometary shape, there is no indication of a velocity shift in any line emission, and thus no sign of ablated surface or photoevaporation flows.

Fig. 11 Same as Fig. 8 (16″ resolution) but for the [C II] in IС 1396D.

Fig. 12 Integrated intensity maps (colors, 16″ resolution) of the different tracers (denoted on the right of each panel) overlaid with contours of [C II] integrated intensity (white lines) in IC 1396D. White and black areas have the same meaning as in Fig. 5. Integrated intensities were calculated over the components detected by the dendrogram analysis (see the main text). Black crosses with numbers mark positions whose spectra are shown in Fig. 13. The bottom panel shows the schematic decomposition as in Fig. 5.

4.3.4 IC 1396E

Figure 14 shows the channel maps of [C II], Fig. 15 shows the integrated intensity maps, Fig. 16 plots the spectra of selected positions marked in Fig. 15. At the bottom of Fig. 15, we again sketch the structures identified in IC 1396E.

(i) The channel maps of [C II] (Fig. 14) show that the rim is clearly outlined at a velocity of about 1 km s⁻¹, while the inside of the BRC or the northern part of the map has somewhat blueshifted emission (− 1 to 0 km s⁻¹). The [C II] integrated intensity traces the edge of the bright rim, illuminated by the exciting sources from the south (Fig. 15). The rim is also visible in ¹³CO at 1–2 km s⁻¹. We discuss the implication of the velocity gradient toward the inside of the BRC in Sect. 5.

(ii) In contrast to the [CII] emission, CO stems mainly from the inside of the BRC. As shown in previous works (see Sect. 2), the ¹²CO lines show a very prominent outflow signature in their spectra, together with a strong self-absorption feature in IC 1396E (Fig. 16). BIMA2, an intermediate-mass protocluster and the driving source of the central outflow, is located west of position 3 in Fig. 15. Source C, the driving source of the northern outflow (Codella et al. 2001), is located east of position 7 (between the two CO peaks).

(iii) The channel maps of CO at a larger scale (Fig. E.5) than the area mapped in [C II] show the two tails with different velocities as already discussed by Patel et al. (1995) (Sect. 1). The eastern tail connects to the rim (i) at its most redshifted velocity.

The effect of outflows on [C II] is not obvious in the individual spectra (Fig. 16), but the analysis of the velocity widths (Fig. 17) shows its influence in the [C II] spectra. We performed the same analysis as in IC 1396A, fitting the spectra with multiple Gaussians. To account for the skewed shape of the spectra, we derived the red and blue velocity widths separately, namely as the width from the peak velocity to the velocity where the intensity drops to 10% of the peak intensity on the respective side (Sicilia-Aguilar et al. 2019). Figure 17 shows the results for ¹³CO(3−2) and [C II]. The ¹³CO(3−2) width shows significant enhancements at the clumps produced by the central outflow. In the [C II] width maps, there are two local peaks in the middle of the BRC, associated with the H₂ knots identified by Beltrán et al. (2009); one peak in the red velocity with the chain of H₂ knots G, and the other in the blue velocity with the H₂ knots E and P. These results suggest that the [C II] emission traces the ambient gas interacting with the outflows.

The width of the [C II] red velocity (second panel of Fig. ) shows that the eastern side of the rim has a broader line in [C II]. It is visible as a wing in the [C II] spectrum at position 4 (Fig. 16). The CO(3−2) peak and the H₂ knots are located about 1′ (0.27 pc) to the north from the rim. At position 4, the peak velocity of [C II] and ¹³CO(3−2) match, and this red wing in [C II] may be a sign of photoevaporation flow (see Sect. 5). On the other hand, at the intensity peak of the rim (position 1), the ¹³CO(3−2) peak velocity and the absorption peak of CO(3−2) correspond to the dip in the [C II] line profile; thus, the [C II] profile is likely affected by self-absorption.

Fig. 13 Spectra at the selected positions in IC 1396D marked in Fig. 12 (16″ resolution).

Fig. 14 Same as Fig. 8 (16″ resolution) but for the [C II] in IC 1396E.

Fig. 15 Integrated intensity maps (colors, 16″ resolution) of the different tracers (denoted on the top of each panel) overlaid with contours of [C II] integrated intensity (white lines) in IC 1396E. White and black areas have the same meaning as in Fig. 5. The bottom right panel shows the schematic decomposition as in Fig. 5. Integrated intensities were calculated by direct integration over the velocity range of (−5, 8) km s−1, except for CO(3−2) with the range of (−10, 10) km s−1. Red pluses are the positions of outflow-driving sources: BIMAl-3 (southern ones, Beltrán et al. 2002) and the source C (northern one, Codella et al. 2001). Black crosses with numbers mark positions whose spectra are shown in Fig. 16.

Fig. 16 Spectra at the selected positions in IC 1396E marked in Fig. 15 (16″ resolution).

Fig. 17 Red wing (left two panels) and blue wing (right two panels) line width at 10% of the peak intensity for 13CO(3−2) and [C II] in IC 1396E at 25″ resolution. Red pluses are the positions of the outflow-driving sources: BIMAl–3 (southern ones, Beltrán et al. 2002) and the source C (northern one, Codella et al. 2001) as in Fig. 15. The crosses on the [C II] panels are the positions of the H2 knots from Beltrán et al. (2009), with the size of the crosses corresponding to the H2 flux. Some of the H2 knots clusters are labeled according to Beltrán et al. (2009).

4.4 [¹³C II] in IC 1396A

We analyzed the [¹³C II] hyperfine emission lines in IC 1396A, where the S/N of the [C II] observations is the highest. The [¹³C II] ²P_2/3-²P_1/2 transition splits into three hyperfine components with relative intensities of 0.625, 0.25 and 0.125 (Ossenkopf et al. 2013). The strongest component (F = 2−1) is 11.2 km s⁻¹ away from the [C II] line and often blended with its wing, but it can be recovered using a polynomial fit to the wing of the [C II] line. The S/N of individual positions in the map is not sufficient for a [¹³C II] analysis; we thus averaged the spectra over larger areas. To determine which area to average, we performed a 2D dendrogram analysis on the [C II] integrated intensity using the same parameters as in the 3D dendrogram analysis (Sect. 4.2) except for the number of pixels for a leaf. A 3D dendrogram is not appropriate for this purpose because the area to average should be common at all velocities, so we change the number of pixels for a leaf to 3², suitable for a 2D dendrogram. Figure E.6 shows the seven clusters identified. We averaged the spectra within each cluster to disentangle the [¹³C II] emission from the noise floor. We fit the baseline around each hyperfine line by a second-order polynomial and subtracted it (see Appendix C for the discussion on the impact on the order of the baseline). We have a clear detection of the [¹³C II] F = 2−1 line in most of the clusters. In a second step we average the three hyperfine lines weighted by their relative intensities as described in Guevara et al. (2020). The averaged [¹³C II] emission scaled up by the local ¹²C/¹³C abundance ratio α = 68 ± 15 (Milam et al. 2005) is shown in Fig. 18. We find that the scaled up [¹³C II] line only slightly overshoots the [C II] line, indicating weak optical depth effects.

We then quantified the velocity-resolved optical depth. Since the [¹³C II] line is weak, it is easily affected by baseline or noise fluctuations. Therefore, we performed a fit to the averaged [¹³C II] line using a superposition of two Gaussian profiles. In some velocity bins the scaled [¹³C II] are weaker than the [C II] intensity, which suggests fully optically thin [C II] within errors. To guide the fit into a physically meaningful solution and to avoid numerical problems when calculating the optical depth, we replaced the scaled [¹³C II] with the [C II] intensity for those data points where the former is smaller than the latter during the fit. This is equivalent to assigning an optical depth of zero. This adjustment was necessary due to significant baseline fluctuations, mostly at the wings of the [¹³C II] line where the S/N is lowest. The fitted [¹³C II] spectra are shown by the orange curve in each panel in Fig. 18. We use the fitted [¹³C II] for the following analysis. We assume that the [C II] emission originates from a layer that can be described by a single excitation temperature. The lack of absorption dips in the [C II] spectra that would result from additional colder material along the line of sight supports our approximation. In this simple scenario, we can determine the optical depth by comparing the intensities of the two isotopes (Ossenkopf et al. 2013): $$\frac{T_{\text{mb}{,}^{12}\text{C}}(v)}{T_{\text{mb}{,}^{13}\text{C}}(v)}=\frac{1-e^{-\tau (v)}}{\tau (v)}\alpha .$$(1)

The corresponding error is derived following the Eq. (2) in Kabanovic et al. (2022). It propagates the uncertainties of [C II] and [¹³C II] intensities, and the ¹²C/¹³C abundance ratio α, using the first-order Taylor expansion of each component. The dominant term in the optical depth error comes from the error of [¹³C II] (>= 75%). The resulting velocity-resolved optical depth is shown by the black dots with error bars in Fig. 18. As mentioned above, we forced the optical depth to zero when the [C II] intensity is equal or exceeds the α-scaled [¹³C II] intensity. Consequently only positive error bars indicate physical meanings, and thus shown. The derived τ values are mostly smaller than 1 as indicated by the dashed cyan line in each panel, and it is compatible with the optically thin [C II] emission within the error bars. While it cannot be completely excluded that the [C II] emission is optically thick in a small area of locally high column density, we conclude that the optical depth does not significantly affect the observed [C II] emission.

4.5 Heating efficiency in IC 1396A

We investigated the heating efficiency by calculating the ([C II] 158 µm + [O I] 63 µm)/FIR in IC 1396A. We fit a gray body with the frequency dependence of the dust absorption coefficient β = 2 to the PACS 70 µm, 160 µm, and SCUBA-2 850 µm and derived the total FIR as the integral of the fitted function. In contrast to the fitted temperature and dust column density, the integrated FIR is insensitive to the assumed β. Figure 19 shows the spatial distribution of ([C II] 158 µm + [O I] 63 µm)/FIR. It has a value of 3–4% in the cavity around LkHα 349c and the northern part of the BRC. When interpreting this quantity as heating efficiency and comparing it with the values in previous works, two things have to be considered. One is that we refer to FIR, the total FIR flux only, while some literature reference is made to the “TIR.” The TIR is defined as the total infrared flux at 3–1100 µm (Dale & Helou 2002) and better traces the total input energy, because it includes the emission from all dust grains, including PAHs and very small grains (VSGs). The ratio of FIR to TIR varies depending on the environment but is typically about 50% (Pabst et al. 2021; Okada et al. 2013). Here we use FIR due to lack of sufficient MIR data. Another difference is that the cooling line flux does include the [O I] emission, which is often ignored due to its rare availability in observations compared to [C II]. In IC 1396A the ratio of [O I] 63 µm/[C II] ranges from 0.3–1.2. Taking these two differences into account, the heating efficiency in IC 1396A is similar to the upper bounds of the data shown, for example, in Pabst et al. (2021).

The spatial distribution of the heating efficiency in Fig. 19 indicates that the variation is not due to a decrease in the photoelectric heating efficiency on PAHs. The photoelectric heating efficiency is theoretically expressed as a function of the charging parameter, γ = G₀T^0.5/n_e, where G₀ is the radiation field strength, T is the gas temperature, and n_e is the electron density (Bakes & Tielens 1994). As discussed in Sect. 4.3.1, the UV radiation in the southern part of the BRC is likely shielded by the dense molecular gas west of the rim. A low UV field, cold and dense gas results in a low value of γ, which should increase the abundance of neutral PAHs and thus increase the heating efficiency. However, Fig. 19 presents a contradictory spatial distribution; the heating efficiency is clearly lower in the southern region than in the northern region. Furthermore, the estimated⁸ γ in IC 1396A is on the order of 100, a value for which the heating efficiency is high and not sensitive to the exact value of γ (Bakes & Tielens 1994; Weingartner & Draine 2001b). This also supports the view that the spatial difference of the measured heating efficiency cannot be attributed to the charging status of PAHs.

One possible reason for different heating efficiencies in the southern and northern regions of the BRC is a difference in the dust properties. Weingartner & Draine (2001b) modeled the photoelectric heating efficiency with different size distributions of dust grains listed in Weingartner & Draine (2001a). Models with strong enhancement at the size of < 10 Å (the C abundance in very small grains is parameterized as b_c) show higher efficiencies, and models with R_V = 5.5 (R_V is the ratio of visual extinction to reddening), representing dense environments, have lower efficiencies compared to models with R_V = 3.1 (see also Fig. 6 in Okada et al. 2013). If the dust in the southern part is composed of fewer small grains relative to the large grains (i.e., larger R_V and smaller b_c), this is consistent with a lower heating efficiency in the southern region compared to the northern region.

Fig. 18 Averaged spectra in each dendrogram cluster (Fig. E.6) in IC 1396A except cluster 2, which does not have a high enough S/N for the [13C II] analysis. Colored stepped data are the [C II] spectra (thick) and the combined and scaled [13C II] spectra (thin). The orange curve is the fit to the [13C II] spectra using two Gaussian profiles. Derived optical depths of [C II] are shown as black data points with error bars. See the main text for details.

Fig. 19 ([C II]+[O I] 63 µm)/FIR integrated intensity ratio in IC 1396A at 25″ resolution. A jump in values along the shape of the rim is reflecting the [O I] 63 µm detection limit (see Fig. 5 for the area where [C II] and [O I] 63 µm are detected).

Fig. 20 Density plot of the correlation between [C II] integrated intensities and IRAC 8 µm (IC 1396A and E) or IRC 9 µm (IC 1396B and D) at 25″ resolution. The solid line is the relation derived in the Orion A region by Pabst et al. (2021), adapted to this study (see the main text). The dotted line indicates the slope of a linear relation (slope of unity in the log scale) with an arbitrary scaling (offset in the log scale).

4.6 [C II]-MIR correlation

In all subregions a clear correlation between the [C II] integrated intensities and the MIR photometric observations (IRAC 8 µm or IRC 9 µm) is observed (Fig. 20). We compare it with the results in the Orion A region (Pabst et al. 2021), shown as a solid line in Fig. 20. We note that their relation is adapted for Fig. 20 because they use the extended calibration factor (0.74) and filter width (2.9 µm) for the IRAC 8 µm flux (Pabst, priv. communication), which is not considered in this study. The correlation in IC 1396 agrees with that in Orion A for the flux range of [C II] (≥ 10⁻⁷ W m⁻² sr⁻¹ = 10⁻⁴ erg s⁻¹ cm⁻² sr⁻¹) in IC 1396A and B, while we clearly see a turnover in fainter areas. This turnover is probably due to emission along the line of sight, because the [C II] observations are velocity-resolved and only the velocity components associated with the individual BRCs are extracted, while the IRAC 8 µm and IRC 9 µm represent the total emission along the line of sight. If we subtract about 5 MJy sr⁻¹ from the IRAC 8 µm or IRC 9 µm data, the full range of the data in Fig. 20 follow the correlation of the black line in IC 1396A, B, and D.

In IC 1396E the [C II] emission is overall about a factor of 2 weaker than in other regions in the relation to the 8 µm or 9 µm. This could be either due to the contribution of very hot thermal dust to 8 µm in IC 1396E with active YSOs and jets inside the BRCs, or due to a larger [O I] contribution to the cooling in a higher density environment in IC 1396E, although the latter remains speculative without available [O I] data.

5 Discussion

5.1 Searching for the signature of photoevaporation

The interpretation of the velocity structure depends on whether the BRCs are located in front of or behind the ionized gas (Okada et al. 2015). As described in Sect. 2, the optical and infrared images suggest that IC 1396A is located in front of the ionized gas and illuminated from the backside along the line of sight. Blueshifted lines (shifted relative to the other spectral lines or other parts of the BRC) trace gas moving faster toward the observer. As shown in Sect. 4.3.1, the peak velocity of [C II] in IC 1396A is overall blueshifted compared to ¹³CO(3−2). This [C II] velocity structure is inconsistent with a photoevaporating flow, which would be a flow toward the ionized gas and thus red-shifted. We also find no other evidence of photoevaporating gas in [C II] in IC 1396A. The redshifted ring in the channel maps of IC 1396A (at −6 km s⁻¹) is not consistent with photoevaporating gas because part of this ring is also visible in ¹³CO(3−2) (position 4 in Fig. 6).

For IC 1396E, comparison with the mean velocity of the ionized gas from Pedlar (1980) (−1.9 km s⁻¹) would indicate that the BRC is behind the expanding ionized gas (Patel et al. 1995). However, the central velocity derived from the H166α observations in Pedlar (1980) actually ranges from −4.5 km s⁻¹ to +0.6 km s⁻¹ at different positions in IC 1396, with an error of ~2 km s⁻¹, and thus the velocity alone does not give a clear picture of the geometry. The sharp boundary in the Hα image (Fig. 14 bottom right panel) rather suggests that IC 1396E is in front of the ionized gas. On the other hand, the contrast between the inside and outside of the BRC in Hα is not as strong as in IC 1396A (Barentsen et al. 2011). If the BRC is completely in front of the ionized gas, Hα should be almost completely extinct, considering that the extinction is estimated to be A_υ ~ 20 mag (Nisini et al. 2001). IC 1396E could be located in the middle of the ionized gas along the line of sight. With this uncertainty in the geometry, it is difficult to conclude whether the redshifted component in the eastern side of the [C II] rim in IC 1396E is a sign of a photoevaporating flow. It potentially traces the photoevaporating gas leaving the BRC on the backside of the rim along the line of sight toward the ionized gas, but it is not excluded that this is an additional local distortion by the outflows.

In both IC 1396A and IC 1396E, there is a gradient of the velocity in [C II] from the rim (most redshifted) to the rear side of the BRC (the side opposite to the exciting sources), and the blue wings at northwest and southwest of IC 1396A are prominent. Our interpretation is that the [C II] gas comes dominantly from the surface of the BRC, and is blown away from the BRC toward the observer, and from the rim to the rear side of the BRC. The curved shape of the p–v diagram in IC 1396A (Fig. E.1) supports this interpretation of ablation flows.

5.2 3D kinematic structures and comparison of different BRCs

Pelayo-Baldárrago et al. (2023) identified seven subclusters in IC 1396 based on parallax and proper motion measurements from Gaia. Subclusters A, B, E, and F have a large number of members to be statistically significant and are relatively free from contamination by sources not associated with IC 1396. As mentioned in Sect. 1, they found that these subclusters are clearly distinct in proper motion and partly in age. Inside the field of our observations, IC 1396A has a clear association with their subcluster A (six Gaia sources in the field of our [C II] map) and IC 1396E with their subcluster B (nine Gaia sources in the field of our [C II] map). IC 1396D has a dispersed population (with each of their subclusters A, E, and F having a single Gaia source within the field of our [C II] map), and IC 1396B has only one Gaia source belonging to subcluster A within the field of our [C II] map.

The distance to the subclusters A and B determined by Pelayo-Baldárrago et al. (2023) is 908 ± 73 pc and 911 ± 75 pc, respectively. The projected diameter ofIC 1396 on the sky (basically the width of Fig. 1) is approximately 45 pc at a distance of 925 pc. Given the considerable uncertainty in the derived distance, it remains unclear whether IC 1396E is located in front or behind the ionized gas. A comparison of the distances of the Gaia sources that overlap with the [C II] map in IC 1396A and IC 1396E does not provide a conclusive result either.

We can, however, combine the mean proper motion of the subclusters A and B (Pelayo-Baldárrago et al. 2023) and the [C II] line-of-sight velocity to derive the relative 3D motion between IC 1396A and E. Table 3 summarizes the velocity calculated by the proper motion and the [C II] line-of-sight velocity. It is evident that IC 1396A and IC 1396E are moving away from each other in 3D (both in the plane of the sky and along the line of sight), which is consistent with an expansion of the entirety of IC 1396 (Patel et al. 1995; Pelayo-Baldárrago et al. 2023).

IC 1396A and E have more, young sources listed by Pelayo-Baldárrago et al. (2023) and embedded YSOs identified in the literature, as mentioned in Sect. 1, while those in IC 1396B and D are less prominent. This difference is not attributable to the current projected distance of each BRC from the exciting sources. It rather must be associated with the presence of high density clumps in IC 1396A and E (Sicilia-Aguilar et al. 2014; Serabyn et al. 1993). The question whether these BRCs have a higher gas density in general will be addressed in a follow-up paper by deriving physical properties using PDR models.

6 Summary

We investigated the dynamical and physical structures of four BRCs in IC 1396 using velocity-resolved observations of [C II] 158 µm, [O I] 63 µm, and 145 µm with (up)GREAT on board SOFIA, combined with low-J (3–2) transitions of CO, ¹³CO, and C¹⁸ O from the JCMT archive. The individual results are as follows:

• In IC 1396A, the peak velocity of [C II] is blueshifted compared to ¹³CO(3–2) and C¹⁸O(3–2), especially on the back side (the side farthest from the exciting sources). Because IC 1396A is likely in front of the ionized gas along the line of sight, the blueshifted [C II] cannot be explained by photoevaporation. It instead indicates that gas was blown off from the rim to the back side of the BRC.

• The spatial distribution and the self-absorption signature of the emission lines suggest that the northern part ofIC 1396A is clumpier and the UV radiation reaches the middle of the BRC, while the southern part is strongly shielded from the UV radiation by a high column of hot and dense gas.

• Our analysis of [¹³C II] in IC 1396A did not reveal any evidence that the [C II] emission has a significant optical depth.

• The heating efficiency measured by the ([C II] + [O I] 63 µm)/FIR is 3–4% in IC 1396A. This is similar to the upper limits from previous studies of other regions, taking the different definitions of efficiency in different studies into account. The southern part of IC 1396A has a lower efficiency than the northern part, which may be due to fewer small grains and a high density.

• In IC 1396B and IC 1396D, distinct velocity components are observed in both the [C II] and CO emission, with no indication that they are physically connected. They are probably independent clouds behind and in front of the ionized gas along the line of sight.

• The width of the [C II] emission in IC 1396E indicates an interaction between outflows and the ambient gas.

• Except for IC 1396E, the [C II] intensity shows a good correlation with the 8 µm or 9 µm emission in an intensity range that is significantly higher than that of the background emission, which is quantitatively consistent with previous studies.

• By combining the line-of-sight [C II] velocity with the proper motion measured by Gaia, we confirm that IC 1396A and IC 1396E are moving away from each other in 3D, which supports the scenario that the entire region is expanding.

The interpretation of the [C II] and [O I] 63 µm spectral shape and intensity strongly depends on their optical depth. This study demonstrates that the combination of [¹³C II] with [C II], and [O I] 145 µm with [O I] 63 µm leads to conclusive results because we can disentangle the effect of optical depth and the intrinsic dynamics in the velocity profile.

Acknowledgements

This work is based in part 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 carried out within the Collaborative Research Centre 956, sub-project A4 and C1, funded by the Deutsche Forschungsgemeinschaft (DFG) - project ID 184018867. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. This work is partly based on observations with AKARI, a JAXA project with the participation of ESA. This research made use of astrodendro, a Python package to compute dendrograms of Astronomical data (http://www.dendrograms.org/). The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. N.S. acknowleges supprot by the BMWI via DLR, project number 50 OR 2217 (FEEDBACK-plus). S.K. acknowledges support by the BMWI via DLR, project number 50 OR 2311 (Orion-Legacy). Y.O. et al. gratefully acknowledge the Collaborative Research Center 1601 (SFB 1601 sub-project A6 and B2) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 500700252.

Appendix A Baseline correction using PCA

Principal component analysis was used to characterize and correct for the instrument-origin baseline structures in the spectra. Its application to different types of observations is described in detail in Buchbender et al. (in prep.). Here, we summarize the procedure applied to the chopped OTF data presented in this paper and present the improvement in the data quality.

The basic idea is to extract principal components from the OFF (sky) measurements, and subtract these components from the ON spectra. This ensures that the corrected baseline components are not affected by the line emission in on-source spectra. The whole procedure is written in python with the pygildas interface, so that the program can be called from GILDAS (Buchbender et al. in prep.). Figure A.1 shows as an example the difference of adjacent OFF measurements (OFF-OFF) in a chopped OTF line for the HFAV pixel 0, clearly showing a common broad baseline pattern. This is inferred to be caused by the gain drift of the instrument, and it also appears in the ON spectra. While the patterns are different for different pixels, the pattern of the same pixel is stable in each observing leg. Therefore, we construct the principal components from the set of OFF-OFF spectra for each pixel in each flight. We used only the velocity range of –57 to 0 km s^–1 for [OI] 63µm (HFA) and –57 to 30 km s^–1 for [C II] (LFA), which covers the expected range of the emission line well and excludes the velocity of the atmospheric line. We resampled each OFF-OFF spectrum in velocity to the goal spectral resolution (0.5 km s^–1), subtracted the mean over this velocity range as an offset from each OFF-OFF spectrum, and normalized them with respect to the standard deviation. For all six flights, the number of OFF-OFF measurements per pixel per flight is more than 1600, which is statistically sufficient to extract the principal components.

Figure A.2 shows examples of the first three principal components (HFAV pixels 0 and 2, flight on March 6, 2020). Mathematically, one can derive as many components as the number of input spectra. We used only those components that show a broad baseline structure by comparing the standard deviation of a narrow and broad velocity range. We measured a running mean and standard deviation of 5 channels in the 0.5 km s^–1 bin, and we used the component if the median of the running standard deviation is less than 0.7 of the standard deviation of the entire velocity range used. In the examples in Fig. A.2, only component 1 is used for HFAV_0, while three components are used for HFAV_2. An alternative criterion could be the percentage of variance explained by each of the selected components. This works well when there is a clear dominant component, but it is difficult to decide on a consistent number as a criterion for differently behaving pixels and receivers.

Fig. A.1 Difference of adjacent OFF measurements in one of the chopped OTF lines in the flight on March 6, 2020, for HFAV pixel 0. Each spectrum is smoothed to the velocity resolution of about 0.5 km s-1. The vertical axis is in Kelvin.

Fig. A.2 Examples of derived principal components (HFAV pixels 0 and 2 for the flight on March 6, 2020).

Fig. A.3 Examples of the baseline rms with and without the PCA (HFAV pixels 0 and 4 for the flight on March 6, 2020).

After selecting the components as above, the contributions of these components to the individual ON spectra are calculated and then subtracted. The ON spectra are also resampled in velocity to the goal spectral resolution before the principal components are subtracted.

Figure A.3 shows typical histograms of the baseline rms over the entire velocity range used. The plot of HFAV_0 represents cases where the PCA clearly improves the baseline quality, where the rms distribution becomes narrower and generally shifts to lower noise after the PCA is applied. The plot of HFAV_4 represents cases where the original data are already good enough that the PCA has no significant effect. This is the case for all LFA data.

Appendix B Comparison of different methods for generating integrated intensity maps

In this appendix we compare the integrated intensity maps obtained using different methods: direct integration of the spectra with a fixed velocity range, 3D dendrograms, and the masked moment method. It should be noted that the purpose of the integrated intensity maps presented in this paper is to identify spatial structures and to compare the spatial distribution of different emission lines. The precise selection of methods and thresholds does not significantly impact the outcome of the analysis in this paper.

Fig. B.1 Integrated intensity maps of IC 1396D derived using different methods and thresholds at 16″ resolution. The upper three panels are from the 3D dendrogram method, with thresholds of 1σ, 2σ, and 3σ from left to right. The middle three panels are from the masked moment method, also with the threshold of 1σ, 2σ, and 3σ from left to right. The bottom panel is the direct integration of the spectra with the velocity range from –12 km s–1 to +10 km s–1. See the main text for details of the different methods. Black crosses with numbers mark positions where the masks along the spectral axis are shown in Fig. B.3.

Fig. B.2 Integrated intensity maps of IC 1396A derived using the 3D dendrogram method, with thresholds of 1σ (left) and 3σ (middle), and the direct integration of the spectra in the velocity range –13 km s–1 to –4 km s–1. The color scale is the same for all panels.

Fig. B.3 Masks along the spectral axis at the positions marked in Fig. B.1 derived using different methods for calculating the integrated intensity. The blue dotted line shows the 1σ level. Note that the mask is calculated by averaging the 3 × 3 spatial pixels, but the spectrum shown in each panel is the one at the central pixel.

When using the 3D dendrogram analysis to create a mask for integrated intensity maps, the parameter of a minimum significance of a leaf is less important because we do not distinguish the individual detected clusters, while the minimum intensity limit directly affects the results because its high value cuts off weak emission. The top three panels in Fig. B.1 compare the derived integrated intensity maps using 1σ, 2σ, and 3σ as the minimum intensity threshold, where σ = 0.6 K is the median rms noise in IC 1396D, and the bottom panel shows the direct integrated intensity using the velocity range from –12 km s^–1 to +10 km s^–1. The 3σ threshold not only loses extended faint emission, but also suppresses the intensities at bright regions.

On the other hand in Fig. B.2 we show the same comparison for IC 1396A, which has an overall higher S/N. This demonstrates that the results are essentially indistinguishable for high S/N data.

We also compare the integrated intensity maps using the dendrogram-based masks with the masked moment method (Adler et al. 1992). This method also creates a 3D mask for integration based on the 3D average of a small volume around each pixel. We used 3 pixels here to be comparable to the minimum cluster size of the dendrogram analysis. For each spatial and spectral bin, we calculated an average over 3 × 3 × 3 (spatial and spectral) pixels, and if the average was greater than a certain threshold, we assigned a value of 1 to the 3D mask of the central pixel. We then derived the integrated intensity map by integrating the 3D spectral cube over this mask, as with the dendrogram-based masks. We also used 1σ, 2σ, and 3σ as threshold values (Fig. B.1 middle row). Overall, the derived integrated intensity maps look very similar between the dendrogram-based mask method and the masked moment method when using the same threshold.

Figure B.3 shows examples of the mask at the selected positions indicated in Fig. B.1. The masks over a spectral axis at these positions are shown as 0-1 step plots. The example of position 1 clearly shows that the 3σ threshold masks too much structure, indicating that a 2σ or 1σ threshold should be used instead. This is also consistent with experience with the masked moment method (e.g., Simon et al. 2001). An advantage of the dendrogram-based method is that a strong outlier does not affect neighboring pixels. Sometimes the masked spectral bins are shifted by one bin compared to the dendrogram-based masks (Fig. B.3). On the other hand, the masked moment method is more robust against small fluctuations around the threshold (see position 2 in Fig. B.3).

Appendix C Impact of different baseline subtractions on the [C II] optical depth estimates

As mentioned in Appendix A, the overall baseline quality for the [C II] spectra in IC 1396A is good enough so that the PCA has no significant effect. However, for the [¹³C II] hyperfine components, especially F = 2 – 1 is affected by a wing of the [C II] profile, and an additional baseline subtraction around each hyperfine component is required. Here we compare the derived [¹³C II] spectra and [C II] optical depths when using the second-order polynomial (Fig. 18) and first-order polynomial (linear; Fig. C.1). Except for clusters 1 and 3, the best fit second polynomial is very close to linear, so there is no difference in both cases (therefore not shown in Fig. C.1). Cluster 1 (left panel of Fig. C.1) has the largest difference depending on the order of the baseline among all clusters. The [¹³C II] spectrum at a velocity of between-14 km s^–1 and-10.5 km s^–1 indicates that the second- order polynomial (blue data in the figure) better characterizes the baseline. If we use a linear baseline, the derived optical depths are overall higher (around unity) than the calculation with the second-order polynomial baseline. However, the velocity dependence of the optical depth is more reasonable with the second-order baseline; the optical depth has a peak at a velocity close to the [¹³C II] peak velocity, instead of having the highest optical depth at the edge of the line profile as shown with orange data points in Fig. C.1 left panel. Therefore, we subtracted the second-order polynomial as the baseline.

Fig. C.1 Comparison of different orders of the baseline used to extract [13C II] hyperfine components for cluster 1 and cluster 3. Colored stepped data are the combined and scaled [13C II] spectra subtracted by the baseline with the first and second-order polynomial. Curves with corresponding colors are the fit to the [13C II] spectra using two Gaussian profiles, and derived optical depths of [C II] are shown as data points with error bars.

Appendix D Atmospheric model fit in the [O I] 145 µm observations

The [O I] 145 µm is located on a shoulder of a significant atmospheric ozone feature in the signal sideband. Figure D.1 shows an example of the atmospheric model fit in the "sky minus hot" spectrum, which is the difference between the OFF and hot-load measurement. Because the precise profile of the ozone feature depends on the pressure-temperature profile of the atmospheric layers, which is not precisely determined during the calibration process, the fit is often suboptimal as illustrated in this figure. This discrepancy results in a significant uncertainty in the absolute intensity calibration of the [O I] 145 µm. However, we focus in this paper on the line profile of [O I] 145 µm within the velocity range that is gray shaded in Fig. D.1. The atmospheric feature in this range has a smooth and monotonic profile, and thus the derived [O I] 145 µm line profile is reliable despite an imperfect fit of the atmospheric model.

Fig. D.1 Example of the "sky minus hot" spectra (the hot-load measurement subtracted from the OFF measurement and the atmospheric model fit in the [O I] 145 µm observations toward IC 1396A (March 6, 2020; LFAV pixel 1). The shaded velocity range is the plotted range in Fig. 6.

Appendix E Supplementary figures

Here we present supplementary figures. Figures E.1–E.4 are the position-velocity (p-v) diagrams of IC 1396A, B, and D. Figure E.2 is the [C I] 492 GHz integrated intensity map, shown in Okada et al. (2012), in order to relate with the positions discussed in this study. Figure E.6 shows the clusters detected in the 2D dendrogram analysis in IC 1396A, where the [C II] optical depth is calculated in Sect. 4.4. Figure E.5 show channel maps of ¹³CO(3–2) covering a larger area than the [C II] observations in this paper.

Fig. E.1 Overlay (center) of the integrated intensity of the [C II] (color) and 13CO(3–2) (white contours) of IC 1396A; the two cuts, (a) and (b) (red lines), mark where the p-v diagrams shown on the left and right were made. In the p-v diagrams the color map is [C II] (scale shown in the color bar) and the contours are 13CO(3–2) (with a spacing of 2.5 K). The spatial positions are vertically aligned between the p-v diagrams and the map. All images have a resolution of 16″.

Fig. E.2 Integrated intensity map of the [C I] 492 GHz (colors) in IC 1396A as presented in Okada et al. (2012), overlaid with contours of [C II] integrated intensity, both at 25″ resolution. Marks have the same meaning as in Fig. 5.

Fig. E.3 Same as Fig. E.1 but for IC 1396B and with white contours representing CO(3–2). In the p-v diagrams the color map is [C II] (scale shown in the color bar) and the contours are CO(3–2) (with a spacing of 10 K, starting at 5 K). The spatial positions are vertically aligned between the p-v diagrams and the drawn cuts. All images have a resolution of 25″.

Fig. E.4 Same as Fig. E.1 but for IC 1396D and with four cuts, (a), (b), (c), and (d). In the p-v diagrams the color map is [C II] (scale shown in the color bar) and contours are CO(3–2) (with a spacing of 10 K, starting at 5 K). The spatial positions are horizontally and vertically aligned between the p-v diagrams and the drawn cuts. All images have a resolution of 25″.

Fig. E.5 Channel maps of 13CO(3–2) (prior to resampling to the same resolution and mapping area as [C II]), overlaid with contours of [C II] in IC 1396E.

Fig. E.6 2D dendrogram clusters on the [C II] integrated intensity map in IC 1396A. We used the 25″ resolution map because the purpose of the dendrogram analysis is to average over a larger area for the [13C II] analysis.

References

All Tables

All Figures

Fig. 1 Overview of IC 1396. In this three-color composite image, blue is DSS-гed, green is the 9 µm map from the MIR all-sky survey by the Infrared Camera (IRC; Onaka et al. 2007; Ishihara et al. 2010) on board AKARI (Murakami et al. 2007), and red is the WIDE-L (140 µm) map from the FIR all-sky survey by the Far-Infrared Surveyor (FIS; Kawada et al. 2007; Doi et al. 2015) on board AKARI. The maps are not convolved to the common spatial resolution, and the coarsest spatial resolution is 88″ (Takita et al. 2015) for FIS (red), which is about a third of the size of the box around at IC 1396A. The fields observed with upGREAT in the [C II] line emission are overlaid as white boxes. The regions are labeled following the notation of Weikard et al. (1996). The position of the exciting source, HD 206267, is marked.

In the text

	Fig. 2 Dust temperature (left) and the continuum optical depth at 100 µm (right) fitted with a single temperature model to the AKARI/FIS data (see the main text). The white cross marks the position of HD 206267, and the white circle in the left panel indicates a 1 deg radius, where the dust temperature is 25 K from a purely geometric estimate (see the main text).
In the text

	Fig. 3 Channel maps of [C II] in IC 1396A (16″ resolution). The color scale is Tmb averaged over a velocity bin centered at the value noted in each panel with a width of 0.5 km s−1. Red crosses mark LkHα 349a (southeast) and c (northwest). HD 206267 is east of the BRC at (RA, Dec)(J 2000) = (21:38:57.619, 57:29:20.54).
In the text

	Fig. 4 Same as Fig. 3 but for the [O I] 63 µm emission (colors) overlaid with the [C II] 158 µm channel map contours (white with a spacing of 5 K) of IC 1396A.
In the text

Fig. 5 Integrated intensity maps (colors, 16″ resolution) of the different tracers (denoted on top of each panel) and overlaid with contours of the [C II] integrated intensity (solid white lines, 16″ resolution) in IС 1396A. Integrated intensities were calculated via a direct integration over the velocity range (−13, −4) km s−1. White areas indicate that a region was not observed, and black areas indicate that no line was detected (positions where the integrated intensities are below 3σ of the baseline are masked). Red crosses mark LkHα 349a (southeast) and с (northwest). Black crosses with numbers mark the positions of the spectra shown in Fig. 6. The bottom-right panel shows the schematic decomposition. Different colors indicate a difference in velocity (red represents the most redshifted component and blue is the most blueshifted), and the Roman numerals correspond to the description in the text. The black box border indicates the size of the images in the other panels.

In the text

Fig. 6 Spectra at selected positions in 1С 1396A, as marked in Fig. 5 (16″ resolution). The [C II] and [O I] spectral baselines are offset for a better visibility of the individual lines, and the vertical dotted lines are for guidance when comparing the peaks between different lines. [O I] 145 µm spectra are affected by the atmospheric ozone feature at velocities of < −10 km s−1 (see the main text).

In the text

Fig. 7 Peak velocity (upper panels) and line width at 10% of the peak intensity (lower panels) in units of km s−1 for 13CO(3−2), C18O(3−2), and [C II] (from left to right) in IС 1396A at 25″ resolution. The contours are the integrated intensity of [C II] at 25″ resolution. Note the different scale of the color wedge used for the [C II] line width, as the spatial distribution of the [C II] width would not be visible at the color scale used for the other two lines.

In the text

	Fig. 8 Channel maps of [С II] in 1С 1396В (16″ resolution). The bottom-right panel shows the Hα image from the archive of the INT (2.5 m Isaac Newton Telescope)/WFC (Wide Field Camera) Photometric Hα Survey of the Northern Galactic Plane (IPHAS; Barentsen et al. 2014) for comparison.
In the text

Fig. 9 Integrated intensity maps (colors, 16″ resolution) of the different tracers (denoted on the top of each panel) overlaid with contours of [C II] integrated intensity (solid white lines) in IС 1396B. White and black areas have the same meaning as in Fig. 5. Integrated intensities were calculated over the components detected by the dendrogram analysis (see the main text). Black crosses with numbers mark positions whose spectra are shown in Fig. 10. The right panel shows the schematic decomposition as in Fig. 5.

In the text

	Fig. 10 Spectra at the selected positions in IС 1396B marked in Fig. 9 (16″ resolution).
In the text

	Fig. 11 Same as Fig. 8 (16″ resolution) but for the [C II] in IС 1396D.
In the text

Fig. 12 Integrated intensity maps (colors, 16″ resolution) of the different tracers (denoted on the right of each panel) overlaid with contours of [C II] integrated intensity (white lines) in IC 1396D. White and black areas have the same meaning as in Fig. 5. Integrated intensities were calculated over the components detected by the dendrogram analysis (see the main text). Black crosses with numbers mark positions whose spectra are shown in Fig. 13. The bottom panel shows the schematic decomposition as in Fig. 5.

In the text

	Fig. 13 Spectra at the selected positions in IC 1396D marked in Fig. 12 (16″ resolution).
In the text

	Fig. 14 Same as Fig. 8 (16″ resolution) but for the [C II] in IC 1396E.
In the text

Fig. 15 Integrated intensity maps (colors, 16″ resolution) of the different tracers (denoted on the top of each panel) overlaid with contours of [C II] integrated intensity (white lines) in IC 1396E. White and black areas have the same meaning as in Fig. 5. The bottom right panel shows the schematic decomposition as in Fig. 5. Integrated intensities were calculated by direct integration over the velocity range of (−5, 8) km s−1, except for CO(3−2) with the range of (−10, 10) km s−1. Red pluses are the positions of outflow-driving sources: BIMAl-3 (southern ones, Beltrán et al. 2002) and the source C (northern one, Codella et al. 2001). Black crosses with numbers mark positions whose spectra are shown in Fig. 16.

In the text

	Fig. 16 Spectra at the selected positions in IC 1396E marked in Fig. 15 (16″ resolution).
In the text

Fig. 17 Red wing (left two panels) and blue wing (right two panels) line width at 10% of the peak intensity for 13CO(3−2) and [C II] in IC 1396E at 25″ resolution. Red pluses are the positions of the outflow-driving sources: BIMAl–3 (southern ones, Beltrán et al. 2002) and the source C (northern one, Codella et al. 2001) as in Fig. 15. The crosses on the [C II] panels are the positions of the H2 knots from Beltrán et al. (2009), with the size of the crosses corresponding to the H2 flux. Some of the H2 knots clusters are labeled according to Beltrán et al. (2009).

In the text

Fig. 18 Averaged spectra in each dendrogram cluster (Fig. E.6) in IC 1396A except cluster 2, which does not have a high enough S/N for the [13C II] analysis. Colored stepped data are the [C II] spectra (thick) and the combined and scaled [13C II] spectra (thin). The orange curve is the fit to the [13C II] spectra using two Gaussian profiles. Derived optical depths of [C II] are shown as black data points with error bars. See the main text for details.

In the text

	Fig. 19 ([C II]+[O I] 63 µm)/FIR integrated intensity ratio in IC 1396A at 25″ resolution. A jump in values along the shape of the rim is reflecting the [O I] 63 µm detection limit (see Fig. 5 for the area where [C II] and [O I] 63 µm are detected).
In the text

Fig. 20 Density plot of the correlation between [C II] integrated intensities and IRAC 8 µm (IC 1396A and E) or IRC 9 µm (IC 1396B and D) at 25″ resolution. The solid line is the relation derived in the Orion A region by Pabst et al. (2021), adapted to this study (see the main text). The dotted line indicates the slope of a linear relation (slope of unity in the log scale) with an arbitrary scaling (offset in the log scale).

In the text

	Fig. A.1 Difference of adjacent OFF measurements in one of the chopped OTF lines in the flight on March 6, 2020, for HFAV pixel 0. Each spectrum is smoothed to the velocity resolution of about 0.5 km s-1. The vertical axis is in Kelvin.
In the text

	Fig. A.2 Examples of derived principal components (HFAV pixels 0 and 2 for the flight on March 6, 2020).
In the text

	Fig. A.3 Examples of the baseline rms with and without the PCA (HFAV pixels 0 and 4 for the flight on March 6, 2020).
In the text

Fig. B.1 Integrated intensity maps of IC 1396D derived using different methods and thresholds at 16″ resolution. The upper three panels are from the 3D dendrogram method, with thresholds of 1σ, 2σ, and 3σ from left to right. The middle three panels are from the masked moment method, also with the threshold of 1σ, 2σ, and 3σ from left to right. The bottom panel is the direct integration of the spectra with the velocity range from –12 km s–1 to +10 km s–1. See the main text for details of the different methods. Black crosses with numbers mark positions where the masks along the spectral axis are shown in Fig. B.3.

In the text

	Fig. B.2 Integrated intensity maps of IC 1396A derived using the 3D dendrogram method, with thresholds of 1σ (left) and 3σ (middle), and the direct integration of the spectra in the velocity range –13 km s–1 to –4 km s–1. The color scale is the same for all panels.
In the text

	Fig. B.3 Masks along the spectral axis at the positions marked in Fig. B.1 derived using different methods for calculating the integrated intensity. The blue dotted line shows the 1σ level. Note that the mask is calculated by averaging the 3 × 3 spatial pixels, but the spectrum shown in each panel is the one at the central pixel.
In the text

Fig. C.1 Comparison of different orders of the baseline used to extract [13C II] hyperfine components for cluster 1 and cluster 3. Colored stepped data are the combined and scaled [13C II] spectra subtracted by the baseline with the first and second-order polynomial. Curves with corresponding colors are the fit to the [13C II] spectra using two Gaussian profiles, and derived optical depths of [C II] are shown as data points with error bars.

In the text

	Fig. D.1 Example of the "sky minus hot" spectra (the hot-load measurement subtracted from the OFF measurement and the atmospheric model fit in the [O I] 145 µm observations toward IC 1396A (March 6, 2020; LFAV pixel 1). The shaded velocity range is the plotted range in Fig. 6.
In the text

Fig. E.1 Overlay (center) of the integrated intensity of the [C II] (color) and 13CO(3–2) (white contours) of IC 1396A; the two cuts, (a) and (b) (red lines), mark where the p-v diagrams shown on the left and right were made. In the p-v diagrams the color map is [C II] (scale shown in the color bar) and the contours are 13CO(3–2) (with a spacing of 2.5 K). The spatial positions are vertically aligned between the p-v diagrams and the map. All images have a resolution of 16″.

In the text

	Fig. E.2 Integrated intensity map of the [C I] 492 GHz (colors) in IC 1396A as presented in Okada et al. (2012), overlaid with contours of [C II] integrated intensity, both at 25″ resolution. Marks have the same meaning as in Fig. 5.
In the text

	Fig. E.3 Same as Fig. E.1 but for IC 1396B and with white contours representing CO(3–2). In the p-v diagrams the color map is [C II] (scale shown in the color bar) and the contours are CO(3–2) (with a spacing of 10 K, starting at 5 K). The spatial positions are vertically aligned between the p-v diagrams and the drawn cuts. All images have a resolution of 25″.
In the text

	Fig. E.4 Same as Fig. E.1 but for IC 1396D and with four cuts, (a), (b), (c), and (d). In the p-v diagrams the color map is [C II] (scale shown in the color bar) and contours are CO(3–2) (with a spacing of 10 K, starting at 5 K). The spatial positions are horizontally and vertically aligned between the p-v diagrams and the drawn cuts. All images have a resolution of 25″.
In the text

	Fig. E.5 Channel maps of 13CO(3–2) (prior to resampling to the same resolution and mapping area as [C II]), overlaid with contours of [C II] in IC 1396E.
In the text

	Fig. E.6 2D dendrogram clusters on the [C II] integrated intensity map in IC 1396A. We used the 25″ resolution map because the purpose of the dendrogram analysis is to average over a larger area for the [13C II] analysis.
In the text