Infrared view of the multiphase ISM in NGC 253

André Beck, Vianney Lebouteiller, Suzanne Madden, Christof Iserlohe, Alfred Krabbe, Lise Ramambason, Christian Fischer, Maja Kaźmierczak-Barthel, Serina Latzko, Juan-Pablo Pérez-Beaupuits

► To cite this version:

André Beck, Vianney Lebouteiller, Suzanne Madden, Christof Iserlohe, Alfred Krabbe, et al. Infrared view of the multiphase ISM in NGC 253: I. Observations and fundamental parameters of the ionised gas. Astronomy and Astrophysics - A&A, 2022, 10.1051/0004-6361/202243822. hal-03712457

HAL Id: hal-03712457

https://hal.science/hal-03712457

Submitted on 28 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Infrared view of the multiphase ISM in NGC 253

## I. Observations and fundamental parameters of the ionised gas

A. Beck<fnref n="1" />,<fnref n="2" />, V. Lebouteiller<fnref n="2" />, S. C. Madden<fnref n="2" />, C. Iserlohe<fnref n="1" />, A. Krabbe<fnref n="1" />, L. Ramambason<fnref n="2" />, C. Fischer<fnref n="1" />, M. Kaźmierczak-Barthel<fnref n="1" />, S. T. Latzko<fnref n="1" />, and J. P. Pérez-Beaupuits<fnref n="3" />,<fnref n="4" />

Received 20 April 2022 / Accepted 28 June 2022

## ABSTRACT

*Context.* Massive star formation leads to enrichment of the interstellar medium with heavy elements. On the other hand, the abundance of heavy elements is a key parameter with which to study the star-formation history of galaxies. Furthermore, the total molecular hydrogen mass, usually determined by converting CO or [C II] 158 µm luminosities, depends on the metallicity as well. However, the excitation of metallicity-sensitive emission lines depends on the gas density of the H II regions where they arise.

*Aims.* We used spectroscopic observations of the nuclear region of the starburst galaxy NGC 253 from SOFIA, Herschel, and Spitzer, 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 whether or not these lines could originate from different regions within the observing beam. The ([Ne II] 13 µm + [Ne III] 16 µm)/Hα line flux ratio on the other hand is independent of the depletion onto dust grains but sensitive to the Ne/H abundance ratio and is 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⁻<fnref n="3" /> and ~170–212 cm⁻<fnref n="3" />) and a metallicity of solar value.

**Key words.** galaxies: starburst – galaxies: individual: NGC 253 – infrared: ISM

## 1. Introduction

NGC 253 is the nearest spiral galaxy showing a nuclear starburst with a star-formation rate (SFR) of $3 M_{\odot} yr^{-1}$ (Radovich et al. 2001). Its proximity makes this galaxy a perfect model with which to study the interstellar medium (ISM) in the extreme environments of massive star formation (see Table 1 for basic properties of NGC 253). Near-infrared (NIR) 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\alpha$ (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_{\odot} yr^{-1}$ was estimated, implying that the nuclear starburst is severely suppressed by the outflow. The existence of a central supermassive black hole (SMBH) has been discussed within the literature (e.g., Vogler & Pietsch 1999; Günthardt et al. 2015). The putative central SMBH 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, 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 reduce the cold molecular gas reservoir, which is needed to sustain the ongoing starburst. However, the nature of this hypothetical BH, 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 these 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 BH or AGN (Fernández-Ontiveros et al. 2009) – the major energy input to the ISM is uncertain, and this series tackles this topic.

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. However, the metallicity is a key parameter when trying to understand the ongoing mechanisms in the ISM. As metals are byproducts of star formation, the metallicity of a galaxy tracks its evolutionary history. Additionally, to estimate the molecular hydrogen gas mass reservoir of galaxies from CO emission, the conversion

<footnote n="1">¹ Deutsches SOFIA Institut, University of Stuttgart, Pfaffenwaldring 29, 70569 Stuttgart
e-mail: abeck@dsi.uni-stuttgart.de</footnote>

<footnote n="2">² Laboratoire AIM, CEA/Service d'Astrophysique, Bât. 709, CEA-Saclay, 91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="3">³ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote n="4">⁴ Centro de Astro-Ingeniería (AIUC), Pontificia Universidad Católica de Chile, Av. Vicuña Mackena, 4860 Macul, Santiago, Chile</footnote>
<table><caption>**Table 1.** General properties of NGC 253.</caption><thead><tr><th>Quantity</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td>RA (J2000)</td><td>0h47m33.1s</td><td>Müller-Sánchez et al. (2010)</td></tr><tr><td>Dec (J2000)</td><td>-25°17′18″.3</td><td>Müller-Sánchez et al. (2010)</td></tr><tr><td>Distance</td><td>3.5 Mpc</td><td>Rekola et al. (2005)</td></tr><tr><td>Image scale</td><td>17 pc/″</td><td></td></tr><tr><td>SFR</td><td>3 M<sub>☉</sub> yr<sup>-1</sup></td><td>Radovich et al. (2001)</td></tr><tr><td>Inclination</td><td>78°</td><td>Pence (1981)</td></tr><tr><td>v<sub>sys</sub></td><td>236 km s<sup>-1</sup></td><td>Pence (1981)</td></tr></tbody></table>

**Notes.** Right ascension (RA) and declination (Dec) of the nucleus. Distance calculated from planetary nebula luminosity function. The SFR was calculated from ISO FIR observations, $v_{sys}$ is the heliocentric systemic velocity of the source.

factor is known to be a function of metallicity (e.g., Bolatto et al. 2013b, and references therein). Recent studies showed that, in low-metallicity environments in particular, a large fraction of H₂ is not traced by CO emission (e.g., Madden et al. 2020, and references within). Hence, knowing the metallicity is the main parameter to estimate the molecular hydrogen mass, which itself is a key ingredient of star formation.

This study uses recent SOFIA observations of [O III]52, 88 µm, [O I]63, 146 µm, and [C II] 158 µm in conjunction with *Spitzer* and *Herschel* MIR and FIR spectroscopy. We gather MIR to FIR observations from the various telescopes, some observations of which will be used in this study. In this paper, we address the properties of the ionised gas in the central region of NGC 253 using Hu $\alpha$ (7–6), [Ne II] 13 µm, [Ne III] 16 µm, [S III] 19 µm, [S III] 33 µm, [O III] 52 µm, [O III] 88 µm, [N II] 122 µm, and [N II] 205 µm. Due to their wide range of critical densities and ionisation potentials, we are able to pin down the densities traced by these different emission lines. Additionally, we quantify the abundances of various species, and thus the overall metallicity of the central region, since ratios like O/H and N/O provide measures for the star-formation history of a galaxy (e.g., Pérez-Montero et al. 2013; Spinoglio et al. 2022). In a followup paper, we present a more complex multiphase model using MULTIGRIS (Lebouteiller & Ramambason 2022; 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 the intensity of the radiation field of the different components of the ISM, the temperature and luminosity of the putative black hole or AGN, and the age of the starburst.

This work is structured as follows: Sect. 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, and determination of 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.

## 2. Observations and data reduction

### 2.1. SOFIA/FIFI-LS

We used the Field-Imaging Far-Infrared Line Spectrometer (FIFI-LS Fischer et al. 2018) on board the Stratospheric

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Photometric image from *Spitzer/IRAC* at $\lambda = 8$ µm in logscale. The white contours show observations from SOFIA/FIFI-LS (solid – [O I] 146 µm and [C II] 158 µm, dashed – [O III] 88 µm, dotted – [O III] 52 µm and [O I] 63 µm), the two black boxes show apertures of *Spitzer/IRS* short- and long high modules. The red box shows the approximate footprint of the *Herschel/PACS* observations.</figcaption></figure>

Observatory For Infrared Astronomy (SOFIA, Young et al. 2012) to 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 $\lambda = 51–120$ µm (blue channel) and in the $\lambda = 115–203$ µm (red channel) wavelength regimes. Each channel projects an array of $5 \times 5$ spatial pixels (spaxels) on the sky, covering a total area of $30'' \times 30''$ and $60'' \times 60''$ in the blue and red channels, 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).

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'' \times 1.5''$ and $3'' \times 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 for more details). For our [O III] 52 µm, [O III] 88 µm, [O I] 146 µm, and [C II] 158 µm observations, we used the data from the pipeline that is PWV calibrated. The [O I] 63 µm lies in a wavelength regime that is contaminated with a telluric absorption feature and is therefore not used in this study.

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
<table><caption>**Table 2.** Summary of our observed lines with FIFI-LS.</caption><thead><tr><th>Species</th><th>Flight dates</th><th>R</th><th>PSF ["]</th></tr></thead><tbody><tr><td>[O III] 52 µm</td><td>2015 Oct. 14, 2018 Nov. 07<br>2018 Nov. 08</td><td>980</td><td>5.5</td></tr><tr><td>[O I] 63 µm</td><td>2015 Oct. 14, 2015 Oct. 23</td><td>1300</td><td>6.0</td></tr><tr><td>[O III] 88 µm</td><td>2015 Oct. 23, 2018 Nov. 07<br>2019 Nov. 02, 2019 Nov. 05</td><td>600</td><td>9.0</td></tr><tr><td>[O I] 146 µm</td><td>2015 Oct. 23, 2018 Nov. 07<br>2018 Nov. 08, 2019 Nov. 05</td><td>1000</td><td>15.0</td></tr><tr><td>[C II] 158 µm</td><td>2015 Oct. 14, 2015 Oct. 23<br>2019 Nov. 02</td><td>1150</td><td>16.0</td></tr></tbody></table>

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).

<figure><img src="image_3.png" /></figure>

**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.

(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, varying the flux density within the normal-distributed flux den-

sity uncertainties in each spectralbin 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 pipeline 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.

## 2.2. Additional data

### 2.2.1. Spitzer/IRS

We complement our data set with single slit spectra obtained with the Infrared Spectrograph (IRS, Houck et al. 2004) on board the *Spitzer* Space Telescope. Data are available at the Combined Atlas of *Spitzer* IRS Spectra (Lebouteiller et al. 2015, CASSIS<fnref n="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, because the low-resolution observations are saturated. Both high-resolution modules have a velocity resolution of $\Delta v = 500 \text{ km s}^{-1}$, comparable to the spectral resolution of our FIFI-LS observations. Slit sizes of the short- and 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. As the nucleus of NGC 253 is partly extended (see Sect. 2.3.2 and Fig. 6), we chose a full extraction procedure as this ensures we retrieve most of the flux of the source. 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–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⁻<fnref n="1" /> with respect to the expected velocity from the source (see Table 1). The minimum line width is 300 km s⁻<fnref n="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 two times the RMS of the flux density of the continuum (in W m⁻² µm⁻<fnref n="1" />) multiplied by the spectral resolution (in µm). In this case, the upper limit is given by two times the rms error of the continuum.

### 2.2.2. Herschel

We add to the FIFI-LS and *Spitzer* data observations with the Photodetector Array Camera and Spectrometer (PACS, Poglitsch et al. 2010) on board the *Herschel* Space

<footnote n="1">¹ https://cassis.sirtf.com/</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Line flux maps from the SOFIA/FIFI-LS observations. Colour maps show the line flux in $10^{-17}$ W m$^{-2}$ on a 3" × 3" (left column) and 1.5" × 1.5" (right column) grid, respectively. Pixels with a S/N 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).</figcaption></figure>

Observatory<fnref n="2" /> (Pilbratt et al. 2010). We used PACSman v3.63 (Lebouteiller et al. 2012) 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 providing the line flux and line flux error in one spaxel, respectively.

We obtained line flux maps and corresponding uncertainties for the [N III] 57 µm, [O I] 63 µm, [O III] 88 µm, [N II] 122 µm, [O I] 146 µm, and [C II] 158 µm transitions. See 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), and therefore we do not reduce these observations and list only the results in Table 3. For details of the data reduction, see Pérez-Beaupuits et al. (2018).

## 2.3. Data homogenisation

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

### 2.3.1. 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 disc. 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 and avoid contamination from background (i.e., disc-) emission, we extracted the line flux and errors from the resampled maps using a circular aperture with a radius of two times 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, because 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 the uncertainties from Herschel/PACS (see Table 3). The resulting spectra from this procedure are shown in Fig. 5.

All line fluxes and corresponding statistical errors from the FIFI-LS observations are listed in the second block of Table 3. The line fluxes listed are in good agreement with the values from Herschel/PACS as explained in Sect. 2.3.3. We compared the [O III] 52 µm line flux with observations from the Kuiper Airborne Observatory (KAO), as this emission line was not observable with PACS. Compared to observations published by Carral et al. (1994), our value for [O III] 52 µm emission is lower by ~50%. However, the beam size of KAO at 52 µm is

<footnote n="2">² Obs. ID 1342210652, 1342212531</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Spectra of the nuclear region observed with the Herschel/PACS spectrometer. Only the spectrum of the central 3 × 3 spaxels is shown here (Sect. 2.3.3). The dashed line shows the redshift-corrected rest wavelength of the emission line.</figcaption></figure>

approximately eight times larger than that of SOFIA at that wavelength, so the KAO observations most likely include some background emission. To confirm that the KAO measurements of [O III] 52 µm include diffuse emission, we extracted the spectrum from a 43" aperture, which is the beam size of the KAO observations 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⁻². This is in good agreement with the findings of Carral et al. (1994), who yielded a line flux of 90 ± 20 W m⁻².

### 2.3.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 over-

lapping 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, and therefore 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<fnref n="3" />) at $\lambda$ = 12.13 µm, 14.65 µm, and 21.3 µm to assert that no wavelength-dependent scaling is necessary. The more recent WISE observations are saturated and therefore cannot be used to verify the scaling. The continuum of the IRS spectra agrees within ~10% with all three MSX bands, and therefore 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.

### 2.3.3. 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

<footnote n="3">³ Available at https://irsa.ipac.caltech.edu/data/MSX/</footnote>
<table><caption>**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), and *Herschel/PACS* and SPIRE.</caption><thead><tr><th>Species</th><th>Transition</th><th>Rest wavelength<br>[µm]</th><th>Flux<br>[10<sup>-16</sup> W m<sup>-2</sup>]</th><th>Critical density<br>[cm<sup>-3</sup>]</th><th>Ionisation energy<br>[eV]</th></tr></thead><tbody><tr><td colspan="6"><strong>Spitzer/IRS</strong></td></tr><tr><td>[S IV]</td><td><sup>2</sup>P<sub>3/2</sub> – <sup>2</sup>P<sub>1/2</sub></td><td>10.511</td><td>1.29 ± 0.21</td><td>5 × 10<sup>4</sup> [e]</td><td>34.79</td></tr><tr><td>H<sub>2</sub></td><td>0 – 0 S(2)</td><td>12.279</td><td>11.6 ± 1.06</td><td>2 × 10<sup>5</sup> [H<sub>2</sub>]</td><td>–</td></tr><tr><td>Hu α</td><td>n = 7 – 6</td><td>12.372</td><td>6.19 ± 1.58</td><td>–</td><td>–</td></tr><tr><td>[Ne II]</td><td><sup>2</sup>P<sub>1/2</sub> – <sup>2</sup>P<sub>3/2</sub></td><td>12.814</td><td>446.40 ± 44.64</td><td>7 × 10<sup>5</sup> [e]</td><td>21.56</td></tr><tr><td>[Ne V]</td><td><sup>3</sup>P<sub>2</sub> – <sup>3</sup>P<sub>1</sub></td><td>14.322</td><td>&lt;0.78</td><td>3 × 10<sup>4</sup> [e]</td><td>97.12</td></tr><tr><td>[Cl II]</td><td><sup>3</sup>P<sub>1</sub> – <sup>3</sup>P<sub>2</sub></td><td>14.368</td><td>5.52 ± 1.06</td><td>4 × 10<sup>4</sup> [e]</td><td>12.97</td></tr><tr><td>[Ne III]</td><td><sup>3</sup>P<sub>1</sub> – <sup>3</sup>P<sub>2</sub></td><td>15.555</td><td>29.91 ± 2.33</td><td>3 × 10<sup>5</sup> [e]</td><td>40.96</td></tr><tr><td>H<sub>2</sub></td><td>0 – 0 S(1)</td><td>17.035</td><td>11.93 ± 1.19</td><td>2 × 10<sup>4</sup> [H<sub>2</sub>]</td><td>–</td></tr><tr><td>[P III]</td><td><sup>2</sup>P<sub>3/2</sub> – <sup>2</sup>P<sub>1/2</sub></td><td>17.885</td><td>&lt;2.02</td><td>–</td><td>19.77</td></tr><tr><td>[Fe II]</td><td><sup>7</sup>D<sub>7/2</sub> – <sup>7</sup>D<sub>9/2</sub></td><td>17.936</td><td>&lt;3.10</td><td>–</td><td>7.90</td></tr><tr><td>[S III]</td><td><sup>3</sup>P<sub>2</sub> – <sup>3</sup>P<sub>1</sub></td><td>18.713</td><td>97.81 ± 3.78</td><td>2 × 10<sup>4</sup> [e]</td><td>23.34</td></tr><tr><td>[Ar III]</td><td><sup>3</sup>P<sub>0</sub> – <sup>3</sup>P<sub>1</sub></td><td>21.831</td><td>&lt;2.00</td><td>4 × 10<sup>4</sup> [e]</td><td>27.63</td></tr><tr><td>[Fe III]</td><td><sup>5</sup>D<sub>3</sub> – <sup>5</sup>D<sub>4</sub></td><td>22.925</td><td>10.30 ± 2.23</td><td>10<sup>5</sup> [e]</td><td>16.19</td></tr><tr><td>[Ne V]</td><td><sup>3</sup>P<sub>1</sub> – <sup>3</sup>P<sub>0</sub></td><td>24.318</td><td>&lt;1.50</td><td>3 × 10<sup>5</sup> [e]</td><td>97.12</td></tr><tr><td>[O IV]</td><td><sup>2</sup>P<sub>3/2</sub> – <sup>1</sup>P<sub>1/2</sub></td><td>25.913</td><td>7.75 ± 4.70</td><td>1 × 10<sup>4</sup> [e]</td><td>54.94</td></tr><tr><td>[Fe II]</td><td><sup>6</sup>D<sub>7/2</sub> – <sup>6</sup>D<sub>9/2</sub></td><td>25.988</td><td>21.00 ± 5.15</td><td>2 × 10<sup>6</sup> [H], 10<sup>4</sup> [e]</td><td>7.90</td></tr><tr><td>H<sub>2</sub></td><td>0 – 0 S(0)</td><td>28.221</td><td>&lt;5.00</td><td>700 [H<sub>2</sub>]</td><td>–</td></tr><tr><td>[S III]</td><td><sup>3</sup>P<sub>1</sub> – <sup>3</sup>P<sub>0</sub></td><td>33.480</td><td>166.00 ± 19.00</td><td>7 × 10<sup>3</sup> [e]</td><td>23.34</td></tr><tr><td>[Si II]</td><td><sup>2</sup>P<sub>3/2</sub> – <sup>2</sup>P<sub>1/2</sub></td><td>34.814</td><td>241.00 ± 22.40</td><td>3 × 10<sup>5</sup> [H], 10<sup>3</sup> [e]</td><td>8.15</td></tr><tr><td>[Ne III]</td><td><sup>3</sup>P<sub>0</sub> – <sup>3</sup>P<sub>1</sub></td><td>36.014</td><td>&lt;20.00</td><td>6 × 10<sup>4</sup> [e]</td><td>40.96</td></tr><tr><td colspan="6"><strong>SOFIA/FIFI-LS</strong></td></tr><tr><td>[O III]</td><td><sup>3</sup>P<sub>2</sub> – <sup>3</sup>P<sub>1</sub></td><td>51.815</td><td>43.34 ± 4.17</td><td>4 × 10<sup>3</sup> [e]</td><td>35.12</td></tr><tr><td>[O III]</td><td><sup>3</sup>P<sub>1</sub> – <sup>3</sup>P<sub>0</sub></td><td>88.356</td><td>87.60 ± 16.91</td><td>510 [e]</td><td>35.12</td></tr><tr><td>[O I]</td><td><sup>3</sup>P<sub>0</sub> – <sup>3</sup>P<sub>1</sub></td><td>145.525</td><td>43.01 ± 3.72</td><td>10<sup>5</sup> [H]</td><td>–</td></tr><tr><td>[C II]</td><td><sup>3</sup>P<sub>3/2</sub> – <sup>3</sup>P<sub>1/2</sub></td><td>157.741</td><td>452.95 ± 13.51</td><td>3 × 10<sup>3</sup> [H], 50 [e]</td><td>11.26</td></tr><tr><td colspan="6"><strong>Herschel/PACS</strong></td></tr><tr><td>[N III]</td><td><sup>2</sup>P<sub>3/2</sub> – <sup>2</sup>P<sub>1/2</sub></td><td>57.32</td><td>57.76 ± 40.73</td><td>3 × 10<sup>3</sup> [e]</td><td>29.60</td></tr><tr><td>[O I]</td><td><sup>3</sup>P<sub>1</sub> – <sup>3</sup>P<sub>2</sub></td><td>63.185</td><td>372.50 ± 57.59</td><td>5 × 10<sup>5</sup> [H]</td><td>–</td></tr><tr><td>[O III]</td><td><sup>3</sup>P<sub>1</sub> – <sup>3</sup>P<sub>0</sub></td><td>88.356</td><td>74.25 ± 17.19</td><td>510 [e]</td><td>35.12</td></tr><tr><td>[N II]</td><td><sup>3</sup>P<sub>2</sub> – <sup>3</sup>P<sub>1</sub></td><td>121.8</td><td>94.16 ± 22.28</td><td>300 [e]</td><td>14.53</td></tr><tr><td>[O I]</td><td><sup>3</sup>P<sub>0</sub> – <sup>3</sup>P<sub>1</sub></td><td>145.525</td><td>43.06 ± 7.104</td><td>10<sup>5</sup> [H]</td><td>–</td></tr><tr><td>[C II]</td><td><sup>3</sup>P<sub>3/2</sub> – <sup>3</sup>P<sub>1/2</sub></td><td>157.741</td><td>380.30 ± 44.61</td><td>3 × 10<sup>3</sup> [H], 50 [e]</td><td>11.26</td></tr><tr><td colspan="6"><strong>Herschel/SPIRE</strong></td></tr><tr><td>[N II]</td><td><sup>3</sup>P<sub>1</sub> – <sup>3</sup>P<sub>0</sub></td><td>205.3</td><td>20.63 ± 2.13</td><td>45 [e]</td><td>14.53</td></tr><tr><td>[C I]</td><td><sup>3</sup>P<sub>2</sub> – <sup>3</sup>P<sub>1</sub></td><td>370.415</td><td>12.25 ± 1.25</td><td>2 × 10<sup>3</sup> [H]</td><td>–</td></tr><tr><td>[C I]</td><td><sup>3</sup>P<sub>1</sub> – <sup>3</sup>P<sub>0</sub></td><td>609.135</td><td>4.93 ± 0.56</td><td>300 [H]</td><td>–</td></tr></tbody></table>

