Mélanie Chevance
email: chevance@uni-heidelberg.de

Suzanne C Madden

Christian Fischer

William D Vacca

Vianney Lebouteiller

Dario Fadda

Frédéric Galliano

Remy Indebetouw

J M Diederik Kruijssen

Min-Young Lee

Albrecht Poglitsch

Fiorella L Polles

Diane Cormier

Sacha Hony

Christof Iserlohe

Alfred Krabbe

Margaret Meixner

Elena Sabbi

Hans Zinnecker

The CO-dark molecular gas mass in 30 Doradus

Keywords: ISM, general -photon-dominated region (PDR) -ISM, clouds -ISM, structure -Magellanic Clouds -ISM, individual objects, LMC-30Doradus

Determining the efficiency with which gas is converted into stars in galaxies requires an accurate determination of the total reservoir of molecular gas mass. However, despite being the most abundant molecule in the Universe, H 2 is challenging to detect through direct observations and indirect methods have to be used to estimate the total molecular gas reservoir. These are often based on scaling relations from tracers such as CO or dust, and are generally calibrated in the Milky Way. Yet, evidence that these scaling relations are environmentally dependent is growing. In particular, the commonly used CO-to-H 2 conversion factor (X CO ) is expected to be higher in metal-poor and/or strongly UV-irradiated environments. We use new SOFIA/FIFI-LS observations of far-infrared fine structure lines from the ionised and neutral gas and the Meudon photodissociation region model to constrain the physical properties and the structure of the gas in the massive star-forming region of 30 Doradus in the Large Magellanic Cloud, and determine the spatially resolved distribution of the total reservoir of molecular gas in the proximity of the young massive cluster R136. We compare this value with the molecular gas mass inferred from ground-based CO observations and dust-based estimates to quantify the impact of this extreme environment on commonly used tracers of the molecular gas. We find that the strong radiation field combined with the half-solar metallicity of the surrounding gas are responsible for a large reservoir of "CO-dark" molecular gas, leaving a large fraction of the total H 2 gas ( 75%) undetected when adopting a standard X CO factor in this massive star-forming region.

INTRODUCTION

The strong correlation between tracers of star formation and tracers of molecular gas observed on kpc scales in galaxies (e.g. Kennicutt 1998;Wong & Blitz 2002;Leroy et al. 2008;Bigiel et al. 2011) suggests that star formation is mostly fueled by molecular gas (e.g. Krumholz & McKee 2005;Elmegreen 2007;McKee & Ostriker 2007;Krumholz et al. 2009;Ostriker et al. 2010). While Glover & Clark (2012a) show that this correlation does not necessarily imply a direct causation, the density and shielding conditions required for star formation are similar to those for the formation of molecular hydrogen, H 2 (e.g. Krumholz & Gnedin 2011;Glover & Clark 2012a;Glover & Clark 2012b). This is the most abundant molecule in the interstellar medium (ISM), and is generally a major constituent of the star-forming gas. However, H 2 is a symmetric molecule, with no electric dipole moment. The lowest rotationvibrational transition is thus the quadrupole transition at ∼28 µm (hν/k ∼ 500 K), which is not excited in the cold gas phase associated with star formation (∼ 10 K). Other molecules, such as carbon monoxide (CO), are easily excited even at low temperatures, and are used as a proxy for the total reservoir of molecular gas, both in local galaxies (e.g. Bigiel et al. 2011;Leroy et al. 2011;Schruba et al. 2012;Saintonge et al. 2017), and at high redshift (e.g. Daddi et al. 2010;Walter et al. 2014;Genzel et al. 2015;Pavesi et al. 2018;Tacconi et al. 2018). The emission of the CO J = 1 -0 transition at 2.6 mm, hereafter CO(1-0), is frequently used as a tracer of the molecular gas mass. A simple CO-to-H 2 conversion factor X CO is used to infer the H 2 column density, N(H 2 ), from the CO(1-0) line intensity, I CO(1-0) , such that N(H 2 ) = X CO × I CO(1-0) . In the Milky Way, the common adopted value for this conversion factor is X CO,MW = 2 × 10 20 cm -2 (K km s -1 ) -1 (see e.g Bolatto et al. 2013, and references therein).

While CO seems to trace the bulk of H 2 in solar metallicity galaxies relatively well (e.g. Bolatto et al. 2013), detection of CO in dwarf galaxies is challenging, even for metallicities only a few times below solar (e.g. Madden et al. 2013;Cormier et al. 2014). In particular, CO emission has been shown to be disproportionately weak in the Large Magellanic Cloud (LMC) and other low-metallicity galaxies (Wong et al. 2011;Hughes et al. 2010) compared to their star formation activities (estimated via [C II] 158 µm emission or far-infrared luminosity; e.g. Madden et al. 2013;Cormier et al. 2014;De Looze et al. 2014;Cormier et al. 2015). By contrast, the star formation activity of some lowmetallicity galaxies as probed by the ratio of the [CII] luminosity over the CO luminosity (e.g. Stacey et al. 1991) can be very high (Madden 2000), as is particularly the case for blue compact dwarf galaxies. As a result, star-forming low-metallicity galaxies are often outliers on the Schmidt-Kennicutt relation (Kennicutt 1989), displaying lower gas surface densities at a given star formation rate density than solar metallicity galaxies (e.g. Taylor et al. 1998;Galametz et al. 2009;Schruba et al. 2012;Cormier et al. 2014). Whether this results from an decreased gas depletion time (and in this case, H 2 is indeed less abundant than expected from solar metallicity galaxy relations) or from the fact that CO fails at being an effective tracer of the total molecular gas in these environments (with the H 2 reservoir and the gas depletion time actually being comparable to normal galaxies) remains unanswered.

Indeed, the X CO factor is expected to vary with environmental properties, and notorious uncertainties exist when inferring the H 2 reservoir from CO, especially in low-metallicity environments. While a constant X CO factor, or a linear scaling of this value with metallicity, is often assumed to derive the total mass of molecular gas, these values remain rough estimates (e.g. Schruba et al. 2012, see also the review by Bolatto et al. 2013). These uncertainties strongly limit our interpretation of the evolution of star formation activity and gas mass at high redshift (e.g. Daddi et al. 2010;Tacconi et al. 2010;Genzel et al. 2012). Indeed, contrarily to H 2 , which is efficiently self-shielded, CO is photodissociated by the ra-diation field at low dust column densities. Therefore, X CO has been suggested to increase in metal-poor and/or highly irradiated environments, where the transition between C + /C 0 /CO is shifted deeper into the cloud, whereas a substantial reservoir of H 2 , called "COdark" molecular gas, can be found outside the CO-emitting regions (where carbon is found in atomic or ionised form; e.g. Israel 1997;Wolfire et al. 2010;Bolatto et al. 2013, Madden et al. in prep.).

