© ESO 2019

1 Introduction

The life cycle of the interstellar medium (ISM) is a central part of understanding star formation and galaxy evolution. The ISM not only influences present and future star formation but it also interacts with recently formed stars and continues evolving. A photon-dominated region (PDR; Tielens & Hollenbach 1985; Sternberg & Dalgarno 1995) is one of the places with such interactions, where the ultraviolet (UV) radiation from stars dominates the physical and chemical conditions of the surrounding ISM. In the Galaxy, detailed analyses of PDRs that take account of the source geometry have been performed (e.g., Andree-Labsch et al. 2017). On the other hand, some of the dominant emission lines such as [C II] 158 μm or [O I] 63 μm are often used to probe star formation in distant galaxies where individual star-forming regions are rarely resolved. In order to bridge our understanding in Galactic regions to high redshift work, studies in nearby galaxies with low (or different from solar) metallicity are necessary. The best local templates are the Large Magellanic Cloud (LMC) and the Small Magellanic Cloud (SMC) since they have significantly sub-solar metallicity (0.5Z_⊙ and 0.2Z_⊙ for carbon (Garnett 1999), where Z_⊙ is the solar abundance by Sofia et al. (2004)). They are also the nearest gas-rich systems (50 and 63 kpc, respectively, where we quote the median distance from the NASA/IPAC Extragalactic Database, NED), enabling excellent spatial resolution compared to other extragalactic targets.

Carbon is the second most abundant metal in the ISM. Dominant carbon-bearing species in PDRs are C⁺, C⁰, and CO. The [C II]/CO ratio in environments with different metallicities has been studied observationally and theoretically (Mochizuki et al. 1994; Poglitsch et al. 1995; Israel et al. 1996; Madden et al. 1997, 2011; Bolatto et al. 1999; Röllig et al. 2006), with nearly all studies finding that ratio increases with decreasing metallicity. This can be understood as a thicker C⁺ layer in regions where the UV radiation penetrates deeper due to lower dust extinction and, to a lesser extent, less self-shielding. By contrast, the transition H⁰ –H₂ does not shift as much as that of C⁺ –C⁰ or C⁺ –CO at lower metallicity because H₂ self-shielding is the dominant factor controlling the H⁰ –H₂ transition and is metallicity independent. An ISM phase where hydrogen is in molecular form but CO is photo-dissociated is called the CO-dark molecular gas. Understanding the mass and spatial distribution of this phase is important for estimates of the total gas mass using CO emission (Bolatto et al. 2013). However, it is not clear whether a similar explanation applies to the [C I]/CO ratio. Bolatto et al. (2000a) show that the gradient of [C I]/CO against metallicity is slightly shallower than the case A model by Bolatto et al. (1999), which assumes that the extent of the C⁰ layer in a clump scales linearly with the inverse of the metallicity (1/Z) and uses a mean column density of C⁺, C⁰, and CO over a clump size distribution. The latest KOSMA-τ PDR model (see Sect. 3) show different dependencies of C⁺ /CO and C⁰ /CO on metallicity; the difference in C⁰/CO between the Galactic model and LMC model is much smaller than for C⁺ /CO, for reasons that are not yet well-understood.

The velocity profile of spectral line emission represents the turbulent motions of emitting material. In our observing beam, different phases of the ISM coexist on spatial scales smaller than the beam and along the line-of-sight. When their motions are significantly different along the line-of-sight, we can distinguish them as separate velocity components in the line profile. This is especially useful for [C II] emission, which can originate from ionized gas, atomic gas, and molecular gas phases (Pérez-Beaupuits et al. 2012; Carlhoff 2013; Röllig et al. 2016; Fahrion et al. 2017). Okada et al. (2015, hereafter Paper I) show that the [C II] emission line profile is substantially wider than that of CO and [C I] in the N159 star-forming region in the LMC. The fraction of the [C II] integrated intensity that cannot be fitted by the CO-defined line profile varies across the map from 20 to 50%. The significant difference in the velocity profile of the [C II] emission observed in many sources is interpreted as an indication of a substantial fraction of the [C II] emitting material being accelerated relative to the quiescent material, for example it undergoes ablation (Dedes et al. 2010; Mookerjea et al. 2012; Okada et al. 2012; Schneider et al. 2012; Simon et al. 2012; Pilleri et al. 2012), or photoevaporated (Sandell et al. 2015). Although an order of a few km s⁻¹ displacement is consistent with modeled photoevaporation in globules (Lefloch & Lazareff 1994), only a few studies to compare the observed line profiles of different emission lines with simulation incorporating full PDR chemistry have been conducted (Bisbas et al. 2018). An interesting question is whether the velocity profile difference is observed only in [C II] or also in [O I], because the [O I] intensities constrain physical conditions of the ablated material. Since 2014 the German REceiver for Astronomy at Terahertz Frequencies (GREAT) on board the Stratospheric Observatory for Infrared Astronomy (SOFIA) enables observations of velocity-resolved [O I] 63 μm emission. These observations often show strong absorption features (Ossenkopf et al. 2015; Leurini et al. 2015). Since the [O I] 145 μm emission is typically optically thin, the velocity profile of this line gives unambiguous information about the gas dynamics with minimal opacity effects on the line profile.

As mentioned above, the LMC provides a good low metallicity template to bridge Galactic and high redshift observations, since it is close enough that star-forming regions are spatially well-resolved, and independent resolution elements do not show significant blending of different velocity components. In the LMC, 30 Doradus (30 Dor) and other high-mass star-forming regions in the molecular ridge south of 30 Dor (see Fig. 1), provide a good opportunity to study star formation and the ISM in a low metallicity environment.

30 Dor is a well-studied, prominent star-forming region excited by the super star cluster R136, where 117 massive stars provide a total ionizing fluxof 4 × 10⁵¹ Ly photons s⁻¹ (Crowther & Dessart 1998). Many observations suggest that the ISM must be highly clumped; Poglitsch et al. (1995) show that the integrated intensity maps of the [C II] and CO(1-0) emission are spatially strongly correlated, suggesting that the molecular and the photodissociated gas appear coextensive over ~30 pc. Kawada et al. (2011) show that the [O III] emission is widely extended spanning more than 150 pc from R136. Indebetouw et al. (2013) resolved clumps with the Atacama Large millimeter/submillimeter Array (ALMA) and show that most of the flux in the CO map is contained within an area of 5–10%. They also show that relative velocities of individual clumps, typically with a line width of 2–3 km s⁻¹, dominate the velocity structure, with an additional overall east–west gradient. Yeh et al. (2015) suggest that the H₂ rotational–vibrational transition in30 Dor is due to fluorescence without evidence of shock excitation. Chevance et al. (2016) model FIR fine-structure lines and determine the three-dimensional structure of the gas by comparing the modeled UV radiation field strength with the emitted radiation by stars.

N158 is located south of 30 Dor with an elongated H II region (Henize 1956). The southern part excited by the OB association L101 (Lucke & Hodge 1970) is identified as N158C or NGC2074. LH101 contains two populations of massive stars (≤2 Myr and 3–6 Myr; Testor & Niemela 1998), creating elongated H II regions in a north–south direction with an extension to the west seen in Hα and the optical [O III] emission (Fleener et al. 2010; Galametz et al. 2013). Fleener et al. (2010) found that the majority of the YSOs and central stars of ultracompact H II regions in NGC2074 are not earlier than the late O-type, whereas there are several evolved early O-type stars, indicating that their formation may have started at a similar time, a few 10⁵ years ago.

N160 is the hottest among the star-forming regions in the molecular ridge in terms of the dust temperature (Galametz et al. 2013) and the CO excitation (Heikkilä et al. 1999; Bolatto et al. 2000b). N160A is the brightest cloud in the radio continuum which has three compact H II regions (Heydari-Malayeri & Testor 1986; Indebetouw et al. 2004). In the southwest of N160A, there is another radio continuum source N160D, where the ionized gas has a shell-like structure with a WR star as well as several OB stars at the center (Heydari-Malayeri & Testor 1986; Martín-Hernández et al. 2005; Fariña et al. 2009), and two compact H II regions are located at the eastern edge of the shell.

N159 is the most southern region among the four regions that we study here. Three CO cores are identified as N159 W, E, and S by Johansson et al. (1998). Recent ALMA observations reveal detailed filamentary structures and converging flows at N159 W and E (Fukui et al. 2015; Saigo et al. 2017). Lee et al. (2016) report that the PDR properties that produce a good fit to the [O I], [C II], and [C I] emission predict too low CO intensities compared to the observations, suggesting that the CO emission may be excited by low-velocity C-type shocks.

A global age gradient from 30 Dor southwards to N159 has been suggested by various studies. Nakajima et al. (2005) found an age gradient of the YSOs from N160 to N159 toward the tip of N159 S and suggest a propagation of triggered star formation from north to south along the molecular ridge. Israel et al. (1996) measure a lower CO luminosity compared to [C II] and far-infrared (FIR) in 30 Dor compared to N159 and N160, presumably reflecting that 30 Dor is the most evolved region. On the other hand, since both N160 and N159 contain young objects as well as evolved objects, Fariña et al. (2009) suggest a common time for the origin of recent star formation in the N159–N160 complex as a whole, while sequential star formation at different rates is probably present in several subregions. The analysis by Gordon et al. (2017) indicated that N159 is indeed younger than N160, but these authors did not conclude whether there is an overall north–south age gradient. de Boer et al. (1998) proposed a scenario in which star formation is triggered at the leading edge due to the bow-shock of the LMC with its motion through the halo of the Milky Way; in this scenario, an age gradient is expected because of the rotation of LMC.

