Testing the Momentum-driven Supernova Feedback Paradigm in M31

Momentum feedback from isolated supernova remnants (SNRs) have been increasingly recognized by modern cosmological simulations as a resolution-independent means to implement the effects of feedback in galaxies, such as turbulence and winds. However, the integrated momentum yield from SNRs is uncertain due to the effects of SN clustering and interstellar medium (ISM) inhomogeneities. In this paper, we use spatially-resolved observations of the prominent 10-kpc star-forming ring of M31 to test models of mass-weighted ISM turbulence driven by momentum feedback from isolated, non-overlapping SNRs. We use a detailed stellar-age distribution (SAD) map from the Panchromatic Hubble Andromeda Treasury (PHAT) survey, observationally-constrained SN delay-time distributions, and maps of the atomic and molecular hydrogen to estimate the mass-weighted velocity dispersion using the Martizzi et al. ISM turbulence model. Our estimates are within a factor of 2 of the observed mass-weighted velocity dispersion in most of the ring, but exceed observations at densities $\lesssim 0.2$ cm$^{-3}$ and SN rates $>2.1\times 10^{-4}$ SN yr$^{-1}$ kpc$^{-2}$, even after accounting for plausible variations in stellar-age distribution models and ISM scale height assumptions. We conclude that at high SN rates the momentum deposited is most likely suppressed by the non-linear effects of SN clustering, while at low densities, SNRs reach pressure equilibrium before the cooling phase. These corrections should be introduced in models of momentum-driven feedback and ISM turbulence.