Other indirect probes of the molecular gas have been used to trace this potentially important reservoir of H 2 . These include γ-ray emission (e.g Grenier et al. 2005), dust emission with an assumed dust-to-gas mass ratio (e.g Planck Collaboration et al. 2011), or emission from the [C I] line (e.g [START_REF] Papadopoulos | [END_REF]Bell et al. 2006Bell et al. , 2007;;Valentino et al. 2018). However, the calibration of these diagnostics as a function of the environment also remains uncertain. For example, the dust-to-gas mass ratio has also been shown to depend strongly on the metallicity (Roman-Duval et al. 2014;Rémy-Ruyer et al. 2014, 2015) and shows significant scatter between galaxies, most likely due to different star formation histories (see Galliano et al. 2018, for a review).

The question remains as to whether there is a true deficit of H 2 in low-metallicity environments relative to their star formation activity, or a lack an accurate measurement of the gas reservoir under lower dust shielding conditions. Low-metallicity galaxies are often characterised by higher and harder radiation fields (e.g. Madden et al. 2006;Cormier et al. 2015;Cormier et al. 2019), and a reduced dust content (e.g. Rémy-Ruyer et al. 2014), which both contribute to changing the structure of the ISM, making it in particular more porous to photons than in higher metallicity galaxies (Cormier et al. 2015;Chevance et al. 2016). These properties affect the star formation process by potentially reducing the reservoir of gas available for star formation, and influence how the stellar radiation interacts with the surrounding medium. These characteristics may resemble the ones of the chemically young medium (∼ 0.4 Z , e.g. Yuan et al. 2013;Ly et al. 2016) found in galaxies at redshift ∼ 2-3. This makes nearby low-metallicity star-forming dwarf galaxies perfect laboratories to investigate these questions at high spatial resolution. This is a critical step in order to obtain a better understanding of the conditions under which most of the stars in the Universe formed (at z ∼ 2 -3, Madau & Dickinson 2014). Eventually, nearby dwarf galaxies can be used to calibrate new diagnostics for the total reservoir of molecular gas as a function of the environment (e.g. [START_REF] Lupi | [END_REF].

In this paper we focus on the extreme environment around the Super Stellar Cluster (SSC) R136 on the 30 Doradus region (hereafter "30Dor") in the LMC. At a distance of 50 kpc [START_REF] Walker | [END_REF]), the LMC is our closest neighbouring galaxy, and has a moderately low metallicity of 0.5 Z (Rolleston et al. 2002;[START_REF] Pagel | CNO Abundances in Dwarf and Spiral Galaxies[END_REF]. 30Dor is the most prominent and most massive star-forming region in the LMC. The SSC R136 in its centre contains a large population of massive stars (see the recent review by Crowther 2019, and references therein), irradiating the surrounding gas. In a previous paper (Chevance et al. 2016), we have combined atomic finestructure line observations of 30Dor from Herschel and Spitzer with ground-based CO data to provide diagnostics on the properties and the structure of the gas, by modeling it with the "Meudon" model of photodissociation regions (PDR) 1 (Le Petit et al. 2006;Bron et al. 2014). In the study by Lee et al. (2019), we expanded the neutral gas tracers to include additional CO lines accessible from the Herschel/SPIRE FTS observations of the 30Dor region.

The observed (or modeled) reservoir of molecular gas in a single region is dependent on the current evolutionary stage of the region (e.g. Kruijssen et al. 2018). A small molecular gas mass can indeed result from H 2 having been consumed by the on-going star formation event, or having been photodissociated by the stellar feedback from young massive stars, while a large H 2 reservoir potentially indicates that the star formation event just started in this region. Measuring the total molecular gas mass in a single region therefore does not allow us to determine whether the star formation is more efficient at low metallicity. However, this is a first step towards a better understanding of how stellar feedback interacts with the surrounding sub-solar metallicity gas.

In the current paper, we determine additional constraints on the physical conditions of the ionised gas with a larger spatial coverage achieved by recent observations with the Stratospheric Observatory for Infrared Astronomy (SOFIA; Young et al. 2012;Krabbe et al. 2013) of FIR fine structure lines from the ionised and neutral gas. We make use of these new observations to predict the distribution of the total H 2 column density and mass surface density in this region. We present the data in Section 2 and discuss empirical diagnostics of the ISM conditions in Section 3. In Section 4, we use the Meudon PDR model to calculate the mass surface density and total mass of H 2 in 30Dor and determine the X CO factor for this region, by comparison to the standard Galactic value X CO,MW . We discuss environmental variations of this measurement and compare it with literature values in Section 5. Key results and conclusions are summarised in Section 6.

OBSERVATIONS

We present here the observational data used in this study to constrain the physical conditions of the gas in 30Dor and predict the total reservoir of molecular gas. We use new observations of FIR fine structure lines from the ionised and neutral gas with SOFIA, in addition to archival spectrometric and photometric observations from Herschel and molecular gas observations from Mopra.

SOFIA/FIFI-LS

We present new observations of the [O III] lines at 52 µm and 88 µm, the [O I] line at 145 µm and the [C II] line at 158 µm, taken with the Far Infrared Field-Imaging Line Spectrometer (FIFI-LS) instrument (Krabbe et al. 2013;[START_REF] Fischer | [END_REF]Colditz et al. 2018) aboard SOFIA. FIFI-LS is an imaging spectrometer with a wavelength dependent spectral resolution R between ∼ 500 and ∼ 2000. The instrument provides two channels for simultaneous observations between 51-120 µm (blue channel) and 115-203 µm (red channel). Each channel has a 5 × 5 pixel projection onto the sky, with a size of 6 ×6 (∼1.5 pc × 1.5 pc) in the blue channel and 12 ×12 (∼2.9 pc × 2.9 pc) in the red channel, yielding field of views of respectively 30 ×30 (blue channel) and 60 ×60 (red channel). Each pixel is a so-called "spaxel", which means that internally the light is dispersed spectrally (using a grating) over 16 pixels for each spaxel, providing an integral-field data cube for each observation, covering a total spectral bandwidth between 1500 and 3000 km s -1 . The detectors consist of stressed, gallium-doped germanium photoconductors.

For both [O III] lines and the [C II] line the data was acquired in asymmetric chop mode with a fast mapping "AABAA" scheme Data reduction and flux calibration. Data reduction was carried out using the FIFI-LS data reduction pipeline [START_REF] Vacca | Proceedings of the ADASS XXIX Valentino F[END_REF]. The chop and nod cycles are first subtracted. Then, each spectrum is resampled to a regular wavelength grid with approximately eight pixels per FWHM of the instrument's spectral profile function at the observing wavelength. The spectra from the 25 spaxels are re-arranged back to their original 5 × 5 positions on the sky which differ slightly from the ideal 12 rectangular raster.

A telluric correction is then applied and all data cubes are combined to get full maps with regular spatial sampling. At constant wavelength, spatial positions are interpolated locally with a polynomial surface fit to a rectilinear, regular spatial grid in right ascension and declination. The spatial sampling is 1 for the blue channel and 2 for the red channel. Flux calibration is established by observing calibration sources such as the planet Mars and measured spectra are compared to theoretical spectra, scaled to the distance between Mars and Earth. A full calibration across the wavelength range is done every flight series. The absolute flux calibration uncertainty is assumed to be 20%, of which 10% is relative uncertainty from FIFI-LS seen between flight series and different calibrators. The remaining uncertainties arise from atmospheric corrections and calibrator modelling uncertainties.

Atmospheric correction. Absorption of far-infrared light in the atmosphere due to water vapour or trace gases like ozone can have a significant impact on determining the true astronomical signal even at typical observing altitudes between 39,000 and 45,000 feet. Thus, telluric correction is critical to measure source flux. The main influencing factors for atmospheric transmission are flight altitude, source elevation and precipitable water vapour (PWV). The PWV is variable and could not be measured on board SOFIA at the time of observation. Predictions of the variability are available from weather forecasts but have not proven to be reliable. The FIFI-LS pipeline uses a modelled transmission spectrum appropriate for the flight altitude and source elevation, with a standard value for water vapor overburden, as calculated with the ATRAN model [START_REF] Lord | A new software tool for computing Earth's atmospheric transmission of near-and far-infrared radiation[END_REF]. Based on a worst-case assumption of a PWV uncertainty of 50%, the transmission uncertainty here is smaller than 5% for all observed lines. After the spectral resampling of each nod subtracted data cube, the spectrum is divided by the convolved transmission spectrum for that cube (total integration time per cube is about 75s). This allows a correction for changing elevation and altitude on the time scale of minutes. This technique can be problematic near narrow atmospheric features since the convolution with the instrument's resolution significantly changes the shape of those features whereas the astronomical signal interacts with the unconvolved transmission profile of the atmosphere. However, there are no narrow lines at the specific observed line locations in 30Dor.

Flux measurement. We use FLUXER 2 to measure the flux from the combined and calibrated cubes. For each pixel, the continuum is modeled with a zeroth order polynomial and the emission line is fitted using a gaussian function. The resulting emission maps for [O III] at 52 and 88 µm, [C II] at 158 µm and [O I] at 145 µm are presented in Figure 1.

Combination of PACS and FIFI-LS observations for [O I] 145 µm. As the [O I] 145 µm fine structure line requires the longest integration time with FIFI-LS, we have maximised the total coverage by not duplicating the PACS observations. We have ensured partial overlap between the maps for cross-calibration purposes. As a result, a small region in the central part of 30Dor was not observed by FIFI-LS in [O I] 145 µm, as visible in Figure 1. To recover the full map, we combine the FIFI-LS and PACS maps as follows. We first check the agreement between the relative flux calibrations of the FIFI-LS and PACS observations for the areas that were observed by both instruments (equivalent to approximately 20 FIFI-LS beams). They agree at the ∼ 10% level on 12 beam size once both maps are convolved to the same spatial resolution. In view of the ∼20% absolute calibration of PACS and FIFI-LS we can then establish that the relative flux calibrations between both sets of observations are in agreement and we therefore simply match the FIFI-LS and PACS maps on the same 2 × 2 pixel grid. We assign the observed FIFI-LS (respectively PACS) flux to pixels that only have flux observation above a signal-to-noise ratio (SNR) of 3 with FIFI-LS (respectively PACS), and the SNR-weighted mean flux of the FIFI-LS and PACS fluxes for pixels that are covered by both instruments.

Archival data

The Herschel/PACS (Poglitsch et al. 2010) spectroscopic data of 30Dor have been presented in Chevance et al. (2016), and we refer 2 http://www.ciserlohe.de/fluxer/fluxer.html to Cormier et al. (2015) for the full description of the PACS observations and data reduction. We use here only the maps of [O I] at 145 µm to complement the field of view covered by FIFI-LS, as described in Section 2.1.

One of the parameter of the PDR model is the total extinction in the V band, A max V , which scales with the depth of the PDR layer. Ratios of lines originating from the same depth in a PDR cloud are not sensitive to the size of the cloud (e.g. [C II] and [O I]; see also Section 4.1). In addition to the PDR tracers observed by FIFI-LS, tracers originating from different depths into the cloud are therefore necessary to constrain the total depth of the cloud A max V . We use observations of the CO(1-0) from Mopra (Wong et al. 2011) at a spatial resolution of 42 , which we also use to estimate the mass surface density of the CO-bright molecular gas (Section 4.2). For the atomic gas, we use H I observations of the LMC from Kim et al. (2003) at a spatial resolution of 1 .

In order to calculate the FIR luminosity, we use the Herschel/PACS and SPIRE maps of the Large Magellanic Cloud at 100, 160 and 250 µm, first published in Meixner et al. (2013) as part of the HERITAGE project. We also make use of the observations of 30Dor obtained as part of the Spitzer (Werner et al. 2004) Legacy program "Surveying the Agents of a Galaxy's Evolution" (SAGE; Meixner et al. 2006). We used the four channels of IRAC (Fazio et al. 2004) at 3.6, 4.5, 5.8 and 8.0 µm and MIPS (Rieke et al. 2004) observations at 24 and 70 µm. As the MIPS 24 µm map is saturated in several pixels, we use the IRS spectra (Indebetouw et al. 2009; see also Sect. 5.2.1) to calculate the 24 µm synthetic photometry in the MIPS 24 bandpass and compare to the original map. These are in excellent agreement in parts where the Spitzer/MIPS map is not saturated. We apply the dust SED model of Galliano et al. (2011, AC composition) to these data, as described in Chevance et al. (2016). Since our analysis here is limited by the spatial resolution of the CO(1-0) observations (i.e. 42 ), we use the photometric bands of SPIRE up to 250 µm to constrain the SED model in order to maximise the accuracy on the dust mass (this contrasts with Chevance et al. (2016), where only the FIR luminosity was used, and obtained at a good accuracy with the photometric bands up to 160 µm, which lead to a resolution of 12 ).

EMPIRICAL DIAGNOSTICS

[C II]/[O III] ratio

In metal-rich environments, [C II]158 µm and [O I]63 µm lines are the dominant FIR emission lines. However, this is often not the case in dwarf galaxies, at lower metallicity (e.g. Hunter et al. 2001;Cormier et al. 2015), where [O III] lines become more prominent.

In 30Dor, we now extend the analysis of Chevance et al. (2016), covering projected distances from the centre of the excitation source R136 only up to ∼40 pc, to the entire FIFI-LS map. We show that the emission from the [O III] lines (both at 52 µm and 88 µm) is relatively extended, such that the [O III]88 µm line is brighter than [C II] in the entire covered area, up to a projected distance of 60 pc from the centre of R136 (see Figure 2). energy photons (the ionisation potential of O + is 35.1 eV) can travel over long distances before interacting with the ISM, revealing the high level of porosity of low-metallicity environments. As a result of this porosity of the ISM, R136 is indeed expected to be the dominant source of ionisation in 30Dor (Kawada et al. 2011;Lopez et al. 2011).

Density of the ionised gas

The ratio of the [O III] fine-structure cooling lines is a good probe of the electron density in the highly ionised gas. The [O III] 52 µm line was unobservable with Herschel in the LMC so this diagnostic could not easily be evaluated for star forming regions of the LMC before SOFIA. The ratio of the [S III] lines at 18 µm and 33 µm, commonly observed by Spitzer, is sensitive to high gas densities (the critical densities of these lines are respectively 1.5 × 10 4 cm -3 and 4.1 × 10 3 cm -3 ), which provides only limited constraints in the (low) electron density regime of 30Dor (Indebetouw et al. 2009 estimate an electron density of 300-600 cm -3 using this ratio; see also Chevance et al. 2016). The critical density of the [O III]-emitting material (n crit,OIII88 = 5.1 × 10 2 cm -3 ; n crit,OIII52 = 3.6 × 10 3 cm -3 ) is lower than those of classical dense H II regions, allowing us to probe the low-density ionised gas.

The ratio map of [O III]52 µm/[O III]88 µm, along with the theoretical variations of this ratio as a function of the electron density and temperature are shown in Figure 3. The theoretical ratio depends mostly on the density and only slightly on the temperature. We choose a typical temperature of 10 000 K for the ionised gas around a star-forming region to determine the density constrained by our observations. In our FIFI-LS observations of 30Dor, the ratio [O III]52 µm/[O III]88 µm ranges between 0.48 and 1.84 at a resolution of ∼ 10 , leading to a maximum electron density of a few 10 2 cm -3 (up to ∼ 500 cm -3 ). Where the ratio is lower than ∼ 0.65 (typically at large distances from the centre of R136), we can only infer an upper limit of ∼ 30 cm -3 for the electron density.

PDR MODELING

In the proximity of a star forming region in low-metallicity environments, the strong FUV radiation from the nearby massive stars combined with the reduced metal and dust abundances results in a shift of the C + /C 0 /CO transition deeper into the gas clouds. This effect induces a smaller CO core, surrounded by a larger shell of C 0 and C + which can potentially encompass an important reservoir of H 2 (e.g. Wolfire et al. 2010;Bolatto et al. 2013). The H 2 molecule is self-shielded at a smaller column density relative to that for the CO molecule, which makes the photodissociation rate of CO larger than that of H 2 at the cloud surface.

In this section, we show how detailed PDR modeling can be used to determine the physical conditions of the neutral gas and we compute the depths at which H 2 and CO form. This model, constrained by observations, allows us to calculate the total mass of molecular gas.

Physical conditions of the PDR

We use the Meudon PDR model (Le Petit et al. 2006;Le Bourlot et al. 2012;Bron et al. 2014) to calculate the mass surface density of molecular gas in 30Dor. This model computes the atomic and molecular structure of interstellar clouds in a 1D stationary plane parallel slab of gas and dust illuminated by a radiation field (from UV to radio) arising from one or both sides. We use a constant pressure equation of state throughout the cloud.

We use the method described in Chevance et al. (2016) to find the PDR model parameters that provide the best fit to the observed line ratios

([O I] 145 + [C II] 158 )/L FIR , [O I] 145 /[C II] 158 and [C II] 158 /CO(1-0)
, at a spatial resolution of 43 (determined by the CO(1-0) observations). The parameters we fit for are the incident radiation field on one side of the cloud, G UV , the pressure in the gas, P, and the total visual extinction of the model, A max V . We fit for these three parameters within each 30 pixel. The resulting spatial distributions of these three parameters are presented in Figure 4. The incident radiation field on the PDR, G UV , ranges between 7 × 10 2 and 7 × 10 3 (in units of the standard radiation field observed in the solar neighborhood, where the integrated energy density between 911.8 Å to 2400 Å is G 0 = 6.8 × 10 -14 erg cm -3 ; Mathis et al. 1983). We note that we observe a minimum of G UV in the direction of R136. This indicates that the surface of the PDR is likely to be far away from the cluster in that region, and therefore the received flux by the PDR is low (see Section 5.3). The measured pressure, P, ranges between 4.3 × 10 5 and 1.6 × 10 6 K cm -3 . A max V ranges between 2 and 3 mag. By comparing the predictions of the constrained model for each pixel to the [C II] observations, we measure the area filling factor, Φ A , similarly to Chevance et al. (2016). We find Φ A ∼ 1 -6, which means that typically several clouds are present along the line of sight.

Given the lower resolution at which this analysis is performed (set by the CO data), the results for G UV , P, and Φ A measured here are consistent with the results of Chevance et al. (2016). Our measurements of G UV and P are also in close agreement with the values constrained at a similar spatial resolution by Lee et al. (2019), using the Meudon PDR model, and using the [C I] line at 370 µm instead of the CO(1-0) line. If we exclude the constraint provided by the ratio

[C II] 158
CO(1-0) , it is possible to constrain the PDR parameters G UV and P at a higher resolution of ∼ 18 , limited by the FIFI-LS [C II] observations, instead of 43 as presented here. In this case, the values and spatial distributions of G UV and P using the FIFI-LS observations are in close agreement with the results of Chevance et al. (2016) using the PACS observations, but this does not allow us to constrain A max V , necessary for the determination of the total molecular gas mass. The reason is that [C II] and CO come from different depths into the PDR cloud, which makes the ratio [C II]/CO sensitive to A max V , contrary to ratios involving only [C II], [O I], and L FIR (originating from similar parts of the PDR cloud). We note that A max V constrained here using the ratio [C II]/CO(1-0) is slightly larger than A max V for an ensemble of clouds is therefore likely to be smaller for [C I] emitting clouds than for CO emitting clouds. We note however that this small uncertainty on A max V only slightly affects the predicted H 2 mass and the results presented in Section 4.2, since the transition between H and H 2 occurs closer to the surface of the cloud (A V 1).

Recent observations of several regions of the LMC, including

Total H 2 reservoir

From the Meudon PDR model, using the constrained parameters P, G UV , A max V , and Φ A from Section 4.1 as inputs, we determine the column density and mass surface density of molecular gas predicted by the model for each 30 ×30 pixel. We estimate the H 2 column density between 1.9×10 21 and 1.3×10 22 cm -2 and the H 2 mass surface density to be approximately between 30 and 220 M pc -2 . The spatial distribution of the H 2 mass surface density is presented in the left-hand panel of Figure 5. The position of the maximum corresponds to the location of the [C II] peak. It is also co-spatial with the CO (2-1) detected with ALMA by Indebetouw et al. (2013). When integrating over the observed field of view, we find a total mass of H 2 of 1.8 × 10 5 M .

We compare the mass surface density of H 2 estimated using the PDR modeling with the mass surface density of H 2 that we can infer from the CO (1-0) emission at 43 resolution, using a standard galactic α CO,MW conversion factor of 4.3 M pc -2 (K km s -1 ) -1 (corresponding to X CO,MW = 2 × 10 20 cm -2 (K km s -1 ) -1 ), which can be approximately applied down to metallicities of ∼ 0.5 Z according to Bolatto et al. (2013). Given the large uncertainties and the large scatter of the α CO measurements, it is indeed not clear whether α CO always scales with metallicity or if there is a plateau above ∼ 0.5 Z . The difference between the modeled H 2 mass surface density from the Meudon PDR model and the CO-bright H 2 mass surface density calculated from the CO (1-0) observations, normalised by the modeled total H 2 mass surface density, is presented in the middle panel of Figure 5. This quantity can be interpreted as the fraction of molecular gas not probed by CO (1-0), i.e. as the fraction of "CO-dark" molecular gas, and its spatial distribution will be discussed in more details in Section 5.3. As shown in the middle panel of Figure 5 we find that 75-97% of the total H 2 mass surface density is associated with CO-dark gas within several tens of parsecs around the SSC R136. A large fraction of the molecular gas mass surface density is hence not probed when using CO (1-0) observations and a standard Milky Way X CO factor in this half-metallicity, highly irradiated environment.

We can then calculate the equivalent X CO conversion factor adapted for this region using:

X CO,30Dor = N(H 2 ) model I CO(1-0) (1)
where N(H 2 ) model is the H 2 column density from the PDR model in cm -2 and I CO(1-0) is the observed CO intensity in K km s -1 . The X CO,30Dor conversion factor thus calculated ranges between 8×10 20 and 4.0×10 22 cm -2 (K km s -1 ) -1 (between ∼ 4 and ∼ 20 times X CO,MW ; right-hand panel of Figure 5). We compare this range of values with previous literature measurements in Section 5.2.

DISCUSSION

Effect of the spatial resolution

At the distance of the Andromeda galaxy (∼ 780 kpc ; Ribas et al. 2005) the entire region of ∼ 44 pc × 65 pc mapped here (limited by the overlap between the [O I] and CO(1-0) field of views) would fall approximately in only one spaxel of the FIFI-LS spectrometer.

To test whether the same analysis presented in this paper based on unresolved observations would lead to similar conclusions, we integrate all quantities over the covered area of 30Dor and repeat our analysis. Indeed, we have shown in Chevance et al. (2016) that the physical conditions derived from integrating over a region are typically biased toward the conditions of the high-pressure, highlyirradiated sub-regions, which might also affect the measured total H 2 reservoir in an unresolved case, where different environments are mixed within a beam. We do not see a significant effect of the resolution on the measured X CO factor or fraction of CO-dark gas (between ∼ 10 and 50 pc scale resolution). The integrated IR tracers over the full map are best reproduced by a PDR model with G UV = 1390, P = 6.18×10 5 K cm -3 and A max V ∼ 3 mag (Chevance et al. 2016). The H 2 mass surface density predicted by a model with these particular input conditions is 194 M pc -2 . It would therefore represent a total H 2 mass of 3.4 × 10 5 M over the same region as mapped in Figure 5. This is less than a factor of two higher than the total value of 1.8 × 10 5 M for the sum of the H 2 mass estimated with the PDR model at 43 from Figure 5. In addition, while the mean X CO,30Dor

factor is ∼ 1.9 × 10 21 cm -2 (K km s -1 ) -1 , the X CO,30Dor from the integrated model would be ∼ 2.3 × 10 21 cm -2 (K km s -1 ) -1 . This weak dependence of our measurement with the resolution could be due to the fact that we are still probing a relatively limited range of environments, under the dominant influence of R136.

Comparison with previous X CO estimations

In this section, we present some of the different methods commonly used to estimate the total mass of molecular gas. We compare our determination of the X CO factor in 30Dor from Section 4.2 to values obtained from these different methods in previous studies. The results are summarised in Table 2.

Direct measurement of H 2 in the warm gas

Observations of the warm H 2 have been carried out with the Infrared Spectrograph (IRS; Houck et al. 2004), targeting H 2 midinfrared (MIR) lines (H 2 0-0 S(1), S(2) and S(3) transitions) in 30Dor at low-and high-resolution (Lebouteiller et al. 2008;Indebetouw et al. 2009). However, the low-resolution observations mostly provide upper limits only and the high-resolution observations are only available for a limited number of pointings. These limitations make it challenging to determine the distribution of the total molecular gas mass using direct H 2 observations from IRS.

More recently, H 2 1-0 S(1) has been observed in the entire 30Dor nebula at 1 resolution with the CTIO 4-meter Blanco Telescope (Yeh et al. 2015). These high-resolution observations of the warm molecular gas in 30Dor reveal that the H 2 emission originates from the PDR, and no evidence for shock excitation is found. The H 2 mass inferred from these observations is obviously a strong lower limit on the total molecular gas mass as it does not include the mass associated with other levels of H 2 and does not trace the reservoir of cold molecular gas associated with star formation.

H 2 gas mass determination by dust

Since direct H 2 observations are challenging, the FIR emission from dust remains a commonly used tracer of the gas mass, through Centre: Fraction of the CO-dark gas calculated by comparing the predicted H 2 column density from the Meudon PDR model and the estimation through the CO(1-0) observations and a galactic X CO,MW factor. The size of the circles indicates the total H 2 mass surface density, between 30 and 220 M pc -2 . Right: Ratio of the conversion factors X CO,30Dor over X CO,MW . The minima of X CO,30Dor correspond to the peaks of 12 CO (1-0) (black contours). X CO,30Dor increases close to the location of R136 (red circle) and in the North-West region, where [C II] is bright and were P increases as well.

fitting of the dust spectral energy distribution (SED; e.g. Thronson 1988;Israel 1997;Dame et al. 2001;Leroy et al. 2007Leroy et al. , 2009Leroy et al. , 2011;;Gratier et al. 2010;Bolatto et al. 2011, in particular in dwarf galaxies). The gas-to-dust mass ratio (G/D) is believed to be fairly constant within a region and can be calculated as the ratio of the neutral gas mass (M H I ) to the dust mass (M dust ) at locations where neutral hydrogen dominates. Alternatively, the Galactic G/D value of ∼ 158 (Zubko et al. 2004) can be scaled with metallicity following Rémy-Ruyer et al. (2014). The H 2 mass is then calculated as:

M H 2 = M dust × (G/D) -M H I (2)
in an atomic or molecular gas dominated medium. This method has the advantage of being independent of CO emission observations, but the determination of the total gas mass from the dust mass is susceptible to variations of G/D. Especially, on small scales, this method relies on the assumption that gas and dust are well mixed.

Even if this assumption is verified, the determination of the total gas mass still depends on the dust model used (e.g. Galliano et al. 2011;Rémy-Ruyer et al. 2014). In addition, while G/D roughly scales with metallicity (e.g. Sandstrom et al. 2013), Rémy-Ruyer et al. (2014) show that there is a large scatter in the G/D-metallicity relation between galaxies and therefore the scaling of G/D with metallicity is uncertain. We apply this method to the 30Dor observations in order to compare the gas mass surface density determined from dust modeling with the results of Section 4.2 obtained with our PDR modeling. We convolve the Spitzer/IRAC (3.6, 4.5, 5.8 and 8.0 µm bands) and MIPS (24 and 70 µm bands) and the Herschel/PACS (100 and 160 µm bands) and SPIRE (250 µm band) photometry data with a Gaussian function in order to degrade the resolution of the observations to 18 (similar to the resolution of SPIRE at 250 µm). We use the phenomenological SED fitting procedure presented in (Galliano et al. 2018, Appendix C) to estimate the amount of dust along the line of sight for each pixel. This method uses a least-square fitter with the AC grain mixture of Galliano et al. (2011) and the diffuse interstellar radiation field intensity distribution from Dale et al. (2001). This model was designed to fit the Herschel broadband photometry of the LMC, and to remain consistent with the ob-served elemental abundances. The free parameters include the total infrared luminosity (L TIR ) between 3 and 1000 µm, the minimum starlight intensity, the difference between the maximum and minimum starlight intensities, the starlight intensity distribution powerlaw index and the PAH-to-total dust mass ratio. The G/D for 30Dor is determined using the broken power-law relationship between G/D and the metallicity defined in Rémy-Ruyer et al. (2014). The abundance ratio in 30Dor is log(n(O)/n(H)) = -3.75 , as measured by Pellegrini et al. (2011). This is in the regime where G/D scales linearly with the metallicity (this relation becomes non-linear for Z 0.1; Rémy-Ruyer et al. 2014) and we derive (G/D) 30Dor = 450. Using this value of G/D and the H I observations by Kim et al. (2003) in equation 2, we find a mass surface density of molecular gas on average about a factor of 3 higher than that measured with our PDR modeling, with similar spatial variations (see Figure 6). This leads to a range of X CO,30Dor factor of 3 -12 × 10 21 . The higher values of molecular gas mass surface density and X CO,30Dor found here in comparison to Section 4.2 can possibly be explained by an effective lower G/D in highly irradiated regions compared to the relationship calibrated by Rémy-Ruyer et al. ( 2014) (here by a factor of ∼ 3), or by the fact that we have considered an atomic or molecular gas dominated medium and we have not removed the possible contribution from the ionised gas. In addition, the above equation 2 only takes into account the scaling of G/D with metallicity, but not with gas density, which results from the accretion of material onto grains. This accretion in the dense ISM leads to an increase of the dust mass (up to a factor of 2), but also an increase of the emissivity of the grains by up to a factor of ∼ 3 (Köhler et al. 2015). This could therefore potentially lead to an overestimate of the dust mass by up to a factor of 6. Despite these caveats, the determination of the total H 2 mass surface density from dust modeling also strongly suggest the presence of a large reservoir of CO-dark molecular gas.

Other literature measurements on local scales

Another classical method to estimate the mass of molecular gas is to apply the virial theorem to resolved CO clouds. This method is limited to the nearest galaxies due to resolution and sensitivity requirements. It has been in particular applied in the LMC (Israel et al. 2003;Fukui et al. 2008;Hughes et al. 2010;Wong et al. 2011;Pineda et al. 2009), leading to X CO ∼ 8 × 10 20 cm -2 (K km s -1 ) -1 , relatively similar to the dust-based estimate from Roman-Duval et al. (2014, see Section 5.2.4). However, this method likely leads to underestimation of the total H 2 mass, and even more so in lower metallicity environments, as the CO core is likely to be smaller than the entire H 2 cloud, due to the presence of a larger PDR envelope. As a result, virial measurements only consider the gas bright in CO and will not capture the potential CO-dark molecular gas reservoir. Besides, this method requires the CO molecular clouds to be resolved, and makes the assumption that the CO cores are in virial equilibrium, which is still debated (e.g. Ballesteros-Paredes 2006). On smaller scales in 30Dor, Indebetouw et al. (2013) observed a ∼ 1 × 1 region in 12 CO(2-1) and 13 CO(2-1) with ALMA, centred on α = 5 h 38 m 48 s and δ = -69 • 0.4 50 and also covered by our PACS and FIFI-LS observations. They estimate the molecular mass from the 12 CO and 13 CO intensities and calculate an α CO factor of 8.9 M pc -2 (K km s -1 ) -1 (corresponding to X CO ∼ 4 × 10 20 cm -2 (K km s -1 ) -1 ) for the mapped region. However, we note that these observations have high spatial resolution (∼ 1 pc), and Indebetouw et al. (2013) measure specifically the H 2 mass associated with the CO clumps. The X CO factor including H 2 possibly existing between the CO clumps is therefore expected to be higher.

Comparison with literature measurements on global scales

On larger scales, Roman-Duval et al. ( 2014) use a single temperature blackbody modified by a broken-law emissivity to calculate the dust mass from the SED over the entire LMC, at a resolution of 15 pc. They find an upper limit of X max CO = 6 × 10 20 cm -2 (K km s -1 ) -1 for the LMC. Our values are about an order of magnitude higher than this upper limit, which can easily be explained by the fact that we are specifically focusing on a few tens of parsecs around R136, where the radiation field is the strongest and where the differential dissociation between CO and H 2 is expected to be the most noticeable. In addition, the particular parameterisation of the dust model of Gordon et al. (2014) used by Roman-Duval et al. (2014) effectively results in an emissivity in the LMC significantly lower than in the Milky Way (i.e. a lower submillimeter excess in the LMC than in the MW), which leads to lower dust masses than measured here. At lower spatial resolution, Israel (1997) presents dust-based estimates of the total molecular gas mass in various starforming regions of the LMC and Small Magellanic Cloud (SMC). Their measured X CO for the entire 30Dor region (over 235 pc) is 8.4×10 21 cm -2 (K km s -1 ) -1 . This is again much higher than what would be expected from a simple scaling of the fiducial Milky Way value by the half-solar metallicity of the LMC and is consistent with our dust-based measurement, and with the upper range of our

The CO-dark molecular gas mass in 30 Doradus 11 PDR-based measurement. The strong variation between the average X CO factor found for the LMC and X CO,30Dor measured here and in Israel (1997) indicates that it is not only the half solar metallicity ISM but the combination with the strong radiation from the SSC that together drives the large fraction of CO-dark molecular gas in our observations.

Another way to constrain the global H 2 mass is by assuming a constant molecular gas depletion time. This method is based on the observed constant scaling between H 2 surface density and the SFR surface density on large scales for nearby spiral galaxies. Schruba et al. (2012) estimate the CO-to-H 2 conversion factor in several dwarf galaxies of the IRAM HERACLES survey (Leroy et al. 2009) making the assumption of a constant depletion time. They find α CO ∼ 100 M pc -2 (K km s -1 ) -1 (i.e. X CO ∼ 4.6 × 10 21 cm -2 (K km s -1 ) -1 ) for the integrated LMC. However, similarly to X CO and G/D, the depletion time might vary with the environment. In addition, this method can only be applied on large scales, where observations average over a sufficiently large number of individual starforming regions. On scales smaller than ∼ 1 kpc, the correlation between H 2 surface density and the SFR surface density breaks down, due to the insufficient statistical sampling of regions following independent evolutionary cycles (e.g. Schruba et al. 2010;Kruijssen & Longmore 2014). Specifically, the ratio between the H 2 mass surface density and the SFR surface density of a single star-forming region (i.e. its instantaneous depletion time) varies a priori between 0, before the onset of star formation, and infinity, after the dispersion of the gas by feedback. This approach can therefore not be applied to determine the H 2 mass in an individual star-forming region.

Environmental variations of the CO-dark gas mass

In this section, we look for environmental factors that could explain the spatial variation of the CO-dark gas throughout 30Dor. We first calculate the average distance of the PDR in each pixel to the centre of R136, by comparing the emitted radiation field from the SSC, G stars , to the incident radiation field on the PDRs, G UV , which we derive from the PDR modeling as described in Section 4.1. We use for this the same approach as in Chevance et al. (2016). We use the brightness temperature and position in the plane of the massive stars in R136 from the literature, and assume that they are randomly distributed along the line of sight, within a sphere of 6 pc radius. We then propagate the ionising flux from each star, using a 1/R 2 relation (with R the physical distance from the centre of the cluster), to calculate the average G stars in a cube centred on R136, assuming no extinction along the line of sight. For each pixel of the map (i.e. each observed line of sight), we then estimate the distance z of each pixel to the plane of R136 by matching G UV to one of the possible values of G stars along this line of sight. Finally, we calculate the 3D distance of each pixel to the centre of R136 as R = z 2 + d 2 , where d is the projected distance between a given pixel and the centre of the source.

One might expect an increase of the CO-dark gas fraction with the incident radiation field G UV (or equivalently a decrease with the physical distance to R136). However, such a correlation is not observed in our results. While this may seem surprising, it can be explained by the observed anti-correlation between the CO-dark gas fraction and A max V , in combination with the observed decrease of A max V with distance to R136 (respectively decreasing G UV ; see Figure 7), such that clouds with smaller shielding columns are observed towards the outskirts of the region. The first of these cor-relations is indeed expected from the model: the deeper the cloud, the better CO will trace the total H 2 reservoir. On the contrary, for a small size cloud with a low A V , the PDR layer including H 2 but no CO represents a non negligible part of the cloud. In this context, the fact that we observe a decrease of A max V with increasing distance to R136 (respectively decreasing G UV ) likely balances the expected increase of the CO-dark gas fraction with the intensity of the radiation field.

Through the strong correlation with A max V as noted above, there is a direct correlation between the CO-dark gas fraction and the ratio [C II]/CO(1-0). This correlation is investigated further in Madden et al. (in prep.) by quantifying the total molecular gas mass on global scales in the Dwarf Galaxy Survey sources. 88 µm, previously inaccessible with Herschel, ranges between 0.48 and 1.84 in 30Dor, and is a good probe of the electron density in the low-density ionised gas. Overall, we measure a maximum density of ∼ 600 cm -3 around R136. This is in agreement with previous measurements of the electron density of the ionised gas in this region (e.g. Indebetouw et al. 2009;Pellegrini et al. 2011;Chevance et al. 2016), although the different line ratios used in these studies might probe different phases of the ionised gas.

Following a similar approach as in Chevance et al. (2016), we have determined the spatial distribution of the pressure and radiation field in 30Dor, by fitting our observations with the Meudon PDR model. We find that the incident radiation field G UV ranges between 7 × 10 2 and 7 × 10 3 (in units of the standard radiation field defined in Mathis et al. 1983) and the pressure P ranges between 4.3 × 10 5 and 1.6 × 10 6 K cm -3 , at a resolution of 43 (∼ 10 pc). By making use of the ratio [C II]/CO(1-0), which is sensitive to the depth of the PDR clouds, we determine a maximum total visual extinction A max V of 2 to 3 mag. In 30Dor, the intense radiation field induced by the massive excitation source R136 on the surrounding half-solar metallicity gas creates an extreme environment. Based on the physical properties of the gas determined through the PDR modeling, and by comparing to CO(1-0) observations, we estimate that between 75% and 97% of the H 2 would be undetected, and therefore can be considered as "CO-dark" molecular gas, if a standard X CO,MW factor of 2 × 10 20 cm -2 (K km s -1 ) -1 would be applied in this region. As a result, a more appropriate CO-to-H 2 conversion factor around R136 is X CO,30Dor ∼ 1-4×10 21 cm -2 (K km s -1 ) -1 , in agreement with previous determinations of this conversion factor in 30Dor and other regions of the LMC using different methods. We note that "CO-dark" gas here does not presume whether the CO lines are simply too faint relatively to the Milky Way environment, or if a significant fraction of the H 2 gas is CO-free and co-spatial with ionised and neutral carbon. This question is investigated in more detail in Madden et al. (in prep.). The spatial variation of the CO-dark fraction throughout the mapped regions is mostly driven by the amount of visual extinction A max V (i.e. the depth of the clouds). This strong dependence seems to balance any expected trend with the incident radiation field or the distance to R136. The influence of the SSC on the neutral and molecular gas is therefore still dominant out to several tens of parsecs away from the centre of the cluster. Future observations of other star-forming (or more diffuse) regions in the LMC and SMC would allow us to disentangle the combined effects of the radiation field and the metallicity on the properties of the ISM (e.g. the pressure, porosity, and size of the PDR clouds) and ultimately probe the environmental dependence of the CO-dark gas fraction.

  The ratio [C II]/[O III]88 µm shows low values (between 10 -2 and ∼ 0.5 at the resolution of ∼ 4 pc achieved by the FIFI-LS observations), with the lowest values reached around the peaks of the ionised gas probed by the [O III] lines. The ratio [C II]/[O III]52 µm ranges between 10 -2 and ∼ 0.3. These low values reached by the ratio [C II]/[O III] (even at large distances from the centre of R136) highlight the fact that high

Figure 1 .

 1 Figure 1. Intensity of the FIR fine structure lines [O III] at 52 µm, [O III] at 88 µm, [O I] at 145 µm and [C II] at 158 µm in Wm -2 sr -1 , observed by SOFIA/FIFI-LS. The red circle has a 10 diameter and indicate the location of R136. For [O I] 145 µm, the map is a combination of the FIFI-LS and PACS observations(see text). The outline of the PACS and FIFI-LS fields of view are indicated in cyan and black respectively. The white pixels in the maps are the ones for which the emission line fit has a a signal-to-noise ratio lower than 3. These are not taken into account in our analysis. We note the layered structure of the gas, with the peak of the ionised gas ([O III]) being closer to R136 than the peak of the neutral gas ([C II] and [O I]). The latter also corresponds to the peak of the molecular gas mass surface density as determined in Section 4.2.

  These values are on average larger than the density range of 10-100 cm -3 probed in the vicinity of the peak of the [O III] lines in Chevance et al. (2016), using the ratio [N II] 122 µm / [N II]

Figure 2 .

 2 Figure 2. Ratio of the FIR fine structure lines [C II] 158 µm / [O III] 88 µm at 18 (∼4 pc) resolution. We note that the color bar is in logarithmic scale. The black contours represent the [C II] line intensity and the white contours represent the [O III] 88 µm line intensity. The yellow circle indicates the location of R136. The ratio [C II]/[O III] 88 µm is below 1 in the entire area covered by FIFI-LS, revealing a porous medium bathed in strongly ionising radiation field tens of parsecs away from the excitation source R136.

V

  constrained byLee et al. (2019) (∼ 1.5-2 mag), relying mostly on the ratio [C II]/ [C I]. This is not surprising given the expected structure of PDR regions, where CO lines are emitted deeper into the clouds than [C I] lines, and as already shown in fig.17ofChevance et al. (2016). The average A max

Figure 3 .

 3 Figure 3. Left: Map of the observed ratio [O III] 52 µm/ [O III] 88 µm (restricted to the pixels with SNR > 5). The black contours represents the [O III] 88 µm intensity. The yellow circle indicates the location of R136. Right: Theoretical ratio [O III] 52 µm/ [O III] 88 µm at the temperatures of 5 000 (dotted line), 10 000 (solid line) and 15 000 K (dashed line) as a function of the electron density n e . The normalised histogram of the observed values of the ratio [O III] 52 µm/ [O III] 88 µm is shown on the right. The normalised histogram of the inferred densities (using the theoretical prediction at 10 000 K) is shown on the top.For the lowest value of the observed ratio ( 0.65), only an upper limit can be determined for the density. The values in the grey area should therefore be considered as upper limits at ∼ 30 cm -3 .

Figure 4 .

 4 Figure 4. Constrained parameter of the Meudon PDR model, using the observed ratios ([O I] 145 + [C II] 158 )/L FIR , [O I] 145 /[C II] 158 and [C II] 158 /CO(1-0), at a resolution of 43 The spatial distributions of the pressure P (left; in K cm -3 ), the incident radiation field (centre; in units of the Mathis field) and the total visual extinction A max V (right; in magnitudes) are presented, along with the location of the SSC R136 (red circle). We overplot the [C II] intensity as black contours in the central panel to guide the eye.

Figure 5 .

 5 Figure 5. Left: Mass surface density of H 2 (in M pc -2 ) predicted by the PDR model. The black contours represent the 12 CO (1-0) emission. The red circle indicate the location of the SSC R136. The maximum of the H 2 mass surface density corresponds to the peak of CO intensity, as well as of [C II] and [O I].Centre: Fraction of the CO-dark gas calculated by comparing the predicted H 2 column density from the Meudon PDR model and the estimation through the CO(1-0) observations and a galactic X CO,MW factor. The size of the circles indicates the total H 2 mass surface density, between 30 and 220 M pc -2 . Right: Ratio of the conversion factors X CO,30Dor over X CO,MW . The minima of X CO,30Dor correspond to the peaks of 12 CO (1-0) (black contours). X CO,30Dor increases close to the location of R136 (red circle) and in the North-West region, where [C II] is bright and were P increases as well.

Figure 6 .

 6 Figure 6. Mass surface density of molecular gas from the dust modeling. In black are the contours of the molecular mass surface density from the PDR modeling at 100-150-200 M pc -2 . The red circle indicates the location of R136.

  paper, we have presented new observations of the FIR lines [O III] 52 µm and 88 µm, [O I] 145 µm and [C II] 158 µm in the massive star forming region 30 Doradus of the LMC, observed during Cycles 4 and 5 of SOFIA with FIFI-LS. These observations extend and complement previous Herschel/PACS observations of this region and allow us to explore the influence of the distance to the SSC R136 on the physical properties of the gas. The spatially extended [O III] emission and the overall low values of the ratio [C II]/[O III] 88 µm reveal a porous medium where high-energy photons emitted by R136 propagate over several tens of pc before interacting with the ISM. The ratio [O III] 52 µm/[O III]

Figure 7 .

 7 Figure 7. Top left: G UV as a function of the physical (deprojected) distance of the PDR clouds to the centre of R136. Clouds closer to the centre of the SSC receive a stronger radiation field. The dotted line represents a 1/D 2 relation for the intensity of a central source with distance D, scaled to the mean emitted flux 30 pc away form R136 (see text and Chevance et al. 2016). The fact that the data points lie slightly above this relation come from the spatial extension of the source, taken into account to infer the distances. Top right: CO-dark mass fraction as a function of the physical distance of the PDR clouds to the centre of R136. No significant correlation is observed. Bottom left: CO-dark mass fraction as a function of the extinction in individual clouds. The observed anti-correlation is expected from the PDR model, as for clouds of smaller depth the (CO-dark) H 2 layer in the PDR represents a larger fraction of the cloud. Bottom right: Extinction in individual clouds as a function of the incident radiation field. The observed correlation between A V and G UV compensates the expected increase of the CO-dark gas fraction with G UV .

Table 2 .

 2 Measured CO-to-H 2 conversion factor in the LMC and in star-forming (SF) regions of the LMC and SMC.

	Reference	Method	X CO [ cm -2 (K km s -1 ) -1 ]	Object	Physical size
		PDR modeling	0.8 -4 × 10 21	30Dor	10 pc
	This study	PDR modeling	2 × 10 21	30Dor	∼ 50 pc
		Dust modeling	3 -12 × 10 21	30Dor	15 pc
	Indebetouw et al. (2013)	Virial equilibrium	4 × 10 20	CO clumps in 30Dor	1 pc
	Fukui et al. (2008)	Virial equilibrium	7 × 10 20	CO clouds in LMC	40 pc
			8.4 × 10 21	30Dor	235 pc
	Israel (1997)	Dust modeling	1.3 × 10 21	LMC SF regions	235 pc
			12 × 10 21	SMC SF regions	270 pc

The Meudon PDR model is public and available online at the following address: http://ism.obspm.fr.MNRAS 000, 1-14 (2020)

MNRAS 000, 1-14 (2020)

ACKNOWLEDGEMENTS

We thank the referee, Neal Evans, for constructive feedback that improved this paper. The authors would like to thank Annie Hughes for providing the CO(1-0) data. MC and JMDK gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through an Emmy Noether Research Group (grant number KR4801/1-1) and the DFG Sachbeihilfe (grant number KR4801/2-1). JMDK gratefully acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme via the ERC Starting Grant MUSTANG (grant agreement number 714907). DC 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). The authors would like to thank the French Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP cofunded by CEA and CNES, for supporting our SOFIA flights and collaboration with SOFIA projects. This research was made possible through the financial support of the Agence Nationale de la Recherche (ANR) through the programme SYMPATICO (Program Blanc Projet NR-11-BS56-0023) and through the EU FP7. This The CO-dark molecular gas mass in 30 Doradus 13 work is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart.