© The Authors 2024

1 Introduction

Molecular clouds (MCs) are the densest and coldest regions of the interstellar medium (ISM) and, as the name suggests, are defined as those regions where hydrogen exists mainly in molecular form. Such environments are also the regions where dense cores and stars form. The formation and evolution of MCs has been studied with numerical simulations in a large number of recent works (see e.g. the reviews by Henshaw et al. 2023; Chevance et al. 2023; Hacar et al. 2023, and references therein). Once star formation begins, it is of particular interest to explore the effects that stellar feedback has on the evolution and possible disruption of MCs. Here, the three most important forms of stellar feedback are supernova explosions, stellar winds, and ionising radiation. In particular, the latter can contribute to cloud dispersal as it decreases the density of the cloud in and around sites of star formation. As a consequence, this also further inhibits star formation, contributing to the low star formation efficiency generally observed (e.g. Krumholz et al. 2014; Peters et al. 2017; Rathjen et al. 2021, 2023). In addition, photons emitted by stars dissociate or ionise the surrounding gas. The first theoretical study of this phenomenon, performed by Strömgren (1939), showed that, under the assumption of uniform density and isotropic conditions, stellar radiation produces a shell where hydrogen is ionised; the radius r of this shell is $r\propto Q_0^{1/3}{n}^{-2/3}$, where Q₀ is the ionising flux and n the gas density. These portions of MCs are called HII regions. Later on, Spitzer (1978) calculated the expansion velocity of the ionisation fronts, finding r ∝ t^4/7 (see also Bisbas et al. 2015; Williams et al. 2018). Following on from these fundamental works, ionising radiation and its role in MCs have been extensively studied in a large number of numerical works (Hosokawa & Inutsuka 2006; Krumholz & Matzner 2009; Fall et al. 2010; Murray et al. 2010; Dale et al. 2012; Walch et al. 2012; Girichidis et al. 2016; Gatto et al. 2017; Haid et al. 2019; Scheuermann et al. 2023, and many more).

The radiation from the stars has a substantial impact on the chemistry of the clouds. At the boundaries of the aforementioned HII regions, the so-called photodissociation regions (PDRs) form, where the chemistry is dominated by far-ultraviolet (FUV) photons and their progressive attenuation (Wolfire et al. 2022, and references therein). In this context, one of the brightest emission lines originating from star-forming regions is the [CII] line at 157.7 µm (e.g. Stacey et al. 1991; Brauher et al. 2008), which is one of the main coolants in the warm ISM (Tielens & Hollenbach 1985; Stacey et al. 1991; Hollenbach & Tielens 1999; Stutzki 2001; Röllig et al. 2006; Wolfire et al. 2003; Ossenkopf et al. 2013; Beuther et al. 2014; Pineda et al. 2013, 2014; Klessen & Glover 2016). The [CII] emission has been extensively observed towards HII regions (Schneider et al. 2012; Goicoechea et al. 2015; Röllig et al. 2016; Pabst et al. 2017, 2020; Schneider et al. 2020; Tiwari et al. 2021; Luisi et al. 2021; Beuther et al. 2022; Kabanovic et al. 2022), and is found to generally highlight a bubble structure with bright rims and very little emission from the inner part. [CII] is also frequently used as a tracer of the star-formation rate in galaxies (e.g. Stacey et al. 1991; Boselli et al. 2002; De Looze et al. 2014; Herrera-Camus et al. 2015; Sutter et al. 2019; Bisbas et al. 2022) as well as a tracer in cloud–cloud interactions (Bisbas et al. 2017; Schneider et al. 2023).

The far-infrared (FIR) dust continuum emission is the thermal emission from dust grains with temperatures in a typical range of 10-50 K. It is usually defined as the integrated emission in the wavelength range from 40 to 500 µm (e.g. Sanders & Mirabel 1996; Goicoechea et al. 2015; Pabst et al. 2021). It is also a meaningful measure for understanding HII regions, as it is intimately associated with star formation, as shown in numerous recent studies (Brauher et al. 2008; Graciá-Carpio et al. 2011; De Looze et al. 2011; Casey et al. 2014; Herrera-Camus et al. 2018; Pabst et al. 2022; Dunne et al. 2022).

The ratio between [CII] and FIR emission is of particular interest for understanding stellar feedback in MCs as both are thought to trace star formation. However, extragalactic observations show that this ratio, which is typically of the order of 10⁻²–10⁻³, decreases with increasing infrared luminosity across galaxies (Malhotra et al. 1997; Luhman et al. 1998; Malhotra et al. 2001; Luhman et al. 2003; Casey et al. 2014; Díaz-Santos et al. 2017; Smith et al. 2017; Herrera-Camus et al. 2018). The origin of this so-called [CII] deficit is still unclear, although different hypotheses have been proposed, including AGN contributions to FIR emission (Sargsyan et al. 2012), further ionisation of C⁺ into C²⁺ (Abel et al. 2005), fine-structure lines overcoming [CII] as coolants (Luhman et al. 2003, for instances), and saturation of the [CII] line (Muñoz & Oh 2016; Rybak et al. 2019; Bisbas et al. 2022). Evidence for a relationship between the [CII] deficit and the star formation rate (SFR) has also been presented, both theoretical (Narayanan & Krumholz 2017) and observational (Stacey et al. 2010; Smith et al. 2017; Hu et al. 2019). To get a better understanding of the origin of the [CII] deficit, resolved studies of individual HII regions and MCs are required. However, so far only a few have been reported in the literature, and these explicitly consider both the [CII] and FIR emission on MC scales (Goicoechea et al. 2015; Pabst et al. 2021).

Here, we extend our previous work on synthetic emission maps from MCs (Ebagezio et al. 2023) and provide a simulation counterpart to the two [CII] and FIR observations on MC scales mentioned above in order to enable a comparison between real and simulated data, and to take advantage of our knowledge of the simulation input to obtain more detailed insights into the physical phenomena that cause the [CII] deficit. To this end, we used the SILCC-Zoom simulations, which include stellar feedback and a chemical network to properly account for the ionised carbon from which the [CII] emission emerges (Seifried et al. 2017; Haid et al. 2019). The paper is organised as follows. In Section 2 we present the numerical simulations of the MCs that serve as a basis for our analysis, as well as the post-processing and analysis techniques that we use. Next, in Section 3 we describe our results: we investigate the dust temperature, and the [CII] and FIR emission maps and their ratio. These results are discussed in Section 4, where we present a more detailed investigation of the causes of the [CII] deficit, and make comparisons to observational data. In Section 5 we summarise our work and draw conclusions.

2 Numerical methods

This work is a further analysis of the synthetic emission maps described in Ebagezio et al. (2023). Therefore, we refer to this paper and to the related SILCC-Zoom simulation papers (Seifried et al. 2017; Haid et al. 2019) for a detailed description of the numerical simulations, the post-processing routines, and the radiative transfer. Here, we briefly summarise the most important concepts of these simulations, and describe the additional processes not incorporated into the aforementioned works.

2.1 SILCC-Zoom simulations

The simulated MCs that we use in this work are part of the SILCC-Zoom project (Seifried et al. 2017; Haid et al. 2019, Walch et al., in prep.). This consists of zoom-in simulations of specific regions of the larger-scale simulations within the SILCC project (see Walch et al. 2015; Girichidis et al. 2016, for details).

The SILCC project aims to simulate the life-cycle of molecular clouds in typical solar neighbourhood conditions. To do so, it models a region of a stratified galactic disc of 500 pc × 500 pc × ±5 kpc in size using the adaptive mesh refinement code FLASH 4.3 (Fryxell et al. 2000; Dubey et al. 2008). The maximum allowed resolution is 3.9 pc. Both hydrodynamical (HD) and magnetohydrodynamical (MHD) runs are simulated within the project. The chemistry is modelled using a chemical network for H⁺, H, H₂, C⁺, O, CO, and e⁻ (Nelson & Langer 1997; Glover & Mac Low 2007a,b; Glover et al. 2010), also including the most important heating and cooling processes of the gas. We assume solar metallicity, corresponding to an elemental abundance of carbon and oxygen of 1.4 × 10⁻⁴ and 3.16 × 10⁻⁴, respectively, relative to total hydrogen (Sembach et al. 2000). The chemical network is applied on-the-fly in the simulations, that is, the chemistry is evolved together with the hydrodynami-cal evolution. In Ebagezio et al. (2023) we show that having an on-the-fly chemical network permits to obtain significantly more accurate chemical abundances with respect to post-processed chemical networks which assume chemical equilibrium, which is in agreement with findings of for example Glover & Mac Low (2007b); Seifried et al. (2017); Hu et al. (2021); Seifried et al. (2022). The interstellar radiation field (ISRF) is set to G₀ = 1.7 in Habing units (Habing 1968) and the cosmic ray ionisation rate is 3 × 10⁻¹⁷ s⁻¹ with respect to atomic hydrogen. We also consider the attenuation of the ISRF using the TREERAY/OPTICALDEPTH module (Wünsch et al. 2018). We consider simulations with and without magnetic fields. The magnetised runs have a magnetic field B initialised to be parallel to the x-direction with a strength of B_x,0 = 3 µG in the galactic mid-plane. The dust-to-gas ratio is fixed to be 1/100 everywhere in the simulations. We refer to Walch et al. (2015) for further details concerning the initial conditions.

At the beginning of the SILCC simulations and up to a given time t₀ (see Table 1), supernova explosions drive turbulence. At t₀, further explosions are stopped and local over-densities are already visible. At this point, we focus on some of these over-dense regions (i.e. those which are part of the SILCC-Zoom project), zooming-in onto them with a resolution as high as 0.12 pc while keeping the larger-scale environment on a lower resolution. In this paper we use two HD and two MHD clouds. We refer to them as MC1-HD, MC2-HD, MC1-MHD, and MC2-MHD, respectively (see Table 1).

All runs include stellar feedback. We use sink particles to model the formation and evolution of stars (see Haid et al. 2018, 2019, for more details). The ionising radiation feedback from each massive star is treated with TreeRay (Wünsch et al. 2021) and is fully coupled to the chemical network (Haid et al. 2018, 2019). We note that currently we do not include non-ionising stellar radiation below 13.6 eV, which would add to the ISRF, that is G₀ (see also Appendix A).

