Supernova feedback efficiency and mass loading in M82

Supernova feedback efficiency and mass loading in the starburst and galactic superwind exemplar M82

David K. Strickland11affiliation: Department of Physics and Astronomy, The Johns Hopkins University, 3400 N. Charles St., Baltimore, MD 21218, USA. and Timothy M. Heckman11affiliation: Department of Physics and Astronomy, The Johns Hopkins University, 3400 N. Charles St., Baltimore, MD 21218, USA.

We measure the net energy efficiency of supernova and stellar wind feedback in the starburst galaxy M82 and the degree of mass-loading of the hot gas piston driving its superwind by comparing a large suite of 1 and 2-dimensional hydrodynamical models to a set of observational constraints derived from hard X-ray observations of the starburst region (the fluxes of the He and Ly-like lines of S, Ar, Ca and Fe, along with the total diffuse – 8 keV X-ray luminosity). These are the first direct measurements of the feedback efficiency and mass-loading of supernova heated and enriched plasma in a starburst galaxy. We consider a broad range of plausible parameters for the M82 starburst, varying the age and mode of star formation, the starburst region size and geometry, and supernova metal yields. Over all these varied input parameters all the models that satisfy the existing observational constraints have medium to high thermalization efficiencies (% %) and the volume-filling wind fluid that flows out of the starburst region is only mildly centrally mass loaded (). These results imply a temperature of the plasma within the starburst region in the range 30 – 80 million Kelvin, a mass flow rate of the wind fluid out of the starburst region of and a terminal velocity of the wind in the range . This velocity is considerably larger than the escape velocity from M82 () and the velocity of the H emitting clumps and filaments within M82’s wind (). Drawing on these results we provide a prescription for implementing starburst-driven superwinds in cosmological models of galaxy formation and evolution that more accurately represents the energetics of the hot metal-enriched phases than the existing recipes do.

Subject headings:
ISM: bubbles — galaxies: individual : NGC 3034 (M82) — galaxies: halos — (galaxies:) intergalactic medium — galaxies: starburst — X-rays: galaxies

1. Introduction

The presence of metals in the intergalactic medium (IGM) at all red-shifts (Songaila, 1997), the present-day galaxy mass-metallicity relationship (Garnett, 2002; Tremonti et al., 2004), and -kpc-scale holes in the IGM at red-shift (Adelberger et al., 2003) have all been interpreted as the consequence of energetic metal-bearing outflows from galaxies.

There is no doubt that at least one form of galactic wind exists: starburst galaxies are commonly observed to host supernova-driven galactic-scale winds (Heckman et al., 1990; Heckman, 2003) in both the local and high redshift Universe (Lehnert & Heckman, 1996a; Pettini et al., 2001; Shapley et al., 2003). These superwinds are phenomenologically complex and are comprised of multiple different gas phases moving at different velocities (Bland-Hawthorn, 1995; Dahlem, 1997; Heckman et al., 2000, 2001; Veilleux et al., 2005). However, the very hot plasma long hypothesized to be the prime mover in superwinds (Chevalier & Clegg, 1985), often called the wind fluid, has only recently been directly observed in the archetypal superwind associated with the starburst galaxy M82 (Griffiths et al., 2000; Strickland & Heckman, 2007).

The most fundamental long-standing uncertainty in our understanding of the nature of superwinds and their larger scale role is the value of the energy per particle of the metal-enriched wind-fluid that drives an individual superwind. Both supernovae (SNe) and stellar winds from massive stars are responsible for supplying the kinetic energy that drive superwinds (Stellar winds only contribute % of the total mechanical power of a starburst, but they may supply up to half of the gas mass returned to the ISM, see Leitherer et al. 1992; Leitherer & Heckman 1995). What fraction of the of kinetic energy released per star that ultimately goes supernova is actually available to drive bulk motions of the ISM? As it is understood that SNe dominate the energy return this is often termed the SN thermalization efficiency or the SN feedback efficiency, even though stellar winds do play some small part energetically. Is the metal-enriched mix of SN and stellar wind ejecta diluted and cooled by mixing with the ambient ISM, i.e. is the wind-fluid mass-loaded? If this ejecta is to escape the galaxy and enter the IGM then its specific energy must exceed its initial gravitational potential energy, and furthermore the total energy in the wind-fluid must be sufficient to perform the mechanical work necessary to move the interstellar and halo gas that lies between the starburst region and the IGM (Mac Low & McCray, 1988; Strickland et al., 2004b).

For moderate to high thermalization efficiency (efficiencies 10%), the wind-fluid is expected to have an initial temperature within the starburst region in the range 10 – 100 million Kelvin (Chevalier & Clegg, 1985), even if it has been lightly mass-loaded with cold ambient gas. Such temperatures are far hotter than the neutral, warm ionized or soft X-ray emitting plasmas that are the most easily observed phases within superwinds (Dahlem, 1997; Veilleux et al., 2005).

However it has only recently become technically possible to detect the diffuse hard X-ray emission that would be associated with a plasma in the temperature range (K) . The starburst regions of even the nearest starburst galaxies only subtend on the sky and are heavily populated with X-ray-luminous compact objects. It was only with the advent of -spatial resolution X-ray astronomy in 1999 following the launch of the Chandra X-ray Observatory that made it possible for Griffiths et al. (2000) to resolve out the point source emission in the archetypal nearby starburst galaxy M82 and reveal residual diffuse hard X-ray emission, although of uncertain origin.

It is important to distinguish between the material that actually drives a superwind, and the other more easily observed phases of swept-up and entrained gas that are part of the wind. We define the wind-fluid as the material whose high thermal and/or ram pressure drives the superwind. The wind fluid certainly comprises a significant fraction of the merged SN ejecta and massive star stellar wind material, probably admixed (mass-loaded111We will use the term mass loading to denote the addition and efficient mixing of ambient ISM into the full volume of the wind-fluid, which is consistent with the terminology in Suchkov et al. (1996). Technically mass loading need not affect the entire flow and can be localized, see e.g. Dyson et al. (1993). Furthermore the form of mass-loading considered in this work occurs only within the starburst region itself (“centralized mass loading”) although ultimately it affects the entire flow. We find it conceptually beneficial to differentiate between mass loading and entrainment in superwinds. Mass loading affects the density of one gaseous phase by mixing in material from another phase, almost always a cooler and denser phase. The origin and fate of the secondary phase is unimportant given the focus on the properties of the primary phase. Entrainment is the incorporation into the superwind of subsidiary gas phases that leads to creating or sustaining a multiphase superwind. The entrained material is usually ambient disk or halo gas over-run or swept-up by the wind fluid, and retains much of its original character, albeit with some modification. Although entrainment processes may lead to mass loading of the wind fluid the two processes are not synonymous. In principle entrainment could occur without leading to significant mass loading of the wind fluid. One can not realistically hope to represent or approximate the totality of a multiphase superwind as a single phase mass loaded flow. Only specific components of the superwind are amenable to such treatment, of which the wind fluid is the clearest example.) with some ambient gas from with starburst region.

Although other gas phases may dominate the total gas mass budget of a superwind, they are less important with respect to the issue of energy and metal transfer into the IGM. Furthermore the ejection of these phases is less than certain, as the measured warm neutral and warm ionized medium (WNM, WIM) gas velocities are typically comparable to the galactic escape velocities (Heckman et al., 2000; Martin, 2005; Rupke et al., 2005). Unfortunately even the most recent attempts to include superwinds in cosmological models of galaxy formation use only simplified single-phase winds, often with energetics and masses based on the observed WNM/WIM properties of superwinds (see e.g.  Springel & Hernquist, 2003a; Bertone et al., 2005). These models do not properly address the physics of hot-phase metal transport in addition to having overly low values of the energy per particle, which casts the reliability of their conclusions into doubt.

Attempts have been made to use measurements of the temperature of the extended soft X-ray emitting plasmas ( K) in starburst galaxies to constrain whether the hot gas will escape the host galaxy (see e.g.  Martin, 1999), but such an approach neglects both the significant kinetic energy expected to be associated with this material and that the soft X-ray emitting plasma is an indirect probe of the properties of the true wind-fluid (Strickland & Stevens, 2000).

The detection of apparently thermal diffuse hard X-ray emission of temperature K in the starburst region of M82 in the 33 ks-long Chandra ACIS-I observations of Griffiths et al. (2000) suggested it was possible to directly observe and constrain the properties of the wind fluid in a starburst superwind.

We recently obtained and analyzed a new 18 ks Chandra ACIS-S observation of M82, and also reprocessed archival Chandra ACIS-I and XMM-Newton observations222Chandra ObsId numbers ADS/Sa.CXO#Obs/2933 (catalog 2933), ADS/Sa.CXO#Obs/361 (catalog 361) and ADS/Sa.CXO#Obs/1302 (catalog 1302). The XMM-Newton data ObsIds are 0112290201 and 0206080101. (total exposure of 48 ks with ACIS-I, 75 ks with XMM-Newton PN and 102 ks with each MOS instrument). We confirmed the presence of both diffuse continuum and diffuse Fe He line emission within the central starburst region (Strickland & Heckman, 2007, hence forth referred to as Paper I). Having assumed that the line-emitting plasma is metal-enriched to , we found that the observed iron line luminosity of was consistent with the properties of the wind fluid expected given M82’s star formation rate, provided that the SN thermalization efficiency was high. However the diffuse hard X-ray continuum luminosity of was too luminous to be the thermal bremsstrahlung associated with the wind fluid given the parameters assumed when interpreting the iron line emission.

