© V. Lebouteiller et al. 2019

1 Introduction

The [C II] 158μm line emits undera variety of conditions in the interstellar medium (ISM) corresponding to the cold and warm neutral medium and warm ionized gas, owing to the relatively low ionization potential to produce C⁺ ions (11.3 eV) and to the relatively low energy of the ²P_3∕2 fine-structure level (91.3 K). In the neutral gas, C⁺ may exist in the H⁰ phase but also in the H₂ phase, in regions where CO is photodissociated while H₂ is self-shielded and shielded by dust, the so-called CO-dark molecular gas (e.g., Madden et al. 1997; Grenier et al. 2005; Wolfire et al. 2010). The [C II] line has been used to trace the star-formation rate (SFR; e.g., De Looze et al. 2014; Pineda et al. 2014), to infer physical conditions in photodissociation regions (PDRs), and to calculate the CO-to-H₂ conversion factor (X_CO; e.g., Jameson et al. 2018; Herrera-Camus et al. 2017; Pineda et al. 2017). The growing number of observations in both Galactic and extragalactic environments (with routine detections at z > 6; e.g., Aravena et al. 2016) has renewed the interest in understanding [C II] as a diagnostic tool.

It is crucial to provide astrophysical experiments that can isolate the phase contributions to the [C II] emission toward well-chosen regions and to examine how the origin of [C II] depends on environmental parameters such as the metallicity or the radiative and mechanical feedback from young stellar clusters. On the one hand, metallicity plays a role through the lower abundance of dust, resulting in less efficient shielding from UV photons and to a larger layer of CO-dark H₂ gas in PDRs (Madden & Cormier 2018, Madden et al., in prep.). On the other hand, young massive stars shape the surrounding ISM, thereby modifying the relative filling factors of warm ionized gas and molecular clouds, resulting in a lower “effective” extinction on average when measured over the scale of a star-forming region. The age of the molecular cloud in which stars will form is another important parameter, notably because of the H₂ formation timescale (e.g., Franeck et al. 2018).

It is often assumed that most of the [C II] emission arises from a given dominant ISM phase or alternatively that the relative contributions from the various ISM phases can be recovered from photoionization and photodissociation models (e.g., Cormier et al. 2015). Another method to disentangle the origin of [C II] is to compare its velocity profile to those of CO and H I 21 cm, which are assumed to trace the “CO-bright” H₂ gas and the H⁰ gas, respectively, with the remaining [C II] emission attributed to CO-dark H₂ gas (the contribution from the ionized gas being usually determined indirectly for lack of reliable velocity-resolved ionized gas tracers). Using velocity profiles, Pineda et al. (2014) studied [C II] in the Milky Way with the Herschel GOT C+ survey and found approximately equal contributions (20− 30%) from dense PDRs, CO-dark H₂ gas, cold atomic gas, and ionized gas. The warm neutral medium does not seem to contribute significantly to the observed [C II] intensity (see also Fahrion et al. 2017). Pineda et al. (2014) also find that the extragalactic SFR relationship with gas surface density can be recovered only by combining [C II] from all the ISM phases, suggesting that SFR determinations using PDR-specific tracers may be difficult to calibrate if the fraction of UV photons absorbed in PDRs varies significantly across objects or within star-forming regions. Langer et al. (2014) calculated that the fraction of CO-dark H₂ gas, f_dark, in the Herschel GOT C+ survey is ≈75% in the diffuse clouds and down to ≈20% in dense clouds. Toward the Galactic star-forming region M 17-SW, Pérez-Beaupuits et al. (2015) found that about half of the [C II] traces the CO-dark H₂ gas while ≈ 36% originates from the ionized gas.

Metallicity effects and massive stellar feedback effects are conveniently probed through observations of nearby star-forming dwarf galaxies. Fahrion et al. (2017) examined the ISM at ≈ 200 pc scales in the nearby star-forming galaxy NGC 4214 and found that about half of the [C II] traces the atomic gas and that most (≈ 80%) of the H₂ gas is CO-dark. The authors also found that f_dark is higher in the low-metallicity diffuse region where a super stellar cluster is located, although the relative influence of metallicity and feedback remains uncertain. Magellanic Clouds allow smaller spatial scales to be attained by resolving star-forming regions. Okada et al. (2015) observed N 159 in the Large Magellanic Cloud (LMC; ≈ 1∕2Z_⊙) and found that the fraction of [C II] that cannot be associated with CO increases from ≈ 20% around the CO clumps to ≈50% in the interclump medium, and that the overall fraction of [C II] associated with the ionized gas is ≤ 15% (see also Okada et al. 2019 for other regions). In LMC-30 Dor, 90% of [C II] originate in PDRs and 80−99% of H₂ is not traced by CO (Chevance et al. 2016; Chevance 2016). Requena-Torres et al. (2016) observed several star-forming regions in the Small Magellanic Cloud (SMC; ≈ 1∕5Z_⊙) and found that most of the H₂ gas is traced by [C II] and not CO. Using 54 lines of sight in the LMC and SMC, Pineda et al. (2017) find that most of the molecular gas is CO-dark. Overall, the finding that CO-dark H₂ gas is predominant in low-metallicity environments is consistent with the picture of CO-emitting regions occupying a smaller filling factor due to the relatively lower dust-to-gas mass ratio, and it is in fact possible to obtain X_CO conversion factors close to the Milky Way value in the Magellanic Clouds if filling factor effects are accounted for, as shown by Pineda et al. (2017). However, the dependence of f_dark gas with metallicity may be indirect, and model results from the Herschel Dwarf Galaxy Survey (DGS; Madden et al. 2013) indicate that f_dark is mostly a function of the effective cloud extinction, with the latter showing some dependence with metallicity (Madden et al., in prep.).

Summarizing these results, [C II] traces a significant fraction of the H₂ gas and the fraction of CO-dark H₂ gas increases from dense CO peaks to the diffuse medium (see also Mookerjea et al. 2016 in M 33), and also increases with lower metallicity, and stronger stellar feedback (radiation and/or dynamical and mechanical). The fraction of [C II] in the ionized gas is usually ≲20% as also found in the global analysis of the DGS (Cormier et al. 2015) and in resolved regions such as LMC-30 Dor (Chevance et al. 2016) and IC 10 (Polles et al. 2019). From a kinematics point of view, the [C II] line is always found to be broader than CO and [C I] but narrower than H I, with a profile agreeing better with CO (see also Braine et al. 2012; de Blok et al. 2016).

In a first publication (Lebouteiller et al. 2012), we investigated N 11B, part of the second largest (≈ 150 pc in diameter) giant H II region N 11 in the LMC after 30 Dor. We found a remarkable correlation between the total cooling rate traced by [C II]+[O I] and the polycyclic aromatic hydrocarbon (PAH) mid-infrared emission, suggesting that [C II] emission is predominantly originating from PDRs with a uniform photoelectric-effect heating efficiency. In the present study, we examine the velocity structure of [C II] in N 11B and other regions within N 11 obtained with the GREAT instrument (Heyminck et al. 2012) onboard the SOFIA telescope (Young et al. 2012). LMC-N 11 has been studied in detail especially at infrared and submillimeter wavelengths (e.g., Israel & Maloney 2011; Herrera et al. 2013; Galametz et al. 2016) and has been fully mapped in CO(1–0) (Israel et al. 2003; Wong et al. 2011) and H I 21 cm (Kim et al. 2003) (Fig. 1). The N 11 region was chosen in order to access various environments (e.g., PDRs, quiescent CO clouds, H II regions, ultracompact H II regions; see e.g., Lebouteiller et al. 2012; Galametz et al. 2016) which were expected to result in distinctive [C II] velocity profiles.

The objectives are (i) to measure the quantity of CO-dark H₂ gas traced by [C II] and the fraction of molecular gas that is CO-dark, (ii) to identify potential [C II] components associated with atomic gas, and (iii) to probe the influence of the environment (in particular stellar feedback). A specific focus is given in the present study on the velocity profile decomposition method. The comparison between the [C II], CO, and H I spectral profiles is complex and methods usually determine average properties of the spectral profile along various lines of sight or else use profile fitting starting with the tracers with relatively simple velocity structure and adding velocity components as needed for the more complex ones. Here we use a statistical approach with as few assumptions as possible on the number of components and their properties. The observations are described in Sect. 2. We derive velocity-integrated properties in Sect. 3. The profile decomposition and the associated results concerning the physical properties of the components are presented in Sect. 4.

Fig. 1 LMC-N 11 in CO, Hα, H I, [C II], [O I], and [O III]. The SOFIA/GREAT pointings are shown as black/white circles and labelled in the CO panel. The black rectangles show the Herschel/PACS coverage. Region names are labelled in the [C II] panel. GREAT pointings #1 (N 11A) and #9 (part of N 11D) were not observed in [O III], and pointing #12 (N 11 “South”) was not observed with PACS.

2 Observations

The observations used in this study are summarized in Table 1. Below we describe the details of each set of observations.

2.1 SOFIA/GREAT

Twelve pointings were observed within LMC-N 11 (Table 2; Fig. 1) with SOFIA/GREAT in [C II] 158μm and [N II] 205μm as part of program 01_0030 (PI: Lebouteiller). Observations were conducted on 2013, July 19 and 28 deploying from Christchurch, New Zealand.

