X-Ray Properties of NGC 253's Starburst-driven Outflow

In order to constrain the properties from the hot gas outflow, spectra were extracted using ciao. We defined several regions along the minor axis of the outflow as a function of distance from the kinematic center of NGC 253 (Müller-Sánchez et al. 2010). The areas where spectra were extracted are shown in Figure 3. From the center to the north of the outflow, the regions were 025 in height and 1' in width, while the southern regions were 05 in height and 1' in width. The region sizes were set to achieve ≳5000 counts per region (in order to reliably constrain the metal abundances) and were placed qualitatively in different sections of the outflow (the starburst nucleus, north, and south).

Zoom In Zoom Out Reset image size

Background regions were defined for each of the observations and subtracted from each of the source regions. The background regions totaled an area of 32 arcmin². These regions were placed far from the disk in order to avoid diffuse emission from NGC 253, while also being located on the same ACIS chip as the source regions. The background regions are ≥25 from the center of the galaxy.

Once extracted and background-subtracted, the spectra were modeled with XSPEC (Arnaud 1996) Version 12.12. The XSPEC model that was used to fit the spectral data differed depending on the region. All the models included a multiplicative constant component (const), two absorption components (phabs*phabs), a power-law component (powerlaw), and at least one optically thin, thermal plasma component (vapec). The const component was allowed to vary and accounted for variations in emission between the observations. The first phabs component accounted for the galactic absorption in the direction of NGC 253, ${N}_{{\rm{H}}}^{\mathrm{MW}}=1.4\times {10}^{20}\ {\mathrm{cm}}^{-2}$ (Dickey & Lockman 1990), and was frozen. The second component ${N}_{{\rm{H}}}^{\mathrm{NGC}\,253}$ represented the intrinsic absorption of NGC 253 and was allowed to vary. Both phabs components have their metal abundance set to solar, consistent with observations of the NGC 253 nuclear region (Mills et al. 2021; Beck et al. 2022). The powerlaw component accounted for the nonthermal X-ray emission from, e.g., X-ray binaries, and the vapec component represented the thermal plasma with variable abundances and temperatures. We adopted solar abundances from Asplund et al. (2009) and photoionization cross sections from Verner et al. (1996).

Several other components were added as necessary to improve fits. The central region required a second vapec component that accounted for a hotter thermal plasma component. The central, S1, and S2 regions also included an AtomDB ¹⁴ CX component vacx that accounted for the emission from electrons captured by ions from neutral atoms (Smith et al. 2012). The vacx and additional vapec components have their metal abundances tied to the first vapec component. Regions N1 and N2 were found to not need the vacx component based on F-tests.

The final model used for regions N1 and N2 was const*phabs*phabs*(vapec+powerlaw). For regions S1 and S2, the model was const*phabs*phabs*(vapec+powerlaw+vacx), and for the central region, the model was const*phabs*phabs*(vapec+vapec+powerlaw+vacx).

X-ray emission was also evident in the disk of NGC 253. To constrain its properties, we extracted and modeled spectra from nine 1' by 1' regions shown in Figure 5. The disk spectral models were similar to those of the outflow except the abundances were not allowed to vary individually because there were not enough counts to reliably constrain them. Instead, the total metallicity of the apec component, Z, was set to solar to be consistent with observations (Mills et al. 2021; Beck et al. 2022). The XSPEC model for the disk was const*phabs*phabs*(apec+powerlaw).

In Figure 3, we present the X-ray spectra extracted from the outflow. We find several emission lines in these regions as identified in the spectrum of the central region. We detect emission lines from Ne, Fe, Mg, Si, S, and Ar. The Fe xxv, Ar xvii, and S xv lines are only apparent in the central region because of the hotter plasma located there.

Using the models described in Section 2 and the extracted spectra, we find the best-fit values of column density ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$, temperature kT, and metal abundances. We report these values in Table 2 and plot them in Figure 4. We find that the models fit the data well, with the reduced χ² values in each region being near one, and the largest being that for the central region—a reduced χ² of 1.28 and a null hypothesis of 7.67 × 10⁻¹³.

