Diagnosing the interstellar medium of galaxies with far-infrared emission lines II. [C ii ], [O i ], [O iii ], [N ii ] and [N iii ] up to z=6

Context. Gas cooling processes in the interstellar medium (ISM) are key to understanding how star-formation processes occur in galaxies. Far-infrared (FIR) ﬁne-structure emission lines can be used as a tool to understand the gas conditions and trace the di ﬀ erent phases of the ISM. Aims. We model the most important far-infrared (FIR) emission lines throughout cosmic time back to z = 6 with cosmological hydrodynamical simulations. We study how di ﬀ erent physical parameters, such as the interstellar radiation ﬁeld (ISRF) and metallicity, impact the ISM phases traced by FIR line luminosities and connect those with the star-formation rate (SFR). Methods. We implement a physically motivated multi-phase model of the ISM by post-processing EAGLE cosmological simulation with C loudy look-up tables. In this model, we assume four phases of the ISM: dense molecular gas, neutral atomic gas, di ﬀ use ionised gas (DIG) and H ii regions. Results. Our model shows good agreement with the observed luminosity–SFR relation up to z = 6 in the FIR emission lines analysed and we also provide linear ﬁts. Our predictions also agree with observations in terms of diagnostic diagrams involving various line ratios. Conclusions. We ﬁnd that [C ii ] is the best SFR tracer of the FIR lines even though it traces multiple ISM phases, while [O iii ] and [N ii ] can be used to understand the DIG-H ii balance in the ionised phase. In addition, line ratios like [C ii ] / [O iii ] and [N ii ] / [O i ] are useful to trace parameters such as ISRF, metallicity and speciﬁc star-formation rate. These results help to interpret observations of FIR line emission from the local Universe to high-z galaxies.


Introduction
After the Universe was largely ionised (the period known as "cosmic dawn" at z 6-8), the composition of the gas and its cooling gradually changed, affecting the star-formation processes in galaxies (Dayal & Ferrara 2018).Since then, star formation and black hole accretion processes co-evolved with cosmic time and shaped the evolution of galaxies (Madau & Dickinson 2014).The study of the interstellar medium (ISM) in the local Universe allows us to comprehend the current star-formation processes, but the ISM gas cooling budget may not be the same at earlier cosmic epochs (Carilli & Walter 2013).Recent observational data have opened new paths for describing and understanding the ISM gas processes of local galaxies (e.g.Malhotra et al. 2001;De Looze et al. 2014;Cormier et al. 2015).However, the complete picture of how the ISM evolves over cosmic time and how its conditions are connected with star formation in galaxies is not well understood.
Far-infrared (FIR) fine-structure emission lines (Table 1) are dominant in the gas cooling of the ISM and can help us to understand the star-formation processes, from a theoretical (e.g.Tielens & Hollenbach 1985;Goldsmith et al. 2012;Wolfire et al. e-mail: ramos@astro.rug.nl2022) and observational perspective (e.g.Ferkinhoff et al. 2010;Carilli & Walter 2013;Cormier et al. 2015;Díaz-Santos et al. 2017;Herrera-Camus et al. 2018b).These lines are less affected by dust extinction than optical lines and, at high redshift (hereafter high-z), are shifted to the (sub-)mm wavelength range accessible to ground-based telescopes (Hodge & da Cunha 2020;Förster Schreiber & Wuyts 2020).The most important ISM emission cooling line is [C ii] at 158 µm.This line traces different phases of the ISM: photo-dissociation regions (PDRs), Hii regions, diffuse ionised gas (DIG; also known as the warm ionised medium (WIM), e.g.Haffner et al. 2009;Kewley et al. 2019), molecular clouds, and the cold and warm neutral media (CNM and WNM, respectively) (e.g.Cormier et al. 2012;Croxall et al. 2017;Abdullah et al. 2017;Abdullah & Tielens 2020).[C ii] is thus a very important cooling line in the range of 20-8000 K due to its low ionisation potential (11.3 eV compared to 13.6 eV for hydrogen: Gong et al. 2012;Goldsmith et al. 2012).Furthermore, it is easily observable, as its luminosity is around 1% of the FIR luminosity of galaxies (e.g.Stacey et al. 1991;Brauher et al. 2008).Other important FIR lines include: i) atomic oxygen ([O i]) at 63 and 145 µm, which traces the denser and warmer neutral ISM environments important for star formation (Malhotra et al. 2001;Goldsmith 2019) (Nagao et al. 2011).Using these lines to understand the ISM evolution of galaxies requires a self-consistent model for all of these FIR lines over the z = 0-6 range.
For many decades, the lack of suitable instruments has hampered observations of these FIR lines in high-z galaxies (Inoue et al. 2014).Fortunately, observations taken with telescopes such as IRAM, CSO, Herschel and ALMA have provided the first high-z detections of lines like [C ii], [N ii], [O i] and [O iii] (e.g.Maiolino et al. 2005;Ferkinhoff et al. 2010Ferkinhoff et al. , 2011;;Inoue et al. 2016;Uzgil et al. 2016).Moreover, recent emission line surveys like ALPINE (Le Fèvre et al. 2020) and REBELS (Bouwens et al. 2021) are gathering data for larger samples of high-z galaxies, which are ideal to diagnose the ISM of galaxies over cosmic time.
With these new observations, different tools can be used to describe the different physical conditions in the ISM of highz galaxies.The most common and accessible tool is the use of emission line ratios that reflect the physical conditions of the ISM.For example, the [C ii]/[N ii] ratio has been used to describe and constrain whether Hii regions and/or PDRs contribute to the ISM phases of galaxies (Decarli et al. 2014), and can be used to estimate the amount of ionised gas in [C ii] (Croxall et al. 2017).Another useful line ratio is [O iii]/[C ii], used to understand ionised and neutral gas in high-z galaxies (e.g.Harikane et al. 2020;Carniani et al. 2020).This ratio has the advantage that [O iii] can be brighter than [C ii] at redshifts around the reionisation epoch (z 7, Inoue et al. 2014Inoue et al. , 2016) ) and can be observed efficiently with ALMA (Bouwens et al. 2021).Finally, the [O iii]/[N iii] ratio is used to estimate the gas metallicities (Nagao et al. 2011;Herrera-Camus et al. 2018a), although other line ratios can also be used for this purpose (Fernández-Ontiveros et al. 2021).
A more sophisticated tool to describe the physical processes of the ISM is the use of line luminosity predictions from simulations or analytical models.Simple models involve ratios between FIR emission lines and/or the total luminosity in the IR to obtain physical conditions such as hydrogen density, FUV radiation flux or stellar temperatures (e.g.Malhotra et al. 2001;Ferkinhoff et al. 2011), while more complex models use magnetoradiation hydrodynamics simulations of the Universe with different radiative transfer codes to predict emission line luminosities (e.g.Olsen et al. 2021;Katz et al. 2022;Pallottini et al. 2022, and references therein).Some of these studies focus on specific emission lines such as [C ii], [O i] or [O iii] in analytic models (e.g.Goldsmith 2019;Ferrara et al. 2019;Yang & Lidz 2020) and simulations (e.g.Moriwaki et al. 2018;Leung et al. 2020;Ramos Padilla et al. 2021), while others study interesting line ratios like [O iii]/[C ii] (Arata et al. 2020;Vallini et al. 2021).A large effort has also been made to model various FIR emission lines in models at different cosmic times (e.g.Vallini et al. 2013Vallini et al. , 2015;;Olsen et al. 2015Olsen et al. , 2017Olsen et al. , 2021;;Popping et al. 2019;Sun et al. 2019;Katz et al. 2019Katz et al. , 2022;;Pallottini et al. 2019;Yang et al. 2021bYang et al. , 2022)).However, these studies do not examine the emission of these FIR lines in a consistent way in terms of their redshift evolution, from the local Universe to the epoch of reionisation, due to computational constraints or focus on certain cosmic epochs (with the exception of Popping et al. (2019) which modelled the [C ii] line).Therefore, we need better models to understand current observations consistently.
With this in mind, we aim to predict luminosities of the main FIR lines in a cosmological context through the use of cosmological hydrodynamical simulations to infer the physical conditions of galaxies across a wide range of redshifts.The goal of this paper is to test the impact of physical parameters on the FIR emission lines tracing different ISM phases in galaxies.We will use these predictions as diagnostic tools, which will be useful for both current and future observations.To do this, we model the emission of FIR lines by post-processing the hydrodynamical simulations of the Evolution and Assembly of GaLaxies and their Environments (EAGLE) project (Schaye et al. 2015;Crain et al. 2015) with a physically motivated model of the ISM presented in Ramos Padilla et al. (2021, hereafter Paper I).We use Cloudy (Ferland et al. 2017) cooling tables (Ploeckinger & Schaye 2020) to predict the emission from different ISM phases in galactic environments.Throughout this paper, we assume the ΛCDM cosmology from the Planck Collaboration et al. ( 2014) results (Ω = 0.307, Ω Λ =0.693, H 0 = 67.7 km s −1 Mpc −1 and σ 8 = 0.8288).
In this paper, we first briefly describe the simulation data and the ISM model that we use to predict FIR emission line luminosities (Sect.2).Next, we present the results of the FIR emission line predictions from the simulations and how they compare with the observations from the local Universe all the way out to z∼6 (Sect.3).In Sect. 4 we evaluate some FIR diagnostic diagrams used in various high-z studies.After that, we discuss the potential systematic uncertainties that can affect the predictions and the comparisons with observations (Sect.5).Finally, we present our conclusions in Sect.6.