We performed velocity-resolved mapping observations of CO, ¹³CO, [C I] ³P₁–³P₀ and ³P₂–³P₁, and [C II] in N159, N160, N158, and 30Dor, and velocity-resolved pointed observations of [O I] 145 and 63 μm at selected positions in N159 and 30 Dor. In Sect. 2 we present the data that we use in this study. A brief summary of the KOSMA-τ PDR model is given in Sect. 3. Section 4 presents all the results and discussion. We start by presenting a line profile analysis of the [O I] emission (Sect. 4.1) and other emission lines (Sect. 4.2). Next we estimate the contribution of the atomic gas (Sect. 4.3) and ionized gas (Sect. 4.4) to the [C II] emission. An overview of the line profile analysis is given in Sect. 4.5. In Sect. 4.6 we present details of the spatial distribution of detected emission lines in individual regions, and discuss the column density of CO, C⁰, and C⁺ in the low metallicity LMC environments in Sect. 4.7. The results of the PDR model fitting are given in Sect. 4.8. Our main findings are summarized in Sect. 5.

Fig. 1 Center panel: IRAC 8 μm map with boxes showing the areas we observed. Four inset color maps with contours: integrated intensity map of [C II] (color) overlaid with contours of CO(4-3) integrated intensity at 16′′ resolution. The red lines outline the area observed in [C II]. Contour levels are at intervals of 5 K km s−1 for N158, 10 K km s−1 for N160 and 30 Dor, and 15 K km s−1 for N159. The unit of color bar is K km s−1. In the maps of N159 and 30 Dor, the positions of the [O I] observations are marked following the numbering in Table 2. In 30 Dor, the position of R136 is shown as a red circle, and the two cuts to make position-velocity diagrams (Fig. 11) are shown as black lines. In N160, two H II regions N160A1 and A2 (Martín-Hernández et al. 2005) are marked.

2 Observations and data reduction

2.1 [C_II] 158 μm observations with GREAT on board SOFIA

[C II] mapping observations were made in four regions using the German REceiver for Astronomy at Terahertz Frequencies (GREAT¹; Heyminck et al. 2012)on board the Stratospheric Observatory for Infrared Astronomy (SOFIA; Young et al. 2012). We observed N159 and 30 Dor in July 2013, as part of the open time and guaranteed time in cycle 1 observations, and N160 and N158 in July 2015, as part of the guaranteed time in cycle 3 observations. The L2 channel was tuned to [C II], and the backend was XFFTS with 2.5 GHz bandwidth in 2013, FFTS4G with 4 GHz bandwidth in 2015. The [N II] 205 μm was observed in parallel only for N159. All mapping observations were made in total-power on-the-fly (OTF) mode with 6′′ step size. Other parameters of observations are summarized in Table 1. See also Paper I for details of the observations in N159. For 30 Dor, we found that the primary OFF position is contaminated in the [C II] emission, and we executed total-power pointed observations of this contaminated OFF position against another more distant OFF position (see Table 1). We averaged 12 spectra with an integration time of 15 s each to obtain the OFF spectrum, and added it to the mapping spectra to correct the OFF contamination. The uncertainty of the OFF spectra was estimated by the standard deviation of the 12 OFF spectra, which is 0.2 K as a median of all velocity bins between 220 and 280 km s⁻¹, and 0.7 K as a maximum.

The data were calibrated by the standard GREAT pipeline (Guan et al. 2012), which converts the observed counts to the main beam temperature(T_mb). One improvement since Paper I is that we interpolated the internal calibration (load) measurement and the OFF measurement in time to better account for the gain drift. We compared this new calibration in N159 with the one used in Paper I and confirmed that the resulting spectra do not have a significant change but we get better signal-to-noise ratio (S/N). We then subtracted second-order polynomial baselines and spectrally resample to 1 km s⁻¹ channel width. We convolved the [C II] map to different spatial resolutions as required for further analysis. For the discussion of the spatial distribution, we convolved it to 16′′ resolution (Sect. 4.6). Since we use a simple Gaussian kernel for the convolution, where we define the width of the Gaussian kernel by the difference between the square of the goal beam size and the square of the original beam size, we cannot avoid losing the spatial resolution slightly compared to the original resolution. To compare the line profiles with the [O I] emissions, we also use the 16′′ resolution. For the analysis of the line profile and the PDR modeling, we use a spatial resolution of 30′′ for all regions in order to make use of CO(2-1) data (see Sect. 2.3). This also gives maps with a better S/N. When comparing the [C II] emission with the H I emission, we convolved the [C II] map to 1′ resolution.

2.2 [O_I] 145 and 63 μm observations with upGREAT on board SOFIA

We observed[O I] 145 and 63 μm at selected positions in N159 and 30 Dor with the upGREAT (Risacher et al. 2016) onboard SOFIA in July 2017, as part of the guaranteed time in cycle 5 observations. The two polarizations of the Low Frequency Array (LFAH and LFAV) were tuned to the [O I] 145 μm line, and the High Frequency Array (HFA) was tuned to the [O I] 63 μm line, and they were operated in parallel. Each array has 7 pixels in a hexagonal configuration and the beam size of the [O I] 145 and 63 μm is 13′′ and 6.3′′, respectively. We selected two positions in each of N159 and 30 Dor, where the line profiles of the CO and [C II] emission are significantly different, and executed pointed single-phase chopped observations. Since the [O I] 63 μm line for the source velocity of the LMC (~250 km s⁻¹) is located at the edge of the atmospheric window and at a difficult tuning range for the quantum cascade laser (QCL) used in the HFA, we optimized the observations to the [O I] 145 μm emission. We chose the rotation of the arrays so that some of the LFA outer pixels observed other positions of interest in parallel. We focused on the HFA at only one position in the 30 Dor region, where we observed four point rasters (central position plus 4′′ east, west, and north) with a step of 4′′ in order to fill a LFA beam, so that we can compare the line profile of the [O I] 63 μm with the [O I] 145 μm and [C II].

The central pixels of LFAH and LFAV are aligned within 2′′, which is small enough compared to the LFA beam size. The central pixel of HFA is about 3′′ away from the LFA central pixels, and outer pixels of HFA and LFA are located at different positions because the size of the array is scaled by the beam size.

Before applying the standard calibration with the GREAT pipeline, we made a frequency correction and a gain correction. The frequency correction is always needed for the HFA, since the QCL locks only to a discrete frequency. We correct the frequency scale scan by scan by using the narrow atmospheric [O I] line. In addition, we had a synthesizer reference problem for the LFA arrays during flight, which results in an imprecise LO frequency and requires a similar frequency correction for the [O I] 145 μm observations. In the observed LFA band, only a wider atmospheric feature is available, thus the frequency uncertainty of the [O I] 145 μm in the present observations is ~5 MHz, which corresponds to 0.7 km s⁻¹. We also applied a gain correction using the total raw count (see Appendix A for details). Then the data were calibrated by the standard GREAT pipeline, as for the mapping observations. In Table 2, we list the positions where at least one of the [O I] lines was detected.

2.3 CO and [C_I] observations with APEX

We carried out complementary observations with the Atacama Pathfinder Experiment (APEX²; Güsten et al. 2006)of CO and [C I] emission lines. The list of observed lines is shown in Table 3. CO(2-1) and its isotopes were observed with The Swedish Heterodyne Facility Instrument (SHeFI; Vassilev et al. 2008) in 2012–2014 (project ID: M-081.F-0034-2008, M-093.F-0036-2014). For N159, only pointed observations toward 12 positions for CO(2-1) and one position for ¹³CO(2-1) were made.CO(3-2) and its isotopes, CO(4-3), and [C I] ³P₁–³P₀ emissions were observed with the MPIfR heterodyne receiver FLASH⁺ (Klein et al. 2014) in 2012–2016 (project ID: M-081.F-0034-2008, M-093.F-0036-2014, M-097.F-0002-2016), and CO(6-5) and [C I] ³P₂–³P₁ with the CarbonHeterodyne Array of the MPIfR (CHAMP⁺; Kasemann et al. 2006)in 2012–2014 (project ID: M-081.F-0034-2008, M-089.F-0008-2012, M-093.F-0036-2014). All observations were made in the OTF mode except for the CO(2-1) and ¹³CO(2-1) in N159. Weused the same OFF position as in the GREAT observations. We reduced the spectra in the same way as GREAT; resampled to 1 km s⁻¹ and spatial resolution of16′′ or 30′′, depending on the analysis.

2.4 H_I data from ATCA and Parkes

To investigate the contribution of H I gas to the [C II] emission, we use the combined H I survey data by Australia Telescope Compact Array (ATCA) and Parkes telescope (Kim et al. 2003). The angular resolution is 1′, corresponding to a linear scale of ~15 pc, and the spectral resolution is 1.65 km s⁻¹. The cube covers heliocentric velocities between − 33 and 627 km s⁻¹, with emission from the LMC arising between 190 and 386 km s⁻¹. The rms of the brightness fluctuations in the final cube is 2.4 K per channel. Integrating over three channels, this corresponds to a column density sensitivity limit of ~ 2.2 × 10¹⁹ cm⁻².

2.5 [O_I] observations with PACS

We also use the [O I] 63 and 145 μm data, observed with the Photoconductor Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) on the Herschel Space Observatory³ (Pilbratt et al. 2010), in order to better constrain the PDR modeling. All these data are presented in Cormier et al. (2015). A detailed analysis of these lines in 30 Dor has been presented by Chevance et al. (2016), and in N159 W by Lee et al. (2016). In order to avoid regridding more than once, we downloaded the level 2 data in the original 5 × 5 spaxels from the HerschelScience Archive and directly extracted the spectra. We assume that the brightness of individual spaxel is uniformly distributed in a 9.4′′ × 9.4′′ square. We defined a geometrical area which corresponds to each pixel in our GREAT and APEX dataset as a circle with diameter of 30′′, and extracted spectra as an average within this circle. As mentioned in Cormier et al. (2015), some observations in 30 Dor have contaminated OFF spectra, where we re-ran the last step of the pipeline to skip the procedure of the OFF subtraction and simply used only ON spectra because we are not interested in the continuum level.

2.6 FIR continuum emission with Spitzer and Herschel