In this paper we consider in greater detail what constraints the observed hard X-ray properties of M82 can place on the properties on the superwind and the wind-fluid, in particular the efficiency of SN heating () and the degree of central mass-loading (). Ultimately these determine the energy per particle of the wind fluid. We restrict analysis to the starburst region within pc of the center of M82. We calculate the X-ray emission using both multi-dimensional hydrodynamical simulations and the 1-dimensional analytical Chevalier & Clegg (1985) model (hence forth referred to as the CC model). By using the simple but relatively accurate CC model we can explore the predicted X-ray properties of the wind fluid in models covering a wider range of parameter space than if we used multi-dimensional hydrodynamical simulations alone.

A variety of recent papers have presented theoretical calculations of the broad-band X-ray emission from individual star-forming clusters, using either the adiabatic CC model, or the non-adiabatic solution of Silich and collaborators (Cantó et al., 2000; Stevens & Hartwell, 2003; Oskinova, 2005; Silich et al., 2005; Ji et al., 2006), with the aim of seeing if the particular model used can explain the observed data. Our motivation differs from these papers, in that our aim is to use theoretical models in combination with observational data on M82 to constrain fundamental properties of the superwind. In addition our approach differs from these previous papers in its use of X-ray line luminosities rather than purely broad-band X-ray luminosities. In reality the broad-band emission might be dominated by non-thermal continuum sources. Furthermore line emission is more sensitive than the thermal continuum luminosity to interesting parameters such as the plasma temperature.

2. Observational constraints on high temperature gas in M82

The observational constraints derived in Paper I on the diffuse hard X-ray emission within pc of the nucleus of M82 from the Chandra and XMM-Newton observations are summarized in Table 1. These comprise a robust 6.69 keV iron line detection ( in significance), an upper limit on the 6.96 keV iron line, and from the Chandra observations the luminosity of the diffuse continuum. For the Chandra-based observations the values refer to the diffuse emission alone, after the exclusion of detected point-like X-ray sources and corrections for undetected point sources based on the observed point source luminosity function. These luminosities also include 32% and 20% upward corrections to the luminosities derived from the ACIS-S and ACIS-I data to account for regions excluded from the spectra due to remove the emission from detected point sources. For the XMM-Newton observation the line luminosities include all emission within pc, but the contribution to the line luminosities from point-like sources is expected to be small in this particular observation.

Instrument ObsID
(1) (2) (3) (4) (5) (6)
XMM-Newton EPIC 0206080101
Chandra ACIS-S 2933
Chandra ACIS-I 361+1302

Note. – This table summarizes the results presented in Strickland & Heckman (2007) on the emission from within pc of the nucleus of M82. Column 2: The ObsID is the identification number allocated to each observation by the Chandra and XMM-Newton Science Centers. Column 3: Logarithm of the total diffuse X-ray luminosity in the – 8 keV energy band, including line emission. Column 4: Logarithm of the keV Fe K line luminosity. Column 5: Logarithm of the keV Fe He line luminosity. Column 6: Logarithm of the keV Fe Ly line luminosity. The quoted uncertainties are 68.3% confidence for one interesting parameter, while the upper limits are 99.0% confidence.

Table 1Diffuse continuum and iron line luminosities from Paper 1.

Note that the region of bright diffuse hard X-ray emission studied in this paper is associated with the nuclear starburst region. As can be seen from Fig. 1 this region is much smaller than the full 12-kpc-scale superwind that is most often studied in optical (Lehnert et al., 1999) or soft X-ray emission (Stevens et al., 2003).

Figure 1.— A three color composite image of M82 showing the region within a 500 pc radius of the nucleus (white circle) in comparison to the galaxy and superwind. Soft X-ray emission in the 0.3–2.8 keV energy band is shown in red, optical R-band emission (starlight) in green, and diffuse hard X-ray emission in the 3 – 7 keV energy band in blue. The X-ray images are adaptively smoothed Chandra ACIS-S images with point source emission removed and interpolated over. The image is ( kpc) on a side.

As discussed in Paper I the nature of the diffuse broad-band continuum is not clear. If the continuum is thermal then the relative strength of the iron line emission to the observed continuum is much lower than expected for a plasma of Solar elemental abundances, let alone the iron-enriched plasma expected from recent core-collapse SNe. We demonstrated that a wind mass-loaded with additional ambient gas could produce the total broad-band X-ray luminosity that we observe, but the required degree of mass loading could not sufficiently reduce the relative strength of the iron line to the continuum. Thus the best fit bremsstrahlung models for the continuum do not necessarily provide meaningful estimates of the temperature of the iron-emitting plasma, and should be interpreted with caution. Power law spectral models provide a marginally better fit to the continuum than do thermal bremsstrahlung models. Thus a non-thermal process such as Inverse Compton X-ray emission might be responsible for the majority of the diffuse hard X-ray continuum. The X-ray luminosity associated with Inverse Compton emission depends sensitively on the unknown magnetic field strength within the starburst region. In Paper I we estimated it might account for % of the observed diffuse X-ray luminosity in the – 8 keV energy band, although other authors have in the past given higher luminosity estimates for Inverse Compton emission from M82 (Schaaf et al., 1989; Moran & Lehnert, 1997).

The upper limits on the iron 6.96 keV / 6.69 keV line flux ratio333Note that the line ratios quoted in Paper I were the ratio of the photon fluxes, whereas in this paper we quote the energy flux ratio. of place upper limits on the temperature of if the plasma is in collisional equilibrium (see Table 2). The Coulomb relaxation and ionization equilibrium timescales presented in Paper I are much shorter than the timescale for plasma to flow out of the starburst region, so the assumption of collisional ionization equilibrium is reasonable.

Basis for temperature estimate Detector (K)
(1) (2) (3)
Fe Ly / Fe He XMM-Newton
Fe Ly / Fe He ACIS-S
Fe Ly / Fe He ACIS-I
Ca He / Fe He XMM-Newton
Ca He / Fe He ACIS-S
Ar He / Fe He XMM-Newton
Ar He / Fe He ACIS-S
S He / Fe He XMM-Newton
S He / Fe He ACIS-S
(Ca Ly + Ca He) / Fe He XMM-Newton
(Ca Ly + Ca He) / Fe He ACIS-S
(Ar Ly + Ar He) / Fe He XMM-Newton
(Ar Ly + Ar He) / Fe He ACIS-S
(S Ly + S He) / Fe He XMM-Newton
(S Ly + S He) / Fe He ACIS-S
S He / Ca He XMM-Newton
S He / Ca He ACIS-S
S He / Ar He XMM-Newton
S He / Ar He ACIS-S
Ar He / Ca He XMM-Newton
Ca Ly / Ca He XMM-Newton
Ca Ly / Ca He ACIS-S
Ar Ly / Ar He XMM-Newton
Ar Ly / Ar He ACIS-S
S Ly / S He XMM-Newton
S Ly / S He ACIS-S
Bremsstrahlung fit to diffuse continuum ACIS-S
Bremsstrahlung fit to diffuse continuum ACIS-I

Note. – Temperature estimates are given only for line ratios where one or more of the lines involved is detected at . The quoted uncertainties are 68.3% confidence for one interesting parameter, while both upper and lower limits are 99.0% confidence.

Table 2Temperature estimates for the diffuse hard X-ray emission.

Other lines in the (keV) hard X-ray band may in principle be used to constrain the plasma temperature in a starburst region. In particular the helium-like and hydrogen-like ions of sulphur (at keV and keV respectively), argon ( keV and keV) and calcium ( keV and keV) can be detected in CCD spectra from the Chandra ACIS and XMM-Newton EPIC instruments. These lines are most intense in plasmas with temperatures in the range – 7.5, somewhat lower than the optimal temperature for the iron lines, but this might be advantageous in diagnosing starbursts where mass loading is significant and/or the efficiency of SN thermalization is low.

Figure 2.— (Upper panel) The XMM-Newton EPIC pn (black data points with error bars), MOS1 (red data points) and MOS2 (green data points) spectra of the central 500 pc of M82. The best-fit spectral model (red solid line) is shown, along with the named X-ray line components of the model (dashed and dotted lines with a color matching the parent data). The He-equivalent lines of S, Ar, Ca and Fe, along with Si Ly and S Ly, are all detected with significances of . The other lines have lower significance. (Lower panel) The residuals from the best-fit model.

In Paper I we specifically excluded data at the expected energies of the S, Ar and Ca lines from the fitting of the M82 nuclear region spectra. Repeating the analysis described Paper I but now fitting all the data in the energy range – 9.0 keV in both the XMM-Newton and Chandra ACIS-S spectra yields the line detection significances and line luminosities given in Table 3. We fit for the intensity of the helium-like and hydrogen-like S, Ar and Ca lines, having fixed the line energies at the expected values. These lines are faint, but in some cases likely to be real. As an example we display the best-fit to the XMM-Newton EPIC pn + MOS spectra in Fig. 2. Detection significance was assessed using Monte Carlo methods with several thousand simulated spectra. Where the line detection significance was we give the 99% upper limit on the line luminosity (and not the best-fit value). The Chandra ACIS-I spectra gave similar results to those from ACIS-S, but as the line detections are less significant than in the XMM-Newton data there is little reason to present the values.