2.2 Chemical post-processing

Before the actual synthetic emission maps are produced, we post-process the data to model physical phenomena which are not included in the simulation itself. In particular, we include: (i) CO freeze-out, (ii) splitting of H₂ into para- and ortho-H₂, (iii) microturbulence, (iv) C⁺ → C²⁺ thermal ionisation, and (v) C⁺ → C²⁺ photo-ionisation in regions affected by stellar feedback as a consequence of stellar radiation. The effects (i) – (iv) are implemented in a pipeline¹ developed by Pierre C. Nürnberger, which couples the FLASH simulation output to RADMC-3D input required for the radiative transfer (Section 2.3). As pointed out in Ebagezio et al. (2023), the photoionisation (v) has a large impact on the abundance of C⁺ within the HII regions and, as a consequence, also on the synthetic [CII] maps of those regions. It is thus essential to perform this postprocessing step (v) in order to accurately model the [CII] within the HII regions.

For this purpose, we sample the gas density, gas temperature, the stellar temperature and bolometric luminosity (the latter is obtained from the ionising radiation flux modelled in the 3D, MHD simulations) over the typical ranges occurring in our simulations, and build a grid of PDR models with CLOUDY (Ferland et al. 2017). For this post-processing step we use the entire frequency range also covering that of the ISRF, that is the local G₀ generated by nearby stars, to accurately calculate C⁺ and C²⁺ abundances. Next, using the information on the stellar content in our simulations, we obtain from this pre-calculated PDR database the updated fractional abundance of C⁺ that takes into account the C⁺ → C²⁺ photo-ionisation. We then use it to replace in each cell the original abundance of C⁺ given by FLASH. This update step is only done in regions where the 3D simulations predict the presence of ionising radiation coming from the stars. More details on this post-processing routine are provided in Ebagezio et al. (2023).

2.3 Radiative transfer

We use the RADMC-3D code (Dullemond et al. 2012) to do the radiative transfer calculations of the simulated clouds and to obtain synthetic emission maps. RADMC-3D is an open-source, 3D radiative transfer code, capable of performing both line and dust continuum radiative transfer calculations. We create synthetic emission maps of the [CII] 158 µm line and of FIR emission between 3 and 1100 µm. We define t as the simulation time since the very beginning of the SILCC simulations and t_evol as t_evol = t − t₀. We consider the evolutionary stages at t_evol = 2, 3, and 4 Myr for the HD simulations, and t_evol = 3–5 Myr for the MHD simulations, as the MHD clouds generally evolve somewhat slower (Seifried et al. 2020; Ebagezio et al. 2023). Table 1 gives an overview of the simulations and the times at which the synthetic emission maps are produced. All emission maps are calculated for lines of sight (LOS) along the x-, y-, and ɀ-axis. Furthermore, for MC1-HD and MC2-HD we produce maps at intervals of 100 kyr along the ɀ-axis. We use these highly time-refined emission maps in Section 3.5 to be able to better study the evolution of the considered characteristics.

2.3.1 Line radiative transfer

In the following, we briefly describe the method to obtain the [CII] emission; for more details we refer to Franeck et al. (2018) and Ebagezio et al. (2023). In order to accurately model the [CII] line emission with RADMC-3D, we include microturbulence assuming that microturbulent broadening is as strong as the thermal broadening. Furthermore, we use the large velocity gradient (LVG) approximation (Ossenkopf 1997; Shetty et al. 2011a,b) to calculate the level population. We consider 201 equally spaced velocity channels for a velocity range of ±20 km s⁻¹ around the [CII] rest frequency, which results in a spectral resolution of dv = 0.2 km s⁻¹. This velocity range and spectral resolution enable us to capture the contribution of Doppler-shifted emission (Franeck et al. 2018). As we use the LVG approximation, we need to specify the collisional rates for C⁺. For this, we rely on the Leiden Atomic and Molecular DAtabase (LAMDA, Schöier et al. 2005). We consider para-H₂, ortho-H₂, H, and e⁻ as collisional partners for C⁺. We also model the effect of the Cosmic Microwave Background by assuming a background black body radiation at T = 2.725 K (Fixsen 2009), which, however, only adds a negligible contribution both to the level population of C⁺ and the obtained [CII] intensity.

2.3.2 Dust continuum radiative transfer

For calculating the dust continuum emission, the dust temperature T_d must be known. The FLASH simulation data contain information about dust temperature set by the ISRF and the collisional coupling of dust and gas. Its calculation does, however, currently not account for the radiative contribution of the stars formed in the simulations (see Klepitko et al. 2023, for our recent implementation to couple this radiation to the dust). Therefore, in a first step we use RADMC-3D to calculate the dust temperature resulting from the present stars using their luminosity and stellar surface temperature. In a second step, we consider for each simulation cell the maximum of the dust temperature given by FLASH and the one given by RADMC-3D. In this way, we properly account for the fact that dust can also be heated via energy transfer from the gas due to collisions as well as heating from the ISRF (which both are accounted for in FLASH). As this effect is not included in RADMC-3D, this could lead to rather low dust temperatures in regions with very low stellar radiation, whereas in regions with a strong stellar radiation, the dust temperature is mainly set by this radiation. Hence, by taking the maximum of the two available dust temperatures, we consider the effect of all dust heating sources in our T_d calculation.

RADMC-3D calculates T_d using a Monte Carlo method described in Bjorkman & Wood (2001) with the continuous absorption method of Lucy (1999). For this purpose, information on absorption and scattering coefficients for different wavelengths are needed. These coefficients depend mainly on the composition and size of the grains. In general, the extinction coefficient κ scales roughly as κ ∝ λ⁻² at long wavelengths, but varies significantly at shorter wavelengths. We use the values from Weingartner & Draine (2001), which aim to model a mixture of carbonaceous and silicate grains typical of MCs under solar neighbourhood conditions (their model B with R_V = 4.0). For the calculation of T_d we use 10¹⁰ photon packages (see Appendix B for a resolution study).

Once the dust temperature is calculated, we use RADMC-3D to calculate the FIR dust continuum emission. We produce maps of 101 wavelengths between 3 and 1100 µm, equally spaced in logarithmic space. We then compute the integrated FIR luminosity by integrating in the range 40–500 µm as this is used in different observational works (e.g. Sanders & Mirabel 1996; Goicoechea et al. 2015; Pabst et al. 2021), allowing for a straightforward comparison with observations.

Fig. 1 Column density maps of MC2-HD from tevol = 2–4 Myr projected along the y-direction (from left to right) showing the effect of star formation (black stars) on the dispersal of the cloud.

3 Results

Before analysing the results concerning the [CII] and FIR emission, we give a brief overview of the general evolution of the clouds. For more details we also refer to Haid et al. (2019). In Fig. 1 we show the evolution of the column density for MC2-HD from t_evol = 2–4 Myr. The snapshots show the filamentary and complex gaseous structure of the MC. Over time, star formation sets in at various places in the cloud, which subsequently disperses the clouds by creating large HII regions which are almost completely devoid of any gas.

In Fig. 2 we show the mass of all massive stars (M > 8 M_⊙) formed for all four clouds considered in this work. In order to better compare different clouds, we here use the time t_Sf, which has passed since the first massive star has formed in a specific cloud at t_Sf,0, that is t_Sf = t − t_Sf,0 (t_Sf,0 is larger than t₀ used to define t_evol, see Table 1). Table 2 shows t_Sf,0 and the relation between t_evol and t_Sf for our clouds². We note, however, that, unless explicitly specified, throughout the paper we use t_evol.

For the clouds MC1-HD, MC2-HD and MC1-MHD there is a phase of rapid star formation within the first ~ 1 Myr followed by a phase with a low star formation rate. Overall, for the clouds including magnetic fields, the star formation rate appears to be lower due to the stabilising effect of the magnetic field against gravity (Girichidis et al. 2018; Seifried et al. 2020; Ebagezio et al. 2023; Ganguly et al. 2023). Particularly, for MC2-MHD, star formation starts with a rather low rate, but continues over a longer time with an indication of acceleration towards late times.

Fig. 2 Time evolution of the mass of all massive stars (M > 8 M⊙) formed in each MC. The clouds without magnetic fields form stars more rapidly.

Fig. 3 Dust temperature PDFs for all four clouds at tevol = 3 (solid lines) and 4 Myr (dashed lines). There is a wide spread of temperatures over about 1 order of magnitude, which increases as star formation proceeds over time.

3.1 Dust temperature

We next analyse the dust temperature, T_d, of our MCs. This is important in order to understand the characteristics of the FIR continuum that we analyse in subsequent steps. First, we consider the distribution of dust temperatures of the four different clouds at t_evol = 3 and 4 Myr (Fig. 3). Overall, there is a wide spread of observed dust temperatures from about 10 to 100 K. Over time the spread increases and the dust becomes on average warmer due to the heating by newly formed stars. This effect is particularly pronounced for the MC1-MHD, where star formation is slightly delayed (see Table 2).

Next, to make predictions comparable to observations, in Fig. 4 we show four projected dust temperature maps, showing the mass-weighted mean dust temperature, 〈T_d〉, averaged along the LOS. Contrary to Fig. 3, we show snapshots for different values of t_evol. In this way, the chosen four snapshots correspond to different stages of star formation in MCs. Moving from top left to bottom right, we show snapshots with progressively more advanced star formation and thus also more and more stars can be found. In general, we find 〈T_d〉 ≃ 15 K in regions not affected by stellar feedback, and 〈T_d〉 ≲ 50 K in regions where stellar feedback is relevant. Feedback-affected regions are confined to a small area of a few pc when star formation has set in recently (t_sf ≤ 1 Myr, top row of Fig. 4), and are extended to a large portion of the cloud (a few 10 pc) when star formation has progressed over a few Myr and several massive stars have been formed (bottom row of Fig. 4).

We emphasise that the LOS averaging reduces the observed range of dust temperatures. The actual dust temperature in the vicinity of stars can be considerably larger than 〈T_d〉, reaching values as high as 100 K and above (see Fig. 3). Overall, the reported dust temperatures match reasonably well with other numerical works. In projected maps, Ali et al. (2018) and Inoguchi et al. (2024) find similar dust temperatures of a few 10 K in HII regions embedded into MCs, also reporting an expansion of regions of warm dust over a few Myr as star formation proceeds. Furthermore, the higher values for T_d in the (unprojected) 3D distribution up to ~100 K (see Fig. 3) are also found in simulations by Kim et al. (2023) and analytical work of Sommovigo et al. (2020), supporting the robustness of our results.