**Notes.** See Sects. 2.1 and 2.2 for data reduction, and Sect. 2.3 for extraction of the nuclear emission. Critical densities from Malhotra et al. (2001), Giveon et al. (2002), Kaufman et al. (2006), Meijerink et al. (2007) and Cormier et al. (2012), calculated at T = 100, 300 and 10 000 K for collisions with H atoms, H<sub>2</sub> molecules and electrons, respectively. The critical densities are listed for the main collisional partners. Uncertainties listed are only statistical errors from the fit. In addition, systematic uncertainties are 15% for IRS (Sloan et al. 2015) and PACS (Pérez-Beaupuits et al. 2018), and 20% for FIFI-LS (FIFI-LS team, priv. comm.). Transitions and wavelengths from NIST (Kramida et al. 2020), ionisation potentials to create the ion from Draine (2011).

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 disc. Line fluxes and errors of the brightest fine-structure lines observed with PACS are listed in Table 3.

### 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.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Spectra from our FIFI-LS observations in the nuclear region of NGC 253. See Sect. 2.3.1 and Fig. 3 for details of the extraction aperture.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Spatial profiles from the cross-dispersion direction of the LH module of IRS (see Lebouteiller et al. 2015, for details). Black shows the model PSF from observations of several point sources. Green diamonds show the spatial profile of NGC 253 in the observations, red is a Gaussian fit with a free FWHM. From the fit, we obtained a broadening of the PSF of 6.68''.</figcaption></figure>

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²⁺. This is particularly important, because the elemental abundances are uncertain for NGC 253.

### 3.1. Optical depth

Luminous IR objects tend to be highly extincted sources, making it sometimes necessary to correct for extinction even at MIR wavelengths, particularly for $A_V \ge 10$ mag (Armus et al. 2007). We therefore 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), to $A_V = 14.0$ mag (Rieke et al. 1980, from $Br\alpha/Bry$), and up to $A_V = 19.1$ mag (Engelbracht et al. 1998, from [Fe II]), depending on the tracers.

To determine the optical depth of the nuclear region of NGC 253, we used the SED fitting tool MAGPHYS (da Cunha et al. 2008). Based on a set of photometric observations between the UV and FIR range as input, MAGPHYS fits the SED of a galaxy. MAGPHYS varies different parameters for the stellar 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 $\chi^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 far-UV (FUV) and near-UV (NUV) bands, 2MASS in the *J*, *H*, and *Ks* bands, Spitzer/IRAC bands 1–4 and Herschel/PACS bands 70, 100, and 160 µm. The MSX filters are not available in MAGPHYS and are therefore 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
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** *Spitzer/IRS* spectrum (blue) and MSX broadband flux density measurements (red dots) of the nuclear region (11.1" × 22.3") of NGC 253. Detected emission lines in the *Spitzer* spectrum are labelled accordingly. MSX filter bandwidths are indicated in dashed lines.</figcaption></figure>

<table><caption>**Table 4.** Photometric bands, isophotal wavelengths $\lambda_{Isophot}$, and flux densities $F$ of the different observatories used for the SED fit with MAGPHYS (Sect. 3.1).</caption><thead><tr><th>Observatory</th><th>Band</th><th>λ<sub>Isophot</sub> [µm]</th><th>F [Jy]</th></tr></thead><tbody><tr><td rowspan="2">GALEX</td><td>FUV</td><td>0.153</td><td>0.03</td></tr><tr><td>NUV</td><td>0.253</td><td>0.10</td></tr><tr><td rowspan="3">2MASS</td><td>J</td><td>1.2</td><td>0.33</td></tr><tr><td>H</td><td>1.65</td><td>0.60</td></tr><tr><td>K<sub>s</sub></td><td>2.2</td><td>0.73</td></tr><tr><td rowspan="4">Spitzer/IRAC</td><td>1</td><td>3.6</td><td>0.78</td></tr><tr><td>2</td><td>4.5</td><td>0.87</td></tr><tr><td>3</td><td>5.8</td><td>3.10</td></tr><tr><td>4</td><td>8.0</td><td>8.83</td></tr><tr><td rowspan="3">Herschel/PACS</td><td>Blue</td><td>70.0</td><td>779.62</td></tr><tr><td>Green</td><td>100.0</td><td>822.10</td></tr><tr><td>Red</td><td>160.0</td><td>551.37</td></tr></tbody></table>

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, such as the bar, is 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 IR luminosity $L_{TIR} = 1.37 \times 10^{10} L_{\odot}$. 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 \pm 0.4) \times 10^9 L_{\odot}$, 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 the findings of a recent study (Pérez-Beaupuits et al. 2018)

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Broadband measurements from GALEX, 2MASS, *Spitzer/IRAC*, and *Herschel/PACS* observations and the resulting SED fit from MAGPHYS (blue). We also include the continuum from the SOFIA/FIFI-LS observations (Sect. 2.1) which are in good agreement with the fitted SED.</figcaption></figure>

