Infrared view of the multiphase ISM in NGC 253 II. Modelling the ionised and neutral atomic gas

Context. Multi-wavelength studies of galaxies and galactic nuclei allow us to build a relatively more complete picture of the interstellar medium (ISM), especially in the dusty regions of starburst galaxies. An understanding of the physical processes in nearby galaxies can assist in the study of more distant sources at higher redshifts, which cannot be resolved. Aims. We aimed to use observations presented in the first part of this series of papers to model the physical conditions of the ISM in the nuclear region of NGC253, in order to obtain primary parameters such as gas densities and metallicities. From the created model we further calculated secondary parameters such as gas masses of the di ff erent phases, and estimated the fraction of [C ii ]158 µ m from the di ff erent phases, which allowed us to probe the nuclear star-formation rate. Methods. To compare theory with our observations we used MULTIGRIS , a probabilistic tool that determines probabilities for certain ISM parameters from a grid of Cloudy models together with a set of spectroscopic lines. Results. We find that the hypothetical active galactic nucleus within NGC253 has only a minor impact compared to the starburst on the heating of the ISM as probed by the observed lines. We characterise the ISM and obtain parameters such as a solar metallicity, a mean density of ∼ 230cm − 3 , an ionisation parameter of log U ≈ − 3, and an age of the nuclear cluster of ∼ 2Myr. Furthermore, we estimate the masses of the ionised (3 . 8 × 10 6 M ⊙ ), neutral atomic (9 . 1 × 10 6 M ⊙ ), and molecular (2 . 0 × 10 8 M ⊙ ) gas phases as well as the dust mass (1 . 8 × 10 6 M ⊙ ) in the nucleus of NGC253.


Introduction
In the extreme environments of galactic nuclei, the various heating and cooling mechanisms at work within the interstellar medium (ISM) are notoriously entangled.Consequently, the physical conditions of the ISM are difficult to constrain.In particular the heating mechanisms are strenuous to unravel, as they are typically and indirectly probed through specific corresponding cooling processes.Multi-wavelength observations and models, ideally with tracers sensitive to the different ISM phases, are therefore mandatory to link heating and cooling processes which regulate matter cycle and star formation.Nearby galaxies such as NGC 253, with a distance of 3.5 Mpc (Rekola et al. 2005), are ideal laboratories to study these effects.Our understanding of the physics in nearby galaxies may then be useful for more distant sources, where spatial resolution is far more limited.Models of the ISM are becoming ever more complex, able to account for an increasing number of mechanisms, encompassing not only stellar emission, but also the effects of active galactic nuclei (AGN), shock, cosmic rays etc.Furthermore, an increasing number of ISM cooling diagnostics can be considered by the models, in-cluding many infrared observables, such as atomic fine-structure emission lines.
In this study we used observations of the nuclear region of NGC 253, presented in Beck et al. (2022) -herafter B22 -from infrared telescopes to model the ionised and neutral atomic gas.Although being one of the archetypical starburst galaxies, the major excitation conditions in the centre of this galaxy are not yet well understood (e.g.Vogler & Pietsch 1999;Günthardt et al. 2015): The nuclear starburst likely plays a crucial role in the gas excitation, however, the central supermassive black hole (BH) and associated AGN may have an impact as well (Fernández-Ontiveros et al. 2009).When constraining the impacts of heating on the ISM, it is also necessary to constrain a variety of relevant physical parameters, such as the metallicity, optical depth, or gas density.For instance, the latter varies widely within the literature for NGC 253, depending on the regions observed and the tracers used (e.g.Puche & Carignan 1991;Wall et al. 1991;Carral et al. 1994;Engelbracht et al. 1998).In B22, we used a subset of the observed emission lines and analytical approach, from which we obtained solar metallicity, an optical depth of 4.35 mag, and a mean density of ∼ 150 cm −3 for the ionised gas within the central ∼ 100 pc.From the more complete picture built by the multiwavelength study it is also possible to derive simultaneously and self-consistently masses of the ionised and neutral atomic hydrogen as well as the dust mass, which are the goals of this work.As of this study, these have only been obtained on much larger scales (Melo et al. 2002;Weiß et al. 2008).
The [C ii] 158 µm line has been recently proposed as a powerful probe to determine not only the local star-formation rate (SFR, Herrera-Camus et al. 2018), but also for the CO-dark gas (Madden et al. 2020).However, the validity of [C ii] 158 µm as a tracer is hindered by its ubiquitous nature, potentially originating from ionised, atomic and molecular phases of the ISM.An objective in this study is to determine masses of the different phases, and calculate the origin of the [C ii] 158 µm from the different phases.Furthermore, from our models we estimate and compare SFRs from different tracers.
This paper is organised as follows: In Sect. 2 we explain the model grid and the code used to infer ISM parameters from this grid.Section 3 continues with creating a model for the ionised gas, by using emission lines that originate primarily in H ii regions.This model serves as a starting point to create a model for the ionised and neutral atomic gas in Sect. 4. The discussion in Sect. 5 presents parameters that can be inferred from the model explained in Sect. 4.