INTRODUCTION
Supernova (SN) feedback plays a critical role in galaxy formation by regulating the phase structure of the interstellar medium (ISM; McKee & Ostriker 1977;Mac Low & Klessen 2004;Joung & Mac Low 2006), launching galactic winds (Strickland & Heckman 2009;Heckman & Thompson 2017;Zhang 2018), accelerating cosmic rays (Drury et al. 1994;Socrates et al. 2008;Caprioli et al. 2010;Girichidis et al. 2016), and enriching the intergalactic and circumgalactic medium with metals Weinberg et al. 2017;Telford et al. 2019). Through a combination of these phenomena, feedback regulates the global star-formation efficiency of galaxies (Ostriker & Shetty 2011;Hopkins et al. 2012). Simulations imply that, without feedback, galaxies would rapidly convert cold gas into stars, resulting in up to a factor of 100 overproduction of stars compared to what is observed (Navarro & Benz 1991;Hopkins et al. 2011).
Unfortunately, current state-of-the-art cosmological simulations that study the evolution of galaxy population over cosmic time cannot resolve the spatial scales on which supernova remnants (SNRs) interact with the ISM. Even modern 'zoom-in' simulation of isolated galaxies can only marginally resolve SNRs (Hopkins et al. 2014(Hopkins et al. , 2018a, and properly resolved SNRs can only be obtained in simulations of smaller regions of the ISM disk (Gatto et al. 2017;Kim & Ostriker 2017;Karpov et al. 2020). This limitation motivated development of subgrid models of SN feedback at the resolution limit of simulations. Initial efforts to implement SN feedback in the form of thermal energy deposition were ineffective due to efficient radiative cooling in high-density starforming regions (Katz 1992). The quest to limit plasma cooling and runaway star-formation spawned a variety of subgrid models which employed techniques like delayed gas-cooling (Stinson et al. 2006;Governato et al. 2007Governato et al. , 2010, stochastic thermal feedback (Dalla Vecchia & Schaye 2012), an effective equation of state for a pressure-supported multi-phase ISM, with hydrodynamically decoupled wind particles (Springel 2000;Springel & Hernquist 2003;Oppenheimer & Davé 2006;Vogelsberger et al. 2014) These techniques ranged from being unphysical in nature to being inaccurate in the details of the SN-ISM coupling (Martizzi et al. 2015;Rosdahl et al. 2017;Smith et al. 2018).
More recent cosmological simulations (e.g., FIRE, Hopkins et al. 2018a,b) have explored subgrid models that deposit momentum, which unlike thermal energy, cannot be radiated away before impacting ambient gas (Murray et al. 2005;Socrates et al. 2008;Agertz et al. 2013). During the Sedov-Taylor phase of SNRs, the blast wave increases its momentum yield by a factor of 10-30 as it sweeps up ambient ISM. Later, it transitions into a cold, dense, momentum-conserving shell that ultimately merges with the ISM (Chevalier 1974;Cioffi et al. 1988; Thornton et al. 1998;Martizzi et al. 2015;Kim & Ostriker 2015;Karpov et al. 2020). This momentum budget per SN has been quantified by several realistic models of the ISM (e.g., Martizzi et al. 2015;Kim & Ostriker 2015;Li et al. 2015;Walch & Naab 2015;Karpov et al. 2020). It has been shown to effectively drive turbulence and winds, and reproduce key features of galaxies such as the Kennicutt-Schmidt relation and galactic morphologies (Martizzi et al. 2016;Smith et al. 2018;Hopkins et al. 2018a). More recent studies however have shown that the momentum deposition per SN depends sensitively on effects like SN clustering (Sharma et al. 2014;Gentry et al. 2017), entrainment of cold clouds (Pittard 2019), the abundance pattern of the ISM (Karpov et al. 2020), thermal conduction (El-Badry et al. 2019), enhanced cooling due to fluid-instability-driven mixing across the contact discontinuity (Gentry et al. 2019), the SN delay-time distribution (DTD) model (Gentry et al. 2017;Keller & Krui-jssen 2020), and pre-SN feedback via winds, photoionization and radiation pressure (Fierlinger et al. 2016;Smith et al. 2021) Observations that are specifically sensitive to SN feedback can identify a reliable subgrid model. Generally, cosmological simulations calibrate subgrid models to reproduce bulk properties of the galaxy population such as the stellar mass function and stellar mass to halo mass relation, but this necessarily limits the predictive power of the simulations (Schaye et al. 2015). Extragalactic multi-wavelength have served as useful references for setting subgrid model components such as SN rates (e.g. Mannucci et al. 2006), the efficiency of SN energy driving ISM turbulence (e.g. Tamburro et al. 2009;Stilp et al. 2013) and mass-loading in supersonic winds (e.g. Martin 1999;Veilleux et al. 2005;Strickland & Heckman 2009). However, the main source of uncertainty in modern subgrid models stem from a poor understanding of the SN-ISM interaction physics that originates on scales of 10-100 pc, which is beyond the reach for most distant surveys. In this respect, the resolved environments of Local Group galaxies provide detailed information available on stellar populations, ISM distribution and kinematics, and SNRs at the highest affordable spatial resolution. They are the ideal testing grounds for SN feedback models.
In this work, we test models of mass-weighted ISM turbulence predicted by SN momentum feedback models against observations of turbulence in M31, with a focus on the long-lived, prominent 10-kpc star-forming ring (Lewis et al. 2015;Williams et al. 2015). The proximity of M31 pushes the frontier of turbulence studies to < 100 pc, where the effects of feedback are spatiallyresolved, complete with detailed maps of the atomic ISM distribution (e.g. Braun 1991;Nieten et al. 2006;Braun et al. 2009), and spatially-resolved stellar age distribution (SAD) measurements with sensitivity down to masses ≈ 1.5 M obtained by the Panchromatic Hubble Andromeda Treasury survey (PHAT; Dalcanton et al. 2012;Williams et al. 2017). We can use these SADs to estimate SN rates by taking into account currently known constraints on the efficiencies of the different progenitor channels of core-collapse and Type Ia SNe, expressed in the form of SN delay-time distributions (DTDs) (Maoz et al. 2014;Zapartas et al. 2017;Eldridge et al. 2017). The SADs of the older stellar populations allow us to quantify the Type Ia SN rate as a function of location, which is important for a 'green-valley' galaxy like M31 (Mutch et al. 2011;Davidge et al. 2012), and lacks correlation with conventional star-formation rate tracers. This paper is organized as follows. In Section 2 we describe our analytical momentum-driven ISM turbulence model, and how we use stellar population and ISM data to constrain ISM densities, SN rates and velocity dispersion in the M31 ring. Section 3 describes the results of our analysis and checks on potential systematics, and in Section 4, we discuss the implications of these results on the assumed subgrid models of feedback used in cosmological simulations.

MODELING ISM VELOCITY DISPERSION IN M31
Here, we compare the observed non-thermal velocity dispersion in M31's neutral (H i) and cold ISM with the predicted turbulent velocity dispersion from the SN momentum-driven ISM turbulence model of Martizzi et al. (2016). Our calculations are supplemented by measured SN rates from the SAD of the PHAT survey and known forms of the SN DTD. We describe these efforts below. For all measurements, we assume that the distance to M31 is 785 kpc (McConnachie et al. 2005), and 1 = 3.78 pc at the distance of M31.
We restrict our analysis to the 10 kpc star-forming ring of M31 (Figure 1). We expect the main source of turbulence here to be star-formation, which we are mainly interested in testing, as opposed to other sources of turbulence observed in galaxies such as galactic spiral arms and magnetorotational instabilities (Tamburro et al. 2009;Koch et al. 2018;Utomo et al. 2019). Both atomic and molecular gas are most abundantly located and detected at high signal-to-noise along the ring (Braun et al. 2009;Nieten et al. 2006). Additionally, since the gas scale height in pressure-supported star-forming disks can vary with radius (Utomo et al. 2019), staying within the ring helps justify the use of a constant scale height in Eq. (6) and (7). We do however assess the impact of variable scale heights later on in Section 3.2.
In the following sub-sections, we describe our methodology for modeling the ISM velocity dispersion using momentum-driven turbulence from SNe, and the observations we use for comparison.

Momentum Injection by Supernovae
Following Martizzi et al. (2015) and Martizzi et al. (2016) -hereafter, M15 and M16 respectively -we assume that the non-thermal velocity dispersion in H i is a result of SNR momentum-driven turbulence driven on spatial scales comparable to the radius at which the SNR merges with the ISM, i.e. the shock velocity becomes of the order of velocity dispersion in the ISM. M15 quantified the final momentum (p fin ) driven by an isolated Figure 1. Map of H i column density derived from 21 cm observations by Braun et al. (2009). Overlaid is the footprint of the PHAT survey, with observation bricks outlined as large red rectangles (Dalcanton et al. 2012). Stellar age distributions are measured in 83 cells within the PHAT area (Williams et al. 2017). The shaded red squares show the locations of cells located between deprojected radii of 10-13 kpc from M31's center; they roughly cover the main starforming ring of M31 (see Section 2.4). We will compare our model (Section 2) with the observed velocity dispersion in this ring.
SNR well past its shell formation stage in a turbulent ISM as where p fin /m * is the momentum deposited per mass of stellar population (we set m * = 100 M per M15), Z is the metallicity and n h is the ISM density 1 . We note that this form of p fin is similar to other independent highresolution studies of momentum deposition by SNRs in an inhomogenous ISM (e.g. Kim & Ostriker 2015).
M16 used this subgrid model of momentum feedback to simulate the SN-driven ISM at 2-4 pc spatial resolution, and showed that the resulting velocity dispersion in a steady-state Milky Way-like ISM can be described by an analytical equation where the energy injection rate of SN momentum-driven turbulence is balanced by its corresponding rate of decay. The resulting mass-weighted velocity dispersion (σ p ) is given by the Eq 22. in M15, which we repeat here for convenience, where ρ is the density of gas,ṅ SN is the SN rate per unit volume, and f is a factor that accounts for momentum cancellation when multiple blast waves interact. We set f = 1 for our fiducial runs, then revisit the issue of f in Section 4. We will use this predicted σ p in different regions of M31's ring, as a function of the measured ρ andṅ SN , for comparison with observations in the subsequent sections.

Measurement of SN rate density (ṅ SN )
We set the SN rate per unit volume,ṅ SN using the detailed SAD maps from the PHAT survey and the known properties of SN DTD. We use the Williams et al. (2017) SAD map of the northern third of the disk of M31, spanning a total area of about 0.8 deg 2 (Figure 1). An SAD is measured in each of 826 spatial cells, each 83 wide. Each cell contains the stellar mass formed per look-back time (M ij , where i is the cell and j is the age bin), esti-mated by comparing resolved color-magnitude diagrams of the stars in the region with stellar isochrone models.
We then convert these SAD maps into maps of the SN rate per cell using observationally constrained DTDs. The DTD is defined as the SN rate versus time elapsed after a hypothetical brief burst of star formation. When convolved with the SAD maps described in the previous section, the DTD provides the current SN rate in each region of the galaxy Maoz et al. 2014) in the following way: where M ij is the stellar mass formed in cell i in the age-interval j given by the SAD map, and Ψ j is the DTD value in the age bin j. We use the form of the corecollapse DTD given Eq A.2 of Zapartas et al. (2017), which accounts for the effects of binary stellar interactions at Z . For Type Ia SNe, we assume the parametric form of the Type Ia DTD from Maoz et al (2012) based on the compilation of all observational constraints to date, where t is the delay-time between star-formation and SN. The form of the Type Ia and core-collapse SN DTDs are shown in Figure 2. The volumetric SN rate in cell i (ṅ SN,i for use in Eq. (2)) can then be estimated from R i asṅ where A i is the cell size of each SAD region (≈ 83 × 83 or 310 × 310 pc 2 ) and z sn is scale height of the vertical distribution of SNe. The factor cos(i) accounts for the extended line of sight through the disk as a result of the galaxy's inclination angle i = 77 • (Corbelli et al. 2010), so z sn → z sn /cos(i). We assume z sn = 150 pc for core-collapse SNe and z sn = 600 pc for Type Ia SNe, as explained in Section 2.5.

Measurements of ISM density and velocity dispersion
Most of the ISM mass in star-forming regions is in the atomic (H i) and molecular phases, so we use maps of the 21 cm line of H i (Braun et al. 2009) and the 115 GHz line 12 CO(J=1-0) (Nieten et al. 2006) in M31. The data cubes of Braun et al. (2009) were obtained using the Westerbork Synthesis Radio Telescope (WSRT) and the Green Bank Telescope (GBT), with a spatial resolution of 30 (or 113 pc at the distance of M31). The H i column density (N HI ) and non-thermal velocity dispersion (σ HI ) were measured by Braun et al. (2009) from the 21 cm emission along each line of sight assuming a model of an isothermal, turbulence-broadened line profile. We note that evidence for opacity-corrected H i features in 21 cm is somewhat inconclusive in more recent observations in M31 and M33 (Koch et al. 2018(Koch et al. , 2021, so we use the opacity-uncorrected map of Braun et al. (2009) (their Fig. 15). The difference in the predicted velocity dispersion from the two different version of the density maps is about ≈ 12%, which does not affect our conclusions later. Molecular hydrogen column densities were obtained from the 12 CO(J=1-0) emission map of Nieten et al. (2006) using the single-dish IRAM 30m telescope. The survey covered 2 • × 0.5 • of the M31 disk, yielding a map of CO-line intensity at a final angular resolution of 23 (spatial resolution ≈ 87 pc at the distance of M31). The CO-line intensities were converted into H 2 column densities (N H2 ) using the conversion factor X CO = 1.9 × 10 20 mol cm −2 (K km s −1 ) −1 assumed by Nieten et al. (2006). The total mass of H 2 in M31 is about 14% that of H i in the M31 ring. For convenience, we use the H i velocity dispersion as a proxy for H 2 velocity dispersion using the radius-independent ratio of σ HI /σ H2 = 1.4 measured by the HERACLES CO and THINGS H i surveys of nearby galaxies (Mogotsi et al. 2016), as well as in M33 (Koch et al. 2019). We note here that our inclusion of H 2 measurements is done to account for the velocity dispersion of the 'mass-weighted' ISM, in order to be consistent with the M15 and M16 models, which also predicts the mass-weighted turbulent velocity dispersion.
We combine the density and velocity dispersion of the atomic and molecular phases into an effective massweighted ISM. The total mass-weighted non-thermal velocity dispersion in the H i and molecular phases is We assume the vertical distribution of ISM in M31 is centered on the disk midplane, and approximately Gaussian for the molecular phase and exponential for H i, consistent with observations of our Galaxy (Dickey & Lockman 1990;Ferrière 2001). Each phase is characterized by an 'effective' scale height, which we discuss further in Section 2.5.
For each SAD cell with H i column density N HI , the H i density along the line of sight z is, As in Eq. 5, the scale height has been corrected for the inclination of M31 with the factor cos(i). Similarly for each SAD cell with H 2 column density N H2 , the corresponding H 2 volume density is, For ease of interpretation (given our simplified ISM model), we will compare the observed velocity dispersions with the minimum velocity dispersion predicted by models. We enforce this by assuming all SNe explode at the mid-plane density, i.e. n h = n h (z = 0). This is approximately a lower-limit on σ p (we call this σ min p ) per SAD cell since SNe exploding away from the midplane but still within the scale height of gas would interact with lower densities than at the midplane, deposit greater momenta (Eq. 1), and contribute to an effectively higher σ p per SAD cell. This lower limit is also a good assumption since we neglect all other sources of stellar feedback (e.g. winds, cosmic rays) that could add to the momentum budget per SAD cell depending on the environment. Comparing this minimum feedback from SNe with observations leads to some interesting insight as we show in Section 4. For each value of (N HI , N H2 ), we derive (n HI , n H2 ) using Eq (6) and (7), convert to a total hydrogen mass density ρ = m p (n HI + 2n H2 )/X H in units of g/cm 3 (where m p = 1.67 × 10 −24 g and X H = 0.76 is the mass fraction of hydrogen), and feed it into Eq. (2). We also take the total number density of hydrogen, in units of atoms cm −3 as n h = n HI + 2n H2 , for use in Eq 1.

Galactocentric radii of SAD cells
We first calculate the deprojected distance of each cell from the center of M31 using the method in Hakobyan et al. (2009). Let (α, δ) be the sky-projected location of each cell centroid, and (∆α, ∆δ) be the skyprojected angular offset from M31 center (located at α M 31 = 00 h 42 m 44.3 s , δ M 31 = +41 • 16 9 ) 2 . Assuming a position angle of M31's disk, θ p = 38 • (Corbelli et al. 2010), the location (u, v) of each SAD cell in M31's coordinate system is u = ∆α sinθ p + ∆δ cosθ p v = ∆α cosθ p − ∆δ sinθ p The radial distance of each cell in the plane of M31 from the M31 center, corrected for M31's inclination (i = 77 • ; Corbelli et al. 2010), is, where d is the angular distance from the center in arcseconds. We identify the "ring" as SAD cells with 10-13 kpc, as shown by the shaded region in Figure 1.

Assumptions about SN and ISM scale heights
In this section, we describe plausible ranges and fiducial values for our free parameters: the atomic and molecular scale heights (z HI and z H2 ) and SN scale heights z sn (here on, we will specify the separate scale heights of SNe Ia and CC as z Ia and z cc respectively).
Core-collapse SNe generally occur at lower effective scale heights than SNe Ia (Hakobyan et al. 2017).
In the Milky Way, open clusters younger than 100 Myrs are all situated within 200 pc of the midplane, with an effective scale height of 60-80 pc (Joshi et al. 2016;Soubiran et al. 2018). Since their age and velocity distribution nicely follows that of field stars (Baumgardt et al. 2013;Soubiran et al. 2018), we can assume that the general population of core-collapse SN progenitors in the Milky Way also has a scale height of 60-80 pc. However, the disk of M31 is kinematically hotter and more extended than the Milky Way (Ivezić et al. 2008;Collins et al. 2011). Based on the ratio of scale heights to scale lengths observed in edge-on disk galaxies (Yoachim & Dalcanton 2006, Collins et al. (2011) proposed that the M31 disk could be 2-3 times thicker than the Milky Way (although this may be an over-estimate as the galaxies in the Yoachim & Dalcanton 2006 sample are different and less massive than M31). We therefore assume that in M31, 60 pc < z cc < 200 pc is a plausible range for the scale height of core-collapse SNe.
Older (∼Gyr) stars are mostly concentrated in the thin and thick disk with measured scale heights in the range of 140-300 pc and 500-1100 pc respectively in the Milky Way (Li et al. 2018;Mateu & Vivas 2018). The thin disk is slightly younger, with ages in the range of 7−9 Gyrs compared to the thick disk's age of ∼ 10 Gyrs (Kilic et al. 2017). The measured shape of the Type Ia SN DTD suggests that progenitors younger than 10 Gyrs will produce the majority of Type Ia SNe (Maoz & Badenes 2010), so we assume Type Ia SN progenitors are roughly distributed at the same scale height as the thin disk, ∼300 pc. This is also consistent with the scale height of SDSS white dwarfs (De Gennaro et al. 2008;Kepler et al. 2017;Gentile Fusillo et al. 2019) and about 4 times the scale height of young core-collapse progenitors, so for simplicity we assume that in M31, z Ia ≈ 4z cc , and z cc is in the range mentioned previously.
H i scale heights in M31 were measured by Braun (1991) in the range of z HI = 275 − 470 pc between radii of 10-13 kpc in M31. We are not aware of any scale height measurements of the molecular phase in M31, but the Milky Way can provide some supplementary information. Studies of the H 2 profiles traced by CO in the Milky Way have measured a half-width at half-maximum scale height of 50-80 pc (consistent with being a bit smaller than z cc ), which is about a factor of 3 lower than the scale height of H i in the Milky Way (Marasco et al. 2017).
Given these constraints, we can assume that z cc is always less than z Ia , z H2 is always less than z HI , z Ia 4z cc and z H2 ≈ z HI /3 . Given the range of values allowed by observations, we first analyze our results for a fiducial model where z cc = 150 pc, and z HI = 350 pc, giving z Ia = 600 pc, and z H2 = 117 pc. We then change the values of these parameters and their ratios within the plausible ranges discussed previously to assess the impact of assumptions in Section 3.2.

RESULTS
In this section, we show the distribution of SN rates across the M31 ring as measured from our SAD map and DTDs, and a comparison of our predicted velocity dispersions predicted by these rates with the observed values along the ring.

Distribution of SN Rates
Figure 3(a) shows our SN rate distribution estimated from the DTDs and SADs as described in in Section 2.2 (Eq 3). The integrated SN rate in the region we identify as M31's 10 kpc ring is 1.74 × 10 −3 SN yr −1 , with roughly 39% contribution from SN Ia and 61% from core-collapse SNe.
The fraction of this SN rate of Type Ia versus corecollapse is shown in Figure 3(b). About 75% of the ring has a higher core-collapse rate than Type Ia. These regions coincide well with young star-forming regions identified in UV and IR images of M31 (Lewis et al. 2017), and are mostly concentrated in the inner parts of the ring. Regions with the highest core-collapse rates, exceeding that of Type Ia by more than a factor of 3, coincide with the well-known star-forming region OB54, with nearly 3.8×10 5 M of stars younger than 300 Myrs (Johnson et al. 2016, also see Figure 7 in this paper). SNe Ia generally dominate the total SN rate near the edges of our ring region, coinciding with the inter-arm region as seen in Figure 3(b), and exceeding the corecollapse rate by up to a factor of 3 in some SAD cells.  As evidence of the high characteristic SN rate of the ring, we also show the distribution of optically-selected SNRs in M31 by Lee & Lee (2014) in Figure 3(b). The majority of SNRs are concentrated along the M31 ring, and particularly associated with regions of higher corecollapse fraction. A more quantitative test of whether the observed SNR distribution is consistent with the SN rates will be the subject of a future paper, since it requires a more rigorous analysis of the poorly understood completeness of SNR catalogs (particularly at optical wavelengths).

Comparison of model and observed velocity dispersion
We compare the observed (σ obs ) versus mininum predicted velocity dispersion (σ min p ) in the mass-weighted ISM in Figure 4. The observed velocity dispersion exhibits a range of values spanning 4-12 km/s, whereas the predicted values extend up to 20 km/s or higher. On average, we find that for our fiducial model described in Section 2.5, the σ min p mostly exceed the observed values σ obs , but within a factor of two for 84% of the SAD pixels in the ring. To understand why our velocity dispersion model over-estimates the observed values in  the column density and SN rate, the two fundamental parameters in our model, in Figure 5. We find hint of a negative correlation in σ min p /σ obs with N H and a positive correlation with SN rate. In particular, SAD pixels with log (N tot /cm −2 ) < 21.3 and SN rate or Log (R i /yr −1 ) > −4.7 mostly have σ min p /σ obs > 1. We examine this more closely in Section 4.
We briefly discuss the impact of model uncertainties which are certain to alter the σ min p /σ obs measurements. The SN rates can vary by an average of 15% (maximum of 50%), depending on the isochrone model used for constructing SADs (Williams et al. 2017), but this has a relatively small impact on our result. For example, using the MIST SAD solutions (which we have been using), we have about 82% SAD pixels with σ min p /σ obs > 1, whereas using the PARSEC SAD solutions results in 74% of SAD cells having σ min p /σ obs > 1, and the correlations with density and SN rates remain. Assumptions about the scale height of gas and SNe directly affect the midplane densities and SN rate densities, which has a larger effect on σ min p /σ obs measurements. We therefore assess the impact of varying scale heights on the σ min p /σ obs values as as shown in Figure 6. For smaller gas scale heights and larger core-collapse scale heights, σ min p /σ obs decreases. This is because smaller gas scale heights imply a higher volume density of ISM for a given column density (Eq. 6, 7), which reduces the momentum deposition and turbulence driving based on Eq. 1. For larger SN scale heights, the SN rate per unit volume is smaller (Eq. 5), which likewise reduces the momentum deposition rate (Eq. 1). Within the plausible range of scale heights discussed in Section 2.5 and marked as a box in Figure 6, about 60% of SAD pixels still have overpredicted velocity dispersion. The part of the parameter space in Figure 6 where the fraction of over-predicted cells are below 10 − 20% involves z cc > z HI , which is unlikely given the close association of core-collapse SNe with gas-rich star-forming regions.

Insight Into Momentum feedback efficiency
Our analysis has shown that simple models of ISM turbulence driven by isolated non-overlapping SNRs are consistent within a factor of 2 of observations for most of the star-forming/ISM environment of the M31 ring. Some of the discrepancy can be explained by variation in model parameters (e.g. scale heights of stars and gas) as explained in Section 3.2, but in this discussion, we particularly focus on cells with Log (N tot /cm −2 ) < 21.3 and Log (R i /yr −1 ) > −4.7, since all the cells in this range have σ min p /σ obs 1 even after accounting for plausible variation in SN rates and scale heights in Section 3.2.
Regions with σ min p /σ obs > 1 values are interesting in the sense that there are multiple sources of stellar feedback such as stellar winds, radiation pressure, photoionization and cosmic rays, but here the hydrodynamical momentum from SN blast-waves alone over-predict the observed ISM turbulence.
One reason behind these over-predicted cells could be that σ obs were underestimated in our maps, but this is unlikely. More recent, sensitive VLA-based H i surveys (e.g. Koch et al. 2018Koch et al. , 2021 suggest that a clean separation of thermal and non-thermal components of the 21 cm line is non-trivial. It is likely that the assumption of Braun et al. (2009) of an isothermal H i component along the line of sight results in some residual thermal contribution to the non-thermal velocity dispersion. Thus the non-thermal velocity dispersion in M31 we are using from Braun et al. (2009) may be an upper-limit to the actual turbulence contribution.
The regions with σ min p /σ obs > 1, especially in the cells with Log (N tot /cm −2 ) < 21.3 and Log (R i /yr −1 ) > −4.7, may therefore indicate a drawback in the SN momentum-driven turbulence model, so we investigate these regions visually in Figure 7. From here on, we couch the column density and SN rate cutoff in terms of a volume density and SN rate surface density cutoff, to be consistent with values used in simulations. The low column density cutoff corresponds to n h < 0.2 cm −3 for our fiducial scale heights, and the SN rate cutoff corresponds to a SN rate surface density Σ SN > 2.1 × 10 −4 SN yr −1 kpc −2 .
The low-density cells are situated at the upper and lower edges of the star-forming ring as shown in Figure 7. Comparison with Figure 3(b) shows that these regions also have a higher rate of SN Ia than CC, with an average SN Ia/CC ratio ≈ 1.67 in these cells. M16 noted that at low densities, SNRs have longer cooling timescales, and may come into pressure equilibrium with the ISM before cooling or depositing a significant amount of momentum (McKee & Ostriker 1977). This missing physics in our model is likely the reason for the over-predicted σ min p . SAD cells with Log (R i /yr −1 ) > −4.7 are spatially correlated with the prominent star-forming region OB54 as mentioned in Section 3.1. One possibility, as raised by M15 and M16, is that in high star-forming regions, overlapping shocks from close-proximity SNRs might cancel some of the outgoing momentum (parameterized by f in Eq 2, which we had set to 1). For example, a reduction of more than a factor of 2 in σ min p is achieved with f < 0.2, which is consistent, though slightly less than ) and to that measured from 21 cm + CO-line observations (σ obs ) in each SAD cell in the M31 ring (see Figure 1) is compared with the observed column density (atomic + molecular; panel a) and the estimated SN rate (panel b) in those cells. Both panels provide the same information, but NH and SN Rate switches between the x-axis and the colorbar. The points represent the lower limit to the ratio as it corresponds to SNe exploding in M31's midplane (see Section 2.5).  by clustered SNe driving a hot, over-pressurized outflow (Sharma et al. 2014;Gentry et al. 2017). This explanation is plausible given the detection of X-ray emission in this region by Kavanagh et al. (2020), and was also given as an explanation by previous energy-balance studies that similarly observed an excess of SN energy over the measured ISM turbulent energy in the central high star-forming regions of galaxies (Tamburro et al. 2009;Stilp et al. 2013;Koch et al. 2018;Utomo et al. 2019). A remaining possibility is that a non-negligible ISM mass in these regions is sustained at warmer phases, invisible to 21 cm or CO-line maps, due to the cumulative heating by SNe and pre-SN processes like winds and photoionization.
The results above indicate that fiducial models of momentum feedback from SNe used by most cosmological simulations, which generally assume non-overlapping, non-clustered SNRs, may require adjustment at low densities and at high SN rates due to aforementioned non-linear effects of clustering and SNR evolution at low-densities. This can be quantified by a suppression factor f (n h , Σ SN ) < 1 for n h 0.2 cm −3 and Σ SN > 2.1 × 10 −4 SN yr −1 kpc −2 , although a more precise form of this relation will be explored in a subsequent paper where we account for the energy and momentum carried away by any high-velocity outflows or warm diffuse gas from these regions.
The result also highlights the role of Type Ia SN feedback in low-density regions of the ISM, where it can exceed core-collapse SN rates (Li et al. 2020a,b). The energetics of Type Ia SNe are particularly pronounced in the central few kpc of M31 (though not explored in this paper), where it is likely responsible for the bright X-ray halo emission and depleted metal abundances in the region (Tang et al. 2009;Telford et al. 2019).

Comparison with previous studies of ISM energy balance
As the molecular ISM in our data is only ∼ 14% of the atomic ISM, our results primarily explore feedback in the atomic ISM, and therefore it is interesting to compare our work with previous studies of energy balance in the ISM traced by atomic hydrogen. Tamburro et al. (2009) showed that SN energy alone can drive turbulence in atomic gas within the optical radius of nearby galaxies, with an approximate coupling efficiency of ∼ 10% (Thornton et al. 1998;Mac Low & Klessen 2004). Similar results were also obtained by Stilp et al. (2013) with globally-averaged H i observations. More recently, Koch et al. (2018) and Utomo et al. (2019) extended these techniques to M33, with the latter study allowing coupling efficiency to vary with radius.
A key difference between our work and previous ones is that we examine spatially-resolved H i line profiles along different lines of sight, as opposed to globally or radially-stacked H i profiles. This allows us to compare the observed ISM turbulence with the local properties of the star-forming and ISM environment.
Phenomenologically, there are a few key differences between our work and previous studies also worth mentioning. Similar to Utomo et al. (2019), we account for turbulence driven by momentum-conserving phase of isolated SNRs, and constrain the efficiency of this momentum-feedback driving the observed turbulence, as opposed to previous studies that considered the efficiency of initial SN energy (= 10 51 ergs) going into the ISM turbulence (the majority of which will be radiated away without impacting the gas). The M15 and M16 models also assume that turbulence is driven at the radius where SNR dissolves into the ISM, which strongly depends on the ISM density (i.e. when v s ∼ σ). This is different from the assumption of constant driving scale (equal to the scale height) or decay timescale in Tamburro et al. (2009) and Koch et al. (2018). These assumptions affect the predicted SN feedback. For example, Utomo et al. (2019) showed that a spatiallyvarying decay timescale allowed SNe to drive turbulence in M33 out to 7 kpc instead of 4 kpc, by which point the star-formation rate and gas densities decrease by an order of magnitude compared to the central region. Bacchini et al. (2020) similarly showed that a variable decay timescale makes SNe efficient enough to drive turbulence in the THINGS galaxies throughout, as opposed to just within the optical radius (Tamburro et al. 2009).
Despite the differences in methodology, our work agrees with previous studies that SN energy driving is inefficient, particularly in regions characterized by high star-formation rates. The higher spatial resolution offered by a more nearby galaxy like M31 reveals that regions where our models disagree with observations also correlate with regions of clustered star-formation, signifying the importance of taking into account clustering effects in SN feedback models. A direct comparison with previous studies is complicated given the differences in methodology, but our work highlights the importance of spatially-resolved observations in Local Group galaxies in the study of SN feedback.

CONCLUSION
In this paper, we have tested the paradigm of SN momentum-driven ISM turbulence developed by recent high-resolution vertical disk simulations. We compare model prescriptions with resolved observations of stellar populations and ISM in the prominent 10-kpc star-forming of M31, where stellar feedback is expected to be the main source of turbulence. The spatially-resolved PHAT stellar photometry in the northern third of the disk provides detailed stellar-age distributions (SADs) in ≈ 310 pc 2 cells, which we convolved with known forms of the SN delay-time distribution to predict the corecollapse and Type Ia rates across M31's ring. We used ISM densities of the neutral atomic gas (traced by 21 cm H i line maps) and molecular gas (traced by 12 CO(1-0) line maps) alongside the SN rates to predict the steadystate mass-weighted turbulent velocity dispersion, using the feedback prescriptions of Martizzi et al. (2015) and Martizzi et al. (2016). We compared these model estimates against the scaled turbulent velocity dispersion obtained from H i and CO maps of M31. We assumed all SNe explode in the galaxy midplane where the lineof-sight density is highest, effectively providing a lowerlimit on the predicted velocity dispersion. We summarize the following key results from our work :-1. We find an integrated rate of ≈ 1.7 × 10 −3 SN yr −1 in the ring covered by PHAT, with 61% contribution from core-collapse SNe. Regions with dominant core-collapse contribution coincide with known star-forming regions as expected, while regions with dominant Type Ia contribution fall near the edges of the ring.
2. We found that the minimum predicted velocity dispersion exceed observed values in 84% of the ring covered by PHAT for our fiducial model within a factor of 2. Some of the discrepancy can be explained by varying the assumptions regarding SADs and ISM/SN scale heights within plausible limits, but for densities 0.2 cm −3 and SN rates > 2.1 × 10 −4 SN yr −1 kpc −2 , the discrepancy appears to increase.
3. SAD cells with SN rates > 2.1 × 10 −4 SN yr −1 kpc −2 where velocity dispersion is over-predicted are spatially correlated with dense concentration of young clusters embedded in a bright thermal X-ray region. This supports the possibility of clustering of SNe in this regime, which is not captured in our momentum feedback model. Clustering of SNe can lower the momentum deposited per SN and mass-weighted turbulence in the ISM as a result of converging shocks from adjacent explosions, mass-loaded outflows, or higher mass fraction in warmer ISM phases due to cumulative action of stellar winds and SNe.
4. The low density ( 0.2 cm −3 ) regions where velocity dispersion is over-predicted coincide with the edges of our ring region where Type Ia SNe dominate the injection rate by nearly a factor of 2. However given the overall low SN rate in these regions, it is likely that the discrepancy could be due to isolated SNRs coming into pressure equilibrium with the ISM before significant amount of cooling and momentum deposition takes place-another effect not included in our models.
Our results provide observational support for including adjustments in fiducial subgrid models of momentum feedback, to account for SNR evolution in clustered and in low-density environments. The work underscores the importance of resolved stellar photometry and cloudscale atomic and molecular ISM observations for assessing feedback models and ISM turbulence. Newer, more sensitive observations at high spectral resolution, such as the VLA maps of Koch et al. (2021) can provide more detailed characterization of turbulence in atomic clouds in M31. Preliminary comparison have shown that the 2nd moment line-widths of the VLA maps are within a factor of 2 of Braun et al. (2009) non-thermal values, although the former do not yet cover the full PHAT area or the M31 ring. We will expand our present analysis in future papers with data from the ongoing Local Group L-Band Survey 3 which will cover all of M31, as well as M33 and four Local Group dwarfs, providing H i maps of unprecedented sensitivity at a wide range of spatial resolution. This extension would allow us to test feedback models with spatially resolved observations across a wide range of star-forming conditions, and empirically obtain corrections to the fiducial models to be included in cosmological simulations. ACKNOWLEDGMENTS S.K.S is grateful to Robert Braun for sharing the WSRT 21 cm maps of H i density and non-thermal velocity dispersion in M31. SKS and LC are grateful for support from NSF grants AST-1412549, AST-1412980 and AST-1907790. Parts of this research were supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. We acknowledge support from the Packard Foundation. E.R-R is supported by the Heising-Simons Foundation and the Danish National Research Foundation (DNRF132).