We also use the Spitzer MIPS 70 μm continuum from the SAGE project (Meixner et al. 2006) and PACS 100 μm, 160 μm, SPIRE 250 μm and 350 μm by the HERITAGE project (Meixner et al. 2013) to construct the dust spectral energy distribution (SED) to be fitted with the PDR model. As for the [O I] line emission, we took a geometrical average within the circle of a 30′′ diameter to make the same grid as the GREAT and APEX dataset. We did not use SPIRE 500 μm because its spatial resolution is worse than 30′′.

3 PDR model

We use the KOSMA-τ PDR model (Störzer et al. 1996; Röllig et al. 2006, 2013) to investigate the effect of metallicity on the PDR structure and chemistry and to investigate physical properties of the four regions. As presented in Röllig et al. (2013), the revised KOSMA-τ model treats the dust-related physics consistently and computes the dust continuum emission together with the line emission of chemical species. This has been now extended to include all dust size distributions tabulated in Weingartner & Draine (2001), including the LMC and SMC dust models. By applying the elemental abundances for the LMC (X(C) = 7.94 × 10⁻⁵, X(O) = 2.51 × 10⁻⁴; Garnett 1999), the model can be compared with the far-infrared continuum SED and the line intensities in the LMC.

The KOSMA-τ model assumes either a single clump filling the beam (hereafter “non-clumpy model”) or an ensemble of clumps with a power law mass spectrum dN∕dM ∝ M^−α where α = 1.8 (hereafter “clumpy model”). Each clump follows the mass-size relation $M_{\text{clump}}\propto R_{\text{clump}}^{\gamma }$ with γ = 2.3, and has a power-law radial density profile outside of the core, where the density is constant (Störzer et al. 1996; Cubick et al. 2008). The line and continuum intensities are calculated by summing the contribution from each clump in a model grid of three independent parameters; total mass (m), mean gas density (n), and far-ultraviolet (FUV; hν = 6–13.6 eV) flux (χ) in units of the Draine field (2.7 × 10⁻⁶ W m⁻²).

Fig. 2 [O I] 145 and 63 μm spectra in N159 at the positions marked in Fig. 1. The CO(4-3) and [C II] spectra were extracted from the maps with a resolution of 16′′, while the [O I] observations are shown with the original spatial resolution (13′′ for the [O I] 145μm and 6.3′′ for the [O I] 63 μm).

4 Results and discussion

The properties of the detected emission lines in each region are summarized in Table 3. Figure 1 shows the integrated intensity maps of [C II] and CO(4-3) with 16′′ angular resolution, together with an overview of IRAC 8 μm emission. In Appendix B, integrated intensity maps of the detected emissions at 30′′ resolution and spectra of selected positions are shown for N160, N158, and 30 Dor. For N159 data we refer to Paper I. The [O I] 145 μm and 63 μm lines at selected positions (Table 2, and marked in Fig. 1) are shown in Figs. 2 and 3.

4.1 The velocity resolved [O_I] emissions in the LMC

We observed the velocity resolved [O I] 63 and 145 μm emission in N159 and 30 Dor for the first time. In N159, the central pixels of LFAH and LFAV detected the [O I] 145 μm at two observed positions. We averaged the LFAH and LFAV spectra using a weight based on the baseline noise. In addition, we detected the line with the LFAV pixel 3 at the first observed position. We used only the LFAV data for this position because the LFAH pixel 3 performed much worse than the LFAV pixel 3. The [O I] 145 μm spectra at these three positions are shown in Fig. 2 and the positions are marked in Fig. 1. The 63 μm emission was detected by the HFAV center pixel and is also shown in Fig. 2. A direct comparison of the line profile is difficult because of the different beam sizes. The line profile of the [O I] 145 μm has a mixture of the characteristics of the CO(4-3) and [C II] profiles. At position 1 in Fig. 2, the blue wing of the [O I] 145 μm line is identical to that of [C II], while the broad red wing and a second velocity component at ~ 240 km s⁻¹ is missing in the [O I] 145 μm. On the other hand, at position 3, the [O I] 145 μm profile follows thatof CO(4-3), which is narrower than the [C II] line. The [O I] 145 μm profile is difficult to interpret at position 2 due to low S/N, but the [O I] 145 μm peak velocity may be between that of [C II] and CO.

In 30 Dor, we detected [O I] emission only at the first pointed positions, where the four point raster observations were carried out (see Sect. 2.2). The six positions listed in Table 2are the positions of different pixels. In Fig. 3, the four raster positions were averaged for the [O I] 63 μm, so that we can compare its line profile with other emission lines. We took a simple mean instead of a weighted mean in order to weight different spatial positions equally. The averaged [O I] 63 μm spectra have line profiles that tend to be more similar to the [C II] profile in terms of the overall line width. At position 1, a dip around ~ 250 km s⁻¹ may be self-absorption, because it corresponds to the peak of CO and [O I] 145 μm. On the other hand, the profiles at positions 4 and 5 seem to show two velocity components around ~ 242 km s⁻¹ and ~ 250 km s⁻¹, where CO(4-3) has only a very weak emission for the former component, and [C II] emission fills the velocity between these two components, although it is not excluded that all of CO(4-3), [C II], and [O I] 63 μm have an absorption feature. The [O I] 145 μm emission behaves as in N159: a mixture of the characteristics of the CO(4-3) profile and [C II] profile. At position 1 it is closer to the CO(4-3) profile, at position 2 to the [C II] profile, and at position 3, it has the blue profile of CO(4-3) and red profile of [C II]. Figure 4 shows that the [O I] 63 μm spectra at individual raster positions are different, indicating that different velocity components often originate from spatially separated components. This is consistent with the ALMA observations by Indebetouw et al. (2013), showing that the relative velocities of individual clumps dominate the velocity structures.

We compare the [O I] integrated intensities measured by upGREAT to those by PACS. We extracted PACS spectra following the method described in Sect. 2.5 but using a circle with diameter of 13′′. For [O I] 145 μm, the integrated intensities match within 20% except for the two faintest positions (N159 position 2 and 30 Dor position 3), where they match within 30%. For [O I] 63 μm, we make comparisons only at positions in 30 Dor where the upGREAT observations consist of the four point raster. For the bright positions (30 Dor position 1 and 5), the discrepancy between PACS and upGREAT (raster positions averaged) is within 20%, which is also consistent with the uncertainties caused by the gain correction (Appendix A). In view of the difficult to assess beam coupling effects for our undersampled data set, and the different beam width of the observations, we consider theagreement encouraging.

Fig. 3 Same as Fig. 2 but for 30 Dor. Here the [O I] 63 μm spectra are the average of the four raster positions (see text).

Fig. 4 [O I] 63 μm spectra of four raster positions and their average (black histogram) in 30 Dor positions 1, 4, 5, and 6.

4.2 Line profile analysis

In Paper I, we show that a combination of Gaussian profiles defined by CO(3-2) typically well reproduces the line profiles of CO, ¹³CO, and [C I] emissions. [C II], by contrast, has a very different line profile that is typically wider than the CO profile, but is often not symmetric and is complex in N159. The other three regions studied here also show the same trend. Appendix B shows spectra of selected positions in each region.

To quantitatively investigate the velocity profiles, we fit one of the CO spectra with Gaussians as a reference, and applied the fitted centervelocity and width in the fit to the other emission lines, fitting only their amplitudes. For the mapping observations, we used the CO(3-2) as a reference emission line and performed fitting at each grid position in a 30′′ resolution map (a detailed description is given in the Paper I). For the pointed observations of [O I] 63 and 145 μm, we used the 16′′ resolution data and used CO(4-3) as a reference line. We fit the reference CO spectra with two Gaussians, except for eye-selected positions where there are clearly three velocity components identified. In the following, we refer to the fitted center velocity and width of those Gaussians as “the CO-defined line profile”, and investigate how well the (amplitude-scaled) CO-defined line profile describes the other emission lines.

We compare the sum of the fitted Gaussians with the integrated intensity for each emission line. For all detected CO and ¹³CO emissions in four regions, the integrated intensity is well reproduced by the sum of Gaussians. As shown in the Paper I, the ¹³CO(3-2) and [C I]³P₁–³P₀ lines in N159 show that the sum of Gaussians yields larger values than the integrated intensity where the emission is strong, which is interpreted as a broadening of the ¹²CO lines by the optical depth effect. The ¹³CO(3-2) in N160 shows the same trend, but it is not the case for 30 Dor and N158. For the [C I]³P₁–³P₀ emission, we do not see the trend other than in N159 because of the large scatter due to the lower S/N in other regions.

What is common across all four regions is that the sum of the Gaussians of the [C II] emission is on average 30% lower than its integrated intensity, which suggests that the CO-defined line profile can reproduce only a fraction of the [C II] emission. As in N159, the [C II] profile is closer to the CO-defined profile toward the CO peaks, and the fraction of the [C II] emission that cannot be fit by the CO-defined Gaussians increases toward the regions between CO peaks, up to 60% in some positions in 30 Dor. Possible contributions from the atomic gas and the ionized gas are discussed in the next subsections.

As shown in Sect. 4.1, the [O I] emission has a profile between CO(4-3) and [C II]. The Gaussian fit result confirms this conclusion: the fraction of the [O I] emission that cannot be fit by the CO-defined Gaussians are clearly correlated to its fraction for the [C II] emission, and somewhat smaller than the [C II] emission. In N159, this fraction for the [O I] 145 μm is about 25% at position 1 and 3, whereas the fraction for the [C II] is 35–45%. In 30Dor, the fraction for [O I] 145 μm is 0–20%, [O I] 63 μm is 0–30%, and [C II] is 15–45%.

Fig. 5 [C II] and H I spectra at a [C II] blob in N159 (05:39:45.9, −69:44:32.2, upper panel)and 30 Dor (05:38:31.6, −69:04:50.0, lower panel) at 1′ spatial resolution.

4.3 Contribution of the atomic gas to the [C_II] emission

