Infrared view of the multiphase ISM in NGC 253 I. Observations and fundamental parameters of the ionised gas

Context. Massive star-formation leads to enrichment with heavy elements of the interstellar medium. On the other hand, the abundance of heavy elements is a key parameter to study the star-formation history of galaxies. Furthermore, the total molecular hydrogen mass, usually determined by converting CO or [C ii] 158 $\mu$m luminosities, depends on the metallicity as well. The excitation of metallicity-sensitive emission lines, however, depends on the gas density of H ii regions, where they arise. Aims. We used spectroscopic observations from SOFIA, Herschel, and Spitzer of the nuclear region of the starburst galaxy NGC 253, as well as photometric observations from GALEX, 2MASS, Spitzer, and Herschel in order to derive physical properties such as the optical depth to correct for extinction, as well as the gas density and metallicity of the central region. Methods. Ratios of the integrated line fluxes of several species were utilised to derive the gas density and metallicity. The [O iii] along with the [S iii] and [N ii] line flux ratios for example, are sensitive to the gas density but nearly independent of the local temperature. As these line ratios trace different gas densities and ionisation states, we examined if these lines may originate from different regions within the observing beam. The ([Ne ii] 13 $\mu$m + [Ne iii] 16 $\mu$m)/Hu $\alpha$ line flux ratio on the other hand, is independent of the depletion onto dust grains but sensitive to the Ne/H abundance ratio and will be used as a tracer for metallicity of the gas. Results. We derived values for gas phase abundances of the most important species, as well as estimates for the optical depth and the gas density of the ionised gas in the nuclear region of NGC 253. We obtained densities of at least two different ionised components $(<84$ cm$^{-3}$ and $\sim 170 - 212$ cm$^{-3})$ and a metallicity of solar value.


Introduction
NGC 253 is the nearest spiral galaxy showing a nuclear starburst with a star-formation rate (SFR) of 3 M yr −1 (Radovich et al. 2001). Its proximity makes this galaxy a perfect model to study the interstellar medium (ISM) in the extreme environments of massive star-formation. Near-infrared observations revealed the presence of a bar which funnels gas and dust from outer regions into the nuclear starburst (Iodice et al. 2014). Besides fuelling of the starburst by the stellar bar, an outflow exiting from the nuclear region has been observed in Hα (Bolatto et al. 2013a), and CO (Leroy et al. 2015;Walter et al. 2017). From these observations, a mass outflow rate of about 3 − 9 M yr −1 has been estimated, implying that the nuclear starburst is severely suppressed by the outflow. The existence of a central supermassive black hole has been discussed within literature (e.g. Vogler & Pietsch 1999;Günthardt et al. 2015). The putative central supermassive black hole and subsequent active galactic nucleus (AGN, Fernández-Ontiveros et al. 2009) should have a noticeable impact on the ISM and the star-formation activity in the nuclear region, by accreting the surrounding stars and ISM at high velocities. The accretion creates highly ionised species such as O 3+ or Ne 4+ , which are also observable via their fine-structure lines in the mid-infrared (MIR, Fernández-Ontiveros et al. 2016, and references therein). The outflows originating from the starburst and the hypothetical AGN reduces the amount of the cold molecular gas reservoir, which is needed to sustain the ongoing starburst. The nature of this hypothetical black hole, however, its signatures in the surrounding ISM and consequences for star-formation and the H 2 gas mass reservoir are still undetermined (Fernández-Ontiveros et al. 2009). Due to the presence of this many different components − bar (Engelbracht et al. 1998), starburst (Rosenberg et al. 2013), outflows (Bolatto et al. 2013a;Krips et al. 2016), and possibly a black hole or AGN (Fernández-Ontiveros et al. 2009) − the major energy input to the ISM is uncertain and will be tackled in this series.
Additionally, there is still little knowledge about the chemical composition of the gas in the nuclear region of NGC 253 (Martín et al. 2021). This is mostly due to difficulties that occur when optical lines like [O ii] 3737Å or [O iii] 5007Å are observed in such a highly extincted starburst galaxy. The metallicity, however, is a key parameter when trying to understand the ongoing mechanisms in the ISM. Since metals are by-products of star-formation, the metallicity of a galaxy tracts its evolutionary history. Additionally, to estimate the molecular hydrogen gas mass reservoir of galaxies from CO emission, the conversion factor is known to be a function of metallicity (e.g. Bolatto Spinoglio et al. 2022). In a followup paper we present a more complex multiphase model using MULTIGRIS Ramambason et al. 2022) a Bayesian method determining probability density distributions of parameters of a multisector Cloudy (Ferland et al. 2017) grid from a set of emission lines. With this approach we will constrain parameters of the ISM such as intensity of the radiation field of different components of the ISM, temperature and luminosity of the putative black hole or AGN, and age of the starburst.
This work is structured as follows: Section 2 describes our observations and the data reduction, as well as archival data from different observatories, and the determination of line fluxes and errors. Section 3 follows with an analysis using some of the observed lines, determining the density and metallicity of the ionised gas, as well as the optical depth of the nuclear region. A summary and outlook of future work concludes this study in Sect. 4.

SOFIA/FIFI-LS
We used the Field-Imaging Far-Infrared Line Spectrometer (FIFI- LS Fischer et al. 2018)   Notes. Neutral and ionised species, observing flight dates, spectral and spatial resolution (FWHM of PSF) at the given wavelength of the observed lines. More details about spectral resolution can be seen in Fischer et al. (2016) obtain emission line data. FIFI-LS is an imaging spectrometer for the FIR providing a wavelength dependent spectral resolution of R = 600 − 2000. This instrument provides two independent observing channels allowing simultaneous observations in the λ = 51 − 120 µm (blue channel) and in the λ = 115 − 203 µm (red channel) wavelength regimes. Each channel projects an array of 5 × 5 spatial pixels (spaxels) on the sky, covering a total area of 30 × 30 and 60 × 60 in the blue and red channel, respectively (Klein et al. 2014;Colditz et al. 2018).
Observations were carried out on SOFIA flights starting from Palmdale, California during cycle 3 in October 2015, cycle 6 in November 2018, and cycle 7 in November 2019 as guaranteed time campaigns. A summary of the observed lines, flight dates, and spatial and spectral resolution is listed in Table 2. Figure 1 shows an outline of our observations (white contours).  Fig. 2. Spectra of the pipeline spaxel (1.5 × 1.5 ) located toward the nucleus. The orange spectrum shows the observed spectrum (including errors from the ramp fit) without correction for atmospheric transmission. The telluric absorption feature at 51.7 µm can clearly be seen in the continuum. The blue spectrum shows the pipeline-corrected spectrum, including the errors determined as described in Sect. 2.1. The grey area shows the spectral region used to calculate flux density uncertainties. Green shows the best fit from the Monte-Carlo method.
The FIFI-LS data were reduced with the FIFI-LS standard data reduction pipeline (Vacca et al. 2020). The pipeline performs a flux calibration, and resamples the data onto a spatial grid of 1.5 × 1.5 and 3 × 3 (pipeline spaxels) in the blue and red channels, respectively. The transmission of the atmosphere and hence the observed line flux at a given wavelength depend mainly on the total upward atmospheric precipitable water vapour (PWV). Inflight PWV values are determined by targeting specific water vapour features which are used in the data reduction pipeline for atmospheric calibration purposes (see Fischer et al. (2021) and Iserlohe et al. (2021)  Flux density errors in the data cubes are the propagated errors from the ramp fit of the readout of the detector, but do not include errors from flux calibration or the PWV correction (FIFI-LS team, private communication). We therefore estimated the statistical flux density uncertainty from the root-mean-square (RMS) error in the continuum region of each pipeline spaxel. Due to the small spectral coverage the continuum is nearly constant and the RMS error is not amplified by a non-constant continuum. For the [O iii] 52 µm, [O iii] 88 µm, [O i] 146 µm, and [C ii] 158 µm observations we excluded spectral regions with the expected line position. Due to a strong absorption feature at 51.7 µm, the flux density < 51.8 µm was also not included for the root-mean-square calculation (see Fig. 2 for an example spectrum and fit).
To determine the line fluxes and their statistical uncertainties in each pipeline spaxel, we used a Monte-Carlo approach by varying the flux density within the normal-distributed flux density uncertainties in each spectral bin randomly over 500 steps. The spectrum obtained from this variation was fitted with a linear continuum and a Gaussian with variable mean, width and amplitude (see Fig. 2). From the integral over the Gaussian component we obtained the line flux in each step − giving a sample of 500 line fluxes in each pipelines spaxel, from which we calculated the mean and standard deviation, which are the line flux and line flux error, respectively. These procedures were applied to the spectra from all pipeline spaxels, yielding the line flux maps shown in Fig. 3.

Spitzer/IRS
We complement our data set by single slit spectra obtained with the Infrared Spectrograph (IRS, Houck et al. 2004) onboard the Spitzer Space Telescope. Data are available at the Combined Atlas of Spitzer IRS Spectra (Lebouteiller et al. 2015, CASSIS 1 ) under AOR key 9072640.
While the nuclear region of NGC 253 was observed with low-and high-resolution modules, we only use the short-high (SH) and long-high (LH) module observations, from 9.6 − 19.6 µm and 18.7 − 27.2 µm, since the low-resolution observations are saturated. Both high-resolution modules have a velocity resolution of ∆v = 500 km s −1 , comparable to the spectral resolution of our FIFI-LS observations. Slit sizes of the shortand long-wavelength module are different, 4.7 × 11.3 and 11.1 × 22.3 , to better match the PSF size of each module.
CASSIS offers two reduction processes, the optimal extraction adapted to point-like sources, and full extraction adapted to extended sources. Since the nucleus of NGC 253 is partly extended (see Sect. 2.3.2 and Fig. 6), we chose full extraction procedure as this ensures we retrieve most of the source's flux. The full extraction accounts for losses of flux density, and integrates the flux density over the full slit size. Although the nucleus of NGC 253 is not uniformly extended over the whole slit, which would be the ideal application for the full extraction, this process yields better results than optimal extraction. For each emission line in the IRS spectrum we fitted a Gaussian line profile and a polynomial of order 0 to 4 for the underlying continuum using lmfit. From the line fits, we obtained the line flux and line flux errors. In cases where several emission lines were close together, we fitted the continuum and all lines simultaneously. The high order for the continuum is necessary to account for the several broad emission and absorption features from, for instance, polycyclic aromatic hydrocarbons and silicates in the MIR. Line centroids are allowed to vary within 400 km s −1 with respect to the expected velocity from the source (see Table 1). The minimum line width is 300 km s −1 , which is about the spectral resolution of IRS. Some lines are only weakly detected. We give an upper limit for the line flux, if the obtained line flux from the fit is smaller than 2× the RMS of the flux density of the continuum (in W/m 2 /µm) times the spectral resolution (in µm). In this case, the upper limit is given by 2× the RMS error of the continuum. W m −2 on a 3 × 3 (left column) and 1.5 × 1.5 (right column) grid, respectively. Pixels with a signal-to-noise ratio smaller than 3 are masked. The black circles depict the aperture used to extract the spectrum from the nuclear region (Sect. 2.3.1). We also show contours from the fitted underlying continuum (grey).

Herschel
We add to the FIFI-LS and Spitzer data observations with the Photodetector Array Camera and Spectrometer (PACS, Poglitsch et al. 2010) onboard the Herschel Space Observatory 2 (Pilbratt et al. 2010). We used PACSman v3.63  for data reduction and analysis of the PACS spectral cubes. PACSman uses the rebinned data cubes from the PACS standard pipeline and fits a Gaussian emission line and a second order polynomial for the underlying continuum in each of the 9.4 × 9.4 spaxels. PACSman provides a Monte-Carlo based method to determine the line flux errors in each spaxel. Similar to the procedure in Sect. 2.1, the flux density is varied within the Gaussian-distributed uncertainties in each of the 3000 steps. Each of the obtained spectra are fitted yielding a line flux in each spaxel. From the sample of 3000 line fluxes we calculated the mean and standard deviation, ultimately provides the line flux and line flux error in one spaxel, respectively.
We obtained line flux maps and corresponding uncertain- Fig. 4 for the observed spectra from the central 3 × 3 spaxels of the field-integral unit (cf. Sect. 2.3.3).
We finally complement our data set with Herschel/SPIRE-FTS observations. The SPIRE line fluxes from Pérez-Beaupuits et al. (2018) are already corrected for the semi-extended nuclear source (see Sect. 2.3), hence, we do not reduce these observations and list only the results in Table 3

Data homogenisation
We focus on emission from the nuclear source in this study. Hence, we have to determine the line fluxes and line flux errors arising only from the nucleus for each instrument.

SOFIA/FIFI-LS
The FIFI-LS line flux maps (Fig. 3) are oversampled, without knowledge of the correlation between the pipeline spaxels. In addition, the line flux maps are much more extended than only the nucleus and show emission from the bar and disk. To extract line fluxes and corresponding errors from only the nucleus, we proceeded as follows: we calculated an effective PSF size by fitting a Gaussian to the spatial profile of the nucleus for each emission line − the ideal PSF is broadened by the non-Nyquist sampling of the observations and the resample procedure in the pipeline. We resampled both, the line flux maps (Fig. 3) and line flux error maps to a spatial grid of the effective PSF size. For the line flux error maps we assumed that all pixels are fully correlated − meaning that the errors are cumulative. To get most of the line flux from the nucleus, not contaminated from background (i.e. disk-) emission, we extracted the line flux and errors from the resampled maps using a circular aperture with a radius of 2× the effective PSF FWHM, assuming that the pixels are uncorrelated. This means that errors add up quadratically. This procedure yields an upper limit for the line flux errors, since the assumption of the first step that all pixels are correlated is not entirely fulfilled. However, it is a good approximation as the errors are of the same order as uncertainties from Herschel/PACS (see Table 3). The resulting spectra from this procedure are shown in   Carral et al. (1994). Applying the same Monte-Carlo method as in Sect. 2.1 we get a line flux of 107.7 ± 30.8 W m −2 . This is in good agreement with Carral et al. (1994), who yielded a line flux of 90 ± 20 W m −2 .

Spitzer/IRS
The spatial profile of the nuclear source is not point-like, but slightly extended, as can be seen from the cross-dispersion direction of the IRS LH slit (see Fig. 6). From the comparison of the observed spread and the ideal PSF, we get a spatial extent of the nucleus of 6.68 , corresponding to 113.6 pc under the assumed distance of 3.5 Mpc (see Table 1). This is in good agreement with the 6 from [Ne ii] 13 µm observations by Engelbracht et al. (1998). We assume that the source size is the same for all MIR and FIR observations. As the spatial extent of the source is of the order of the size of the SH aperture size, the flux density in the SH observations is naturally underestimated due to a loss of flux outside the aperture. This leads to an offset between the continuum flux density of both, SH and LH modules. To account for any lost flux density in the smaller SH aperture, we stitched the overlapping continua of both modules. We fit the spectrum of both modules in the overlapping continuum region (∼ 18.7 − 19.6 µm) linearly. From the fit we get an offset of the two continua of 1.55, hence, we scale the SH spectrum by a factor of 1.55. We assume that the continuum and line emission arise from the same region, which is valid because the MIR continuum emission is mostly from hot dust in H ii regions. We used broad band photometric images from the Midcourse Space Experiment (MSX 3 ) at λ = 12.13 µm, 14.65 µm, 21.3 µm to assert that no wavelengthdependent scaling is necessary. The more recent WISE observations are saturated and thus cannot be used to verify the scaling. The continuum of the IRS spectra agrees within ∼ 10% with all three MSX bands, thus no wavelength-dependent scaling is needed. The resulting spectrum and MSX measurements are shown in Fig. 7. Footprints of SH and LH observations can be seen in Fig. 1.
Line fluxes from both modules and other parameters of the observed lines are listed in Table 3. Our values are in good agreement with values from Bernard- Salas et al. (2009) who used the same data set.

Herschel/PACS
To get the integrated line flux from the nuclear source from the PACS line flux maps created in Sect. 2.2.2, we summed up the line flux from the central 3 × 3 spaxels (28.2 × 28.2 ) to limit background emission and noise from the weakly illuminated outer pixels, and to only extract the nuclear source flux. A comparison of our line fluxes with Pérez-Beaupuits et al. (2018) shows consistency within the error bars for all lines except [C ii] 158 µm, which is ∼ 20% weaker in our case. However, when comparing the [C ii] 158 µm line flux from all 5 × 5 spaxels, the results are again in good agreement. The large difference between these two values is due to a significant contribution of [C ii] 158 µm from the underlying disk. Line fluxes and errors of the brightest fine-structure lines observed with PACS are listed in Table 3.

Analysis
We used ratios of our MIR and FIR integrated line fluxes to determine parameters of the ISM. These fine-structure lines are in most cases not affected by extinction or self-absorption. Furthermore, when studying line ratios from the same species (e.g. [O iii] 52/88 µm) the ratio is independent of the ionic abundance of O 2+ . This is particularly important, since the elemental abundances are uncertain for NGC 253.

Optical depth
Luminous infrared objects tend to be highly extincted sources, making it sometimes necessary to correct for extinction even at MIR wavelengths, particularly for A V ≥ 10 mag (Armus et al. 2007). Hence, we have to correct the MIR line fluxes for extinction or make sure that such correction is not needed for our modelling approach. The optical depth determinations for the nuclear region in NGC 253 vary within the literature from A V = 5.5 mag (Pérez-Beaupuits et al. 2018, from dust SED fit), A V = 14.0 mag (Rieke et al. 1980, from Brα/Brγ)   emission (e.g. SFR and initial mass function), and dust attenuation and emission (e.g. optical depth and dust mass) and fits the SED by minimising the χ 2 between observations and modelled SED. It is not possible to fix any of the parameters, all variables are free. As input data we used photometric observations from the archives of different observatories, covering the spectral range from the UV to the FIR. We used data from GALEX in the FUV and NUV bands, 2MASS in the J, H and K s bands, Spitzer/IRAC bands 1 − 4 and Herschel/PACS bands 70, 100, and 160 µm. The MSX filters are not available in MAGPHYS and hence are not used to model the SED. Details about the different photometric bands are reported in Table 4.
We obtained the flux density maps from the respective archives of the different observatories. To get only the flux density from the nucleus, we convolved the source size of 6.68 (see Fig. 6) with the resolution of the different observatories at the given wavelength and extracted the flux density from an aperture of the resulting diameter. We assume that the nuclear emission is dominant in this region and that emission from other structures like the bar are negligible. The obtained measurements are listed in Table 4 and shown in Fig. 8. Using these flux densities as input, MAGPHYS determines the best-fitting SED and returns the fitted parameters. From the fit, we obtained a total infrared luminosity L TIR = 1.37 × 10 10 L . We also used the prescription given in Galametz et al. (2013) (see their Eq. 7 and Table  3) to calculate L TIR from FIR data only. From this prescription -using the Herschel/PACS photometric data only -we obtain L TIR = (9.2 ± 0.4) × 10 9 L , which is in good agreement with the result from MAGPHYS. MAGPHYS also returns the optical depth, which is A V = 4.35 mag. This again is in good agreement with a recent study (Pérez-Beaupuits et al. 2018) who obtained A v = 5.5 ± 2.5 mag from the 100 µm dust continuum emission. The fitted SED is shown in Fig. 8.
Using the models from Weingartner & Draine (2001) with R V = 3.1 and the optical depth we obtained (A V = 4.35 mag) we calculated optical depths at other wavelengths. We yield an optical depth at λ = 20 µm of A 20 µm = 0.13 mag, which translates into an extinction correction for the line fluxes of ∼ 12%. The wavelengths of the emission lines used for the line flux ratios are close, thus, the extinction affects the line ratios even less. In case of the ([Ne ii] 13 µm + [Ne iii] 16 µm)/Hu α, the extinction correction is 3.5%. We find that the extinction correction has little effect for the set of lines used in our study.

Ionised gas density
Densities of the nuclear region of NGC 253 have been reported in the literature to range from ∼ 4.3 × 10 2 cm −3 (Carral et al. 1994, from [O iii]) up to 5 × 10 3 cm −3 (Engelbracht et al. 1998, from [Fe ii]) for the ionised gas. The results depend on the tracers A&A proofs: manuscript no. NGC253 Table 3. Line fluxes and uncertainties for the nuclear region, as well as properties of the emission lines observed with Spitzer/IRS (Fig. 7), SOFIA/FIFI-LS (Fig. 3)  used which could, in principle, apply to different regions within the observing beam or could be associated with different regions due to poorer spatial resolution. of the ionised gas, since they are independent of the ionic abundance and weakly dependent on the electron temperature. Each of these line flux ratios are not only ideal tracers of the density, due to their different critical densities (see Table 3), they potentially also trace different regions of the ionised ISM, as shown in Fig. 9. Owing to the critical densities, the line flux ratio [N ii] 122 µm/[N ii] 205 µm is a good tracer for the diffuse outskirts of H ii regions (Oberst et al. 2011  thermore, the [Ne v] 14 µm/[Ne v] 24 µm traces the gas phase whose conditions are dominated by the AGN, since Ne 4+ is not produced by stars but by the accretion disk from an AGN or fast radiative shocks (Izotov et al. 2012). We used a Monte-Carlo method to obtain a density from these line flux ratios and also to propagate uncertainties. With the Monte-Carlo method we were also able to investigate if the influence of the electron temperature is indeed negligible.
To determine the hydrogen density, we used the Python package PyNeb (Luridiana et al. 2013). PyNeb calculates the density of the ionised gas, if two density dependent emission lines and the temperature of the gas are provided, by solving the equilibrium equations of the n-level atom and determining the level populations. We assumed a normal distribution for the line fluxes, where the mean of the distribution is given by the line flux and the standard deviation corresponds to the statistical uncertainties given in Table 3. For the calculations, the temperature is uniformly distributed between 7000 K and 13000 K, which are typical values for the gas in H ii regions (Osterbrock & Ferland 2006). In each of the 10 5 steps, a temperature and line flux for each of the ionised species of a pair of lines were randomly picked and a density was calculated. The results of these calculations are shown in the left panel of Fig. 10  within a range where the line flux ratio is not sensitive to the density, hence, we only give an upper limit for the density. For the latter line flux ratio, we used both lines from the FIFI-LS observations to reduce systematic effects.
The densities derived from all line flux ratios are significantly lower than results from other studies (e.g. Engelbracht et al. 1998). For the [N ii] line flux ratio, this is most likely due to the low critical densities of both transitions. Both, the [N ii] 122 µm and [N ii] 205 µm lines are hence associated with diffuse, low density regions and do not trace high density regimes. The [S iii] line flux ratio traces higher density regions (see Fig. 9), but the S + ion has a higher ionisation potential (see Table 3) compared to, for example, Fe + . Fe + , however was used in Engelbracht et al. (1998). The different excitation potential and critical density of the species could explain the discrepancy between their and our study. The ratio [O iii] 52 µm/[O iii] 88 µm was also used by Carral et al. (1994). The different results from Carral et al. (1994) and this study is most likely due to the different sizes of observing beams, leading to a higher [O iii] 52 µm line flux and, hence, to a higher density, as discussed in Sect. 2.3.1. We conclude that, although all densities are in moderate agreement with each other, to appropriately model the ionised gas in the nuclear region at least a two phase model would be needed, due to the different critical densities and ionisation potentials of the used species.
We calculated the metallicity with the available emission lines, self-consistently. The [Ne ii] 13 µm, [Ne iii] 16 µm, and Hu α together provide excellent constraints for this purpose, by using the ([Ne ii] 13 µm+[Ne iii] 16 µm)/Hu α line flux ratio, which is sensitive to the Ne/H abundance in the gas (Bernard Salas et al. 2001;Lebouteiller et al. 2008). The lines arise from the ionised gas and Ne is not depleted on dust grains, hence, the Ne/H is an ideal tracer for the metallicity of the gas.
To calculate the Ne/H abundance ratio, we again used the Python tool PyNeb. For a given temperature T , density n, and line flux ratio [Ne ii] 13 µm/Hu α, PyNeb predicts the Ne + /H abundance ratio by solving the equilibrium equations and level populations for the Ne and H atom. The procedure is the same for the line flux ratio [  1.00 ± 0.24. The similar results suggest that the density does not seem to play a major role in the metallicity determination in the examined density regimes. From the Ne/H abundance, we use equations given by Izotov & Thuan (1999); Izotov et al. (2006); Nicholls et al. (2017), who studied elemental abundances in galaxies, to deduce abundances for C, N, O, Si, S, and Fe. The abundances of these elements are given in Table 6.
As discussed in, for instance, Madden et al. (2020), the metallicity is of great importance for determining the CO-to-H 2 conversion factor α CO . Using the results from previous studies, α CO (in M pc −2 K km s −1 −1 ) ranges from ∼ 1 up to ∼ 40, covering more than 1.5 orders of magnitude. From our obtained metallicity and Eq. (6) in Madden et al. (2020), we yield a more accurate α CO of 3.8 +5.8 −2.0 .

Summary and outlook
In this first of a series of papers we presented a wealth of MIR and FIR data of the nuclear region of the starburst galaxy NGC 253 from the SOFIA, Herschel, and Spitzer telescopes. We described the data reduction steps for each observatory and noted caveats that appeared during these processes. The full list of observed lines, the line fluxes and errors from the nuclear region are listed in Table 3, the extended line flux maps we obtained from SOFIA/FIFI-LS observations are shown in Fig. 3. Furthermore, we used the obtained line fluxes to derive fundamental parameters of the local ISM, such as density n, metallicity Z, and optical depth A V . We showed that for our observed MIR emission lines no extinction correction is needed. We used the line flux ratios of [S iii], [O iii], and [N ii], which are almost independent of the local temperature and mostly depend on the gas density or, more precisely, on the local density of electrons, which are the main collisional partners to excite these states. We conclude that there are at least two ionised gas phases in the central region of this galaxy. From the calculated density and assumed temperature, we used the line flux ratio of ([Ne ii] 13 µm + [Ne iii] 16 µm)/Hu α which is an ideal tracer for the Ne/H abundance, to calculate the local metallicity, for which we obtain solar value. In a forthcoming paper we will use all observed lines to create a more sophisticated multi-phase model using MULTIGRIS, a Bayesian approach to determine the probability density distribution within a set of Cloudy models Ramambason et al. 2022) of ISM parameters such as density, intensity of radiation field and age of the starburst.