Observational constraints
In B22 we presented far-infrared observations from SOFIA/FIFI-LS, complemented by mid-infrared and farinfrared archival data from Herschel/PACS, Herschel/SPIRE, and Spitzer/IRS.Apertures are shown in Fig. 1.We obtained line fluxes and corresponding line flux errors of 30 emission lines, covering different species (C, N, Ne, O, Si, . . . ) and ionisation states of the nuclear region of NGC 253.The observed emission lines are sensitive to various density regimes, due to their wide range of critical densities.For instance the critical densities for [Ne ii] 13 µm and [N ii] 205 µm are 7 × 10 5 cm −3 and 45 cm −3 , respectively.The different ionisation states with ionisation potentials up to 97.12 eV to create Ne 4+ also show the variety of physical conditions in the centre of NGC 253.Owing to the large range of ionisation potentials and density regimes to excite these infrared emission lines, we are now able to draw a more complete picture of the embedded regions in the nucleus of NGC 253 and estimate the physical conditions there, by modelling all lines simultaneously.
In addition to the line fluxes listed in Table 3 of B22 we also used observations of the CO rotational transitions CO(1 − 0) to CO(3 − 2) from Israel et al. (1995) and CO(4 − 3) to CO(8 − 7) from Pérez-Beaupuits et al. (2018) to compare a posteriori.Higher-J CO emission lines are also available in Pérez-Beaupuits et al. (2018).However, these lines are not available in the model grid (Sect.2.2) and hence we could not use them as constraints nor for a comparison a posteriori.

Model grid
The grid of models used in this study is Star-Forming Galaxy with an X-ray component (SFGX, Ramambason et al. 2022) 2018).In addition to the star cluster, an X-ray source with inner disk temperature T X2 and luminosity L X is modelled by a multicolour black body (Mitsuda et al. 1984).In this study we refer to this combination of star cluster and X-ray source, which are assumed to be co-spatial, as cluster.A cluster in this context is described by a parameter set (L ⋆ , t, L X , T X ).The X-ray source in these Cloudy models is particularly important in our study, because of the presence of highly ionised emission lines such as [Ne v] 14 µm and [O iv] 26 µm.Ions like Ne 4+ and O 3+ with ionisation potentials > 54.9 eV (or λ = 22.6 nm) to create the ion, are unlikely to be created by stars and likely imply the presence of an AGN or a different source of X-ray emission.
An ISM component in the grid is characterised by the metallicity Z, ionisation parameter U, hydrogen volume density at the illuminated edge of the cloud n, and the depth.We emphasise that the metallicity refers to the O/H abundance.Other elements (such as N, Ne, C, Ne, Si) scale with the O/H abundance -see Ramambason et al. (2022) for the respective relations.Since the luminosity is already included in the cluster properties, the ionisation parameter is defined by the distance r between the cloud and the cluster.SFGX was created to investigate physical conditions in low-metallicity galaxies.Since the optical depth A V as a function of depth varies with metallicity (i.e. two clouds with the same physical size but different metallicities have different A V ), a more general definition of the cloud depth was defined by a "cut" parameter ξ as shown in Table 1.
Notes. (a) In some models -in particular ones with large X-ray luminosities or high CRIR that keep the ISM warm-the CO photodissociation front is not reached within the Cloudy simulation.
In the following a parameter set (Z, U, n, ξ) will be called a component.Note that in case of modelling more than one component, the metallicity is the same for all components.The SFGX grid contains a total of 32000 models, which are divided into 544000 sub-models by introducing 17 cut values between ξ = 0 and ξ = 4.We refer to Ramambason et al. (2022) for more details about the grid of models and an application to the Herschel Dwarf Galaxy Survey (Madden et al. 2013).

Inference method
We used MULTIGRIS for the modelling approach (Lebouteiller & Ramambason 2022).MULTIGRIS uses a Bayesian approach, designed to calculate posterior probability density distributions of parameters of a given grid of models, by using priors (e.g.line fluxes and errors) as constraints.Upper and lower limits for line fluxes can be used as well.
MULTIGRIS allows a linear combination of components, or continuous distributions for the parameters.This is particularly important, as the ISM in galaxies is not homogeneous on (kilo-)parsec-scales, but instead heterogeneous and a mix of several components (e.g.Lacy et al. 1982;Burton et al. 1990;Snow & McCall 2006;Cormier et al. 2019).By considering several components or using, for instance, power-laws for certain parameters, various geometries that are relevant to studies of complex objects like galaxies can be explored.Hence, we investigated two different model architectures: a one cluster and two component architecture (1C2S), and an architecture where U, n, ξ, and t are continuously distributed following a power-law (PLaw).Both architectures have been used to model the ISM in galaxies, see for instance Péquignot (2008); Polles et al. (2019) approaches with a discrete sampling and Baldwin et al. (1995); Richardson et al. (2016); Ramambason et al. (2023) for continuous distribution of ISM parameters.A power-law approach, for instance for the depth and density takes the porosity and clumpiness of the ISM into account, by enabling a larger number of low density (diffuse) clouds and only a few clouds of higher density and large depth, or vice versa.Each parameter that is described by a power-law introduces three variables: the slope α of the power-law, and the lower and upper boundary, between which the power-law is valid.
The choice for a probabilistic approach is justified by the complexity of the problem, for example: -The combination of two or more components (i.e.stellar populations or gas phases) yields a large number of free parameters making it difficult to gauge degenerate solutions.
-Upper limits or asymmetric uncertainties of line fluxes are difficult to handle with deterministic approaches.-Known parameters (e.g. the metallicity) and their uncertainty cannot be used as priors in a deterministic approach.
Probabilistic tools such as MULTIGRIS are able to overcome all the mentioned issues of deterministic techniques.For more details about the methodology, sampling techniques and applications of MULTIGRIS see Lebouteiller & Ramambason (2022).We note that in B22 we listed the obtained line fluxes and statistical uncertainties in Table 3.The uncertainties did not include systematical or calibration errors.We included these uncertainties by defining line sets for each instrument.Each line set includes the emission lines observed by one instrument.We assume that systematic and calibration uncertainties affect the observed emission lines of an instrument all in the same way.For each line set, we introduced a scaling parameter to account for these systematic and calibration uncertainties.This also allowed us to consider potential offsets due to the different spatial resolution and wavelength coverage and hence the different size of the nucleus.In particular we also divided the Spitzer/IRS emission lines in two observation sets since this instrument consists of two different modules with different aperture sizes, namely short-high and long-high modules.See the discussion in Sect.2.3.2 of B22 for details.