In order to estimate the contribution of the atomic gas to the [C II] emission, we compared the [C II] spectra with the H I spectra at 1′ resolution. In general, the H I emission has much broader line profiles than those of [C II] in N159, N160, N158, and regions not associated with CO peaks in 30 Dor. Towards CO peaks in 30 Dor, the H I shows absorption features. Figure 5 shows the [C II] and H I spectra at two positions in N159 and 30 Dor, where the [C II] emission dominates over the CO and [C I] emission. At the position in 30 Dor, the [C II] velocity component at ~ 255 km s⁻¹ has a counterpart in the CO emission, but the component at ~ 265 km s⁻¹ is only seen in [C II] and not in CO, and the H I emission shows an absorption feature at this velocity. This is consistent with a picture of the H I self-absorption by cold molecular clouds with a mixture of atomic hydrogen (Li & Goldsmith 2003; Gibson et al. 2005; Kavars et al. 2005; Klaassen et al. 2005; Tang et al. 2016), although we do not resolve individual clouds in the LMC and the absorption feature is not as narrow as specified in their studies.

We estimate the H I gas contribution to the [C II] emission using the line profile. We scaled the H I spectra to fit the wing of the [C II] emission and attributed it as a contribution of the H I gas to the [C II] emission. In most regions in N158, N159 and N160, the estimated fraction of the H I gas contribution to the [C II] emission is 15% or less, except for N159 E and the western edge of N160, where the fraction is higher. In 30 Dor, we could not fit the wing well because the velocity profile is more complex than in other regions. Since the H I profile is much broader than the [C II] line profile, subtracting this contribution does not make the [C II] line profile significantly narrower, and it still cannot be explained by the CO-defined line profile.

We then estimate the thermal pressure (or density) of the H I gas that is required to emit the fitted wings of the [C II] intensity assuming a temperature (Pineda et al. 2017). Using the equation of optically thin subthermal [C II] emission from Goldsmith et al. (2012) and the gas temperature of 100 K, the derived thermal pressure is 4 × 10³–10⁵ K cm⁻³ in three regions except for 30 Dor. They are higher than the standard Galactic ISM (e.g., the cold neutral medium with T_k = 100 K and n = 30 cm⁻³; Draine 2011), and consistent with previous studies suggesting higher thermal pressure in the LMC (Welty et al. 2016; Pineda et al. 2017). Regions with higher pressure ≥ 3 × 10⁴ K cm⁻³ occur around N159 E, the [C II] peak in N159 W, the N160 peak and its southern and western part. This is consistent with Welty et al. (2016), who find higher pressure toward complex regions including H II regions, molecular clumps and supernova remnants, where the pressures may be enhanced by energetic activities. The derived pressure is not sensitive to the assumed temperature for T > 80 K (Pineda et al. 2017).

With the above method, we may underestimate the [C II] integrated intensity that can be attributed to the atomic phase because we assume that the H I line profile is optically thin. Braun et al. (2009) and Braun (2012) show that flat-topped H I emission profiles can be modeled by an opacity effect, and the opacity correction of the H I column density can be an order of magnitude in our observed regions. In the above method, we do not use the H I column density itself but use the H I profile to fit the [C II] emission. A flat-topped profile and steeper line wings of the H I emission lead to an underestimate of the corresponding [C II] emission, if the [C II] emission is not as optically thick as the H I emission.

4.4 Contribution of the ionized gas to the [C_II] emission

As shown for the case of N159 in Paper I, the velocity-resolved [N II] spectra are a useful tool to investigate the contribution of the ionized gas to the [C II] emission. However, the [N II] observations were made only in N159. Here we examine a possible ionized gas contribution through the spatial distribution of the residual of the [C II] emission after subtracting the CO-defined profile. In contrast to Galactic-scale analysis (e.g., Pineda et al. 2014), the contribution of the ionized gas is expected to vary locally in spatially resolved massive star-forming regions. As described in Sect. 4.2, we fit the [CII] spectra with Gaussians defined by the CO profile. We refer to the residual at greater velocities than the peak of the Gaussian (if more than one Gaussians are used in the fit, the peak with the largest velocity) as the red-side residual, and the residual at smaller velocities than the bluest Gaussian peak as the blue-side residual.

The spatial distributions of the integrated intensity of the red-side residual and blue-side residual are clearly different in all regions. The left panel of Fig. 6 shows that the integrated intensity of the red-side residual in N159 correlates with the Hα emission, which confirms the result of Paper I, meaning that the ionized gas makes some contribution to red wings of the [C II] emission. On the other hand, the blue residual is strong in N159 E and the [C II] blob (see Paper I). In N159 E, ALMA observations resolved a complex velocity structure consist of a few colliding filaments (Saigo et al. 2017). The [C II] blue residual in this study may be a gas component at peculiar velocities in this kinematically complex region with physical conditions that emits dominantly in [C II], rather than CO.

In N160, neither the red-side nor blue-side residual shows a correlation with the Hα emission. The red-side residual is strong towards the extended [C II] emission in the northwest (Fig. 1), which may be the ionized gas contribution because the radio continuum emission extends weakly there (Mills & Turtle 1984; Israel et al. 1996). In other regions, the spatial distribution of the red-side and blue-side residuals do not give a strong indication of an ionized gas contribution. In 30 Dor, the region ~ 1′ northeast of the [C II] peak at 30 Dor-10 shows a strong blue-side residual, which may be due to a gas stream (Sect. 4.6).

Fig. 6 Integrated intensity maps of the red-side (left panel) and blue-side (right panel) residual of the [C II] emission after subtracting the CO-defined line profile (see text) in N159 (color, at 30′′ resolution) overlaid with contours of Hα (Chen et al. 2010).

4.5 Discussion of the line profile shapes

The pointed observations of [O I] 145 μm and 63 μm at N159 and 30 Dor suggest that neither the CO nor the [C II] line profiles simply represent the [O I] profile. Since 16′′ corresponds to ~ 4 pc at the distance of LMC, the observed line profile must be a composite of different cloud components. This interpretation is consistent with the fact that the [O I] 63 μm profile varies among the four raster positions with a spacing of 4′′ (Fig. 4), and that individual clumps detected by ALMA in 30 Dor typically have linewidths of 2–3 km s⁻¹ (Indebetouw et al. 2013), which is narrower than the results presented here within a 30′′ beam. We propose the following interpretation for the observed different line profiles among CO, [C II] and [O I] in the LMC. In addition to a possible ionized gas contribution to the [C II] emission, we consider that our beam includes several PDR components that are spatially separated and/or are in different physical phases, and each component contributes to a certain velocity range in the observed line profiles depending on their dynamics. In this case, line ratios of individual velocity components can vary depending on the physical properties of the corresponding gas component. When a PDR component emits dominantly the [C II] line and the intensity of the CO emission is below our detection limit (for example, a low density case), we detect it as a [C II] component which cannot be reproduced by the CO-defined velocity profile. As far as the contribution of the H I gas and ionized gas to the [C II] emission is excluded, we could call it CO-dark molecular gas. We should, however, keep in mind that the derived fraction depends on the detection limit of the CO emission.

4.6 Spatial distributions in individual regions

Figure 1 shows the integrated intensity map of CO(4-3) and [C II] in the observed four regions, and Figs. 7–9 show the channel maps (N159 is not shown). As discussed in Paper I, the spatial distribution of the integrated intensity of [C II] and CO(4-3) emission is not the same in N159. This is also the case for N158 to some degree, whereas in 30 Dor and N160, they are more similar. The integrated intensity maps of all detected emissions converted to 30′′ are shown in Appendix B. In the following, we describe the detailed morphology in each region except for N159, for which we refer the reader to Paper I.

4.6.1 N160

The peak of the [C II] integrated emission is located at N160A, which is the brightest region in the radio continuum and Hα emission (Heydari-Malayeri & Testor 1986; Martín-Hernández et al. 2005). N160A consist of subregions A1, A2, and another cluster in between (Heydari-Malayeri et al. 2002). Both CO(4-3) and [C II] emission show an elongated shape tracing those subregions (Fig. 1). In radio continuum emission, another H II region is identified as N160D southwest of N160A, with the ionized gas having a shell-like structure with a WR star as well as several OB stars at the center (Heydari-Malayeri & Testor 1986; Martín-Hernández et al. 2005; Fariña et al. 2009). Two compact H II regions are located at the eastern edge of the shell, where our CO(4-3) has the strongest peak and the [C II] has a peak as well. It corresponds to N160-4 in Johansson et al. (1998). The shell structure can be traced in the [C II] map in the bottom panels of Figs. 7 and 10, where positions of OB stars and compact H II regions are marked. OB stars in Fariña et al. (2009) are distributed not only in N160A and D, but continue in southeast direction of the observed region toward another radio continuum peak (Mills & Turtle 1984). The [C II] emission extends northwards from N160A, which can be also seen in the [C II] map by Israel et al. (1996). There is one late O-star to the north of this extended distribution (Fariña et al. 2009). To the east of N160A, we detected two clouds in CO(4-3), where the [C II] emission is very weak (N160-1 and N160-3 in Johansson et al. 1998). Figure B.1 shows that these clouds are bright in [C I] ³P₁–³P₀ and low-J CO, while CO(6-5) is not detected, indicating a cold and quiescent environment. However, they must have exciting UV sources because the [C II] emission is weakly detected and the IRAC 8 μm map shows diffuse structure as well as cores.

The [C II] and CO(4-3) channel maps (Fig. 7) show an overall good correlation, with a global velocity gradient from the northeast to the southwest. The shell structure in N160D also follows the gradient, but the [C II] emission at the eastern edge of the shell extends in velocity up to 244 km s⁻¹ (it is also seen as a [C II] wing in the spectrum of position 4 in Fig. B.1). The [C II] emission to the north of N160A appears at a velocity of 235–236 km s⁻¹, and it connects to clouds at N160A, while the CO(4-3) emission is detected at 233-235 km s⁻¹ and is somewhatisolated from N160A.

Fig. 7 1 km s−1 wide channel maps of the [C II] emission (color) overlaid with contours of CO(4-3) emission in N160. The contours start at 1.5 K with an interval 1.5 K. The central velocity is written in each panel.