Fig. 4 Maps of the mass-weighted mean dust temperature, 〈Td〉, for four snapshots of the simulated clouds. We show, from top left to bottom right, MCs with progressively more advanced star formation. We overplot the stars present in the simulations, where the size and colour indicate the age and temperature of the stars, respectively. In general, regions not affected by stellar feedback have 〈Td〉 ≃ 15 K, and regions where stellar feedback is important have 〈Td〉 ≃ 50 K. These regions are small (a few pc) and limited to the direct vicinity of the stars in the early stages of star formation (e.g. MC2-MHD, top left), but involve larger parts of the cloud (a few 10 pc) at later stages (e.g. MC2-HD, bottom right).

3.2 [CII], FIR, and [CII]/FIR maps

Next, we analyse the [CII], FIR, and [CII]/FIR ratio maps. In short, we show that the [CII]/FIR ratio is particularly low in regions affected by recent massive star formation, and highlight two different mechanisms that can explain it.

To disentangle these two mechanisms, we show in Fig. 5 maps of [CII], FIR and [CII]/FIR for two different snapshots. The left column shows the [CII] moment 0 maps, I_[CII], the central column gives the distribution of the FIR intensity, I_FIR, and the right column shows I_[CII]/I_FIR. The top row corresponds to MC2-MHD at t_evol = 5 Myr and the bottom row corresponds to MC2-HD at t_evol = 4 Myr. These two snapshots were chosen as they depict the two typical reasons, which lead to a low [CII]/FIR ratio.

For discussing these two reasons, we focus on specific regions in these maps. First, we consider the region in MC2-MHD enclosed by the green square in the top row of Fig. 5. It contains a relatively young star (age of 0.33 Myr), which did not yet evacuate the gas in its surroundings, that is it did not yet create a low-density HII region. As a consequence, the heating due to its radiation causes a centrally concentrated peak of FIR emission with I_FIR ≃ 5 erg s⁻¹ cm⁻² sr⁻¹. This is approximately three orders of magnitude above the values obtained for the surrounding cloud. Conversely, I_[CII] in that area is only larger by a factor of a few compared to the wider-spread, high-intensity region extending to the left, as well as down- and upwards from the centre. As a result, I_[CII]/I_FIR drops in the vicinity of the central star.

A similar behaviour is found in the proximity of the stars on the left side of MC2-HD (bottom row, green square), where again I_[CII] is comparable to the emission coming from the surroundings, whereas I_FIR is significantly higher in the green box compared to its surroundings. Again, this leads to a lower I_[CII]/I_FIR in that area than in its surroundings. We note that overall in MC2-HD I_[CII] is roughly an order of magnitude higher than in MC2-MHD. This is likely due to the stronger excitation of C⁺ caused by more free electrons acting as collisional partners and higher temperatures, which in turn are a consequence of the overall larger amount of massive stars in MC2-HD (see Fig. 2).

Beside the situation found in the two green squares in Fig. 5, there is another reason leading to a low [CII]/FIR ratio. This is indicated in the central part of MC2-HD (bottom row), enclosed by a black rectangle. The region corresponds to an evolved HII region, showing a large extent with low I_[CII]/I_FIR values. Unlike for the other two regions, here this low ratio seems to be rather due to the lack of [CII] emission coming from the HII region, and less due to the enhancement of I_FIR: there is only a moderately enhanced FIR emission, of the order of a few 10⁻¹ erg s⁻¹ cm⁻² sr⁻¹, in the direct vicinity of stars, whereas the rest of the highlighted region is not distinguishable in FIR with respect to the other parts of the cloud. This is due to the dynamic effect of the ionising radiation coming from the more evolved, bright stars in the central part of the highlighted area. It sweeps up the surrounding gas, thus reducing the total column density in the HII region. Hence, even though 〈T_d〉 is larger in that area (see Fig. 4, bottom-right plot), I_FIR is relatively low. Furthermore, I_[CII] is very faint in this area (≃10⁻⁶ erg s⁻¹ cm⁻² sr⁻¹), comparable to the emission coming from very outskirts of the cloud. This is not only due to the aforementioned gas sweep-up, but mainly due to the second photo-ionisation C⁺ → C²⁺ caused by the stellar radiation (see Section 2.2 and Ebagezio et al. 2023, for details on the numerical implementation).

This is demonstrated in Fig. 6, where we show the (projected) fraction of C⁺ with respect to the sum of C⁺ and C²⁺, $$f_{{\text{C}}^{+}}=\frac{N_{{\text{C}}^{+}}}{N_{{\text{C}}^{+}}+N_{{\text{C}}^{2+}}},$$(1)

where N_i is the column density of the corresponding carbon ion. As can be seen, at early stages there is almost no conversion of C⁺ into C²⁺, meaning $f_{{\text{C}}^{+}}\simeq 1$. Contrary to that, at later stages a significant part of the C⁺ has become further ionised. We find that due to the further photo-ionisation, the C⁺ content is in parts reduced down to a few 1% of the overall carbon content, in particular in the HII region marked with the black rectangle, where the [CII] intensity is particularly low (compare with Fig. 5). This clearly shows that a significant part of the reduction in I_[CII] is due to the further photo-ionisation of C⁺ into C²⁺.

We emphasise that the [CII] intensity in this region would be significantly higher without the further photo-ionisation of C⁺. Not accounting for this would thus also affect the [CII]/FIR ratio, resulting in values which could be too high by up to a factor of ~2.5 (see Figures 6 and A1 in Ebagezio et al. 2023).

Hence, we speculate that we identify a general evolution of the [CII]/FIR ratio as a function of the age of the associated HII region:

1. At early stages (which we marked with the green squares in Fig. 5), the reduction of the [CII]/FIR ratio occurs on small scales, driven by a significant enhancement of I_FIR in regions with newly formed stars.

2. At later stages (black rectangle in Fig. 5), and thus in more evolved HII regions, the reduction of the [CII]/FIR ratio involves larger scales and is caused mainly by the further photo-ionisation of C⁺.

In order to investigate this more in detail, in Section 3.4 we systematically analyse the I_[CII]/I_FIR ratio in all HII regions that we have in our simulated clouds.

Fig. 5 Emission maps of I[CII] and IFIR, and maps of I[CII]/IFIR, for MC2-MHD, tevol = 5 Myr (top), and MC2-HD, tevol = 4 Myr (bottom). In the very young HII regions highlighted with a green rectangle, the low [CII]/FIR ratio is mostly due to the high IFIR. Conversely, the black rectangle highlights a more developed HII region, where the [CII] deficit is mainly caused by the C+ → C2+ ionisation due to stellar radiation.

Fig. 6 Fraction of C+ (with respect to C+ and C2+) for MC2-HD at different times projected along the y-direction. As time evolves, C+ gets progressively converted into C2+ inside the HI regions due to further photo-ionisation by the radiation of the newly formed stars. This leads to a drop in the [CII] intensity, in particular in the region enclosed with the black rectangle (compare with Fig. 5).

3.3 Projection effects

Before we further investigate the [CII]/FIR ratio, we first explore the importance of projections effects, which can be substantial due to the turbulent (sub) structure of MCs. In general, MCs are often observed to have filamentary (Bally et al. 1987; André et al. 2014; Rezaei Kh. et al. 2020; Hacar et al. 2023) and/or sheetlike substructures (Zucker et al. 2021; Rezaei Kh. & Kainulainen 2022). This is verified for the simulations analysed in this paper, too (Ganguly et al. 2023, 2024). For this reason, the aspect ratio of a cloud is strongly dependent on the chosen LOS along which it is observed, which also affects the (apparent) shape of HII regions. For instance, a star-forming region in a sheet-like structure exhibits a ring shape when observed face-on, with the formed stars approximately in the centre. If the same HII region is observed edge-on, it may even be not observable at all due to the attenuation by the gas in the sheet (sitting in the foreground).

The impact of projection effects is evident in Fig. 7, where we plot the I_[CII] emission maps of MC2-HD along the three principle axes. With green and light cyan circles, we highlight the stars responsible for the creation of two HII regions, “bubble 1” and “bubble 2”. Indeed, the shape of the two bubbles is significantly different along the three LOS. In particular, bubble 2 cannot be seen at all in the x − y plane, but is clearly evident in the x − ɀ and y − ɀ planes. The bubble has a diameter of ~40 pc and is devoid of [CII] emission inside, due to the large amount of radiation coming from several hot stars inside (T_star ≲ 5 × 10⁴ K), which further photo-ionises C⁺ to C²⁺. Bubble 1 (associated with the stars in the green circle) exhibits different shapes when observed from different LOS as well, even though there is no LOS where the bubble is not visible at all. Nevertheless, the apparent size of the bubble changes significantly for different LOS.

Projection effects can also cause apparent bubbles, similar to those associated with expanding HII regions, but without any stars inside. Such an apparent bubble could be ‘created’ by a dense foreground or background layer that separates the part of the real HII region containing the driving stars from a more distant region of the same region. A possible example of such an apparent bubble is shown in Fig. 8, a zoom-in of the right panel of Fig. 7. In the centre of the figure there is a region with very low I_[CII] surrounded by a much brighter rim, which in combination resembles the structure of actual HII regions. However, there are no stars embedded in this structure. Hence, some caution is advised when analysing bubble structures in order to properly account for projections effects. We note that apparent bubbles with no stellar counterparts are also found in real observations of MCs and HII regions: for instance, Beuther et al. (2022) find several of such bubbles in [CII] emission maps of NGC 7538.

3.4 The [CII]-FIR ratio within HII regions

In this section we study I_[CII], I_FIR, and I_[CII]/I_FIR as a function of the distance from the centre of the HII region (see our definition below), in order to analyse radial trends of the [CII]-FIR ratio and the possible dependence on the age of the HII regions and on the kind of stellar population which creates them.