Line Ion Energy Detector Significance
(1) (2) (3) (4) (5) (6)
He S XV 2.45 XMM-Newton
S Ly S XVI 2.62 XMM-Newton
Ar He Ar XVII 3.13 XMM-Newton
Ar Ly Ar XVIII 3.32 XMM-Newton
Ca He Ca XIX 3.90 XMM-Newton
Ca Ly Ca XX 4.11 XMM-Newton
S He S XV 2.45 ACIS-S
S Ly S XVI 2.62 ACIS-S
Ar He Ar XVII 3.13 ACIS-S
Ar Ly Ar XVIII 3.32 ACIS-S
Ca He Ca XIX 3.90 ACIS-S
Ca Ly Ca XX 4.11 ACIS-S

Note. – The luminosities are given as of the line luminosity in units of , assuming a distance to M82 of Mpc and zero intrinsic absorption. Line energies are in units of keV. Line significances are the values equivalent to the enclosed probability for a two-tailed Gaussian, where is the probability of the line being spurious. In the case of highly significant lines we can only estimate a lower limit to the significance of the line. The quoted uncertainties on the luminosities are 68.3% confidence, while upper limits on are provided at 99.0% confidence if the line significance is .

Table 3Hard X-ray S, Ar and Ca line detections and luminosities in M82.

In the XMM-Newton EPIC pn + MOS spectra of M82’s nucleus the equivalent of the He line from highly ionized S, Ar, Ca and Fe are all detected with significances , while in the shorter (and hence lower S/N) Chandra ACIS-S spectrum of the nuclear diffuse emission the S He and Fe He are detected at while the Ar He and Ca He lines are only detected at . In both instruments Si Ly is also detected at high significance, but we shall ignore this lower ionization state line in the analysis to follow. The only Ly-like line detected at is S Ly in the XMM-Newton EPIC spectra.

Using the robust line detections (i.e. those with significance ) we can estimate the plasma temperature from the ratio of S Ly to S He, and from between different elements from the ratio of the various He lines (assuming their relative elemental abundances are Solar). These temperature estimates are presented in Table 2. The S ions tightly constrain the ion temperature to , while the ratio of the S, Ar or Ca He line to the Fe He line yields . The ratio of S He to either Ar He or Ca He yields a temperature estimate of .

This pattern of the inferred temperature increasing as higher ionization state ions are used (which are of course more sensitive to higher temperatures) indicates that either the plasma is multi-phase, and/or the the plasma is not in collisional ionization equilibrium (CIE). Of these explanations we favor the former, as multi-dimensional numerical hydrodynamical simulations predict a broad range of gas temperatures below (see e.g.  Strickland & Stevens, 2000). Furthermore we discussed in Paper I reasons to expect the Fe-emitting plasma to be in or close to CIE.

Unfortunately these results do not strongly constrain the temperature of the Fe-emitting plasma. One possibility is that the iron line emission arises in the very hottest phases in the superwind, i.e. at , and that most of the S, Ar and Ca emission comes from somewhat cooler components associated with the wind/ISM interaction and the larger scale soft X-ray emission in the M82 superwind. In this case, given the possible soft X-ray contribution, the S, Ar and Ca observed luminosities must be treated as upper limits to the S, Ar and Ca line emission coming from the volume filling wind-fluid. Alternatively the wind fluid is relatively cool, with a temperature corresponding to that measured from the Ar and Ca lines, i.e. , and the Fe lines are produced by some unexplained non-thermal mechanism.

In either case a successful model of M82’s nuclear emission must reproduce the observed hard X-ray iron line emission without over-producing the broad-band continuum and S, Ar and Ca line emission.

3. Theoretical and Numerical Methodology

Hard X-ray emission in multi-dimensional hydrodynamical models of superwinds is dominated by the thermal emission from starburst region itself (Tomisaka & Bregman, 1993; Suchkov et al., 1994; Strickland & Stevens, 2000).

If the starburst region is spherical then the plasma properties of this region obtained from hydrodynamical simulations agree very well with the simple one-dimensional radially-symmetric analytical CC model, in which the mechanical energy and returned gaseous ejecta from core-collapse supernovae are injected at a uniform rate only within a radius , and radiative energy losses are assumed to be negligible. Gravity and any surrounding ambient gaseous medium are ignored.

As the feedback parameters we wish to constrain are related to the observable X-ray emission via the plasma properties in the starburst region it is worthwhile for us to explicitly reproduce the analytical solution obtained by Chevalier & Clegg (1985). We consider the conditions under which the adiabatic assumption inherent in the CC model may be invalid in Appendix A, but we will show that for most reasonable parameters associated with the M82 starburst (described below) the radiative energy losses are not significant.

In real galaxies starburst regions are better described as disks or ring of star-formation rather than spherical regions. To explore the predicted X-ray emission for our default model of the M82’s disk-like starburst we begin by using 2-dimensional hydrodynamical simulations in cylindrical symmetry to predict the plasma properties within the a radius of 500 pc of the center of the galaxy.

However, with appropriate scaling (see Appendix B), it is possible to use the 1-dimensional CC model with a spherical starburst to predict the plasma properties of the wind fluid within a non-spherical starburst region to high accuracy, and to predict the hard X-ray emission from even larger regions with an accuracy no worse than the uncertainties in the observed X-ray luminosities.

These scaled CC models are considerably faster to calculate than a full hydrodynamical model, so we will use them to explore how the constraints on the thermalization efficiency and mass loading change if we adopt a different range of parameters for the M82 starburst.

This approach, in particular the use of the scaled CC model, is applicable only to the study of the wind fluid in the vicinity of the starburst region. It is not intended for use as a model of an entire superwind as it neglects the interaction between the wind fluid and the multiphase ambient ISM. The simplification of neglecting the ambient ISM works because the only significant source of thermal hard X-ray emission in a superwind is the wind fluid itself, in particular the wind fluid within the starburst region. The interactions between the wind fluid and the multiphase ISM are important in the large scale wind and for its emission at soft X-ray and longer wavelengths (e.g.  Strickland & Stevens, 2000; Cooper et al., 2008), but with the exception of altering the net efficiency of SN+SW heating and any mass-loading of the wind fluid these interactions are not expected to alter the properties of the wind fluid within the starburst region to any significant degree given its energetic dominance.

3.1. 2-Dimensional Hydrodynamical Modeling

The hydrodynamical simulations were performed using the multi-dimensional numerical hydrodynamics code VH-1 (Blondin et al., 1990). Simulations were performed in 2-dimensional cylindrical symmetry with a cell size of 2 pc over 500 by 500 cell computational grid. Reflecting boundary conditions were imposed along the and axes, with inflow/outflow boundary conditions are the other edges of the computational grid.

Energy and mass are injected within the starburst region in the manner described in Strickland & Stevens (2000). In general the starburst region was modeled as disk-like region of the size described in § 3.3. No gravitational forces or radiative cooling are considered for the purposes of consistency with the analytical CC model described below.

The computational grid was initially populated with a tenuous ISM of density and temperature K. The simulations were run until the wind had blown the initial ISM off the grid and reached a steady-state solution, which typically occurred – 2 Myr after the onset of the simulation. Gas properties and X-ray luminosities are measured at a time of 3 Myr within a spherical region of radius of 500 pc to closely match the observed region discussed in § 2.

3.2. The 1-Dimensional Analytical CC Solution

For completeness we present the analytical flow solution derived by CC, and how we use it to derive the fluid properties as a function of radial position. More complicated, but also more generalized, methods of deriving the solution are presented in Cantó et al. (2000) and Silich et al. (2004). The basic model is completely specified by only three parameters: the energy and mass injection rates (); and the radius of the starburst region () within which the energy and mass are injected.

CC present the stationary solution to the hydrodynamical equations as a function of the Mach number , where is the local speed of the flow and the sound speed is . The pressure , density and adiabatic index have their normal meanings. At radii mass and energy are injected at a uniform constant rate per unit volume and , where the starburst volume , and the solution for the Mach number is


This approach to the energy and mass injection implicitly assumes that the exact manner in which stars return mechanical energy to the ISM can be ignored due to rapid local thermalization and mixing, i.e. individual stellar winds and SN blast waves collide and merge on spatial scales small compared to the size of the starburst region. If radiative energy losses are to be small then these collisions must occur with a time shorter than the cooling time scale for stellar wind bubbles or supernova remnants (SNRs). Theoretical investigations of some of these issues may be found in Larson (1974), Heckman et al. (1990), Cantó et al. (2000), Melioli & de Gouveia Dal Pino (2004) and Melioli et al. (2005).

At larger radii then , and the solution is


For any given fractional radius the Mach number is found using a numerical root-finding technique (in our case the Newton-Raphson method) on the appropriate equation.

To obtain the density at any radius we integrate Equation 1 of CC to obtain


Within the starburst region the gas density, pressure and temperature vary only slowly with radius, with maxima at . The central density is effectively determined by the rate at which the injected material can flow out of the starburst region, where is the terminal velocity of the wind at large radius, and the central temperature is set by the energy injected per particle. More specifically, the central temperature and density are and for .

At radii the mass flow rate through the spherical surface of radius is equal to the total mass injection rate within the starburst region, so the density is


To derive the local velocity of the flow we must first determine the speed of sound. If we make the assumption that the flow is adiabatic then we can obtain the speed of sound at any point ( or ), as we know the Mach number and the initial energy per particle is conserved. Specifically we solve Equation 3 of CC and substitute for using Equation 3, to obtain