4.6.2 N158

The observed region is N158C. The southeast and the most northwestern CO peaks (Fig. 1) correspond to N158-2 and N158-1 respectively in Johansson et al. (1998). The channel maps (Fig. 8) indicate three different cloud components: north-south structures at right ascensions of 05:39:17 and 05:39:07, and an east-west structure with a declination of − 69:30:00. The H II region traced by Hα is also elongated in the north–south direction (slightly inclined from northeast to southwest; Fleener et al. 2010), and is located between the two north–south structures seen in [C II] emission. The eastern CO(4-3) peak corresponds to the prominent dust lane in the Hα image (Fleener et al. 2010). The strongest peak of the CO(4-3) and [C II] toward the southeast is shifted, and two compact H II regions identified by Indebetouw et al. (2004) are located between the CO(4-3) peak and the [C II] peak. The ionized gas extends also in an east-west direction at the northern edge of our observed region. This clearly shows that the [C II] emission is distributed dominantly outside the ionization front, together with the CO emission. On the other hand, the [C II] distribution is more or less continuous spatially and spectrally, while the CO(4-3) emission shows more clearly separated clouds.

Fig. 8 Same as Fig. 7 but for N158.

4.6.3 30 Dor

There are two main clouds of CO and [C II] emission in our observed regions (Fig. 1): northeast (30 Dor-10) and southwest ofR136 (30 Dor-12, Johansson et al. 1998). The brightest part of 30 Dor-10 is closer to R136 with a relatively sharp boundary on the side facing R136, while both CO(4-3) and [C II] are extended towards the northeast. Indebetouw et al. (2013) identify about one hundred clumps with the ALMA observations at the central part of 30 Dor-10. Their clumps 52 and 72 are the main contributors to our brightest peak and the neighboring peak towards the east.

In the channel maps, the peak velocities of 30 Dor-10 and 30 Dor-12 appear similar, while the emission at > 255 km s⁻¹ appears only in 30 Dor-10 (Fig. 9). In 30 Dor-10, we see four different velocity streams, connecting at the CO(4-3) and [C II] emission peaks: northeast-southwest direction at ~ 244 km s⁻¹, north-southdirection at ~ 247 km s⁻¹, east-west direction at ~ 256 km s⁻¹, and south-north direction at ~ 260 km s⁻¹. The CO(4-3) and [C II] emission show basically similar structures, except for the following few differences. At ~ 243 km s⁻¹, the [C II] emission has an additional northwest-southeast stream, which is not seen in CO(4-3). The CO(4-3) has a peak northeast of the brightest peak (see the channel map at ~ 251 km s⁻¹), which does not have a corresponding [C II] peak; [C II] instead shows a hole there (Fig. 1). At ~ 260 km s⁻¹, the CO(4-3) shows a cloud south of 30 Dor-10, while the [C II] emission extends continuously towards the south and it extends further in the velocity domain up to ~ 266 km s⁻¹.

Position-velocity diagrams along two cuts (Fig. 1) are shown in Fig. 11. Cut (a) goes through the so-called “stapler nebula”, where Kalari et al. (2018) found three molecular clouds with CO(2-1). These authors suggested that we see the tails of pillar-like structures whose ionized heads are pointing toward R136, and that the large observed velocity dispersion can be explained by the velocity difference between the head and tail of a photo evaporating cloud. In our observations, the [C II] velocity is shifted compared to the CO(4-3) velocity (Fig. 11a, see also positions 10 and 11 in Fig. B.6), which is consistent with their picture since the contribution from the head should be more dominant in [C II] emission than the CO(4-3) emission. Cut (b) goes through two weak [C II] blobs in the integrated intensity map, which appear as a velocity component at ~265 km s⁻¹. Positions 10 and 11 in Fig. B.6 also show this velocity component at 30′′ resolution. The p–v diagrams show that there is no clear connection between different velocity components among these two cuts. As discussed in Sect. 4.3, the H I spectra show an absorption feature at ~265 km s⁻¹. Unfortunately, the spatial resolution of the H I observations is insufficient to determine whether the absorption is spatially associated with the [C II] blobs.

4.7 Column densities

We derived column densities of C⁺, C, and CO at each position in our four regions. The detailed derivation is presented in Paper I. Here we describe it briefly. We assume local thermodynamic equilibrium (LTE) for each species, and a uniform excitation temperature (T_ex) for different velocity bins at each position. Assuming that CO and ¹³CO have the same T_ex, beam (area) filling factor (η), and an isotope ratio ¹² C/¹³C of 49 (Wang et al. 2009), we derived the optical depth of ¹³CO(3-2) (τ₁₃) from the intensity ratio CO(3-2)/¹³CO(3-2). With the absolute intensity of ¹³CO(3-2) and τ₁₂ we can estimate T_ex, and then the CO column density for a given beam filling factor (η). To estimate the C⁰ column density, we defined T_ex of the [C I] emission as follows. In N159 and 30 Dor, we derived T_ex by the ratio of the two [C I] intensities where both lines are detected and took its median. We used this median T_ex as a constant T_ex over the whole map in each region. For the other two regions, where we do not have the [C I] ³P₂–³P₁ line, we used the T_ex obtained in N159. To estimate the C⁺ column density, we have to assume T_ex of [C II]. Here we also used a constant T_ex over the map in each region, which is determined so that it is above the T_ex lower limit at any position of the map (i.e., for a τ ≫ 1 case).

Figure 12 shows the relation between the column density N(CO), N(C), and N(C⁺) integrated over the whole velocity range. Here we present beam averaged column densities, namely the column densities presented in the Paper I multiplied by the beam filling factor. In most of the positions N(C⁺) is in the range of 1–10 times N(CO), and N(C⁰) is 0.1–1 times N(CO). There is some hint that N159 has a lower N(C⁺)/N(CO) ratio compared to other regions, but the trend is weak. N160 appears to have two different populations; for the two clouds in the northeast, where the [C II] emission is very weak and the [C I] emission is strong, there are several positions with lower N(C⁺) and higher N(C). This may be due to a lower radiation field in this region. Toward the main cloud of N160 and in N158, the N(C)/N(CO) ratio is smaller than in N159 and 30 Dor, while the N(C⁺)/N(CO) ratio is in a similar range as N159 and 30 Dor. This result is not very sensitive to a different T_ex assumption for the [C I] emission in N160 and in N158, for example using T_ex of CO, since it is above the energy of the [C I]³P₁–³P₀ transition (24 K) at most positions. It is most likely not a radiation field effect since the representative radiation field strength is χ = 60–220 in N159 and N160 (Israel et al. 1996; Pineda et al. 2008; Lee et al. 2016) and 10³ –10⁴ in 30 Dor (Pineda et al. 2012; Chevance et al. 2016). It is also inconsistent with a scenario that a more evolved region has less CO, since 30 Dor is the most evolved region and N159 is likely to be one of the youngest regions in our analysis. A comparison with the SMC data is also puzzling. In Fig. 12 the data of N66 (plume and ridge), N25, and N88 by Requena-Torres et al. (2016) are plotted. Their N(C⁺)/N(CO) ratio falls within the range of values obtained for LMC regions, while they find a larger N(C)/N(CO) ratio. As mentioned in Sect. 1, N(C⁺)/N(CO) is expected to be large in a lower metallicity environment, while N(C)/N(CO) is expected to be less dependent on metallicity. We do not see this effect. For N25 and N88, where the radiation field is χ =20–80 (Israel & Maloney 2011), the high N(C)/N(CO) may be dueto the low radiation field. However, the data point with the highest N(C)/N(CO) comes from N66, where the radiation field is estimated as 80–570 (Israel & Maloney 2011), which is not very different to the LMC regions. The result here indicates that N(C⁺)/N(CO) and N(C)/N(CO) in resolved regions are affected significantly by local conditions and are not a simple function of metallicity.

Fig. 9 Same as Fig. 7 but for 30 Dor and contour interval of 1.0 K.

Fig. 10 [C II] integrated intensity map of N160 (color scale is the same in Fig. 1). Open blue squares show the position of stars with a spectral type of O7 or earlier, while red triangles indicate stars with spectral type O7 to B1 (Fariña et al. 2009). The blue asterisks are compact H II regions (Indebetouw et al. 2004), and the red diamond marks the position of an H2 O master (Lazendic et al. 2002).

Fig. 11 Position–velocity diagrams of [C II] (color) and CO(4-3) (contour) emission along the cuts shown in Fig. 1. The contour spacing is 0.6 K.

4.8 PDR modeling

4.8.1 Assumptions and fitting procedure

We fit the KOSMA-τ PDR model (Sect. 3) to the absolute intensities of the observed line emission and continuum SED at each position of the 30′′ resolution maps. We assume uncertainties of emission lines to be 10% for the APEX and continuum observations (Meixner et al. 2013) and 15% for GREAT and PACS observations (PACS Spectroscopy performance and calibration document; Cormier et al. 2015). In case the baseline noise is bigger than these uncertainties, we took the baseline noise as an uncertainty. We used the clump ensemble model with a lower mass limit of 10⁻³M_⊙ and an upper mass limit of 10³M_⊙. The upper mass limit is consistent with the ALMA result for the LMC (Indebetouw et al. 2013). For the dust model, Weingartner & Draine (2001) provide six different models for the LMC; two types of environment characterized by different extinction curves (LMC average and LMC 2) with three different abundances of very small grains each. LMC 2 is a supergiant shell southeast of 30 Dor that partially overlaps our observed regions (Misselt et al. 1999). We selected the LMC average environment with the abundance of very small grains of 10⁵ b_c = 2.0 (model 28; numbering by the rows of Weingartner & Draine 2001 Tables 1 and 3), which is consistent with the polycyclic aromatic hydrocarbon (PAH) abundance quoted by Galametz et al. (2013). We used a sum of the Gaussians defined by the CO velocity profiles (Sect. 4.2) as input intensities except for the [O I] emission, where we used the integrated intensity from the velocity-unresolved PACS data. In this way, the derived physical properties are the average of different gas components that are detected in CO(3-2). The effects of different model assumptions and input intensities are discussed in Sect. 4.8.3.

We selected the emission to be used in the model fit as follows. We consider two cases, where we fit only the line emission or both the line and continuum emission. We consider five continuum data (Sect. 2.6) as independent data points and treat them as if they were individual line data in the χ² fit. For the selection of line emission to be included in the fit, we have four cases: using all lines, using all lines except for the [O I] lines, using all lines except for the [C I] lines, and using only optically thin lines (¹³CO, [C I], and [O I] 145 μm). For every selection, we consider the map positions where at least four lines are detected in order to derive three model parameters (n, m, and χ, see Sect. 3 for their definitions). The exception is the case for the optically thin case fitted together with the continuum. In this case, we accept positions with three lines since otherwise most of the map positions do not fulfill the criteria, and it gives still a reasonable fit because the number of the total data points is more than four together with the continuum.

The original model grid covers log(n) = 2–7 in density, log(m) = −3 to 3 in total mass, and log(χ) = 0–6 in the UV field strength, all with a logarithmic step of 1. We first extended the mass grid to log (m) = 6. In a clumpy model, no extrapolation is needed to extend the mass range because it is just a scaling of the total intensity. We interpolated the model to a 0.1 step in the logarithmic scale of n, m, and χ. Then we calculated the reduced χ-square (χ²) of all the grid points and take the grid point with the minimum χ² as a solution. In this way, we can avoid to fall into a local minimum.

Fig. 12 Column density N(C) (upper panels) or N(C+) (lower panels) against N(CO) for the beamfilling factor of 0.1 (left panels) and 0.5 (right panels). The SMC data are taken from Requena-Torres et al. (2016). The black line indicates that the two column densities arethe same.

Fig. 13 Results of fitting the line and continuum emission at the CO peaks in N159 W (upper two panels), in N160A (middle panels), and in 30 Dor-10 (lower panels). Black points are the observed data and colored points represent the models. The left panels are the line SED of CO(2-1) to (6-5), [C I]3 P1–3P0, [C I]3P2–3P1, [C II], [O I] 63 μm, and [O I] 145 μm. 12 CO and 13CO are plotted together. For N159, CO(2-1) and 13CO(2-1) are not used in the fit because it is a pointed observation and the beam sizes are slightly different. The red model uses all line emissions but no continuum in the fit. The blue model excludes [O I] from the fit but includes the continuum emission. The green model fits all the lines and the continuum emission. The purple model fits only optically thin lines and the continuum. In 30 Dor the purple and green models give the same fit results (Table 5). Input line intensities are the sum of Gaussians except for the [O I] lines, where the integrated intensities are used (see Table 4).

4.8.2 Result

We calculated the best-fitting model at each position of the 30′′ resolution maps of the four star-forming regions. Figure 13 and Table 5 show an example fit (the coordinates and input intensities are listed in Table 4). In most of the observed positions within the four star-forming regions (except for the CO peak at N159 W), we find the following overall trend. When fitting with only the line emission (red lines in Fig. 13), the lines are relatively well fit, with some cases the model underestimating the [C I] emission. However the continuum emission is clearly underestimated especially at shorter wavelengths, indicating that the dust temperaturein the best-fitting model is too low. It is consistent with the fact that the CO-ladder starts to be flat already around J = 5 or 6, resulting in an underestimate of CO(6-5). When we exclude the [O I] lines and include the continuum in the fit (blue lines), the estimated UV field is much stronger. The continuum SED as well as the rising trend from CO(4-3) to CO(6-5) are relatively well reproduced, but the [O I] lines are heavily overestimated and the [C I] emission is even more severely underestimated than in the first model. When we fit everything (green lines), the fit of [O I] improves with compromising other parts. In most cases, fitting only optically thin lines (purple lines) or fitting without [C I] lines (not shown) does not change the result much, or makes the fit worse.

The fundamental problem for fitting the line emission is that the CO and [C I] have two solutions at low n – low χ and at high n – high χ range for their given intensities (Fig. 14). The [C II] intensity contours are also aligned from low n – low χ to high n – high χ. In principle, the [O I] emission lines have a different dependency on n and χ, and hence including the [O I] emission in the fit should provide a good constraint. However, for most positions that we study here, the solutions of n and χ indicated by the [O I] emission do not overlap with those that other emissions suggest (Fig. 14a). This is why it looks like that the fit excluding[O I] (blue lines in Fig. 13) look the most reasonable for lines other than the [O I]. However, in this case the low n – low χ solution and high n – high χ solution are degenerate (Fig. 14b), and very different physical conditions are obtained for spatially adjacent pixels. It is also not justified to simply ignore the [O I] emission, because it is difficult to explain the overestimate of both [O I] lines. For these reasons, we consider that fitting with all the lines and the continuum (green model) gives the most stable solution for representative physical properties in the sense that it excludes physical properties that would yield completely different line or continuum intensities compared to what is observed.

Since we fit the absolute intensities, the model results include the area filling factor as well (Table 5). However, it is sensitive to the derived density and it ranges from the order of 0.01 to 100 in four regions even with the same model assumption, which is again an effect of the degeneracy of the low n – low χ solution and high n – high χ solution.

The CO peak of N159 W shows a different result (Fig. 13, top panel). There all models that fail to reproduce the [C II] emission and the [O I] 63 to 145 μm ratio (Lee et al. 2016). When we fit only the optically thin lines, the model overestimates all the CO, [C II], and [O I] 63 μm line intensities, which is qualitatively consistent with an argument that it is an optical depth effect among clumps. Quantitatively, we assume that the N_cloud clouds are aligned along the line of sight, with each having a line intensity of I_line and the optical depth of τ_line. The model predicts the line intensity of I_line × N_cloud, while the observed line intensity is $I_{\text{line}}\times {\sum }_{i=0}^{N_{\text{cloud}}-1}{e}^{-i{\tau }_{\text{line}}}$ (Okada et al. 2003). To obtain one order of magnitude difference between the model and observations as seen in Fig. 13, we need N_cloud = 10 for a τ_line = ∞ case and N_cloud = 15 for a τ_line = 1 case. The fitted clump ensemble model when using only optically thin lines gives an area filling factor of > 20 for any combination of model assumptions, which supports the idea of heavily overlapping clumps, although the filling factor is sensitive to the derived density as mentioned above.

Lee et al. (2016) modeled [O I] 145 μm, [C I]³P₂–³P₁, [C II] and FIR with the Meudon PDR code and derived P/k = 10⁶ K cm⁻³ and G₀ = 100. Their model underestimates the CO emission, especially high-J CO, and they concluded that CO is heated by low-velocity shocks. When we use a non-clumpy model of the KOSMA-τ model and exclude the CO emissions from the fit, we reproduce qualitatively their results. In Fig. 15, we compare our clumpy PDR model with the high-J CO intensities observed with the Spectral and Photometric Imaging REceiver (SPIRE) on Herschel provided by Lee et al. (2016). The gap between CO(6-5) and CO(7-6) is partially because the SPIRE data are convolved to the 42′′ beam size. The original beam size of CO(7-6) is 33′′ (Lee et al. 2016), close to our 30′′, and its value is 15% higher. The figure shows that the clumpy model provides a reasonable fit to the CO up to around J = 10. The very shallow slope of the higher-J CO-ladder is consistent with the shock contribution, while low- and mid-J CO can be also explained by a PDR model when we introduce clumpiness. In fact, the very similar line profiles of the CO and [C I] is consistent with the picture that they are both emitted by a PDR gas, although the proposed shock in Lee et al. (2016) has a velocity of only ~10 km s⁻¹, which may be hard to see in an emission line with a width of ~8 km s⁻¹ (Paper I).

For the continuum SED, the model prediction of the SED is typically broader than the observed SED, independent of which model we adopt. It is typically characterized by a strong drop of 70 μm continuum (e.g., N159 in Fig. 13). When we choose a dust model with fewer small grains (see Sect. 4.8.3), the model SED is slightly narrower because of weaker 70 μm emission, but the fit to the observed SED is not significantly better. This is consistent with the model SED by Bron (2014) where the contribution of PAHs to the 70 μm emission is a few percent for χ = 1–1000. At the column density peak, it is possible that the optically thin assumption for the 70 μm continuum emission does not fully hold; we estimated the total hydrogen column density from N(CO), N(C⁰), and N(C⁺) in Sect. 4.7 with the carbon abundance of 7.9 × 10⁻⁵ (Garnett 1999) and converted it to a dust opacity using the extinction cross section (Weingartner & Draine 2001), which is at maximum 0.15 at the peak of 30 Dor-10 in case of η = 0.1. This adds ~10% uncertainty to a few positions around the CO peak, but does not have a significant impact on the continuum fit, although it is not excluded that the beam filling factor is less than 0.1 and the 70 μm continuum is marginally optically thick. Chevance et al. (2016) suggest that the FIR continuum has a significant contribution from the ionized gas outside of the molecular clouds in 30 Dor by analyzing its correlation with PAH and [O III] emission. This may cause an overestimate of the UV field in the fit with the continuum SED, but it does not explain the observed narrow peak in the continuum SED and its contribution should be higher outside of molecular clouds.

Fig. 14 Contours of χ2 at the 30 Dor-10 CO peak position for the green model (fit all line and continuum emission; panel a) and the blue model (excluding [O I]; panel b) in the lower panel of Fig. 13.

Fig. 15 CO ladder including the high-J CO data by Lee et al. (2016) at the CO peak in N159 W. Black data points and the model lines are the same as in the upper panel of Fig. 13. Blue data points are SPIRE observations at 42′′ resolution (Lee et al. 2016).

4.8.3 Comparison of different inputs and model assumptions