Zoom In Zoom Out Reset image size

Table 2. Spectral Fit Results ^a

Region	r	${N}_{{\rm{H}}}^{\mathrm{NGC}253}$	kT	O/O☉	Ne/Ne☉	Mg/Mg☉	Si/Si☉	S/S☉	Fe/Fe☉	χ2/dof
	(kpc)	(×1022 cm−2)	(keV)
N2	0.52	0.17 ± 0.09	0.71 ± 0.08	<0.96	<0.51	${0.45}_{-0.18}^{+0.45}$	${1.0}_{-0.27}^{+0.58}$	${1.0}_{-0.28}^{+0.76}$	${0.09}_{-0.03}^{+0.06}$	240/227
N1	0.26	${0.66}_{-0.09}^{+0.12}$	0.73 ± 0.04	1	1	${1.0}_{-0.17}^{+0.20}$	${1.4}_{-0.20}^{+0.23}$	${1.4}_{-0.39}^{+0.44}$	${0.05}_{-0.03}^{+0.04}$	551/498
Center b	0	0.86 ± 0.09	0.98 ± 0.02	${1.3}_{-0.57}^{+0.97}$	${0.89}_{-0.30}^{+0.50}$	${1.8}_{-0.55}^{+0.92}$	${4.0}_{-1.1}^{+1.9}$	${6.6}_{-1.8}^{+3.2}$	${2.7}_{-0.74}^{+1.4}$	1979/1552
S1	−0.39	0.15 ± 0.04	0.65 ± 0.03	${0.20}_{-0.08}^{+0.12}$	${0.69}_{-0.20}^{+0.27}$	${0.49}_{-0.12}^{+0.15}$	${0.56}_{-0.11}^{+0.14}$	${0.66}_{-0.28}^{+0.32}$	0.26 ± 0.04	752/636
S2	−0.92	${0.18}_{-0.09}^{+0.10}$	0.37 ± 0.04	${0.24}_{-0.10}^{+0.15}$	${0.60}_{-0.21}^{+0.31}$	${0.35}_{-0.13}^{+0.19}$	1	1	${0.21}_{-0.06}^{+0.08}$	272/258

Notes.

^a Abundances with values of 1 are frozen to solar values in the fits. ^b The temperature for the second, hotter thermal plasma component is ${{kT}}_{2}={5.5}_{-0.8}^{+1.2}$ keV.

Download table as:  ASCIITypeset image

All the best-fit values, with the exception of Ne, have distributions that peak in the center and decrease outward along the outflow to the north and south. For ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ the gradient is similar in the northern and southern outflows, with a peak of ${N}_{{\rm{H}}}^{\mathrm{NGC}253}=(8.6\pm 0.9)\times {10}^{21}$ cm⁻² in the center and decreasing to ${N}_{{\rm{H}}}^{\mathrm{NGC}253}\approx (1\mbox{--}2)\times {10}^{21}$ cm⁻². kT peaks at the center, with kT = 0.98 ± 0.02 keV and decreases outward along the outflow axis. We note that the second, hotter component in the center had a best-fit value of ${{kT}}_{2}={5.5}_{-0.8}^{+1.2}$ keV.

The metal abundances also peak in the center, with S having the highest value of ${6.6}_{-1.8}^{+3.2}$ for S/S_⊙ and Ne having the smallest peak with a relatively flat distribution within the uncertainties. We note in Section 4.2 that the high S value may be due to contributions from the second temperature component. The profiles show that there is little enhancement of the metals beyond the central region, with the outer regions containing metal abundances of around the solar value or lower. The Fe abundance is particularly low, with all the outflow regions below ≲0.30Fe/Fe_⊙. The decrease in abundances along the outflow may suggest mass loading or indicate that other physics should be considered, such as nonequilibrium ionization.