The speed of the flow at any point is then simply . For the flow is subsonic, with . The sonic point is reached at , beyond which the flow velocity rapidly tends toward the terminal velocity .

With the density and the speed of sound determined at every point within the flow the thermal pressure is also now known. The kinetic temperature , where is the mean mass per particle for a fully ionized plasma of Solar composition.

Numerically we derive the gas density, temperature and velocity as a function of radius, iterating outward in thin radial shells. The numerical implementation we use is based on code originally created by Ian Stevens (Stevens & Hartwell, 2003). This semi-analytical code reproduces the density and pressure to within 0.1% of the values predicted by the full 2-D hydrodynamical simulation using the VH-1 code for a spherical starburst region. The scaling used to approximate a non-spherical starburst region using the spherical CC model is described in Appendix B.

Figure 3 presents the radial variation in several fluid variables and associated derived quantities for one particular set of input parameters, as predicted using our 1-dimensional CC model code.

Figure 3.— (a) An example of the radial CC solution described in § 3.2. Note the relatively uniform number density, thermal pressure and temperature within the assumed starburst region of radius pc, and rapid drop in those parameters outside the starburst region. The net energy and mass injection rates assumed are and . (b) Local wind velocity and sound speed . (c) Logarithms of the ionization-state-related parameter (Equation 14) and cumulative volume emission integral . The plasma is in ionization equilibrium when . Values are given in c.g.s units unless it is stated otherwise.

3.3. Parameters of the M82 nuclear starburst

3.3.1 The starburst age and magnitude

The flow times for the region of interest ( pc) are short (a few Myr) compared to the lifetime of the massive stars powering the flow ( Myr). The properties of the starburst region will very rapidly approach those of the steady-state CC solution after any change in the energy and mass injection rate.

Nevertheless, in order to apply the CC model to M82 we must determine the history (i.e. age and star formation rate as a function of time) of the starburst event in order to derive the current mechanical energy and mass injection rates and the elemental abundances in the injected gas. Furthermore, determining the degree to which the M82 starburst event can be treated as one the simple limiting cases of either a single instantaneous burst (SIB) or as continuous star formation (CSF) is important, in particular for comparison to Starburst99 models (Leitherer et al., 1999).

Recent observational studies have greatly expanded our knowledge of M82’s global star formation history since its close encounter with M81 approximately 0.6 — 1 Gyr ago (see e.g.  de Grijs, 2001). Nevertheless significant gaps still remain in our understanding of the history of the nuclear starburst.

Exterior to the central kiloparsec star formation appears to have peaked Myr ago in a burst of comparable intensity to the current nuclear starburst event, but star formation effectively ceased in the disk – 30 Myr ago (de Grijs et al., 2001; Mayya et al., 2006).

Within the central kiloparsec (projected separations from the nucleus of pc) the majority of recent studies have concentrated on the 200 compact star clusters found there (O’Connell et al., 1995; Melo et al., 2005). The brightest star clusters typically appear to have ages Myr (Satyapal et al., 1997; McCrady et al., 2003), although two massive clusters (M82-F and M82-L) at the edge of this region have ages of Myr (Smith & Gallagher, 2001; Smith et al., 2006a). Förster Schreiber et al. (2003a, b) present a detailed analysis of the NIR emission (1 — 45) from a variety of regions within pc, and find that two bursts (of ages and Myr respectively) are required in each region to simultaneously explain the strength of both nebular emission features and CO absorption from Red Super Giant (RSG) stars. This work, and the estimated ages of the brightest star clusters, indicates that the nuclear star formation is not well represented as a single instantaneous starburst event.

Note that the cluster-based studies do not provide a complete view of the starburst, as they are mainly observed in regions of lower-than-average extinction444The typical measured extinction towards the clusters is — 6 magnitudes, see Melo et al. (2005). toward the nucleus (Smith et al., 2006a). Nor does this work strongly constrain the star formation rate between 10 — 60 Myr ago. Old clusters in this age range, of similar mass to the bright clusters, would be considerably fainter and thus harder to detect than the Myr-old clusters555Note that this statement, and the previously mentioned cluster ages, is based on single star evolutionary models. Mas-Hesse & Cerviño (1999) argue that if binary systems are taken into account the strongly ionizing phase of a single coeval cluster can last up to Myr after formation, and that the cluster can simultaneously produce nebular emission and have a substantial number of RSGs. If this is indeed the case then the clusters identified as having ages of Myr may be as old as Myr., due to the rapid decrease of the ionizing and bolometric output at ages Myr. UV spectroscopic studies of nearby low-obscuration starbursts have shown that the massive star population within the clusters is typically younger the the field massive star population (see e.g.  Tremonti et al., 2001; Chandar et al., 2005). The probable cause for this effect is the gradual dissolution of the massive star clusters on times scales of — 30 Myr (for a recent review see de Grijs & Parmentier, 2007, and references therein). We are not aware of any similar UV study of the field stellar population in M82’s nucleus, so there is no direct evidence for the presence of star formation intermediate in age between the 10 Myr-old clusters and the Myr-old clusters. Nevertheless it remains very plausible that significant star formation continued in between the two periods.

Another reason to suspect that the nuclear burst has been active for longer than Myr is the dynamical age of the H cap, a large nebular feature associated with the superwind at a height kpc above the mid-plane of the M82 (Devine & Bally, 1999; Lehnert et al., 1999). Its true space velocity is unknown, but assuming a mean velocity of equivalent to that of the H filaments in the lower wind (McKeith et al., 1995; Shopbell & Bland-Hawthorn, 1998) it would have taken Myr to reach its current location if its origin was within the disk of the galaxy.

We thus choose to treat M82’s nuclear starburst as having experienced a uniform rate of star formation since the onset of the burst years ago, with negligible star formation before that time. We will treat as a model parameter in the range — 50 Myr.

In Table 4 we present the bolometric luminosity, and total SN plus stellar wind energy and mass injection rates (all per unit star formation rate) for continuous star formation models calculated using version 5.1 of Starburst99 population synthesis code (Leitherer et al., 1999). The star formation rate is for a Salpeter IMF between the mass limits of 1 — . For a Kroupa (2001) IMF between the mass limits 0.01 — producing the same total mass of massive stars (mass range 8 — ) as in the Salpeter IMF model, the star formation rate is a factor 1.76 times larger. The Geneva stellar evolutionary tracks for an initial metallicity equal to Solar, and with high stellar wind mass loss rates, were used.

(Myr) () () () ( yr) () ( yr)
(1) (2) (3) (4) (5) (6) (7)
10 0.19 6.1 1.1
20 0.26 5.1 1.3
30 0.29 4.7 1.4
40 0.31 4.5 1.4
50 0.31 4.4 1.3

Note. – Column 1: Age since the onset of continuous star formation. Column 2: Starburst radiated bolometric luminosity per unit star formation rate. Column 3: Mechanical energy injection rate due to SNe and stellar winds per unit star formation rate. Column 4: Mass injection rate due to SNe and stellar winds per unit star formation rate. Note that the values in Columns 2 – 4 are specific to a Salpeter IMF between the lower and upper mass limits of 1 — . Column 5: Star formation rate derived from , where is the observed IR luminosity of M82. Column 6: The mechanical energy injection rate associated with this SFR. Column 7: The mass injection rate associated with this SFR. See § 3.3.1 for details on the stellar evolutionary tracks and stellar wind models used.

Table 4Continuous star formation models for M82.

The total 8 — IR luminosity of M82 is , based on the latest re-calibrated 12 – IRAS fluxes (Sanders et al., 2003). IR emission dominates M82’s spectral energy distribution, so by equating to the starburst bolometric luminosity for each assumed age yields star formation rates between 4 — (Salpeter IMF). These values and the associated SN plus stellar wind energy and mass injection rates are tabulated in Table 4. By way of comparison the formulae in Kennicutt (1998), which assume CSF with a large age, give a star formation rate of for M82 based on , when corrected to the 1 — mass range Salpeter IMF we use.

These star formation rates are not unreasonable. Over the time scale of interest the net mass of stars formed is at maximum % of the dynamical mass of within pc (Götz et al., 1990), even if we assume a Kroupa IMF extending down to a lower mass limit of .

3.3.2 The total size of the starburst region

The brightest region of diffuse non-thermal 5 GHz radio emission in M82 subtends an angular extent of (), the same region occupied by the compact radio sources presumed to be young SN remnants (Kronberg et al., 1985; Muxlow et al., 1994). This is similar to the size of the region occupied by the brightest star clusters seen in the NIR, (, McCrady et al. 2003). These major-axis angular diameters closely match the base of the superwind as measured by the region of split optical emission lines, which has a diameter of Götz et al. (1990).

If we approximate the starburst region as a flat cylinder of diameter and total height then the starburst volume and surface area in each case is pc, pc (5 GHz) and pc , pc (NIR). The equivalent spherical CC model should then have a starburst region radius of pc or pc (Appendix B).

The starburst in M82 is possibly concentrated along the stellar bar (see e.g.  Achtermann & Lacy, 1995; Westmoquette et al., 2007), and hence is more geometrically complex than a uniform disk. Furthermore star formation may not have occurred simultaneously throughout the full starburst region, but evidence for inwardly or outwardly propagating star formation within pc of the nucleus remains inconclusive (Förster Schreiber et al., 2003a). Exploring the full significance of these effects requires 3-dimensional hydrodynamical simulations that are outside the scope of this paper.