Methodology
In this section, we first describe the sets of simulations that we use in this work (Sect.2.1), then we briefly explain the initial model used to characterise the structure of the ISM (Sect.2.2).Finally, we present in detail the addition of Hii regions as a new ISM phase (Sect.2.3) in our model.
For this work, based on the results from Paper I, we used two simulations from the EAGLE suite: Ref-L0100N1504 and Recal-L0025N0752, as described in Table 2.The main differences between the two simulations are the box-size of the simulation (100 and 25 cMpc (comoving Mpc), respectively), the mass resolution (∼10 6 M and ∼10 5 M , respectively), and the calibration of the physical parameters of the subgrid routines to reproduce the galaxy stellar mass function (GSMF; Schaye et al. 2015).Both simulations are similar in terms of "weak convergence", which means numerical results converge in different simulations after re-calibrating the sub-grid parameters (Furlong et al. 2015;Schaye et al. 2015).
In this study, we retrieve simulated galaxies from the SPH data (The EAGLE team 2017) by using textscFoF (Friends-of-Friends) and SUBFIND algorithms (Springel et al. 2001;Dolag et al. 2009) in the dark matter halos.With these algorithms, the sub-halos containing the particle with the lowest value of the gravitational potential are called "central" galaxies.We focus on these "central" galaxies to estimate line luminosities.We use galaxies with at least 300 star particles (i.e.stellar masses higher than ∼10 8 M and ∼10 8.5 M for Recal-L0025N0752 and Ref-L0100N1504, respectively) within 30 pkpc (proper kpc) from the centre of the potential.We selected our sample of galaxies from the EAGLE database (McAlpine et al. 2016) in redshifts (the closest snapshot) between z = 0 and z = 6 in steps of ∆z = 1, using a single snapshot at each redshift (the snapshot closest in time; e.g. at z = 6, we use the z = 5.97 snapshot from EAGLE), where the total number of galaxies depends on the simulation used.The z = 6 cutoff was selected due to the availability of observational data and the number of galaxies recovered in EAGLE required to compare them statistically.For Recal-L0025N0752, we select all retrieved galaxies, while for Ref-L0100N1504, we randomly select up to 1 000 galaxies that fulfil the previous conditions per redshift.In the last columns of Table 2, we present the total number of galaxies used per redshift slice in each of the simulations.
We selected a total sample of 8 277 galaxies simulated with EAGLE at redshifts between z = 0 and z = 6.In each of those galaxies, we model the emission coming from the eight FIR emission lines (Table 1) that trace different phases of the ISM.

The multi-phase ISM model
To predict these emission lines, we used the ISM model presented in Paper I, where the luminosity estimations from the [C ii] emission line at 158 µm in the local Universe (z = 0) showed good agreement with observations.However, this model must be improved if we want to properly account for other FIR lines such as [O iii] and [N iii], which probe denser ionised regimes.Therefore, we add Hii regions as a new phase in our ISM model.In Fig. 1, we illustrate the path from the EAGLE simulations (Sect.2.1) to the total luminosity of the lines in the current model.The main difference between Paper I and this work is the addition of the contribution to the line luminosities from Hii regions (Sect.2.3).In this subsection, we briefly explain the model presented in Paper I.
After selecting the sample of galaxies for which we want to calculate the line luminosity, we retrieve and post-process the gas and star particle data of the galaxies in the simulation.Physical properties such as total hydrogen number density were estimated for all gas particles using the information available in the SPH particle data.For example, we calculated the total hydrogen number density as n(H) = ρX H m H , with m H the hydrogen mass, ρ the density and X H the SPH weighted hydrogen abundance.In addition, the fraction of neutral hydrogen was estimated for all gas particles according to Rahmati et al. (2013), following ionisation equilibrium as described in Sect.2.2.1 in Paper I. We calculated the background interstellar radiation field (ISRF) from the star formation rate (SFR) surface density of the gas particles and the local ISRF from the star particles as described by Olsen et al. (2017).The local ISRF is estimated with starburst99 models (Leitherer et al. 2014) for star particles with a distance below the smoothing length of a given gas particle (as described in Sect.2.2.1 in Paper I).The sum of the background ISRF and the local ISRF defines the total ISRF impinging on the neutral cloud.With the fraction of neutral hydrogen, we split the gas into neutral and ionised components.The neutral gas mass and the ISRF are used to estimate the sizes of the neutral clouds, while the ionised mass is used to describe the diffuse ionised gas (DIG).
The neutral clouds are defined as concentric spheres with densities following a Plummer profile.Each ISM phase in the neutral cloud is defined by a different radius as described in Sect. 2. in Paper I. The transition between atomic and molecular hydrogen defines the limit between the neutral atomic gas and the dense molecular gas, while the transition between ionised to neutral carbon defines the limit for the inner core region of the cloud that is completely shielded from FUV radiation.This definition assumes that the inner core can only be traced by neutral species such as CO, and the inner core is therefore ignored in the estimation of the total luminosity of the FIR emission lines1 .
On the other hand, the volumetric structure of the DIG is assumed to be spherical with a radius drawn from a smoothed broken power-law distribution.The parameters of the powerlaw distribution are derived based on calibrating the [N ii] line (at 122 and 205 µm) predictions of 492 galaxies from the Ref-L0100N1504 and 200 galaxies from Recal-L0025N0752 simulations to the observational dataset from Spinoglio et al. (2015), Fernández-Ontiveros et al. (2016) and Cormier et al. (2019) (containing 8, 53 and 2 galaxies, respectively, of a sample of 63 galaxies).This led us to DIG clouds with average sizes of ∼900 pc for Recal-L0025N0752 galaxies and ∼2 kpc for Ref-L0100N1504, which is around 3 times the maximum softening length of the original simulations.
These structures define the phases of the ISM that contribute to the total luminosity of the emission line: DIG, neutral atomic and dense molecular gas.The estimation of the luminosities for each phase is obtained by using Cloudy (v17.01, Ferland et al. 2017) cooling tables of shielded and optically thin gas (Ploeckinger & Schaye 2020).The sum of those phases then gives us the total luminosity for a FIR emission line.For a complete description of the model, we refer the reader to Sect.2.2 and 2.3 in Paper I.
Table 2. EAGLE simulations used in this work.The box-size, number of particles and gas particle mass define the resolution of the simulation.The right columns show the number of galaxies used in this work for a given simulation at each redshift.

Name in
Box-size Number of SPH Gas mass Number of galaxies Schaye et al. (2015) (cMpc) particles (M )

Hii regions as a new ISM phase
The luminosity estimations for FIR lines such as [O iii] and [N iii], which require high ionisation potentials (35.1 and 29.6 eV, respectively), can only be found in Hii regions and other dense, ionised regimes.The addition of this phase as an ISM component to the model allows us to compare the predicted luminosities of these lines with observations.Therefore, we update the model described in Paper I to simulate the Hii regions production of the most prominent FIR lines, including as shown in Table 1, in order to infer the ISM conditions in galaxies from the local Universe out to z ≈ 6.
In Paper I, we use starburst99 to calculate the ISRF from the stars close to gas particles (the distance between the particles is less than one smoothing length).Now, we also use the spectrum from starburst99 to calculate the emission coming from Hii regions.To generate this spectrum, we adopted the Geneva stellar models (Schaller et al. 1992) with standard mass loss for five metallicities (Z = 0.001, 0.004, 0.008, 0.020, and 0.040).We split starburst99 grid values in young (≤ 100Myr) and old ages (> 100Myr), to improve our estimations at younger ages.For young ages, we estimate the parameters every 1 Myr, while for old ages we calculate 100 steps on a logarithmic scale up to 10 Gyr.We assume a total stellar mass of 10 4 M , with a Kroupa initial mass function (IMF).The SPH star particles are divided to match the stellar mass from the starburst99 models, assuming a random exponential distribution of the original ages of the SPH particles.By doing this we try to avoid the poor sampling of star formation that can affect luminosity estimates, especially in Hii regions.
We use a a photoionisation model from Cloudy to simulate the line emissitivites of Hii regions based on starburst99 spectra of their underlying stellar populations.We calculate the stellar atmospheres, or spectral energy distributions (SEDs), in the [Gyr] −3.0 0.9 0.1 photoionisation models for a given age, metallicity and the ionising photon flux (Q).In this way, the stellar mass and metallicity from the star SPH particles are used to obtain Q coming from a stellar cluster from the starburst99 grids.Thus, we use these three physical parameters (age, metallicity and Q) to construct spherical clouds where the emissivities depend only on them.
The range of values for these parameters is presented in Table 3, totalling 25 200 grid points for all redshifts in this work.
In terms of the structure, we assume a fixed density of ∼30 cm −3 to resemble classical Hii regions similar to the Strömgren sphere, where densities are in the range of 10-100 cm −3 (Draine 2011).Choosing a higher or lower density affects the Hii luminosity of some lines as we describe in Sect. 5.In ionisation balance the Strömgren radius is with Q the ionising photon flux in units of s −1 , n(H) the total hydrogen number density in units of cm −3 and α B is the Case B recombination rate coefficient.These Hii regions are radiation bound in terms of ionised hydrogen.Once the ratio of ionised hydrogen to the total hydrogen density drops below 5%, the calculation stops, which defines the outer radius.The inner radius is set to 1% of the expected Strömgren sphere radius (Eq.1), as in Yang & Lidz (2020), allowing us to obtain "ultracompact" Hii regions, which can have sizes of ∼0.03 pc (Kurtz 2005).We calculate the emissivities without further iteration on the assumed density, which is adequate for clouds with densities typical of Hii regions (Ferland 1996).For the rest of the parameters in Cloudy we choose a similar configuration as the one used by Ploeckinger & Schaye (2020).
The starburst99 outputs are also used to calculate the L FUV and Q for each star SPH particle.Thus, the look-up luminosity tables from Cloudy are used to interpolate the star SPH particles in terms of age, metallicity and Q.As a result, we obtain the respective Hii luminosity for each star particle for the FIR emission lines.This luminosity is then used as part of the total contribution of the ISM phases (DIG + neutral atomic gas + dense molecular gas + Hii regions) for a given simulated galaxy.