In Fig. 9 we show an example of the procedure we follow. The top row shows, from left to right, the [CII] emission map of MC2-HD in the x − ɀ plane at t_evol = 2, 3, and 4 Myr, respectively. By visual inspection of all snapshots for all times projected from the x-, y- and ɀ-direction, we identify for each simulation various HII regions (or bubbles), which we list in Table 3. We also identify their associated stars; for example in Fig. 9, we highlight with a green circle the stars which are responsible for the HII region of interest (namely bubble 2 of MC2-HD, see Table 3).

Next, in the middle row of Fig. 9 we show, from left to right, the azimuthal average of I_[CII], I_FIR, and I_[CII]/I_FIR, taken in radial bins at a distance d, as a function of the distance from the centre of the bubble. We here define the centre of the bubble as the position of the “main” star of the HII region. It is defined as the hottest star in the oldest sink particle associated with the HII region (see Table 3). This definition is motivated by the consideration that the size of the bubble is mainly determined by the time available for the expansion (this is the reason why we choose the oldest sink) and that the hottest star produces most of the radiation. We point out, however, that such a ‘main’ star may not necessarily be the star which constitutes the main source of ionising radiation for all cells in the HII regions. This can happen when a bright star forms later on in a different sink particle. Hence, our definition of a main star is mainly used to define a geometric centre of the associated HII region for further investigation.

The bottom row shows the same quantities as the middle row, but for each distance d we now average over all the pixels, which have a distance D smaller than the given distance d: $$I(D<d)=\frac{1}{N}{\displaystyle \sum _{i=1}^N{I}_i}.$$(2)

Solid, dashed, and dotted lines represent the calculation performed for the x − y, x − ɀ, and y − ɀ plane, respectively.

By looking at both, the emission maps and I_[CII](d) (Fig. 9, top and middle row), we note that at an early evolutionary stage I_[CII] decreases with increasing d. When going to later evolutionary stages, I_[CII] is low at small d (up to a few pc), then rapidly increases and eventually slightly decreases again at large d. This is a direct consequence of the expansion of the HII regions with time, as well as of the further ionisation of C⁺ inside such regions (lowering I_[CII]) and its high excitation in the rims (enhancing I_[CII]). The distance within which such a lack of [CII] is observed, depends not only on the age of the HII region but also on the LOS considered (see e.g. the various blue lines for t_evol = 4 Myr), which can be attributed to the aforementioned projection effects (Section 3.3).

The FIR intensity (central panel of the middle row) is very high at early evolutionary stages, reaching ~0.1 erg s⁻¹ cm⁻² sr⁻¹ in the proximity of the formed stars and then rapidly decreases with increasing d due to the drop in both density and dust temperature as we go further away from the stars. At later stages, the dust density is less peaked around the star and the dust is also warmer at large distances due to the presence of other stars: this results in a less pronounced (but still existing) decreasing trend of I_FIR with increasing d. The additional peaks in I_FIR (those at d = 5–7 pc) are due to other neighbouring stars.

The features of I_[CII] and I_FIR together lead to an I_[CII]/I_FIR ratio which progressively increases with d (right panel in the middle row of Fig. 9). The increase is sharp at early stages, with I_[CII]/I_FIR reaching a roughy constant value of ≃10⁻² at d ≃ 5 pc. At later stages, there is a central region where I_[CII]/I_FIR ranges from ~10⁻⁶ to ≲10⁻³, depending on the snapshot and the LOS, followed by an increase up to values of ≃10⁻³. This is due to the aforementioned lack of [CII] emission inside the HII region and due to the enhanced I_FIR emission in the part of the cloud surrounding it.

The profiles of I_[CII], I_FIR, and I_[CII]/I_FIR averaged over areas with D ≤ d (bottom row of Fig. 9) follow approximately the same trend as those averaged in radial bins. This is due to the fact that, in the averaging process, the outer bins contribute with more points than the inner bins as they contain more pixels. However, high values of I in the inner bins can still bring I(D < d) to values above I(d) at high d (this is the case, for instance, for I_FIR at t_evol = 2 Myr).

The behaviour of I_[CII], I_FIR, and I_[CII]/I_FIR described above refers to a specific HII region within MC2-HD. However, the same physical phenomena, namely the dust heating, the clearing of the HII region due to ionising radiation, and the C⁺ excitation and ionisation, produce similar trends in the other HII regions simulated here as well. However, differences in the stellar population and the surrounding environment can produce quantitative differences between those HII regions. In Fig. 10 we show I_[CII]/I_FIR as a function of d for all HII regions identified in our simulations. For each bubble, we identify the main star of the HII region with the procedure described above. We list the main star’s properties for each bubble and time in Table 3. We note that the stellar temperature, T_★, evolves according to the model described in Gatto et al. (2017), Peters et al. (2017) and Haid et al. (2019). Comparing this table with the I_[CII]/I_FIR plots in Fig. 10, we identify the following results:

• At small distances (d < 2 pc), I_[CII]/I_FIR drops to very small values as anticipated in Fig. 9. The lowest values reached are, however, independent of the main star temperature: for instance, for MC1-MHD at t_evol = 4 Myr, the main stars in the bubbles 2 and 3 have very similar temperatures, but the lowest values of I_[CII]/I_FIR are 10⁻⁴ and 10⁻⁶, respectively.

• At intermediate distances (2 < d < 5–10 pc) the [CII]/FIR ratio quickly increases to values around ~10⁻² to ~10⁻³. A [CII]/FIR ratio < 10⁻³ extending up to radii of 5–10 pc is only clearly visible in the bubbles 1 and 2 of MC2-HD and, to a smaller extent, in bubble 1 of MC1-HD. These three bubbles are characterised by an extremely hot star (45 000–50 000 K) and generally a large number of stars inside them. This causes an extended HII region with a lack of [CII] emission, which is the reason of the low I_[CII]/I_FIR ratio (Section 3.2). However, the geometry and the LOS from which we observe have a great impact on whether this feature is observable or not (see Section 3.3).

• At large distances (d > 5–10 pc), I_[CII]/I_FIR(d) is rather flat. Its value, however, decreases within approximately the first Myr of the HII region’s lifetime from ~ 10⁻² to ~ 10⁻³. This is a consequence not only of the ionising radiation of the main star, but also of the subsequent birth of numerous other stars in the HII region itself, as well as in other parts of the cloud (see top row of Fig. 9). The progress in star formation thus leads to an increased dust heating, increasing the FIR emission (central panel of Fig. 9) and thus lowering the I_[CII]/I_FIR ratio.

• Despite the qualitative trends of I_[CII]/I_FIR with distance and time, we were not able to identify clear quantitative trends. In general, we do not find a relation of the absolute value of I_[CII]/I_FIR with neither the stellar age t_★, its temperature or luminosity, nor the local column density evaluated within a radius of 2 pc (a proxy for the density in the surroundings of the star). We also did not find a relation of I_[CII]/I_FIR with a combination of two out of the four aforementioned quantities. We speculate that the main reason for this is the fact that the I_[CII]/I_FIR ratio is strongly influenced by a nonlinear interaction of all four quantities, projection effects as well as the exact density profile, which is not necessarily spherically symmetric.

Motivated by the analysis of Orion A presented by Pabst et al. (2021), in Fig. 11 we show the distribution of I_[CII] and I_FIR for all pixels of one snapshot, colour-coded according to the distance to the main star. In all three panels we can identify a “main branch”, characterised by a high I_[CII] for a given I_FIR. At early stages (t_evol = 2 Myr), the pixels at small distance constitute a portion of this branch. However, as time evolves, the pixels at small distance are not on the “main branch” any more: they seem to form a separate branch which is located below the main one. This appearance of a second branch corresponds to the drop of the I_[CII]/I_FIR ratio over time discussed before. The reasons for this are, as stated before, the stellar ionising radiation and the subsequent second ionisation of carbon.

The ionising radiation disperses the cloud: therefore, the density of both dust and C⁺ drop at small d. In combination with the second ionisation of C⁺ into C²⁺ (see Fig. 6), this leads to the observed drop of I_[CII]. For the FIR emission, the reduction of the density is partially compensated by the increase in dust temperature (see e.g. Fig. 4), resulting in a smaller decrease in I_FIR. Overall, this creates the separated branch, characterised by a lower I_[CII] for a given I_FIR. For other bubbles and/or LOS we reproduce these two branches only partially. This is mostly due to projection effects (see also Section 3.3).

In Fig. 11 we also overplot data for M42 in Orion A and data for Orion B obtained by Pabst et al. (2021, see their Figure 5³). The overall trend is similar showing an increase of I_[CII] with I_FIR and larger values being found preferentially at smaller distances⁴, although the observational data tend to lie slightly above our data points. This might be attributed to the potentially higher radiation field strength and density in this very compact region, which would lead to elevated intensities. For example, the luminosity of the binary system θ¹ Ori C in Orion A sums up to more than 3 × 10⁵ L_⊙ (GRAVITY Collaboration 2018) compared to less than 2 × 10⁵ L_⊙ in the region simulated by us. In addition, projection effects and the generally highly complex density structures could cause deviations. Hence, given these potential differences, we consider the agreement between the simulated and observed data as remarkable. In particular the slope for the high-intensity regions is similar, specifically for our early time at 2 Myr.

Finally, given that we see a slight decrease of I_[CII] at a given I_FIR over time in our simulations, this could mean that observationally found relations lying at lower I_[CII] belong to further evolved HII regions. However, given that (i) the spread in our simulated data is already rather large and (ii) the aforementioned differences in the properties of the MC/HII region could potentially cause intrinsic variations, we consider the quantitative predictive power of such a time evolution as rather limited.

Fig. 7 [CII] emission map of MC2-HD at tevol = 4 Myr for three different LOSs. The green and cyan circles highlight the stars (see Fig. 4 for a legend of the star symbols) responsible for the formation of ‘bubble 1’ and ‘bubble 2’, respectively. Because of the complex structure of the cloud, the shape of such HII regions is strongly dependent on the chosen LOS. As an example, “bubble 2” is not visible at all when looking along the z-axis (left panel). This explains the variations of the [CII] deficit as a function of the bubble radius for different LOS (Fig. 10).