From the best-fit values, we calculate several other properties of the outflow for each region and present them in Table 3. The norm is defined as norm = (10⁻¹⁴EM)/4π D², where the emission measure is EM = ∫n_e n_H dV. By setting n_e = 1.2n_H, we calculate n_e as ${n}_{{\rm{e}}}={(1.5\,\times \,{10}^{15}\mathrm{norm}{D}^{2}/{fV})}^{1/2}$ where f is the filling factor and V is the volume. In computing V, we assume that each region has a cylindrical volume of radius R that we measure using broadband, X-ray brightness profiles. We make 2'' slits along the major axis to create the brightness profile, and we define R as the scale that encompasses 95% of the X-ray brightness (for reference, we also list the values corresponding to the 68% X-ray brightness in Table 3 as well). We note that the measurement of R is influenced by our choice of regions over which we extracted the profiles. We found that if we defined regions larger than the ones used in Section 2.1, the 95% emission radii were overestimated due to more background emission being included.

Table 3. Physical Parameters of the Central and Outflow Regions

Region	Norm	LX	LCX/LX	EM	Ra	Va	ne a	P/ka	tcool a
		(×1038 erg s−1)		(×1062 cm−3)	(×1021 cm)	(×1063 cm3)	(cm−3)	(×106 K cm−3)	(Myr)
N2	1.5 × 10−4	1.1	...	0.32	1.5 (0.68)	5.3 (3.0)	${0.07}_{-0.02}^{+0.04}$ (${0.09}_{-0.03}^{+0.05}$)	${1.2}_{-0.33}^{+0.67}$ (${1.5}_{-0.43}^{+0.88}$)	${48}_{-13}^{+28}$ (${36}_{-10}^{+21}$)
N1	6.3 × 10−4	5.2	...	1.6	1.3 (0.58)	3.9 (2.2)	0.17 ± 0.02 (0.22 ± 0.03)	${2.8}_{-0.33}^{+0.41}$ (${3.8}_{-0.44}^{+0.55}$)	${21}_{-2.4}^{+3.0}$ (${15}_{-1.8}^{+2.3}$)
Center b	6.6 × 10−4	47	0.42	1.4	0.73 (0.26)	1.3 (0.59)	0.30 ± 0.02 (${0.44}_{-0.03}^{+0.04}$)	${6.7}_{-0.52}^{+0.56}$ (${10}_{-0.78}^{+0.84}$)	${17}_{-1.4}^{+1.5}$ (${12}_{-0.91}^{+0.97}$)
S1	7.1 × 10−4	7.7	0.20	1.5	1.3 (0.63)	7.8 (4.7)	${0.13}_{-0.01}^{+0.02}$ (0.16 ± 0.02)	${1.9}_{-0.20}^{+0.24}$ (${2.5}_{-0.26}^{+0.30}$)	${25}_{-2.6}^{+3.1}$ (${19}_{-2.0}^{+2.4}$)
S2	2.5 × 10−4	3.3	0.33	0.53	1.4 (0.73)	9.1 (6.0)	${0.07}_{-0.01}^{+0.02}$ (0.09 ± 0.02)	${0.60}_{-0.12}^{+0.16}$ (${0.74}_{-0.15}^{+0.20}$)	${27}_{-5.2}^{+7.1}$ (${21}_{-4.2}^{+5.7}$)

Notes.

^a The values listed for R, V, n_e, P/k, and t_cool are for the 95% emission radius. The values in parentheses are for the 68% emission radius (see Section 3.1 for details of these measurements). ^b The parameters of the hotter component kT₂ in the central region are: norm₂ = 2.7 × 10⁻⁴, n_e,2 = 0.19 ± 0.02 cm⁻³, P₂/k = (2.4 ± 0.2) × 10⁷ K cm⁻³, and ${t}_{\mathrm{cool},2}={211}_{-19}^{+21}$ Myr. For the 68% emission radius, the values are: n_e,2 = 0.28 ± 0.03 cm⁻³, ${P}_{2}/k={3.6}_{-0.34}^{+0.36}\times {10}^{7}$ K cm⁻³, and ${t}_{\mathrm{cool},2}={140}_{-13}^{+14}$ Myr.

Download table as:  ASCIITypeset image

