Physical conditions in the gas phases of the giant HII region LMC-N11: II. Origin of [CII] and fraction of CO-dark gas

(abridged) The ambiguous origin of [CII] 158um in the interstellar medium complicates its use for diagnostics concerning the star-formation rate and physical conditions in photodissociation regions (PDRs). We observed the giant HII region N11 in the Large Magellanic Cloud with SOFIA/GREAT in order to investigate the origin of [CII] to obtain the total H2 gas content, the fraction of CO-dark H2 gas, and the influence of environmental effects such as stellar feedback. We present an innovative spectral decomposition method that allows statistical trends to be derived. The [CII] line is resolved in velocity and compared to HI and CO, using a Bayesian approach to decompose the profiles. A simple model accounting for collisions in the neutral atomic and molecular gas was used in order to derive the H2 column density traced by C+. The profile of [CII] most closely resembles that of CO, but the integrated [CII] line width lies between that of CO and that of HI. Using various methods, we find that [CII] mostly originates from the neutral gas. We show that [CII] mostly traces the CO-dark H2 gas but there is evidence of a weak contribution from neutral atomic gas preferentially in the faintest components. Most of the molecular gas is CO-dark. The fraction of CO-dark H2 gas decreases with increasing CO column density, with a slope that seems to depend on the impinging radiation field from nearby massive stars. Finally we extend previous measurements of the photoelectric-effect heating efficiency, which we find is constant across regions probed with Herschel, with [CII] and [OI] being the main coolants in faint and diffuse, and bright and compact regions, respectively, and with PAH emission tracing the CO-dark H2 gas heating where [CII] and [OI] emit. Our study highlights the importance of velocity-resolved PDR diagnostics and higher spatial resolution for HI observations.


Introduction
The [C II] 158 µm line emits under a variety of conditions in the interstellar medium (ISM) corresponding to the cold and warm The reduced spectra are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc. u-strasbg.fr/viz-bin/cat/J/A+A/632/A106 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 2 P 3/2 fine-structure level (91.3 K). In the neutral gas, C + may exist in the H 0 phase but also in the H 2 phase, in regions where CO is photodissociated while H 2 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 2 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 wellchosen 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 2 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 2 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 2 gas and the H 0 gas, respectively, with the remaining [C II] emission attributed to CO-dark H 2 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 2 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 2 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 2 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 2 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/2 Z ) 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 2 is not traced by CO Chevance 2016). Requena-Torres et al. (2016) observed several star-forming regions in the Small Magellanic Cloud (SMC; ≈1/5 Z ) and found that most of the H 2 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 COdark. Overall, the finding that CO-dark H 2 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 2 gas and the fraction of CO-dark H 2 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 ) 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 photoelectriceffect 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 2 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.

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

SOFIA/GREAT
Twelve pointings were observed within LMC-N 11 (  (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 Notes. (a) Line spread function full width at half maximum (spectral resolution). (b) Point spread function full width at half maximum (spatial resolution). (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 1 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 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 −1 for passband L#2 and 0.20 km s −1 for L#1, but spectra were rebinned to obtain a channel width of 1.2 km s −1 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 (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 −1 ; 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 θ 2 , 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% 3 (see also Fahrion et al. 2017).

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 −1 . The column density is calculated using a conversion factor for optically thin gas 1.823 × 10 18 K km s −1 cm −2 (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 0 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 0 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 selfabsorption) and to the fact that multiple components exist along  the line of sight with different opacities. For the H 0 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).

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 −1 . The column density is calculated using a conversion factor 4 4 We use the notation X to distinguish it from the factor X including the CO-dark H 2 gas contribution. X CO = 2 × 10 20 cm −2 (K km s −1 ) −1 , which is the fiducial value for resolved CO clumps without the contribution of the CO-dark H 2 gas between CO clumps as discussed in Roman-Duval et al. (2014). In Table 4 we provide the CO intensity, the H 2 column density using X CO , N(H 2 |CO), and the molecular gas fraction , ignoring the CO-dark H 2 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 −1 . The line width ranges between ≈2 and ≈5 km s −1 . 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 #4 is 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    (45 )  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.

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 −1 ), and [Ne III] 3868 Å observed with the HR2 grating with R = 19 600 (≈15 km s −1 ). 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.).

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.

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;  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 −1 ) 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 −1 . 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 −1 , while it is ≈30−40 km s −1 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