Fig. 8 Zoom-in on an apparent cavity in the [CII] emission from the full [CII] emission map of MC2-HD at tevol = 4 Myr shown in the right panel of Fig. 7. This apparent cavity does not contain any stars (see Fig. 4 for a legend of the star symbols) and could be a potential consequence of projection effects, which can create apparent bubbles similar to expanding HII regions.

Fig. 9 [CII] intensity maps and profiles for MC2-HD at different times. Top row: maps of I[CII] for MC2-HD (x − ɀ plane) at tevol = 2–4 Myr (from left to right). Stars are represented using the same symbols as in Fig. 4. The green circles show the stars which are responsible for creating the HII region of interest, i.e. bubble 2 (see Table 3), inspected in the middle and bottom row. Middle row: average I[CII], IFIR, and I[CII]/IFIR at a given distance d from the main star within the HII region highlighted by the green circles in the upper row. Solid, dashed, and dotted lines represent the values for the x − y, x − ɀ, and y − ɀ plane, respectively. Bottom row: same as in middle row, but now for an average over all points with a distance D < d.

Fig. 10 I[CII]/IFIR as a function of the distance d from the main star of all HII regions found in the investigated clouds. At small d, I[CII]/IFIR drops to 10−6 − 10−4, depending on the regions and snapshots. At larger radii, the ratio is ~10−2 at early stages, and ~10−3 at later stages. The drop of the ratio at later stages is mostly due to the overall dust heating in the entire cloud caused by the various stars which have already formed.

Fig. 11 Distribution of I[CII] as a function of IFIR for all pixels within 10 pc from the main star of MC2-HD, bubble 2, in the x − ɀ plane. The colour bar indicates the distance from the main star. To guide the readers eye, we show with the black line a constant ratio of I[CII]/IFIR = 10−2. For pixels at close-by distances, the values of I[CII] and IFIR decrease over time due to the evacuation of the bubble, whereas at larger distances the values remain in a similar range. This leads to a separation into a “main branch”, formed by distant pixels, and a branch of close-by pixels, located below the main one. The squares show observational data from Orion A & B taken from Pabst et al. (2021, binned for the sake of readability) for comparison.

3.5 The [CII] deficit on molecular clouds scales

3.5.1 Time evolution of the [CII] deficit

Next, we investigate the [CII] and FIR luminosities (and their ratios) integrated over the entire cloud scale (≲100 pc) to mimic unresolved MCs in observations and study how they evolve over time. In Fig. 12 we show the evolution of L_[CII], L_FIR, and L_[CII]/L_FIR with time, using t_sf as a measure, hence referring to the point when the first massive star forms (see Table 2). We plot the results for all three LOS; note that differences among the different LOS are rather minor and hard to discern in the plot. For MC1-HD and MC2-HD, we add snapshots for the ɀ-direction, which are separated by 0.1 Myr only, allowing us a more detailed time analysis.

The [CII] luminosity generally increases with time, as already pointed out in Ebagezio et al. (2023). This is consistent with the common usage of L_[CII] as a star formation tracer (e.g. Herrera-Camus et al. 2015; Sutter et al. 2019; Bigiel et al. 2020; Pabst et al. 2021; Liang et al. 2024). For MC1-HD, the increase in L_[CII] slows down after t_sf ≃ 1.5 Myr. For MC2-HD, L_[CII] reaches a plateau at t_sf - 1 Myr, then slowly decreases up to t_sf ≃ 2 Myr. This general behaviour matches well with the evolution of the mass of all massive stars in the simulations (compare with Fig. 2). In addition, we see a second increase of L_[CII] for MC2-HD after t_sf ≃ 2 Myr. We attribute it to the stellar radiation permeating now larger portions of the cloud due to its fractal structure including large, low-density areas (compare figure 20 in Seifried et al. 2020) and partially overlapping HII regions, leading to a larger amount of exited C⁺. In addition, the stellar mass evolution in Fig. 2 also indicates a new star formation event in MC2-HD at late times. For the two MHD clouds, the increase in L[CII] is somewhat slower, which we attribute to the lower star formation rate in the beginning (Fig. 2).

Concerning the evolution of L_FIR, we first note that for the two HD clouds, L_FIR increases rapidly, then reaches a maximum and finally starts decreasing again. The maximum occurs at t_sf ≃ 2 Myr for MC1-HD and at t_sf ≃ 1 Myr for MC2-HD. The initial increase is caused by the bright FIR emission stemming from the area in the vicinity of the newly formed massive stars (indicated by the areas in the green rectangles in Fig. 5). The further evolution is indicated in the same figure, considering the bubble contained in the black rectangle. It shows a more evolved HII region, where there is more widespread FIR emission, but at a much lower intensity than at earlier stages.

We investigate the detailed evolution leading to this behaviour in Fig. 13, where we show the LOS-averaged, mass-weighted dust temperature 〈T_d〉, the dust column density N_d, and I_FIR for MC2-HD, bubble 2 (see Table 3) at t_evol = 2, 3, and 4 Myr. Over time, there is a progressively larger region with contains relatively warm dust, namely with 〈T_d〉 ≳ 20 K. However, the warm region becomes less dense (middle row), as the stellar feedback pushes the dust towards the edge, forming a bubble. In the shell of the bubble, 〈T_d〉 is measurably lower than inside.

Hence, there are two phenomena which counteract each other: (i) the enlargement of the warm regions leads to an increase of L_FIR, and (ii) the redistribution of dust into the cooler shell leads to a decrease. The net result of these two processes is shown in the bottom row: there is a first phase when L_FIR increases, as the enlargement of the hot areas overcomes the decrease of N_d , and a second phase when LFIR decreases again, as the redistribution (and cooling) of dust now takes over.

This explanation so far refers to a single bubble. The evolution of L_FIR of an entire cloud can then be explained as the combined effect of several bubbles evolving in the manner described above. Specifically, for MC2-HD in Fig. 12, we can see the phenomenon (i) followed by phenomenon (ii), which in turn is followed by phenomenon (i) attributed to a second phase of star formation (Fig. 2).

Finally, we consider the L_[CII]/L_FIR ratio. Due to the stronger increase of L_FIR compared to that of L[CII], the ratio decreases initially with t_sf, down to a roughly constant value of ~10⁻³ with a scatter of a factor of a few. In MC2-HD, at t_sf ≳ 1.5 Myr, L_[CII]/L_FIR increases again due to the aforementioned sharp increase of L_[CII]. Overall, we can thus explain the origin of the [CII] deficit by the different qualitative evolution of the [CII] and FIR luminosity.

We note that the apparent slower increase of L_[CII] and L_FIR for the MHD clouds with respect to the HD clouds is due to the delayed evolution that the MHD clouds experience: magnetic fields slow down the star formation process as they act against gravitational collapse leading to a lower star formation rate (Girichidis et al. 2018; Seifried et al. 2020; Ebagezio et al. 2023; Ganguly et al. 2023).

Fig. 12 Evolution of L[CII], LFIR, and L[CII]/LFIR (from left to right) as a function of tsf. The four clouds are represented with different colours and the three LOS with different symbols; the highly time-resolved datapoints along the x − y-projection for MC1-HD and MC2-HD are connected with a line. Overall, we observe a general tendency of a decreasing L[CII]/LFIR with increasing tsf. There is, however, some scatter in the values for different clouds at a given tsf, which we attribute to different star formation activities. We note the logarithmic y-axis scale in the right panel.

Fig. 13 Mass-weighted dust temperature (〈Td〉, top row), dust column density (Nd, middle row), and FIR intensity (IFIR, bottom row) along the y-axis of a portion of MC2-HD centred around bubble 2 (see Table 3) at tevol = 2–4 Myr (from left to right). The expansion of the bubble progressively enlarges the region with warm dust, characterised by 〈Td〉 ≳ 20 K. However, Nd in the warm region is reduced by the expansion of the bubble, which pushes the dust into the much cooler shell of the bubble. As a consequence, IFIR arises from a progressively larger, but less FIR-bright area. This is also reflected in the total LFIR, reported for each snapshot in the IFIR maps.

Fig. 14 Same as Fig. 12, but as a function of $x_{{\text{H}}^{+}}$. We note a general increase of L[CII] with increasing $x_{{\text{H}}^{+}}$ and a maximum in LFIR at $x_{{\text{H}}^{+}}\simeq {10}^{-1}$ in the two HD clouds. As a consequence, L[CII]/LFIR decreases for increasing $x_{{\text{H}}^{+}}$ up to this point, and then increases again.

Fig. 15 Same as Fig. 12, but as a function of L★,tot. We fitted a power-law (Eq. (4), green line) for the relation of L[CII]/LFIR with L★,tot. The relation holds as long as LFIR increases (middle panel). After reaching its maximum, the decrease of LFIR, which has no counterpart in L[CII], breaks this relation, which we therefore do not consider in the fit.

3.5.2 Dependence on H+ and stellar luminosity

Due to the creation of H+ by stellar ionising radiation, we find a rough correlation between the mass of H⁺ and t_sf in our simulations. For this reason, we next investigate the dependence on L_[CII]/L_FIR on the total mass fraction of ionised hydrogen. This is defined as $$x_{{\text{H}}^{+}}=\frac{M_{{\text{H}}^{+}}}{M_{\text{H},\text{tot}}},$$(3)

where $M_{{\text{H}}^{+}}$ is the total mass of H⁺ in the cloud, and M_H,tot is the total hydrogen mass.

The dependence of L_[CII], L_FIR, and L_[CII]/L_FIR on $x_{{\text{H}}^{+}}$ is shown in Fig. 14. Symbols and colours have the same meaning as in Fig. 12. Indeed, for each cloud the relations of L_[CII], L_FIR, and L_[CII]/L_FIR with $x_{{\text{H}}^{+}}$ are qualitatively similar to those with t_sf. However, the differences among the clouds are smaller when $x_{{\text{H}}^{+}}$ instead of t_sf is used. This suggests that $x_{{\text{H}}^{+}}$ might be a better parameter to describe the evolution of L_[CII], L_FIR and L_[CII]/L_FIR than t_sf. However, the correlation of L_[CII] with $M_{{\text{H}}^{+}}$ is only of indirect nature: first, inside an HII region both the amount of C⁺ (influencing L_[CII]) and H⁺ are mainly set by a common factor, that is the amount of the stellar radiation. Second, in an HII region the excitation of C⁺ is mainly due to free electrons, whose abundance in turn is directly proportional to the abundance of H⁺, which is the main source of free electrons in HII regions. Also for L_FIR the coupling with $M_{{\text{H}}^{+}}$ occurs only indirectly via the amount of stellar radiation, which both sets $M_{{\text{H}}^{+}}$ and heats up the dust grains, thereby leading to an increase in L_FIR.