We assume a filling factor of f = 1, though we note that due to processes such as mass loading, it may be lower. The other measured quantities are the thermal pressure P/k = 2n_e T and the cooling time t_cool = 3kT/Λn_e. We find the radiative cooling function Λ using CHIANTI (Dere et al. 1997) at solar metallicity assuming a thin, thermal plasma.

The results in Table 3 show that R, and consequently the volume V, increases with distance from the center of the galaxy. n_e peaks in the central region with a value of n_e = 0.30 cm⁻³ and decreases along the outflow. The pressure also peaks in the central region with a value of 6.7 × 10⁶ K cm⁻³. The cooling time has the shortest value of t_cool = 17 Myr in the central region, and the longest value of 48 Myr in region N2. S1 and S2 have comparable cooling times of 25–27 Myr.

As noted in Section 4.4, when using smaller regions, the cooling times are much shorter (the shortest time is t_cool =1.26 Myr). The cooling times of larger regions appear to overestimate those of the smaller regions.

We also calculate the same physical quantities shown in Table 3 for the hotter component, kT₂, which is only present in the central region. We find that the hotter component has a lower n_e of 0.19 cm⁻³ than the cooler component. The hotter component has a higher pressure of 2.4 × 10⁷ K cm⁻³ and a longer cooling time of 211 Myr than the respective values in the cooler component.

Using the spectral model described in Section 2 we measure kT and ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ for the diffuse X-ray emission in the disk (Table 4). We find that the kT values are in the range 0.18–0.30 keV and the ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ values are in the range (0.31–0.72) ×10²² cm⁻². As shown in the right panel of Figure 5, we find the highest kT in Region 1 and the lowest in Region 4.

Zoom In Zoom Out Reset image size

Table 4. Disk Spectral Results

Region	${N}_{{\rm{H}}}^{\mathrm{NGC}253}$	kT	χ2/dof
	(×1022 cm−2)	(keV)
1	${0.40}_{-0.14}^{+0.19}$	${0.30}_{-0.06}^{+0.04}$	212/185
2	${0.37}_{-0.15}^{+0.16}$	${0.24}_{-0.02}^{+0.05}$	186/178
3	${0.72}_{-0.19}^{+0.08}$	${0.20}_{-0.02}^{+0.05}$	206/199
4	${0.71}_{-0.11}^{+0.10}$	${0.18}_{-0.01}^{+0.02}$	191/181
5	${0.54}_{-0.10}^{+0.09}$	${0.24}_{-0.01}^{+0.02}$	281/237
6	<0.53	${0.25}_{-0.03}^{+0.11}$	279/273
7	${0.69}_{-0.06}^{+0.05}$	${0.20}_{-0.01}^{+0.01}$	396/335
8	${0.31}_{-0.20}^{+0.29}$	${0.24}_{-0.05}^{+0.03}$	225/243
9	${0.64}_{-0.07}^{+0.07}$	${0.20}_{-0.01}^{+0.02}$	349/303

Download table as:  ASCIITypeset image

Previous X-ray studies of NGC 253 have been performed and are complementary to our analysis. Strickland et al. (2000) used 15 ks of Chandra X-ray observations to study NGC 253's outflow, much shallower than the 365 ks data considered in this work. In their work, X-ray spectra were extracted along the starburst outflow from regions spanning a similar area to those in this paper. There were five regions, each of which measured about 033 in height and 08 in width. The spectra were fitted using a single-temperature plasma model, where kT, ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$, and the normalization were allowed to vary. In their analysis, the model was only fit to the soft X-rays (0.3–3 keV) due to near-zero count rates at harder X-ray energies in most regions.

Strickland et al. (2000) found that the ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ values in the southern hemisphere range from 0.034 × 10²² to 0.056 ×10²² cm⁻² and kT values range from 0.46 to 0.63 keV. In the northern outflow, their ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ values range from 0.473 × 10²² to 0.922 × 10²² cm⁻² and their temperatures range from 0.59 to 0.66 keV.