Model for the ionised gas
In this first step we investigated model results using the emission lines arising from H ii regions only, meaning all lines with an ionisation potential ⩾ 13.6 eV.We note that lines with a lower ionisation potential (for example [C ii] 158 µm and [Si ii] 35 µm) may also arise in H ii regions (e.g.Abel et al. 2005;Chevance et al. 2016).However, we did not use these lines as constraints in the first step, because of their ambiguous origin.But we let the model predict which fraction of [C ii] 158 µm comes from the ionised gas only.This is done by comparing the cumulative line flux of [C ii] 158 µm at the ionisation front and at the solution for the cut-value ξ.Lines with lower ionisation potentials will be included in Sect. 4. Also we note that the [Fe iii] 23 µm emission line is not included -due to the strong dependence of Fe emission lines on the dust-to-gas ratio, these emission lines will be handled separately (see Sect. 4.2).
Recently, Behrens et al. (2022) reported that the cosmic ray ionisation rate (CRIR) in selected clouds in the nuclear region of NGC 253 is about 10 4 × larger than the average Galactic value ζ 0 = 2 × 10 −16 s −1 (Indriolo et al. 2007).The value for the CRIR used in SFGX, however, is only 3ζ 0 (Ramambason et al. 2022).We will show in Sect.3.4 that cosmic rays only mildly affect the ionised gas and hence the effect on the model results are negligible.The effect of the increased CRIR on the neutral atomic and molecular gas, however, is significant and will be explored on a smaller sub-grid (see Sect. 4).

Systematic offsets within the dataset
In B22 we showed the data reduction and cross-calibration of observations of NGC 253 with SOFIA, Spitzer, and Herschel.
To account for potential offsets between the line sets (i.e. the emission lines as observed by one instrument) due to calibration or aperture effects, we performed MULTIGRIS runs where the different line sets were successively added (i.e.first run with IRS short-high (SH) and long-high (LH), second run adding the Notes. (a) Edge length of the rectangular aperture (Fig. 1). (b) Radius of the circular aperture.See B22 for the exact SOFIA apertures and their determination.For comparison, the effective radius of the nucleus in the near-infrared bulge is 9.1 ′′ or 150 pc (Iodice et al. 2014).
PACS lines etc.).MULTIGRIS determines the -small -systematic offsets between datasets so that observations match better with models.The LH module of Spitzer/IRS served as standard (i.e.offset LH ≡ 1).As the aperture size of this instrument matches the best with the angular size of the nucleus (6.68 ′′ ), we believe that this observation contains low background emission and we therefore expect the LH line set to be of the best quality.In this way, systematic offsets between line sets should be the dominant source of uncertainties.To simultaneously investigate the effect of different model architectures, we performed the same runs for a 1C2S and for a PLaw architecture.These architectures fulfil the request for more than one component as claimed in B22, but are still as simple as possible.Due to the limited number of constraints, a simple model with fewer parameters preferably avoids overfitting.Table 2 shows the obtained offsets between observation sets that will be further used in this work.Line fluxes shown in Table 3 of B22 were multiplied by these factors to account for the offsets introduced by calibration uncertainties and PSF effects between the instruments.The systematic offsets of all line sets are smaller than 1.This could be due to different beam sizes, meaning that the Herschel and SOFIA observations include more background emission because of their lower spatial resolution.This would also explain that the PACS, SPIRE, and FIFI-LS offsets are all in good agreement.The short-high observations of Spitzer/IRS are also < 1.Since this observation does not fully sample the nucleus due to a small aperture size, we performed a scaling correction of 1.55 based on the continua of SH and LH in B22.The lines may not arise from the exact same region as the thermal emission at ∼ 20 µm, which could lead to an over-correction of the line fluxes and thus a systematic offset < 1 in the SH emission lines.

Impacts of the model configuration
To evaluate the performance, we used the marginal likelihood as well as "pnσ" values, both of which are computed by MULTIGRIS.The marginal likelihood measures the probability that the model grid including any prior reproduces the observations and will reach a maximum for the most likely set of parameters under the assumption that the priors do not change.Too few parameters are not able to reproduce all the emission lines and will therefore have a low marginal likelihood, while too many parameters will cause overfitting, that is that the model contains more parameters than can be justified by the constraints.The pnσ values on the other hand describe how many draws of the posterior distribution agree within nσ of the constraints or observables.The higher the pnσ values, the better the model is able to reproduce the observations within their uncertainties.
Generally, both architectures show high pnσ values with p3σ = 99.19% and p3σ = 96.44%for the 1C2S and PLaw architecture, respectively.Moreover, both are able to reproduce all emission lines (with the exception of [O iv] 26 µm) within the uncertainties (see Fig. 2).The log marginal likelihood is better for the 1C2S architecture (−31.8 compared to −37.3 for the PLaw architecture).However, this does not mean that the 1C2S architecture is a more realistic or "better" one -see for instance Lebouteiller & Richardson (2023) for a discussion on model architectures.To conclude, the model architecture does not seem to play a significant role for the modelling of the ionised gas.