Individual line luminosities
In this section, we present the predictions of individual lines from our model and compare them with observations.For each line, we first examine the relation between the line luminosity and the SFR of the galaxy at different cosmic epochs.Our predictions of the line luminosity-SFR relation are also compared with observational measurements collected and homogenised from published work (see Appendix A).Then we check the contribution to the line luminosities from each of the ISM phases, i.e., the DIG, Hii regions, neutral atomic and dense molecular gas.We mainly discuss the results of the five FIR emission lines: and [N iii] at 57 µm.The other three emission lines listed in Table 1 behave similarly to their aforementioned pairs (e.g.[O iii] at 52 µm is similar to [O iii] at 88 µm).We release all the data in a Zenodo repository as described in Appendix B.

[C
The most important and brightest of these FIR lines is [C ii] at 158 µm.This line follows a clear trend with SFR (e.g.Stacey et al. 1991;Brauher et al. 2008;Stacey et al. 2010), and it is used as a SFR tracer at different redshifts (e.g.De Looze et al. 2014;Herrera-Camus et al. 2015;Schaerer et al. 2020).In Fig. 2, we show the relationships between SFR and [C ii] luminosity (L [C ii] ) for the seven redshift slices analysed in this work.We compare the predictions from the EAGLE simulations Ref-L100N1504 and Recal-L0025N0752 with predictions from other simulations (Olsen et al. 2015(Olsen et al. , 2018;;Katz et al. 2022), semi-analytic models (Lagache et al. 2018;Popping et al. 2019) and linear relationships derived from observations (De Looze et al. 2014;Schaerer et al. 2020).We also plot the linear relationships that we infer from our model (see Appendix C) and extrapolate them outside the dynamic range covered by the simulations but within 10 −3.5 M yr −1 < SFR < 10 3.5 M yr −1 .With this extrapolation, we can compare the observations with high SFR values that the simulations do not cover.This is necessary as high-z observations generally only include galaxies with very high SFR due to sensitivity limits.In general, the agreement between our models, observations and other models is good within the typical scatter ∼0.4 dex (De Looze et al. 2014), especially at z = 0 and z = 6, where there are more observational constraints.
In Ramos Padilla et al. (2021) we showed that the SFR-L [C ii] relationship at z = 0 could be reproduced with a model similar to that implemented in this work.Therefore, it is not surprising that the current model still reproduces this SFR-L [C ii] relationship.Compared to Popping et al. (2019), the only other model that covers the same redshift range as this work for [C ii], we find a flatter slope in the SFR-L [C ii] relationship.Especially at z = 1-3, the differences in the slopes can lead to up to ∼1.8 dex change in L [C ii] at SFRs around 1000 M yr −1 , but the relationships are more similar at other redshifts.The reason for these discrepancies may reside in the different galaxy formation physics in EA-GLE and the Santa Cruz SAM used in Popping et al. (2019).Unfortunately, the comparison between the two galaxy formation models is out of the scope of this work.
At z = 2, the linear relationship of Olsen et al. (2015) predicted from seven simulated galaxies shows a behaviour similar to that of Popping et al. (2019) and agrees with the estimated scatter of our linear regression (0.2 dex) at SFRs around 1 M yr −1 .Over the redshift range z = 1-3, the extrapolation of our relation show a potential small offset compared with observations.However, this discrepancy is not significant, taking into account the small sample size, large scatter (around 0.5 dex), potential bias towards galaxies with high line luminosities, and systematics in deriving luminosities and SFR.Furthermore, altering the assumptions in our modelling process could affect our predictions, as we show in Sect. 5.
At z = 4-6, most of the models and observations match well with the linear relationships derived from our predictions.Similar predictions from Vallini et al. (2015), Leung et al. (2020) and Carniani et al. (2020), which are not shown in the plot, also agree at z = 6, indicating that the SFR-L [C ii] relationship can be tight at higher redshift, even if there is some scatter in the data related to observational errors (Schaerer et al. 2020).This demonstrates that our physically motivated model of the ISM is valid not only for estimating the luminosity of [C ii] in the local Universe but also up to z = 6.

[C ii] deficit
Although a linear SFR-L [C ii] relation is commonly used to assume that L [C ii] is a good SFR tracer, observational data of local galaxies show a decrease of L [C ii] at IR luminosities above 10 12 L , known as the "[C ii] deficit" (e.g.Díaz-Santos et al. 2017;Herrera-Camus et al. 2018b).However, recent observations of z > 4 galaxies show no evidence of such a deficit (Matthee et al. 2019;Carniani et al. 2020;Schaerer et al. 2020, e.g.).In Ramos Padilla et al. (2021), we examined this deficit by comparing the SFR-L [C ii] relationship with deviations from the star-forming main-sequence (MS), following the suggestion of Herrera-Camus et al. (2018b).In this work, we compare the L [C ii] /SFR ratio with the offset from the MS (∆MS) in Fig. 3.We estimate the specific SFR (sSFR = SFR/M ) of our galaxies and normalise this by the derived sSFR for the MS from Speagle et al. (2014), which gives us the ∆MS.
We find that the L [C ii] /SFR ratio almost always decreases with ∆MS.If the decrease in the L [C ii] /SFR ratio extends to a higher ∆MS, it supports the idea that starburst galaxies may not follow the SFR-L [C ii] relationship shown in Fig. 2 (Herrera-Camus et al. 2018b;Ferrara et al. 2019).Therefore, L [C ii] may not be a good tracer of SFR for starburst galaxies, as we discuss in Sect.3.6.

Contribution to L [C ii] from each ISM phase
We present the average contribution of the ISM phases to L [C ii] at each redshift as a function of SFR for the galaxies in the Recal-L0025N0752 simulation in Fig. 4. We examine the ISM phase contributions of the Recal-L0025N0752 simulation because we expect this simulation to behave more like observed local galaxies, as shown in Paper I. In general, we find that most of the [C ii] emission comes mainly from the atomic phase, especially at z > 2, in agreement with the general assumption that the neutral gas is the dominant ISM phase contributing to L [C ii] as estimated from observations (e.g.Croxall et al. 2017;Cormier et al. 2019) and suggested by models (e.g.Olsen et al. 2015Olsen et al. , 2018;;Lagache et al. 2018).However, the contribution by neutral atomic gas to L [C ii] changes with redshift, from ∼20-40% at z ≤ 1, to ∼70-90% at z ≥ 3. The most important reason for these differences is the contribution from the DIG.At z = 0, the DIG dominates (the contribution is greater than 50%) in most of the galaxies with a SFR< 0.1 M yr −1 , then the DIG contribution reduces to 30% at z = 2, and finally it is negligible at z = 6.This negligible contribution of the DIG at z = 6 does not agree with the estimated contribution of 44% by Olsen et al. (2018).This result is expected as we estimate the size of the DIG using a physical assumption instead of using the smoothing length as in Olsen et al. (2018), which leads to a more compact DIG.However, as our modelled SFR-L [C ii] relation shows a better agreement with observations, this may imply that a higher contribution from the atomic gas is needed to match the observed galaxies at z = 6.Therefore, our estimations seem to favour the atomic gas as the main responsible for the L [C ii] at z = 6.
On the other hand, the fractional contribution of Hii regions to L [C ii] reaches its maximum at z = 2, where they contribute up to 80% of the luminosity.This trend is expected as Hii regions trace young stars, and it is well known that the co-moving star formation rate density reaches its peak value at z ≈ 2 (Madau & Dickinson 2014).In general, the average contribution by Hii regions is ∼20% of L [C ii] over all z < 4. For molecular gas, the contribution is on average 10% with a maximum at z = 0.In both the molecular and Hii regions phases we obtain higher fractional contributions to L [C ii] at higher SFR.This confirms the results found by Olsen et al. (2015) and in Paper I that the contribution of a given phase to L [C ii] depends on the global SFR of the galaxy.
These predictions are not in agreement with recent results from Tarantino et al. (2021) in resolved regions of two local galaxies (M 101 and NGC 6946) where the ionised gas contribution to L [C ii] is negligible (average upper limit of 12%).This disagreement in the ISM phases may be related to the spatial constraints of the observations, as they focus on the arms regions of galaxies where denser neutral gas is expected.For example, Pineda et al. (2018) estimated the L [C ii] contribution coming from different regions inside M 51.They found that the region between the arms in M 51, where the diffuse ionised gas is expected be located, is approximately 20% of the total L [C ii] , compared with ∼80% that comes from the arms and the nucleus.In our model, we calculate the global properties of galaxies (in a 30 pkpc aperture), and thus we expect a higher contribution from the diffuse gas, as this ISM phase can cover more extended regions throughout a galaxy.The information from the diffuse phases may therefore be missing when regions close to the arms of disk galaxies are observed.In general, contributions from different phases will depend on the scales over which we observe the galaxies (Tarantino et al. 2021).
It is important to note that our model seems to overpredict the contribution of the ionised phases (DIG + Hii regions) for the Recal-L0025N0752 simulation.In Fig. 5, we examine the relation between the ionised phases to [C ii] as a function of metallicity (12+log (O/H)).We compare the median predictions of Ref-L100N1504 and Recal-L0025N0752 with the relationship obtained by Cormier et al. (2019), based on the data from their work and Croxall et al. (2017), where at a higher metallicity there is a higher contribution to L [C ii] from the ionised phase.For Ref-L100N1504, we estimate fractional contributions ∼10% higher than those inferred from the Cormier et al. (2019) fitting function, although at z = 0 our predictions can be 40% lower for higher metallicities.In Recal-L0025N0752 we see that the fractional contributions are always above the empirical fitting function, especially at z = 1 and z = 2, by around 40%.This means that there may be an overestimation in the ionised component of the Recal-L0025N0752 luminosities.
The simulations and observations presented in Fig. 5 are in any case difficult to compare.The simulated galaxies may not have the same metallicity calibration as the observed metallicities in Croxall et al. (2017) and Cormier et al. (2019), and the method by which fractional contributions of the ionised gas are calculated may be different.Nonetheless, in general, we find that the contribution of the ionised phase to [C ii] increases with increasing metallicity, as observed by Cormier et al. (2019).

The SFR-L [N ii] relationship
The emission lines at 122 and 205 µm of [N ii] are commonly used to trace ionised gas phases around neutral clouds, because its ionisation potential is only slightly above that of hydrogen.These lines are also used to disentangle the ionised gas contribu- Observational data Upp.limits obs.data Lensed galaxies Upp.limits lensed gal.

Ref-L100N1504
Recal-L0025N0752 Estimated lin.reg.and Recal-L0025N752 (right).We show the median values of the different redshifts (solid lines) and their 25th and 75th percentiles (dotted lines).We only show the bins with more than 3% of the total sample for the respective simulation.
tion to the [C ii] luminosity (e.g.Goldsmith et al. 2015;Ferkinhoff et al. 2015;Pavesi et al. 2016;Croxall et al. 2017;Cormier et al. 2019;Langer et al. 2021), as we do in this work.The relationship of these lines with SFR has been explored in the local Universe (e.g.Herrera-Camus et al. 2016;Zhao et al. 2016a); however, at higher redshifts we have very little observational data, especially for the 122 µm line.Due to this lack of data at z > 0, we focus here on the et al. (2016, eq. 10), we assume the values of the collisional excitation coefficients from Tayal (2011), and abundances close to solar.For the relation of Mordini et al. (2021), we use the sample of AGN galaxies assuming the conversion from infrared luminosities (L IR ) to SFR of Kennicutt & Evans (2012).We use the AGN sample in Mordini et al. (2021) as those galaxies also follow the relation between SFR and the luminosity of the PAH feature at 6.2 µm (see their Fig. 6).At z = 0, the luminosity predictions of [N ii] follow a similar relationship to the observational relations of Zhao et al. (2016a) and Mordini et al. (2021).A potential reason for the difference of one to two orders of magnitude with respect to the Herrera-Camus et al. ( 2016) relation is that the assumption of solar abundances is incorrect2 .If we assume an abundance below solar in the Herrera-Camus et al. ( 2016) relation, the relation is closer to our model.At z = 1 and z = 2 our model is consistent with the upper limits of the observational data.At higher redshifts, 3 < z < 6, our models agree with the observations, although the range of L [N ii] /SFR in the observations is larger than the estimations (almost 0.5 dex at z = 4).Most of the data at these redshifts come from the work of Cunningham et al. (2020), who observed 40 gravitationally lensed galaxies with the Morita Atacama Compact Array (ACA) of ALMA.The luminosities of these galaxies depend strongly on the lensing magnification factor, which could create a large scatter in the inferred luminosities (due to uncertainties in the assumed lensing model), as observed in the orange data points.Our model shows a good agreement with the observations, including the lensed galaxies after correction for magnification.

Contribution to L [N ii] from each ISM phase
In Fig. 7, we show the contribution of the different phases of the ISM to the [N ii] line at 205 µm.The luminosity of this line comes mainly from the two ionised phases (DIG and Hii regions), as expected from observations (e.g.Langer et al. 2021).Interestingly, the relative dominance of these two phases seems to change with redshift and SFR.At higher SFR most of the luminosity comes from Hii regions, while at low SFR, most of the luminosity comes from the DIG phase.The DIG can contribute significantly to the luminosity, which is clearer at lower SFR when there are not so many Hii regions.In the local Universe, the contribution from the DIG dominates (∼80%) over the contribution from Hii regions.At higher redshifts, 1 < z < 4, the contribution of the two phases is split relatively evenly.At z > 5, Hii regions dominate the line emission.In this redshift range, at the highest SFRs, Hii regions contribute significantly more than the DIG.The transition point between the phases (the SFR where DIG contributes less than 50% and Hii regions more than 50%) is around 1 M yr −1 at z = 0 and decreases to 0.1 M yr −1 at z = 6.There is also a very small contribution from the atomic phase (<4%) at some redshifts, but this is negligible compared to the ionised phases.i] 63 and 145 µm   3

[O
The [O i] emission lines at 63 and 145 µm trace warm gas in neutral clouds and are commonly detected in galaxies in the local Universe (Malhotra et al. 2001).However, the line at 145 µm is fainter than the 63 µm line due to its lower spontaneous decay rate and higher upper level energy (Goldsmith 2019).As we did with [N ii], here we focus on the results of the brighter line of [O i], i.e., the emission line at 63 µm.
In Fig. 8 we present predictions of [O i] 63 µm luminosity as a function of SFR at different redshifts.Our model at z = 0 agrees with the local Universe relationships of De Looze et al. ( 2014) and Mordini et al. (2021).For the relation of Mordini et al. (2021), we use the sample of star-forming galaxies assuming the conversion from L IR to SFR of Kennicutt & Evans (2012).At z > 1, extrapolation of the linear fits to our model predictions are consistent with the observational upper limits.There are only a few detections at z = 2, z = 3 and z = 6 that can be directly compared with our predictions.There are seven detections in the z = 2 bin, of which five are > 0.5 dex from the extrapolated linear fit.However, all of these galaxies are gravitationally lensed (Brisbin et al. 2015), and no corrections were applied to correct for lensing because the magnification factors were unknown (see also Zanella et al. 2018).Therefore, it is reasonable to assume that these five measurements should be treated as upper limits.In the other redshift bins, we observe a very good agreement with the linear relationship because the lensing magnification has been taken into account (Rigopoulou et al. 2018;Zhang et al. 2018;Rybak et al. 2020).
Other numerical models, not shown in the plots, also predict the [O i] 63 µm line and show similar results.For example, Olsen et al. (2018) predict the emission of 30 simulated galaxies at z = 6 in the MUFASA cosmological simulation.All of their galaxies have SFR∼10 M yr −1 and L  2021) at z = 0, with a difference of ∼0.5 dex at log(SFR[M yr −1 ]) = 1.We conclude that the model predictions presented in this work provide a better match with observations across a wide range of redshifts than previous models that only focus on a single redshift slice.

Contribution to L [O i] from each ISM phase
Contributions to the [O i] 63 µm line come mainly from neutral clouds, i.e., neutral atomic gas and dense molecular gas, as shown in Fig. 9.At z = 0, the contribution to L [O i] from molec-  ular gas is ∼40% while for atomic gas it is ∼39%.These percentages change with redshift: the contribution from molecular gas decreases and from atomic gas increases with increasing redshift.At z = 2 the percentages are ∼7% and ∼85%, respectively, while at z = 4 they are ∼1% and ∼98%.The contributions from the other phases are negligible, especially for the Hii regions.On average the contribution from the DIG is less than 10%; however, it can be very high (> 80%) in galaxies with very low SFR (< 10 −2 M yr −1 ) in the local Universe.At z ≥ 3, the molecular fraction, which is calculated from the line luminosity in the region defined by the radius at which the transition from atomic to molecular H occurs, in the Plummer profile (Equation 31 in Paper I) is very low.At those redshifts, even though the average density of the neutral cloud is higher than in the local Universe, the ISRF is also high, which causes the dominant emission to come from the atomic gas instead of the molecular gas.These results support the understanding that the [O i] 63 µm line originates in warm neutral environments (Malhotra et al. 2001; Goldsmith 2019) even in high-z galaxies (z∼6).(Olsen et al. 2018;Katz et al. 2022;Kannan et al. 2021).

[O
Our predicted SFR-L [O iii] relationships largely agree with observations, although it seems that our predicted relationship is steeper than some relations found in the literature.At z = 0, we compare our model with the predictions of De Looze et al. ( 2014) and Mordini et al. (2021).For the relation of Mordini et al. (2021), we use the sample of star-forming galaxies assuming the conversion from L IR to SFR of Kennicutt & Evans (2012).Both observational relationships appear to have slopes flatter than our model, separated by about an order of magnitude at SFR= 10 M yr −1 .The reason for this difference is that  at low SFRs, our luminosity predictions coming from Hii regions drop off sharply for both simulations (as we will discuss in the next subsection), leading to a steep slope even though our predicted L [O iii] values agree with the observational data at SFR= 1 M yr −1 .At z = 2-4, due to the small sample size in the observational datasets, it is very difficult to properly assess the level of agreement between the observations and our model predictions.In some cases, the observations agree well with the extrapolation of our predicted SFR-L [O iii] relations, with a few exceptions.For example, the galaxy located 1.1 dex below the extrapolation at z = 4 is AzTEC 1 (Tadaki et al. 2019).The discrepancy in this galaxy is not related to gravitational lensing, but rather this galaxy has a very low [O iii]/[C ii] ratio compared with SPT-S J041839-4751.8, a galaxy that closely follows our relation (De Breuck et al. 2019).This means that this galaxy may be an outlier and very different physical conditions can shift its position in the SFR-L [O iii] relation, as we discuss in Sect. 4.
At z = 5, we compare our results with those of Katz et al. (2022), which show a similar slope but lower L [O iii] by ∼1.3 dex.The same model was used to calculate L [O iii] at z = 6, showing that their predictions, together with the models of Olsen et al. (2018), underestimate L [O iii] compared with observations by 0.7 dex at SFR∼100 M yr −1 .This is not the case with our model: it agrees very well with the observations and with the linear relation from Harikane et al. (2020), with a difference of 0.2 dex at SFR∼100 M yr −1 .It is also interesting to note that the numerical results of Kannan et al. (2021) are very similar to those presented in this work, with the difference of a slightly higher slope in their work, which gives a difference of 0.4 dex at SFR∼10 M yr −1 .

Contribution to L [O iii] from each ISM phase
In Fig. 11 we show the contribution of the ISM phases to L [O iii] , of which the ionised ISM is the only contributor, as expected.The dominant contributor to L [O iii] are the Hii regions at most redshifts.However at z = 0, the emission coming from the Hii regions drops sharply when the SFR decreases.The reason for this sharp drop is the low ionising photon flux coming from the star SPH particles in the simulated galaxies at low SFR.This sets the [O iii] line luminosities from the Hii regions to almost negligible values compared to the DIG, which explains the trend observed in z = 0.At z = 1, the lack of ionising photon flux affects galaxies less than at z = 0, and at these redshifts the Hii regions dominate (72%) over the DIG (28%).For redshifts from z = 2 to z = 6 the contribution from the Hii regions changes from 85% to 99%.The [N iii] emission line at 57 µm is very similar to the [O iii] 88 µm emission line, due to its excitation properties (Table 1).
Both trace Hii regions, with the difference that [N iii] is fainter in galaxies with metallicities below solar (e.g.Nagao et al. 2011;Rigopoulou et al. 2018).The relationship between the SFR and L [N iii] is less well known than the other FIR emission lines presented in this work.In Fig. 12, we present the SFR-L [N iii] relationship.We compare predictions from the EAGLE simulations with the observational sample and the linear relation at z = 0 by Mordini et al. (2021) in AGN galaxies assuming the conversion from L IR to SFR of Kennicutt & Evans (2012), as we did in the case of [N ii].
At z = 0, our model predictions agree with the observations, with a mean offset of 0.2 dex, as they are inside the observational scatter for the overlapping SFR range.However, the L [N iii] predictions have the same problem as the L [O iii] predictions: the low ionising photon flux coming from SPH star particles in simulated galaxies at low SFR.This leads to a steeper relation of L [N iii] with SFR than that of Mordini et al. (2021), which follows the other observational results.At z ≥ 1, even though there is not much information, we find that the extrapolation of our linear fit is consistent with the upper limits and detections in two galaxies at z ≈ 2 (H-ATLAS J091043.1-000321,Lamarche et al. 2018) and z ≈ 3 (HERMES J105751.1+573027,Rigopoulou et al. 2018).These two galaxies are magnified by around an order of magnitude (magnification factors of 11.5 and 10.9, respectively), therefore the agreement with our model, after applying the magnification correction, shows that the slope of the linear relationship may be correct for these redshifts.

Contribution to L [N iii] from each ISM phase
The contributions of the ISM phases to the [N iii] 57 µm emission line are very similar to the contributions of the phases to the [O iii] 88 µm line (Fig. 11), as shown in Fig. 13.The main difference resides in the exact percentages between the two dominant ionised ISM phases (DIG and Hii regions).At z = 0, the contributions of the DIG and Hii regions to L [N iii] are 38% and 62%, respectively.At higher redshifts these percentages change to 15% and 85% at z = 2, 8% and 91% at z = 4, respectively.Finally, at z = 6 the Hii regions are responsible for almost all the

Summary of FIR line luminosities
The relationships of the SFR with the different FIR line luminosities presented in Sects.3.1-3.5 depend on redshift (Figs. 2,6,8,10 and 12).We assume a potential redshift dependency in the linear fits of our line predictions (Eq.C.1), which shows a good agreement with the observations.Nonetheless, we want to quantify how these FIR lines and their dependence on SFR evolve with redshift.To do this, we plot in Fig. 14 the ratio between luminosity and SFR (L/SFR) versus redshift for the eight lines modelled in this work.First, we compare the results of the two sets of simulations Ref-L100N1504 and Recal-L0025N0752.Taking into account the scatter of the predictions, we find a good agreement between them, even though some of the physical properties of galaxies in the simulations are different, as we discuss in Sect. 5.
The evolution of the L/SFR ratio is almost flat for most FIR lines at z ≥ 4. At lower redshifts, 1 < z < 3, however, there are drastic changes in the L/SFR ratio.For example, the [O iii] lines at 52 and 88 µm have higher L/SFR values at z = 2 and then decrease sharply towards z = 0, by almost an order of magnitude.The opposite occurs for [C ii], where at z = 2 the L/SFR ratio has lower values than at other redshifts, although the difference is less than 0.5 dex.
These changes are related to the effects present in galaxies during "cosmic noon", where the cosmic star-formation history reaches its peak value at z ≈ 2 and around half of the stellar mass of the local Universe was formed, affecting the different phases of the ISM traced by the FIR lines in Fig. 14 (Madau & Dickinson 2014;Förster Schreiber & Wuyts 2020).This result is expected because the EAGLE galaxies reproduce the observed trend in the SFR density and stellar mass assembly across cosmic time (Furlong et al. 2015) and so the ISM phases evolve accordingly with redshift (Figs. 4,7,9,11 and 13).
The two [O iii] lines have a similar shape in Fig. 14, which explains why both lines can be used to constrain gas density and metallicity at high-z (Yang et al. 2021a).The [N ii] pair shows a similar behaviour, where the main difference resides in the value of the L/SFR ratio at z = 0, although the scatter is especially large at this redshift.This could be important when estimating electron densities from observations from the ratio of these two [N ii] lines (Croxall et al. 2017;Langer et al. 2021).Finally, for both [O i] lines, we see a clear difference between the estimated luminosities of the lines of around 1.15 dex over most of the redshift range.This difference of more than an order of magnitude is expected from theoretical models.If the difference in [O i] luminosities is more than 1.15 dex, it may indicate higher kinetic temperatures (> 400 K) and/or lower gas densities ( 10 cm −3 ) in observations (Goldsmith 2019).
With this information, we can ask which FIR emission line is the best SFR tracer across cosmic time.The FIR line showing the least variation with z is [C ii].However, this tracer may not be ideal in some cases: observations, analytical models and simulations (e.g.Schaerer et al. 2020;Popping et al. 2019;Carniani et al. 2020) suggest that there might be a weak evolution for the SFR-L [C ii] relation with respect to redshift.Our model suggests that there is a slight evolution, although less than for the other lines, as shown by the coefficient c 2 in Table C.1 and Equation C.1.The luminosity evolution may be related to the active star-formation processes that occur in starburst galaxies, as shown in Fig. 3, and in galaxies during "cosmic noon" at z ≈ 2. The regulation of these star-formation processes is reflected in the changes of the ISM phases (Fig. 4), which tend to be more stable in the case of Another possibility is to use L [O iii] as an SFR tracer, since the [O iii] lines tend to be equal to or brighter than L [C ii] in some redshift ranges (e.g.Arata et al. 2020;Carniani et al. 2020;Vallini et al. 2021;Bouwens et al. 2021;Katz et al. 2022).However, our L [O iii] predictions at z = 0 may be underestimated for lower SFR (see Sect. 3.4.2),which may explain the decrease in the median value of the L/SFR ratio for these lines.It is also possible to use L [O iii] together with L [O i] at 63 µm as an SFR tracer (Mordini et al. 2021) to balance the neutral and ionised components of the ISM in these lines.Nevertheless, this use would require having access to both luminosities (L [O iii] and L [O i] ) at z < 4 to confirm the trends presented in Fig. 14.The [O i] 63 µm line alone could also be used to trace the SFR, as this line has similar properties to the [C ii] line (Katz et al. 2022).For example, in the predictions of Olsen et al. (2018), [O i] is brighter than [C ii] at z = 6, and therefore is more easily observable.The trend of [O i] being brighter than [C ii] is also predicted in this work, confirming the results from Olsen et al. (2018).Therefore, L [C ii] may be the best SFR tracer across the entire redshift range of this work (z = 0-6), but at high-z other FIR lines such as [O iii] and [O i] are also very useful.

Diagnostic diagrams using FIR lines
We now examine our predictions of FIR emission line strengths from the EAGLE simulations in the context of diagnostic diagrams.Typically, diagnostic diagrams use emission line ratios that reflect the physical conditions of the ISM.In this work, we focus on two diagnostic diagrams: one that normalises the emission line luminosity with SFR for [C ii] and [O iii] at 88 µm (Harikane et al. 2020), and the other that uses the ratios between lines at 88 µm, 205 µm and 63 µm, respectively.With these two diagnostic diagrams we investigate whether these ratios trace physical quantities related to the ISM, such as radiation field and density.
Other line ratios in the FIR are also of great interest for different types of studies (e.g.Sect.3.3 of Cormier et al. 2015).We therefore make our model predictions for different emission lines of these simulated EAGLE galaxies publicly available, as described in Appendix B. Similar FIR diagnostic diagrams, such as those presented by De Breuck et al. ( 2019) and Li et al. (2020) are not discussed in this work, but plots are provided as supplementary material3 .

Comparison with observations
We compare our model predictions with the observational dataset of Appendix A and the results from the latest version of the SIGAME framework presented by Olsen et al. (2021).The latest version of the SIGAME framework uses the SIMBA simulations (Davé et al. 2019) with the SKIRT (Camps & Baes 2020) radiative transfer code and Cloudy (Ferland et al. 2017) to predict line luminosities, similar to this work.The SIGAME predictions come from galaxies in two simulated boxes of 25 cMpc and 100 cMpc volumes in the local Universe, similar to the EAGLE box sizes used in this work.Olsen et al. (2021) find that their L [C ii] predictions appear to be an extension to higher SFRs of the model predictions presented in Paper I. However, SIGAME tends to have higher line luminosities relative to the SFR-L [C ii] relationship (their Fig. 11).In addition, their 25 cMpc box returns even higher line luminosities up to ∼0.5-1.5 dex of the 100 cMpc box.
In Fig. 15, we compare the L [O iii] /SFR and L [C ii] /SFR ratios of our model predictions with the observational data and the SIGAME predictions.For the local Universe, our predictions share a similar range of values with the observational data, which tends to be in the range between 6.5 < log(L [C ii] /SFR) < 7.6 and 5.8 < log(L [O iii] /SFR) < 7.6.However, most of the simulated galaxies tend to be above (log(L [O iii] /SFR) > 6.8) or below (log(L [O iii] /SFR) < 6.0) the observational data.In contrast, SIGAME predictions tend to have very high [C ii] luminosities, which shifts most of the SIGAME simulated galaxies to log(L [C ii] /SFR) > 7.5 values, with the SIGAME values peaking an order of magnitude higher than the observed galaxies and our simulated ratios.This difference is expected from the comparisons presented by Olsen et al. (2021).As noted above, the SIGAME predictions of L [C ii] /SFR are higher than those described in Paper I, which are similar to those in this work.The SIGAME L [O iii] /SFR values seem to be similar to our predictions, with values between 6.0 < log(L [O iii] /SFR) < 8.
At z > 0, there are very few observed galaxies with which we can compare our results.At z = 2, we only have two measurements for one galaxy: H-ATLAS J091043.0-000322.These two measurements come from Herschel observations following different data reduction methods (i.e.pipeline versions) in Zhang et al. (2018) and Lamarche et al. (2018).ALMA observations presented by Lamarche et al. (2018) show that the [C ii] luminosity is around half of the Herschel measurement (right data point on the panel), but this difference cannot be fully explained.Therefore, we assume that H-ATLAS J091043.0-000322lies somewhere between the two measurements presented in the z = 2 panel, agreeing with the predictions from Recal-L0025N0752.
At z = 3, we have three measurements for SDP 81 (Valtchanov et al. 2011;De Looze et al. 2014;Zhang et al. 2018) and one for HLock01 (Rigopoulou et al. 2018) Zhang et al. (2018), which is the leftmost point in the z = 3 panel.The main difference between the values of Valtchanov et al. (2011) and Zhang et al. (2018) is the different data reduction methods, and therefore we rely on the more recent results of Zhang et al. (2018).Both SDP 81 and HLock01 are close to the model predictions of Recal-L0025N0752.
At z = 4, we have measurements for two galaxies: SPT-S J041839-4751.8(De Breuck et al. 2019) and AzTEC 1 (Tadaki et al. 2019).Our results agree with the location of both objects, a bit to the left of the results from the z = 0 galaxies.
Finally, at z = 6 we have measurements for two galaxies: [DWV2017b] CFHQ J2100-1715 companion (Walter et al. 2018) and [MOK2016b] HSC J121137.10-011816.4 (Harikane et al. 2020).[DWV2017b] CFHQ J2100-1715 companion is in the same region where most of our model predictions are at z = 6.In contrast, [MOK2016b] HSC J121137.10-011816.4 falls in the same upper region as SPT-S J041839-4751.8 in the z = 4 panel.The latter two galaxies have SFRs around 100 M yr −1 and have L [O iii] higher than their respective companions in their panels.These high SFRs can explain their positions in the diagnostic diagrams.Harikane et al. (2020) presented two other galaxies in the same redshift range for the ratios presented in Fig. 15 besides [MOK2016b] HSC J121137.10-011816.4.Those two galaxies have been removed from the comparison because they have been identified as QSOs.
For the other diagnostic diagram, in Fig. 16, we compare the [C ii]/[O iii] and [N ii]/[O i] ratios of our predictions with the observational data and SIGAME predictions.For the local Universe (z = 0 panel), we observe that our predictions, observational data and SIGAME predictions tend to be grouped in the region between −0.5 The SIGAME predictions match most of the observational data, especially for the 100 cMpc box, although they do not reach the values close to log( The difference between the SIGAME 100 cMpc and 25 cMpc boxes comes from the high L [C ii] that galaxies in the latter box can have, something that we also commented on at the beginning of this section.Our model behaves similarly to Fig. 15, where most of the observations lie between the two concentration regions for Recal-L0025N0752, and the Ref-L100N1504 agree with the observations.Our estimates with a simple model of the ISM are in the same parameter space as observations.Unfortunately, a completely fair comparison cannot be made because of selection bias in both the observational and simulated galaxies. Unfortunately, as we are using four FIR emission lines, there are few observational data points to which we can compare at higher redshifts.The only galaxy with all four FIR lines is SPT-S J041839-4751.8 at z = 4 (De Breuck et al. 2019).Our model matches the position of this galaxy in this diagram, similar to the results presented in Fig. 15.The location of this galaxy and most of our predictions at high-z in Fig. 16 seem to coincide with some of the observational data at z = 0, which may imply that some physical parameters (e.g.sSFR, metallicity and/or density) in these galaxies will be similar at different redshifts.
The comparison presented in Figs. 15 and 16 show that our model predictions largely match the parameter space of the observational data in the diagnostic diagrams.We conclude that our model predictions can be used to interpret the physical parameters of observed galaxies.

Physical parameters in diagnostic diagrams
Now, we use the diagnostic diagrams to infer the ISM physical conditions in galaxies.We use the simulation data of all modelled galaxies and compare the sensitivity of the line luminosityto-SFR ratio to eight physical parameters as estimated in Paper I: gas mass (M gas ), stellar mass (M * ), metallicity (Z/Z ), specific star formation rate (sSFR), interstellar radiation field (ISRF), total hydrogen number density in the neutral clouds (n(H) cloud ), external pressure (P ext ) and radius of the neutral clouds (R cloud ).We divide these physical parameters into seven ranges to compare them with the observational dataset of Appendix A and, as a reference, our mock data from Recal-L0025N0752 at z = 0.
We begin by presenting the physical parameters in the 17.We note that the impact of almost all the physical parameters is perpendicular to the observational sample.This effect arises because our mock data also tend to be perpendicular to the observational data, especially at z = 0. From these physical parameters, we see that most of the predicted galaxies at the upper-left boundary of the observational contour tend to have higher M gas , M * , sSFR, ISRF, n(H) cloud ) and P ext , and lower R cloud .In addition, low values of M * do not reach the bottom-right limit of the observational contour, while the metallicity spans all over the observational contours and does not have a clear trend.Harikane et al. (2020) used the /SFR diagnostic diagram to explain the physical properties of galaxies at z = 6-9 compared with the local Universe.Using simple Cloudy grids, they found that one of the probable reasons for the location of some of their galaxies in the upper-right region of the diagnostic diagrams was a high ionisation parameter, which is proportional to the ISRF.Their result is similar to what we find from Fig. 17; however, it is also important to compare the diagnostic diagram with our metallicity and density results.In terms of metallicity, Harikane et al. (2020) find that L [O iii] /SFR decreases with metallicity while L [C ii] /SFR does not change.We find the same result for L [C ii] /SFR, but the range of change of the L [O iii] /SFR ratio is not directly dependent on metallicity.In terms of density, Harikane et al. (2020) find that both ratios decrease with density, which we also find for our predictions of L [C ii] /SFR -but not for L [O iii] /SFR.A simple reason for this discrepancy is  that L [C ii] /SFR depends mainly on density while L [O iii] /SFR is more dependent on other physical parameters: namely metallicity, sSFR and gas mass.Figs. 15 and 17 highlight the importance of the [C ii] and [O iii] emission lines in understanding the physical conditions of gas of galaxies, especially by using their ratio, as other recent studies have done (e.g.Harikane et al. 2020;Arata et al. 2020;Carniani et al. 2020;Vallini et al. 2021;Bouwens et al. 2021) and as we show below.We check its correlation with other FIR lines that trace the neutral and ionised gas components, such as the ratios in Fig. 18.In this diagnostic diagram, we note how the physical parameters cross the observational data contours in different ways.Interestingly, the physical parameters in the bottom row panels show a "spoon-like" shape showing that most of the observed galaxies have low to moderate values for the ISRF, n(H) cloud ), and P ext , and moderate to high values for the R cloud .The opposite occurs in the region with low [C ii]/[O iii], (i.e.high ISRF, n(H) cloud and P ext , and low R cloud ), which coincides with our predictions for high-z galaxies (z > 3, in Fig. 16).In terms of mass, M gas is more extended but follows a similar trend to the parameters mentioned above, while M * does not show a clear trend.The most intriguing trends in the diagram come from metallicity and sSFR.In terms of metallicity, the [C ii]/[O iii] ratio seems to be a good tracer for metallicities close to solar, while the [N ii]/[O i] ratio is a good tracer for metallicities below log(Z/Z ) 0.5, which agrees with some results for high-z galaxies (Arata et al. 2020).In terms of sSFR, both ratios do a good job of separating high and low sSFR of galaxies in a zigzag pattern across the observational sample region.This supports the idea that the [C ii]/[O iii] ratio can be used for starbursting systems (Vallini et al. 2021).
We have shown that both diagnostic diagrams track the behaviour of the physical parameters presented in the simulated galaxies using the luminosities of the main FIR lines.Our model agrees with the observational data in most of the parameter space, and in some cases with other simulations (e.g.SIGAME).Therefore, it is reasonable to expect that the physical model could be used in inferring the physical parameters in FIR diagnostic diagrams to trace different physical parameters.
Comparing our model predictions with new observational data will give an idea of what kind of physical parameters are expected in those galaxies.At the same time, our model can constrain the expected luminosities of FIR lines when no other measurement is available.However, we emphasise that the modelled FIR lines can also be used to study other types of problems.For example, ratios such as [N ii]/[C ii] or other diagnostic diagrams can be used to characterise different types of galaxies (e.g.luminous infrared galaxies (ULIRGS, where L IR > 10 12 L )) in the local Universe (e.g.Farrah et al. 2013) and at high-z (e.g.Cunningham et al. 2020).Another example is that the ratio between the two [O lines can be used to test the efficiency of the black hole feedback in galaxies, as was done with the Il-lustrisTNG simulations (Inoue et al. 2021).Similarly, ratios like can be used as a metallicity indicator (e.g.Nagao et al. 2011;Rigopoulou et al. 2018;Fernández-Ontiveros et al. 2021), which will be explored in future studies.

Model assumptions
Although our model emission-line luminosities and ratios are in good agreement with observational results and other simulations, our physically motivated model is based on a number of assumptions.The effect of some of these assumptions is generally not visible in these types of predictions due to uncertainty in the observations and our current understanding of some of the physical processes involved.However, these assumptions may be important for models that try to predict the FIR line emission of galaxies, especially at high-z.In this section we highlight the most important assumptions.
For neutral clouds and Hii regions, we assume static spherical geometries to describe the densities and temperatures within those environments.This assumption is made for simplicity, although we know that these structures may not be spherical.Physical processes such as radiation destroy spherical geometries (Deharveng et al. 2010;Peñaloza et al. 2018), which may lead to rough or incorrect line luminosity predictions (Decataldo et al. 2020).However, approximations using mass distributions (Eq.15 in Paper I) may smooth out these differences, since cloud masses follow scaling relations that seem to be valid for observations at different redshifts (Dessauges-Zavadsky et al. 2019).
A problem with our modelled line luminosities from Hii regions is that we assume a fixed density (∼30 cm −3 ).We can in-  crease the luminosity by increasing this density and vice-versa.The use of different densities could be important when comparing the contributions of DIG and Hii regions, especially in lines such as [O iii] and [N iii] (see Figs. 11 and 13).This DIG-Hii region balance is still unclear in ionised emission lines, and although some estimates exist from optical wavelengths (e.g.Poetrodjojo et al. 2019), a change in this balance can lead to different metallicities, which may affect high-z studies (Sanders et al. 2017).In this work, we calibrate the DIG with observational data at z = 0, which may also bias the balance between these two ISM phases.Fortunately, the results presented in this work show that these assumptions seem to be in agreement with the observations, which may represent a likely first step in understanding the DIG-Hii region balance.Finally, our predictions depend on Cloudy lookup tables, which can give different emissivities depending on the assumed initial abundances or dust configurations (Ploeckinger & Schaye 2020).Different photoionisation models could lead to different interpretations of physical parameters coming from line ratios, such as metallicity or ionisation parameter (e.g.Ji & Yan 2022).Furthermore the intrinsic thermodynamics may not be the same in terms of cooling and heating functions in cosmological simulations and photoionisation models like Cloudy (Robinson et al. 2021).Therefore, some care must be taken when interpreting the line luminosities predicted by our model.

The use of EAGLE
In this work we use the EAGLE simulation as a proxy of what we would expect to see in the Universe.However, by using a cosmological simulation like EAGLE, our line luminosity predictions of the ISM model are expected to inherit the limitations of the simulation.An example of these limitations is the lack of starburst-like galaxies within EAGLE (Wang et al. 2019b).As we discussed in Paper I, this restricts our comparison at z = 0 mainly to SFR below 10 M yr −1 but also limits the comparisons at other redshifts (Katsianis et al. 2017).To compare our predictions with high SFR observations, we extrapolate the linear relations in the range of −3.5 < log(SFR) < 3.5.However, care must be taken when using this extrapolation.
Another important physical property within EAGLE that affects our predictions is the gas metallicity, which is usually studied through the gas-phase mass-metallicity relation (MZR).Bellstedt et al. (2021) show that the MZR in EAGLE galaxies, as measured by Zenocratti et al. (2020) at z = 0, does not behave similarly to other cosmological simulations or semi-analytical models.As shown in Fig. 4 of Bellstedt et al. (2021), the metallicities in EAGLE have values around 12+log(O/H)∼9 for stellar masses between 9 < log(M * [M ]) < 11, which is high compared to observations with metallicities going from 12 + log(O/H)∼8.5 to 12 + log(O/H)∼9.2 in the same stellar mass range.This may affect the metallicity that can be recovered from FIR lines if only z = 0 is used.However, MZR in EAGLE depends on resolution and box-size used due to the assumed AGN and star formation feedback processes (De Rossi et al. 2017).In Figs. 17 and 18, we present the predicted luminosities in different boxes and redshifts studied in this work and find that some of the FIR line ratios can be useful to infer the metallicity.Therefore, it will be important to compare the FIR line predictions that trace metallicity with observations in the future more consistently.
Although we base our predictions on EAGLE, we expect that similar physical models and/or cosmological simulations can be used to understand the gas properties with FIR emission lines.For example, Olsen et al. (2021) show that these kinds of predictions can be obtained with a different gas fragmentation scheme using SIGAME, which are similar to observations and our predictions.However, most of the galaxies in SIGAME have higher SFR than those studied in this work.This difference in the sample of simulated galaxies comes from choosing a different simulation (SIMBA instead of EAGLE), which allows for the forma-tion of starburst-like systems.Therefore, although some physical assumptions are different for each model and are limited by the simulation used (e.g.Vallini et al. 2015;Lagache et al. 2018;Popping et al. 2019;Pallottini et al. 2019;Leung et al. 2020), we expect predictions to behave similarly to those presented in this work.In the future, we will compare results from different simulations, such as SIMBA and IllustrisTNG in an efficient way to reduce the bias that the initial assumptions of the simulations may introduce.

Observational data from samples
For our comparison between observations and our model line luminosities, we have collected a heterogeneous sample of observed galaxies with FIR emission line information.We have transformed the luminosities of the lines to the same reference cosmology (Planck Collaboration et al. 2014) and we have estimated the SFR in most of the cases following the L IR -SFR relation described by Kennicutt & Evans (2012).However, there are some other important issues that may affect this heterogeneous sample.
The most important issue may come from the initial mass function (IMF) assumed to estimate the SFRs.Unfortunately, IMF information is not always present in papers.Assuming a star-formation law that takes into account a standard IMF is a possible solution, as we have tried in this work.However, this may add additional uncertainties.For example, EAGLE uses an instantaneous SFR, which is different from the SFR averaged over the last 10 or 100 Myr typically used in observations.Another problem is that some of the SFR (or L IR ) may come from spectral energy distribution (SED) models rather than empirical laws using FIR photometric bands, like IRAS, PACS or SPIRE.In addition, the IR luminosity definition can have different flavours (e.g.L IR , L FIR or L TIR ) that use different wavelengths ranges to estimate luminosities, adding to the spread of SFR estimates (Förster Schreiber & Wuyts 2020) and leading to differences of factors between 1.1 and 1.7 (De Looze et al. 2014;Lagache et al. 2018).Another possibility is to use other SFR tracers, such as Hα-based, UV-based, or radio-based SFRs as shown in many studies (e.g.Wang et al. 2016Wang et al. , 2019a)).However, in some cases, such as in high-z galaxies, it can be difficult to obtain such SFRs estimates and other lines like [C ii] come into play as a SFR tracer (Le Fèvre et al. 2020), as we have also shown in this work.
In addition to theoretical considerations, some of the FIR line measurements have systematic uncertainties.For example, calibration pipelines or the use of different instruments may affect the comparison, as shown by Lamarche et al. (2018) in the case of H-ATLAS J091043.0-000322,discussed in Sect. 4.1. Furthermore, apertures in the observational sample may affect the analysis of the fixed-size aperture we selected for the simulated galaxies (30 pkpc).For example, the balance between DIG and Hii regions as main contributors to the ionised lines could depend on the selected aperture, as indicated by Mannucci et al. (2021).Finally, gravitational lensing can introduce a large uncertainty in the luminosities, which can be reduced by a factor of 30 to 40 in some galaxies when corrected by the magnification factor (e.g. in Eyelash and SPT-S J041839-4751.8,Zhang et al. 2018;Rizzo et al. 2020, respectively).In addition, these magnification factors may change depending on the observed region of the galaxy (Lamarche et al. 2018).In this work, we highlight those galaxies to warn about these potential complications, but in general, they show good agreement with our predictions.

Summary and Conclusions
We have modelled FIR emission lines by post-processing EA-GLE cosmological simulations using Cloudy.We have predicted the luminosities of the eight most important lines in the FIR up to z = 6, using a physically motivated model that traces four different ISM phases: dense molecular gas, neutral atomic gas, diffuse ionised gas and Hii regions.We have also collected a sample of observed galaxies from the literature with which to compare with our predictions.Our main conclusions are as follows.
1. Predictions from our model replicate observed galaxies in the SFR-FIR line luminosity relationship to an average degree of 0.5 dex over the range z = 0-6, which is reasonable considering the observational measurements errors.We also compare with different models showing similar level of agreement.We model the SFR-FIR line luminosity relationship for each of the eight lines with a linear relation, each of which shows a slight evolution with redshift.2. We have presented the expected contributions of each ISM phase to each FIR line.These contributions change as a function of SFR.For the [C ii] 158 µm line, the main contributor is the neutral atomic gas, with considerable contributions from Hii regions at z = 1-4 and the DIG at z < 2, which may be related to metallicity.In the future we expect to use this model and its predictions to understand the effect that AGN can have on the ISM, as well as the physical parameters traced by these lines and their ratios.We make our model predictions and collected observational sample publicly available to allow potential users to compare with their work and/or interpret new observations.We envisage that our predictions will also be useful in planning for future FIR missions.
Acknowledgements.The authors thank Karen Pardos Olsen for providing the data from the SIGAME framework.This research has made use of the 3Mdb database 4 (Morisset et al. 2015) from which initial experiments to model Hii regions were performed.This research made use of Astropy, 5 a communitydeveloped core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018).This research made use of several Python packages, among them: numpy (Harris et al. 2020), pandas (Wes McKinney 2010) and matplotlib (Hunter 2007).This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.This research has made use of NASA's Astrophysics Data System Bibliographic Services.The EAGLE simulations were performed using the DiRAC-2 facility at Durham, managed by the ICC, and the PRACE facility Curie based in France at TGCC, CEA, Bruyères-le-Châtel.We would like to thank the Center for Information Technology of the University of Groningen for their support and for providing access to the Peregrine high performance computing cluster.