Our kT values in the southern outflow are in a similar range (0.37–0.65 keV); but our ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ values differ greatly ((0.15–0.18) × 10²² cm⁻²). We find the opposite when comparing our northern outflow values, where the ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ values are comparable to ours ((0.17–0.66) × 10²² cm⁻²), but the kT results are discrepant (0.71–0.72 keV).

The difference between our derived values and those of Strickland et al. (2000) may result from their limited energy range (0.3–3 keV) and differences in the spectral model they adopted. As a result of Strickland et al. (2000) using a softer X-ray range, their derived temperatures are lower in regions of higher ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ (the northern outflow), where softer X-rays are attenuated. By contrast, in the southern outflow where ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ is lower, we find comparable kT values. If we adopt the energy range of Strickland et al. (2000), then we find lower temperatures in the northern outflow (0.66 ± 0.04 and ${0.71}_{-0.09}^{+0.08}$ keV), consistent with Strickland et al. (2000)'s values within the range of the errors.

As for the spectral fitting, Strickland et al. (2000) adopted a variable-abundance, single-temperature, MEKAL hot plasma model for the combined southern regions. Once they found best-fit metallicities, they adopted and froze these values when fitting the other regions while letting kT, ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$, and the normalization vary. Thus, by assuming a constant value for the metal abundances, Strickland et al. (2000) did not constrain the metallicity gradients along the outflow as we did in this work.

In Bauer et al. (2007), X-ray spectra from the XMM-Newton Reflection Grating Spectrometer (RGS) were extracted along the starburst outflow. Four of their five regions overlap with ours: two from the southern outflow, one from the center, and one from the northern outflow. They estimated the kT of each region using line ratios of Si, Mg, Ne, and O. They did not measure other parameters (e.g., ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ or metal abundances) using spectral fits because of the limited statistics of the data. In the northern outflow, they found that the temperatures range from 0.21 to 0.61 keV, in the center they range from 0.22 to 0.79 keV, and in the southern regions they range from 0.21 to 0.46 keV and 0.25–0.31 keV. Our temperatures are about a factor of 1.4 larger in the first southern region and about a factor of 1.2 larger in all the other comparable regions.

The discrepancies between Bauer et al. (2007) and this work may be a result of the different analysis performed. Rather than modeling the spectra of an entire region to find a best-fit temperature, Bauer et al. (2007) derived a temperature from the ratios of the same elements in different ionization states. They found that even within the same regions, the ratios yielded different temperatures. They noted that the different temperatures may be a result of the elements sampling different parts of the gas along the line of sight. If this is the case, then finding a best-fit temperature for the whole region would indeed result in a different value than temperatures from line ratios that sample different parts of the gas. Another physical explanation may be that there is a second temperature component in the southern regions. While we did not find a second temperature component to be required to improve fits there, it may still be physically present. If that is the case, then the second temperature component could be in better agreement with the results from Bauer et al. (2007).

Another reason for the discrepancy between our work and Bauer may be the different instruments used. As mentioned before, Bauer et al. (2007) used X-ray spectra from the XMM-Newton RGS. This instrument only covers the 0.33–2.5 keV energy range, ¹⁵ in contrast to our Chandra data that are in the range 0.5–7 keV. Their lower temperatures may be a result of them probing the softer X-ray regime. Another explanation is that Chandra is currently less sensitive to softer X-rays than it was when it originally launched due to contamination in the detectors (Plucinsky et al. 2003, 2018). About 60% of our data were taken well after the contamination began affecting Chandra's soft X-ray sensitivity. As a result, our temperature values may be skewed upward.

Mitsuishi et al. (2013) used X-ray spectra from XMM-Newton and Suzaku to derive kT values and metallicities from defined disk, superwind, and halo regions in NGC 253. Their superwind region was approximately the size and location of both of our southern regions, measuring 1.4 kpc × 0.9 kpc. They used an XSPEC model similar to ours, with two absorption components to account for the Milky Way's absorption and the intrinsic NGC 253 absorption. Different from our model, they opted for a two-temperature plasma (two vapec components) and incorporated a zbremms component to include contributions from point-source emission, rather than a power-law component, which is preferred for X-ray binaries (Lehmer et al. 2013). For their two temperatures, they found a warm component of 0.21 keV and a warm–hot component of 0.62 keV. These values are comparable to the range of our best-fit kT values. Their metallicities are also comparable to ours in the southern hemisphere within the range of the errors. However, their ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ measurement is lower than ours with an upper limit of 0.05 × 10²² cm⁻², while our value is (0.15–0.18) × 10²² cm⁻².

