The origin of the [CII]-deficit in a simulated dwarf galaxies starburst

We present [CII] synthetic observations of smoothed particle hydrodynamics (SPH) simulations of a dwarf galaxy merger. The merging process varies the star-formation rate by more than three orders of magnitude. Several star clusters are formed, the feedback of which disperses and unbinds the dense gas through expanding HII regions and supernova (SN) explosions. For galaxies with properties similar to the modelled ones, we find that the [CII] emission remains optically thin throughout the merging process. We identify the Warm Neutral Medium ($3<\log T_{\rm gas}<4$ with $\chi_{\rm HI}>2\chi_{\rm H2}$) to be the primary source of [CII] emission ($\sim58\%$ contribution), although at stages when the HII regions are young and dense (during star cluster formation or SNe in the form of ionized bubbles) they can contribute $\gtrsim50\%$ to the total [CII] emission. We find that the [CII]/FIR ratio decreases due to thermal saturation of the [CII] emission caused by strong FUV radiation fields emitted by the massive star clusters, leading to a [CII]-deficit medium. We investigate the [CII]-SFR relation and find an approximately linear correlation which agrees well with observations, particularly those from the Dwarf Galaxy Survey. Our simulation reproduces the observed trends of [CII]/FIR versus $\Sigma_{\rm SFR}$ and $\Sigma_{\rm FIR}$, and it agrees well with the Kennicutt relation of SFR-FIR luminosity. We propose that local peaks of [CII] in resolved observations may provide evidence for ongoing massive cluster formation.