In this context we note that H_α emission could further help to disentangle the apparent correlation between $M_{{\text{H}}^{+}}$ and L[CII]. If L[CII] is dominated by emission from within the HII region, this might lead to a good spatial correlation between the H_α. and [CII] emission, whereas for the case that the [CII] emission comes mainly from the photodissociation region in the (neutral) shell, H_α and [CII] could be spatially disjunct. Overall, this might depend crucially on the evolutionary stage (Fig. 5) and/or the considered direction (Fig. 7), which is why we defer a thorough investigation of this question to a subsequent work.

Next, we investigate the relation of L_[CII], L_FIR, and L_[CII]/L_FIR with the stellar luminosity, L_★,tot, as the impact of stellar feedback is determined, to some extent, by the stellar luminosity. We define L_★,tot as the sum of the bolometric luminosity of all stars with a mass M_★ 8 M_⊙⁵. The relation of L_[CII] and L_FIR with L_★,tot (Fig. 15) is qualitatively similar to the relations with t_sf and $x_{{\text{H}}^{+}}$, showing an initial increase and, in the case of MC2-HD, a subsequent slight decrease and a last sharp increase. However, the values for different clouds are more scattered, for a given L_★,tot, than in the case of $x_{{\text{H}}^{+}}$.

Similar results were also found by Vallini et al. (2017), who modelled in a simplified manner the photoevaporation of MCs by stellar UV radiation and the resulting [CII] luminosity. Like here, the authors report a general increase of L_[CII] with G₀ (their parametrisation of the stellar luminosity). Also quantitatively our simulations and their models seem to reasonably match. Assuming a typical G₀ of ~100 caused by the stellar radiation (see Fig. A.1), the authors obtain [CII] luminosities of a few 10³ L_⊙ in good agreement with our findings. Furthermore, also the specific [CII] luminosities, L_[CII]/M, of 10^{−(1.5–1.8)} L_⊙/M_⊙ reported by the authors for these G₀ values roughly agree with our findings of L_[CII]/M ≃ 10^−(1.5–2). Overall, the comparability of the findings obtained via these two different modelling approaches thus supports the robustness of the results reported here.

Lastly, the relation between L_[CII]/L_FIR and L_★,tot shows a decrease of L_[CII]/L_FIR for increasing L_★,tot for all clouds. For MC2-HD, it is followed by an increase at L_★,tot ≳ 4 × 10⁶ L_⊙, which corresponds to the new star formation phase (see Fig. 2) and the corresponding additional C+ excitement already mentioned before. Remarkably, apart from the data points associated with this phase, the scatter of L_[CII]/L_FIR as a function of L_★,tot is minimal compared to the cases of $x_{{\text{H}}^{+}}$ and t_sf. This relation can be well modelled with a power law, as shown in the right panel of Fig. 15. We fit a power law to all snapshots of MC1-HD and MC2-HD apart from those at very late stages, that is, those after L_FIR reaches its maximum (t_sf = 1 Myr for MC2-HD and t_sf = 2 Myr for MC1-HD). For this, we use the logarithmised values of L_[CII] and L_★,tot and fit a linear relation with the POLYFIT routine from NUMPY using equal weights for each data point and a χ² minimisation. Translating this back into a power-law relation (i.e. not-logarithmic data) yields the following best fit (green line in Fig. 15): $$L_{[\text{CII]}}/L_{\text{FIR}}=4.15\times {\left(\frac{L_{\star ,\text{\hspace{0.17em}tot\hspace{0.17em}}}}{L_{\odot }}\right)}^{-0.56}.$$(4)

This fit provides a potential method to assess the total stellar luminosity of a cloud by using [CII] and FIR luminosity measurements. This relation holds up to when L_FIR increases. As we discussed in detail in Section 3.5.1, this corresponds to the phase of life of MCs, when stellar feedback heats the dust in the dense regions close to the star formation sites, but has not yet completely dispersed these regions and compressed the dust and gas into cool, dense shells surrounding the HII regions.

4 Discussion

4.1 Saturation of [CII] in HII regions

The origin of the [CII] deficit has been often attributed to the saturation of [CII] caused by stellar radiation, which prevents the [CII] emissivity from increasing above a certain amount, no matter how strong the radiation field becomes (Kaufman et al. 1999; Muñoz & Oh 2016; Rybak et al. 2019; Bisbas et al. 2022). Also Vallini et al. (2017) report in their idealised models of photoe-vaporating MCs a saturation of L_[CII] for radiation field strengths corresponding to G₀ ≳ 10³.

For the purpose of testing this hypothesis, we use RADMC-3D to calculate the level population of the C⁺ ions (after correction for possible ionisation into C²⁺) and analyse its dependence on the dust temperature, in order to link it to the FIR emission. More specifically, let n_u be the number density of C⁺ ions in the upper level, ²P_3/2 , and n_l the number density of those in the lower level,² P_1/2. We then investigate the fraction of ions in the upper level, n_u/(n_u + n_l) as a function of T_d. In Fig. 16 we show this relation for MC2-HD at t_evol = 4 Myr in form of a mass-weighted phase diagram. We overplot the mean and the median (solid and dashed red line, respectively). We do not consider cells with T_d < 18 K, as here the main heating source is the ISRF, in which we are not interested in here. Furthermore, as for T_d > 55 K, the number of cells is very low, which makes the mean and median very noisy, we also do not consider them here.

The fraction of ions in the excited state remains low (n_u/(n_u + n_l) ≃ 0.03) up to T_d ≃ 23 K, then slowly increases until T_d ≃ 40 K, where the distribution appears to flatten around ~0.6. We point out that the theoretical maximum of (n_u/n_l) is 2, because of the statistical weights of the upper and lower level, and therefore max(n_u/(n_u + n_l)) = 2/3. This value is only reached for very high gas temperatures and when the density is above the critical density for C+ of around 10³ cm⁻³. Hence, we do not expect (n_u/(n_u + n_l)) to exactly reach the maximum of 2/3 in our simulations. Nonetheless, the clear sign of flattening of (n_u/(n_u + n_l)) for T_d ≥ 40 K strongly indicates that the [CII] line indeed saturates in this regime. Since the regions with high T_d are those close to bright stars (Figs. 4 and 13, top row), this shows that in the vicinity of bright stars the [CII] line is indeed saturated.

Hence, overall the saturation of the [CII] line contributes to the rather low L_[CII]/L_FIR ratio, as most of the FIR emission comes from the vicinity of the stars – where [CII] is saturated. In addition, as stated before, also the conversion of C⁺ into C²⁺ contributes to the L_[CII]/L_FIR low ratio, lowering L_[CII] by up to a factor of ~2.5 compared to the case when no further ionisation would be considered (Ebagezio et al. 2023).

Fig. 16 Mass-weighted 2D histogram of the fraction of C+ ions in the 2P3/2 level (nu/(nu + nl)) as a function of the dust temperature Td. The solid red line represents the mass-weighted average of nu/(nu + nl) for a given Td bin, and the dashed line represents the median. The thin, dotted black line represents nu/(nu + nl) = 2/3, that is, the maximum possible value. We show the result for MC2-HD at tevol = 4 Myr, with 18 K ≤ Td ≤ 55 K. On average, C+ reaches rough saturation at Td ≃ 40 K. This confirms that, in vicinity of stars, the [CII] line is saturated, causing, in part, the [CII] deficit.

4.2 The [CII] deficit and its link to evolutionary stages of star formation

In the following we summarise our results discussed so far concerning the origin and evolution of the [CII] deficit. Particularly, we link it to the different stages of MCs and star formation. We suggest the following correlation between star formation and the [CII] deficit:

1. Phase 1: pristine MCs: before star formation sets in, we see widespread and diffuse [CII] and FIR emission with typical intensity ratios of ≳10⁻² (regions outside the coloured boxes in Fig. 5). This results in an overall L_[CII]/L_FIR ratio averaged over the entire clouds of ≳10⁻².

2. Phase 2: onset of star formation and young HII regions: directly after new stars were formed, the dust temperature and density in their surroundings are large, causing a high I_FIR in these regions. Simultaneously, the increase in I_[CII] is limited in this phase, due to the saturation of [CII] and the second ionisation of C⁺ into C²⁺. In combination, this causes the I_[CII]/I_FIR ratio to be low close to the newly formed stars (see green enclosed regions in Fig. 5), whereas in the surrounding regions the I_[CII]/I_FIR ratio still remains around 10⁻² . Over time, the regions affected by the radiative feedback expand leading to a global decline of the L_[CII]/L_FIR ratio in MCs to values around ~10⁻³.

3. Phase 3: evolved HII regions: in this phase, L_FIR remains constant or even decreases as the dust is pushed into cold shells while the HII regions are evacuated (region encapsulated with the black rectangle in Fig. 5). The elevated excitation of C⁺ by free electrons and higher temperatures created by the stellar radiation, which now permeates the entire cloud, compensates the decrease of I_[CII] inside the evacuated HII region. Overall, this leads to a phase with roughly constant values of L_[CII]/L_FIR around ~10⁻³.

4. Phase 4: secondary star formation: the potential onset of a secondary star formation phase can lead to a subsequent increase of both L_[FIR] and L_[CII], where the latter appears to be more pronounced. This results in a slight increase of L_[CII]/L_FIR as seen for MC2-HD (see orange dots in Fig. 15).

We emphasise that throughout phase 2 to 4, the L_[CII]/L_FIR ratio is lower than in a cloud not affected by stellar radiation (phase 1). From Fig. 12 it can be seen that this stage is reached within short timescales of ≲1 Myr, which is typically significantly shorter than the duration of the star forming phase of clouds, during which they get disrupted. As shown in Haid et al. (2019), for the simulations underlying this publication, the ionising feedback disrupts the clouds within a few Myr (see also Walch et al. 2012). These timescales for cloud disruption are in line with a number of other numerical as well observational results, which are discussed in detail in the review of Chevance et al. (2023, see their Chapters 4.2 and 5 for a summary and an extensive list of literature). Hence, during most of their star forming phase, MCs will exhibit a rather low value of I_[CII]/I_FIR.