The analysis in this paper is similar to that performed on M82 in Lopez et al. (2020). As the two nearest starburst galaxies where this analysis has been completed, it is worthwhile to compare their results for a more complete picture of the hot phase in galactic outflows. Differences in their metal enhancement trends, for example, would motivate investigations into how the differences arise and whether they depend on galaxy properties.

We detect the same metals in NGC 253 that were detected in M82, and we find that their distributions vary across the outflow. Ne is elevated in M82's northern hemisphere. However, for NGC 253, Ne has a flatter distribution with neither the northern nor southern outflow containing an enhancement. Mg and Fe have flat distributions in M82, while in NGC 253 they peak at the center and decrease with distance from the disk. Both M82 and NGC 253 have elevated amounts of S and Si in the central region that decrease with distance. Also, for both galaxies, the O abundance is not well constrained in regions of high column density because soft X-rays are attenuated. Both M82 and NGC 253 have an enhanced α/Fe ratio in the outflow regions.

S and Si abundances peak in the central region of both galaxies and decrease with distance. However, in M82 the decrease in S and Si abundance is gradual, while in NGC 253 the decrease is sudden. As reported in Lopez et al. (2020), the S and Si line fluxes originate predominantly from the hot component. In our analysis of NGC 253, only the central region has this hotter component, which is the likely explanation for why the S and Si peaks are so prominent in that location.

Das et al. (2019, 2021) found α-enhancement in their analysis of the Milky Way's CGM. Gupta et al. (2022) also found supersolar Ne/O in the Galactic eROSITA bubbles. The origin of the α-enhancement or nonsolar abundance ratio was unclear, but it was speculated that the enhancement may be a product of galactic outflows.

In our analysis of the galactic wind of NGC 253, we find α-enhancement already in the outflow regions. We also find supersolar Ne/O ratios in the southern outflow regions. These elevated α/Fe ratios and supersolar Ne/O in the wind may provide an origin for the observed metal enrichment of the Milky Way CGM.

In order to provide a more detailed comparison to a theoretical galactic wind model, we take advantage of the smaller count threshold necessary to accurately constrain kT and ${N}_{{\rm{H}}}^{\mathrm{NGC}253}$ (compared to the signal necessary to measure the metal abundances). Specifically, we create 42 regions along the outflow for spectral extraction that are 25 in height along the minor axis and 1' in width along the major axis. They span the same distance from the starburst as the larger regions we describe in Section 2. We adopt an XSPEC model of const*phabs*phabs*(vapec+powerlaw) to fit the spectra associated with these regions. We freeze the metallicities in the vapec component to the best-fit values found in the associated larger regions. Additionally, we look at lateral surface brightness profiles perpendicular to the wind axis for each region and find 68% and 95% emission radii (as in Section 3.1). From these radii, we calculate n_e for each region.

We first compare the best-fit kT and n_e values to those expected in the galactic wind model of Chevalier & Clegg (1985, CC85). The CC85 model describes a supernova-powered galactic wind that assumes spherical outflow geometry and adiabatic expansion. This model provides a good starting point for comparison because of its simplicity.