Here we discuss the effect of different model assumptions and inputs on the obtained physical properties. First, we compare the results between the model with an upper mass limit of 10³ and 10¹M_⊙. In the model, a smaller clump has a higher mean density (see the mass-size relation in Sect. 3), and removing clumps with the mass range of 10¹–10³ increases the mean ensemble density by a factor of 3.7. The derived densities using an upper mass limit of 10¹ in the fit are typically higher than the results with an upper mass limit of 10³ by a factor of 1–4. On the other hand, the choice of the lower mass limit between 10⁻³ and 10⁻² does not make a significant difference. This is because clumps with a mass of 10⁻³ are dominated by C⁺ and do not contribute much to the CO and [C I] emissions in most of our density and UV field strength parameter space; for example, with n =10⁴ and χ =10², a clump with m = 10⁻³ gives line ratios of [C II]/CO(3-2) ~ 5 × 10⁵ and [C II]/[C I]³P₁–³P₀ ~ 7 × 10⁴.

We also ran the KOSMA-τ model with the smallest number of very small grains with the extinction curve of LMC average (model 26) and LMC 2 (model 29) in Weingartner & Draine (2001). The extinction curve in the LMC 2 shows a weaker 2175 Å bump than that of the LMC average. Thus, the model 29 has a size distribution with the fewest small grains in Weingartner & Draine (2001). The smaller amount of small grains results in a lower heating efficiency, and consequently a lower temperature, which changes the continuum SED shape. The most affected emission line is [O I] because the [O I] lines have the highest level energies. When we compare the fitting result of model 26 and model 28 (same extinction curve, but a different amount of very small grains), the derived density and UV field of model 26 are higher than those of model 28, because it has a lower heating efficiency and requires a larger density and/or stronger UV field to produce the same amount of the line emission. The derived mass is not significantly different, because the mass is well constrained by the continuum emission, which does not depend significantly on small grains. The density and UV field derived with the model 29 is even slightly higher, but the difference between model 26 and model 29 (different extinction curve) is not as significant as the difference between model 28 and 26.

As input intensities, we examine the following cases; (1) integrated intensities over the whole velocity range, (2) a sum of the Gaussians defined by the CO velocity profiles together with the integrated intensity of [O I] emissions, (3) same as (2) but scale the PACS [O I] intensities using the fraction of the sum of Gaussians to the integrated intensity for the [C II]. The default input used above is case (2). Case (2) corresponds to the assumption that the [O I] emission has the same line profile as CO lines, and case (3) corresponds to the assumption that the [O I] emission has the same line profile as [C II]. The pointed observations at a few positions in N159 and 30Dor (Sect. 4.1) show that the real case is somewhere between those two assumptions. The obtained mass is stable against the choice of the input intensities. The obtained density and UV field shows a significant scatter (≤ 1 dex) but no systematic trend.

4.8.4 Spatial distribution of density, mass, and UV field

Figures 16–19 show the spatial distribution of density, mass, and UV field in our four regions obtained by the model fit using all line and continuum emissions. Areas where at least one of the [O I] lines is detected are indicated in the figures, since the fitting with and without [O I] emission makes a difference to the obtained range of density and the UV field strength (see Sect. 4.8.2). On the other hand, the derived mass distribution is stable against different fitting and model assumptions, tracing the CO clouds.

When we only look at areas where [O I] emission detected, the UV field and mass distribution appear anti-correlated, indicating that the areas outside of giant molecular clouds are more excited by OB stars that are not embedded by molecular clouds. Surprisingly, the distribution of the density is quite similar to that of the UV field. This may mean that only dense clumps survive in higher UV field regions. However, it contradicts the interferometry observations, where high-density tracers (HCO⁺ and HCN) are detected toward the center of molecular clouds in N159 (Seale et al. 2012) and the detection of NH₃ toward N159 W (Ott et al. 2010). In N160, there are patches with almost two orders of magnitude stronger UV field than surrounding pixels, which are not fully consistent with the positions of early type stars (Fig. 10). Also, the derived density at these positions are 10⁶ –10⁷ cm⁻³, which suggests a very high pressure. It is unlikely that this jump in the UV field is real. Instead, we think that it is due to the degeneracy problem between a high n–high χ solution and a low n–low χ solution (see Sect. 4.8.2). In 30 Dor, the derived UV field strength close to R136 is somewhat too high because it isslightly higher than the radiation field derived by the star radiation in the plane of R136 (Chevance et al. 2016).

4.8.5 Discussion

Overall, the results of the PDR model fit suggest that a model with one clump ensemble component cannot reproduce all of the line and continuum emissions consistently at most of the observed positions. A natural next step would be to consider two clump ensembles with different physical properties. With the velocity-resolved spectra, it is well justified; instead of using the sum of Gaussians that conform to the CO line profiles, we would assign each velocity component to a separate clump ensemble. They should be fitted simultaneously, with free parameters to distribute the continuum emission to each component. We did not pursue this strategy in this paper since we do not have enough data points to fit additional free parameters. A very useful future addition to the present data would be velocity-resolved high-J CO spectra and mid- (to high-J) ¹³CO spectra, to better constrain the density and UV field, and investigate the origin of low-J and high-J CO emissions.

The results of the PDR model fit using different tracers provide a warning in interpreting the “best” fit of the model. When excluding one line or continuum emission from the fit, the physical properties of the best fit model vary significantly, and can jump by a few orders of magnitude in some cases. This clearly indicates that we should use as many tracers as possible to obtain a consistent picture.

5 Summary

We mapped 3–13 arcmin² areas in 30 Dor, N158, N160, and N159 in [C II] 158 μm with GREAT on board SOFIA, as well as CO(2-1) to (6-5), ¹³CO(2-1) and (3-2), [C I] ³P₁–³P₀ and [C I] ³P₂–³P₁ with APEX. We also observed the velocity-resolved [O I] 145 and 63 μm at selected positions in N159 and 30 Dor for the first time. The results are as follows.

• In all four star-forming regions, the line profiles of the CO, ¹³CO, and [C I] emission are similar while [C II] typically has a wider line profile or an additional velocity component. On average, 30% of the emission in [C II] cannot be reproduced by the CO-defined line profile. This fraction is lower toward molecular clouds.

• The [O I] line profiles match those of CO at some positions, but are more similar to the [C II] profiles at other positions. This indicates that we cannot simply assume that the velocity components of the [O I] emission are the same as either CO or [C II].

• Using the H I and [C II] line profiles, we estimated that contribution of atomic gas to the [C II] emission is 15% orless. The thermal pressure that is required to emit the [C II] associated with the atomic gas is 4 × 10³–10⁵ K cm⁻³, which ishigher than the standard Galactic ISM value and consistent with previous studies. For some positions in 30 Dor, the H I absorption coincides with a velocity component in [C II], which may be because of a cold molecular cloud with a mixture of atomic hydrogen.

• We interpret the different line profiles among CO, [C II], and [O I] in the LMC as contributions from spatially separated clouds and/or clouds in different phases, which give different line ratios depending on their physical properties.

• We investigate channel maps in individual regions and distinguish clouds.

• We derived the column density of CO, C⁰, and C⁺ and compared them to those of SMC regions measured by Requena-Torres et al. (2016). We do not see a clear correlation between N(C⁺)/N(CO) and metallicity. The trend cannot be fully explained by a different UV field strength and metallicity.

• We modeled the line and continuum emission using the latest KOSMA-τ PDR model, which treats the dust-related physics consistently and computes the dust continuum SED as well as the line emission. At most positions, a single clump ensemble does not satisfactorily reproduce all the observed emissionlines and continuum. Toward the CO peak at N159 W, we propose a scenario in which the CO, [C II] and [O I] 63 μm emission are affected by mutual shielding between clumps.

• We show that the best fit model results depend sensitively on which combination of continuum and emission lines are used in the fit.

Fig. 16 Derived density, mass, UV field, and the χ2 of the fit (from left to right panels) for the N159 region. The results are obtained using a fit to all emission lines and the continuum (green lines in Fig. 13). Areas surrounded by black lines are where at least one of the [O I] lines isdetected.

Fig. 17 Same as Fig. 16 but for N160.

Fig. 18 Same as Fig. 16 but for N158.

Fig. 19 Same as Fig. 16 but for 30 Dor.

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 A3 and C1, funded by the Deutsche Forschungsgemeinschaft (DFG). We thank M.-Y. Lee for fruitful discussions and providing the data of Herschel/SPIRE at N159 W. We thank R. Braun for providing the opacity-corrected H I column density in the LMC.

Appendix A: Gain correction using the total raw count

In the HFA data, we often see that the total raw count is drifting (Fig. A.1). Since the raw count of the load measurements follow the same drifting trend in general, this is most likely determined by a gain drift of the receiver rather than a change of the sky transmission. The lower panel of Fig. A.1 shows the sky−hot count divided by the gain (see Guan et al. 2012, for the detailed formula). When the gain is calculated from the last load measurement before the sky observations, the sky−hot shows a strong drift between two load measurements (black points), which leads a corresponding variation of the fitted atmospheric transmission. The first method to correct it is to scale the raw count of individual observations to the one right after each load measurement, assuming that the total count variation between two load measurements is caused by the gain drift. The resultant sky−hot is presented as green points in the lower panel of Fig. A.1. The drift between two load measurements is corrected, but a new load measurement often does not provide the same level of sky−hot; it is jumping around. For our paper, relative intensities between different measurements are more important than the absolute intensity, because firstly we would like to achieve a good S/N by averaging measurements over 10–20 min of observing time and secondly we would like to compare four raster points in 30 Dor which are spread among the whole leg. Therefore, we took a somewhat radical approach: we scaled all hot measurement to the first hot measurement, all cold measurement to the first cold measurement, and all sky measurement to the first sky measurement for each pixel of each source. The sky−hot after this gain correction is shown as red triangles in lower panel of Fig. A.1. This assumes that the sky opacity is constant over the observing time of one source, and any change in the total raw count is caused by a gain variation. This correction implies that all measurements have the same atmospheric opacity correction, but the absolutecorrection is uncertain since the choice of the measurement to be scaled to is arbitrary. We note that the lower panel of Fig. A.1 indicates that taking the first measurement is not a bad choice because it has an intermediate gain value. This method is justified by following reasons. The variation in sky−hot among different load measurement blocks (green points) is not synchronized among different pixels. If this variation is because of the variation of the sky opacity, all pixels should show the same trend. The second reason is that for similar observations of another source (NGC2024; Graf et al., in prep.), strong baseline structure is suppressed significantly by our adopted scaling method.