Results for the ionised gas
We started by narrowing down the parameter space of our model grid.This is necessary since a new model grid with a larger CRIR is required but the full set of Cloudy models would take a long time to run.To do so, we investigated parameters that can already be determined from the solution for the ionised gas.For instance, emission lines from species with high ionisation potentials such as [O iv] 26 µm (55 eV) and [Ne v] 14 µm (97 eV) are almost exclusively created by X-rays (Abel & Satyapal 2008).Hence, the parameters of the X-ray source are already constrained when using these emission lines.Another parameter that can already be determined is the metallicity Z, provided we assume that the metallicity does not vary significantly within the nuclear region.We further estimate the stellar luminosity L ⋆ from the ionised gas model, assuming that this parameter should not vary greatly once the neutral atomic gas lines are taken into account.We focus on the PLaw architecture in this section and will compare the results with those from the 1C2S architecture at the end of this section.Figure 3 shows the probability density functions (PDF) of the parameters determined in this section.Mejía-Narváez et al. (2020) interpret the PDFs, for instance for the metallicity, as a physical distribution within the ISM, however, we think that the PDF in our case is dominated by the uncertainty rather than a physical metallicity distribution function.Table 3 shows the resulting median and standard deviations for all variables of the model.For the power-law distributed parameters an average and standard deviation is given as well.
Table 3 and Fig. 3 show that the inferred metallicity from the inference is around the solar value.This is in good agreement with results from our analytic approach in B22, where we calculated Z = 1.0 Z ⊙ using the ([Ne ii] 13 µm + [Ne iii] 16 µm)/Hu α line flux ratio.The probability is little outside the given uncertainties.
With 7.70 +8.58 −5.25 × 10 8 L ⊙ we obtained a somewhat lower stellar luminosity L ⋆ than previous studies, although with a large uncertainty and broad PDF as can be seen in row 1 of Fig. 3. Watson et al. (1996) estimated a luminosity from young stars of 1.5 × 10 9 L ⊙ .Radovich et al. (2001) calculated 1.5 × 10 10 L ⊙ , although with a much poorer spatial resolution (∼ 2 ′ ) and a higher extinction correction of 11 mag, compared to our result of A V = 4.35 mag in B22.
The mean age of the nuclear stellar cluster that we obtained from our solution is 2.54 +1.04 −1.01 Myr.This is slightly younger than previous estimates such as Kornei & McCrady (2009) respectively, perhaps due to a different choice of initial-massfunctions (IMF).We will investigate effects of different IMFs in Sect.5.3.
One major question that is tackled with this study is the characterisation of the central BH, because the physical properties and the impact of the BH on the ISM are not completely understood.The total X-ray luminosity obtained from our models, which we assume to originate from the central BH or AGN is 5.9 × 10 39 erg s −1 , which is in good agreement with Chandra observations (4.7 × 10 39 erg s −1 , Lopez et al. 2023) and Fernández-Ontiveros et al. (2009).The X-ray luminosity shows large uncertainties, but the PDF shows only low probabilities for luminosities below 10 39 erg s −1 and above 4 × 10 40 erg s −1 (see in Fig. 3).To confirm that an X-ray component is important in this model, we made runs with no X-ray source.For both architectures, the pnσ values and marginal likelihood drop by several percents, and the highly ionised emission lines are not well reproduced.The X-ray component is indeed necessary to model the observed emission lines, in particular to recover the [O iv] and [Ne v] emission.For comparison, Sgr A * has a similar luminosity of L X = 4 × 10 39 erg s −1 (Kaneda et al. 1997), suggesting that the central BH of NGC 253 resembles a low-luminosity AGN rather than a typical extragalactic AGN which have luminosities in the range of L X = 10 40 − 10 43 erg s −1 (e.g.Fornasini et al. 2018).
The power-law distributions of the density n, ionisation parameter U, and cut parameter ξ show that a model with low density (i.e.diffuse), low depth, and low ionisation parameter clouds is preferred.However, the range of these three parameters is small (see lower and upper boundaries in Table 3), suggesting that most of the ionised gas emission arises from a mostly diffuse gas component.The average gas density in the power-law model is 135 cm −3 , which is in good agreement (although with large uncertainties) with our analytic results in B22, where we used the [O iii] 52/88 µm, [S iii] 19/33 µm, and [N ii] 122/205 µm line flux ratios and obtained 84 ≲ n ≲ 212 cm −3 .As expected, the range for cut-values ξ is somewhat beyond the ionisation front ξ = 1.0.A fraction of the emission from species with ionisation potentials near 13.6 eV (e.g.[N ii]) possibly originates beyond the ionisation front, which is why ξ ≡ 1.0 or below is not an expected or reasonable solution.
In Table 3 we also show the resulting parameters for the 1C2S configurations.Both, the PLaw and 1C2S architecture are in good agreement within the uncertainties, showing again that the choice of the configurations seems to be only a second order effect when considering only the ionised gas emission lines.Notes. (a) Fixed in the respective grid.