However we can investigate the significance of concentrating the supernova activity within a smaller disk-like region than the two considered above. We adopt the starburst disk geometry used in Strickland & Stevens (2000) which has pc and pc. This produces significantly lower starburst volumes and surface areas: pc and pc. The equivalent spherical starburst radius is pc.

For the majority of the models we consider we adopt the first starburst region geometry described above: a disk of diameter pc and total height pc (or pc for the scaled CC models). The effect of the size of the starburst region is explored using the scaled CC model with radii pc and pc.

3.3.3 The effective energy and mass injection rates

The most significant contributors to the energy and mass return rates in a starburst event are core-collapse SNe and the stellar winds from the massive stars. Stellar winds are thought to contribute an order of magnitude less mechanical energy than the SNe (assuming equal thermalization efficiencies), but they may contribute up to % of the total mass returned to the ISM (Leitherer, Robert, & Drissen, 1992; Leitherer & Heckman, 1995). As described above, we used the Starburst99 population synthesis code to find the combined SNe plus stellar wind energy and mass injections rates for a continuous burst of star formation as a function of age. Note that the values of and given in Table 4 have not been corrected for any initial mechanical energy losses due to radiation before thermalization occurs (the thermalization efficiency ) or for mass loading by the addition of cold ambient ISM (with a rate ).

The thermalization efficiency and degree of mass loading are expected to be a function of the local gaseous environment the massive stars find themselves in, as well as the local supernova rate per unit volume and the density of massive stars. Low ambient gas density and/or high SN rate per unit volume should lead to a higher hot gas filling factor, which in turns increases the thermalization efficiency by reducing the initial radiative energy losses (Heckman et al., 1990). However, at very high stellar densities the SN and stellar wind ejecta may itself become dense enough to cause drastic energy losses, e.g. within the central regions of super star clusters (Silich et al., 2004; Tenorio-Tagle et al., 2007).

To account for these factors we assume that only a fraction of all the SNe and stellar winds occur in regions where conditions lead to non-negligible thermalization, and the mean thermalization efficiency in these regions is what we define as . Essentially is a modification to the total star formation rate of the starburst to account for localized spatial variations within the starburst of the SN and stellar wind feedback. We will refer to as the participation fraction, as it is effectively the fraction of the total star formation that is actually involved in driving the superwind. For the remaining fraction 1- of the SNe and stellar winds we assume that radiative energy losses are so significant that no mechanical energy is supplied to the starburst region, and that the ejecta never becomes part of the outflow. We treat as a model variable, and refer the reader to Appendix C for a discussion of plausible values for given the gas content of the starburst and the star cluster mass function.

The net energy injection rate into the hot volume-filling gas in the starburst region is


The net mass injection rate is


where the mass loading factor is a parameter commonly used to indicate the relative amount of mass loading. If then the wind is not mass loaded. In this paper we consider only mass loading that operates within the starburst region, i.e. central mass loading in the terminology adopted by Suchkov et al. (1996) and Strickland & Stevens (2000).

Written as a function of the thermalization efficiency, the mass loading factor and participation factor, the central temperature, pressure and density are




The central temperature is thus linearly proportional to the ratio of the mechanical energy thermalization efficiency to the mass loading factor, and independent of the size of the starburst region or participation factor . Thus it is an excellent diagnostic tool for quantifying the effects of feedback, provided that this temperature can be accurately measured.

The terminal velocity of the wind fluid


is also a function of . An observational determination of is beyond the capabilities of current X-ray observatories, but the X-ray Microcalorimeter Spectrometer (XMS) instrument on the future International X-ray Observatory (IXO) will provide the high resolution non-dispersive X-ray spectroscopy necessary to directly measure the velocity of the hot phases in superwinds.

3.3.4 The metal abundance within the wind fluid

The mixed SN ejecta and stellar wind material will be enriched in heavy elements compared to the cooler ambient ISM within the starburst region. As the emissivity of the thermal X-ray emission from a hot plasma depends strongly on the metal abundance of the hot plasma it is necessary to consider both SN-related enrichment and dilution via mass loading.

If we assume mass-loading is associated with efficient mixing, then the elemental abundance of an element in the wind fluid is


where is the abundance of the element in the SN plus stellar wind ejecta and is the abundance in the ambient ISM, and we have made use of Equation 7 which defines .

The primary elemental coolants for plasmas with temperature K are oxygen, neon and iron (e.g. see Sutherland & Dopita, 1993). Iron is of particular importance for this work given the keV and keV lines from helium-like and hydrogen-like iron ions. Elements such as S, Ar and Ca also produce weaker X-ray lines at energies keV that may also serve as a diagnostic of the wind-fluid.

Limongi & Chieffi (2007) present the IMF-averaged SN ejecta (both SN II and SN Ib/c) element abundances with respect to Solar666All abundances quoted in this paper are relative to the scale defined in Anders & Grevesse (1989).. Their stellar evolution models include mass loss due to stellar winds, but the quoted yields do not include the hydrogen and lighter elements returned to the ISM by the winds. They find , , , , and . Note the low O/Fe abundance ratio compared to the IMF-averaged yields presented in Gibson et al. (1997), which is most probably a consequence of wind mass loss in the stellar evolution models used by Limongi & Chieffi (2007).

The time-averaged nature of an IMF-averaged yield is appropriate for continuous star formation with an age of Myr, but the abundances in the ejecta for a younger or instantaneous burst may differ from that in Limongi & Chieffi (2007). To investigate this, and the effect of using a different set of SN yields, we computed the abundances of the combined SN ejecta and stellar wind material from the yields given by Starburst99. Starburst99 includes wind yields by combining the stellar evolutionary wind mass loss rates with observed stellar surface abundances. SN yields are based on Woosley & Weaver (1995) for SN type II, but note that the the stellar evolution models that Woosley & Weaver used did not in general include wind mass loss. Starburst99 also assumes that stars that become Wolf-Rayet stars (initial mass ) are assumed to explode as SN type Ib/c and contribute no metals to the final yields. For a continuous star formation model the SN plus stellar wind ejecta abundances predicted by Starburst99 are relative constant for ages Myr after the onset of the burst, at levels of , , and (Starburst99 does not calculate the yields of Ne, Ar and Ca).

Given the uncertainties we will adopt for these models. For computational convenience we will further assume the same abundance relative to Solar for the other elements.

The stellar and H ii region O and Ne abundances in M82 are Solar or higher (Achtermann & Lacy, 1995; Origlia et al., 2004; Smith et al., 2006b). The are very few measurements of iron abundances, the only one we are aware of is that of Origlia et al. (2004) who derive a stellar iron abundance of from NIR spectroscopy. The total baryonic mass of M82 is based on its peak rotational velocity (see Strickland et al., 2004a). Galaxies of this mass have a mean nebular oxygen abundance of (on the Anders & Grevesse 1989 scale), based on the galaxy mass-metallicity relationship (Tremonti et al., 2004). Thus the measured metallicity of the stellar population and the warm ISM in M82 are not unusual given its mass, nor do they suggest that the ambient ISM has been significantly enriched with nucleosynthetic products from the current starburst. We adopt for our models, again assuming the same value for all other elements.

Figure 4.— (a) The broad band cooling functions used in the model calculations. The total cooling due to metals for unit metal abundance with respect to Solar abundance (1) and cooling from hydrogen and helium only (2) for non-equilibrium cooling are both from Sutherland & Dopita (1993). The APEC plasma code (Smith et al., 2001) was used to calculate the emission in the hard X-ray energy band ( – 8 keV) from metals (3) and hydrogen plus helium (4), assuming collisional ionization equilibrium. (b) X-ray cooling curves per unit metal abundance for helium-like and hydrogen-like ions of S, Ar, Ca and Fe in comparison to the 2–8 keV cooling curve for all metals (3).

3.4. Broad-band cooling and X-ray line emission

The luminosity per unit volume () due to the thermal emission from an element of gas of density , temperature and metal abundance is


where is the mass fraction of hydrogen for a Solar abundance plasma, is the mass of a hydrogen atom and is the cooling function for the specific energy band or line transition being considered.

For broad-band X-ray emission, or the total cooling rate, the cooling function for a hot optically thin plasma may be split into a metallicity-independent component and metallicity-dependent component (evaluated for Solar abundance), such that , where is the metal abundance with respect to the Solar. This assumes that the relative metal abundances are Solar, as otherwise the contribution from each element must be considered separately.

To assess the total cooling rate we use the zero-field non-equilibrium cooling functions from Sutherland & Dopita (1993). We calculated the metallicity-dependent cooling from their published Solar and zero metallicity cooling functions. At temperatures above K (typical of the vast majority of the models we will consider) the difference between equilibrium and non-equilibrium cooling functions is negligible.

To calculate the cooling functions for the — 8 keV hard X-ray band and the He and Ly like lines of S, Ar, Ca and Fe we used version 1.3.1 of APEC hot plasma code (which assumes collision ionization equilibrium), accessed using the Xspec spectral fitting program (Arnaud, 1996; Smith et al., 2001). For the line emission from any specific element the cooling function is an exactly linear function of the metallicity , i.e. .

In Fig. 4 we plot the total and 2 – 8 keV cooling functions (split into their metallicity-independent and metallicity-dependent components) that we use in our calculations. The line-specific cooling functions for hydrogen-like and helium-like sulphur, argon, calcium and iron for unit metal abundance are also shown.

The X-ray line emissivities we use apply to hot gas in collisional ionization equilibrium (CIE). In our the calculations using the scaled CC model we evaluate the following expression in each radial shell as we numerically integrate outwards:


If the ionization parameter then the plasma is likely to be in CIE (Masai et al., 2002). If s then the outflowing gas is over-ionized and in recombining non-equilibrium. Masai et al. show that the broad band X-ray emissivity of such a recombining plasma is less than that of a CIE plasma. The ionization parameter scales with thermalization efficiency and mass loading proportionally as , as increasing the mass loading or decreasing the thermalization efficiency acts both to increase the plasma density and decrease the outflow velocity.

For parameters appropriate to M82 we find that the plasma is typically in CIE at radii even when . Even if we integrate the emission out to radii the total luminosity is dominated by emission from within , justifying our use of CIE emissivities. The radial variation in and the cumulative emission integral shown in Fig 3c demonstrate this for one set of possible model parameters. This result is consistent with the alternative ionization equilibrium and Coulomb relaxation timescale calculations presented in Paper I, which also indicated that the hot plasma within the starburst region should be in ionization equilibrium.

For a single temperature plasma the luminosity is the product of the emission integral with the cooling function . The cumulative volume emission integral is defined as . For a CC-style radial wind this is dominated by material within the relatively dense starburst region, as shown in Fig 3c. Thus the total emission integral is proportional to


which again increases as either the mass loading factor is increased, or as the thermalization efficiency is decreased. However such variations in and/or are also associated with a decrease in the central temperature, and for increased mass loading a decrease in the mean metal abundance, so that the net luminosity can ultimately decrease even as increasing or decreasing monotonically increase the net emission integral.

3.5. Observational constraints

The primary aim of this project is to assess the range of starburst properties that are consistent with the hard X-ray Fe-line-emitting plasma in M82, and which models are deemed acceptable depends on which constraints are used — this is unavoidable. These are not an arbitrary set of observational and conceptual constraints but have been chosen based on realistic physical expectations.

Based on the results and discussion presented in § 2 we adopt the following as the full observational constraints on the properties of the wind fluid within a radius of 500 pc of the center of M82.

  • The total X-ray luminosity of the wind fluid in the – 8 keV energy band must be less than or equal to the broad band luminosity estimated from the ACIS-S observation for the diffuse hard X-ray emission (Table 1), i.e. .

  • The Fe He line luminosity must be consistent with the observed value given in Table 1. In order to be conservative in our estimation of the constraints the observations place on the thermalization efficiency and mass loading, i.e. to accept the largest number of models matching the observational constraints, we adopt the weakest constraints on the allowed value of , i.e. using the ACIS-S observation rather than the better constrained values from the ACIS-I or XMM-Newton data sets. Models producing Fe line luminosities within 68.3% confidence region from the ACIS-S observation, , will be judged as matching the constraint on the keV line emission. This measurement of is slightly larger than, but still consistent with the line luminosity estimated from the ACIS-I and XMM-Newton observations.

  • The Fe Ly luminosity must be less than or equal to the observationally-determined upper limits, of which the best determined value is from the XMM-Newton observation. We specify that should be for a model to be judged to be consistent with the observational constraints (Table 1).

  • The luminosity of the predicted He and Ly-like lines of S, Ar and Ca must be less than or equal to the observed values and/or 99% upper limits. Here the best observational constraints are from the XMM-Newton data (Table 3). Acceptable models must then satisfy all of the following: , , , , , and . Upper limits are used as the S, Ar and Ca-line-emitting plasma is not necessarily the same phase as the Fe-line-emitting plasma.

Models that (a) have total radiative losses less than 30% of the net energy injection rates (see Appendix A) and (b) predict X-ray luminosities consistent with the full set of constraints given above are considered as having successfully matching the full set of constraints.

We will consider the effect of using a different set of observational constraints. These alternatives deliberately ignore what we know about the iron line emission (the iron lines are not used as constraints) and instead attempt to fully account for all of either the broad-band hard X-ray luminosity or the detected S, Ca and Ar line luminosities:

  • The total X-ray luminosity of the wind fluid in the – 8 keV energy band must be within the 3 envelope given in Table 1 for the model to considered successful, i.e. . In addition the the total radiative losses must be less than 30% of the net energy injection rate. We shall refer to these constraints as the broad-band-only constraints.

  • The predicted line luminosities for the S He, S Ly, Ar He and Ca He must be within of the values given in Table 3 for the best fit to the XMM-Newton spectrum. The allowed limits on the Ar Ly and Ca Ly line luminosities are the same as given for the full set of constraints. In addition the model must predict broad band and total radiative energy losses less than the upper limits also used in the full set of constraints: and total radiative losses less than 30% of the net energy injection rate. We shall refer to these constraints as the SArCa constraints.

We allow larger ranges in the broad-band or line luminosity () as we found that in practice that it is relatively hard to find regions of parameter space where models can satisfy the broad-band-only or SArCa constraints.

4. Simulations and the results thereof

Figure 5.— Predicted luminosities and temperatures as a function of the input thermalization efficiency and mass loading factor , as predicted using the 2-dimensional hydrodynamical simulations with parameter set A (see § 4.1). (a) The logarithm of the Fe He luminosity. The grey-shaded region labeled F shows the region of parameter space that matches the full set of observational constraints discussed in § 3.5. The region labeled BB denotes the region satisfying the broadband-only constraints, while the region labeled SArCa satisfies the S, Ar, and Ca line constraints. The region labeled NA denotes the region of parameter space where the adiabatic assumption ceases to the valid. (b) As (a), except the labeled contours show the logarithm of the luminosity in the – 8 keV energy band. (c) As (a), except the labeled contours show the logarithm of the S Ly line luminosity. (d) As (a), except the contours show the logarithm of the central plasma temperature in Kelvin.

4.1. Acceptable values of and .

Mass loading factors of or thermalization efficiencies of mark the practical upper and lower limits on these parameters for a M82-like starburst. Each of these values on its own will reduce the central temperature to K, a temperature at or below which the emissivity of the Fe He emission drops dramatically (Fig. 4).

To explore this region of the mass loading and thermalization efficiency parameter space we adopt the following values for the other important model parameters: , a disk-like starburst with pc and pc ( pc), and . We treat the starburst as continuous star formation at a rate of beginning 30 Myr ago, which leads to and . We shall refer to this combination of input model parameters as parameter set A (see Table 5).

Parameter Set SF Mode/Age
() ( yr) (pc) ()
(1) (2) (3) (4) (5) (6) (7)
A CSF / 30 Myr 1.4 300.0 5.0 1.00
B CSF / 10 Myr 1.1 300.0 5.0 1.00
C SIB / 10 Myr 2.2 300.0 5.0 1.00
D CSF / 30 Myr 1.4 125.5 5.0 1.00
E CSF / 30 Myr 1.4 404.9 5.0 1.00
F CSF / 30 Myr 1.4 300.0 2.5 1.00
G CSF / 30 Myr 1.4 300.0 7.5 1.00
H CSF / 30 Myr 1.4 300.0 5.0 0.71
I CSF / 30 Myr 1.4 300.0 5.0 0.50

Note. – Column 1: Adopted name of combinations of model parameters. Column 2: The mode of star formation assumed, either continuous star formation (CSF) or an single instantaneous burst (SIB) beginning or occurring the specified time before the present. Columns 3 and 4: The mechanical energy and mass injection rates (See § 3.3.1) associated with the assumed star burst. Column 5: The effective starburst radius adopted in calculations using scaled CC model. Column 6: The assumed metal abundance of heavy elements relative to Solar (see § 3.3.4). Column 7: The assumed participation fraction (see § 3.3.3 or Appendix C).

Table 5Combinations of starburst model parameters explored in the simulations.

We initially computed models covering the full range of parameter space. To prevent the total number of simulations from becoming excessive we considered only a limited number of distinct values for and , specifically and , where N ranges from 0 to 3 in increments of 0.25, resulting in the computation of 169 models using the 2-dimensional hydrodynamical code described in § 3.1.

We find that only limited regions of parameter space yield models that satisfy any one of the sets of observational constraints we chose to apply (Fig. 5). Furthermore, the regions in : parameter space satisfying the different sets of constraints are well separated from one another.

Models that satisfy the full set of observational constraints (labeled F in Fig. 5) have moderately high thermalization efficiencies () and moderately low but non-negligible mass loading (). These models reproduce the observed Fe line fluxes and limits without overproducing S, Ar or Ca line emission, but in doing so they can only account for a small fraction (%) of the total observed diffuse luminosity of in the – 8 keV energy band (Fig. 5b).

The range of central thermal pressures found in these matching models, K cm (given in Table 6) is consistent with estimates of the pressure within the starburst region obtained from optical spectroscopy: K cm (Heckman et al., 1990; Smith et al., 2006b; Westmoquette et al., 2007).

To match the observed S, Ar and Ca line fluxes using a wind model (the region labeled SArCa), models require a lower thermalization efficiency and a comparable or slightly higher mass loading factor (, ).

To reproduce the observed diffuse X-ray luminosity in the – 8 keV energy band with thermal emission alone requires both high thermalization efficiency and high mass-loading factors: and (the region labeled BB in Fig. 5). Models in this region of parameter space fail to match the full set of observational constraints because they typically predict more S, Ar or Ca emission (the S Ly luminosity is shown in Fig. 5c) than is observed when the Fe line luminosity is correctly reproduced, or too little Fe He emission when they predict less than or equal to the amount of S, Ar or Ca line emission observed.