[C II] line width
The expected thermal broadening 5 for C + in the ionized gas is ≈6 km s −1 (FWHM for a temperature of 10 000 K). Macroturbulent motions (typically a few km s −1 ) may increase this value to ≈10 km s −1 , 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 −1 . The observed FWHM for the main [C II] velocity component is smaller than ≈6 km s −1 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 −1 for the low-and high-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  (Lebouteiller et al. 2012). The comparison with the theoretical ratio in the ionized gas for typical densities between 10 1−3 cm −3 (Fig. 5

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 farinfrared 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 line emission, 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  We assume that the velocity structure reflects the presence of several individual components, which we refer to as "clouds",

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).

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 contributing to 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 −1 . 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 −1 ; 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 −1 is motivated by the typical turbulence measured for interstellar clouds (e.g., Welty et al. 1994Welty et al. , 1996. -The minimum component separation is set to 0 km s −1 (components with same velocity but potentially different line width), 1 and 2 km s −1 (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 A106, page 11 of 30 A&A 632, A106 (2019) Intensity (Gaussian area) Free, uniform >0 Line width Free, half-normal >σ min Separation between components Free, half-uniform >∆v min 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 burnin 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.

Model-independent quantities
For each velocity component we compile the [C II] line intensity, the H 0 column density, the H 2 column density derived from CO, N(H 2 |CO), and the resulting CO-traced molecular gas fraction f (H 2 |CO) = 2N(H 2 |CO)/(2N(H 2 |CO) + N(H 0 )). The CO-dark H 2 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 2 |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 2 |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 2 |CO) 10% or f (H 2 |CO) 60%, suggesting a sharp transition between CO-bright H 2 gas and either CO-dark H 2 or atomic gas. Components with a low f (H 2 |CO) are our best candidates for evidence of significant CO-dark H 2 gas amount.

Selection of components
Our models account for the neutral gas only (atomic and molecular). The ionized gas contribution to the integrated line emission   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 2 gas may be related to a CO clump in which case the cloud emits both in [C II] and CO. The CO-bright H 2 gas is not considered in the model because we only calculate the properties of [C II]-emitting gas, and as such, the H 2 column density associated with CO is simply calculated using a fiducial X CO factor (Sect. 2.4). The strategy goes as follows (illustrated in Fig. 10 (1), the number density n(H 0 ) in the atomic gas is calculated using the H 0 column density and assuming that the individual cloud size along the line of sight lies in the range 1−10 pc 6 . 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 0 ) therefore ranges from a few cm −3 to ∼10 3 cm −3 . Pop. level = f(n H2 , T H2 )

Molecular gas Neutral atomic gas
Pop

3) I([CII]) (at.) vs. I([CII]) (mol.) → f coll,H2 ([CII])
T H2 : fixed 4) n H2 : pressure eq. / overdense → N(H 2 |C + ) I([CII]) (mol.) 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 H 0 or H 2 and the fraction of C into C + in each phase (see Lebouteiller et al. 2013 for details).
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. 2013Pineda et al. , 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 −3 , 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 −4 , 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 2 column density traced by C + , N(H 2 |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 2 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 metalpoor environments due to the larger photoelectric-effect heating efficiency (itself due to the low dust content), to the lack of metal A106, page 14 of 30 V. Lebouteiller et al.: Origin of [C II] and fraction of CO-dark gas in  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 2 gas is therefore twice larger than in the H 0 gas, reaching up to ≈10 3 cm −3 , 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 2 gas may differ significantly.
From these results, we can then calculate the following quantities for each velocity component: -The total H 2 column density, -The fraction of CO-dark H 2 gas, -The mass fraction of molecular gas, We note that the mass of CO-dark H 2 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 2 (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 0 ), N(H 2 |CO), and I([CII]). The correlation between f dark and [C II]/CO is tighter for low densities, corresponding to low H 0 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 7 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 contribution of 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 2 gas lie around 10 1.5−3.5 cm −3 , corresponding to a pressure of 10 3.5−5 K cm −3 (Fig. 13). The total H 2 column density N(H 2 ) measured this way varies across pointings from ≈10 21 cm −2 (in #5) to ≈2 × 10 22 cm −2 (in #10).    13. Histogram of the modeled pressure in the neutral atomic (blue) and CO-dark H 2 gas (red) for all the velocity components. Results from all decomposition methods are combined (see Table 6).

[C II] in the neutral atomic gas
The [C II] components with low f (H 2 |CO) values (Fig. 9) could a priori be either from CO-dark H 2 gas or atomic gas. When including the contribution from the CO-dark H 2 gas calculated by the models, most velocity components reach a large total (including CO-dark H 2 ) molecular gas fraction f (H 2 ) 80% (Fig. 14),  Globally, we find that >95% of the [C II] emission across all pointings can be attributed to CO-dark H 2 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) Fig. 14. Same as Fig. 9 but with the total molecular gas fraction (i.e., including the CO-dark H 2 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 summary, [C II] mostly traces the CO-dark H 2 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 2 gas lies around 10 3−5 K cm −3 (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 2 gas probed by our [C II] observations corresponds to the diffuse ISM (median density ∼200 cm −3 for all pointings) filling the space between CO clumps rather than thin CO-dark H 2 layers around CO clumps.

Fraction of CO-dark H 2 gas
We show in Fig. 17 the distribution of the fraction of CO-dark H 2 gas, f dark , versus [C II]/CO for each pointing. The fraction of CO-dark H 2 gas is also a proxy for the fraction of DNM (i.e., CO-dark H 2 gas and optically-thick H I) since the CO-dark H 2 gas traced by [C II] is the main contribution ( 95%) to the DNM according 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 2 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-togas 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 2 mass is not traced by CO. Other studies in the Magellanic Clouds also show significant CO-dark H 2 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 2 gas. Chevance (2016) observed LMC-30 Dor and determined from PDR models that 80% of the molecular gas is CO-dark H 2 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 2 gas in the model is by construction traced by C + and because the 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 no correction 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 0 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 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 2 |CO) 10 20 cm −2 ) as compared to CO peaks (∼70−90% for N(H 2 |CO) ∼ 10 21 cm −2 ). 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 2 gas lies in the range 10 21−22 (K km s −1 ) −1 for most of the bright velocity components (Fig. B.2), in good agreement with values obtained in Israel (1997), Galliano et al. (2011), andChevance (2016).

Physical parameters controlling f dark
We find that f dark is somewhat larger 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 2 |CO) distribution, from pointings with a high and almost flat f dark distribution until column densities of ≈10 21.3 cm −2 to pointings showing a steep decrease of f dark for column densities above ≈10 20 cm −2 . The difference between the various pointings is most pronounced for large column densities N(H 2 |CO) > 10 21 cm −2 . 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 2 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  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 2 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 2 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 2 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 0 , for lack of other transitions or dust measurements corresponding to individual velocity components, but we may hope to disentangle both parameters from N(H 2 |CO) or [C II]/CO. Since N(H 2 |CO) is proportional to the CO intensity which correlates with extinction on parsec scales , 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 2 |CO) 10 20 cm −2 would be due to the increase of CO column density with the cloud size while the mass in the CO-dark H 2 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 2 |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 2 |CO) does not accurately trace A V (for instance because our observations combine multiple clouds).

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.2 Z ) as compared to the LMC (≈0.5 Z ). 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 2 conversion factor for the SMC is taken as X CO = 1 × 10 21 cm −2 (K km s −1 ) −1 (Roman-Duval et al. 2014). We selected data points with detections in H I, [C II], and CO, and with I([CII]) > 2 × 10 −8 W m −2 sr −1 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 −8 W m −2 sr −1 ), and the fraction of CO-dark H 2 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 2 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 2 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 2 gas with [C II] or [O I] in fact improves greatly with the ISRF strength G 0 (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.).