INTRODUCTION
The star-formation rate (SFR) is of fundamental importance for understanding the cyclic process of global star formation in galaxies across the epochs (see Madau & Dickinson 2014, for a review). Measuring it can reveal the properties and the evolutionary stages of the observed interstellar medium (ISM). Schmidt (1959) and Kennicutt (1998) were the first to find a strong correlation between the SFR per unit area and the gas surface density, a relation frequently referred to as "Schmidt-Kennicutt relation". Since then, various methods based on continuum bands and optical/nearinfrared (IR) emission lines have been used to measure SFR in different systems (see Kennicutt & Evans 2012, for a review). Recent attempts using fine-structure lines such as [Oi] at 63 µm and [Oiii] at 88 µm (e.g. Hunter et al. 2001;Brauher et al. 2008;De Looze et al. 2014;Olsen et al. 2017), as well as Hα, UV and IR (e.g. Shivaei et al. 2015) have shown good correlations with the SFR. The far-IR fine-structure transition of [Cii] 2 P 3/2 − 2 P 1/2 at a rest frame wavelength of 157.7 µm (hereafter referred to simply as [Cii]) is also widely used as a promising diagnostic of the SFR (see Stacey et al. 1991;Boselli et al. 2002, for early attempts). The Atacama Millimeter Array (ALMA) is able to provide unprecedented resolution of high-redshift (z > 1) observations in this line (using Bands 5 − 10, depending on the redshift), opening an entirely new window in the study of the Early Universe ISM.
[Cii] is one of the brightest lines originating from starforming galaxies (e.g. Stacey et al. 1991;Brauher et al. 2008). The ionization potential of atomic carbon is 11.3 eV, slightly lower than the ionization potential of hydrogen (13.6 eV). Under typical ISM environmental conditions, the [Cii] emission line is a result of the interaction between the ISM gas and FUV photons. In general, the emission of [Cii] represents ∼ 10 −3 of the total far-infrared (FIR) continuum emission. Furthermore, it is also found to be associated with the outer shells of H 2 -rich clouds, where star-formation takes place. Thus, [Cii] plays a very important role in photodissociation regions (PDRs) as a coolant, particularly at low visual extinctions (Hollenbach & Tielens 1999;Wolfire et al. 2003). It has an upper-state energy of hν/k B ∼ 91 K and its critical density spans approximately three orders of magnitude depending on the temperature (Goldsmith et al. 2012).
[Cii] may, therefore, arise from different phases of the ISM, such as Hii-regions, PDRs and cold molecular gas (Velusamy & Langer 2014;Abdullah et al. 2017;Croxall et al. 2017;Accurso et al. 2017;Lagache et al. 2018;Ferrara et al. 2019;Cormier et al. 2019) depending on the environmental parameters, such as the intensity of the FUV radiation, metallicity, the cosmic-ray ionization rate (Bisbas et al. 2015b(Bisbas et al. , 2019(Bisbas et al. , 2021 and the intensity of X-rays (Mackey et al. 2019). Interestingly, the studies of Velusamy & Langer (2014) and Accurso et al. (2017) found that ∼ 60 − 85% of the total [Cii] emission arises from the molecular gas phase, thus naturally explaining its correlation with SFR (see also Madden et al. 2020). However, numerical simulations of Franeck et al. (2018) show that if the cloud is young enough, its emission in [Cii] arising from the molecular phase may be smaller than 20 per cent. In this regard, this fine-structure line may not be a good tracer for the CO-dark 1 H 2 -gas.
The [Cii]/FIR luminosity ratio is observed to decrease with increasing infrared luminosity (Malhotra et al. 1997(Malhotra et al. , 2001Luhman et al. 1998Luhman et al. , 2003Casey et al. 2014). The origin of the so called '[Cii]-deficit' is still be-1 The term 'CO-dark' gas has been introduced by van Dishoeck (1992) ing investigated despite numerous efforts proposing a variety of mechanisms behind it (e.g. Malhotra et al. 2001;Luhman et al. 2003;Stacey et al. 2010;Graciá-Carpio et al. 2011;Sargsyan et al. 2012). Suggestions include optically thick [Cii] emission in large columns of dust, conversion of singly (C + ) 2 to doubly (C 2+ ) ionized carbon, and fine-structure lines (e.g. [Oi]) overcoming [Cii] as coolants. Muñoz & Oh (2016) studied how strong FUV radiation fields can drive a low ratio of [Cii]/FIR due to thermal saturation of the [Cii] emission (see also Kaufman et al. 1999). This mechanism has been recently confirmed observationally by Rybak et al. (2019) and, as we will see later, it is also in accordance with our simulations for which we find a decreasing [Cii]/FIR as the surface densities of FIR and SFR increase. Narayanan & Krumholz (2017) provided a theoretical model suggesting that the cloud structure in galaxies with increasing SFRs, and hence increasing gas surface densities, is responsible for the [Cii]/FIR ratio decrease. Using a large sample of ∼ 15, 000 resolved regions, Smith et al. (2017) was able to show that even extragalactic regions of a few hundred of parsecs in size appear to be [Cii]-deficient. However, the effect is more prominent in the high-redshift Universe where distant and, thus, [Cii]-faint sources may emit only ∼ 10% of the expected [Cii] based on their observed FIR luminosity.
In a series of hydrodynamical simulations with a resolution of 4 M per gas particle, Hu et al. (2016Hu et al. ( , 2017Hu et al. ( , 2019 examined the global star formation process and how supernovae (SNe) affect the SFR in dwarf galaxies, as well as the underlying ISM microphysics including heating/cooling mechanisms and dust sputtering. Follow-up work by the griffin 3 collaboration (Lahén et al. 2020;hereafer 'L20') proposed that dwarf galaxy mergers may result in a significant population of star clusters. In particular, during the merging process the SFR may increase up to three orders of magnitude and can form clusters in the range of 10 2 − 10 6 M . Such a large variation of SFRs in the dynamical evolution provides an interesting set of three-dimensional morphological ISM distributions, including feedback, which can in turn reveal insights on the origin of the [Cii]-SFR correlation and the [Cii]-deficit.
The focus of this work is to perform [Cii] synthetic observations (see Haworth et al. 2018, for a review) of the griffin dwarf galaxy merger simulations presented 2 The [Cii] notation refers to the emission of the line whereas the C + notation to the actual species and/or its abundance.  (De Looze et al. 2011Pineda et al. 2014;Herrera-Camus et al. 2015, 2018Zanella et al. 2018;Sutter et al. 2019) and numerical (Olsen et al. 2015(Olsen et al. , 2017Vallini et al. 2015;Lupi & Bovino 2020) studies, find a [Cii]-SFR relation of the form: where 0.8 m 1.2 and −8 b −5 depending on the type and redshift of the galaxy. When correlating the star-formation rate surface density (Σ SFR ) with the surface [Cii] luminosity (Σ CII ), the above relation takes the form: (2) Table 1 provides a summary of the m, b, m Σ and b Σ values used in this work. This paper is organized as follows. Section 2 gives a description of the selected snapshots from the griffin SPH simulations and the post-processing strategy. Section 3 discusses how the Cii luminosity and the SFR vary in time throughout the merging process. Section 4 studies the origin of the [Cii] emission and Section 5 its relation with the FIR emission and how the feedback from massive clusters leads to a [Cii]-deficient medium. Finally, in Section 6 we compare our results with observations and examine the relation between SFR and FIR as well as the [Cii] and SFR. We conclude in Section 7.

DESCRIPTION OF SIMULATIONS AND THE POST-PROCESSING TECHNIQUE
The hydrodynamical simulations of the dwarf galaxy merger are fully described in L20. These are Smoothed Particle Hydrodynamics simulations using the SPHGal code presented in Hu et al. (2014Hu et al. ( , 2016Hu et al. ( , 2017, which is a modified version of gadget-3 (Springel 2005) adopting a modern formulation of SPH that includes timedependent artificial diffusion (viscosity and conduction) to improve SPH fluid-mixing behaviour. The calculations include a non-equilibrium model for cooling and chemistry that directly integrates the rate equations of H 2 , H + and CO, while obtaining non-equilibrium abundances for H, C + , O, and free electrons from residual conservation laws (Nelson & Langer 1997;Glover & Mac Low 2007). They also take into account metal line cooling from 12 different metal species (H, He, N, C, O, Si, Mg, Fe, S, Ca, Ne and Zn) based on the cooling tables of Wiersma et al. (2009) as implemented in Aumer et al. (2013). For reference and for the Z = 0.1 Z environmental condition, the initial abundances of C + and O relative to hydrogen were 2.46 × 10 −5 and 4.9 × 10 −5 , respectively (see L20 for details).
In addition, routines treating the interstellar radiation field, as well as stellar feedback in terms of photoionization, photoelectric heating and SNe are included. Ionization due to collisions is also treated as described in Hu et al. (2017). The treatment of Hii regions has been described in Hu et al. (2017) and it is a Strömgren-type approximation for photoionization as discussed in Hopkins et al. (2012). This approach was chosen so as to reduce the computational cost as opposed to a more detailed radiative transfer approach. We note that the adopted chemistry model does not account for the conversion of Cii to Ciii which may somewhat overestimate the calculated Cii luminosity. However, as we will see later, the contribution of Hii regions to the total Cii luminosity is small and thus a more detailed approach shall not alter the results presented in this work. The dynamical expansion of Hii regions has been also benchmarked against the results of the STARBENCH workshop (Bisbas et al. 2015a) with reasonable agreement.
With 4 M per SPH particle, the simulation resolves the Sedov-Taylor stage of individual SN-remnants for > 90% of the ambient SN-densities (Hu et al. 2017(Hu et al. , 2019Steinwandel et al. 2020). Prior to the collision, the two dwarf galaxies are identical with virial masses of M vir = 2 × 10 10 M , virial radii of r vir = 44 kpc and are composed of a dark matter halo and a gas-rich disk with a rotationally supported (old) stellar population. The two galaxies are set on parabolic orbits with a peri-centric distance of 1.46 kpc and an initial separation of 5 kpc. The collision is not edge-on but includes an inclination similar to the Antennae galaxies merger, as described in Lahén et al. (2019).
Under the optically thin assumption (see §3), the total Cii luminosity is given by the expression where Λ CII is the Cii cooling function of the i−th SPH particle 4 (in units of erg s −1 cm −3 ), m SPH is its corresponding mass, m p is the proton mass and the summation is over all particles within the volume of interest. We construct luminosity maps at a resolution of 1024 2 uniform pixels, where we project the SPH particles. The luminosity of each pixel is then a direct summation of the SPH particles along the line-of-sight using the above equation.
We calculate the total far-infrared (FIR) emissivity by adopting the dust cooling rate which is given by the expression: where ρ is the gas density, B ν (T d ) is the Planck function at frequency ν, T d is the dust temperature and κ ν is the dust opacity. Glover & Clark (2012) fit this relation using the expression (see also Hu et al. 2017): where D is the dust-to-gas mass ratio relative to the solar value 5 . Here, we set D = 0.1 since the modelled dwarf galaxies have a metallicity of Z = 0.1 Z . For the purposes of this work this linear relation between dust-to-gas ratio and metallicity is generally a good assumption, although metal-poor systems with metallicities lower than the one examined here may not follow such a relation (Herrera-Camus et al. 2012;Rémy-Ruyer et al. 2014). Equation 5 is valid for 5< T d < 100 K. The ∝ T 6 d dependency arises from the Stefan-Boltzmann law and the opacity term. The total FIR luminosity as well as the FIR luminosity maps, are then constructed as described above for the Cii luminosity.
In case we consider velocity-resolved emission properties and unless stated otherwise, we impose a lower observational limit for the [Cii] luminosity of L CII = 0.5 L , corresponding to a velocity integrated emission of W(CII) 0.6 K km/s (see Eqn. A1 for a resolution of 1024 2 and Appendix E for the effect of using a different lower limit). This is a reasonable assumption considering the sensitivity of instruments (see also Franeck et al. 2018). The FIR luminosity considered in this analysis, corresponds to the pixels that satisfy the aforementioned Cii observational criterion. In a similar way, the surface, Σ, is estimated by the area covered from the above number of pixels.
There are two merging events; a first passage encounter occurring at t ∼ 80 Myr and a second encounter leading to a final coalescence occurring at t ∼ 170 Myr. We post-process snapshots from t = 10 Myr to t = 390 Myr with a 10 Myr step. We therefore post-process a total of 39 snapshots. The SFR is calculated from the simulation snapshots following the methodology of Hu et al. (2016). This methodology is a stochastic star formation approach where the local SFR is sf ρ/t ff , where ρ is the gas density, t ff is the free-fall time, and sf = 0.02 is the star formation efficiency. The SFR quantity used in the present work is an average over the past 1 Myr.
The top row of Fig. 1 shows the total gas column density, N tot , the middle row the corresponding [Cii] luminosity and the bottom row the FIR luminosity in four different snapshots. These are (from left to right) during the first encounter at t = 70 Myr, during the second encounter at t = 160 and t = 170 Myr, and after the gas settling in the central part at t = 280 Myr. During the first encounter, the bar-like structure becomes bright in [Cii] only in its central part where the gas surface density becomes high enough N tot 10 21 cm −2 . During the second encounter, clusters are formed in the dense parts. They produce ionizing radiation creating Hiiregions. These, in turn, are responsible for the bubblelike features seen at t = 170 Myr. The positions of the three most massive clusters formed are indicated in these snapshots. The masses of the clusters are 1.6, 1.2, and 7.9 × 10 5 M for the first, second and third most massive cluster, respectively (L20). In the final phase shown in Fig. 1, the two disks merge and enhance feedback from Hii-regions. Subsequent SNe disperse and unbind the gas, creating the irregular shape seen at t = 280 Myr.

TIME EVOLUTION OF IONIZED CARBON LUMINOSITY AND STAR FORMATION RATE
The time evolution of the SFR, L CII and L CII /SFR for all gas mass is shown in Fig. 2. The top left panel shows the SFR calculated for the entire computational domain (thin light blue line) as well as for the inner 1 kpc (thick dark blue line). The latter SFR is the one we will use throughout this work. The two SFRs are in excellent agreement when the dwarf galaxies experience an encounter. We note that for the calculation of SFR in this panel, we use the outputs of the simulation in time intervals of ∆t = 1 Myr.
As described in L20 (see their Fig. 2), the first pericentric passage occurs at ∼ 50 Myr and the first apocenter at ∼ 80 Myr. During that period of time, a tidal bridge forms connecting the two galaxies which results in an approximately two orders of magnitude increase of the SFR, from ∼ 10 −4 to ∼ 10 −2 M yr −1 . The second and much stronger encounter occurs between ∼ 150 to ∼ 180 Myr with the SFR peaking at ∼ 160 Myr. This is the starburst phase where multiple clumpy star formation regions exist. During this second period, the SFR reaches values as high as 0.2 − 0.3 M yr −1 , corresponding to a mini-starburst. Earlier works by Hu et al. (2016Hu et al. ( , 2017 showed that in an isolated dwarf galaxy with properties similar to those modelled here, the SFR is approximately 10 −4 −10 −3 M yr −1 and relatively constant. This in turn means that throughout the merging process of such dwarf galaxies, the SFR can vary between two and three orders of magnitude depending on the evolutionary stage. It is interesting to explore whether or not a potential assumption of optically thin [Cii] emission is valid. We do this since [ 13 Cii] emission of the Large Magellanic Cloud studied by Okada et al. (2019) shows that [Cii] may become optically thick, having implications for the extragalactic observations of this line and, thus, the obtained [Cii]-SFR relation. For this investigation, we perform additional calculations with the radiative transfer code radmc-3d (see Appendix A). In the top right panel of Fig. 2, we plot with black lines the L CII calculated with radmc-3d versus time for three different  viewing angles (x−y, x−z, y−z planes). On top of these three lines, we plot with red solid line the corresponding radmc-3d calculations in the optically thin limit. As can be seen, all aforementioned lines are indistinguishable showing that in these simulations [Cii] can be very well approximated as optically thin. Furthermore, we plot with dashed magenta line the L CII calculated using Eqn. 3 which is in excellent agreement with the radmc-3d results. Finally, in this panel we also plot with a thin solid purple line the L CII value given from Eqn. 1 for m = 0.80 and b = − 5.73 corresponding to the DGS of De Looze et al. (2014). As we describe in §6.1, our simulations show an excellent agreement with these observations which best represent our modelled galaxies.
The modelled galaxies are low-mass with low metallicity so any opacity effects are negligible. We actually verified this by performing radmc-3d calculations. However, Franeck et al. (2018) showed that during molecular cloud formation at solar metallicity, [Cii] can become quickly optically thick. On the other hand, Bisbas et al. (2021) showed that metal-poor clouds may remain optically thin for H 2 column densities up to ∼10 23 cm −2 while the presence of strong FUV intensities can positively contribute to the increment of the [Cii] optical depth. Such conditions are exceptional for the modelled dwarf galaxies, thus [Cii] may remain always optically thin. We argue that the [Cii] emission of systems with similar properties to the simulated galaxies is in general optically thin. However, we cannot unambiguously demonstrate that larger systems and especially galaxies with metallicities close to solar will remain [Cii] optically thin.
The first encounter results in an increase of L CII spanning approximately two orders of magnitude, reaching ∼ 2×10 4 L at ∼ 70 Myr. When the dwarf galaxies reach the apocenter at ∼ 80 Myr, the column density decreases resulting in a decrease in SFR and L CII . The second, stronger, encounter results in a much more prominent increase in L CII , reaching values ∼ 5 × 10 5 L . As described in L20, SN-feedback from the clusters disperses the dense distribution of gas. This leads to a decrease of L CII at ∼ 180 Myr. However, for times > 200 Myr, L CII fluctuates following the trend of SFR due to the settling of the gas in the central region and the associated feedback.
The bottom left panel of Fig. 2 shows how the L CII /SFR ratio evolves in time. In this panel, we have considered an average SFR value over the preceding 5 Myr, every 10 Myr. As can be seen, this ratio remains approximately constant at a value of ∼ 3 − 6 × 10 6 L /M yr −1 for t > 200 Myr. Overall, the ratio does not strongly vary when compared to the fluctuations of both L CII and SFR that span more than three orders of magnitude throughout the evolution. The L CII /SFR ratio decreases during the first and especially during the second encounter. This relatively small fluctuation of this ratio compared to the corresponding one observed in both SFR and L CII individually, indicates that L CII is a good tracer for estimating the SFR.
The bottom right panel of Fig. 2 shows the total gas mass versus time at an interval of ∆t = 10 Myr. After t > 50 Myr, the average gas mass is ∼ 2×10 7 M , which is ∼ 1/4 of the total gas mass of the simulation setup described in L20. Figure 3 shows the time evolution of the [Cii] velocity integrated emission for the regions within which the three most massive clusters form. The plotted velocity integrated emission (here produced with radmc-3d; see Appendix A) is an average over 25 pixels which are centered around each cluster position shown in Fig. 1. The linear size of each of these pixels is approximately 7.8 pc, thus covering an area of ∼ 61 pc 2 . The 25 pixel-sized region, therefore, corresponds to an area of ∼ 1.5 × 10 3 pc 2 . We explore the behaviour of the [Cii] emission for three different viewing angles and find that the trends remain unchanged. We note that the differences observed in each viewing angle are not arising from optical depth effects but rather due to the different projections of mass covered by the aforementioned area. We find that the optically thin emission holds for each different line-of-sight in these areas, as well. The most prominent feature is observed for Cluster-3. As can be seen in the second and third panel of Fig. 1, this cluster is formed during the second encounter (at t∼ 170 Myr) and creates an Hii region, which eventually removes the ISM that satisfies the [Cii]-bright observational criterion (see §2). Thus, the [Cii] emission of that region decreases, reflecting the trend shown in the top panel of Fig. 3. This indicates that local peaks of W CII in resolved observations may provide evidence for ongoing massive cluster formation.

ORIGIN OF THE [CII] EMISSION
The interesting question about the origin of the [Cii] emission has been explored by various groups both numerically and observationally. Here, we analyse the simulation outputs and we study the contribution of [Cii] emission arising from the different ISM phases to the total emission. Each ISM phase (ionized, atomic, molecular) is identified according to the relative abundances (χ) of H + , H and H 2 . In particular, the photoionized ISM (Hii regions in which the energy of photons exceed the 13.6eV ionization potential of hydrogen) has [Cii] originates mainly from the WNM component at a percentage of ∼58%. CNM contributes an approximately constant ∼18% at all times. A similar contribution (∼14%) arises from the ionized material. Hii regions contribute in general a low percentage to the total emission (∼10%), although there are certain times where they are in an early dense evolutionary stage, thus dense, in which their contribution dominates all phases. Finally, the [Cii] emission originating from molecular gas is always negligible. a fixed χ H + = 0.9998 and a gas temperature in the range 10 4 < T gas < 1.3 × 10 4 K. The ionized ISM (resulting from both photoionization and collisional ionization) is defined as the gas with T gas > 10 4 K minus the aforementioned Hii contribution. The atomic ISM is defined as χ HI >2χ H2 and the molecular ISM is defined as χ HI <2χ H2 . We additionally divide the atomic medium in Warm Neutral Medium (WNM; 3 < log T gas < 4), and Cold Neutral Medium (CNM; log T gas < 3) (e.g. Wolfire et al. 2003). Each of the L CII of the aforementioned four ISM phases is then compared to the total Cii luminosity. Figure 4 shows this contribution throughout the duration of the simulation. The shaded region marks the duration of the second encounter. The emission of [Cii] originating from WNM dominates over the corresponding emission of all other ISM phases. In particular, WNM contributes an average of ∼58% in agreement with previous works (e.g. Hu et al. 2017), whereas the contribution of the CNM is ∼18%. The emission of the ionized gas (photoionized and collisionally ionized combined) has an average contribution of ∼24%. Throughout the simulation, the gas that is collisionally ionized remains as the main contributor of the [Cii] emission at this phase. Interestingly, the emission originating from Hii regions varies substantially throughout the duration of the simulation, showing that it depends strongly on its evolutionary stage. On average, the contribution remains quite low (∼5%). However, there are certain times e.g. at t = 160, 250, and 320 Myr where the [Cii] emission from Hii regions dominate over all different ISM phases, with a contribution as high as ∼50 − 60%. At early times (e.g. t < 200 Myr), this sudden increase in L CII is due to the newly formed Hii regions which, in their early evolutionary stages, are dense and very bright in [Cii]. Notably, such a high contribution has been recently observed in ionized regions of the inner Galaxy (Langer et al. 2021). At later times after the main encounter (e.g. t > 200 Myr), Hii regions are mainly formed due to supernova explosions which create bubbles of ionized material. The contribution of L CII originating from the molecular gas is negligible (∼0.02%) at all times. Figure 5 shows phase-plots (2D histograms) of gas temperature versus total number density, weighted with [Cii] luminosity. We consider six snapshots at the times of t = 40, 70, 120, 160, 250 and 390 Myr. Comparing these panels with Fig. 4, it can be seen that the upper bright part of the phase-plots (log n H ∼ 0 − 2, log T gas >3), which is the WNM component of the ISM, plays the most dominant role in the origin of [Cii]. It is interesting to note that at t = 70 Myr the bright curved rim of the WNM (starting at log n H ∼ 0, log T gas ∼ 4.0 with a declining trend as log n H increases) is a result of strong cooling in relatively dense and warm regions with log T gas <4, which is the temperature of the ionized gas in an Hii-region. Such strong cooling is associated with PDRs located ahead of the ionization front of the newly formed Hii regions. This bright rim is also seen at all times during and after the second encounter, thus making PDRs a considerable contributor to the origin of [Cii] emission.
Before the encounter at t = 40 Myr (upper left panel of Fig. 5), the density of WNM component is mainly in the range of −1 < log n H < 1. As the simulation progresses, the density of this ISM component increases and at the particular t = 170 Myr time (bottom left panel of Fig. 5), the above range extends up to log n H ∼ 4. Such densities are much higher than those found locally in the Milky Way (e.g. Wolfire et al. 1995). There are two main reasons that cause this; i) low metallicities shift the equilibrium curve to higher densities (Hu et al. 2016), and ii) the high FUV intensities due to feedback from cluster formation as well as supernova feedback, shift the equilibrium curve even further (Hu et al. 2017). In Figure 5. Phase-plots (Tgas versus total H-nucleus, nH, number density) weighted with LCII for snapshots at t = 40, 70, 120, 170, 250 and 390 Myr. Before the second encounter (top row), LCII originates from the WNM (see Fig. 4). Once the second encounter occurs (bottom row), Hii regions form and their ionized gas takes over as the main contributor to the total LCII. The horizontal straight lines at log Tgas = 4 is the gas temperature of the interior of Hii regions. Figure 6. Phase-plots of Tgas versus the number density of the three C + colliding partners for the t = 170 Myr snapshot (see also Fig. 5). Here, LCII is weighted with the abundance of each colliding partner. From left-to-right, each panel shows the number densities of atomic hydrogen (Hi), molecular hydrogen (H2) and electrons (e). The black solid lines correspond to the critical density of each partner versus Tgas (Goldsmith et al. 2012). We find that collisional de-excitation is negligible.
Appendix B, we additionally show mass-weighted phaseplots for the aforementioned snapshots. Figure 6 shows phase-plots of the t = 170 Myr snapshot for the three colliding partners (Hi, H 2 and e). The Cii luminosity is weighted with the corresponding relative abundance of each of the aforementioned colliding partners. The solid line in each panel shows the critical density of each partner as calculated by Goldsmith et al. (2012). Gas that falls in the right-hand part of each critical density relation is collisionally de-excited. We find that at all times, L CII associated with collisional de-excitation due to Hi and H 2 , is negligible. The same finding applies for electrons as collision partner, except for a ∼20 Myr period during the second encounter (t∼160 − 180 Myr; see Fig. 6) where the [Cii] luminosity arising from gas with n e >n crit,e is ∼30 − 50%. This is because in this short period, compact and dense Hii regions form, containing considerable amount of dense ionized gas. However, even during that period, collisional de-excitation still plays a minor role to the total L CII meaning that photoelectric heating is the dominant source of [Cii] emission at all times. Figure 7, shows 2D histograms of the Λ CII cooling function versus n H at t = 170 Myr for the ISM gas at the inner 1 kpc from Cluster-3 (see Fig. 1). At this time, the emission of [Cii] originates ∼65% from WNM, ∼15% from CNM, ∼20% from ionized gas, while Hii regions and molecular gas have negligible contributions ( 0.02%). As can be seen, the majority of Λ CII is associated with gas with n 10 3 cm −3 , which is approximately the critical density for collisions with Hi. For densities lower than the aforementioned, Λ CII scales as ∝ n 2 H while for higher ones it scales as ∝ n H (line is thermalized). It is therefore evident that collisional deexcitation plays a minor role. In Appendix C we present a theoretical approach as to how Λ CII builds as a function of n H and in Appendix D we show 2D histograms of Λ CII for the inner 0.1, 0.2 and 0.3 kpc regions from Cluster-3.
Breaking down the ISM components, it can be seen that main contributors (WNM, ionized and CNM) have a significant fraction of their Λ CII associated with n H 10 3 cm −3 gas. Similarly, the molecular component, although playing a small role, follows the same trend. It is worth mentioning that the Hii region gas is entirely thermalized meaning that the Cii emission that arises from this component, is a result of collisional de-excitations. This is in good agreement with the observational results of Sutter et al. (2021). Hii regions are the places where C + is ionized to form C 2+ . Given that the contribution of this ISM component is negligible, accounting for the transition between the two aforementioned ionization states of carbon can be excluded from our analysis.
In previous numerical studies, Accurso et al. (2017) found that ∼75% of the [Cii] emission in Milky Way as well as ∼60 − 80% in galaxies of the Herschel Reference Survey, arises from their molecular regions. In molecular cloud simulations, Franeck et al. (2018) found that [Cii] is primarily emitted from the cold (T gas ∼ 40− 65 K) neutral medium with densities n H ∼ 50−500 cm −3 . Yet, these simulations did not include star formation and stellar feedback. In isolated dwarf galaxy simulations, Lupi & Bovino (2020) identified the diffuse (n H 100 cm −3 ) neutral gas to contribute most of the [Cii] emission, while only a small fraction of [Cii] originates from higher density gas associated with dense PDRs. They do not, however, actually show the temperature of the [Cii] emitting gas, so the explanation that the CNM is the dominant component is only based on the density criterion. Interestingly, they do show that the typical densities of [Cii]-bright gas are higher for lower metallicity environments. Here, we do not base our ISM definition on density cuts as the full density and temperature evolution of the gas is available to us. In this way, we identify the WNM to be the dominant source of [Cii] emission. Following Wolfire et al. (2003), our definition of WNM is based on the gas temperature. Therefore the WNM can have higher or lower densities, depending on the local balance of heating and cooling terms. The phase-plots of Fig. 5 indicate that a T gas -based criterion is more appropriate in the case of vastly varying star formation rates, since there is warm but dense gas during the merging process, such that the WNM phase 6 shifts to higher gas densities. Our results agree with Hu et al. (2017) who also examined isolated dwarf galaxies and identified WNM to be the most [Cii]bright ISM phase.

RELATION BETWEEN THE [CII] LINE EMISSION AND THE FIR EMISSION
The left panel of Fig. 8 shows how the [Cii]/FIR ratio relates with Σ FIR , in which the observational criterion of L CII > 0.5 L has been applied. As can be seen, the ratio [Cii]/FIR decreases with increasing Σ FIR . For high L FIR , this makes the ISM gas to emit more brightly in FIR in relation to [Cii], which results in the known '[Cii]-deficit' (Malhotra et al. 1997(Malhotra et al. , 2001Luhman et al. 1998Luhman et al. , 2003Combes 2018). We plot (black squares) observations of 52 nearby galaxies (z< 0.2) from the shining 7 sample (Herrera-Camus et al. 2018). These are not dwarf galaxies and have higher metallicities. However, we find that the [Cii]/FIR ratio neither depends strongly on the metallicity nor does it depend strongly on the density distribution. Thus, this ratio can be compared to objects that may not necessarily satisfy the properties of a dwarf galaxy. On the other hand, Σ FIR strongly depends on metallicity (assuming a linear relation between dust-to-gas and metallicity). This in turn results in a shift of the SHINING galaxies to higher Σ FIR as can be seen in the left panel of Fig. 8. Had the modelled galaxies been at solar metallicity, it would have increased the derived Σ FIR by approximately one order of magnitude thereby matching with the lower end of the SHINING sample. Nevertheless, given that a decreasing [Cii]/FIR ratio with a comparable slope is observed in our simulations, it is interesting to explore and understand its origin.
While many different mechanisms leading to a [Cii]deficit medium have been proposed, we emphasize here 6 Note also that Fig. 13 of Hu et al. (2017) shows how sensitive the L CII cumulative functions are versus density and versus Tgas.
7 Survey with Herschel of the Interstellar medium in Nearby Infrared Galaxies on the effect of thermal saturation of [Cii] emission. The effect of thermal saturation of [Cii] leading to a [Cii]-deficit medium was suggested by Kaufman et al. (1999) and studied in detail theoretically by Muñoz & Oh (2016), with Rybak et al. (2019) to provide followup observational evidence for its existence in dusty starforming galaxies (with masses of ∼ 10 10 M ) at a redshift of z ∼ 3. As Muñoz & Oh (2016) explain, the thermal saturation of [Cii] is a direct quantum mechanical consequence of the saturation of the upper fine-structure energy state when the gas temperature exceeds the Cii excitation temperature of 91 K. Once the latter occurs, the population of the upper state cannot increase further leading to an approximately constant emissivity while the FIR dust emissivity is free to increase more. Their theoretical models lead to the expression: where f CII is the fraction of total gas traced by [Cii]. As described in Muñoz & Oh (2016), the value of f CII = 0.13 (which is also adopted in this work) is a good estimate based on observations of Milky Way clouds and various extragalactic sources. The above relation is shown in solid black line in the left panel of Fig. 8

and its power-law is in agreement with the Herrera-Camus et al. (2018) observations as well as our simulations.
In general and throughout the duration of the simulation, for SFR > 10 −2 M yr −1 it is found that the gas is so warm that [Cii] becomes thermally saturated. This increase in temperature is a direct consequence of the increase in FUV photoelectric heating as a result of the high star formation activity, which eventually leads to the birth of massive star clusters. High FUV intensities are to be expected in galaxy mergers. For instance, the PDR study of Bisbas et al. (2014) find an average of G 0 10 2.5 in the Antennae merging system. In our simulations, the consequence of high FUV intensities for high SFRs is demonstrated with the colour bar of Fig. 8 which shows that low [Cii]/FIR ratios are tightly connected with high values of star formation rate.
The right panel of Fig. 8 shows how the [Cii]/FIR ratio relates with Σ SFR . For the latter quantity, we use the same surface, Σ, obtained from the observed area of L CII > 0.5 L . The simulation points here are also colour-coded with SFR. As expected from the above discussion, the [Cii]/FIR ratio decreases with increasing Σ SFR . We compare our results with observations of the Kingfish 8 program presented in Smith et al. (2017).  These are spatially resolved observations of 54 nearby galaxies. Based on a galaxy sample with higher Σ SFR than examined here, Smith et al. (2017) find a bestfitting relation of the form: [CII]/FIR 10 −3 Σ SFR 12.7 −1/4.7 .
As with the [Cii]/FIR vs Σ FIR relation, our simulations have a similar slope with the observations and the above best-fit relation. We note that Eqn. 7 represents a single power-law fit to both local and high-redshift sources and that it can be applied when young stars provide the dominant energy source on scales greater than a few hundred parsecs . Thus, deviations of our simulations from this power-law are to be expected. Figure 9 shows a zoom-in of the central region at t = 170 Myr, where three massive clusters have been formed (see also Fig. 1). In particular, the [Cii] emis-sion, FIR emission and their ratio are shown. Here, we only highlight Cluster-3 which is responsible for the strong [Cii] and FIR emission in its immediate surroundings. In the right panel of Fig. 9 showing the [Cii]/FIR ratio, the dashed circles centered on Cluster-3 show radial distances with steps of 0.1 kpc. As can be seen, the innermost part has a very low [Cii]/FIR ratio, of the order of 10 −4 − 10 −3 and it is thus '[Cii]-deficit' when compared to the outer regions which have a [Cii]/FIR ratio of 10 −3 − 10 −2 . Such a [Cii]/FIR mapping has been observed in the central region of Orion Molecular Cloud 1 by Goicoechea et al. (2015) as well as in the wider Orion Nebula complex, recently, by Pabst et al. (2021).
The above correlation can be also seen in Fig. 10 in which the density-weighted dust and gas temperature as well as the luminosities of [Cii] and FIR are plotted versus the radial distance from Cluster-3. The aforementioned quantities are averaged over shells of 0.05 kpc thickness. For R < 0.1 kpc, the dust temperature is high with T d > 50 K which is a consequence of the very strong FUV radiation field emitted from the massive star cluster. This results in a high FIR luminosity (see Eqn. 5), with values up to ∼ 1.7 × 10 9 L in the R < 0.1 kpc. Similarly, the high gas temperatures of T gas > 3 × 10 3 K in that region results in a high luminosity of [Cii] ∼ 2.7 × 10 5 L . This leads to a [Cii]/FIR ratio of ∼ 1.5 × 10 −4 and thus a [Cii]-deficit gas. In the outer regions e.g. in the shell of 0.4 < R < 0.5 kpc, the dust temperature is ∼ 40 K which reduces the FIR luminosity an order of magnitude i.e. ∼ 1.5 × 10 8 L . On the other hand, the gas temperature, although it is also reduced, remains much higher than the 91 K excitation temperature of [Cii] i.e. ∼ 500 K. This makes the [Cii] luminosity to decrease by a factor of ∼ 3, thus leading to a higher [Cii]/FIR ratio. As shown in Appendix D, the ISM gas immediately around the cluster is thermalized and thus grows with ∝ n H . This growth cannot compensate with the ∝ T 6 dust correlation of Eqn. 5, leading to a [Cii]-deficit medium.
In these hydrodynamical simulations, DGR is constant in space and time. However, dust could be destroyed due to high FUV radiation fields or strong shocks (e.g. Draine & Salpeter 1979;Jones et al. 1994;Zhukovska et al. 2016). By performing the first hydrodynamical multiphase ISM simulations including dust sputtering due to SNe, Hu et al. (2019) showed that DGR can decrease by ∼30% in the volume filling warm gas compared to that in the dense clouds. We expect that such a decrease in DGR would locally result in a lower FIR emission in regions of very high FUV intensity. At the same time, the strength of the FUV In the third panel, the FIR luminosity is displaced downwards by a factor of 5 × 10 −4 to ease the comparison with [Cii] field is expected to be somehow more extended since the attenuation due to dust will be smaller and thus, G 0 will decrease primarily due to geometric dilution following a ∼r −2 law. Considering all the above, we expect [Cii]/FIR to locally decrease, which could result in a "less [Cii]-deficit" ISM gas, but the effect may be small compared to the [Cii]/FIR value obtained from the entire simulation.
We also note that while we do not include the conversion of Cii to Ciii in our chemical network as mentioned in §2, we do expect the derived [Cii]/FIR ratio Figure 11. PDR simulation of a nH = 300 cm −3 total number density interacting with various FUV intensities in the range G0 = 1 − 10 6 . Top left: emissivity of [Cii] versus the visual extinction, AV . As G0 increases, [Cii] emissivity increases from ∼ 6.5 × 10 −24 to ∼ 1.7 × 10 −23 erg s −1 cm −3 at which point it saturates, approaching asymptotically a maximum value. Bottom left: emissivity of FIR calculated assuming that the total dust cooling is equal to radiative dust heating (Eqn. 8), versus AV . Under the optically thin assumption, the FIR emission is given by integrating Eqn. (8) along the line-of-sight. Top right: gas temperature versus AV when thermal balance has been reached. The temperature at the surface of the PDR increases from ∼ 120 K to ∼ 1.6 × 10 3 K as G0 increases. Bottom right: [Cii]/FIR versus ΣFIR assuming optically thin emission. Due to the thermal saturation of the [Cii] emissivity, the ratio [Cii]/FIR decreases for G0 > 10 in these simulations, leading to a [Cii]-deficit medium. The gray triangles represent the simulation data as discussed in Fig. 8. to decrease if the higher ionization states of carbon were taken into account, thereby again enhancing the [Cii]deficit.

Photodissociation region calculations
To explore the decrease in the [Cii]/FIR ratio in greater detail, we perform PDR calculations using the publicly available 3d-pdr photodissociation region code 9 (Bisbas et al. 2012). The code uses the UMIST2012 database of reaction rates (McElroy et al. 2013) and performs iterations over thermal balance by taking into account various heating and cooling processes. It calculates the abundances of species, the gas temperatures as well as the emissivities of various 9 https://uclchem.github.io/3dpdr coolants using the Large Velocity Gradient approximation (Sobolev 1960;Castor 1970;de Jong et al. 1975). The dust temperature due to FUV heating is calculated using the treatment of Hollenbach et al. (1991) in which T d ∝ G 0.2 0 . In these PDR calculations, we explore the response of the [Cii] and FIR emissivities in a one-dimensional uniform density cloud with a total H-nucleus number density of n H = 300 cm −3 , as it interacts with various FUV intensities in the range G 0 = 1 − 10 6 , normalized to the spectral shape of Draine (1978). We use a subset of the UMIST2012 chemical network which contains 33 species (including e − ). For the purposes of this test, we also assume a cosmic-ray ionization rate of ζ CR = 3 × 10 −18 s −1 and metallicity of Z = 0.1 Z to imitate as closely as possible the adopted ISM con-ditions of L20. The cloud has a visual extinction of A V = 10 mag which is related to the total column density, N tot , as A V = A V,0 N tot (Z/Z ), where A V,0 = 6.3×10 −22 mag cm 2 (Weingartner & Draine 2001;Röllig et al. 2007). The size of the cloud is therefore taken to be L 170 pc. Figure 11 illustrates the results from the PDR simulations described above. The top-left panel shows how the [Cii] emissivity, which represents the [Cii] cooling rate, increases for increasing G 0 per cloud depth. For G 0 = 1, the emissivity at the surface of the cloud is ∼ 6.5×10 −24 erg s −1 cm −3 . As G 0 increases, the emissivity increases but becomes thermally saturated for G 0 > 10 5 at which point 10 it is ∼ 1.7 × 10 −23 erg s −1 cm −3 as seen in the top left panel of Fig. 11. Higher FUV intensities would increase the emissivity asymptotically to a maximum value, close enough to the aforementioned saturated value. On the other hand, for the assumed density of n H = 300 cm −3 , the local dust cooling, corresponding to the FIR emissivity, is approximately equal to the dust heating rate due to radiation. The latter is given by the expression (Glover & Clark 2012) where G 0 is the local (attenuated) FUV intensity. Therefore, the FIR emission is given by integrating the above expression along the line-of-sight, thus Λ dust = Γdr. As can be seen, the FIR emission scales linearly with the FUV intensity. The bottom-left panel shows how Λ dust relates with G 0 per cloud depth. High FUV intensities heat up the gas, as can be seen in the top-right panel. The thermal balance calculations performed, show that for G 0 = 1 the gas temperature at the surface of the cloud is T gas ∼ 120 K while for G 0 = 10 6 it is ∼ 1.6 × 10 3 K. Assuming optically thin emission for both [Cii] and FIR in this example, we integrate along the line of sight to obtain the corresponding integrated emission. This is shown in the bottom-right panel of Fig. 11, which correlates [Cii]/FIR with Σ FIR . As Σ FIR increases, [Cii]/FIR decreases leading to a [Cii]-deficit medium. Assuming a linear relation between dust-to-gas and metallicity, higher metallicities would drift the plotted curve in this 10 This value can be analytically calculated using the expression frequency, β ij = 1 the escape probability at the edge of the cloud, S ij the source function and B ij the black-body function for the 2.7 K background emission. For the simulation parameters with G 0 = 10 6 , 3d-pdr outputs n i ∼6 × 10 −4 cm −3 and n j ∼2 × 10 −3 cm −3 (j < i). By replacing these values and calculating S ij and B ij accordingly, we obtain Λ CII ∼1.7 × 10 −23 erg s −1 cm −3 . panel rightwards. Here, the simulation data are also shown with gray triangle. The PDR simulation and the simulation data are in excellent agreement. 6. DISCUSSION 6.1. The star-formation rate and far-infrared luminosity relation Kennicutt (1998) has calibrated SFR with FIR luminosity in dusty circumnuclear starbursts, providing the following relation: where C 1.87 × 10 −10 L accounts for the total IR luminosity covering the wavelength range of 3 − 110µm (Kennicutt & Evans 2012). In the above relation, L FIR corresponds to the total bolometric luminosity with the assumption that all FIR will emerge from dust grains heated by the interstellar FUV radiation field. In our simulations, dust heating is tightly connected with the increase of FUV radiation due to the formation of clusters and SN-feedback, so Eqn. 9 can be directly applied (c.f. Rieke et al. 2009, for applying this relation to observations). Figure 12 shows the SFR-FIR correlation for our simulations, colour coded with L CII luminosity. As expected, the luminosity of [Cii] increases with SFR and L FIR . The black solid line corresponds to Eqn. 9. We find that our simulations are in very good agreement with the Kennicutt (1998) calibration for a broad range of L FIR and SFR values, each one spanning approximately four orders of magnitude. Interestingly however, the agreement appears to break for lower values of SFR ( 3 × 10 −3 M yr −1 ). Such an effect was seen also in the recent simulations of Lahén et al. (2021). We speculate that when the UV radiation is low, there are not enough re-processed photons to produce the IR fluxes which would in turn provide reasonable estimates of the SFR predicted by Eqn. 9. In this regard, Lahén et al. (2021) further finds that the best agreement with the true SFR is reached with the 24µm corrected UV tracers.

The star-formation rate and [Cii] luminosity relation
We now compare the resultant [Cii]-SFR relation against observations found in the literature (see Appendix E for the corresponding Σ CII − Σ SFR relation). The comparison is illustrated in Fig. 13. The simulation points are colour-coded with the [Cii]/FIR ratio. As can be seen, the ratio decreases as both SFR and [Cii] increase, in accordance to the discussion in §5. In the Figure 12. Relation of SFR with LFIR, colour-coded with [Cii] luminosity. As expected, LCII increases as SFR and LFIR increase. The solid line corresponds to the Kennicutt (1998) relation (Eqn. 9). The agreement between the simulations and the latter relation is very good.
We plot the best-fitting relations from the following four observational works. From Herrera-Camus et al. (2015), who study a sample of 46 nearby star-forming galaxies from the Herschel Kingfish survey in the absence of strong Active Galactic Nuclei (AGN). From Pineda et al. (2014), who used the Herschel Galactic Observations of Terahertz C + (GOT C+) to study velocity resolved Milky Way clouds found in the Galactic plane. From Sutter et al. (2019), who studied nearby ( 30 Mpc) normal star-forming galaxies with no LIRGs included from Kingfish and BtP. We also plot the De Looze et al. (2014) relation of 42 dwarf galaxies from the DGS sample of Madden et al. (2013). Furthermore, we add the two L CII −SFR scalings discussed in Herrera-Camus et al. (2018) considering the star formation efficiency (SFE = L FIR /M mol , where M mol is the molecular mass). As described in Herrera-Camus et al. (2018), main-sequence, star-forming galaxies and AGNs have scalings similar to the normal SFE (blue dashed line) of the shining survey while LINERs and (U)LIRGs have scalings similar to the high SFE (orange dashed line). During the second encounter of the collision, our simulation has a better agreement with the high SFE slope thus mimicking -even for a short period of time-the average conditions found in more massive and starburst galaxies.  Table 1). In addition, our results compare well with the individual observations presented of DGS by Cormier et al. (2015Cormier et al. ( , 2019. Furthermore, Olsen et al. (2017) using cosmological zoom-in simulations, presented a [Cii]-SFR relation from 30 main-sequence galaxies at a redshift of z∼ 6. These galaxies are of low metallicity (Z = 0.1 − 0.4 Z ), matching our resolved dwarf galaxy simulations, albeit the Olsen et al. (2017) models exhibit a higher SFR. The best-fit relation of Olsen et al. (2017) is shown with the green solid line (for Z = 0.1 Z ). Overall, the [Cii] emission from the dwarf galaxy merger simulations of L20 and their corresponding SFR values are in very good agreement with observational trends and particularly with the DGS survey (Cormier et al. 2015(Cormier et al. , 2019. Notably, high-redshift galaxies with z∼5 have been observed to satisfy the [Cii]-SFR relation as local (z∼0) starbursts do (Herrera-Camus et al. 2021).

CONCLUSIONS
We perform [Cii] synthetic observations in SPH simulations of low metallicity (Z = 0.1 Z ) dwarf galaxy mergers, focusing on the inner 1 kpc radius where star formation is taking place. Over time, the SFR spans more than three orders of magnitude, thus providing a useful collection of [Cii]-SFR and FIR-SFR pairs for comparison against observations. In our analysis, we consider a lower observational limit of L CII = 0.5 L , which corresponds to W(CII) ∼ 0.6 K km s −1 , for a uniform 2D-grid resolution of 1024 2 . We find the following results: 1. For systems with properties similar to the modeled ones, the emission of [Cii] is optically thin. L CII increases during the two merging stages, following the trend of SFR.
2. The simulation is in very good agreement with the Kennicutt (1998)  Further investigations of similar models under similar resolution will help understand the correlation of [Cii] emission with SFR as well as with the global ISM conditions in extragalactic objects with properties similar to the simulated dwarf galaxies. In addition, different parameters in the galaxy formation and evolution model can lead to significant changes in the properties of the ISM and the star cluster formation (e.g. Hopkins et al. 2012;Buck et al. 2019;Li et al. 2020;Hislop et al. 2022). These can all in turn affect the star-formation rate and also the Cii luminosity. Thus more simulations may be needed in order to have a deeper understanding of the results presented here.  (Harris et al. 2020), Matplotlib (Hunter 2007), RADMC-3D (Dullemond et al. 2012), SPHGal (Hu et al. 2014, Wolfire, M. G., McKee, C. F., Hollenbach, D., et al. 2003, ApJ, 587, 278 Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976Zhukovska, S., Dobbs, C., Jenkins, E. B., et al. 2016 APPENDIX A. RADMC-3D CALCULATIONS We perform radiative transfer calculations in selected snapshots, using the publicly available code radmc-3d 11 (Dullemond et al. 2012) and adopting the Large Velocity Gradient approximation (Shetty et al. 2011). The abundances of C + , H, and H 2 as well as the gas temperatures and the gas velocities are taken directly from the hydrodynamical simulation. The rate coefficients for the excitation of C + and its collisions with ortho-H 2 , para-H 2 , H and e − are taken from the Leiden Atomic and Molecular Database 12 (LAMDA; Schöier et al. 2005). We considered a uniform three dimensional grid with a resolution 256 3 . The output spectra cubes have 201 channels and span ±200 km s −1 , giving a spectral resolution of dv = 2 km s −1 . The Doppler-catching switch is considered to account for velocity jumps between cells. We assume that the line is broadened thermally and due to microturbulence with equal contributions. To obtain the brightness temperature, we convert the radmc-3d line intensity using the Planck function in the Rayleigh-Jeans limit.
The computational box used in radmc-3d has a volume of (2 kpc) 3 , containing the ISM of the inner 1 kpc and centered on the merging site. For each snapshot, we perform radiative transfer calculations along three different lines of sight (along x−, y−, and z− axis) to account for the effects due to viewing angle. For each viewing angle, we convert the velocity integrated emission calculated with radmc-3d to Cii luminosity, L CII , using the expression: where k B is Boltzmann's constant, ν the rest frequency of [Cii], c the speed of light, W CII the emission of the i-th pixel and A i its area. Each 2 × 2 kpc 2 map in the radmc-3d calculations contains 256 2 pixels covering equal areas.
B. MASS-WEIGHTED PHASE-PLOTS Figure 14 shows mass-weighted plots for snapshots at t = 40 and t = 170 Myr. As can be seen, the t = 40 Myr indicates a density range of the WNM component similar to that reported for Milky Way (e.g. Wolfire et al. 1995Wolfire et al. , 2003. This is also in agreement with the phase-plot presented in Lahén et al. (2019) (see their Figure 1 covering a much lower density range and a much higher gas temperature range). For the t = 170 Myr snapshot, it can be seen that the origin of most of WNM mass is for densities −1 log n H < 3 in the 3 < log T gas < 4 in temperature range.
Evidently, the emission of Cii originates from this ISM component which, especially during the merger, contains higher densities than expected from Milky Way observations. C. ANALYTICAL SCHEME TO APPROXIMATE THE CII COOLING FUNCTION Here, we outline how to calculate the Cii cooling rate analytically. The outlined rates are applicable for high gas temperatures. The rates of collisional de-excitation with e − , H and H 2 as colliding partners are as follows: The rates in the above relations are measured in units of cm 3 s −1 . The total de-excitation (R c,DEX ) and excitation (R c,EX ) rates are given respectively by the expressions:  The ∝ nH and ∝ n 2 H relations are plotted for comparison. Right panel: Theoretical vs simulation CII cooling. The red solid line is the y = x function to guide the eye. As can be seen the majority of ΛCII function is well reproduced following the analytical expressions discussed in Appendix C.
The above rates are in units of s −1 .
The excitation temperature of Cii is T CII = hν/k B = 91.25 K. In case the CMB temperature is higher than T CII /5, there will be some contribution due to stimulated emission. The contribution is negligible for the work presented in this paper. For the spontaneous emission, we make an escape probability ansatz and use the LVG approximation. Assuming small optical depth of Cii, the rate for spontaneous emission reduces to the corresponding Einstein-A coefficient R s = A CII = 2.291 × 10 −6 s −1 . (C7) Collisional de-excitation and spontaneous emission rates are added to get the total emission rate, while the collisional excitation rate is the only contribution to the total excitation rate in case the stimulated emission is negligible. Hence, we get R tot,excite = R c,EX R tot,emit = R c,DEX + R s .
From that we define asĖ tot,excite = R tot,excite R tot,excite + R tot,emit × R s × E CII (C10) E tot,emit = R tot,emit R tot,excite + R tot,emit × R CMB,EX × E CII ∼ 0, using E CII = k b T CII = hν CII = 1.25988 × 10 −14 erg and assuming that T CMB T CII . The total cooling rate Λ CII is then Λ CII = (Ė tot,excite −Ė tot,emit ) × χ CII × n TOT ∼Ė tot,excite × χ CII × n TOT , in units of erg cm −3 s −1 , where χ CII × n TOT is the number density of C + particles in the volume of interest. The left panel of Fig. 15 plots the above equation versus the local number density, while the right one shows how it compares with the simulation result. The simulation Λ CII data are taken from the snapshot at t = 170 Myr and within 1 kpc from Cluster-3. D. Λ CII COOLING FUNCTION AROUND CLUSTER-3 Figure 16 shows 2D histograms of Λ CII versus n H within 0.1 kpc, 0.2 kpc and 0.3 kpc from Cluster-3 (see Fig. 9 for a visualization of the region). The ISM gas that is within 0.1 kpc has a considerable amount thermalized and thus collisionally de-excited. Since this part grows ∝ n H , the emission cannot compensate with the ∝ T 6 dust growth of FIR luminosity (see Eqn. 5), thus decreasing the [Cii]/FIR ratio leading to a [Cii]-deficit gas. As we increase in radial distance from Cluster-3, Λ CII comes primarily form the lower density medium which grows ∝ n 2 H , thus increasing the [Cii]/FIR ratio. Figure 16. The Cii cooling function versus nH of the ISM gas within 0.1 kpc (left panel), 0.2 kpc (middle panel) and 0.3 kpc (right panel) around Cluster-3. Dot-dashed line is the ∝ n 2 H and dashed line the ∝ nH relations to guide the eye. As can be seen, the ISM close Cluster-3 is by a considerable amount thermalized and thus collisionally de-excited. This causes the gas to be [Cii]-deficit in the vicinity of Cluster-3.

E. EFFECT OF USING A DIFFERENT LOWER OBSERVATIONAL LIMIT FOR L CII
Throughout the paper, we have assumed 0.5 L as a lower observational limit for the Cii luminosity. Based on this assumption, the observational surface (Σ) has been estimated which was used to calculate the Σ CII , Σ FIR and Σ SFR quantities. Here, we explore the response of the aforementioned variables if a different lower limit was adopted. In particular, we explore the cases of L CII > 0 L (all material capable of emitting [Cii]), L CII > 0.1 L and > 1 L . The corresponding results are shown in Fig. 17. The top left panel shows the time evolution of the observational surface when using the different L CII limitations. The top right panel shows the response in the Σ CII − Σ SFR plane. Similarly, the bottom panels show the response in the [Cii]/FIR -Σ FIR and Σ SFR planes. As can be seen, in all cases the trends and the [Cii]/FIR ratio remain unaffected. As the lower L CII limit increases, the observational surface decreases leading to a higher Σ SFR , Σ FIR and Σ CII values. This makes our results in the corresponding panels to drift rightwards. In all panels, the blue colour is for LCII > 0 L , orange for LCII > 0.1 L , green for LCII > 0.5 L (the one we consider in the main text) and red for LCII > 1 L . As can be seen, the observational limit does not affect the trends and the overall results presented.