Thus we can only account for % of the observed – 8 keV diffuse X-ray luminosity with pc with the sum of the predicted thermal emission from the starburst region and the Inverse Compton X-ray emission (Paper I). We conclude that either our estimate of the Inverse Compton luminosity is in error, or there must be an additional non-thermal or quasi-thermal source of hard X-ray continuum emission within the starburst region (e.g.  Dogiel et al., 2002; Masai et al., 2002). Solving this enigma will require higher resolution spectroscopy to establish the shape of the continuum and the ionization state of the line emission.

Models where the total radiative losses are higher than 30% of the energy injection rate, and the adiabatic assumption is no longer valid, are labeled NA in Fig. 5. Note that for these models the temperatures and X-ray luminosities we predict are not valid. The location of this region in space in consistent with the discussion in Appendix A. For M82’s star formation rate and starburst region size we find that as long as then the adiabatic assumption is applicable to either a disk-like or spherical starburst.

Parameters that give rise to highly radiative winds can be rejected even though we can not predict their exact X-ray luminosities. Cooling is so efficient under such conditions that the central temperature K and there is no possibility of Fe He emission (e.g.  Tenorio-Tagle et al., 2007).

Figure 6.— As Fig. 5, except using a finer grid of models distributed linearly in a smaller region of parameter space. The numbered contour levels are of log (a), log (b), log (c) and log (d).

To further refine our knowledge of the exact region of the allowed parameter space for M82’s starburst we ran a further set of full hydrodynamical simulations a subset of the full grid described above. We used a linearly-spaced grid, varying between 0.1 and 1.0 in steps of and between 1.0 and 3.0 in steps of (a total of of 494 models).

The predicted regions of allowed parameter space for parameter set A are shown in Fig. 6. This figure also demonstrates that the exact shape of the allowed region using the full set of constraints is most-strongly influenced by the combination of the constraints on the Fe He line and the upper limit on the S Ly line777For this reason we have not provided plots of the S He, Ar or Ca line luminosities as a function of and .. For a given value of the thermalization efficiency the lower boundary of acceptable mass loading factors is set by requirement to achieve sufficient Fe He emission (the necessary condition is ), while the maximum is set by the requirement not to exceed the observed S Ly luminosity ().

For a given value of the highest allowed is a less likely solution than a lower allowed value of . The highest allowed values of produce the entire observed S Ly flux in the wind fluid, without leaving leeway for S Ly emission from the soft X-ray-emitting plasma within the starburst region.

Were we to change the S Ly constraint from an upper limit to a specific range in values around the observed line luminosity the resulting allowed region would be a narrow strip along what is currently the upper margin of the full allowed region. Adopting such a constraint would be somewhat ad-hoc, given that the S He and Ar He lines are equally well detected as the S Ly line (Table 3). Furthermore, an attempt to match the S He, Ly and Ar He line fluxes, rather than treating them as upper limits, simultaneously with the Fe He line flux leads to no acceptable models being found. Ultimately this is because the S, Ar and Ca line fluxes and upper limits are indicative of a plasma temperature (see Table 2 or the SArCa region on Fig. 6d), yet such a low temperature plasma produces too little Fe He emission to match the observations.

Figure 7.— (a) A comparison between the region of allowed parameter space for M82 predicted by the scaled CC model (shaded regions) in comparison to the regions predicted by the full two-dimensional hydrodynamical simulations (regions enclosed by dot-dashed lines, previously shown in Fig. 5), using parameter set A. (b) As in (a), except the comparison uses the finer grid of models distributed linearly in a smaller region of parameter space.

As shown in Fig. 7, these same general regions of acceptable parameter space are reproduced in simulations making use of the scaled CC model, although in detail there are minor difference caused by the % errors in the luminosity predicted by the scaled CC model compared to the multi-dimensional hydrodynamical simulations (Appendix B).

4.2. Other input model parameters.

Given the computational cost of the full hydrodynamical simulations it is preferable to use the scaled CC model to explore the effect of varying other model parameters such as starburst region size, metal abundance and participation factor. Our experience with parameter set A indicates that the scaled CC model reproduces the allowed values of and to an accuracy of for the full set of observational constraints that include the Fe lines.

Figure 8.— (a) A comparison between the region of allowed parameter space for M82 predicted by the scaled CC model using parameter set B (shaded regions) in comparison to the equivalent regions for parameter set A (dashed lines). (b) As in (a), except the comparison is between parameter set C and parameter set A. Parameter sets B and C differ from A in the assumed starburst star formation mode and age.
Figure 9.— (a) A comparison between the region of allowed parameter space for M82 predicted by the scaled CC model using parameter set D (shaded regions) in comparison to the equivalent regions for parameter set A (dashed lines). (b) As in (a), except the comparison is between parameter set E and parameter set A. Parameter sets D and E differ from A in adopting respectively smaller and larger starburst effective radii .
Figure 10.— (a) A comparison between the region of allowed parameter space for M82 predicted by the scaled CC model using parameter set F (shaded regions) in comparison to the equivalent regions for parameter set A (dashed lines). (b) As in (a), except the comparison is between parameter set G and parameter set A. Parameter sets F and G differ from A in adopting respectively smaller and larger supernova and stellar wind metal abundances .
Figure 11.— (a) A comparison between the region of allowed parameter space for M82 predicted by the scaled CC model using parameter set H (shaded regions) in comparison to the equivalent regions for parameter set A (dashed lines). (b) As in (a), except the comparison is between parameter set I and parameter set A. Parameter sets H and I differ from A in adopting lower values of the participation factor and respectively, i.e. effectively lower star formation rates.

Parameter sets B through I are shown in Table 5. These alternate parameter sets explore the influence of the chosen star formation mode and age, starburst region size, supernova and stellar wind metal abundance, and participation factor on the region of allowed parameter space. All of these calculations were performed using the scaled CC model alone. Predicted central density, pressure and temperature, along with the terminal velocity, total cooling, broad-band hard X-ray ( – 8 keV band) luminosity and Fe He and S Ly line luminosities for and are shown in Table 7.

Parameter Set Range in parameter consistent with full constraints
(1) (2) (3) (4) (5) (6) (7) (8)
AaaCalculated using 2-dimensional hydrodynamical simulations. All other parameter sets were calculated using the scaled 1-dimensional Chevalier & Clegg model. 52 0.55 – 1.0 1.4 – 2.2 7.56 – 7.82 7.15 – 7.38 34.31 – 34.54 1555 – 2097
A 75 0.40 – 1.0 1.2 – 2.3 7.51 – 7.84 7.08 – 7.40 34.23 – 34.55 1470 – 2164
B 26 0.75 – 1.0 1.8 – 2.3 7.49 – 7.66 7.13 – 7.24 34.29 – 34.39 1434 – 1745
C 32 0.45 – 1.0 1.0 – 1.4 7.64 – 7.87 7.22 – 7.47 34.37 – 34.62 1709 – 2239
D 21 0.45 – 1.0 1.0 – 1.4 7.67 – 7.87 7.80 – 8.05 34.20 – 34.44 1778 – 2240
E 60 0.55 – 1.0 1.7 – 2.6 7.48 – 7.76 6.94 – 7.17 34.35 – 34.58 1426 – 1976
F 47 0.60 – 1.0 1.9 – 2.6 7.50 – 7.72 7.26 – 7.43 34.41 – 34.58 1452 – 1874
G 76 0.30 – 1.0 1.0 – 2.0 7.50 – 7.87 6.96 – 7.37 34.11 – 34.52 1452 – 2240
H 43 0.65 – 1.0 2.1 – 2.8 7.47 – 7.70 7.15 – 7.29 34.30 – 34.44 1409 – 1829
I 0

Note. – Column 1: Parameter set name. Column 2: is the number of models that match the full set of observational constraints, out of the 494 models calculated for each parameter set. Columns 3 and 4: The range of thermalization efficiencies and mass loading factors found to match the full set of observational constraints. The remaining columns show the range in a variety of physical properties of the M82 superwind found in the models that match the full set of observational constraints: central temperature (column 5, in Kelvin), central pressure (column 6, in units of K cm), rate of change of superwind momentum (column 7, units g cm s) and wind terminal velocity (column 8, units ).

Table 6Superwind properties that match the full set of observational constraints.

Parameter set B differs from the previously discussed parameter set A only in that is assumes that ongoing star formation began 10 Myr ago, rather than 30 Myr ago. As the younger stellar population is intrinsically more luminous, less star formation is required to match the IR luminosity of M82 and consequently the mechanical energy and mass injection rates in parameter set B are lower than in parameter set A. Furthermore the ratio of energy to mass injection is lower, resulting in a lower central temperature of K for unit and compared to K in parameter set A, and a lower temperature than parameter set A for any given value of and . Cumulatively these differences result in a significant shift in the location of the the models in parameter space that match either the full or SArCa constraints (see Fig. 8a).

Generation of plots equivalent to Fig. 6 (not shown) reveals that the region matching the full set of constraints is again determined by the intersection of the Fe He and S Ly observational constraints. Furthermore, at a given value of and the predicted broad band and (in particular) line luminosities are lower for parameter set B than for parameter set A, with the least difference in luminosity occurring at .