Influence of cosmic rays
Using data from the ALCHEMI spectral survey (Martín et al. 2021), Behrens et al. (2022) showed that the CRIR in the centre is three to four orders of magnitudes larger than the average Milky Way value ζ 0 , while the value used in the SFGX grid is fixed to only 3ζ 0 .Such high values imply a considerable increase of heating to the ISM and will in particular affect emission lines from the neutral atomic and molecular gas (e.g.Sternberg & Dalgarno 1995;Goldsmith 2001;Bisbas et al. 2021).Here we investigated if the higher CRIR also affects emission from the ionised gas.For this purpose we created new Cloudy models using the best parameter set obtained from the model results in Sect.3.3 and changed the CRIR from 6 × 10 −16 s −1 to 10 −13 s −1 as reported in Behrens et al. (2022).The Cloudy models were combined according to their weights or covering factors W i as shown in Table 3.As expected, emission lines from the ionised gas are hardly affected by a change in the CRIR (see Fig.  5 show that the difference in these lines is of the order of one magnitude or even more.

New model grid and sanity checks
To take the higher CRIR and their effect on the neutral atomic (and molecular) gas into account, we created a new sub-grid of SFGX, with an increased CRIR to 10 −13 s −1 .To reduce the computing time for calculating the Cloudy models and for the MULTIGRIS runs, we decreased the parameter space as listed in Table 4.The resulting new grid contains a total of 10880 models.We combined this new grid with the corresponding sub-grid from SFGX, so that we now have a sampling for the CRIR.We will call this new grid SFNX (for star-forming nucleus with an X-ray source) in the following3 .In Sect.3.4 we showed that an increased CRIR only has a minor impact on the line fluxes of ionised gas lines.To confirm that the effect is negligible, we proved that MULTIGRIS finds similar parameter sets in SFNX and in SFGX.We carried out another run with a two component configuration, using only emission lines originating in the ionised gas as done in Sect.3.2.Table 3 shows the resulting mean values and confidence intervals for both runs, which are in good agreement.The only exception is L X , where the results from the new grid are lower but still within uncertainties.This is most likely due to a degeneracy effect between L X and ζ (see Sect. 5.1 for a more detailed discussion) and the higher values chosen for L ⋆ = 10 9 L ⊙ in order to be consistent with the SFGX grid.Since the new grid is able to reproduce the results from SFGX, but also takes the much higher CRIR into account, we were now able to investigate solutions for the neutral atomic gas.

Results for the ionised and neutral atomic gas
In the next step we added emission lines from the neutral atomic gas with the systematic offsets as determined in Sect.3.1 remaining the same.As mentioned earlier, we had not yet included the [Fe iii] 23 µm emission line due to the dependence on the dust- to-gas ratio.Adding emission lines from the neutral atomic gas such as [Fe ii] 18 µm and [Fe ii] 26 µm as constraints, we now enable a small systematic offset for the three Fe lines.We assume that the line fluxes scale linearly with the iron abundance (within only a small offset), and account for small variations in the dustto-gas ratio.We obtain an offset of ∼ 0.3 for both configurations, suggesting that the model would under-predict the Fe lines without the scaling.Since H 2 rotational transitions arise to a significant amount from the warm neutral atomic gas (e.g.Roussel et al. 2007;Togi & Smith 2016) we also include these lines in these new runs.Although we include emission lines from molecular hydrogen, we do not claim that our model is a proper solution for the molecular gas in the centre of NGC 253.
Most of the parameters obtained for the ionised gas model (Sect.3.3 and Table 3) are consistent with the model of the ionised and neutral atomic gas as shown in Table 5.This shows that the approach to start with only ionised gas lines and determining systematic offsets between instruments (Sect.3.1) is valid.Generally, the preferred model in the PLaw architecture has a high abundance of low density, low ionisation parameter and low depth (i.e.diffuse gas) clouds.All three parameters (n, U, and ξ) have a negative slope of the power-law, meaning that higher density clouds with larger depths and higher ionisation parameter are less abundant.In the 1C2S architecture, this is reflected by a smaller covering factor W 2 for the second component.However, some parameters show differences compared to 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Cut the results from the ionised gas.Obviously, the cut parameter ξ is larger than in the model for the ionised gas.To reproduce emission lines from the neutral atomic gas, the model has to reach values beyond the ionisation front (ξ = 1) and even the H/H 2 photo-dissociation front (ξ = 2).Since we do not include CO or other emission lines from the molecular core, larger depths (ξ ⩾ 3) are not needed nor found.The obtained average ξ ≈ 2.3 (cf.Table 1) translates into an A V of 4.5 which is in good agreement with our results in B22 (see also Fig. A.2). Furthermore, the density -in particular the upper limit n upper for the PLaw architecture and the density of the second component n 2 for the 1C2S architecture-is significantly higher than in the ionised gas model.This is expected, since the model needs to account for higher critical densities of many of the emission lines from the warm neutral gas such as The additional parameter introduced for the CRIR, ζ, is ∼ 2 × 10 −15 s −1 .This is significantly lower than predicted by Behrens et al. (2022) (∼ 10 −13 s −1 ), which is expected for two reasons.First, the region observed is much larger in our case (40 ′′ compared to 1.6 ′′ in Behrens et al. 2022), so the locally in-creased CRIR from Behrens et al. (2022) is smeared out within our observing beam.Second, although many lines are sensitive to the CRIR as shown in Sect.3.4, higher fluxes for these lines can be modelled for instance by a higher density, a stronger radiation field and/or a higher X-ray luminosity (see Sect. 5.1).It is possible that there is some degeneracy between these parameters, which may also be reflected by the large uncertainties.

Cosmic ray ionisation rate and X-ray luminosity
Table 3 shows that the X-ray luminosity drops significantly from ∼ 5 × 10 39 erg s −1 to ∼ 10 37 erg s −1 , once the CRIR is increased.This somewhat contradicts our finding that X-rays are needed to reproduce the highest ionised emission lines in our sample.Hence, we performed another MULTIGRIS run, forcing a higher X-ray luminosity of L X = 5 × 10 39 erg s −1 .The resulting parameter sets for both architectures (PLaw and 1C2S) are similar, with the exception of the CRIR, which drops by ∼ 25% in the higher X-ray runs, consistent with findings from Lebouteiller et al. (2017).Cosmic rays affect mostly the neutral atomic and molecular gas, while the ionised gas remains unchanged.X-rays, however, affect the ionised gas -and high ionisation states in particular-but also heat the neutral atomic gas, which creates a degeneracy between the luminosity of the X-ray source and the CRIR.Such difficulties between the discrimination of X-rays and cosmic rays has already been shown, for instance by Meijerink et al. (2006).Observations of the CO + molecule could help to break this degeneracy, as done in M82 (Spaans & Mei-jerink 2007).In conclusion, this implies that both, the CRIR and X-ray luminosity obtained in this study (Tables 3 and 5) should rather be considered as upper limits.

Predictions for the molecular gas
From the PDFs of the model parameters, MULTIGRIS is able to predict the luminosity of other emission lines.We predicted line fluxes for several CO emission lines (J = 1 → 0 to J = 8 → 7) that arise from the molecular gas.We compared the predicted line fluxes with those observed by Herschel/PACS and Herschel/SPIRE (see Pérez-Beaupuits et al. 2018) but did not use them as constraints.Figure 6 compares the observed and predicted CO line fluxes and the corresponding uncertainties in column 4 (grey background).With the exception of CO(1−0) in the 1C2S architecture, all CO emission lines are under-predicted by up to two orders of magnitude.The model clearly fails to account for the molecular gas.We investigated potential causes for the under-prediction of the CO emission.However, for a more detailed model of the molecular gas we refer to a forthcoming paper.Due to the in- creased CRIR, the ISM is much warmer and the H/H 2 and C/CO photo-dissociation front are shifted to larger depths with increasing ζ.Since the Cloudy models, however, stop at a fixed A V = 10 mag, models with a higher CRIR contain less H 2 and therefore also less CO.This, in turn, results in significantly lower CO line fluxes simply due to the low abundance of CO in the whole Cloudy model.To overcome the low CO abundance, we investigated which A V is needed to predict the observed CO line fluxes, by running new Cloudy models with parameters from Table 5 but with an arbitrary high A V of 30 mag.In fact, such a model would be able to reproduce the CO line fluxes as observed by Herschel/SPIRE and JCMT, but the increase of depth obviously affects other lines in the model as well.In particular emission lines like the [C i] and [O i] lines originate partly in the molecular gas and become brighter with larger depth.See also Fig. 4 showing that the cumulative fluxes of these lines are growing even at larger ξ values.Furthermore, we calculated the total molecular hydrogen mass M(H 2 ) from these high A V models (see Sect. 5.4), which is two orders of magnitude larger than typical hydrogen masses in the nuclei of galaxies.To conclude, simply increasing the depth of the models cannot overcome the under-predicted CO emission for which other emission lines would be affected as well as an absurdly high molecular gas mass would be obtained.Furthermore, it would be in contradiction with our results for an A V of 4.35 from the SED fit in B22.
Besides a higher A V , shocks could play a role in the excitation of CO emission (e.g.Lesaffre et al. 2013;Pon et al. 2016;Kamenetzky et al. 2018).Furthermore, Hao et al. (2009) showed that shocks can possibly excite [O iv] 26 µm emission, which is significantly under-predicted in our models as well.Another mechanism not taken into account is the time-dependence of the chemistry.Cloudy assumes a static chemistry, however, simula-tions have shown that a time-dependent chemistry (in particular in nuclei of active galaxies) can result in several orders of magnitude brighter CO line fluxes (Meijerink et al. 2013).

Effects of the initial mass function
To determine the emitted spectrum of an interstellar cloud, Cloudy requires an input spectrum, in our case the spectrum of a starburst cluster.The distribution of the overall mass among the stars within the cluster (the so-called initial-mass-function, IMF) can have a significant effect on the output spectrum of the cluster -and hence also on the physical conditions and spectrum of the illuminated cloud.However, the exact shape of the IMF and whether it universal, is still uncertain and debated (e.g.Hopkins 2018, and references therein).We analysed the impact of a different IMF, by creating a new input spectrum -to eliminate potential uncertainties by using different codes (and therefore different assumptions such as evolutionary tracks, etc.), we used BPASS to predict the stellar emission.While the other parameters remain the same (Eldridge et al. 2017;Ramambason et al. 2022), we change the IMF from a Kroupa-like IMF (Kroupa et al. 1993) to a Salpeter IMF (Salpeter 1955).Thereafter, we ran two Cloudy simulations with the parameters from Table 5, one using the spectrum from a Kroupa IMF, and one using the spectrum from a Salpeter IMF.
Figure 7 shows the resulting line fluxes from these Cloudy simulations.The line fluxes from the two models are generally in good agreement, in particular emission lines originating in the neutral atomic (e.g.[O i] and [C i]) and molecular gas (CO and H 2 ), where they agree within 1%.Emission lines, that to some fraction come from the ionised gas, are slightly more affected.The largest deviations occur for emission lines that purely arise from the ionised gas, with the Kroupa IMF yielding fainter line fluxes than the Salpeter IMF in almost all cases.For instance, the two [S iii] lines are brighter by ∼ 10%, the two [Ne iii] lines even by ∼ 20%.However, the ratio of two emission lines from the triplet of a species (e.g.[N ii], [S iii], [O iii]) are constant for one IMF.This indicates that the reason for the differences is primarily due to different abundances of the respective ions.Indeed, the abundance of higher ionic states (which is one output of Cloudy) are 5 − 10% higher for the Kroupa IMF at depths lower than the ionisation front.We conclude that the particular shape of the chosen IMF can have a significant impact on the resulting line fluxes, in particular those from emission lines originating in the ionised gas.However, by selecting a state-of-the-art IMF and model BPASS, we assume that these effects are minor and that our results are valid.

Secondary parameters obtained from the model
In addition to predicting line fluxes of other emission lines, MULTIGRIS also allows us to predict secondary parameters in post-processing runs.These secondary parameters must be stored in a post-processing grid, which is associated with the main grid, in our case SFGX and SFNX, respectively.The values listed throughout this section are results from the PLaw architecture.However, similar to the results of the primary parameters in Table 5, the post-processing parameters are comparable and agree within uncertainties.Assuming a spherical geometry, masses of dust and the different hydrogen phases were calculated from the Cloudy output.Using our final model (see Table 5) we are then able to estimate the mass of ionised and neutral atomic hydrogen in the centre of NGC 253.In principle, we can also obtain a mass for the molecular hydrogen.However, due to the limited validity of our model regarding the molecular gas, results for this phase have to be taken with caution, as noted by the large uncertainties associated with these parameters.The resulting masses are shown in Table 6.Both, the mass for the ionised hydrogen and for neutral atomic hydrogen are in good agreement with estimates from KAO observations (Carral et al. 1994), who obtained 1 × 10 6 M ⊙ and 5 × 10 6 M ⊙ from the Tielens & Hollenbach (1985) PDR models for the ionised and neutral atomic hydrogen, respectively.We can also estimate the fraction of [C ii] 158 µm associated with the ionised, neutral atomic, and molecular gas.Recently, [C ii] has been increasingly considered as a probe for the CO dark gas or even the total molecular gas mass (e.g.Madden et al. 2020, and references within).It is necessary, however, to correct for emission coming from H ii regions that do not contain any H 2 .In the case of NGC 253 with solar metallicity, the fraction of [C ii] 158 µm from the ionised gas phase is ∼ 26% (see Table 6).This is somewhat lower than the findings in Cormier et al. (2019) (using their Eq.(2) yields 45% of [C ii] 158 µm from the ionised gas), which could be due to different CRIR assumed in our study or the fact that solar metallicities are not covered in Cormier et al. (2019).

Star-formation rates
In Sect.3.3 we showed that the potential AGN in the centre of NGC 253 has little effect on its environment as probed by the observed lines.What really drives the heating in the nucleus seems to be the star-formation activity, usually quantified by the star-formation rate (SFR).Within the last few decades, a number of tracers for the SFR have been proposed, either from theoretical considerations or empirical calibrations (see Kennicutt & Evans 2012, for a review).Each of these tracers has its (dis-)advantages, in particular regarding where the respective probe arises, and if it is affected by extinction.In this section we will compare three different SFR-tracers, namely the Hα, [C ii] 158 µm, and L TIR luminosities.Throughout we assume a distance of 3.5 Mpc (Rekola et al. 2005) to convert fluxes to luminosities.
One of the most important tracers of the SFR is the luminosity of the Hα line, since it is easily observable with optical telescopes and is directly emitted by young stars.From our latest model presented in Sect.4, MULTIGRIS is able to infer the intrinsic (extinction-free) luminosity of this emission line (Table 6).Using the empirical calibration from Calzetti et al. (2007), which is also tested on smaller, ∼ 100 pc scales (Pessa et al. 2021;Belfiore et al. 2023), we obtain SFR = 2.3 +3.1 −0.4 M ⊙ yr −1 .One disadvantage of the Hα emission line is that that it that it typically suffers from extinction in extra-galactic sources, which is why we corrected the line flux from Table 6 assuming a mixed extinction model and the optical depth of 4.5 mag as determined in Sect.4.2.
Another frequently employed probe for the SFR is the total infrared luminosity L TIR .It assumes that in thermal equilibrium, most of the stellar radiation is reprocessed by dust and radiated within the infrared spectral range (i.e. between 3 µm and 1000 µm).We determined the L TIR already in B22 in two different ways, confirmed by our MULTIGRIS approach.The calibration shown in Kennicutt (1998) yields SFR = 1.6 ± 0.4 M ⊙ yr −1 , which is in good agreement with the SFR determined from Hα. Stacey et al. (1991) proposed the [C ii] 158 µm emission line as a probe for the SFR.Since we directly measured this quantity and did not calculate or predict it, we believe that this is the most reliable measurement of the SFR in this study -not taking any calibration uncertainties into account.Using the more recent calibration from Herrera-Camus et al. (2018), we obtain SFR = 1.7 ± 0.5 M ⊙ yr −1 .
All the SFRs we obtained from the different probes are in good agreement with each other.They are also comparable with the results from ISO observations in Radovich et al. (2001), who obtained SFR = 2.1 M ⊙ yr −1 for the nucleus, however, with a larger beam of the ISO telescope compared to our observations.The larger observing beam could result in contamination of infrared emission from the disk and therefore in a slightly higher SFR.Yet, the result from Radovich et al. (2001) is within the uncertainties of our solution.

Comments and caveats of the model grids
Although the model grids SFGX and SFNX that we use cover a wide range of ISM conditions, we note that some mechanisms are not taken into account which may affect the model results, and in particular might lead to the under-prediction of the CO emission lines.These aspects shall be shortly discussed and will be further investigated in future studies.
According to Sánchez (2020) or Sánchez et al. (2021), one potential major contribution to the radiation field can be postasymptotic-giant-branch (AGB) stars.These objects typically have a hard but weak ionising effect on the ISM, which could cascade down into the CO emitting regions.AGB stars are included in the BPASS models, however, the post-AGB phase is not well understood and needs improvements in the model (Eldridge et al. 2017).While the works of Sánchez (2020); Sánchez et al. (2021) show that ionisation occurs on local scales, the limited spatial resolution of our observations in this work leads to the fact that the ionising sources are not resolved.
Another mechanism not considered in the model grids are shocks, for instance from supernovae or galactic outflows.Such outflows have been observed in CO, Hα, and X-ray emission (Bolatto et al. 2013).They are well known to contribute to ISM heating and therefore boost emission lines.
Lastly, the relative abundances of elements varies not only within galaxies but even when comparing H ii regions within one galaxy (García-Rojas & Esteban 2007).The model grids take into account that the abundance, for example of N and C, varies with the O/H ratio.However, it might also be the case that the N and C relative abundance changes as well.A deviation in the relative abundances of elements would directly lead to a change of the relative line fluxes of different species.

Summary
In this study we used a combined set of 30 emission lines from SOFIA, Herschel, and Spitzer observations of the nuclear region of the starburst galaxy NGC 253 to investigate the physical conditions of the ISM on a 100 pc scale.Using MULTIGRIS, a Bayesian code to probabilistically investigate a set of Cloudy models, we constrained parameters of the ionised and neutral atomic gas.After eliminating systematic offsets between the different telescopes, instruments, and modules, we first determined parameters of the ionised gas.The created model was able to reproduce all observed emission lines with the exception of [O iv] 26 µm.We found that the metallicity and density that we calculated from analytic prescriptions in B22 are in good agreement with our probabilistic results.Furthermore, we inferred that the hypothetical AGN in the nucleus of NGC 253 has a minor impact on the heating of the ISM, with luminosities ≲ 6 × 10 39 erg s −1 or ≲ 1.5 × 10 6 L ⊙ .After modelling the ionised gas, we added emission lines originating in the neutral atomic gas, with an increased CRIR as shown by Behrens et al. (2022).We showed that the higher CRIR has little to no effect on the results for the ionised gas.Again, the extended model is able to reproduce most (24 of 30) of the emission lines, and the obtained optical depth is in good agreement with our results from B22.However, the model fails to reproduce most of the CO emission -we refer to a future study to further investigate the molecular gas properties.From our model we were able to calculate gas and dust masses, and determined the fraction of [C ii] 158 µm emission from the different phases (Table 6).Since the nuclear starburst seems to have the major heating impact on the ISM, we calculated the star-formation rate from different tracers which are all in good agreement (0.6 − 1.7 M ⊙ yr −1 ).

Fig. 1 .
Fig. 1.Optical image from Hubble/WFC3 observations of the central ∼ 2 ′ of NGC 253.Overlays show the apertures used to extract the line fluxes from the different observatories.The apertures from Spitzer/IRS are grey (short-high) and black (long-high), green shows the footprint of the Herschel/PACS integral-field-unit.Blue circles are the apertures for the SOFIA/FIFI-LS observations of [O iii] 52 µm (solid), [O iii] 88 µm (dashed), [O i] 146 µm (dashdot), and [C ii] 158 µm (dotted) as described in B22.See also Table2.

Fig. 2 .
Fig. 2. Comparison of modelled/predicted and observed line fluxes.Emission lines from the ionised gas are used as constraints.The abscissa is normalised to the observed line flux (≡ 0) with errorbars (blue).Red and green show the resulting line fluxes from the PLaw and 1C2S architecture respectively.For display reasons, the abscissa in column one are compressed so that the PLaw results for the [Ne v] lines are not entirely visible.Blue arrows denote upper limits on line fluxes.
B.1).The difference in cumulative line flux slightly beyond the ionisation front (i.e.ξ = 1.25) is lower than 10% for most of the lines.Only Hu α, [Ne ii] 13 µm, and [S iii] 33 µm have a larger but still small difference of ≲ 20%.Emission lines from the neutral atomic and molecular gas however, are heavily affected by a change of the CRIR.Figures 4 and