Conclusions
We present a study of the CO-dark H 2 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 2 gas, and to probe the influence of the 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 2 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 0 , H 2 , 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 0 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 0 gas to [C II] is never dominant. -There is a sharp transition between CO-bright and CO-dark H 2 gas, with the latter quickly becoming the dominant H 2 reservoir. -Overall (combining all pointings and all velocity components), more than 90% of the [C II] emission arises in the (CO-dark) H 2 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 2 gas traced by [C II] is rather diffuse with ∼200 cm −3 on average for all pointings. We identify in particular a specific [C II] velocity component toward #5 with a density around ∼100 cm −3 . -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 2 versus [C II]/CO, with some deviations driven by the gas density. -The effective X CO factor including the CO-dark H 2 gas lies in the range 10 21−22 (K km s −1 ) −1 for most of the bright velocity components. -The fraction of CO-dark H 2 gas decreases with increasing CO column density, but while it is rather constant for CO column density <10 20.5 cm −2 , it shows a large dispersion above this value. We argue that, for a given CO column density, the larger fractions of CO-dark H 2 gas are found toward CO clouds at the interface with Hα-bright regions. It is plausible that stellar feedback (either through radiation or A106, page 23 of 30 A&A 632, A106 (2019) 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 2 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 2 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 2 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.