Hence, the (relatively seen) long timespan, at which star forming clouds exhibit a low I_[CII]/I_FIR ratio around ~ 10⁻³, could thus explain the observed [CII] deficit (Malhotra et al. 1997; Luhman et al. 1998; Malhotra et al. 2001; Luhman et al. 2003; Stacey et al. 2010; Casey et al. 2014; Smith et al. 2017; Hu et al. 2019).

4.3 Comparison with [CII] deficit observations

We next compare our results with actual observations. Starting on scales of (and below) individual clouds, a well-studied MC in the solar neighbourhood is the Orion Molecular Cloud 1 (OMC 1). From the Sun, it is located beyond the Orion Nebula cluster at a distance of ≃414 pc (Genzel & Stutzki 1989). It is an actively star forming cloud, similar to the ones we simulate here, despite being more massive. A study of the [CII] deficit in OMC 1 has been performed by Goicoechea et al. (2015), which used Herschel HIFI data to obtain a I_[CII]/I_FIR map. The authors measure the radial change of I_[CII]/I_FIR with distance to their reference centre (the HII region around the Trapezium cluster). On the smallest scale of 0.1 pc, the authors obtain a value of 9.3 × 10⁻⁴, which increases to 3.8 × 10⁻³ when averaged over the entire cloud (see their Table 2). Similar values are found by Jackson et al. (2020) in their analysis of 4 star-forming clumps in the Milky Way (scales of ~ 1 pc), reporting on average L_[CII]/L_FIR ≃ 10⁻⁴−10⁻³. Generally, these observational values agree well with our results shown in Fig. 10, despite a certain underlying scatter. As discussed in Section 3.4, this scatter in the both observations and our simulations can in parts be attributed to different stellar populations, projection effects and global characteristics of the clouds. Furthermore, our findings of an increasing I_[CII]/I_FIR with larger portions of clouds considered is in good agreement with the results for OMC 1 by Goicoechea et al. (2015).

Also the global evolution of L_[CII]/L_FIR reported in our work (Fig. 12, scales of ≲ 100 pc) is in good agreement with observations. Suzuki et al. (2021) find values ranging between 10⁻³ and 10⁻² in the star-forming region RCW 36. In particular, the authors observe a decreasing L_[CII]/L_FIR ratio for increasing L_FIR, which is in excellent agreement with the evolution for L_[CII]/L_FIR found by us: in our simulations L_[CII]/L_FIR decreases with time, while L_FIR increases (Fig. 12), which corresponds to a decrease in L_[CII]/L_FIR with increasing L_FIR as well.

The relation between L_[CII] and L_FIR has also been analysed for extragalactic observations. For M51, Pineda et al. (2018) find L_[CII]/L_FIR ≃ 0.003 in the centre of the galaxy, and somewhat larger values in the inter-arm regions, generally matching the values we present here. Observations of samples of galaxies (Díaz-Santos et al. 2017; Smith et al. 2017; Herrera-Camus et al. 2018) show L_[CII]/L_FIR values between a few 10⁻⁴ and 10⁻² as well. Specifically, the observations also show a decrease of L_[CII]/L_FIR with increasing L_FIR. As noted before, this qualitative behaviour matches well our findings. Sutter et al. (2021) studied 18 Local Volume galaxies and found that HII regions play a significant role in the overall [CII] deficit in their observed sample, in agreement with our findings. Although the authors state that thermalisation, that is collisional de-excitation, of C+ plays a significant role in setting the L_[CII]/L_FIR ratio, they also argue that it cannot be the only effect. In the extragalactic context, the recent galaxy-scale hydrochemical models of Liang et al. (2024) find that the origin of the [CII] deficit can be associated with lower metallicities, an effect not studied here. However, the astrochemical models of Bisbas et al. (2024) using more realistic C/O ratios in low-metallicity environments, find only a weak correlation with the [CII] deficit, in agreement with the high-ɀ observations of Harikane et al. (2020).

Furthermore, numerical simulations (Smith et al. 2017; Bisbas et al. 2022) fit the relation between L_[CII]/L_FIR and the star formation rate surface density Σ_SFR by means of a power-law with an exponent of −1/4.7. In our work we refrained from investigating this correlation, as on cloud scales – contrary to the galactic scales considered by the aforementioned authors – the momentary SFR can vary significantly. However, we do find L_[CII]/L_FIR to be correlated with time (and thus increasing stellar mass) and in particular with the stellar luminosity (see Fig. 15). Specifically, we also find that the dependence of L_[CII]/L_FIR on L_★,tot can be well described by a power-law relation (Eq. (4)).

5 Conclusions

In this paper, we present an analysis of simulations of four star-forming MCs within the SILCC-Zoom project (Seifried et al. 2017; Haid et al. 2019), including stellar feedback and an on-the-fly chemical network. In particular, we analyse the [CII] 158 µm line and the FIR dust continuum emission. The aim of this paper is to investigate the origin and evolution of the so-called [CII] deficit, which is the observed drop in the L_[CII]/L_FIR ratio with increasing SFR and infrared luminosity.

We performed two analyses, considering (i) spatially resolved emission inside MCs and HII regions and (ii) averages over the entire cloud scale. First, concerning the resolved [CII] and FIR emission in and around individual HII regions, we report the following results:

• Due to the heating of dust by stellar radiation, I_FIR is initially high in the vicinity of newly born stars, and then moderately decreases over time in the centre of the HII regions as the HII regions get evacuated and the gas gets compressed into dense and cool shells.

• On the other hand, I_CII only experiences a moderate increase after the onset of star formation and then quickly drops over time within the HII region. In addition to the intrinsic dilution of the gas inside the HII regions, we identify the further photo-ionisation of C⁺ into C² + as a major contributor to this decrease, reducing the C⁺ content to a few percent of its original value.

• In combination, the two effects lead to a strong drop in I_[CII]/I_FIR in the vicinity of newly born stars, that is, in the HII regions surrounding them. Specifically, we find that I_[CII]/I_FIR decreases from values of 10⁻³−10⁻² at scales above 10 pc to values of around 10⁻⁶−10⁻⁴ at scales below 2 pc.

• Projection effects can significantly change the radial profile of I_[CII] and I_FIR, and their ratio. Depending on the geometry of the cloud and on the viewing angle, an HII region may not be visible at all. Conversely, foreground material can cause apparent bubbles that do not host any stars.

Considering the evolution of total luminosities and the origin of the [CII] deficit on larger scales (~ 100 pc, corresponding to extragalactic observations), we report the following findings:

• The L_[CII]/L_FIR ratio decreases from values of ≳10⁻² in MCs without star formation to values of around ~10⁻³ in MCs with star formation within ~1 Myr.

• We attribute this decrease and thus the origin of the [CII] deficit to two main contributors: (i) the saturation of the [CII] line and (ii) the conversion of C⁺ into C²⁺ by stellar radiation.

• The aforementioned drop in the L_[CII]/L_FIR ratio can be divided into two phases after the onset of star formation: during the early evolution of HII regions, the saturation of [CII] and the further ionisation of C⁺ limit the increase in L_[CII] as star formation progresses, while L_FIR increases rapidly, leading to an initial decline of L_[CII]/L_FIR.

• In more evolved HII regions, L_CII stagnates and even partially drops over time due to the aforementioned secondary ionisation of C⁺ and the saturation of [CII]. L_FIR also stagnates as the gas gets pushed into the cooler shells surrounding the HII region. In combination, this keeps the global L_[CII]/L_FIR ratio at low values of ~10⁻³.

Overall, our work shows that the evolution of the [CII] deficit is linked to the various evolutionary stages of MCs and their embedded star formation process. During this process different mechanisms contribute to varying degrees to the [CII] deficit. Our work can thus aid in interpreting observations of [CII] emission from scales of individual HII regions to galactic scales.

Acknowledgements

The authors thank the anonymous referee for the constructive report, which helped to increase the quality of the paper. The authors also thank C. Pabst for providing the observational data reused here in Fig. 11. S.E., D.S., and S.W. thank the Deutsche Forschungsgemeinschaft (DFG) for funding through the SFB 956 “Conditions and Impact of Star Formation” (projects C5 and C6) and SFB 1601 “Habitats of massive stars across cosmic time” (sub-projects B1, B4 and B6). Furthermore, the project is receiving funding from the programme “Profilbildung 2020”, an initiative of the Ministry of Culture and Science of the State of Northrhine Westphalia. T.G.B. acknowledges support from DFG Grant No. 424563772, and the Leading Innovation and Entrepreneur-ship Team of Zhejiang Province of China (Grant No. 2023R01008). The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. We particularly thank the Regional Computing Center Cologne of the University of Cologne for providing the computational facilities for this project by hosting our computing cluster “Odin”.

Appendix A The radiation field below 13.6 eV

As stated in Section 2.1, we do not include the effect of the stars on the radiation field at energies below 13.6 eV in the simulations. Doing so would have the same effect as locally increasing the G₀ factor for the ISRF. In an effort to get a rough estimation of how large this increase of G₀ would be, we present a here a first-order estimate. First, given that for each massive star we know its luminosity and surface temperature T_★, we calculate the total luminosity of all stars in a MC in the range between 6 and 13.6 eV, L_★,ISRF (corresponding to the photon energy range of the ISRF, Draine 1978). Next, we make the assumption that all stars in a given MC are located at the same position. With this assumption, we can now convert L_★,ISRF into an effective G_{0, eff}, which would be present at an assumed distance d from the stars (via geometrical attenuation ∝ d⁻²). Given the typical extent of the clouds and HII regions presented in this work, we choose d = 10, 20 and 30 pc. We note that in this first-order estimate we neglect attenuation by dust.