Appendix A: Observational sample
We have collected measurements from the literature of the FIR emission lines predicted in this work.The observational sample is a heterogeneous selection of galaxies that covers the redshift range between 0 and 6.5.We present references of the works used in this sample of galaxies together with the number of measurements available per line in Table A.1.When possible we recalculated the luminosities and SFRs onto the ΛCDM cosmology used in this work (Planck Collaboration et al. 2014).When the cosmology is not explicitly mentioned, we assume it is the same as used in this work, so no corrections are applied.We flag those measurements where a magnification factor is involved in estimating a luminosity due to gravitational lensing.We estimate SFR in most of the literature samples from the infrared luminosity (L IR ) using the relation of Kennicutt & Evans (2012), which assumes a similar IMF (Kroupa) to the one used in EAGLE (Chabrier).Unless stated otherwise, we use the same infrared luminosities as published in the respective works.In cases where L IR is unavailable or is unreliable due to measurement error, we use SFR estimates from other works.Additionally, we remove strong AGN galaxies (i.e.QSO and Blazar).This observational sample is available as a supplementary material in a Zenodo repository, at https://doi.org/10.5281/zenodo.6576202.

Fig. 1 .
Fig. 1.Flowchart of the sub-grid procedures applied to the SPH to simulate FIR line emission in post-processing.This flowchart is similar to the one presented in Paper I with the main difference being the added Hii regions as a new ISM phase.The dashed lines connect the gas and star environments to the ingredients of the model.