who obtained $A_V = 5.5 \pm 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 $\lambda = 20$ µm of $A_{20\mu 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, and therefore the extinction affects the line ratios even less. In the case of the ([Ne II] 13 µm + [Ne III] 16 µm)/H$\alpha$, the extinction correction is 3.5%. We find that the extinction correction has little effect on the set of lines used in our study.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Density dependence (logscale) of the [O III] 52 µm/[O III] 88 µm, [N II] 122 µm/[N II] 205 µm, and [S III] 19 µm/[S III] 33 µm line flux ratios for different temperatures. The various line flux ratios are sensitive to different density regimes. The results from [N II] 122 µm/[N II] 205 µm and [S III] 19 µm/[S III] 33 µm are indicated in red and green dots, respectively, and for [O III] 52 µm/[O III] 88 µm we show the upper density limit (blue dots).</figcaption></figure>

## 3.2. Ionised gas density

Densities of the nuclear region of NGC 253 are reported in the literature to range from ~4.3×10² cm⁻³ (Carral et al. 1994, from [O III]) up to 5×10³ cm⁻³ (Engelbracht et al. 1998, from [Fe II]) for the ionised gas. The results depend on the tracers 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.

Several observed line flux ratios are sensitive to the density of the H II region. Well-known examples are the [O III] 52 µm/[O III] 88 µm and [N II] 122 µm/[N II] 205 µm ratios (e.g., Carral et al. 1994; Oberst et al. 2011; Herrera-Camus et al. 2016). Additionally, the MIR line flux ratios [Ne III] 16 µm/[Ne III] 36 µm, [Ne V] 14 µm/[Ne V] 24 µm, and [S III] 19 µm/[S III] 33 µm are excellent tracers for the density of the ionised gas because they are independent of the ionic abundance and are weakly dependent on the electron temperature. These line flux ratios are not only ideal tracers of the density – because of their different critical densities (see Table 3) –, but 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), the ratios [S III] 19 µm/[S III] 33 µm and [Ne III] 16 µm/[Ne III] 36 µm trace a higher density phase. Furthermore, the [Ne V] 14 µm/[Ne V] 24 µm traces the gas phase, the conditions of which are dominated by the AGN because Ne⁴⁺ is not produced by stars but by the accretion disc 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 whether or not 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⁵ 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 for the [S III] 19 µm/[S III] 33 µm line flux ratio.

We applied this Monte Carlo method to all line flux ratios listed above. However, from the [Ne III] 16 µm/[Ne III] 36 µm and [Ne V] 14 µm/[Ne V] 24 µm we did not obtain any results as these line flux ratios did not allow *PyNeb* to solve the level populations. We derived the densities listed in Table 5 from the line flux ratios.

The obtained densities from [S III] 19 µm/[S III] 33 µm and [N II] 122 µm/[N II] 205 µm are in good agreement with each other despite the different critical densities. In contrast, the [O III] 52 µm/[O III] 88 µm is lower by a factor of 2. The density obtained from the ratio [O III] 52 µm/[O III] 88 µm lies within a range where the line flux ratio is not sensitive to the density, and therefore we only give an upper limit for the density. For the
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Results from the Monte Carlo simulation to derive the hydrogen density (Sect. 3.2) and the Ne/H abundance ratio (Sect. 3.3). **Left:** main plot showing the derived density in 10 cm⁻³ bins. Inset plots show histograms of the randomly picked temperature (uniform distribution, 1000 K bins), and integrated line fluxes for [S III] 19 µm and [S III] 33 µm. **Right:** inset plots again showing histograms of the randomly picked integrated line fluxes of [Ne II] 13 µm, [Ne III] 16 µm, and Hu α. From the line fluxes, temperature, and density, we obtained the Ne/H abundance ratio in solar units, which is shown in the main plot (0.1 bins).</figcaption></figure>

<table><caption>**Table 5.** Densities derived from the line flux ratios as discussed in Sect. 3.2.</caption><thead><tr><td>Lines</td><td>Line flux ratio</td><td>Density [cm<sup>-3</sup>]</td></tr></thead><tbody><tr><td>[S III] 19/33 µm</td><td>0.59</td><td>170.0 ± 83.6</td></tr><tr><td>[O III] 52/88 µm</td><td>0.49</td><td>&lt;83.6</td></tr><tr><td>[N II] 122/205 µm</td><td>4.56</td><td>211.7 ± 105.6</td></tr></tbody></table>

**Notes.** The [O III] yields densities in the low-density regime, where the line ratio is not sensitive to the density, and so we give only an upper limit. The [Ne V] 14 µm, [Ne V] 24 µm, and [Ne III] 36 µm are upper limits and are therefore not suitable to determine a density.

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 therefore 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⁺. Nevertheless, Fe⁺ was used in Engelbracht et al. (1998). The different excitation potential and critical density of the species could explain the discrepancy between the results of our study and this latter. The ratio [O III] 52 µm/[O III] 88 µm was also used by Carral et al. (1994). The difference between the results of Carral et al. (1994) and those presented here 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, because of the different critical densities and ionisation potentials of the used species.