Fig. A.1 Upper panel: total counts averaged over 5000 channels with a good sensitivity and outside of the source emission versus time in the 30 Dor leg for the pixel HFAV_PX00. The red and blue points indicate hot-load and cold-load measurements, respectively, and black points are the sky observations. Lower panel: sky count minus hot count divided by the gain, averaged over the same channels. See Guan et al. (2012) for the detailed formula. Black diamonds show the case without gain correction, green asterisks are after the gain correction within observations that share the same load measurement, red triangles represent the data after the final grain correction (see text).

In terms of the absolute intensity scale, we estimate the uncertainty by the variation of sky−hot among different load measurement blocks (green points), which is translated to the variation of the atmospherictransmission. Among the data used in this paper, the uncertainty of the absolute intensity is ≤ 66% for HFAV_PX00 and ≤33% for other HFAV pixels. We did the same analysis for LFA channels, although they were much more stable. The maximum uncertainty of the absolute intensity for LFA pixels is 20%.

Appendix B: Integrated intensity maps and selected spectra in N160, N158, and 30 Dor

Here we show the integrated intensity maps and spectra at selected positions of N160, N158, and 30 Dor as shown in Paper I for N159.

Fig. B.1 Integrated intensity maps (colors, 30′′ resolution) overlaid with contours of the IRAC 8 μm emission in N160. The red lines outline the observed area. Blue asterisks in CO(3-2) and [C II] maps mark the positions where the spectra shown in Fig. B.2 have been extracted.

Fig. B.1 continued.

Fig. B.2 Spectra at selected positions in N160, marked in the CO(3-2) and [C II] panels of Fig. B.1.

Fig. B.3 Integrated intensity maps (colors, 30′′ resolution) overlaid with contours of the IRAC 8 μm emission in N158. The[C I]3P1–3P0 is not shown because it is detected only at a few positions around the two CO peaks. Blue asterisks in CO(3-2) and [C II] maps mark the positions mark the positions where the spectra shown in Fig. B.4 have been extracted.

Fig. B.4 Spectra at selected positions in N158, marked in the CO(3-2) and [C II] panels of Fig. B.3.

Fig. B.5 Integrated intensity maps (colors, 30′′ resolution) overlaid with contours of the IRAC 8 μm emission in 30 Dor. The red lines outline the observed area. Blue asterisks in CO(3-2) and [C II] maps mark the positions where the spectra shown in Fig. B.6 have been extracted.

Fig. B.5 continued

Fig. B.6 Spectra at selected positions in 30 Dor, marked in the CO(3-2) and [C II] panels of Fig. B.5.

References

All Tables

All Figures

Fig. 1 Center panel: IRAC 8 μm map with boxes showing the areas we observed. Four inset color maps with contours: integrated intensity map of [C II] (color) overlaid with contours of CO(4-3) integrated intensity at 16′′ resolution. The red lines outline the area observed in [C II]. Contour levels are at intervals of 5 K km s−1 for N158, 10 K km s−1 for N160 and 30 Dor, and 15 K km s−1 for N159. The unit of color bar is K km s−1. In the maps of N159 and 30 Dor, the positions of the [O I] observations are marked following the numbering in Table 2. In 30 Dor, the position of R136 is shown as a red circle, and the two cuts to make position-velocity diagrams (Fig. 11) are shown as black lines. In N160, two H II regions N160A1 and A2 (Martín-Hernández et al. 2005) are marked.

In the text

	Fig. 2 [O I] 145 and 63 μm spectra in N159 at the positions marked in Fig. 1. The CO(4-3) and [C II] spectra were extracted from the maps with a resolution of 16′′, while the [O I] observations are shown with the original spatial resolution (13′′ for the [O I] 145μm and 6.3′′ for the [O I] 63 μm).
In the text

	Fig. 3 Same as Fig. 2 but for 30 Dor. Here the [O I] 63 μm spectra are the average of the four raster positions (see text).
In the text

	Fig. 4 [O I] 63 μm spectra of four raster positions and their average (black histogram) in 30 Dor positions 1, 4, 5, and 6.
In the text

	Fig. 5 [C II] and H I spectra at a [C II] blob in N159 (05:39:45.9, −69:44:32.2, upper panel)and 30 Dor (05:38:31.6, −69:04:50.0, lower panel) at 1′ spatial resolution.
In the text

	Fig. 6 Integrated intensity maps of the red-side (left panel) and blue-side (right panel) residual of the [C II] emission after subtracting the CO-defined line profile (see text) in N159 (color, at 30′′ resolution) overlaid with contours of Hα (Chen et al. 2010).
In the text

	Fig. 7 1 km s−1 wide channel maps of the [C II] emission (color) overlaid with contours of CO(4-3) emission in N160. The contours start at 1.5 K with an interval 1.5 K. The central velocity is written in each panel.
In the text

	Fig. 8 Same as Fig. 7 but for N158.
In the text

	Fig. 9 Same as Fig. 7 but for 30 Dor and contour interval of 1.0 K.
In the text

Fig. 10 [C II] integrated intensity map of N160 (color scale is the same in Fig. 1). Open blue squares show the position of stars with a spectral type of O7 or earlier, while red triangles indicate stars with spectral type O7 to B1 (Fariña et al. 2009). The blue asterisks are compact H II regions (Indebetouw et al. 2004), and the red diamond marks the position of an H2 O master (Lazendic et al. 2002).

In the text

	Fig. 11 Position–velocity diagrams of [C II] (color) and CO(4-3) (contour) emission along the cuts shown in Fig. 1. The contour spacing is 0.6 K.
In the text

	Fig. 12 Column density N(C) (upper panels) or N(C+) (lower panels) against N(CO) for the beamfilling factor of 0.1 (left panels) and 0.5 (right panels). The SMC data are taken from Requena-Torres et al. (2016). The black line indicates that the two column densities arethe same.
In the text

Fig. 13 Results of fitting the line and continuum emission at the CO peaks in N159 W (upper two panels), in N160A (middle panels), and in 30 Dor-10 (lower panels). Black points are the observed data and colored points represent the models. The left panels are the line SED of CO(2-1) to (6-5), [C I]3 P1–3P0, [C I]3P2–3P1, [C II], [O I] 63 μm, and [O I] 145 μm. 12 CO and 13CO are plotted together. For N159, CO(2-1) and 13CO(2-1) are not used in the fit because it is a pointed observation and the beam sizes are slightly different. The red model uses all line emissions but no continuum in the fit. The blue model excludes [O I] from the fit but includes the continuum emission. The green model fits all the lines and the continuum emission. The purple model fits only optically thin lines and the continuum. In 30 Dor the purple and green models give the same fit results (Table 5). Input line intensities are the sum of Gaussians except for the [O I] lines, where the integrated intensities are used (see Table 4).

In the text

	Fig. 14 Contours of χ2 at the 30 Dor-10 CO peak position for the green model (fit all line and continuum emission; panel a) and the blue model (excluding [O I]; panel b) in the lower panel of Fig. 13.
In the text

	Fig. 15 CO ladder including the high-J CO data by Lee et al. (2016) at the CO peak in N159 W. Black data points and the model lines are the same as in the upper panel of Fig. 13. Blue data points are SPIRE observations at 42′′ resolution (Lee et al. 2016).
In the text

	Fig. 16 Derived density, mass, UV field, and the χ2 of the fit (from left to right panels) for the N159 region. The results are obtained using a fit to all emission lines and the continuum (green lines in Fig. 13). Areas surrounded by black lines are where at least one of the [O I] lines isdetected.
In the text

	Fig. 17 Same as Fig. 16 but for N160.
In the text

	Fig. 18 Same as Fig. 16 but for N158.
In the text

	Fig. 19 Same as Fig. 16 but for 30 Dor.
In the text

Fig. A.1 Upper panel: total counts averaged over 5000 channels with a good sensitivity and outside of the source emission versus time in the 30 Dor leg for the pixel HFAV_PX00. The red and blue points indicate hot-load and cold-load measurements, respectively, and black points are the sky observations. Lower panel: sky count minus hot count divided by the gain, averaged over the same channels. See Guan et al. (2012) for the detailed formula. Black diamonds show the case without gain correction, green asterisks are after the gain correction within observations that share the same load measurement, red triangles represent the data after the final grain correction (see text).

In the text

	Fig. B.1 Integrated intensity maps (colors, 30′′ resolution) overlaid with contours of the IRAC 8 μm emission in N160. The red lines outline the observed area. Blue asterisks in CO(3-2) and [C II] maps mark the positions where the spectra shown in Fig. B.2 have been extracted.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.2 Spectra at selected positions in N160, marked in the CO(3-2) and [C II] panels of Fig. B.1.
In the text

	Fig. B.3 Integrated intensity maps (colors, 30′′ resolution) overlaid with contours of the IRAC 8 μm emission in N158. The[C I]3P1–3P0 is not shown because it is detected only at a few positions around the two CO peaks. Blue asterisks in CO(3-2) and [C II] maps mark the positions mark the positions where the spectra shown in Fig. B.4 have been extracted.
In the text

	Fig. B.4 Spectra at selected positions in N158, marked in the CO(3-2) and [C II] panels of Fig. B.3.
In the text

	Fig. B.5 Integrated intensity maps (colors, 30′′ resolution) overlaid with contours of the IRAC 8 μm emission in 30 Dor. The red lines outline the observed area. Blue asterisks in CO(3-2) and [C II] maps mark the positions where the spectra shown in Fig. B.6 have been extracted.
In the text

	Fig. B.5 continued
In the text

	Fig. B.6 Spectra at selected positions in 30 Dor, marked in the CO(3-2) and [C II] panels of Fig. B.5.
In the text