Fig. 2 .Fig. 3 .
Fig. 2. SFR-L [C ii] relation for all redshift slices used in this work.We compare the obtained relations from the EAGLE simulations Ref-L100N1504and Recal-L0025N0752 with predictions from other simulations(Olsen et al. 2015(Olsen et al. , 2018;;Katz et al. 2022), semi-analytic models(Lagache et al. 2018;Popping et al. 2019) and linear relations derived from observations(De Looze et al. 2014;Schaerer et al. 2020).Linear relations inferred from our models are shown as black solid lines over the dynamic range covered by the simulations and extrapolated to lower and higher SFRs as black dotted lines, with the grey shaded area representing the 1 σ error.Collections of observational data (Appendix A) are plotted as grey squares for detections and as grey triangles for upper limits.For lensed galaxies, the markers are plotted in and the upper limits that affect both the SFR and luminosity are plotted as diamonds.

Fig. 4 .
Fig. 4. Contribution from the different ISM phases for the [C ii] emission line in Recal-L0025N0752.The regions define the DIG (blue), neutral atomic gas (orange), dense molecular gas (green) and Hii regions (red).

Fig. 5 .
Fig. 5. Contribution of the ionised gas phase to [C ii].We show the median values of the different redshifts (solid lines) and their 25th and 75th percentiles (dotted lines) with respect to the fit to the observational data made by Cormier et al. (2019, black dashed line).In the left panel, we show the predictions of Ref-L100N1504 while in the right panel we show the predictions of Recal-L0025N0752.We only show the bins with more than 5% of the total sample of galaxies for the respective simulation.
very similar to the estimations of the Ref-L100N1504 simulation.However in an updated model, Olsen et al. (2021) use the SIMBA cosmological simulation at z = 0, and their predicted [O i] 63 µm luminosities are 1.2 dex above the De Looze et al. (2014) relation shown in Fig. 8.This contrasts with our model, which exhibits better agreement with the linear fits of De Looze et al. (2014) and Mordini et al. (

Fig. 11 .
Fig. 10.As Fig. 2 for the [O iii] 88 µm line.We compare the obtained relations from the EAGLE simulations with relations derived from simulations (Olsen et al. 2018; Katz et al. 2022; Kannan et al. 2021) and observations (De Looze et al. 2014; Harikane et al. 2020; Mordini et al. 2021).Linear relations inferred from our models using only Hii regions are shown as indigo solid lines over the dynamic range covered by the simulations and extrapolated to lower and higher SFRs as indigo dotted lines, with the shaded area representing the 1σ error (see Appendix C).

Fig. 14 .
Fig. 14.Evolution of the line luminosity-SFR ratio with redshift for the main FIR lines modelled in this work.We show the median values from Recal-L0025N0752 (dotted lines) and Ref-L100N1504 (solid lines).The shaded regions correspond to the 25th and 75th percentiles.

Fig. 15 .
Fig. 15.Diagnostic diagram for the L [O iii] /SFR and L [C ii]/SFR ratios, similar to that presented byHarikane et al. (2020).Cyan and olive contours show the model predictions from Ref-L100N1504 (cyan) and Recal-L0025N0752 (olive).We compare with observational data (black squares) and SIGAME predictions(Olsen et al. 2021) for the local Universe in 25 Mpc and 100 Mpc simulation boxes (purple and chocolate contours).All panels with redshifts above zero show the z = 0 Recal-L0025N0752 estimations as grey dashed contours.

Fig. 16 .
Fig. 16.Diagnostic diagram using four different FIR emission lines comparing the [C ii]/[O iii] ratio against the [N ii]/[O i] ratio.Colour-codes are the sames as in Fig. 15.

Fig. 17 .
Fig. 17.Physical parameters in the L [O iii]88 /SFR-L [C ii] /SFR diagnostic plot (see Fig.15).All panels show the z = 0 Recal-L0025N0752 model predictions as grey dashed contours and the observational data as green dashed contours.In each panel, colour-coded rectangles represent the 25th and 75th percentiles for a given parameter from low (blue) to high (red) values.The solid black lines connect the median values of each rectangle.The scales for all physical parameters are logarithmic.

Fig. 18 .
Fig. 18.Physical parameters in the diagnostic plot comparing the [C ii]/[O iii] ratio against the [N ii]/[O i] ratio.Colour-codes are the sames as in Fig. 17.
Fig. 18.Physical parameters in the diagnostic plot comparing the [C ii]/[O iii] ratio against the [N ii]/[O i] ratio.Colour-codes are the sames as in Fig. 17.
Tadaki et al. (2019): We assume the cosmology used is a ΛCDM cosmology with Ω = 0.27 and H 0 = 70 km s −1 Mpc −1 .Neeleman et al. (2019): For [NBW2019] J0842+1218C2 we assume that the upper limit for SFR is 100 M yr −1 .Lee et al. (2019): We use the IR luminosity values and [C ii] measurements from Iono et al. (2006).Hashimoto et al. (2019):The reported values for the SFR are assumed to be upper limits, as those come from QSOs.Béthermin et al. ( ; ii) ionised nitrogen ([N ii]) at which traces the ionised medium from DIG and Hii regions (in the local Universe, Cormier et al. 2012; Goldsmith et al. 2015; Zhao et al. 2016a; Croxall et al. 2017; Langer et al. 2021, and at high-z, e.g.Pavesi et al. 2016); iii) [O iii] at 52 and 88 µm, which traces Hii regions around young stars (at high-z Ferkinhoff et al. 2010; Inoue et al. 2014); and iv) [N iii] at 57 µm, which also traces Hii regions