The controlling parameters of the CC85 model are the starburst radius R and the total energy and mass-loading rates, ${\dot{E}}_{{\rm{T}}}=\alpha \times 3.1\times {10}^{41}\times ({\dot{M}}_{\mathrm{SFR}}/{M}_{\odot }\,{\mathrm{yr}}^{-1})\,[\mathrm{erg}\,{{\rm{s}}}^{-1}]$ and ${\dot{M}}_{{\rm{T}}}=\beta \times {\dot{M}}_{\mathrm{SFR}}$, respectively, where ${\dot{M}}_{\mathrm{SFR}}$ is the star formation rate, and α and β are dimensionless parameters. We take R = 200 pc (Strickland et al. 2000; Bauer et al. 2007) and ${\dot{M}}_{{\rm{SFR}}}=2.5\,{M}_{\odot }\,{{\rm{yr}}}^{-1}$ (Ott et al. 2005; Bendo et al. 2015; Leroy et al. 2015). For the 68th percentile contours, we find α₆₈ ≃ 0.32 and β₆₈ ≃ 0.55. For the 95th percentile contours these parameters are α₉₅ ≃ 0.21 and β₉₅ ≃ 0.36.

In Figure 6, we plot the observed kT and n_e values (black and purple dashed lines with "o" markers) with those of the spherical adiabatic CC85 model (gray lines). The disagreement between the CC85 model and the best-fit observed values primarily widens when the flow leaves the wind-driving region (i.e., when r > R). We see that the adiabatic spherical model produces temperature and density profiles that fall off faster than the observed values.

Zoom In Zoom Out Reset image size

This may be a result of missing processes such as nonspherical areal divergence (Nguyen & Thompson 2021). In NGC 253, the super star clusters within the nuclear core are likely arranged in a ring-like geometry (Levy et al. 2022). Nguyen & Thompson (2022) recently showed that energy and mass injection in a ring-like geometry naturally produces a biconical outflow, where the hot phase is cylindrically collimated. When we include the effect of nonspherical flow geometry into the wind model (red and blue lines in Figure 6), we find better agreement with the profiles.

The enclosed 68th and 95th percentile volumes, which were used to determine n_e, can also be used to define an outflow geometry A(s) = π r(s)², where r(s) is the radius of the 68th or 95th percentile volume at height s perpendicular to the wind axis. In Figure 7, we plot A(s) and $d\mathrm{ln}A/{ds}$, where the latter term describes the areal expansion rate. In the left panel of Figure 7, we see that the derived A(s) is highly nonspherical. In the right panel, we see that the derived A(s) (red and blue lines) undergoes an areal expansion rate slower than that of spherical (gray line). These profiles demonstrate that the wind does not have the geometry assumed by the CC85 model.

Zoom In Zoom Out Reset image size

In Figure 8 we show the calculated cooling times (t_cool) for 68% and 95% emission radii, along with the advection times (t_adv) from both the spherical (orange lines) and nonspherical (green lines) wind models described above. The advection times are calculated as t_adv = r/v, where r is the distance along the minor axis to the central starburst and v is the velocity of the wind as it advects, predicted from the models. As noted in Section 3, we assume a filling factor of f = 1, therefore the cooling times derived are upper limits. Using the 68% emission results, we find that the measured t_cool values are of the order of (but somewhat larger than) the predicted advection times from both the CC85 and A(s) models. Similarly, using the derived values of the thermalization efficiency (α) and mass-loading parameter (β) in the core (using either the 68% or 95% contours), we estimate an advection timescale of the order of ≃0.5–1 Myr on ≃0.2–0.5 kpc scales along the wind axis.

Zoom In Zoom Out Reset image size

As shown by Thompson et al. (2016), when the cooling time is less than a few times longer than the advection time, strong bulk radiative cooling may result. Thus, on the basis of the comparison in Figure 8, it is possible that the hot outflow from NGC 253 goes through bulk radiative cooling, perhaps changing the character of both the X-ray emission and the Hα emissions along the wind axis. Bulk cooling from the hot phase is important because it provides a physical origin for high-velocity cool-phase outflows (Wang 1995; Silich et al. 2004).

