X-ray Properties of NGC 253's Starburst-Driven Outflow

We analyze image and spectral data from $\approx$365~ks of observations from the {\it Chandra} X-ray Observatory of the nearby, edge-on starburst galaxy NGC 253 to constrain properties of the hot phase of the outflow. We focus our analysis on the $-$1.1 to $+$0.63 kpc region of the outflow and define several regions for spectral extraction where we determine best-fit temperatures and metal abundances. We find that the temperatures and electron densities peak in the central $\sim$250 pc region of the outflow and decrease with distance. These temperature and density profiles are in disagreement with an adiabatic spherically expanding starburst wind model and suggest the presence of additional physics such as mass loading and non-spherical outflow geometry. Our derived temperatures and densities yield few-Myr cooling times in the nuclear region, which may imply that the hot gas can undergo bulk radiative cooling as it escapes along the minor axis. Our metal abundances of O, Ne, Mg, Si, S, and Fe all peak in the central region and decrease with distance along the outflow, with the exception of Ne which maintains a flat distribution. The metal abundances indicate significant dilution outside of the starburst region. We also find estimates on the mass outflow rates which are $2.8\:M_{\odot}/\rm{yr}$ in the northern outflow and $3.2\:M_{\odot}/\rm{yr}$ in the southern outflow. Additionally, we detect emission from charge exchange and find it has a significant contribution ($20-42$\%) to the total broad-band ($0.5-7$ keV) X-ray emission in the central and southern regions of the outflow.

1. INTRODUCTION Starburst galaxies are characterized by their prominent multiphase outflows known as galactic winds (Veilleux et al. 2005;Rubin et al. 2014;Veilleux et al. 2020) and are the result of an intense period of star formation. These galactic winds have important effects on their host galaxies where they alter the metal content of the disk (Peeples & Shankar 2011;Finlator & Davé 2008) and are able to enrich the surround-Previous work on the hot gas component has found temperature, column densities, and metallicity constraints of NGC 253's outflow. Strickland et al. (2000) measured temperatures and column densities in five different regions along the outflow by fitting a single temperature plasma model using data from Chandra. Bauer et al. (2007) measured temperatures along the outflow using emission line ratios from XMM-Newton spectra. Mitsuishi et al. (2013) measured both metallicities and temperatures of the disk, superwind, and halo region using XMM-Newton and Suzaku data. Missing from previous work is an X-ray study of the outflow that incorporates both the coverage and resolution along the outflow of Strickland et al. (2000) and Bauer et al. (2007), with the more sophisticated model and metallicity constraints of Mitsuishi et al. (2013).
Also missing from previous X-ray analyses of NGC 253's outflow is the inclusion of charge exchange (CX) in spectral models used to constrain the temperatures and metal abundances. CX is the stripping of an electron from a neutral atom by an ion. CX has been found to produce X-ray emission in galactic winds Liu et al. 2012) and recently, Lopez et al. (2020) found that charge exchange contributes between 8%-25% of the total absorption-corrected, broadband (0.5-7 keV) X-ray flux in the starburst M82. As a result, the absence of CX in spectral models can lead to inaccuracies in constraints on the outflow properties.
In this paper we follow the methodology employed in Lopez et al. (2020) to constrain temperature, metal abundances, column density, and number density of NGC 253's outflow using Chandra images and spectra. We employ recent spectral models that include CX to account for its contribution to the line emission. In Section 2, we describe the Chandra observations used and provide details for the procedures used to produce Xray images and spectra of NGC 253. In Section 3, we present the results of our spectral fitting along the outflow as well as across the disk. In Section 4, we compare our work to previous X-ray studies of NGC 253, to the results of Lopez et al. (2020) on M82, and to the predictions of galactic wind models. In this section, we also derive mass outflow rate estimates and quantify the contribution of CX to the broadband emission. In Section 5, we summarize our findings and outline future work to be done to constrain the hot phase in outflows.
Throughout the paper we assume a distance of 3.5 Mpc to NGC 253 (Rekola et al. 2005), where 1 =1.0 kpc.
2. METHODS NGC 253 was observed seven times with Chandra from 1999 to 2018, as detailed in Table 1. The first three observations used the ACIS-S array, while the rest used ACIS-I.
Data were downloaded from the archive and reduced using the Chandra Interactive Analysis of Observations  ciao version 4.14 (Fruscione et al. 2006). Using the merge obs function, the observations were combined, and the wavdetect function identified point sources in the image that were removed with the dmfilth command. The detection sensitivity for the removed point sources was 7.14 × 10 −16 ergs/cm 2 /s. The sensitivity was calculated using WebPIMMS 1 , assuming a point source with a power-law spectrum of photon index 1.7. Each of the Figure 2. Three color image of NGC 253, where blue is broad-band (0.5-7 keV) Chandra X-ray, green is Hα (Lehnert & Heckman 1995), and red is CO (2 − 1) (Leroy et al. 2021) emission. The 1 label corresponds to 1.0 kpc, and in the image, North is up and East is left.
observations were also checked for background flares and none were found. The final, broad-band (0.5 − 7.0 keV) X-ray image of the diffuse emission is shown in Figure 1. From this image, the outflow extends ≈0.63 to the north and ≈1.1 south of the starburst. Additionally, the galactic disk is evident in the diffuse X-rays as well, following the spiral structure as seen in the three-color image with CO (Leroy et al. 2021) and Hα (Lehnert & Heckman 1995) shown in Figure 2.