### 3.3. Chemical composition

The determinations of metallicity in the nuclear region of NGC 253 vary widely within the literature. The metallicity ranges from ~0.5 Z⊙ (Carral et al. 1994; Ptak et al. 1997), to ≤1 Z⊙ (Engelbracht et al. 1998), and up to ≤1.5 Z⊙ (Webster & Smith 1983).

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, and we can use 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; the Ne/H is therefore 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 [Ne III] 16 µm/Hu α, yielding the abundance ratio Ne²⁺/H. We assumed that Ne⁺ and Ne²⁺ are the dominant ionisation stages, that is Ne/H = Ne⁺/H + Ne²⁺/H. Under this assumption, we can use the densities derived in Sect. 3.2.
<table><caption>**Table 6.** Gas-phase elemental abundances relative to the hydrogen abundance for most of the species that were observed.</caption><thead><tr><th>Species</th><th>Abundance (log)</th><th>Ref.</th></tr></thead><tbody><tr><td>Ne/H</td><td>-4.000</td><td></td></tr><tr><td>C/H</td><td>-3.696</td><td>(3)</td></tr><tr><td>N/H</td><td>-4.341</td><td>(3)</td></tr><tr><td>O/H</td><td>-3.313</td><td>(2)</td></tr><tr><td>Si/H</td><td>-4.763</td><td>(1)</td></tr><tr><td>S/H</td><td>-5.053</td><td>(2)</td></tr><tr><td>Fe/H</td><td>-5.583</td><td>(2)</td></tr></tbody></table>

**Notes.** Abundances were calculated from the Ne/H abundance (see Sect. 3.3) and by the use of equations from (1) Izotov & Thuan (1999), (2) Izotov et al. (2006), and (3) Nicholls et al. (2017).

Again, we used the Monte Carlo method to propagate uncertainties in the derived Ne/H abundance ratio. We calculated the Ne/H abundance ratio assuming the different densities we obtained from the [O III], [N II], and [S III] line flux ratios in Sect. 3.2. To do so, we extended our Monte Carlo approach from Sect. 3.2. After calculating the density from the randomly drawn line fluxes (e.g., for [S III] 19 µm and [S III] 33 µm), line fluxes for [Ne II] 13 µm, [Ne III] 16 µm, and Hu α were drawn; the distributions for these lines are again normal distributions. With these three emission lines, as well as the temperature and density from the first step (Sect. 3.2), we calculated the Ne/H abundance ratio for 10⁵ steps. From this sample, we calculated the mean and standard deviation.

The results are shown in the right panels of Fig. 10. We obtain a Ne/H abundance ratio of nearly solar value, $[Ne/H]_{[O III]} = 1.08 \pm 0.41$, $[Ne/H]_{[N II]} = 1.00 \pm 0.24$, and $[Ne/H]_{[S III]} = 1.00 \pm 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), and 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 by Madden et al. (2020) for example, the metallicity is of great importance for determining the CO-to-H₂ conversion factor $\alpha_{CO}$. Using the results from previous studies, $\alpha_{CO}$ (in $M_{\odot} 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 $\alpha_{CO}$ of $3.8^{+5.8}_{-2.0}$.

## 4. Summary and outlook

In this first of a series of papers, we present 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 describe the data-reduction steps for each observatory and note 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, and 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_{\nu}$. We show 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 \mu m + [Ne III] 16 \mu m)/Hu \alpha$ – 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 (Lebouteiller & Ramambason 2022; Ramambason et al. 2022) of ISM parameters such as density, intensity of radiation field, and age of the starburst.

**Acknowledgements.** The authors thank the anonymous referee for useful comments that improved the quality and clarity of the paper. SOFIA, the Stratospheric Observatory For Infrared Astronomy, is a joint project of the Deutsche Raumfahrtagentur im Deutschen Zentrum für Luft- und Raumfahrt e.V. (German Space Agency at DLR, grant: FKZ 500K2002) and the National Aeronautics and Space Administration (NASA). It is funded on behalf of the German Space Agency at DLR by the Federal Ministry for Economic Affairs and Climate Action based on legislation by the German Parliament, the state of Baden-Württemberg and the University of Stuttgart. Scientific operation for Germany is coordinated by the Deutsches SOFIA Institut (DSI) at the University of Stuttgart. Software: Astropy (Astropy Collaboration 2013, 2018), Photutils (Bradley et al. 2021), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020)

## References

* Armus, L., Charmandaris, V., Bernard-Salas, J., et al. 2007, ApJ, 656, 148

* Astropy Collaboration (Robitaille, T. P., et al. 2013, A&A, 558, A33

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

* Bernard Salas, J., Pottasch, S. R., Beintema, D. A., & Wesselius, P. R. 2001, A&A, 367, 949

* Bernard-Salas, J., Spoon, H. W. W., Charmandaris, V., et al. 2009, ApJS, 184, 230

* Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013a, Nature, 499, 450

* Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013b, ARA&A, 51, 207

* Bradley, L., Sipőcz, B., Robitaille, T., et al. 2021, https://doi.org/10.5281/zenodo.596036

* Carral, P., Hollenbach, D. J., Lord, S. D., et al. 1994, ApJ, 423, 223

* Colditz, S., Beckmann, S., Bryant, A., et al. 2018, J. Astron. Instrum., 7, 1840004

* Cormier, D., Lebouteiller, V., Madden, S. C., et al. 2012, A&A, 548, A20

* da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

* Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton: Princeton University Press)