We calculate an estimate of the outflow mass as M = n_e m_H V/1.2f^1/2 where the factor of 1.2 is from our assumed relation of n_e = 1.2n_H, f is the filling factor, which we assume to be 1, and V is the cylindrical volume using a 95% emission radius. We also calculate a mass outflow rate as $\dot{M}={n}_{e}{m}_{{\rm{H}}}{rv}/1.2{f}^{1/2}$ where r is the vertical distance traveled and v is the outflow velocity. For r we use 0.43 kpc in the northern and 0.95 kpc in the southern outflow, where we ignore the 200 pc central starburst region in order to solely capture the outflowing gas. For v we assume a velocity of 1000 km s⁻¹ based on the models described in Section 4.4. We find masses of 1.2 × 10⁶ M_⊙ in the northern outflow and 3.0 × 10⁶ M_⊙ in the southern outflow. For the mass outflow rates with the assumptions from above, we find a rate of 2.8 M_⊙ yr⁻¹ in the northern outflow and 3.2 M_⊙ yr⁻¹ in the southern outflow.

We compare our estimates of mass outflow rate with those from Strickland et al. (2000) and find that they are similar. Strickland et al. (2000) estimated an upper limit on the mass outflow rate of 2.2v₃ M_⊙ yr⁻¹ where v₃ is the outflow velocity in units of 3000 km s⁻¹. Although the mass outflow rates are similar, the assumptions used in the estimates are different and warrant a discussion.

Strickland et al. 2000 used a velocity of 3000 km s⁻¹ based on standard mass and energy injection rates from Leitherer & Heckman (1995). However, we find in our models described in Section 4.4 that the velocity profile produces much lower values. Based on our derived CC85 parameters of α₉₅ ≃ 0.21 and β₉₅ ≃ 0.36, we find that the velocity of the wind is closer to 1000 km s⁻¹. Strickland et al. (2000) also used a filling factor of 0.4. They asserted than on scales of less than about a few hundred parsecs, the pressures of the X-ray- and Hα-emitting phases are similar. By equating their derived pressure (which they write in terms of the filling factor) to the Hα pressure from McCarthy et al. (1987), they were able to derive a filling factor of 0.4.

If we use the assumptions of Strickland et al. (2000) with our formulae, we find that the derived mass outflow rates are too large when compared to the cooler phases. Assuming f = 0.4 and v = 3000 km s⁻¹, we find a mass outflow rate in the hot phase of 13.5 M_⊙ yr⁻¹ in the northern outflow and 15.3 M_⊙ yr⁻¹ in the southern outflow. The estimate for the mass outflow rate is 14–39 M_⊙ yr⁻¹ in the cool molecular phase (Krieger et al. 2019), while it is ≈4 M_⊙ yr⁻¹ in the warm ionized phase (Westmoquette et al. 2011; Krieger et al. 2019). The assumptions of Strickland et al. (2000) would result in an outflow rate in the hot phase on par with that on the molecular phase and greater than that in the warm phase. The results from the assumptions of Strickland et al. (2000) lead to mass outflow rates dominated by the hot phase.

In order to acquire the most accurate measurements of outflow properties, the spectral model needs to include all the components that contribute significantly to the X-ray emission. Previous work has shown that charge exchange, the stripping of an electron from a neutral atom by an ion, is one of these components. In particular, Lopez et al. (2020) found that when including CX in spectral models of the outflows of M82, CX contributes a significant portion (up to 25%) of the broadband (0.5–7 keV) X-ray emission. Similarly, for several nearby star-forming galaxies (including NGC 253), Liu et al. (2012) showed that the Kα triplet of He-like ions required a CX component to account for the observed line ratios. In the context of galactic winds, CX represents the emission of the hot phase interacting with the cooler phases where the hot ionized gas strips electrons from the cooler neutral gas.

In our analysis, we find that three regions' spectral fits were statistically significantly better, based on F-tests, with the inclusion of a CX component: the central and the two southern regions. At these locations, CX accounts for a significant fraction of the emission. As an example, in Figure 9 we show the spectra for Observation 20343 in the central region and note that CX makes up a significant fraction of the spectra. In the central region, the CX contribution is 42%. For region S1, it is 20%, and for region S2, it is 33% of the total broadband X-ray emission. We note that the lack of a CX contribution in the northern outflow may be because of the high column density there, obscuring the soft X-rays necessary to distinguish the CX emission.

Zoom In Zoom Out Reset image size