Outflow Spectral Analysis
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). In Figure 3, the areas where spectra were extracted are shown. From the center to the north of the outflow, the regions were 0.25 in height and 1 in width, while the southern regions were 0.5 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).
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 2 . 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 ≥2.5 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 powerlaw 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 MW H = 1.4 × 10 20 cm −2 (Dickey & Lockman 1990) and was frozen. The second component N NGC 253 H 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 non-thermal 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 center region required a second vapec component that accounted for a hotter thermal plasma component. The center, S1, and S2 regions also included an AtomDB 2 charge-exchange (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.

Disk Spectral Analysis
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 the abundances. Instead, the total metallicity of the apec component, Z, was set to so-  lar 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).

Outflow Spectral Results
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 center region spectrum. 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 NGC253 H , 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 χ 2 values in each region being near one, and the largest being the center region with a reduced χ 2 of 1.28 and a null hypothesis of 7.67 × 10 −13 .
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 NGC253 H the gradient is similar in the northern and southern outflow, with a peak of N NGC253 H = (8.6±0.9)× 10 21 cm −2 in the center and decreasing to N NGC253 H ≈ (1 − 2) × 10 21 cm −2 . For kT , it peaks at the center, with kT = 0.98 ± 0.02 keV and decreases going 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 +1.2 −0.8 keV. The metal abundances also peak in the center, with S having the highest value of 6.6 +3.2 −1.8 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 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 −14 EM )/4πD 2 , where the emission measure is

272/258
a Abundances with values of 1 are frozen to solar values in the fits.
b The temperature for the second hotter thermal plasma components is kT2 = 5.5 +1.2 −0.8 keV. H peak in the center and decrease with distance. With the exception of Ne, the metal abundances also peak in the center and decrease with distance. The metal distributions indicate much dilution of the gas outside of the starburst region.  EM = n e n H dV . By setting n e = 1.2n H , we calculate n e as n e = (1.5 × 10 15 normD 2 /f V ) 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 broad-band, 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 define regions larger than the ones used in Section 2.1, the 95% emission radii were overestimated due to more background emission being included. We assume a filling factor of f = 1, though we note that due to processes like mass loading, the filling factor 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 center region with a value of n e = 0.30 cm −3 and decreases along the outflow. The pressure also peaks in the center region with a value of 6.7 × 10 6 K cm −3 . The cooling time has the lowest value of t cool = 17 Myr in the central region, and the largest 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 larger region cooling times appear to overestimate those of the smaller regions.
We also calculate the same physical quantities shown in Table 3 for the hotter temperature component, kT 2 , that is only present in the center region. We find that the hotter component has a lower n e of 0.19 cm −3 than the cooler component. The hotter component has a higher pressure of 2.4 × 10 7 K cm −3 and a longer cooling time of 211 Myr than the respective values in the cooler component.