We plot the resulting G_{0, eff} for all four clouds considered in this work in Fig. A.1. As can be seen, the effective G_{0, eff} created by the stars increases over time from a few 10 to ~ 100. We again stress that these values present only a very rough estimate, the actual increase in G₀ might vary significantly depending on the chosen location in the cloud. Nonetheless, the substantial increase of G₀ with respect to our fiducial value of 1.7 could have some effect on the temperature profile of the HII regions. In particular in the surrounding shells, where the ionising radiation is already extincted – and thus also the Cloudy post-processing step described in Section 2.2 is not applied – we do not use the G_{0, eff} produced by the nearby stars. This could slightly alter the temperature as well as the CO and thus C⁺ abundance in the PDR region. However, in the volume-wise larger interior of the HII regions, the effect would presumable be only of minor importance as here the heating of the gas would mainly be due to the ionising radiation, which is taken into account in the simulations. Hence, we expect that taking into account radiation below 13.6 eV would have only a moderate effect in a spatially very confined region.

Fig. A.1 Effective G0, eff created by the massive stars in the different clouds (see text). The solid line represents G0, eff at an assumed distance of 20 pc, the shaded area the range for distances between 10 pc and 30 pc. Typical values are of the order of a few 10 - 100.

Appendix B Monte-Carlo resolution study

We perform a convergence study to determine a sufficient number of RADMC-3D photon packages for the calculation of T_d. To do so, we consider MC2-HD at t_evol = 4 Myr, as this is the snapshot which contains the largest number of stars. As RADMC-3D splits up the available photon packages to the sources (i.e. stars), the number of packages per star is the lowest for this snapshot. Hence, convergence for this snapshot should also imply convergence for all other snapshots. In Fig. B.1 we show the distribution of T_d (calculated by RADMC-3D) for different numbers of photon packages. We note that with a low number of photon packages (10⁵ – 10⁶) the distribution is clearly not yet converged, whereas the distributions for 10⁹ and 10¹⁰ photon packages are almost identical (they vary for less than 1 per cent). We thus choose to use 10¹⁰ photon packages in the Monte Carlo dust temperature calculation.

Fig. B.1 Dust temperature distribution calculated with RADMC-3D of MC2-HD at tevol = 4 Myr for different numbers of photon packages. Differences between the distribution for 109 and 1010 photon packages are less than 1 per cent for all dust temperature bins, and therefore we chose to use 1010 photon packages for all dust temperature calculations, as the computational costs for this highest value are still manageable.

References

All Tables

All Figures

	Fig. 1 Column density maps of MC2-HD from tevol = 2–4 Myr projected along the y-direction (from left to right) showing the effect of star formation (black stars) on the dispersal of the cloud.
In the text

	Fig. 2 Time evolution of the mass of all massive stars (M > 8 M⊙) formed in each MC. The clouds without magnetic fields form stars more rapidly.
In the text

	Fig. 3 Dust temperature PDFs for all four clouds at tevol = 3 (solid lines) and 4 Myr (dashed lines). There is a wide spread of temperatures over about 1 order of magnitude, which increases as star formation proceeds over time.
In the text

Fig. 4 Maps of the mass-weighted mean dust temperature, 〈Td〉, for four snapshots of the simulated clouds. We show, from top left to bottom right, MCs with progressively more advanced star formation. We overplot the stars present in the simulations, where the size and colour indicate the age and temperature of the stars, respectively. In general, regions not affected by stellar feedback have 〈Td〉 ≃ 15 K, and regions where stellar feedback is important have 〈Td〉 ≃ 50 K. These regions are small (a few pc) and limited to the direct vicinity of the stars in the early stages of star formation (e.g. MC2-MHD, top left), but involve larger parts of the cloud (a few 10 pc) at later stages (e.g. MC2-HD, bottom right).

In the text

Fig. 5 Emission maps of I[CII] and IFIR, and maps of I[CII]/IFIR, for MC2-MHD, tevol = 5 Myr (top), and MC2-HD, tevol = 4 Myr (bottom). In the very young HII regions highlighted with a green rectangle, the low [CII]/FIR ratio is mostly due to the high IFIR. Conversely, the black rectangle highlights a more developed HII region, where the [CII] deficit is mainly caused by the C+ → C2+ ionisation due to stellar radiation.

In the text

	Fig. 6 Fraction of C+ (with respect to C+ and C2+) for MC2-HD at different times projected along the y-direction. As time evolves, C+ gets progressively converted into C2+ inside the HI regions due to further photo-ionisation by the radiation of the newly formed stars. This leads to a drop in the [CII] intensity, in particular in the region enclosed with the black rectangle (compare with Fig. 5).
In the text

Fig. 7 [CII] emission map of MC2-HD at tevol = 4 Myr for three different LOSs. The green and cyan circles highlight the stars (see Fig. 4 for a legend of the star symbols) responsible for the formation of ‘bubble 1’ and ‘bubble 2’, respectively. Because of the complex structure of the cloud, the shape of such HII regions is strongly dependent on the chosen LOS. As an example, “bubble 2” is not visible at all when looking along the z-axis (left panel). This explains the variations of the [CII] deficit as a function of the bubble radius for different LOS (Fig. 10).

In the text

	Fig. 8 Zoom-in on an apparent cavity in the [CII] emission from the full [CII] emission map of MC2-HD at tevol = 4 Myr shown in the right panel of Fig. 7. This apparent cavity does not contain any stars (see Fig. 4 for a legend of the star symbols) and could be a potential consequence of projection effects, which can create apparent bubbles similar to expanding HII regions.
In the text

Fig. 9 [CII] intensity maps and profiles for MC2-HD at different times. Top row: maps of I[CII] for MC2-HD (x − ɀ plane) at tevol = 2–4 Myr (from left to right). Stars are represented using the same symbols as in Fig. 4. The green circles show the stars which are responsible for creating the HII region of interest, i.e. bubble 2 (see Table 3), inspected in the middle and bottom row. Middle row: average I[CII], IFIR, and I[CII]/IFIR at a given distance d from the main star within the HII region highlighted by the green circles in the upper row. Solid, dashed, and dotted lines represent the values for the x − y, x − ɀ, and y − ɀ plane, respectively. Bottom row: same as in middle row, but now for an average over all points with a distance D < d.

In the text

Fig. 10 I[CII]/IFIR as a function of the distance d from the main star of all HII regions found in the investigated clouds. At small d, I[CII]/IFIR drops to 10−6 − 10−4, depending on the regions and snapshots. At larger radii, the ratio is ~10−2 at early stages, and ~10−3 at later stages. The drop of the ratio at later stages is mostly due to the overall dust heating in the entire cloud caused by the various stars which have already formed.

In the text

Fig. 11 Distribution of I[CII] as a function of IFIR for all pixels within 10 pc from the main star of MC2-HD, bubble 2, in the x − ɀ plane. The colour bar indicates the distance from the main star. To guide the readers eye, we show with the black line a constant ratio of I[CII]/IFIR = 10−2. For pixels at close-by distances, the values of I[CII] and IFIR decrease over time due to the evacuation of the bubble, whereas at larger distances the values remain in a similar range. This leads to a separation into a “main branch”, formed by distant pixels, and a branch of close-by pixels, located below the main one. The squares show observational data from Orion A & B taken from Pabst et al. (2021, binned for the sake of readability) for comparison.

In the text

Fig. 12 Evolution of L[CII], LFIR, and L[CII]/LFIR (from left to right) as a function of tsf. The four clouds are represented with different colours and the three LOS with different symbols; the highly time-resolved datapoints along the x − y-projection for MC1-HD and MC2-HD are connected with a line. Overall, we observe a general tendency of a decreasing L[CII]/LFIR with increasing tsf. There is, however, some scatter in the values for different clouds at a given tsf, which we attribute to different star formation activities. We note the logarithmic y-axis scale in the right panel.

In the text

Fig. 13 Mass-weighted dust temperature (〈Td〉, top row), dust column density (Nd, middle row), and FIR intensity (IFIR, bottom row) along the y-axis of a portion of MC2-HD centred around bubble 2 (see Table 3) at tevol = 2–4 Myr (from left to right). The expansion of the bubble progressively enlarges the region with warm dust, characterised by 〈Td〉 ≳ 20 K. However, Nd in the warm region is reduced by the expansion of the bubble, which pushes the dust into the much cooler shell of the bubble. As a consequence, IFIR arises from a progressively larger, but less FIR-bright area. This is also reflected in the total LFIR, reported for each snapshot in the IFIR maps.

In the text

	Fig. 14 Same as Fig. 12, but as a function of $x_{{\text{H}}^{+}}$. We note a general increase of L[CII] with increasing $x_{{\text{H}}^{+}}$ and a maximum in LFIR at $x_{{\text{H}}^{+}}\simeq {10}^{-1}$ in the two HD clouds. As a consequence, L[CII]/LFIR decreases for increasing $x_{{\text{H}}^{+}}$ up to this point, and then increases again.
In the text

	Fig. 15 Same as Fig. 12, but as a function of L★,tot. We fitted a power-law (Eq. (4), green line) for the relation of L[CII]/LFIR with L★,tot. The relation holds as long as LFIR increases (middle panel). After reaching its maximum, the decrease of LFIR, which has no counterpart in L[CII], breaks this relation, which we therefore do not consider in the fit.
In the text

Fig. 16 Mass-weighted 2D histogram of the fraction of C+ ions in the 2P3/2 level (nu/(nu + nl)) as a function of the dust temperature Td. The solid red line represents the mass-weighted average of nu/(nu + nl) for a given Td bin, and the dashed line represents the median. The thin, dotted black line represents nu/(nu + nl) = 2/3, that is, the maximum possible value. We show the result for MC2-HD at tevol = 4 Myr, with 18 K ≤ Td ≤ 55 K. On average, C+ reaches rough saturation at Td ≃ 40 K. This confirms that, in vicinity of stars, the [CII] line is saturated, causing, in part, the [CII] deficit.

In the text

	Fig. A.1 Effective G0, eff created by the massive stars in the different clouds (see text). The solid line represents G0, eff at an assumed distance of 20 pc, the shaded area the range for distances between 10 pc and 30 pc. Typical values are of the order of a few 10 - 100.
In the text

Fig. B.1 Dust temperature distribution calculated with RADMC-3D of MC2-HD at tevol = 4 Myr for different numbers of photon packages. Differences between the distribution for 109 and 1010 photon packages are less than 1 per cent for all dust temperature bins, and therefore we chose to use 1010 photon packages for all dust temperature calculations, as the computational costs for this highest value are still manageable.

In the text