Most pointings were previously identified as CO-bright peaks in the Magellanic MOPRA Assessment (MAGMA) survey (Wong et al. 2011) or as [C II] peaks in the Herschel/PACS maps (Fig. 1). A few additional pointings were included, notably toward the young stellar cluster LH 10 in N 11B where the gas is mostly ionized, resulting in bright 24 μm and [O III] 88 μm (Lebouteiller et al. 2012). Overall, the pointings span various environments including PDRs (e.g., #2), a quiescent CO cloud (#12), an ultracompact H II region (#3), stellar clusters (#5), among others (Table 2).

The [C II] line was observed in the GREAT L#2 channel while [N II] was observed in the L#1 channel. We used X(A)FFT spectrometer back-ends. The [C II] line was detected toward all pointings. The [N II] line was observed only toward pointings #1 to #6 but was not detected. Some scans were affected by emission in the offset field¹ for the observation of N 11 I (#11). The chopper amplitude was changed during the observations and the contaminated scans were not used for the final release.

The half power beam width (HPBW) is around ≈ 14−15′′ for the [C II] observation and ≈ 19−20′′ for [N II]. The exact value of the GREAT HPBW during the observation is unfortunately unknown. We use θ =14.4′′ as a tentative HPBW, following the SOFIA/GREAT Observation Planning (Version 9, April 29, 2016²).

The data were calibrated following Heyminck et al. (2012) and Guan et al. (2012). The reduced spectra shown in this study are in local standard of rest (LSR) velocity (Fig. 2). Spectra were generated from the level 3 GREAT product using the GILDAS/CLASS package. A baseline of first order was subtracted. Calibration from observed counts to main beam temperature T_mb was calculated with a beam efficiency 0.67 for L1 and 0.65 for L2. The spectral resolution is 0.15 km s⁻¹ for passband L#2 and 0.20 km s⁻¹ for L#1, but spectra were rebinned to obtain a channel width of 1.2 km s⁻¹ to increase the signal-to-noise ratio (S/N). For the conversion to Janskys we use the following equation: S[Jy] = 721 T_mb[K] (GREAT science team; private communication). The flux calibration uncertainty is about 10% (Guan et al. 2012). The final spectra are shown in Fig. 2.

2.2 Herschel/PACS

Herschel/PACS observations were taken as part of the SHINING key program (KPGT_esturm_1, PI: E. Sturm). Small maps were performed with the unchopped scan mode, for regions LMC-N 11 A, B, C, D, and I in [C II], [O I] 63μm, and [O III] 88μm (Cormier et al. 2015). The region N 11 B was also observed in [O I] 145μm, [N II] 122μm, [N II] 205μm, and [N III] 57μm (Lebouteiller et al. 2012). The region N 11 “South” is the only GREAT pointing (#12) with no PACS spectroscopy. All PACS lines appear unresolved given the spectral resolution (55–320 km s⁻¹; Lebouteiller et al. 2012).

The comparison between the [C II] fluxes measured from GREAT individual pointings and from the PACS map is contingent upon the knowledge of the GREAT beam profile and the flux calibration for both instruments. We use θ =14.4′′ for the GREAT HPBW (see Sect. 2.1). For a Gaussian beam profile, the beam solid angle is then 1.13 θ², corresponding to a ≈ 17′′ diameter. Therefore, we convolve the PACS map to a resolution of 14.4′′ and integrate the flux in an aperture of 17′′ diameter. The same steps are performed for the [N II] observations.

The fluxes measured with PACS and GREAT are shown in Table 3 for [C II] and [N II] 205μm. The agreement is good overall for [C II], although the PACS fluxes seem to be lower by a factor of 1.5 on average (see also Fahrion et al. 2017; Schneider et al. 2018). A better agreement could be reached if the effective GREAT HPBW were somewhat larger (≈ 20′′ solid angle). Furthermore, the PACS spectrometer data are calibrated for extended emission, so if the emission is indeed extended, the derived flux will scale with the aperture size; if the emission is point-like however, the total flux of the point source is enclosed in an area somewhat larger than the FWHM of the instrument, and the enclosed energy could be typically underestimated by 10− 15%³ (see also Fahrion et al. 2017).

2.3 H I 21 cm

The H I 21 cm data cube is taken from the ATCA+Parkes observations by Kim et al. (2003), which is, at the time of publication, the highest spatial resolution (60′′) H I survey of the LMC. The pixel size is 40′′ and the velocity channels are 1.6 km s⁻¹. The column density is calculated using a conversion factor for optically thin gas 1.823 × 10¹⁸ K km s⁻¹ cm⁻² (Dickey & Lockman 1990).

Optically thick H I may be significant and an important contributor to the dark neutral medium (DNM) gas as proposed by some authors (e.g., Fukui et al. 2015). We account for optically thick H I by using the results of Braun (2012) for the LMC. The opacity correction is based on a method described in Braun et al. (2009) which computes the temperature as well as the turbulent broadening on scales of 100 pc assuming a spatially resolved isothermal feature (i.e., single temperature along the line of sight). The total column density is then calculated using the profile fit parameters and a correction for the residual emission assumed to be optically thin. To estimate the factor by which we should correct the H⁰ column density from Kim et al. (2003), we use the smooth variation of the ratio between the velocity-integrated opacity-corrected and -uncorrected (optically thin assumption) H⁰ column densities. Figure 3 shows this ratio for the N 11 data points. The dispersion in the ratio is partly due to the opacity correction method (e.g., not considering self-absorption) and to the fact that multiple components exist along the line of sight with different opacities. For the H⁰ column densities measured in individual components (Sect. 4), the correction due to optically thick H I is less than a factor of two (Fig. 3). We consider this correction as an upper limit since other methods to infer the fraction of optically thick H I can result in much lower correction factors (e.g., Lee et al. 2015; Nguyen et al. 2018; Murray et al. 2018).

Fig. 2 Spectra of the main spectral lines used in this study for each GREAT pointing. From top to bottom in each panel: Hα, [Ne III] 3868Å, H I 21 cm, [N II] 205μm, [C II] 158μm, and CO(1–0). The histogram shows the data and the red curve shows the Gaussian fit. For this plot, only one component is considered except for pointing #5. Red dashed vertical lines show the Gaussian fit centroids while black dashed vertical lines show the zero velocity. The spectra are shifted arbitrarily along the y-axis for display purposes.

Fig. 3 Ratio of the opacity-corrected (Braun 2012) to observed H0 column densities across the N 11 region. The observed H0 column density for individual velocity components is always <2 × 1021 cm−2 (see Sect. 4.2).

2.4 CO (1–0)

The CO(1–0) observations are taken from the MAGMA survey described in Wong et al. (2011). Updated reduction is explained in Wong et al. (2017). The original spatial resolution is 45′′ and we use a 15′′ pixel size grid. Velocity channels are 0.53 km s⁻¹. The column density is calculated using a conversion factor⁴ cm⁻² (K km s⁻¹)⁻¹, which is the fiducial value for resolved CO clumps without the contribution of the CO-dark H₂ gas between CO clumps as discussed in Roman-Duval et al. (2014). In Table 4 we provide the CO intensity, the H₂ column density using , N(H₂|CO), and the molecular gas fraction f(H₂|CO) defined as 2 × N(H₂|CO)∕[N(H⁰) + 2 × N(H₂|CO)] (i.e., ignoring the CO-dark H₂ gas).

We also obtained ALMA band-3 observations of N 11B (programs 2012.1.00532.S and 2013.1.00556.S; PI: Lebouteiller). The mosaic covers the area observed with Herschel/PACS map (Sect. 2.2). The spatial resolution is 2.2′′, with a maximum recoverable scale of 12′′ and 47′′ for the 12- and 7-m arrays respectively. The velocity resolution is 0.64 km s⁻¹. The line width ranges between ≈2 and ≈ 5 km s⁻¹. The 7-m, 12-m, and Total power observations were combined for the final data cube. The velocity profiles for the broadest lines are asymmetric, suggesting the presence of multiple components (Fig. 2).

Figure 4 shows the original ALMA map and the projections on the Herschel/PACS [C II] and MAGMA CO(1–0) grids for comparison. While pointings #2 and #3 coincide well with compact CO peaks observed with ALMA, pointing #4is somewhat offset by a few arcseconds with respect to the ALMA clump. Given the GREAT beam size (14.4′′; Sect. 2.1), part of the spatially offset CO cloud emission is expected to contribute to the GREAT pointing #4.

Figure 4 also shows the comparison between MAGMA and ALMA spectra. Overall, there is good agreement between the two datasets over a 45′′ resolution, with the ALMA spectra showing a drastic S/N improvement. We also show the ALMA spectra calculated for a 14.4′′ resolution (i.e., similar to GREAT). On the one hand, ALMA spectra at 14.4′′ and 45′′ resolution are similar for pointing #2, confirming that most of the emission originates from a cloud smaller than 14.4′′ with no significant contamination from nearby clouds. On the other hand, the ALMA spectra toward pointings #3 and #4 show some differences with spatial resolution, indicating that nearby clouds with different properties (intensity, central velocity, line width) contribute to the spectrum calculated over a 45′′ beam for these pointings. Pointing #5 shows no CO emission with ALMA over 14.4′′ resolution but emission is seen for the 45′′ resolution which is likely arising from the main PDR complex corresponding to pointing #2. In the following, the ALMA spectra are used when available in order to compare to the [C II] spectral profiles. The detailed analysis of the ALMA map is deferred to a future work.

Fig. 4 Comparison between ALMA, MAGMA, and Herschel data. Left column, from top to bottom: Herschel/PACS [C II] map of N 11B, the ALMA CO(1–0) map (original resolution, convolved to 14.4′′ and projected on PACS [C II] grid, and convolved to 45′′ and projected on MAGMA grid), and the ratio (log units) of CO MAGMA/ALMA. The GREAT pointings are overlaid as 14.4′′ diameter circles. Right column: MAGMA spectrum (blue histogram) and the ALMA spectra calculated for a 45′′ resolution (red with shade) and 17′′ resolution (black) for the four pointings in N 11B.

2.5 VLT/GIRAFFE

Optical spectroscopy is used to examine tracers of the ionized gas. We use archival ESO/VLT/U2 GIRAFFE observations with the MEDUSA fiber system feeding component. GIRAFFE is a medium-high (R = 5500−65 100) resolution spectrograph in the 3700−9500Å wavelength range. The MEDUSA fibers allow up to 132 separate objects (including sky fibres) to be observed in one go. Each fiber has an aperture of 1.2′′ on the sky.

The sky-subtracted spectra from program 171.D-0237 were downloaded from the GEPI GIRAFFE archive. We examined in particular Hα observed with the HR14a grating with R = 18 000 (≈ 17 km s⁻¹), and [Ne III] 3868Å observed with the HR2 grating with R = 19 600 (≈ 15 km s⁻¹). Some details on the dataset can be found in Torres-Flores et al. (2015).

Most GREAT positions could be associated with a nearby GIRAFFE fiber position, except #12. When several field datasets were available for a given GREAT pointing, we selected only those with the best seeing and S/N. The [Ne III] and Hα spectra are shown in Fig. 2. The Hα line is well detected toward all pointings, while [Ne III] is detected everywhere except in N 11 D and I (GREAT pointings #9 and #10-11 resp.).

3 Integrated properties along the lines of sight

In this section we examine the kinematic properties (radial velocity and line width) integrated along the lines of sight. Results using individual velocity components are discussed in Sect. 4.

3.1 Comparison of the spectral profiles

The CO(1–0) MAGMA spectra can all be fitted reasonably well with a single resolved (i.e., wider than the spectral resolution; Fig. 2) component, with FWHM ≈ 3–8 km s⁻¹ for the lines with the best S/N (see Table 5). The CO(1–0) line was also observed with ALMA in the N 11B region (pointings #2 through #5), which once convolved to 45′′ resolution shows good agreement with the MAGMA data (Sect. 2.4) but also highlights some asymmetry in the line profile (Fig. 4).

Most [C II] GREAT spectra can also be fitted with a single resolved component, although a few pointings show evidence of multiple components (most notably #3, #4, #5, #6, and #8) or an asymmetric profile (e.g., #7, #9). The FWHM of [C II] measured with a single wide component (≈ 4–10 km s⁻¹) is always either similar to or larger than the CO FWHM. The H I profile is always much broader than the CO and [C II] profiles, with FWHM in the range ≈ 16–40 km s⁻¹. Similar results were found in various star-forming regions within nearby galaxies (Braine et al. 2012; de Blok et al. 2016; Requena-Torres et al. 2016; Okada et al. 2015; Fahrion et al. 2017).

We verified that the difference in the profile line width between H I and CO is not due to the spatial resolution by matching both datasets. However, the wider H I profile as compared to [C II] is likely driven by the difference in spatial resolution. The H I spectra, taken with 60′′ resolution (Sect. 2.3), include the emission of clouds outside the GREAT beam. The profile decomposition will allow us to mitigate the biases due to different beam sizes (Sect. 4). Whether from the spectral profiles (Fig. 2) or from the images (Fig. 4), there is no evidence of CO emission not associated with [C II].

Since [N II] was not detected with GREAT (Sect. 2.1), we attempted to use optical tracers (Sect. 2.5) to describe the kinematics of the ionized gas. The observed FWHM of [Ne III] is ≈ 22 km s⁻¹, while it is ≈ 30−40 km s⁻¹ for Hα. The difference between the two line widths is mostly due to the thermal broadening. The thermal broadening (FWHM) in the ionized gas is (1)

with T the temperature and m the ionic mass. For a temperature of 10 000 K (Toribio San Cipriano et al. 2017), we find ≈ 22 km s⁻¹ for H⁺ and ≈ 4.8 km s⁻¹ for Ne²⁺. Convolutions with the instrumental resolution (Sect. 2.5) and with macro-turbulent motions (typically a few km s⁻¹) yield ≈ 30 km s⁻¹ for Hα and ≈ 16 km s⁻¹ for [Ne III]. Both [Ne III] and Hα thus appear somewhat resolved, with a kinematical, nonthermal width of ≈ 15−25 km s⁻¹, which likely reflects the dispersion of several individual components along the line of sight. The only GIRAFFE pointings showing a noticeable velocity structure in [Ne III] are the ones associated with pointings #2, #4, and #5 in N 11 B. Overall, the spectral resolution in the optical tracers is unfortunately not sufficient to decompose the velocity profiles but the fact that several velocity components may contribute to the observed ionized gas tracer profiles provides a useful constraint for the origin of [C II] based on its line width (Sect. 3.2.1).

3.2 Contribution from the ionized gas to [C II]

The [C II] line may originate in the neutral (atomic or CO-dark H₂) or ionized gas. Unfortunately, no direct comparison of the velocity structure in [C II] and ionized gas tracers can be performed because the [N II] 205μm line was not detected with GREAT (Sect. 2.1) and because the spectral resolution for optical tracers is not high enough (Sect. 3.1).

3.2.1 [C II] line width

The expected thermal broadening⁵ for C⁺ in the ionized gas is ≈ 6 km s⁻¹ (FWHM for a temperature of 10 000 K). Macro-turbulent motions (typically a few km s⁻¹) may increase this value to ≈ 10 km s⁻¹, but we also keep in mind that the temperature in the WIM, where a significant fraction of C⁺ could exist, may be colder than 10 000 K. Therefore, we expect [C II] line widths in the ionized gas for a single velocity component to be in the range ≈ 6−10 km s⁻¹.

The observed FWHM for the main [C II] velocity component is smaller than ≈ 6 km s⁻¹ toward pointings #3, #10, #11, #12, and for one of the two components toward pointing #5 (Table 5). We conclude that [C II] is necessarily associated with the neutral gas for these components. This conclusion is strengthened if assuming the main [C II] velocity component in reality corresponds to multiple, narrower components. For some of the other pointings, the [C II] profile indeed shows multiple blended components (#4, #6) or asymmetric profiles (#7, #8, #9), suggesting that the presence of multiple components drives the large total width (Sect. 3.1). The presence of multiple components in the ionized gas is also suggested by the wide [Ne III] profiles in the optical (Sect. 3.1).

Pointing #5 is particularly interesting as it shows two well separated [C II] components, with FWHM of 4.0 and 7.0 km s⁻¹ for the low- andhigh-velocity components respectively. While the FWHM of the low-velocity component is compatible with an origin of [C II] in the neutral gas, the high-velocity component may come from the ionized gas provided it corresponds to a single cloud.

In summary, the observed FWHMs for [C II], together with the fact that [C II] in the ionized gas may originate from more than one component, suggests that the observed [C II] components are unlikely to arise from the ionized gas, except maybe for one component toward pointing #5, and perhaps pointing #1. The line width of all the other [C II] components is compatible with an origin in the neutral gas, where macro-turbulence competes with the thermal broadening. Other complementary methods are examined in the following.

3.2.2 Ionized gas traced by [N II] 122μm and 205μm (pointings #1 through #6)

Owing to the energy range required to produce N⁺ and C⁺ ions (14.5−29.6 and 11.3−24.4 eV respectively) and owing to similar critical densities for [C II] and [N II] for collisions with e⁻, [N II] lines are often used to trace [C II] emission in the ionized gas (e.g., Oberst et al. 2006). We can then compare the observed [N II]/[C II] ratio to the theoretical value in the ionized gas, with any deviation indicating a contribution from [C II] in theneutral gas. In order to compute the theoretical ratio, we first calculate the ionic abundance ratio N⁺ /C⁺ using the solar abundance ratio logN/C = −0.6 (Asplund et al. 2009; also valid for LMC-30 Dor Pellegrini et al. 2011) together with the ratio [(N⁺/N)/(C⁺/C)] = 1.26 computed from MAPPINGS III photoionization grids (Sutherland et al. 2013). We then need to account for the gas density. While the 205μm line has a critical density (180 cm⁻³ for collisions with e⁻) close to that of [C II] 158μm (≈30−50 cm⁻³; Goldsmith et al. 2012), the 122μm line has a larger critical density (400 cm⁻³), meaning that the [N II] 122μm/[C II] ratio depends somewhat more on density. Figure 5 shows the theoretical [N II]/[C II] ratio as a function of density including the ionic abundance correction factor.

Using the upper limits on the [N II] 122 and 205μm lines observed towards N 11B with Herschel/PACS (corresponding to pointings #2, #3, #4, and #5; Sect. 2.2), we find [N II] 122μm/[C II] ≲ 0.05 and [N II] 205μm/[C II] ≲ 0.01 except toward #5 where [N II] 205μm/[C II] < 0.06 (Lebouteiller et al. 2012). The comparison with the theoretical ratio in the ionized gas for typical densities between 10¹⁻³ cm⁻³ (Fig. 5) indicates that the fraction of velocity-integrated [C II] originating in the neutral gas is f_n ([CII]) ≳ 70% using [N II] 122μm and ≳ 90% using [N II] 205μm (except for pointing #5 with ≳55%). Therefore, for pointings #2, #3, and #4, the fraction of [C II] in the ionized gas is not significant. For pointing #5(made of two distinct velocity components), the most stringent constraint is obtained with [N II] 122μm/[C II], with f_n([CII]) ≳ 70% for the sum of both components. Since the two [C II] components toward #5 have similar fluxes, this suggests that the fraction of [C II] originating in the ionized gas must be relatively small for both. The [N II] 205μm upper limits calculated from the GREAT spectra (available for #1 through #6) do not allow us to improve further the determination of f_n ([CII]).

We can also model the velocity profile of [N II] 205μm with GREAT assuming [C II] originates fully from the ionized gas and assuming the same velocity structure (i.e., ignoring potential differences in individual component line widths). Our results are illustrated in Fig. 6 in which we also show the brightest ionized gas [C II] emission possible using the [N II] 205μm upper limits. Thedensity is chosen as 100 cm⁻³ for this test. It is clear that [C II] originates overall from the neutral gas. For pointing #5, the simulated [N II] profile is compatible with the observed profile within uncertainties but we recall that the results using velocity-integrated measurements with Herschel/PACS suggest that the [C II] contribution from ionized gas is not dominant. We can also conclude that relatively faint [C II] components observed toward #3, #4, or #6 cannot be explained by [C II] emission in the ionized gas.

Fig. 5 Theoretical [N II]/[C II] ratios as a function of density. The ratio values include ionic abundance ionization corrections (see text). The red curve indicates ratio values with [C II] originating fully in the ionized gas. The solid black lines show the observed upper limits for all pointings observed with Herschel/PACS (except #5 for [N II] 205μm/[C II] shown as a dashed line).

3.2.3 Photoelectric heating efficiency proxy

Finally, the origin of [C II] can be examined indirectly by comparing the neutral atomic gas cooling traced by the [C II] and [O I] lines to the gas heating. The latter can be traced by far-infrared emission but contamination by warm dust in the ionized phase can become an issue (Lebouteiller et al. 2012). The PAH emission is another tracer of the neutral gas heating, and it is plausible that PAHs actually dominate the gas heating as compared to very small grains. In LMC-N 11B, Lebouteiller et al. (2012) found that the ratio ([C II]+[O I])/PAH is uniform, indicating that the photoelectric heating efficiency is fairly constant (see also Helou et al. 2001; Croxall et al. 2012; Okada et al. 2013). The ratio remains constant even toward the stellar cluster LH 10 (corresponding to pointing #5) where ionized gas dominates the infrared lineemission, suggesting that [C II] and [O I] still trace the neutral atomic gas in the foreground and background.

We now revisit this finding by extending the measurement of ([C II]+[O I])/PAH to all PACS maps in LMC-N 11. Since there is no mid-IR spectra available for all pointings in N 11, we use the IRAC photometry bands to evaluate the PAH emission. We refer to Lebouteiller et al. (2012) for the method. In short, we use the Spitzer/IRAC photometry data points to identify and flag out stellar emission, to estimate and correct for the dust continuum, in order to provide the emission of PAHs either in the IRAC 5.8μm band (corresponding to the PAH 6.3μm feature) or in the IRAC 8.0μm band (corresponding to the PAH 7.7 + 8.6μm features).

Results are shown in Fig. 7, where it can be seen that the cooling rates provided by [C II] and [O I] compensate, with [C II] dominating in the faintest, likely more diffuse, regions, and with [O I] dominating in the brightest, more compact, regions. The sum [C II]+[O I] is however proportional to the PAH emission with no obvious trend with PAH brightness. The remarkably small scatter in ([C II]+[O I])/PAH in many regions including PDRs (#2) suggests that there is no significant extra [C II] emission arising from the ionized phase. Potential PAH emission in the ionized gas is ignored as it is plausible that they are photodestroyed (see e.g., Madden et al. 2006; Lebouteiller et al. 2007; Chastenet et al. 2019).

As for the other diagnostics using velocity-integrated values, this result does not preclude the presence of relatively faint [C II] components that could be associated with the ionized gas. Overall, the results obtained with ([C II]+[O I])/PAH, together with those obtained in Sects. 3.2.1 and 3.2.2, suggest that most of the [C II] arises in the neutral gas. We find no velocity component that can be unambiguously associated with ionized gas.

Fig. 6 Comparison of the GREAT [N II] and [C II] spectral profiles for the six pointings with [N II] observations. The red solid curve shows for both lines the Gaussian fit to the main velocity components as in Fig. 2. For the [N II] spectrum the black curve corresponds to the rms flux while the red dotted curve shows the expected [N II] profile assuming [C II] arises only from the ionized gas, with the uncertainty shown as the gray area (see text). For the [C II] spectrum the red dotted curve and shaded area show the profile of [C II] in the ionized gas calculated using the [N II] rms (i.e., the brightest possible [C II] component arising from the ionized gas).

Fig. 7 Variation of [O I]/PAH (bottom), [C II]/PAH (middle), and ([C II]+[O I])/PAH (top) across the regions in LMC-N 11 mapped with PACS. The size of the symbol scales with the signal-to-noise ratio. Black circles correspond to GREAT pointings. The solid line shows the linear regression for the corresponding ratio while the dashed lines show the regression obtained for the two other ratios.

4 Individual component analysis

In this section we attempt to decompose the CO, [C II], and H I profiles toward the GREAT pointings in order to derive the physical conditions associated with each individual component, with the help of simple two-phase models (neutral atomic and molecular). We assume that the ionized contribution to [C II] is negligible (Sect. 3.2).

4.1 Profile decomposition

We assume that the velocity structure reflects the presence of several individual components, which we refer to as “clouds”, defined by their radial velocity. We ignore clouds with potentially significant velocity gradient along the line of sight. In other words, all tracers (CO, [C II], H I) are assumed to peak at the same velocity for any given cloud. Some caveats pertain to the different spatial resolutions of the observations. The H I emission in particular may arise from clouds that are not contributingto the GREAT beam (Sect. 3.1). However, this confusion is mitigated when individual velocity components are considered, since these are more likely to originate from a relatively more spatially confined cloud.

The velocity profiles of CO, [C II], and H I can be decomposed in various ways, for instance by fitting Gaussian components for the tracer with the narrowest profile, typically CO, and using these centroids for the adjustment of other tracers with other velocity components added as needed (e.g., Okada et al. 2019). We decided to use a statistical approach instead, by adjusting the various profiles simultaneously with velocity components not being fixed or inferred from any other specific profile. Any given component is defined by its velocity and line width, while the component intensity is determined independently for each CO, [C II], and H I. The sum of the components reproduces the global CO, [C II], and H I profiles. We also performed decompositions allowing potentially different line widths for each tracer for any given velocity component, but the main results shown in the following remain unchanged.

We use a Bayesian approach with the Markov Chain Monte Carlo (MCMC) PyMC3 code (Salvatier et al. 2016) along with the Metropolis-Hastings sampling algorithm. Parameters are described in Table 6. We identify three parameters that control the number of possible solutions, namely the number of components, the minimum line width, and the minimum separation between components. Although we could narrow down the number of solutions by introducing these parameters within the Bayesian inference (eventually finding a unique converged solution), other degenerate solutions may also be acceptable. Since we consider that we do not have enough observational constraints and priors for these parameters, we chose to test different fixed values:

• The number of components is set to 10 (minimum to reproduce the H I profile) and 15 (to allow potential significant blends).

• The minimum line width is set to 1 and 2 km s⁻¹. The velocity profiles in Fig. 2 show that, when the [C II] profile is visibly made of more than one component (e.g., #3, #4, #6), the [C II] component associated with CO emission shows the same line width as the CO component, with FWHM values ≳3 km s⁻¹; Table 5). This suggests that macro-turbulence dominates the [C II] and CO line width. Therefore, we impose the same line width for a given velocity component in each tracer and ignore line width differences due to mean molecular weight. The lowest value of 1 km s⁻¹ is motivated by the typical turbulence measured for interstellar clouds (e.g., Welty et al. 1994, 1996).

• The minimum component separation is set to 0 km s⁻¹ (components with same velocity but potentially different line width), 1 and 2 km s⁻¹ (larger values resulting in unsatisfactory fits).

We note that the converged solution for the decomposition, like any other approach, is not claimed to represent the actual velocity structure, but we hope to infer statistically representative trends by making use of the dispersion in the Markov chain (related to flux uncertainties in the spectra), by considering different models (i.e., with different input parameters), and by considering the various pointings. The MCMC method draws a sequence of random variables corresponding to each parameter in the decomposition. Therefore, there are some excursions, especially in the burn-in phase (first iterations, or “states”, in the chain). Once the burn-in phase is finished, the chain enters a high-probability region where the result has converged. The chain still contains a random distribution though, with excursions according to the uncertainties in the parameter. The standard deviation of the Markov chain therefore provides a probability density function. Figure 8 shows an illustration of the simultaneous fit of all the components. The profile decomposition for the other pointings can be found in Appendix A.

4.2 Model-independent quantities

For each velocity component we compile the [C II] line intensity, the H⁰ column density, the H₂ column density derived from CO, N(H₂|CO), and the resulting CO-traced molecular gas fraction f(H₂|CO) = 2N(H₂|CO)∕(2N(H₂|CO) + N(H⁰)). The CO-dark H₂ gas will be accounted for later with the model (Sect. 4.3).

We wish to stress that beam dilution may affect H I and CO observations differently if the respective beam filling factors are different. The beam filling factor for CO is likely much smaller than for H I, and as a result the value of f(H₂|CO) should correspond to the beam filling factor of fully molecular clouds embedded in a mostly atomic medium rather than to the actual molecular gas fraction of a single cloud filling the beam.

The decomposition results are illustrated in Fig. 9, where we show the bivariate kernel density estimate (non-parametric probability density function) of the [C II] emission and f(H₂|CO). The kernel density estimation makes use of the data cloud corresponding to every velocity component of every pointing for all the elements in the Markov chain, and for all model input parameters (number of components, minimum component line width, and minimum component separation). Hence the kernel density estimate conveniently indicates whether different solutions to the profile decomposition cluster around similar locii in the parameter space. Most of the components in Fig. 9 seem to lie either at f(H₂|CO) ≲ 10% or f(H₂|CO) ≳ 60%, suggesting a sharp transition between CO-bright H₂ gas and either CO-dark H₂ or atomic gas. Components with a low f(H₂|CO) are our best candidates for evidence of significant CO-dark H₂ gas amount.

4.3 CO-dark gas properties

4.3.1 Selection of components

Our models account for the neutral gas only (atomic and molecular). The ionized gas contribution to the integrated line emission is negligible but its contribution to specific velocity components, in particular faint ones, is more difficult to assess (Sect. 3.2). An upper limit can be calculated for [C II] in the ionized gas for the pointings with velocity-resolved [N II] 205μm observations (Sect. 3.2.2; Fig. 6). In the following, we ignore all the [C II] components below 2 × 10⁻⁸ W m⁻² sr⁻¹ (corresponding to ≈2σ) to ensure we select [C II] components arising in the neutral gas. Although such [C II] components are faint, and although they do not contribute much to the kernel density estimate plots, we prefer to ignore them for clarity and robustness of the final results. This threshold also prevents over-interpretation of potential [C II] upper limits corresponding to H I-identified components that do not contribute to the GREAT beam.

Fig. 8 Illustration of the profile decomposition method for pointing #2 for one of the model input parameter sets (Nc = 10; σmin = 1 km s−1; Δvmin = 1 km s−1). The gray histogram shows the data. The red, green, and dark blue curves correspond to different components and the dotted vertical line shows the velocity for each component. The black solid curve shows the total inferred profile.

Fig. 9 Bivariate kernel density estimate (non-parametric probability density function) of I([CII]) vs. the molecular gas fraction ignoring CO-dark H2 gas, f(H2|CO), for all pointings. The shade scales with the density of points.

4.3.2 Modeling strategy

While the fraction of [C II] associated with the atomic or molecular phase could be derived from the general profile shapes (e.g., Okada et al. 2019), we choose here a different approach making use of the CO and H⁰ column densities derived for each velocity component. Since we lack additional constraints (e.g., [O I], [C I], or the infrared luminosity) for the individual velocity components, we cannot derive the physical parameters from PDR models such as the impinging radiation field intensity or the cloud extinction for each individual component. Indirect measurements of such parameters making use of 2D projected quantities are themselves subject to caution (see Seifried et al. 2019 for the extinction). Here we rely on a simple agnostic approach to calculate the level population of C⁺ accounting for collisions with H⁰, H₂, and e⁻ as a function of gas temperature and density (see details in Lebouteiller et al. 2013). The clouds that are identified thanks to the velocity decomposition methods are expected to contain both atomic gas, molecular gas traced by CO, and molecular gas traced by C⁺. For a given velocity component, the CO-dark H₂ gas may be related to a CO clump in which case the cloud emits both in [C II] and CO. The CO-bright H₂ gas is not considered in the model because we only calculate the properties of [C II]-emitting gas, and as such, the H₂ column density associated with CO is simply calculated using a fiducial factor (Sect. 2.4).

The strategy goes as follows (illustrated in Fig. 10):

• Step 1: for each component we calculate a range for the atomic gas number density and the [C II] emission is estimated in the atomic phase using the atomic gas number density and column density.

• Step 2: the [C II] intensity from the molecular gas is inferred from the difference of the estimated [C II] emission in the atomic phase with the observed [C II].

• Step 3: we compare the [C II] emission in each phase to compute f_coll,H2([C II]), which is the fraction of [C II] associated with gas where collisions with H₂ dominate, that is, the CO-dark H₂ gas.

• Step 4: from the [C II] intensity from the molecular gas, we calculate the H₂ column density traced by C⁺, N(H₂|C⁺).

For step (1), the number density n(H⁰) in the atomic gas is calculated using the H⁰ column density and assuming that the individual cloud size along the line of sight lies in the range 1− 10 pc⁶. Constraints on the cloud size are mostly given by the H I morphology and typical cloud sizes observed with ALMA in 30 Dor (Indebetouw et al. 2013) or N 11 (this study). The inferred density n(H⁰) therefore ranges from a few cm⁻³ to ~ 10³ cm⁻³. The temperature in the atomic phase is fixed assuming that all H I components that can be associated in velocity with [C II] trace the cold rather than the warm neutral medium (CNM, WNM resp.; see also Pineda et al. 2013, 2017). Accordingly, we estimate a temperature of 100 K from [O I] 63μm/[C II] (Appendix C). The resulting pressure is then about 10^2.5−5 K cm⁻³, which is similar to what is found in predominantly atomic gas within Herschel KINGFISH galaxies (Herrera-Camus et al. 2017) and to other regions in the LMC (Okada et al. 2019) using far-infrared lines. The ionization fraction is fixed to a value of n_e ∕n_H = 10⁻⁴, that is, a typical value for a UV-illuminated diffuse gas in which free electrons are provided by ionization of species with an ionization potential below 13.6 eV (i.e., no significant ionization from cosmic rays or X-rays).

For step (3), the main free parameter is the H₂ column density traced by C⁺, N(H₂|C⁺). We use two constraints for the molecular component: (i) there is no lower limit on the molecular cloud size, but we use the same upper limit as for the atomic gas emission (10 pc), and (ii) the number density in the molecular phase scales by default with that in the atomic phase. For the second hypothesis, the scaling assumes by default thermal pressure equilibrium. The expected temperature of the CNM and of the CO-dark H₂ gas is about 10−100 K in Milky Way conditions (Glover & Clark 2016; Tang et al. 2016; Seifried et al. 2019). A somewhat warmer temperature is expected in metal-poor environments due to the larger photoelectric-effect heating efficiency (itself due to the low dust content), to the lack of metal coolants, and to the harder interstellar radiation field (ISRF). We assume in the following a temperature of 50 K in the molecular phase but our results are not changed significantly if we take 100 K. At thermal pressure equilibrium, the density in the H₂ gas is therefore twice larger than in the H⁰ gas, reaching up to ≈ 10³ cm⁻³, on the low end of values found in Pineda et al. (2017) for a sample of other LMC star-forming regions. We allow the molecular gas to be denser if no solution can be found with a thermal pressure equilibrium hypothesis. This is motivated by the fact that pressure equilibrium might not always be satisfied (e.g., Vázquez-Semadeni et al. 2000; Scoville 2013) and that thermal pressure may depend on the radiation field intensity and star-formation activity (e.g., Ostriker et al. 2010). Furthermore, the estimated densities correspond to averages along the line of sight while the density structure of the neutral atomic phase and CO-dark H₂ gas may differ significantly.

From these results, we can then calculate the following quantities for each velocity component:

• The total H₂ column density,

(2)

• The fraction of CO-dark H₂ gas,

(3)

• The mass fraction of molecular gas,

(4)

We note that the mass of CO-dark H₂ gas we measure is a lower limit because we consider only the gas traced by C⁺. In principle however, other tracers with possibly different properties may also correspond to CO-dark H₂ (Glover & Clark 2016; Clark et al. 2019). In Fig. 11 we show the possible model values for f_dark as a function of [C II]/CO for the observed ranges of N(H⁰), N(H₂|CO), and I([CII]). The correlation between f_dark and [C II]/CO is tighter for low densities, corresponding to low H⁰ column densities (due to the assumed cloud size), and therefore to a low fraction of [C II] in the neutral atomic gas. For such conditions, the [C II]/CO ratio is proportional to first order⁷ to the column density ratio N(C⁺)∕N(CO) (e.g., Crawford et al. 1985). As the atomic hydrogen column density increases, so does the possible contributionof the neutral atomic gas to [C II] (shifting the observed [C II]/CO to larger values) and so does the number density in the atomic and consequently the molecular gas (leading to lower f_dark values). The dependency of f_dark with temperature in the atomic gas is relatively small unless the temperature is much larger than ~ 500 K. Even then, if we were to choose a temperature of several thousand K (WNM conditions), f_dark values would be lower by less than about 10% in value.

An illustration of the model calculation is shown in Fig. 12. In practice, densities in the CO-dark H₂ gas lie around 10^1.5−3.5 cm⁻³, corresponding to a pressure of 10^3.5−5 K cm⁻³ (Fig. 13). The total H₂ column density N(H₂) measured this way varies across pointings from ≈ 10²¹ cm⁻² (in #5) to ≈2 × 10²² cm⁻² (in #10).

Fig. 10 Illustration of the model strategy. Two phases are considered for which the level population of C+ is calculated as a function of density n, temperature T, and the ionization fraction x in the atomic medium. Collisional rates with each partner also depend on temperature. The [C II] intensity is calculated using the level population and the column density of H0 or H2 and the fraction of C into C+ in each phase (see Lebouteiller et al. 2013 for details).

Fig. 11 Possible model values for fdark as a function of [C II]/CO and volume density in the neutral atomic phase (color bar). The models shown are those corresponding to the observed ranges for N(H0) = [1018, 1021.3] cm−2, N(H2 |CO) = [1018, 1021.5] cm−2, and I([CII]) = [10−8, 10−6] W m−2 sr−1. Model results are shown for a temperature in the neutral atomic phase of 100 and 50 K in the molecular phase. The striping is due to the incomplete coverage in the model results.

Fig. 12 Illustration of model results for one components in pointing #1 (top) and in #2(bottom). The [C II] surface brightness is calculated as a function of the H2 column density associated with C+. The horizontal red line shows the observed [C II] surface brightness. The purple and pink shaded areas show the range of model results assuming thermal pressure equilibrium and overdense molecular gas respectively (see text). The bottom and top boundaries are set by the H0 column density together with the cloud size constraint. The right boundary is set by the maximum allowed molecular cloud size. In the top panel the model with thermal pressure equilibrium can reproduce the observations while in the bottom panel, a model with an overdense molecular gas is required.

Fig. 13 Histogram of the modeled pressure in the neutral atomic (blue) and CO-dark H2 gas (red) forall the velocity components. Results from all decomposition methods are combined (see Table 6).

4.4 Results

4.4.1 [C II] in the neutral atomic gas

The [C II] components with low f(H₂|CO) values (Fig. 9) could a priori be either from CO-dark H₂ gas or atomic gas. When including the contribution from the CO-dark H₂ gas calculated by the models, most velocity components reach a large total (including CO-dark H₂) molecular gas fraction f(H₂) ≳ 80% (Fig. 14), except for the faintest [C II] components for which only an upper limit on f(H₂) is often available. This indicates that [C II]-bright regions are dominated by CO-dark H₂ gas while [C II]-faint regions may include a significant contribution from neutral atomic gas.

The weak correlation between f(H₂) and I([CII]) in Fig. 14 indeed suggests that the contribution of the neutral atomic gas to the [C II] emission is more significant toward faint [C II] components, which is confirmed by Fig. B.1 which shows the distribution of the fraction of [C II] tracing CO-dark H₂ gas, f_coll,H2([CII]), versus I([CII]), in particular for #9 and #11.

Figure 15 shows that the [C II] contribution from the neutral atomic gas does not occur preferentially toward components with low [C II]/CO. In this figure one can also clearly identify the two components in #5, which are both dominated by CO-dark H₂ gas. Most pointings seem to show either a wide distribution of [C II]/CO values or two main peaks, with the main peak at low [C II]/CO being due to the velocity component bright in CO and [C II], while the other peaks with large [C II]/CO values correspond to [C II] components with little associated CO emission.

Globally, we find that >95% of the [C II] emission across all pointings can be attributed to CO-dark H₂ gas (i.e., when f_coll,H2([CII]) > 50%), as shown in Fig. 16. This result does not depend on the CO column density and holds in particular for the brightest [C II] components associated with CO. The [C II] components that have the largest contribution from the neutral atomic gas (low f_coll,H2([CII])), correspond preferentially to the components with the faintest [C II] surface brightness and not to components with low CO column density (Fig. 16). Even then, less than ~ 30% of the [C II] emission in these low-[C II]-surface-brightness components may be arising from gas that is mostly atomic (i.e., f_coll,H2([CII]) < 50%). Okada et al. (2019) find similar results toward several star-forming regions in the LMC, with < 15% of [C II] being attributed to neutral atomic gas. We wish to emphasize that the larger fraction of [C II] in the neutral atomic gas in faint [C II] components is driven by the variation of the [C II] flux associated with the molecular phase.

In summary, [C II] mostly traces the CO-dark H₂ gas but there is evidence of a weak contribution from the neutral atomic gas in the faintest [C II] components. The thermal pressure in the CO-dark H₂ gas lies around 10³⁻⁵ K cm⁻³ (Fig. 13), which is similar to the range derived by Pineda et al. (2017) for several diffuse (with no CO detection) LMC regions. Considering the large filling factor of [C II] emission observed in N 11B with Herschel/PACS (Fig. 4), it is therefore plausible that the CO-dark H₂ gas probed by our [C II] observations corresponds to the diffuse ISM (median density ~ 200 cm⁻³ for all pointings) filling the space between CO clumps rather than thin CO-dark H₂ layers around CO clumps.

Fig. 14 Same as Fig. 9 but with the total molecular gas fraction (i.e., including the CO-dark H2 gas). The vertical striping is due to the wide ranges for the molecular gas fraction determination for the faintest [C II] components. The shade scales with the density of points.

Fig. 15 Bivariate kernel density estimate of the fraction of [C II] tracing CO-dark H2 gas vs. [C II]/CO for all pointings. The shade scales with the density of points.

Fig. 16 Cumulative fraction of [C II] plotted against the fraction of [C II] emission arising in the CO-dark H2 gas. Values are calculated with all pointings combined. Plots are shown for all components (red), for components with low CO column density (<2 × 1020 cm−2; dark blue), for components with high CO column density (>1 × 1021 cm−2; blue), and for components with weak [C II] (<4 × 10−8 W m−2 sr−1; green).

4.4.2 Fraction of CO-dark H₂ gas

We show in Fig. 17 the distribution of the fraction of CO-dark H₂ gas, f_dark, versus [C II]/CO for each pointing. The fraction of CO-dark H₂ gas is also a proxy for the fraction of DNM (i.e., CO-dark H₂ gas and optically-thick H I) since the CO-dark H₂ gas traced by [C II] is the main contribution (≳95%) to the DNMaccording to our models. The weak contribution of optically thick H I to the DNM is also found in the local and diffuse ISM (Murray et al. 2018; Liszt et al. 2018).

Most pointings show a well-defined peak in Fig. 17 corresponding to the brightest [C II] component. Fainter [C II] components often result in a wide range of f_dark values because they are more likely to arise from the neutral atomic phase (with relatively large uncertainty; Sect. 4.4.1). From Fig. 17 we see that most of the molecular gas is CO-dark overall, with f_dark ≳ 60%. Since most of the molecular gas is CO-dark, we conclude that the velocity-integrated ([C II]+[O I])/PAH ratio, which suggests a constant photoelectric-effect heating efficiency across N 11 (Sect. 3.2.3), corresponds in fact to the physical conditions of the CO-dark H₂ gas rather than the neutral atomic medium.

The f_dark values we find are in good agreement with the estimates of Galametz et al. (2016) in N 11 using the dust-to-gas mass ratio. Our results are also in general agreement with other studies in nearby low-metallicity environments. Fahrion et al. (2017) used SOFIA/GREAT observations the dwarf galaxy NGC 4214 at ≈ 200 pc spatial scale and found that only ≈ 10% of [C II] could be attributed to the CNM and that ≈80% of the H₂ mass is not traced by CO. Other studies in the Magellanic Clouds also show significant CO-dark H₂ gas fractions. Requena-Torres et al. (2016) examined several star-forming regions in the SMC, and found that most of the [C II] emission originates from CO-dark H₂ gas. Chevance (2016) observed LMC-30 Dor and determined from PDR models that ≳ 80% of the molecular gas is CO-dark H₂ gas.

In N 11, there is a clear correlation between f_dark and [C II]/CO (Fig. 17), which is not surprising because the CO-dark H₂ gas in the model is by constructiontraced by C⁺ and becausethe contribution of [C II] in the neutral atomic gas is relatively small (Sect. 4.3.2). For this reason, the [C II]/CO ratio we use in Fig. 17 is the observed value, i.e., with nocorrection for the [C II] in the atomic gas. For a given [C II]/CO ratio, the range of f_dark values is driven by the H⁰ column density which sets the density in the atomic gas and consequently in the molecular gas (Sect. 4.3.2).

Most points seem to follow the same relationship between f_dark and [C II]/CO (dashed curve in Fig. 17). For [C II]/CO ≳ 10^4.5, this is because f_dark depends little on the gas density (Sect. 4.3.2). For lower [C II]/CO values, f_dark depends more on density and components with similar [C II]/CO ratios may have significantly different f_dark values (e.g., from ≈75% for #9 to ≈90% for #5).

We also find that f_dark is anti-correlated with the CO column density (Fig. 18). Most of the molecular gas is therefore CO-dark, but there is more CO-dark gas in CO-faint regions (~100% for N(H₂|CO) ≲ 10²⁰ cm⁻²) as comparedto CO peaks (~70−90% for N(H₂|CO) ~ 10²¹ cm⁻²). This result is reminiscent of the findings of Okada et al. (2015) that the amount of [C II] that cannot be attributed to the gas traced by CO is larger between CO peaks in another LMC star-forming region, N 159.

The effective X_CO factor including the contribution of the CO-dark H₂ gas lies in the range 10²¹⁻²² (K km s⁻¹)⁻¹ for most of the bright velocity components (Fig. B.2), in good agreement with values obtained in Israel (1997), Galliano et al. (2011), and Chevance (2016).

4.4.3 Physical parameters controlling f_dark

We find that f_dark is somewhatlarger for components with a large [C II]/CO ratio (Fig. 17), which should correspond to molecular clouds illuminated by the UV radiation from massive stars, providing some evidence that stellar feedback may play a role. In Fig. 18, the pointings are ordered based on the shape of the f_dark with N(H₂|CO) distribution, from pointings with a high and almost flat f_dark distribution until column densities of ≈ 10^21.3 cm⁻² to pointings showing a steep decrease of f_dark for column densities above ≈ 10²⁰ cm⁻². The difference between the various pointings is most pronounced for large column densities N(H₂|CO) > 10²¹ cm⁻². We note that the larger f_dark values are not a direct consequence of the fact that a relatively dense molecular phase might be required, since in such cases a somewhat lower column density of CO-dark H₂ gas will be determined (Fig. 10). Interestingly, the sequence of pointings in Fig. 18 also correlates with the presence of bright Hα and 24μm emission near molecular clouds (Fig. 1). For instance pointings #9, #11, and #12 show particularly low f_dark values and high CO column densities and they are also the pointings with the faintest Hα (or 24μm) emission. This result is also shown in Fig. 19 where we report the 24μm emission for each peak in the distribution of each pointing.

Stellar feedback could have an impact on the fraction of CO-dark H₂ gas, either through intense radiation field (ionization, radiation pressure) or through dynamical/mechanical effects from winds and supernovae shocks resulting in disruption/dispersal of molecular clouds (e.g., Dale et al. 2012) and therefore to a lower beam-averaged extinction A_V. The average extinction toward a given pointing is driven by the selective disruption and photodissociation of the most diffuse clouds (thereby selecting out the clouds with the highest extinction) and the overall eroding of all clouds. From a theoretical perspective, Wolfire et al. (2010) find that for fixed A_V the fraction of CO-dark H₂ gas in a molecular cloud is insensitive to the incident radiation field, although with the latter not exceeding 30 times the local Galactic radiation field. In contrast, these latter authors find that f_dark increases with decreasing A_V. On the other hand, kiloparsec-scale simulations of spiral galaxies show that f_dark is a function of the radiation field and of the surface density (Smith et al. 2014). Fahrion et al. (2017) modeled the star-forming dwarf galaxy NGC 4214 and found that the fraction of CO-dark H₂ gas mass depends on the evolutionary stage of the star-forming regions, with a larger fraction towards the naked cluster as compared to compact embedded regions. Madden et al., (in prep.) modeled PDRs in galaxies from the Herschel Dwarf Galaxy Survey (DGS) and find that the effective extinction derived in the model is the main parameter controlling f_dark.

Unfortunately, our models do not allow us to directly examine parameters such as A_V or the radiation field strength G₀, for lack of other transitions or dust measurements corresponding to individual velocity components, but we may hope to disentangle both parameters from N(H₂|CO) or [C II]/CO. Since N(H₂|CO) is proportional to the CO intensity which correlates with extinction on parsec scales (Lee et al. 2018), the trends in Fig. 18 can be understood to first order as the variations of f_dark with A_V for a single cloud hypothesis. With this in mind, the decrease of f_dark for N(H₂|CO) ≳ 10²⁰ cm⁻² would be due to the increase of CO column density with the cloud size while the mass in the CO-dark H₂ gas layer saturates, and A_V should then be the main parameter controlling f_dark. However, the fact that we observe different f_dark values for a given N(H₂|CO) value implies either that a parameter other than A_V controls f_dark (a parameter correlating with [C II]/CO since f_dark correlates best with [C II]/CO; Fig. 17) or that N(H₂|CO) does not accurately trace A_V (for instance because our observations combine multiple clouds).

Fig. 17 Bivariate kernel density estimate of the fraction of CO-dark H2 gas fdark vs. the [C II]/CO ratio for all pointings. The shade scales with the density of points. The dashed curve is an attempt toconnect most of the distribution peaks from different pointings (with the plots ordered according the location of probability peaks along the curve) while the solid curve shows the expected relation between fdark and [C II]/CO at low density (~30 cm−3; Fig. 11).

Fig. 18 Bivariate kernel density estimate of the fraction of CO-dark H2 gas vs. the H2 column density measured from CO (i.e., ignoring CO-dark H2 gas) for all pointings. The shade scales with the density of points. The black curves serve as guides to compare the location of the probability peaks in different pointings. Plots are ordered according to the location of the main peak with respect to the black curves.

Fig. 19 Fraction of CO-dark H2 gas plotted against the CO column density for the peaks identified by eye in Fig. 18. The shade scales with the Spitzer/MIPS 24μm surface brightness.

4.5 Influence of metallicity

In this section we apply our models to the data from Pineda et al. (2017) for LMC and SMC pointings (lines of sight toward H I-bright, CO-bright, and/or 160μm-bright regions). The SMC data points are useful probes of a lower-metallicity environment (≈ 0.2Z_⊙) as compared to the LMC (≈ 0.5Z_⊙). The spectral profiles of H I and CO were used in Pineda et al. (2017) but individual velocity components were not adjusted to mitigate the different angular resolution in each tracer. We use the H I CNM value (integrated over the [C II] line FWHM) in Pineda et al. (2017). For our models, the CO-to-H₂ conversion factor for the SMC is taken as cm⁻² (K km s⁻¹)⁻¹ (Roman-Duval et al. 2014). We selected data points with detections in H I, [C II], and CO, and with I([CII]) > 2 × 10⁻⁸ W m⁻² sr⁻¹ in order to be coherent with the present study (Sect. 4.3.1).

Results are shown in Fig. 20. Overall, the results agree with our study for LMC data points, that is, bright [C II] components have a molecular gas fraction near unity, the contribution of neutral atomic gas to the [C II] emission is larger for faint [C II] components (≲ 6 × 10⁻⁸ W m⁻² sr⁻¹), and the fraction of CO-dark H₂ gas decreases with the CO column density. The comparison with our results illustrates some advantages related to the profile decomposition. First, while our study was performed with fewer pointings, by accounting for multiple velocity components corresponding to various physical conditions, statistically significant results could be obtained. Furthermore, we were able to some extent to probe components with large [C II]/CO values which may correspond to relatively quiescent clouds (Fig. 17).

Toward SMC regions, the fraction of [C II] emission corresponding to CO-dark H₂ gas, f_coll,H2([CII]), is somewhat lower, while f_dark is somewhat larger for a given CO column density as compared to the LMC. In other words, more [C II] can be explained by atomic gas in the SMC, but the fraction of CO-dark H₂ gas is larger. It is possible that the larger f_dark values obtained for the SMC are related to a larger incident radiation field but it is also possible that the gas temperature we assume for the SMC(the same as for the LMC) is not applicable. As noted by Glover & Clark (2016), the detectability of the CO-dark H₂ gas with [C II] or [O I] in fact improves greatly with the ISRF strength G₀ (which increases the photoelectric-effect heating rate and hence the gas temperature). While we have chosen the same temperatures for the LMC and SMC regions, f_dark in the SMC would be compatible with the LMC values if the temperature was scaled up by a factor of only two to four. Results for the DGS using integrated measurements show that there is no simple relationship between metallicity and f_dark and that the latter mostly depends on extinction A_V (Madden et al., in prep.).

5 Conclusions

We present a study of the CO-dark H₂ gas in the giant H II region N 11 in the LMC, with a half-solar metallicity. Twelve pointings were observed with SOFIA/GREAT in [C II]. The [C II] velocity profile is compared to that of CO(1–0) observed with MOPRA and ALMA, and to that of H I observed with ATCA+Parkes. The objectives are to measure the total molecular gas content, the fraction of CO-dark H₂ gas, and to probe the influence of the environment (in particular stellar feedback). First we investigated velocity-integrated properties. The main results obtained are as follows.

Fig. 20 Results from Pineda et al. (2017) for LMC (blue) and SMC (red) regions. The symbol size is inversely proportional to the error bar of the y-axis parameter. The curves are the same as in Figures. 11 and 18.

• The [C II] emission originates mostly from the neutral gas, except for one velocity component toward one of the pointings (#5) for which the result is unclear. This result is found through various methods using the [C II] line width, the [N II]/[C II] ratio in the ionized gas, the [N II] 205μm velocity profile observed with GREAT, and the photoelectric-effect heating efficiency proxy ([C II]+[O I])/PAH.

• The photoelectric-effect heating efficiency proxy ([C II]+[O I])/PAH is found to be constant throughout the regions in N 11 observed with Herschel/PACS, generalizing somewhat the results previously obtained in N 11B by Lebouteiller et al. (2012) that the sum [C II]+[O I] traces the total cooling, that PAH emission traces the gas heating, and that [C II] mostly originates from the neutral gas. Our results suggest that this gas is CO-dark H₂ and not atomic.

• The total profile width of [C II] is found to be between that of CO(1–0) and H I, but the [C II] profile resembles more that of CO(1–0).

We then decomposed the profiles using a Bayesian method and a statistical approach that makes use of the many pointings together with a range of input parameters concerning the number of velocity components, the minimum individual component line width, and the minimum separation between components. A simple model was used to compute the [C II] line intensity accounting for collisions with H⁰, H₂, and e⁻ in a two-phase medium (neutral atomic and molecular) as a function of gas temperature and density. The main results obtained are as follows:

• The variations of I([CII]) are driven by the [C II] emission in the molecular phase. As a result, the [C II] components with the largest contribution from H⁰ gas are preferentially those with a low [C II] surface brightness rather than those with low [C II]/CO or low CO column density. However, the contribution from H⁰ gas to [C II] is never dominant.

• There is a sharp transition between CO-bright and CO-dark H₂ gas, withthe latter quickly becoming the dominant H₂ reservoir.

• Overall (combining all pointings and all velocity components), more than 90% of the [C II] emission arises in the (CO-dark) H₂ gas.

• Most of the molecular gas is CO-dark (fraction between 40 and 100%), in particular toward the brightest [C II] components.

• The CO-dark H₂ gas traced by [C II] is rather diffuse with ~200 cm⁻³ on average for all pointings. We identify in particular a specific [C II] velocity component toward #5 with a density around ~100 cm⁻³.

• The contribution of optically thick H I to the dark neutral medium is not significant.

• Most components follow the same trend of the fraction of CO-dark H₂ versus [C II]/CO, with some deviations driven by the gas density.

• The effective X_CO factor including the CO-dark H₂ gas lies in the range 10²¹⁻²² (K km s⁻¹)⁻¹ for most of the bright velocity components.

• The fraction of CO-dark H₂ gas decreases with increasing CO column density, but while it is rather constant for CO column density <10^20.5 cm⁻², it shows a large dispersion above this value. We argue that, for a given CO column density, the larger fractions of CO-dark H₂ gas are found toward CO clouds at the interface with Hα-bright regions. It is plausible that stellar feedback (either through radiation or dynamical/mechanical effects from stellar winds and supernovae shocks) results in the disruption and dispersal of molecular clouds and in turn a lower extinction on average.

• Our simple models were applied to the LMC and SMC pointings in Pineda et al. (2017). We find circumstancial evidence that the fraction of CO-dark H₂ gas is larger for a given CO column density in the SMC as compared to the LMC, but this conclusion is weakened by the uncertain gas temperature.

The main caveat concerns the derived column density and number density in the atomic medium, with limited spatial resolution, even though the velocity decomposition somewhat mitigates the correspondence between components in the different tracers. Further observations at larger spatial resolution in H I would greatly improve our knowledge of the origin of [C II] and its performance as a CO-dark H₂ gas tracer. Moreover, observations of [O I] would shed light on the physical conditions of the few components where [C II] arises in the neutral atomic medium. In particular, more constraints are needed to examine the incident radiation field and extinction in individual clouds, which would then help the understanding of the nature of stellar feedback responsible for the variations of the fraction of CO-dark H₂ gas. Finally, we are aware that the assumption that the velocity components from different tracers arise from a given cloud is increasingly problematic at increasingly large scales, but we cannot unfortunately thoroughly test this hypothesis with the present dataset.

Acknowledgements

Based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00532.S and ADS/JAO.ALMA#2013.1.00556.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatoryis operated by ESO, AUI/NRAO and NAOJ. V. Lebouteiller wishes to thank S. Kannappan and C. Iliadis for making it possible to do part of this work at UNC Chapel Hill. M.C. gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) through an Emmy Noether Research Group, grant number KR4801/1-1. M.-Y.L. was partially funded through the sub-project A6 of the Collaborative Research Council 956, funded by the Deutsche Forschungsgemeinschaft (DFG). D.C. is supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 702622. FLP is supported by the ANR grant LYRICS (ANR-16-CE31-0011). This work was supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES.

Appendix A Profile decomposition for each pointing

The spectral decomposition of the H I CO, and [C II] spectra obtained for the 12 GREAT pointings is shown in Figs. A.1–A.4. The decomposition is only shown for the model with ten components, a minimum line width of 1 km s⁻¹, and a minimum component separation of 1 km s⁻¹.

Fig. A.1 Profile decomposition for pointings #1, #2, and #3.

Fig. A.2 Profile decomposition for pointings #4, #5, and #6.

Fig. A.3 Profile decomposition for pointings #7, #8, and #9.

Fig. A.4 Profile decomposition for pointings #10, #11, and #12.

Appendix B Kernel density estimates

The kernel density estimate of the fraction of [C II] tracing CO-dark H₂ gas, f_coll,H2([CII]), and I([CII]) is shown in Fig. B.1. The distribution is globally similar to that of f(H₂) versus I([CII]) (Fig. 14) because most H₂ is CO-dark and because the CO-dark H₂ gas is traced by [C II].

Figure B.2 shows the kernel density estimate of the effective X_CO conversion factor and the [C II]/CO ratio. Since the X_CO conversion factor includes the contribution of the CO-dark H₂ gas, it is larger than the fiducial value cm⁻² (K km s⁻¹)⁻¹ (Sect. 2.4). We find an effective X_CO in the range 10²¹⁻²² (K km s⁻¹)⁻¹ for most of the bright velocity components, in good agreement with values obtained in Israel (1997), Galliano et al. (2011), Chevance (2016).

Fig. B.1 Bivariate kernel density estimate of the fraction of [C II] tracing CO-dark H2 gas vs. I([CII]) for all pointings. The shade scales with the density of points.

Fig. B.2 Bivariate kernel density estimate of XCO vs. [C II]/CO for all pointings. The shade scales with the density of points.

Appendix C Gas temperature determination

The temperature in the neutral atomic gas can in principle be estimated from the [O I] 63μm/[C II] ratio observed with Herschel/PACS (Lebouteiller et al. 2012). Apart from optical depth effects for the [O I] 63μm line (see Lebouteiller et al. 2012), another problem resides in the fact that [C II], and presumably [O I], mostly trace the molecular phase rather than the neutral atomic phase, as our present results suggest. For this reason we use the extended emission seen in the PACS maps and assume it to be dominated by atomic gas. For extended emission, the [O I] 63μm/[C II] ratio is in the range ~0.2−0.5 (Fig. C.1 left), corresponding to a temperature of a few hundred Kelvin for densities ≲ 10³ cm⁻³ (Fig. C.1 right), while somewhat lower temperature values would be found if the medium were partly molecular.

Fig. C.1 Left: map of [O I] 63μm/[C II] with PACS with values between ≈−0.7 and ≈ 0.4 in log unit. Right: theoretical ratio [O I]/[C II] (in log units) as a function of the gas density and temperature for a purely atomic gas. We use the N 11B chemical abundances from Toribio San Cipriano et al. (2017).

References

All Tables

All Figures

	Fig. 1 LMC-N 11 in CO, Hα, H I, [C II], [O I], and [O III]. The SOFIA/GREAT pointings are shown as black/white circles and labelled in the CO panel. The black rectangles show the Herschel/PACS coverage. Region names are labelled in the [C II] panel. GREAT pointings #1 (N 11A) and #9 (part of N 11D) were not observed in [O III], and pointing #12 (N 11 “South”) was not observed with PACS.
In the text

Fig. 2 Spectra of the main spectral lines used in this study for each GREAT pointing. From top to bottom in each panel: Hα, [Ne III] 3868Å, H I 21 cm, [N II] 205μm, [C II] 158μm, and CO(1–0). The histogram shows the data and the red curve shows the Gaussian fit. For this plot, only one component is considered except for pointing #5. Red dashed vertical lines show the Gaussian fit centroids while black dashed vertical lines show the zero velocity. The spectra are shifted arbitrarily along the y-axis for display purposes.

In the text

	Fig. 3 Ratio of the opacity-corrected (Braun 2012) to observed H0 column densities across the N 11 region. The observed H0 column density for individual velocity components is always <2 × 1021 cm−2 (see Sect. 4.2).
In the text

Fig. 4 Comparison between ALMA, MAGMA, and Herschel data. Left column, from top to bottom: Herschel/PACS [C II] map of N 11B, the ALMA CO(1–0) map (original resolution, convolved to 14.4′′ and projected on PACS [C II] grid, and convolved to 45′′ and projected on MAGMA grid), and the ratio (log units) of CO MAGMA/ALMA. The GREAT pointings are overlaid as 14.4′′ diameter circles. Right column: MAGMA spectrum (blue histogram) and the ALMA spectra calculated for a 45′′ resolution (red with shade) and 17′′ resolution (black) for the four pointings in N 11B.

In the text

	Fig. 5 Theoretical [N II]/[C II] ratios as a function of density. The ratio values include ionic abundance ionization corrections (see text). The red curve indicates ratio values with [C II] originating fully in the ionized gas. The solid black lines show the observed upper limits for all pointings observed with Herschel/PACS (except #5 for [N II] 205μm/[C II] shown as a dashed line).
In the text

Fig. 6 Comparison of the GREAT [N II] and [C II] spectral profiles for the six pointings with [N II] observations. The red solid curve shows for both lines the Gaussian fit to the main velocity components as in Fig. 2. For the [N II] spectrum the black curve corresponds to the rms flux while the red dotted curve shows the expected [N II] profile assuming [C II] arises only from the ionized gas, with the uncertainty shown as the gray area (see text). For the [C II] spectrum the red dotted curve and shaded area show the profile of [C II] in the ionized gas calculated using the [N II] rms (i.e., the brightest possible [C II] component arising from the ionized gas).

In the text

	Fig. 7 Variation of [O I]/PAH (bottom), [C II]/PAH (middle), and ([C II]+[O I])/PAH (top) across the regions in LMC-N 11 mapped with PACS. The size of the symbol scales with the signal-to-noise ratio. Black circles correspond to GREAT pointings. The solid line shows the linear regression for the corresponding ratio while the dashed lines show the regression obtained for the two other ratios.
In the text

	Fig. 8 Illustration of the profile decomposition method for pointing #2 for one of the model input parameter sets (Nc = 10; σmin = 1 km s−1; Δvmin = 1 km s−1). The gray histogram shows the data. The red, green, and dark blue curves correspond to different components and the dotted vertical line shows the velocity for each component. The black solid curve shows the total inferred profile.
In the text

	Fig. 9 Bivariate kernel density estimate (non-parametric probability density function) of I([CII]) vs. the molecular gas fraction ignoring CO-dark H2 gas, f(H2|CO), for all pointings. The shade scales with the density of points.
In the text

Fig. 10 Illustration of the model strategy. Two phases are considered for which the level population of C+ is calculated as a function of density n, temperature T, and the ionization fraction x in the atomic medium. Collisional rates with each partner also depend on temperature. The [C II] intensity is calculated using the level population and the column density of H0 or H2 and the fraction of C into C+ in each phase (see Lebouteiller et al. 2013 for details).

In the text

Fig. 11 Possible model values for fdark as a function of [C II]/CO and volume density in the neutral atomic phase (color bar). The models shown are those corresponding to the observed ranges for N(H0) = [1018, 1021.3] cm−2, N(H2 |CO) = [1018, 1021.5] cm−2, and I([CII]) = [10−8, 10−6] W m−2 sr−1. Model results are shown for a temperature in the neutral atomic phase of 100 and 50 K in the molecular phase. The striping is due to the incomplete coverage in the model results.

In the text

Fig. 12 Illustration of model results for one components in pointing #1 (top) and in #2(bottom). The [C II] surface brightness is calculated as a function of the H2 column density associated with C+. The horizontal red line shows the observed [C II] surface brightness. The purple and pink shaded areas show the range of model results assuming thermal pressure equilibrium and overdense molecular gas respectively (see text). The bottom and top boundaries are set by the H0 column density together with the cloud size constraint. The right boundary is set by the maximum allowed molecular cloud size. In the top panel the model with thermal pressure equilibrium can reproduce the observations while in the bottom panel, a model with an overdense molecular gas is required.

In the text

	Fig. 13 Histogram of the modeled pressure in the neutral atomic (blue) and CO-dark H2 gas (red) forall the velocity components. Results from all decomposition methods are combined (see Table 6).
In the text

	Fig. 14 Same as Fig. 9 but with the total molecular gas fraction (i.e., including the CO-dark H2 gas). The vertical striping is due to the wide ranges for the molecular gas fraction determination for the faintest [C II] components. The shade scales with the density of points.
In the text

	Fig. 15 Bivariate kernel density estimate of the fraction of [C II] tracing CO-dark H2 gas vs. [C II]/CO for all pointings. The shade scales with the density of points.
In the text

Fig. 16 Cumulative fraction of [C II] plotted against the fraction of [C II] emission arising in the CO-dark H2 gas. Values are calculated with all pointings combined. Plots are shown for all components (red), for components with low CO column density (<2 × 1020 cm−2; dark blue), for components with high CO column density (>1 × 1021 cm−2; blue), and for components with weak [C II] (<4 × 10−8 W m−2 sr−1; green).

In the text

Fig. 17 Bivariate kernel density estimate of the fraction of CO-dark H2 gas fdark vs. the [C II]/CO ratio for all pointings. The shade scales with the density of points. The dashed curve is an attempt toconnect most of the distribution peaks from different pointings (with the plots ordered according the location of probability peaks along the curve) while the solid curve shows the expected relation between fdark and [C II]/CO at low density (~30 cm−3; Fig. 11).

In the text

Fig. 18 Bivariate kernel density estimate of the fraction of CO-dark H2 gas vs. the H2 column density measured from CO (i.e., ignoring CO-dark H2 gas) for all pointings. The shade scales with the density of points. The black curves serve as guides to compare the location of the probability peaks in different pointings. Plots are ordered according to the location of the main peak with respect to the black curves.

In the text

	Fig. 19 Fraction of CO-dark H2 gas plotted against the CO column density for the peaks identified by eye in Fig. 18. The shade scales with the Spitzer/MIPS 24μm surface brightness.
In the text

	Fig. 20 Results from Pineda et al. (2017) for LMC (blue) and SMC (red) regions. The symbol size is inversely proportional to the error bar of the y-axis parameter. The curves are the same as in Figures. 11 and 18.
In the text

	Fig. A.1 Profile decomposition for pointings #1, #2, and #3.
In the text

	Fig. A.2 Profile decomposition for pointings #4, #5, and #6.
In the text

	Fig. A.3 Profile decomposition for pointings #7, #8, and #9.
In the text

	Fig. A.4 Profile decomposition for pointings #10, #11, and #12.
In the text

	Fig. B.1 Bivariate kernel density estimate of the fraction of [C II] tracing CO-dark H2 gas vs. I([CII]) for all pointings. The shade scales with the density of points.
In the text

	Fig. B.2 Bivariate kernel density estimate of XCO vs. [C II]/CO for all pointings. The shade scales with the density of points.
In the text

	Fig. C.1 Left: map of [O I] 63μm/[C II] with PACS with values between ≈−0.7 and ≈ 0.4 in log unit. Right: theoretical ratio [O I]/[C II] (in log units) as a function of the gas density and temperature for a purely atomic gas. We use the N 11B chemical abundances from Toribio San Cipriano et al. (2017).
In the text