Disk Spectral Results
Using the spectral model described in Section 2 we measure kT and N NGC253 H for the diffuse X-ray emission in the disk (Table 4). We find that the kT values are in the range of 0.18 − 0.30 keV, and the N NGC253 H values are in the range of (0.31 -0.72)×10 22 cm −2 . As shown in the right panel of Figure 5, we find the highest kT in Region 1 and the lowest in Region 4. Previous X-ray studies of NGC 253 have been performed and are complementary to our analysis. In Strickland et al. (2000), 15 ks of Chandra X-ray observations were used 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 as the ones in this paper. There were five regions each of which measured about 0.33 in height and 0.8 in width. The spectra were fitted using a singletemperature plasma model, where the kT , N NGC253 H , and 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 NGC253 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 NGC253 H (the northern outflow), where softer X-rays are attenuated. By contrast, in the southern outflow where N NGC253 H 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 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 NGC253 H , and 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, including 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 NGC253 H , 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 − 0.61 keV, in the center they range from 0.22 − 0.79 keV, and in the southern regions they range from 0.21 − 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 a factor of about 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 Reflection Grating Spectrometer (RGS). This instrument only covers the 0.33 − 2.5 keV energy range 3 in contrast to our Chandra data that is in the range of 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(Plucinsky et al. , 2018. About 60% of our data was 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 Xray 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 NGC253 H measurement is lower than ours with an upper limit of 0.05 × 10 22 cm −2 , while our value is (0.15 − 0.18) × 10 22 cm −2 .

Comparison to M82 metal abundances and their
distributions 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 if 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 and decrease with distance. Also, for both galaxies, the O abundance is not well constrained in regions of large 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 center 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. (2019Das et al. ( , 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 non-solar abundance ratio was unclear, but it was speculated that the enhancement may be a product of galactic outflows.

Connection to the CGM
In our analysis of the galactic wind of NGC 253, we find α-enhancement already in the outflow regions. We also find super solar 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.

Comparison to wind models
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 NGC253 H (compared to the signal necessary to measure the metal abundances). Specifically, we create 42 regions along the outflow for spectral extraction that are 2.5 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 a 68% and 95% emission radius (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 Chevalier & Clegg (1985) (CC85) galactic wind model. The CC85 model describes a SNpowered galactic wind that assumes spherical outflow geometry and adiabatic expansion. This model provides a good starting point for comparison because of its simplicity.
In Figure 6, we plot the observed kT and n e values (black and purple dashed lines with 'o' markers) with that of the spherical adiabatic CC85 model (gray lines). The disagreement between the CC85 model with that of 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 that of the observed values.
This may be a result of missing processes like 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 The gray dashed line also shows the CC85 kT predictions but for a 68% emission radius. The red line shows kT predictions from a non-spherical wind model using the 95% emission radius. The blue line is the same as the red but for 68% emission radius. Right: Same models as the right except now ne values are plotted for a 95% (black dashed line with 'o' markers) and for a 68% (purple dashed line with 'o' markers) emission radius. We find better agreement with the non-spherical models and find the CC85 model does not characterize the kT and ne profiles well. (Levy et al. 2022). Nguyen & Thompson (2022) recently showed that energy and mass injection in a ring-like geometry naturally produces a bi-conical outflow, where the hot phase is cylindrically collimated. When we include the effect of non-spherical flow geometry into the wind model (red and blue lines in Figure 6), we find better agreement with the profiles.
The enclosed 68%-ile and 95%-ile volumes, which were used to determine n e , can also be used to define an outflow geometry A(s) = πr(s) 2 , where r(s) is the radius of the 68%-ile or 95%-ile volume at height s perpendicular to the wind axis. In Figure 7, we plot A(s) and d 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 non-spherical. 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.
In Figure 8 we show the calculated cooling times (t cool ) for a 68% and 95% emission radius, along with the advection times (t adv ) from both the spherical (orange lines) and non-spherical (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 order (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 order 0.5 − 1 Myr on 0.2 − 0.5 kpc scales along the wind axis.
As shown by Thompson et al. (2016), when the cooling time is less than a few times longer then 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).

Mass Outflow Rates
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 1, and V is the cylindrical volume using a 95% emission radius. We also calculate a mass outflow rate asṀ = n e m H rv/1.2f 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 outflow- Figure 7. Left: The non-spherical area function for NGC 253's hot wind derived from the observed X-ray (0.5-7 keV) surface brightness contours perpendicular to the minor axis as a function of height s above and below the starburst disk at s = 0. The black circles and diamonds show the lateral size of the contour enclosing 68 % and 95% of the emission, respectively, as a function of distance along the minor axis. The blue and red lines are a cubic spline, fit to both the southern and northern sides independently, for the 68% and 95% contours, respectively. Right: The areal expansion rate for the cubic spline fits to A(s) and for spherical expansion (gray line). These profiles demonstrate that the wind does not have the geometry assumed by the CC85 model. ing gas. For v we assume a velocity of 1000 km/s based on the models described in Section 4.4. We find a mass in the northern outflow of 1.2 × 10 6 M and a mass of 3.0 × 10 6 M in the southern outflow. For the mass outflow rates with the assumptions from above, we find a rate of 2.8M /yr in the northern outflow and 3.2M /yr in the southern outflow. We compare our mass outflow rate estimates with those from Strickland et al. (2000) and find that they are similar. Strickland et al. (2000) estimated a mass outflow rate upper limit of 2.2v 3 M /yr where v 3 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 α 95 0.21 and β 95 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 a few hundred parsecs, the pressures of the X-ray and Hα emitting phase 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 Strickland et al. (2000) assumptions 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 −1 , we find a hot phase mass outflow rate of 13.5 M /yr in the northern outflow and 15.3 M /yr in the southern outflow. In the cool molecular phase, the estimate for the mass outflow rate is 14 − 39 M /yr (Krieger et al. 2019), while in the warm ionized phase, the mass outflow rate is ≈ 4 M /yr (Westmoquette et al. 2011;Krieger et al. 2019). The Strickland et al. (2000) assumptions would result in a hot phase outflow rate on par with the molecular phase rate and greater than the warm phase. The results from the Strickland et al. (2000) assumptions lead to mass outflow rates dominated by the hot phase.

Charge Exchange contribution
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 (CX), 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 NCC 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 Ftests, with the inclusion of a CX component: the center 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 center region, the CX contribution is 42%. For region S1, it is 20%, and for region S2, it is 33% of the total broad-band 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.

CONCLUSIONS
In this paper, we analyze ≈365 ks of archival Chandra data from NGC 253 to produce constraints on its starburst driven outflow that spans 1.1 kpc south and to 0.6 kpc north from the galaxy's center. We extract spectra and use models to find the best-fit parameters (kT , N NGC253 H , and metal abundances) in five regions along the outflow (Figure 3) and in nine regions of the galaxy's disk ( Figure 5).
With the exception of Ne which has a nearly flat distribution, the metal abundances peak in the galactic disk and decrease sharply along the minor axis ( Figure 4 and Table 2). The distributions indicate significant dilution outside of the starburst region. For kT and N NGC253 H , the values also peak in the galactic disk and decrease along the minor axis. The decrease in abundances along the outflow may suggest mass loading or indicate other physics should be considered, such as non-equilibrium ionization.
We constrain the extent of the outflow in the regions and calculate n e , pressure, and cooling times (Table 3). With these measurements we provide estimates of the hot phase mass outflow rates in Section 4.5, finding them to be 2.8 M /yr in the northern outflow and 3.2 M /yr in the southern outflow. We also find that the measured cooling times are of order the advection times predicted from models potentially indicating the presence of bulk radiative cooling in the outflow in Section 4.4. The models also show that incorporating non-spherical geometry provides a better match to observed kT and n e profiles. We also find that the central region and southern outflow have significant contribution of 20 − 42% to the broad-band X-ray emission from a CX component in Section 4.6.
In the future, expanding this approach to other starburst galaxies will reveal whether outflow properties vary with host galaxy properties, such as stellar mass. Chisholm et al. (2018) already showed that there is an inverse relationship between the metal loading factor and galaxy stellar mass in warm gas outflows. It would be worth investigating whether this relationship also holds in the hot phase, where most of the metals are being driven out.
Additionally, understanding of the hot gas in starburst-driven outflows will advance tremendously with the launch of XRISM (Tashiro et al. 2018). The mission will have higher spectral resolution, facilitating measurement of the hot wind velocity and enabling better constraints on the wind energy (Hodges-Kluck et al. 2019) and improving our hot phase mass outflow rate estimates.

ACKNOWLEDGMENTS
We thank Ann Hornschemeier for providing us with the NGC 253 Chandra data and allowing us to use it for this work. We also thank Adam Leroy and the OSU Galaxy/ISM Meeting for useful discussions. SL and LAL were supported by NASA's Astrophysics Data Analysis Program under grant number 80NSSC22K0496. LAL also acknowledges support by the Simons Foundation, the Heising-Simons Foundation, and a Cottrell Scholar Award from the Research Corporation for Science Advancement. ADB acknowledges support from NSF-AST2108140. This work was performed in part at the Simons Foundation Flatiron Institute's Center for Computational Astrophysics during LAL's tenure as an IDEA Scholar. SM gratefully acknowledges the support provided by the National Aeronautics and Space Administration through Chandra Award Number AR2-23014X issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. DDN and TAT are supported by NSF grant 1516967, NASA ATP 80NSSC18K0526, and NASA 21-ASTRO21-0174. AS is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1903834. Software: CIAO (v4.14;Fruscione et al. 2006), XSPEC (v12.12.1;Arnaud 1996)