Fig. 3 .
Fig.3.Probability density distributions for the stellar luminosity L ⋆ , Xray luminosity L X , and metallicity Z for the PLaw architecture (seeSect.3.3).Red shows the median and confidence intervals of the posterior distribution.

Fig. 5 .
Fig.5.Same as Fig.4but for emission lines from the molecular gas and the total infrared luminosity L TIR .

Fig. 6 .
Fig.6.Same as Fig.2.Emission lines from the ionised and neutral atomic gas are used as constraints (white background) with CO emission lines predicted by MULTIGRIS (grey background).

Fig. 7 .
Fig. 7. Comparison of predicted line fluxes from Cloudy runs, using a Kroupa (blue) and Salpeter (orange) IMF input spectrum, respectively.For representation reasons the [Ne v] emission lines are not shown but are both of the order of 10 −28 W m −2 , i.e. within the upper limits of our observations. 1.

Table 2 .
Ferland et al. 2017loudy (a 1D radiative transfer code, seeFerland et al. 2017) models assuming a cloud of spherical geometry.The modelled cloud is illuminated by a representative stellar cluster of a given age t and luminosity L ⋆ .The spectrum of this star cluster is created with BPASS(Stanway & Eldridge

Table 1 .
Definition of the different values for the cut parameter ξ in the SFGX grid.There are tabulated decimal values between the ξ values.

Table 2 .
Mean systematic line flux offsets due to calibration or aperture effects between observation blocks determined from MULTIGRIS runs in linear scale.Only ionised gas lines were used as constraints.The aperture sizes (in ′′ and pc) are added for clarity.

Table 3 .
Resulting mean values and confidence intervals from the inference for the PLaw model from SFGX (see Sect. 3.3) and SFNX (see Sect. 4).Both models results from runs where only emission lines from the ionised gas are considered as constraints.

Table 4 .
Input parameters of Cloudy models in the new sub-grid of SFGX and parameter space of the original SFGX grid.

Table 5 .
Resulting parameters from an inference using emission lines from the ionised and neutral atomic gas as constraints, for the two different architectures discussed in this study.The model grid used for the inference was SFNX.

Table 6 .
Resulting post-processing parameters for the nucleus of NGC 253.