Application of Equation 15 indicates that parameter set B yields emission integrals consistently dex lower than parameter set A for all values of and . For the Fe and S lines the lower temperature of parameter set B when results in a higher cooling function that offsets the reduced emission integral because the lower temperature is closer to the peak of the line-specific cooling functions, such that the predicted values of and are comparable to those of parameter set A (Table 7). But over much of the parameter space the lower temperatures in parameter set B result in lower line emissivities that in combination with the lower emission integral result in lower line luminosities than in parameter set A. Ultimately the regions in parameter space that match the full and SArCa constraints in parameter set B require higher mass loading and higher thermalization efficiency than in parameter set A in order to increase the emission integral without further decreasing the temperature. Increasing strongly increases the net emission integral as , although it slightly decreases the plasma temperature (which must then be offset by increasing the thermalization efficiency). In contrast, reducing the thermalization efficiency is counter productive as the increase in is counteracted by the decrease in temperature that consequently reduces the line emissivities.

Parameter set C represents the star formation in the nucleus M82 as a single instantaneous burst occurring 10 Myr ago. Compared to parameter set A this results in a higher mechanical energy and mass injection rate, although the ratio of energy to mass injection is much closer to that in parameter set A than the lower values of parameter set B.

The higher mass injection for parameter set C leads to higher predicted X-ray luminosities for a given value of and than for parameter set A, which has a significant effect on the regions of parameter space that match the various constraints (Fig. 8b). The region matching the – 8 keV broad band luminosity appears at much lower than in parameter set A, although in terms of absolute mass injection the location of this region is consistent between parameter set A and parameter set C.

Similarly, lower values of are required to match the full set of constraints, the allowed region extending between and . Models with are again excluded as possible matches to the observed properties of M82, but for different reasons than in the case of parameter sets A and B. For parameter set C models with and generate sufficient Fe He luminosity to match the required constraint, but these models predict a Fe Ly/He line ratio greater than the upper limit of 0.34 (i.e.  is too high), and hence do not match the full set of constraints. A combination of constraints shapes upper edge of the allowed region. In addition to the requirement not to overproduce the S Ly emission that we encountered with parameter sets A and B, the upper edge of the allowed region in parameter set C is also set where the Fe He luminosity exceeds the observational limits.

Parameter set D adopts the compact starburst region discussed in § 3.3.2. In terms of the effective radius of the starburst region parameter set D uses pc instead of pc, but is otherwise identical to parameter set A. Consequently the gas density within the starburst region () is considerably higher than in parameter set A () or even parameter set C (), and results in a significant increase in the X-ray emission despite the reduction in starburst region volume (see Equ. 15). The regions of parameter space over which parameter set D matches the various sets of constraints is very similar to that parameter set C (Fig. 9a) for the same reasons we discussed above. The region of parameter space that matches the full set of constraints is smaller than in parameter set C because of the greater increase in luminosity caused by the higher gas density, and no region of matches the SArCa constraints. It is possible that with tighter observational constraints on and there would be no region of parameter space that would match the full set of observational constraints.

The lower volume of the starburst region adopted in parameter set D also leads to higher central pressures. Models that match the full set of observational constraints have , values that are 2 – 4 times higher than the highest independent observational estimates of the thermal pressure in the core of M82. This leads us to suspect that parameter set D is not a valid physical representation on the M82 starburst.

Parameter set E adopts the largest of the possible starburst region geometries discussed in § 3.3.2, for which the effective radius is pc and the resulting central number density is . The region that matches the full set of observational constraints now moves toward higher and higher with respect to the same region in parameter set A (Fig. 9b), for similar reasons to those we discussed regarding parameter set B. Increased mass loading is required to make up for the reduction in gas density due to the larger starburst, but the thermalization efficiency must also increase so as to maintain a central temperature that maximizes the Fe He emissivity (, for parameter sets A, B and E. See Table 6).

In § 3.3.4 we discussed the expected metal abundance of the merged SN ejecta and stellar wind ejecta. The calculations of Limongi & Chieffi (2007) predict a higher iron abundance for the SN ejecta () than the Starburst99 calculations that combine Woosley & Weaver (1995) SN yields with stellar wind ejecta (). Furthermore the Limongi & Chieffi calculations predict that sulphur, argon and calcium are less enriched than iron, while Starburst99 predicts that sulphur is more strongly enriched than iron. In all of the calculations we have discussed until now we adopted a uniform relative enrichment of all elements, . Parameter sets F and G explore the influence of the adopted metal abundance on the region of parameter space that match the observational constraints: parameter set F adopts while parameter set G adopts . Given the already large uncertainties in the yields we feel that exploring the effect of adopting different yields for different elements is beyond the scope of this paper.

The lower metal abundance in parameter set F leads to lower lower S, Ar, Ca and Fe line luminosities than for parameter set A. For the line luminosities are exactly half those in parameter set A, but the fractional difference becomes less as increases and the influence of the mass-loaded material on the net metal abundance of the wind fluid becomes larger (Equ. 12). The result in terms of the region of parameter space that matches the observational constraints is shown in Fig. 10a, and is very similar to that in parameter sets B and E that also generate lower line luminosities for different reasons.

Conversely the higher metal abundance in parameter set G leads to higher line luminosities than parameter set A, although they are generally less than those in parameter sets C and D (see e.g. Table 7). The region of matches is shown in Fig. 10b. The upper edge of the region that matches the full set of observational constraints is set by the upper limit on the S Ly luminosity, while the lower edge of this region is set by a combination of the Fe Ly/He line ratio constraint and the minimum acceptable Fe He line luminosity.

Parameter sets H and I explore the effect of varying the participation fraction (see § 3.3.3). Reducing has the effect of reducing the effective star formation rate without altering the energy per particle, and the emission integral and luminosities for a given value of and are equal to those from parameter set A except multiplied by a factor . For parameter set H, where , the factor 2 reduction in all luminosities moves the regions that match both the SArCa and full set of constraints to higher and higher (Fig. 11a). Parameter set I has and the reduction in luminosity is so significant that the region matching the SArCa constraints has moved off of the region shown in Fig. 11b to , . More significantly, there are no values of or that can match the full set of observational constraints when using parameter set I.

Parameter Set
() (pc)
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
A -0.80 7.22 8.02 2651 278 38.67 38.27 37.46 36.22
B -0.85 7.06 7.91 2341 278 38.56 38.17 37.45 36.24
C -0.59 7.40 7.98 2548 285 39.09 38.70 37.92 36.69
D -0.04 7.98 8.02 2651 121 39.17 38.75 37.93 36.74
E -1.06 6.96 8.02 2651 369 38.49 38.09 37.27 35.98
F -0.80 7.22 8.02 2651 278 38.51 38.11 37.16 35.92
G -0.80 7.22 8.02 2651 278 38.80 38.39 37.64 36.40
H -0.95 7.07 8.02 2651 272 38.37 37.97 37.16 35.92
I -1.10 6.92 8.02 2651 266 38.07 37.67 36.86 35.62

Note. – Column 1: parameter set name. Columns 2, 3 and 4: Logarithms of the central number density (), the central pressure divided by Boltzmann’s constant, and the central temperature (K) respectively for models where and . Column 5: Wind terminal velocity. Column 6: The radius at which the ionization parameter (see Equation 14) drops below . Columns 7, 8, 9 and 10: Logarithms of the plasma luminosity () within a radius of 500 pc. Columns are the total cooling luminosity, the – 8 keV energy band X-ray luminosity, the Fe He line luminosity and the S Ly line luminosity respectively.

Table 7Starburst region properties for and .

5. Discussion

In Strickland & Heckman (2007) we concluded that the diffuse Fe He emission from the central pc of M82 could have been produced by thermalized SN ejecta and stellar wind material from the starburst based on relatively simple estimates of the plasma properties of the emitting material.

The more rigorous approach to the interpretation of the diffuse hard X-ray emission taken in this paper strengthens that conclusion: the plasma responsible for the hard X-ray iron-line emission is the same gas that drives the larger-scale superwind in this starburst galaxy. Furthermore our method allows us to make quantitative statements regarding supernova feedback and the properties of M82’s starburst-driven superwind.

5.1. Supernova Thermalization (Heating) Efficiencies and Mass Loading

The primary aim of this paper was to establish what constraints could be made on the efficiency of supernova heating and the net rate of mass injection within a starburst region.

It is generally accepted that in more normal, or quiescent, star forming galaxies such as the Milky Way on average only % of the erg of initial kinetic energy per supernova ends up heating or moving the ISM, with the remainder being lost as radiation (see e.g.  Thornton et al., 1998; Efstathiou, 2000). Although it had long been suspected that the efficiency of supernova heating is a function of the local interstellar and star forming conditions there was no direct method by which to measure this thermalization efficiency in a starburst galaxy until, as described in § 1, the launch of the Chandra X-ray Observatory.

Theoretical arguments have been made for either lower or higher SN heating efficiencies in more actively star-forming galaxies. Steinmetz & Navarro (1999) argued that radiative energy losses must always be severe because massive stars form in regions of high gas density. In contrast Heckman et al. (1990) argue that thermalization must be efficient in starbursts because a high supernova rate per unit volume leads to an increased filling factor of the hot low density ISM, which in turn reduces the radiative energy losses experienced by subsequent SNe. Thus the SN thermalization efficiency depends on the SN-modulated porosity of the ISM (e.g.  McKee & Ostriker, 1977). Recent numerical simulations by Melioli & de Gouveia Dal Pino (2004) appear to qualitatively demonstrate such a time-dependent non-linear thermalization efficiency. Intriguingly they find that the thermalization efficiency (which they term the heating efficiency) can only be low (%) or very high (%) for any sig