Table 3 .
Sampling of the properties of Hii regions in the Cloudy grid used in this work.The resulting number of grid points is 3600 per redshift.
(De Looze et al. 2014;Harikane et al. 2020;Mordini et al. 2021in previous sections, we focus here on the results of only one of the [O iii] lines, the [O iii] 88 µm line, as we have more observational data for this line.In Fig.10we present the [O iii] 88 µm luminosities as a function of SFR.We compare the predictions from the EA-GLE simulations Ref-L100N1504 and Recal-L0025N0752 with observations of individual galaxies and linear relationships between [O iii] and SFR derived from observations(De Looze et al. 2014;Harikane et al. 2020;Mordini et al. 2021) and numerical models Ferkinhoff et al. 2010;De Looze et al. 2014;Inoue et al. 2014;Harikane et al. 2020;Yang & Lidz 2020;Yang et al. 2021a)f ionised gas in the FIR.These lines may be used as SFR tracers as they come mainly from young stars (e.gFerkinhoff et al. 2010;De Looze et al. 2014;Inoue et al. 2014;Harikane et al. 2020;Yang & Lidz 2020;Yang et al. 2021a).The [O iii] 88 µm line has become very important as it may be brighter than the [C ii] 158 µm line in galaxies close to the reionisation epoch (z 7, . For SDP 81, two of the measurements coincide (De Looze et al. 2014 use the values from Valtchanov et al. 2011): the empty squares with log(L [C ii] /SFR)∼7.6.The other measurement for SDP 81 comes from For the [N ii] lines at 122 and 205 µm, the DIG contributes more than Hii regions in the local Universe, but the opposite is true at high-z, where Hii regions seem to dominate over the DIG.For the [O i] lines at 63 and 145 µm, the contribution of dense molecular gas is important in the local Universe.However, the atomic gas is dominant at high-z.Finally, for the [O iii] lines at 52 and 88 µm and [N iii] at 57 µm, we show that Hii regions dominate, with important contribution from the DIG at low SFR in the local Universe.3. Our predictions indicate that [C ii] may not be a good SFR tracer for starburst galaxies, since the [C ii]/SFR ratio seems to decrease as a function of the offset from the star-forming main-sequence.However, compared to the other FIR lines, [C ii] seems to be the best SFR tracer due to its weak redshift evolution.[O iii] and [O i] may also be good SFR tracers.Nonetheless, our predictions of [O iii] at z = 0 may be underestimated, and more observations of [O i] are necessary at z < 4 to confirm our predictions.4. We compare our predictions in two diagnostic diagrams, and we find reasonable agreement with observations.We compare L [O iii] /SFR and L [C ii] /SFR ratios and find that mock galaxies at high-z tend to have higher L [O iii] /SFR ratios and slightly lower L [C ii] /SFR ratios than galaxies in the local Universe.We also compare the [C ii]/[O iii] and [N ii]/[O i] ratios and find that mock galaxies at high-z tend to have lower [C ii]/[O iii] and [N ii]/[O i] ratios than galaxies in the local Universe.5. Finally, we have examined the impact of physical parameters on these diagnostic diagrams.When we compare physical parameters to line luminosities, we find that L [C ii] /SFR and L [O iii] /SFR ratios trace hydrogen density and ISRF well in the mock galaxies.However, these ratios are not good metallicity tracers, because L [O iii] /SFR does not evolve linearly with metallicity and L [C ii] /SFR does not change with metallicity.Furthermore, we find that [C ii]/[O iii] and [N ii]/[O i] ratios can be good metallicity and sSFR tracers.For example, [C ii]/[O iii] can trace metallicities close to solar and [N ii]/[O i] below solar.On the other hand, we can identify systems with different sSFRs by means of both [C ii]/[O iii] and [N ii]/[O i] ratios, which can be very useful for improving calibrations of [C ii] as a SFR tracer.
for Mrk 33 and discard the information from the LMC and SMC.For 2MASX J12390403+3920437 we use the SFR reported by Duarte Puertas et al. (2017).De Breuck et al. (2019): We use the magnification factor of 32.7 by Spilker et al. (2016).