* Engelbracht, C. W., Rieke, M. J., Rieke, G. H., Kelly, D. M., & Achtermann, J. M. 1998, ApJ, 505, 639

* Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. Astron., 53, 385

* Fernández-Ontiveros, J. A., Prieto, M. A., & Acosta-Pulido, J. A. 2009, MNRAS, 392, L16

* Fernández-Ontiveros, J. A., Spinoglio, L., Pereira-Santaella, M., et al. 2016, ApJS, 226, 19

* Fischer, C., Bryant, A., Beckmann, S., et al. 2016, Proc. SPIE, 9910, 991027

* Fischer, C., Beckmann, S., Bryant, A., et al. 2018, J. Astron. Instrum., 7, 1840003

* Fischer, C., Iserlohe, C., Vacca, W., et al. 2021, PASP, 133, 055001

* Galametz, M., Kennicutt, R. C., Calzetti, D., et al. 2013, MNRAS, 431, 1956

* Giveon, U., Sternberg, A., Lutz, D., Feuchtgruber, H., & Pauldrach, A. W. A. 2002, ApJ, 566, 880

* Günhardt, G. I., Agüero, M. P., Camperi, J. A., et al. 2015, AJ, 150, 139

* Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

* Herrera-Camus, R., Bolatto, A., Smith, J. D., et al. 2016, ApJ, 826, 175

* Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18

* Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

* Iodice, E., Arnaboldi, M., Rejkuba, M., et al. 2014, A&A, 567, A86
Iserlohe, C., Fischer, C., Vacca, W. D., et al. 2021, PASP, 133, 055002

Izotov, Y. I., & Thuan, T. X. 1999, ApJ, 511, 639

Izotov, Y. I., Stasińska, G., Meynet, G., Guseva, N. G., & Thuan, T. X. 2006, A&A, 448, 955

Izotov, Y. I., Thuan, T. X., & Privon, G. 2012, MNRAS, 427, 1229

Kaufman, M. J., Wolfire, M. G., & Hollenbach, D. J. 2006, ApJ, 644, 283

Klein, R., Beckmann, S., Bryant, A., et al. 2014, Proc. SPIE, 9147, 91472X

Kramida, A., Ralchenko, Yu., Reader, J., & NIST ASD Team 2020, NIST Atomic Spectra Database (ver. 5.8), https://physics.nist.gov/asd

Krips, M., Martín, S., Sakamoto, K., et al. 2016, A&A, 592, L3

Lebouteiller, V., & Ramambason, L. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202243865

Lebouteiller, V., Bernard-Salas, J., Brandl, B., et al. 2008, ApJ, 680, 398

Lebouteiller, V., Cormier, D., Madden, S. C., et al. 2012, A&A, 548, A91

Lebouteiller, V., Barry, D. J., Goes, C., et al. 2015, ApJS, 218, 21

Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015, ApJ, 801, 25

Luridiana, V., Morisset, C., & Shaw, R. A. 2013, Astrophysics Source Code Library [record ascl:1304.021]

Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141

Malhotra, S., Kaufman, M. J., Hollenbach, D., et al. 2001, ApJ, 561, 766

Martín, S., Mangum, J. G., Harada, N., et al. 2021, A&A, 656, A46

Meijerink, R., Spaans, M., & Israel, F. P. 2007, A&A, 461, 793

Müller-Sánchez, F., González-Martín, O., Fernández-Ontiveros, J. A., Acosta-Pulido, J. A., & Prieto, M. A. 2010, ApJ, 716, 1166

Nicholls, D. C., Sutherland, R. S., Dopita, M. A., Kewley, L. J., & Groves, B. A. 2017, MNRAS, 466, 4403

Oberst, T. E., Parshley, S. C., Nikola, T., et al. 2011, ApJ, 739, 100

Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (University Science Books)

Pence, W. D. 1981, ApJ, 247, 473

Pérez-Beaupuits, J. P., Güsten, R., Harris, A., et al. 2018, ApJ, 860, 23

Pérez-Montero, E., Contini, T., Lamareille, F., et al. 2013, A&A, 549, A25

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Ptak, A., Serlemitsos, P., Yaqoob, T., Mushotzky, R., & Tsuru, T. 1997, AJ, 113, 1286

Radovich, M., Kahanpää, J., & Lemke, D. 2001, A&A, 377, 73

Ramambason, L., Lebouteiller, V., Bik, A., et al. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202243866

Rekola, R., Richer, M. G., McCall, M. L., et al. 2005, MNRAS, 361, 330

Rieke, G. H., Lebofsky, M. J., Thompson, R. I., Low, F. J., & Tokunaga, A. T. 1980, ApJ, 238, 24

Rosenberg, M. J. F., van der Werf, P. P., & Israel, F. P. 2013, A&A, 550, A12

Sloan, G. C., Herter, T. L., Charmandaris, V., et al. 2015, AJ, 149, 11

Spinoglio, L., Fernández-Ontiveros, J. A., Malkan, M. A., et al. 2022, ApJ, 926, 55

Vacca, W., Clarke, M., Perera, D., Fadda, D., & Holt, J. 2020, in ASP Conf. Ser., eds. R. Pizzo, E. R. Deul, J. D. Mol, J. de Plaa, & H. Verkouter, 527, 547

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Vogler, A., & Pietsch, W. 1999, A&A, 342, 101

Walter, F., Bolatto, A. D., Leroy, A. K., et al. 2017, ApJ, 835, 265

Webster, B. L., & Smith, M. G. 1983, MNRAS, 204, 743

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Young, E. T., Becklin, E. E., Marcum, P. M., et al. 2012, ApJ, 749, L17