Deep Narrowband Photometry of the M101 Group: Strong-line Abundances of 720 H ii Regions

We present deep, narrowband imaging of the nearby spiral galaxy M101 and its satellites to analyze the oxygen abundances of their H ii regions. Using Case Western Reserve University’s Burrell Schmidt telescope, we add to the narrowband data set of the M101 Group, consisting of Hα, Hβ, and [O iii] emission lines and the blue [O ii] λ3727 emission line for the first time. This allows for complete spatial coverage of the oxygen abundance of the entire M101 Group. We used the strong-line ratio R 23 to estimate oxygen abundances for the H ii regions in our sample, utilizing three different calibration techniques to provide a baseline estimate of the oxygen abundances. This results in ∼650 H ii regions for M101, 10 H ii regions for NGC 5477, and ∼60 H ii regions for NGC 5474, the largest sample for this Group to date. M101 shows a strong abundance gradient, while the satellite galaxies present little or no gradient. There is some evidence for a flattening of the gradient in M101 beyond R ∼ 14 kpc. Additionally, M101 shows signs of azimuthal abundance variations to the west and southwest. The radial and azimuthal abundance variations in M101 are likely explained by an interaction it had with its most massive satellite, NGC 5474, ∼300 Myr ago combined with internal dynamical effects such as corotation.


INTRODUCTION
Star forming regions provide a tracer of galactic chemical evolution through the gas-phase metallicity, typically measured in H II regions by the ratio of oxygen to hydrogen atoms.This oxygen is synthesized in highmass stars and released to the interstellar medium (ISM) through stellar winds or supernovae.It is then mixed throughout a galaxy through a variety of mixing processes, such as mergers, turbulence, and instabilities in the disk (see Roy & Kunth 1995).
In the past, most studies of the oxygen abundances1 of H II regions have revealed the presence of monotonically decreasing radial gradients, typically around −0.05 dex kpc −1 (e.g.Vila-Costas & Edmunds 1992;Zaritsky et al. 1994;van Zee et al. 1998), i.e., the inner regions of galaxies are more metal-rich than the outskirts.These gradients were easily explained as a consequence of the inside-out galaxy formation scenario (e.g.Scannapieco et al. 2009).However, the presence of breaks in the oxygen abundance gradients have been observed in numerous galaxies, including M101 (Vila-Costas & Edmunds 1992;Zaritsky 1992).
Traditionally, oxygen abundance breaks, beyond which the abundance gradient of a galaxy flattens, have only been detected by using the so-called "strong-line" methods, those methods that relate the oxygen abundance to the strong recombination and collisionallyexcited (forbidden) lines (e.g., [O II]λ3727, Hβ, [O III]λλ4959, 5007, Hα).These methods are used when auroral lines, like [O III]λ4363, are undetectable and don't allow for the electron temperature of a gas to be estimated, called the "direct" method (Dinerstein 1990;Skillman 1998;Stasińska 2007).Although the strongline methods have gained widespread use with their easily measurable lines, they come with their own host of issues (see, e.g., Pilyugin 2003;Kewley & Ellison 2008, for a discussion).It was these defects that were blamed for producing "spurious" breaks where none existed.
However, with the creation of newer strong-line abundance methods, breaks in oxygen abundance gradients have been reintroduced to the literature.Initially, only galaxies with extended-UV (XUV) emission (Thilker et al. 2007) were found to have flattened gradients, such as M83 (Bresolin et al. 2009) or NGC 4625 (Goddard et al. 2011).While these results were produced using strong-line abundances, the introduction of integral field spectroscopy allowed direct methods to be applied across entire galaxy disks showing that many non-XUV galaxies have flattened radial abundances beyond their R 25 (Sánchez et al. 2014;Sánchez-Menguiano et al. 2016).Even our own Galaxy has some evidence of abundance flattening, measured not from H II regions but from the [Fe/H] ratio of open clusters (Lépine et al. 2011) and the [O/H] ratio of OB stars (Daflon & Cunha 2004).
Numerous physical mechanisms have been proposed to explain the flattening of abundance gradients at large radii.Some attribute it to dynamical effects, such as the barrier effect of corotation, which isolates the inner and outer regions of a galaxy disk from each other (Lépine et al. 2011;Scarano & Lépine 2013), or other resonances with spiral density waves (Sellwood & Binney 2002).Others attribute it to an environmental cause such as minor mergers or satellite interactions increasing the gas content in the outer disk (Bird et al. 2012;López-Sánchez et al. 2015).There is also the possibility of a different star formation efficiency in the outer disk (Bigiel et al. 2010;Bresolin et al. 2012) as would be predicted in a nonlinear Kennicutt-Schmidt law (Esteban et al. 2013).Whatever the cause, it is clear that the chemical evolution of a galaxy is tied to its dynamical evolution.
In addition to breaks in radial abundance gradients, there are also some suggestions that there should be azimuthal variations across galaxy disks as well.This is predicted to be caused by the different scales on which chemical enrichment operates (Roy & Kunth 1995).Given the large difference in timescales necessary to produce oxygen (OB star lifetime; <10 Myr) and mixing on sub-kiloparsec scales (10-100 Myr; Roy & Kunth 1995), there should be some azimuthal inhomogeneity of oxygen abundance.Tentative evidence has been found in the Milky Way via the [Fe/H] ratio in Cepheids (Pedicelli et al. 2009;Lépine et al. 2011;Genovali et al. 2014), as well as in M101 via oxygen abundances (Kennicutt & Garnett 1996;Li et al. 2013).Despite this, detecting azimuthal variations in other galaxies has been difficult because spectroscopic studies are usually limited to bright H II regions resulting in biased coverage.
In an effort to add to the growing evidence of radial breaks and azimuthal variations in the oxygen abundances of galaxies, we have used Case Western Reserve University's (CWRU's) Burrell Schmidt 24/36-inch telescope to add to the deep, wide-field, multiline, narrowband observations of the nearby spiral galaxy M101 (NGC 5457) and its group environment presented in previous papers (Watkins et al. 2017;Garner et al. 2021).In addition to the previous narrowband images targeting Hα, Hβ, and [O III] (Watkins et al. 2017;Garner et al. 2021), we now add deep narrowband images targeting Note-Hubble types are taken from the CVRHS system of Buta et al. (2015) and the sizes of NGC 5477 and NGC 5474 are taken from the RC3 (de Vaucouleurs et al. 1991), while the size of M101 is taken from Mihos et al. (2013).The physical distance to M101 is adopted from Matheson et al. (2012), and references therein.Rproj is the projected distance from M101.
[O II] for the first time.This allows us to sample the entire disks of M101 and its major satellites, NGC 5477 and NGC 5474 (see Table 1), with our large survey area (∼6 deg 2 ) and our photometric depth which allows us to analyze faint inter-arm and outer disk H II regions.M101 was chosen for this survey because its nearby distance (D = 6.9 Mpc; see Matheson et al. 2012 and references therein) enables its properties to be studied in detail.Notably, given its close distance and subsequently high spatial resolution, even the faintest lowluminosity H II regions can be resolved.M101 is also currently interacting with its satellite population, likely the massive satellite NGC 5474 (Beale & Davies 1969;Rownd et al. 1994;Waller et al. 1997;Linden & Mihos 2022).Given that interacting systems have systematically lower oxygen abundances than isolated galaxies and sometimes have flattened gradients (Kewley et al. 2006), and that M101 has been heavily studied by spectroscopic surveys (e.g.Kennicutt & Garnett 1996;Kennicutt et al. 2003;Li et al. 2013;Croxall et al. 2016;Hu et al. 2018), the M101 Group is an excellent group for studying the chemical evolution of interacting galaxies.

OBSERVATIONS AND DATA REDUCTION
The narrowband imaging data for this project was taken using CWRU Astronomy's 24/36-inch Burrell Schmidt telescope located at Kitt Peak Observatory.Full details of our narrowband imaging techniques are given in Watkins et al. (2017), and summarized briefly here.Quantitative information about the narrowband filters and final image stacks for our imaging dataset is given in Table 2.
To measure adjacent stellar continuum for each line, we also observed the galaxy in narrow off-band filters shifted in wavelength by ≈100-150 A from each on-band filter.The Burrell Schmidt images a 1. °65 × 1. °65 field of view onto a 4096 × 4096 back-illuminated CCD, yielding a pixel scale of 1.45 arcsec pixel −1 .In each filter, we took 50-70 1200 s images of the galaxy, randomly dithering the telescope by 10-30 arcmin between exposures.We observed on photometric nights with no moon, yielding night sky levels of 50-100 ADU pixel −1 .We built flat fields in each filter from a combination of twilight flats and stacked night sky exposures (see Watkins et al. 2017), and constructed a spatially dependent scattered light model for the degree-scale point spread function (PSF) of bright stars using deep 1200 s observations of the bright stars Regulus and Arcturus in each filter.Finally, to assist in flux calibration, we also observed spectrophotometric standards (Massey et al. 1988) throughout the course of each observing run.
During data reduction, we first subtracted from each image a median-combined nightly master bias frame, then flat-fielded the images using the composite twilight/dark-sky flat field.We then removed scattered light from bright stars in the field using our degreescale PSF model (following the technique of Slater et al. 2009), after which we model and subtract the sky background using a plane fit to blank sky pixels in each image.We then register and median-stack the images to create final on-and off-band composite images in each filter, with total exposure times varying from 18-24 h each.
We flux calibrate the final image stacks using a variety of techniques, intercomparing the results to estimate our final flux zeropoint uncertainties.We calibrate by (1) deriving a photometric solution from observations of Massey et al. (1988) spectrophotometric standard stars, (2) measuring zeropoints from ugr magnitudes of the ∼150 stars in the M101 field from the Sloan Digital Sky Survey (SDSS; and including a synthesized broadbandnarrowband color term for each filter), and (3) synthesizing narrowband magnitudes using SDSS spectroscopy of ∼100 point sources in the M101 field.These different techniques yielded flux zeropoints that agreed to within ±5 %, which we take to be our absolute flux uncertainties in the imaging data.
Table 2 also details the 1σ limiting depth of the images, measured in two ways.First, the per-pixel noise level (σ pix ) in each image is measured using the standard deviation of pixel intensities measured within regions of blank sky around M101.Second, the limiting depth over larger scales is measured from the scatter in mean flux measured with 30 × 30 (20 × 20 pixel) boxes around M101 (σ 30 ), and traces our sensitivity to diffuse extended narrowband emission.
We stress that as noted in Table 1, throughout this paper we use a different value for M101's R 25 = 8 (Mihos et al. 2013) than the usual RC3 value of R 25 = 14. 4 (de Vaucouleurs et al. 1991).As noted by Mihos et al. (2013), a small patch of starlight to the north of M101's main disk and at the end of the galaxy's northeastern arm is likely to blame for the large radius measured by de Vaucouleurs et al. (1991) despite most of the galaxy not extending that far.We use the areal-weighted R 25 measured by Mihos et al. (2013) as a more robust estimate of the galaxy's size.When comparing our data to spectroscopic data, we will convert their scaled radii to our scale.This has the effect of roughly doubling the scaled radii and halving the abundance gradients reported in the original papers.

H II Region Detection
Our final imaging dataset from the Burrell Schmidt consists of the narrowband imaging described in Section 2, as well as deep broadband imaging in (modified) Johnson B and Washington M filters from Mihos et al. (2013).The broadband photometry has been transformed to standard Johnson B and V as described in Mihos et al. (2013).Our general technique for detecting H II regions in M101 and its satellites is much the same as utilized in Garner et al. (2021) with some key differences described below to account for region detection within a galaxy rather than interspersed throughout the M101 Group.
However, as described in detail by Garner et al. (2021) and noted elsewhere (e.g.Kellar et al. 2012;Watkins et al. 2017), stars will broadly mimic Hα emission as our off-band filters can sample absorption features in the stellar continuum.In order to mask the contribution of these stars, we utilized two surveys targeting both bright and faint stars.We masked those stars in the Tycho-2 Catalog (Høg et al. 2000) brighter than B T = 12.5.All of these stars were found outside the disks of the galaxies and do not affect our analysis.For fainter stars not in the Tycho-2 Catalog, we used entries from the Note-ZP (flux) converts 1 ADU to erg s −1 cm −2 in the master images, while ZP (AB) converts to AB magnitudes.σpix is the per-pixel standard deviation measured in erg s −1 cm −2 arcsec −2 , while σ30 is the 1σ variation in mean intensity measured in 30 ×30 blank sky boxes, in the same units.
Two Micron All Sky Survey (2MASS) All-Sky Catalog of Point Sources (Skrutskie et al. 2006), queried through the VizieR interface of astropy's Astroquery package (Ginsburg et al. 2019).Within a 30 × 30 cutout centered on M101, this resulted in 421 point sources being detected.In Garner et al. (2021) we used Gaia data to remove stars from our imaging of the intragroup field around M101, but here we are working within the bright disk of M101 where such data is not available, thus the necessity of using 2MASS photometry instead.Given this list of point sources, we then make photometric cuts to determine whether a particular point source is a star or an H II region.Ideally, we could use the 2MASS JHK s photometry provided in the Point Source Catalog to determine which sources are stars as they should follow a particular path through a J −H vs. H − K s color-color plot (Straižys & Lazauskaitė 2009).Since many of these point sources lie within or very close to M101, the 2MASS photometry has significant uncertainties due to the high and variable background.Since 2MASS does not provide a local background estimate, we instead perform aperture photometry on our B, V , and Hα on-and off-band images using 8. 7 apertures, or six times the image pixel scale.
Figure 1 shows the color-magnitude diagram for the 2MASS point sources with the marker color indicating their Hα equivalent width (EW) as measured from our narrowband imaging.Star-forming regions should be photometrically blue and have high Hα EW (although older H II regions or regions with a strong underlying stellar population can have slightly lower EWs), while stars should have progressively higher Hα EWs as a function of B − V color, with red stars having the highest Hα EW (see Fig. 2 in Garner et al. 2021).Using these photometric arguments, we define stars as having  (B − V > 0.75 OR B ≥ 17) AND Hα EW < 65 A as illustrated by the box outline in Figure 1.
Having masked both bright stars in the Tycho-2 Catalog around M101 and dim stars in the 2MASS Point Source Catalog in front of M101, we create a twodimensional background object to calculate the background sky level and its uncertainty.As in Garner et al. (2021), all galaxies (both background and M101 Group members) were masked and the sky level was estimated in boxes of 100×100 pixels (145 ×145 ) with filter sizes of 10 × 10 pixels (14. 5 × 14. 5).
Then, we detected sources on the continuumsubtracted Hα image using astropy's PhotUtil package and its segmentation module (Bradley et al. 2021).This program detects sources as objects that have some minimum number of connected pixels, each greater than a background threshold value.In this case, the threshold value of 10σ above the background level on the continuum-subtracted Hα image described above after a two-pixel Gaussian smoothing was applied was chosen.This high threshold value was chosen in order to detect bona fide star-forming regions within the galaxies in the M101 Group.As such, we do not mask the galaxies in our images when detecting these regions.
This process resulted in 514 sources in M101 initially.We then deblended close or overlapping sources using 32 multi-thresholding levels and a contrast parameter of 0.015, resulting in a final source list for M101 of 853 sources.We repeat this same process for both of the satellites in our sample, resulting in 6 and 51 sources before deblending and 11 and 71 sources after deblending for NGC 5477 and NGC 5474, respectively.
With the segmentation maps shown in Figure 2 defining the H II regions in each galaxy, we used the PhotUtil segmentation module to calculate photometric and structural quantities for each region in each of the narrowband images as well as in the broadband imaging of Mihos et al. (2013).Many of these were default calculations for segmentation, including positions and fluxes.We also calculated photometric errors, signal-to-noise, and AB magnitudes for each object in each filter.Additionally, we calculated emission line fluxes and EWs for each narrowband imaging pair (Hα, Hβ, [O III], and [O II]), after correcting for a variety of photometric effects as detailed in the next section.

Photometric Corrections
In narrowband photometry, since the off-band filter measures the strength of the continuum near each emission line captured in the on-band filters, theoretically, the flux differences calculated above are equivalent to each emission line's flux.However, there are several physical processes hindering that simple assumption.Figure 3 shows a synthetic, representative spectrum of a star-forming region generated with the Python Code Investigating GALaxy Emission (cigale; Noll et al. 2009;Boquien et al. 2019) with particular emission and absorption lines marked.We plot our narrowband filter curves on top of the spectrum, showing that three of the four sets of narrowband filters have different problems that we need to correct: the Hβ emission line strength is diminished by the effects of Balmer absorption in the stellar continuum; the [O II] off-band filter does not sam-ple a smooth continuum, instead sitting on high-order Balmer absorption; and the Hα on-band filter captures Hα+[N II] emission while the Hα off-band filter captures [S II] emission.Additionally, but not evident in Figure 3, we must make a correction for any internal reddening.In this section, we describe the steps we take to correct for these various processes.
In each H II region, we correct the observed Hβ EWs for underlying stellar absorption using the following technique.Each H II region contains light from the young star-forming population along with an older underlying background stellar population.For the young population, we assume a stellar Hβ EW width of 5 A, based on models and observations of H II regions (e.g.González Delgado et al. 1999;Gavazzi et al. 2004;Moustakas & Kennicutt 2006).For the underlying old stellar population, we adopt the Hβ EW measured in an annulus around each region with a size five to eight times the H II region size.We then calculate a weighted average of these two EW values, weighting by the off-band (i.e., stellar) flux ratio in the H II region and surrounding annulus.This process yields net stellar absorption EW corrections which fall in the range predicted by stellar population synthesis codes (∼2-12 A; González Delgado et al. 1999).In a small number of cases, typically in low signal-to-noise regions, the resulting EW correction fell well outside the accepted ranges, and in these regions we simply assign an Hβ EW correction typical of other well-measured regions at similar radii in the galaxy.This process results in an average Hβ EW correction of 4.0 A, 2.2 A, and 4.1 A for M101, NGC 5477, and NGC 5474, respectively.
Higher order Balmer lines are increasingly sensitive to stellar absorption which gives rise to the absorption features seen in the [O II] off-band filter in the top panel of Figure 3. Here, the [O II] off-band filter sits on Hη absorption and adjacent to Hζ and Hθ absorption.González Delgado et al. (1999) used synthetic spectra to model the behavior of the Hη absorption line as a function of age and star formation history.They found that the Hη EW ranges from 2-10 A with an average of 5 A. However, we cannot make an EW correction as we did to the Hβ EWs due to the complicated spectrum: the [O II] off-band filter is sampling multiple emission and absorption lines and there is a rapidly changing continuum level between the on-and off-band filters.
We use cigale to model a range of star-forming histories for a short burst of star formation as well as those for an older, pre-existing stellar population.For the young population, we find typical on/off flux ratios of 1.00 ± 0.02, while for the older, pre-existing populations we find a lower ratio of 0.88 ± 0.02.For each region, we

M101
NGC 5477 NGC 5474 Finally, we must make a correction for the fact that our Hα on-band filter measures Hα+[N II] emission while the Hα off-band filter measures [S II] emission.The relative [N II]/Hα and [S II]/Hα line ratios vary with metallicity, electron temperature, ionization state, etc., and not necessarily in the same way (Burbidge & Burbidge 1962;Baldwin et al. 1981;Osterbrock 1989).Our previous work (Garner et al. 2021) corrected the Hα emission using a flux ratio of [N II]/Hα = 0.33 (Kennicutt 1983;Jansen et al. 2000).However, James et al. (2005) demonstrated that this ratio works best for [N II]selected H II regions, and overestimates the flux ratios of other regions and galaxies as a whole.
Instead, we utilize published spectroscopic line ratios for H II regions in M101 to determine the relation be- Therefore, we adopt the following Hα correction: where f Hα,true is the true, corrected Hα flux, F Hα,diff is the measured Hα flux difference, i.e.F Hα,on − F Hα,off , and f [S ii] /f Hα ≈ 10 −0.8 (CHAOS).In the above equation, we use the radial gradients in Figure 4 to determine the [N II]/[S II] ratio, where (2) This correction factor has the effect of reducing the observed Hα flux by ∼7-25 % in M101.There is not enough data for the two satellite galaxies to perform a similar analysis.Instead, we employ a single correction   to the data for each satellite.For NGC 5477, we calculate the correction based on the measured line strengths by Berg et al. (2012), while for NGC 5474, where no spectroscopic data is available, we use the mean correction for M101.These corrections amount to only a few percent in each case.
The last correction we make is the empirical extinction correction following the relation found in Calzetti et al. (1994) and assuming the reddening curve found in Calzetti et al. (2000).This correction uses the Balmer decrement, Hα/Hβ, to derive the extinction at Hβ assuming Case B recombination (T = 10 4 K, n e = 10 2 cm −3 ) so the intrinsic Balmer decrement is Hα/Hβ= 2.86 (Osterbrock 1989).Given the measured extinction, applying this correction is straightforward and applied to the line fluxes and magnitudes in each band in each galaxy.
Figure 5 shows radial profiles for each of the four (dereddened) line fluxes (Hα, Hβ, [O III], [O II]) in each of the galaxies in our sample.H II regions in M101 (circles) and NGC 5474 (plus signs) are colored by density of points, while those in NGC 5477 are simply marked by blue stars.The orange lines shows the median values for M101 in bins of 0.1R 25 out to a radius of 2R 25 and the dashed gray lines mark M101's outer disk, which we define at >3 times the azimuthally averaged disk scale length (430 , 14.5 kpc; Mihos et al. 2013;Watkins et al. 2017).Largely, we recover the trends shown in Watkins et al. (2017).Both NGC 5477's and NGC 5474's H II regions span the same range of flux as those in M101's outer disk.Additionally, the weak radial gradients in median Hα and Hβ fluxes of the inner disk, and, to a lesser extent, [O III] and [O II] fluxes, is most likely a demonstration of the so-called Kennicutt-Schmidt law: molecular gas density in M101 declines exponentially with radius (e.g.Kenney et al. 1991).However, the slight flattening we see in the outer disk could be evidence of different star formation efficiency; indeed, Bigiel et al. (2010) showed that the star formation efficiency is considerably flatter in outer disks, between one and two isophotal radii.
Before estimating abundances and presenting the abundance gradients of the three galaxies in our sample, we briefly provide an overview of the numbers and photometric depth of our sample.In M101 (746 regions), we have reached a minimum photometric Hα flux of 1.6 × 10 −16 erg s −1 cm −2 , which corresponds to an Hα luminosity of 9.1 × 10 35 erg s −1 .Using the SFR−L(Hα) calibration of Kennicutt & Evans (2012), this corresponds to an star formation rate (SFR) of 4.9 × 10 −6 M yr −1 .In NGC 5477 (10 regions) and NGC 5474 (65 regions), the minimum photometric Hα flux is ∼10 −15 erg s −1 cm −2 , which corresponds to an Hα luminosity of ∼10 37 erg s −1 and an SFR of ∼10 −5 M yr −1 .We also quantify our lower photometric limits by calculating the total number of ionizing photons, Q 0 , from Osterbrock (1989,

CROSS-COMPARISON AND OXYGEN ABUNDANCES
Having now corrected our data for the effects of underlying stellar absorption and contaminating emission lines, in this section we investigate the emission-line properties of our sample and calculate strong-line abundances.In Section 4.1, we compare our narrowband photometry to that of high-quality spectroscopy per-formed by the CHAOS group.Finding our data matching the spectroscopic data within reasonable limits, in Section 4.2, we investigate the radial trends between three strong-line ratios: the R 23 parameter (Pagel et al. 1979), the ionization state given by O 32 (e.g.McGaugh 1991), and the excitation state given by P (Pilyugin 2000(Pilyugin , 2001)).Finally, in Section 4.4, we derive oxygen abundances from our narrowband photometry using three different methods.

CHAOS Comparison
To gauge the accuracy of our data reduction techniques described above in Section 3.2, we compare our photometric data with the spectroscopic data of the CHAOS group (Croxall et al. 2016).We perform aperture photometry on the regions targeted by CHAOS in all of our narrowband images.Although the spectroscopic data was collected using slits with 1 width, we use 3 photometric apertures in our data, given its lower image resolution (2. 2 FWHM and 1. 45 pixel −1 ).This leads to a systematic bias in the comparison: we expect our fluxes to be brighter with respect to the CHAOS data, and the line ratios may show systematic differences as well due to gradients between neighboring individual H II regions.
With this noted, we perform standard aperture photometry on each of the 77 regions selected by CHAOS. Figure 6 shows the comparison of our photometric data with the spectroscopic data from CHAOS.Reported are each narrowband flux as well as four line ratios, [O III]/Hβ, [O II]/Hβ, R 23 and O 32 (see Section 4.2).The diagonal dashed lines in each plot is the 1 : 1 line.Reported for each comparison is the scatter from the 1 : 1 line, σ.Although there is some scatter in the line fluxes, the data is consistent with a linear trend.However, there is a slight vertical offset from the 1 : 1 line in each of the line fluxes where CHAOS systemically measures lower line emission than in our photometry.This is likely caused by the effect mentioned above where our larger aperture size will measure more flux.

Strong-Line Ratios
As mentioned in the introduction, the most direct and physically motivated method to determine the oxygen abundance of an H II region or star-forming galaxy is to measure the electron temperature, T e , of the ionized gas using the intensity of one or more temperaturesensitive auroral lines (Dinerstein 1990;Skillman 1998;Stasińska 2007).Unfortunately, these auroral lines are intrinsically faint, making them difficult to measure, and none of our filters target these lines.Therefore, we utilize strong-line abundance calibrations to estimate the metallicity of the ionized gas.The caveat is that strongline methods are indirect and usually model dependent.
Modern strong-line calibrations fall into two broad categories: empirical and theoretical.Empirical methods are calibrated against high-quality observations of individual H II regions with measured direct (i.e., T ebased) oxygen abundances (e.g., Pilyugin 2000Pilyugin , 2001;;Pilyugin & Thuan 2005;Peimbert et al. 2007).However, these are limited as they are only strictly applicable to regions similar to those used to make the calibration and are lacking high-excitation H II regions, especially in the metal-rich regime.Theoretical calibrations get around this issue by using photoionization model calcu-lations across a wide range of nebular conditions (e.g., McGaugh 1991;Kewley & Dopita 2002;Kobulnicky & Kewley 2004).These are not without their problems as well, including oversimplified geometries and the unconstrained role that dust plays in the depletion of metals.
The individual limitations of these various calibrations is compounded by the fact that there exists large, poorly understood systematic discrepancies, in the sense that strong-line abundances generally deviate from T e -based abundances.Furthermore, empirical calibrations and theoretical calibrations disagree with each other.Moustakas et al. (2010) found that empirical calibrations underestimate the "true," T e -based metallicity by ∼0.2-0.3 dex, while the theoretical calibrations yield abundances that are too high by the same amount.Consequently, the absolute uncertainty in the nebular abundance scale is a factor of ∼5 (Kewley & Ellison 2008).Unfortunately, the physical origin of this systematic discrepancy remains unsolved.As such, when comparing trends between multiple calibrations, qualitative trends should be compared, not quantitative differences between the calibrations.
While there are numerous strong-line calibrations in the literature with which to estimate abundances, our narrowband filters limit us to those calibrations that use the metallicity-sensitive R 23 parameter (Pagel et al. 1979): The advantage of R 23 as an oxygen abundance diagnostic is that it is directly proportional to both principal ionization states of oxygen, unlike other diagnostics that have a second-order dependence on the abundance of other elements like nitrogen or sulfur (e.g., Pettini & Pagel 2004;Pilyugin et al. 2010).The disadvantage of R 23 for our purposes is that it must be corrected for stellar absorption and dust attenuation.The most serious complication is that the relation between R 23 and metallicity is famously double-valued.Metal-rich objects lie on the upper R 23 branch, while metal-poor objects lie on the lower branch; the transition between the upper and lower branches happens around a metallicity of 8.4-8.5 dex and is called the turn-around region.
To add to the difficulty of using R 23 as an abundance indicator, there also exists a dependence on an ionization/excitation parameter.The ionization parameter, q, is defined as where S H 0 is the ionizing photon flux per unit area and n is the local number density of hydrogen atoms (Kewley & Dopita 2002).The ionization parameter is usually measured through the ratio of the optical oxygen lines, Some calibrations have attempted to take this into account (e.g.McGaugh 1991; Kobulnicky & Kewley 2004), while others have not (e.g.Zaritsky et al. 1994).The excitation parameter, P , was defined by Pilyugin (2000) in an analogous way to the ionization parameter: Before we estimate oxygen abundances, we investigate the behavior of these strong-line ratios.Figure 7 shows the strong-line ratios R 23 and O 32 as a function of radius for each of the galaxies in our sample.For comparison, the spectroscopic data from CHAOS is plotted as well.The regions in M101 have an R 23 that strongly varies with radius: starting out low for the inner regions, reaching a peak at about 1.4R 25 , and decreasing further out.This likely reflects the non-monotonic relation between R 23 and metallicity, coupled with M101's radial metallicity gradient.This will aid in determining which branch a given H II region is on when we estimate metallicities.Interestingly, the satellites do not have a strong radial gradient in R 23 , which remains roughly constant throughout their populations, equivalent to the regions at intermediate radius of M101.
The ionization state of the H II regions in our sample, as measured by O 32 , shows more scatter than does R 23 .Noticeably, we measure more low ionization regions at larger radii in M101 than CHAOS does.The CHAOS data requires measurable auroral [O III]λ4363, which limits that dataset to only the brightest regions at large radii.Since our detections are determined by measured stronger, collisionally-excited lines, we are able to measure fainter, and thus more H II regions at any radius.Similar statements can be made about the satellites.

Abundance Calibrations & Gradient Shape
As alluded to in Section 4.2, there are poorly understood systematic discrepancies among strong-line calibrations: empirical calibrations generally yield oxygen abundances that are factors of 1.5-5 times lower than abundances derived using theoretical calibrations (Kennicutt et al. 2003;Garnett et al. 2004;Bresolin et al. 2004Bresolin et al. , 2005;;Shi et al. 2006;Nagao et al. 2006;Kewley & Ellison 2008;Moustakas et al. 2010).This being an unsolved issue, we have chosen to estimate our metallicities using three strong-line calibrations: the theoretical calibrations of Kobulnicky & Kewley (2004, hereafter KK04) and McGaugh (1991, hereafter M91); and the empirical calibration of Pilyugin & Thuan (2005, hereafter PT05).This allows us to bracket the range of oxygen abundances one would derive with other existing strong-line calibrations.
All three of these calibrations rely on the R 23 parameter to estimate metallicities and thus have two solutions, one for the upper branch (u) and another for the lower branch ( ).For the KK04 calibration, we have and where x = log(R 23 ).The ionization parameter q in cm s −1 is given by log(q) = 32.81− 1.153y 2 + z(−3.396− 0.025y + 0.1444y 2 ) × [4.603 − 0.3119y − 0.163y 2 + z(−0.48 + 0.0271y + 0.02037y where z = 12+log(O/H) and y = log(O 32 ), where O 32 is given by Equation 5. Note that Equations 7-9 must be solved iteratively for both the ionization parameter and the oxygen abundance; convergence is usually achieved in ∼3 iterations.
For the M91 calibration, we use the functional form as parameterized by Kuzio de Naray et al. ( 2004): and where P is given by Equation 6.Given these calibrations, we now ask what shape the oxygen abundance gradient should take.As mentioned in the introduction, either simple exponential models or broken exponential models are plausible descriptions for metallicity profiles in spiral galaxies.An investigation of Figure 7 shows that the H II regions in M101 have R 23 ratios that vary strongly with radius.The rapidly increasing R 23 values in M101's inner disk are well explained by the fact that inner regions of bright spiral galaxies have high metallicities and negative gradients (e.g.Croxall et al. 2016).However, beyond 1.3R 25 the double-valued nature of the R 23 -O/H relation makes it uncertain whether the data is better explained by a simple exponential model or a broken exponential model.Rather than assigning a branch and deriving metallicities, we use a Bayesian approach to investigate which model best describes the observed R 23 data.In this Bayesian analysis, we first adopt an abundance model: a simple exponential with a gradient and zero point or a broken exponential with a gradient and zero point for the inner region of the galaxy and a break radius followed by a flat outer region.Additionally, both models have a parameter that gives the abundance scatter at a given radius, resulting in a total of three parameters for the simple exponential and four parameters for the broken exponential model.For a given model, we then adopt a calibration that maps oxygen abundance to R 23 , i.e., the calibration of KK04.We adopt uniform priors on the values of each parameter, using the Python package emcee (Foreman-Mackey et al. 2013) to solve for the fit parameters of each model.We report our priors and best fit parameters in Table 3.
Using the best fit parameters from the Bayesian analysis, we generate probability density contours for the R 23 distribution in each model and overlay our observed data for M101 in Figure 8. Comparing the two model contours to the data, we see that both models fit the R 23 data of the inner portion of M101 relatively well, but they differ in the outskirts.In the simple exponential model, the continuing decrease in metallicity results in probability contours which bend downwards in the outskirts, while in the broken model the flat abundances in the outskirts leads to a flattened R 23 profile over that same region.Neither of these trends are precisely described by the observed R 23 datapoints, leaving ambiguity as to the choice of the best model.Indeed, in reality, the proper model may lie in between these two extremes: a weak but still negative abundance gradient in the galaxy outskirts, similar to what has been found for other galaxies (M83, Bresolin et al. 2009;NGC 4625, Goddard et al. 2011;NGC 1058, Bresolin 2019;).
For a more quantitative analysis, we calculated the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) for both models (Akaike 1974;Schwarz 1978; but see Hastie et al. 2016 for an introduction).The AIC penalizes complex models less, while the BIC penalizes complex models more, so an agreement of the two quantities is a good statement of which model truly fits the data best.Both the AIC and BIC indicate that the simple exponential model fits the data slightly better than the broken exponential model.However, even with deep photometry, the lack of points at large radii makes the "goodness of fit" hard to pin down.This motivates a detailed spectroscopic follow up for the outer disk H II regions of M101 to determine their R 23 ratios and oxygen abundances more accurately.In the following, we will calculate abundances for and discuss possible physical origins for both types of gradients in M101.

Oxygen Abundances
Given the calibrations in Section 4.3, we now estimate the metallicities and reasonable uncertainties for each H II region in our sample.As such, we follow the procedure outlined in Moustakas et al. (2010).The goal of this procedure is to estimate the oxygen abundances and uncertainties for regions that lie directly on the R 23 -O/H relation as well as those that lie off of the relation but are statistically consistent with being on the relation given measurement uncertainties.We refer the reader to Moustakas et al. (2010, in particular Fig. 6) for the details, but provide a brief summary here.
For each H II region in our sample, we calculate a metallicity for both the upper and lower branch solutions of a given calibration.We then repeat this for 500 trials in a Monte Carlo algorithm sampling within the given the measurement uncertainties, then no solution exists and these objects are rejected.
Having estimated metallicities and their uncertainties for each region and each calibration, we have to make a decision to which branch of the R 23 -O/H relation each region belongs.A common criterion is to utilize the ratios [N II]/Hα and [N II]/[O II] to break the degeneracy on a region-by-region basis (Contini et al. 2002;Kewley & Ellison 2008), but the [N II] lines are unavailable to us.Instead, we can only make general statements about the populations as a whole.As we have seen in Section 4.3, there are two general shapes that the abundance gradient can take which requires two different assumptions for the regions in the outskirts.

The Gradients of NGC 5477 and NGC 5474
Given the trends evident in Figure 7 for the two satellite galaxies, it is relatively straightforward to assign an R 23 branch to each galaxy's H II regions.Both satellite galaxies' H II regions have R 23 values consistent with being on the upper branch (see Figure 7).We perform a Bayesian linear fit on each radial abundance gradient using the Python package emcee assuming flat priors and a Gaussian likelihood function.The radial abundance gradients are shown in Figure 9 and the fit parameters are reported in Table 4. calibration given the small amount of regions dominated by scatter; the average abundance gradient for NGC 5477 and NGC 5474 are (0.04 ± 0.02) dex kpc −1 and (0.01 ± 0.00) dex kpc −1 , respectively.This is understandable given the small sizes of the satellites, ∼2 kpc for NGC 5477 and ∼6 kpc for NGC 5474, being that both galaxies are likely well-mixed.The relatively flat abundance gradient of NGC 5474 could be surprising given the likely interaction it had with M101 (e.g.Waller et al. 1997;Linden & Mihos 2022).However, given the relatively small size of NGC 5474, gas mixing timescales over sub-kiloparsec scales being on the order of 10-100 Myr (Roy & Kunth 1995), and the estimated age of the interaction being 300 Myr (Mihos et al. 2018;Linden & Mihos 2022), it would be reasonable to assume that NGC 5474 had enough time to smooth out any largescale variations.
Despite the difficulty in investigating dwarf galaxies' oxygen abundances due to their small size and typically low luminosities, a few studies have been able to incorporate dwarfs into their analyses.The CALIFA survey analyzed galaxies down to a luminosity limit of M z < −17 (Sánchez et al. 2014), which includes galaxies like NGC 5474 but excludes NGC 5477.They found that a subsample of galaxies that are undergoing interactions have shallower abundance gradients, regardless of the stage of interaction, and independent of any other quality (morphology, mass, etc.).Those galaxies had an average gradient of ∼ −0.025 ± 0.035 dex R −1 25 (Sánchez et al. 2014), consistent with the gradients we measure for these two satellites.This seems to support the scenario that NGC 5474 has interacted with M101.
However, caution must be taken when overinterpreting these trends.While it does appear that interacting galaxies do tend to have flatter abundance gradients than non-interacting galaxies (e.g.Kewley et al. 2010;Rich et al. 2012;Sánchez et al. 2014), this is not as striking a difference in dwarf galaxies.Numerous studies have found that regardless of the abundance estimator used, there is strong scatter within dwarf galaxies but a lack of any radial abundance gradient within the uncertainties of the abundance estimator (e.g.Pagel et al. 1978;Roy et al. 1996;Hunter & Hoffman 1999;Kniazev et al. 2005;van Zee & Haynes 2006;Lee et al. 2007, but see Pilyugin et al. 2015 for an alternative view).Therefore it is more likely that NGC 5474 and NGC 5477 have abundance gradients consistent with what is expected for dwarf galaxies.

The Gradient(s) of M101
While the satellite galaxies' abundance determinations was relatively straightforward, we now turn to the more complex situation of M101, where we consider both the simple and broken exponential models.Figure 10 shows the abundance gradient for M101.Beyond 1.3R 25 , both upper and lower branch solutions are shown for each individual H II region connected by a dotted line to represent the two general abundance shapes assumed.If we assume that those H II regions are located on the more metal-poor, lower branch, then we recover a simple exponential gradient across all calibrations consistent with spectroscopic studies (see Table 4 for the fit parameters).For comparison, CHAOS reported a gradient of −0.462 ± 0.024 dex R −1 25 (−0.029 ± 0.001 dex kpc −1 ), corrected for our difference in R 25 and distance to M101.Adopting this simple exponential fit would imply that M101's abundance gradient is consistent with a simple application of the insideout galaxy growth scenario.
To recover the broken exponential shape, we assign the outer disk H II regions to the more metal-rich, upper branch.As mentioned before, this decision is mo-

Single Exponential Broken Exponential
Figure 10.The KK04 oxygen abundances of M101 fitted with a simple exponential (solid black line) and a broken exponential line (dashed black line).Regions in the turn-around region (orange circles), upper branch (blue triangles), and lower branch (green inverted triangles) are distinguished.Outside 1.3R25, we display both the upper and lower branch solutions connected by a dotted line for each individual H II region.Characteristic error bars are ±0.02dex for the upper branch, ±0.04 dex for the lower branch, and ±0.2 dex for the turn-around region.As indicated in the plot, a break in the abundance gradient was detected at (15.4 ± 0.5) kpc.
tivated by the radial trend of R 23 for M101 shown in Figure 7 and based on evidence from other galaxies in the literature (see Section 1).To identify the location of the break in the radial abundance gradient of M101, we fit the abundances with two exponential fits joined by a break (a segmented linear regression).The break position was estimated iteratively following the method outlined in Muggeo (2003) and implemented by the Python package piecewise-regression (Pilgrim 2021).Starting with a model of a line with a term that incorporates a change in gradient between two segments of a piecewise function and given an initial guess of the break position, this model is fitted to the data with ordinary linear regression and bootstrap restarting (Wood 2001) is used to find a new break position.The process is iterated until the position of the break converges.In this way an estimate on the uncertainty of the break position is also returned.This process is conceptually similar to that used in Scarano & Lépine (2013).
Figure 10 shows the result of this process for the KK04 calibration and Table 5 records the results of the fits to all of the calibrations.We find a radial break in the KK04 and M91 calibrations at 0.96 ± 0.03 R 25 and 0.89 ± 0.03 R 25 (or 15.4 ± 0.5 kpc and 14.4 ± 0.4 kpc), respectively.We further distinguish between the "inner" and "outer" portions of M101 as defined by those H II regions lying within that calibration's radial break and those lying without.While the fitted gradient in the outskirts of M101 shows a slightly positive slope, given both the large scatter in derived abundances and the large uncertainties in abundances near the R 23 turnaround region, this is likely more consistent with a simple flattening of the gradient.
Figure 11 compares our fitted gradients using the KK04 abundances to other abundance gradients found in the literature.As mentioned in Section 4.2, there is an unexplained deviation between T e -based abundances and strong-line abundances of ∼0.5 dex (e.g.Kewley & Figure 11.A comparison of strong-line (blue lines) and Tebased (orange lines) abundances for M101.Included are our simple and broken exponential fits for the KK04 abundance estimator, as well as the strong-line abundance gradient from Hu et al. (2018) who also used the KK04 abundance estimator.The Te-based abundance gradients are from Croxall et al. (2016), Li et al. (2013), andKennicutt et al. (2003).In all cases from the literature, their gradients have been converted to our choice of R25 and our assumed distance.Ellison 2008;Moustakas et al. 2010).Therefore, the discrepancy between the zero points of the spectroscopic abundances (orange lines) and the strong-line abundances (blue lines) is explained by that difference inherent to abundance estimations.Correcting for our choice of R 25 and our assumed distance, the slope of our simple exponential and the slope of the inner portion of the broken exponential is slightly steeper than that of the spectroscopic studies (Kennicutt et al. 2003;Li et al. 2013;Croxall et al. 2016), but consistent with the other strong-line abundances (Hu et al. 2018).The position of the radial break is also consistent with previous studies (∼18 kpc; Hu et al. 2018).
For the PT05 calibration, we note that this calibration is empirical and can strictly only be applied to regions that are similar to those that were used to make the calibration.For PT05, that restricts H II regions to those with P ≥ 0.4 (Pilyugin & Thuan 2005;Moustakas et al. 2010).Implementing this restriction to our data results in a good fit and a break location that is consistent with the theoretical calibrations of 0.65 ± 0.04 R 25 (10.5 ± 0.7 kpc).The fit with the constraint on the excitation parameter is what is reported in Table 5.We will explore the physical meaning of the break in Section 6.

AZIMUTHAL ASYMMETRIES
As our data completely samples the disk of M101, an important test to undertake is to search for any largescale azimuthal variations.The fundamental processes that govern chemical enrichment in a galaxy operate on different scales.For instance, oxygen is first produced in high-mass stars and then blown to parsec scales through winds or supernovae and then dispersed to kiloparsec scales by various mixing processes (Roy & Kunth 1995).These should be represented in azimuthal deviations from a simple radial gradient since the timescale for galactic differential rotation to homogenize the interstellar medium ( 1 Gyr) is longer than both the oxygen production timescale (< 10 Myr) and mixing timescales (10-100 Myr; Roy & Kunth 1995).
In Figure 12, we plot the spatial distribution of H II regions in M101, NGC 5477, and NGC 5474.Each region is colored by its KK04 metallicity.An inspection of Figure 12 shows no clear morphology to the abundance patterns beyond a radial gradient in all three galaxies.The lack of azimuthal abundance variations in the satellite galaxies is likely evidence of galactic homogenization for the same reasons as for their flat radial abundance gra- dients (Section 4.5).There are minor variations along the spiral arms of M101, possibly similar in nature to the detached outer arm at ∼R 25 (position angles 230°t o 290°) referred to as "arc A" by Li et al. (2013).Relative to the abundance scatter outside of "arc A," the abundance scatter they detected is at the 2σ level.Additionally, there does not appear to be any major difference between H II regions lying in a spiral arm and those lying between arms.As a first test towards a quantitative understanding of azimuthal variations, we investigate the reported asymmetry in abundances of M101 by Kennicutt & Garnett (1996).They suggested that H II regions in the SE have lower oxygen abundances compared to H II regions in the NW as obtained from the use of the R 23 parameter.Later studies reported that there did not appear to be a large-scale azimuthal difference between the SE and NW halves of M101 (Kennicutt et al. 2003;Li et al. 2013).To test this with our expanded sample of H II regions, we divided M101 into two halves with respect to the major axis: a SE part (position angles 37°to 217°, including 270 regions) and a NW part (complementary position angles, with 390 regions).We investigated the radial dependence of R 23 as well as the KK04 abundances and found that for both quantities, the distribution of data points is virtually the same between the two subsamples.
However, given M101's strongly distorted disk, it is possible that there exists azimuthal differences on smaller angular scales or even between different halves of the galaxy.To test this, we measure the abundance gradient in M101's outer disk using azimuthal quadrants advanced incrementally by 5°.The results for each calibration are plotted in Figure 13.We see that regardless of the calibration, the western and southwestern regions of the galaxy have significantly steeper slopes than the rest of the galaxy, a difference of about 0.4 dex.As a further test, we varied the size of the wedge used, ranging from halves to octants, and found the same general trends.If there does exist an azimuthal bifurcation in M101, Figure 13 seems to suggest that the split is almost east-west rather than SE/NW.

DISCUSSION
In the previous sections, we have estimated abundances for ∼700 H II regions throughout the M101 Group.When fitting a profile to the radial abundance gradient of M101, we find ambiguity between two possible models; both a simple exponential and a broken exponential gradient fits the data reasonably well depending on the method used.Additionally, we have also presented evidence of azimuthal abundance variations between the northeastern and southwestern halves of M101.In the following, we will explore these results in the context of M101's asymmetry and interaction history.Galaxy interactions and mergers are fundamental to galaxy formation and evolution and likely lead to strong mixing of abundances in the gas-phase metallicity.As laid out by Toomre & Toomre (1972) and summarized in Toomre (1977), two interacting galaxies will tidally interact during even a minor flyby.Tidal interactions and their associated shocks are thought to trigger gas flows throughout the disk (e.g.Bushouse 1987;Kennicutt et al. 1987;Barnes & Hernquist 1996;Mihos & Hernquist 1996).These gas flows should not only cause bursts of star formation, but it is reasonable to assume that radial mixing of metals in a galaxy should occur as well (e.g.Lacey & Fall 1985;Schönrich & Binney 2009;Bresolin et al. 2009;Rupke et al. 2010;Bird et al. 2012).Therefore, any pre-interaction gas-phase metallicity that exists should effectively be averaged out and imbued with new metals from post-interaction starbursts or accreted gas from the companion galaxy.
An example of where these mixing processes as a result of a flyby occurred is found in the M101 Group, where evidence of a flyby abounds throughout the group.The asymmetric disk of M101 has long been believed to have arisen from an interaction (Beale & Davies 1969;Rownd et al. 1994;Waller et al. 1997).In recent years, low surface brightness (LSB) tidal features have been detected in the outskirts of M101 (Mihos et al. 2013) with colors and stellar populations consistent with a burst of star formation ∼300-400 Myr ago (Mihos et al. 2018).Linden & Mihos (2022, hereafter the Linden model) recently modeled the system with a grazing retrograde encounter between M101 and NGC 5474 with closest approach occurring roughly 200 Myr ago.Their simulation accurately reproduces the tidal morphology of M101 and shows the large radial motion imprinted in the galaxy's disk.Of particular importance, the simulation shows that the young starburst population seen in the NE Plume came from star formation triggered at the contact point when the galaxies collided.
Given the strong observational evidence cited above for an interaction between M101 and NGC 5474 joined with the Linden model, what does this mean for the shape of the radial abundance gradient of M101?Since we have two possible models for the abundance gradient of M101, we will investigate physical mechanisms for producing each gradient.6.2.A Simple Exponential Gradient. . .A simple exponential abundance gradient is predicted by the inside-out growth scenario for disk galaxies (e.g.Scannapieco et al. 2009), and it is likely that M101 followed suit.Since the metallicity of the stars and gas in a galaxy are decoupled, a hint of the original abundance gradient should be present in the old stellar populations.Using HST photometry, Mihos et al. (2018) were able to estimate the metallicities of old RGB stars in the NE Plume and stellar halo at projected distances of 36 kpc and 47 kpc, respectively.Assuming a solar oxygen abundance of 12 + log(O/H) = 8.69 (Asplund et al. 2009), their metallicities correspond to an oxygen abundance of 7.39 and 6.99 for the NE Plume and halo, respectively.If instead of measuring the stellar halo, they measured the extreme outer disk, the slope between these two stellar measurements (−0.036 dex kpc −1 ) is consistent with the slopes of the inner disk found with the strong-line calibrations, hinting that the original abundance gradient of M101 was a simple exponential gradient as expected.
However M101, like most galaxies, has not evolved in isolation.Although the M101 Group has been called one of the poorest groups in the Local Volume (Bremnes et al. 1999), it being dominated by M101 with a scattering of low-mass companions, its galaxies likely have interacted with each other.As described in Section 6.1, M101 and its most massive companion, NGC 5474, likely interacted between 200-400 Myr ago (Mihos et al. 2018;Linden & Mihos 2022).This interaction was strong enough to influence the shape of M101's disk, but was it strong enough to influence its oxygen abundance?
Numerous studies in the literature have investigated the integrated gas-phase metallicity of local galaxies and whether there is a dependence on environment (e.g.Shields et al. 1991;Skillman et al. 1996;Mouhcine et al. 2007;Ellison et al. 2009;Kewley et al. 2010;Peng & Maiolino 2014).While galaxies that reside in denser environments are found to have higher metallicities, this effect is small; the abundance difference between galax- ies in dense and underdense regions is <0.1 dex (Peng & Maiolino 2014;Pilyugin et al. 2017).Relatively little attention has been paid to the abundance distribution within galaxies, i.e., their gradients, and any environmental dependence.Lian et al. (2019) found that there was no change to the shape of the abundance gradient with respect to environment for massive galaxies which confirms earlier work (Kewley et al. 2006;Kacprzak et al. 2015;Pilyugin et al. 2017).It seems that even in a strongly interacting environment like a galaxy group or cluster, the abundance profiles of galaxies do not deviate strongly from the simple exponential model.Perhaps M101 would still retain its original simple exponential abundance gradient despite the interaction with NGC 5474.
However we must express caution at accepting that conclusion for an asymmetric galaxy like M101.The strong asymmetries present in M101's disk will cause its properties to vary as a function of azimuth.For example, Mihos et al. (2013) constructed surface brightness profiles for M101, both azimuthally averaged and in octants.While the azimuthally averaged profile was a Type I profile (Pohlen & Trujillo 2006), i.e., fitted by a simple exponential gradient, the octants revealed the asymmetry.The northeastern portion of the disk was consistent with a Type III (upbending) profile and the southwestern portion of the disk was consistent with a Type II (downbending) profile.
Not only are the two halves with different surface brightness profiles the same halves that have different radial abundance slopes (Figure 13) in M101, but this link between the abundance gradient and surface brightness profile is found in other galaxies as well (e.g.Webster & Smith 1983;Edmunds & Pagel 1984;Vila-Costas & Edmunds 1992;Sánchez et al. 2014;Pilyugin et al. 2014;Bresolin & Kennicutt 2015).Importantly, Type III surface brightness profiles, such as that found in the northeastern portion of M101's disk, are connected with abundance flattening at large radii (Marino et al. 2016).Galaxies with such profiles are thought to have experienced episodes of inside-out growth in their outer disks in the past (Bresolin et al. 2012;Marino et al. 2016) likely caused by tidal interactions with nearby galaxies (Watkins et al. 2019).

. . . Or A Broken Gradient?
Given the asymmetries present in M101's disk, let us now investigate the other abundance profile presented in the current study: the broken exponential gradient.As mentioned before (Section 4.3), abundance profiles with breaks at ∼R 25 beyond which the abundances plateaus have been found in many galaxies (e.g.Bresolin et al. 2009;Goddard et al. 2011;Bresolin et al. 2012;Sánchez et al. 2014;Sánchez-Menguiano et al. 2016).However, the physical nature of this flattening is still an open question since the ability to measure the metal content of the outermost regions of spiral galaxies is relatively recent.Numerous physical processes have been proposed, ranging from radial gas flows (e.g.Lacey & Fall 1985;Goetz & Koeppen 1992;Schönrich & Binney 2009), resonance scattering with spiral density waves (Sellwood & Binney 2002), perturbations by satellite galaxies (Quillen et al. 2009;Bird et al. 2012), or highly efficient star formation in outer disks (Esteban et al. 2013).Likely it will be some combination of these effects that produce the observed flattening at large radii.
In the specific case of M101, we believe that the interaction with NGC 5474 created a burst of star formation in the outer disk of M101 (Mihos et al. 2013(Mihos et al. , 2018)).The burst would have imbued the ISM with oxygen and other metals created in the cores of short-lived massive stars, having the effect of raising the oxygen abundance in M101's disk.However, if the induced gas flows from the interaction were radially restricted, then instead of the entire disk rising in abundance, only the outer disk would increase in abundance, resulting in a flattening of the radial abundance gradient.We propose that this barrier is caused by the corotation radius.As predicted by theory (Mishurov et al. 2009) and seen in simulations (Lépine et al. 2001) and observations (e.g.Marochnik et al. 1972;Mishurov & Zenina 1999;Dias & Lépine 2005;Amôres et al. 2009;Elmegreen et al. 2009;Lépine et al. 2011), corotation has a "pumping out" effect: the corotation barrier produces gas flows in opposite directions on the two sides of corotation, with radial gas flows inward inside corotation and outward on the other side.Scarano & Lépine (2013) compiled results from the literature of the corotation radii and metallicity breaks for 16 galaxies and found a strong correlation.Using the M101 rotation curve from Roberts et al. (1975), Scarano & Lépine (2013) derive a corotation radius for M101 of (15.6 ± 2.2) kpc, adjusted for our assumed M101 distance.This is consistent with the breaks in radial abundance gradients we have found for all three metallicity calibrations we used.Therefore, it is likely that the corotation of M101 is responsible for creating and maintaining the abundance gradient break after the interaction with NGC 5474.
Important to this theory is the longevity of the corotation radius.If the corotation radius was not stationary on long timescales, i.e., if spiral structure was transient (Sellwood & Carlberg 1984;Sellwood 2011), then any break in radial abundance would be smoothed out.Lépine et al. (2011) estimated the minimum lifetime of the Galactic corotation to be about 3 Gyr.Based on the same argument, this can be considered an estimate for the lifetime of corotation for M101 as well (Scarano & Lépine 2013).This timescale could give the inner and outer disks of M101 enough time to chemically evolve differently before and after the interaction with NGC 5474.However, it should be mentioned that the dynamics of the interaction could reasonably change the stationary nature of corotation on short timescales.
Once again, we must express caution at readily accepting this scenario given our data analysis techniques.Pilyugin (2003) noted that abundance breaks found using the R 23 strong-line method may not be real.This is caused by a systematic error depending on the excitation parameter P (Eq. 6) where the method overestimates oxygen abundances in low excitation regions (Kinkel & Rosa 1994;Castellanos et al. 2002;Pilyugin 2003).
To combat this excitation-dependence, we used the so-called P -method presented by Pilyugin (2003) and still found a radial abundance break at 0.65 ± 0.04 R 25 (10.5 ± 0.7 kpc), albeit more interior than the breaks found with the M91 and KK04 calibrations.

CONCLUSIONS
We have conducted a narrowband emission-line survey to measure the oxygen abundance of the entire M101 Group, including M101 and its two major satellites, NGC 5477 and NGC 5474.Using narrowband filters that target Hα, Hβ, [O III], and [O II], we have detected a total of ∼930 H II regions across these three galaxies.Our multi-band detection scheme lets us reject contamination due to foreground and background objects.After correcting for stellar absorption features, [N II] and [S II] emission, and internal extinction, our data shows a good correspondence to spectroscopic flux measurements (Croxall et al. 2016).Additionally, our sample extends down to fainter H II regions in all three galaxies across a wide range of ionization states compared to spectroscopic surveys.
We use our emission line data to derive H II region oxygen abundances for the M101 Group using the R 23 ratio (Pagel et al. 1979).Specifically, we estimated the oxygen abundances with three strong-line calibrations: two theoretical calibrations (McGaugh 1991;Kobulnicky & Kewley 2004) and one empirical calibration (Pilyugin & Thuan 2005).While these methods have an inherent uncertainty due to the double-valued nature of the R 23 -O/H relation at low abundances, we present a variety of models that span the range of abundance pat-terns consistent with the data.Our main conclusions are summarized below.
1. Of 853 H II regions detected in M101, 11 H II regions detected in NGC 5477, and 71 H II regions detected in NGC 5474, we measured oxygen abundances for roughly 75 % of the H II regions in the M101 Group giving us a total of ∼720 H II regions to trace the abundance patterns in the M101 Group.
2. While the two satellite galaxies' oxygen abundances were best fitted with a simple flat gradient, M101's metallicity gradient required more detailed analysis.We compared two models for M101's gradient, a simple exponential model motivated by spectroscopic analyses of M101 and a broken exponential model motivated by spectroscopic analyses of a number of galaxies.Both models provide a plausible description of the data, with some suggestion that an intermediate model of a radial break to a slightly shallower gradient in the outskirts rather than a true flattening might be warranted.
3. Quantitatively, the simple exponential model for M101 has an exponential decline of −0.03 dex kpc −1 , while the satellites were consistent with having shallow or flat gradients within the uncertainties.The broken exponential model for M101 has an exponential decline of −0.04 dex kpc −1 inside a break radius of 14 kpc beyond which the gradient plateaus at 8.5 dex.
4. We also searched for any azimuthal abundance variations across M101's disk.We did not find any strong variations along the spiral arms nor between arm/inter-arm regions.However, we did find that the west and southwest regions of M101 have significantly steeper radial gradients than elsewhere in the galaxy by ∼0.4 dex R −1 25 .
5. We discussed both models for M101's abundance gradient in the context of its interaction history and dynamics.simple exponential gradients are believed to be a feature of inside-out growth and independent of environment.Furthermore, previous spectroscopic analyses of M101 have supported a simple exponential gradient.However, given the strong NE/SW asymmetry of M101's disk, we encourage spectroscopic studies that sample H II regions throughout the outer disk to confirm the azimuthal dependency of the radial gradient.6. Broken exponential gradients are supported by recent studies of other galaxies, both individually and in large surveys (e.g.Bresolin et al. 2009;Sánchez et al. 2014).While numerous physical processes have been proposed to explain the presence of an abundance break followed by flattening, we argue that the possibility of a radial break in M101's disk could be a result of the dynamical corotation barrier (Lépine et al. 2011;Scarano & Lépine 2013) which results in the inner and outer portions of the disk evolving independently from one another.Furthermore, the flattening could be caused by a burst of star formation from M101's interaction with NGC 5474 enriching the outer disk.However, caution must be expressed here as well, since abundances estimated with the R 23 ratio have been known to produce false breaks (Pilyugin 2003).
This work was supported in part by a Towson Memorial Scholarship to R.G. and by grants to J.C.M. from the National Science Foundation (award 1108964) and the Mt.Cuba Astronomical Foundation.A.E.W acknowledges support from the STFC [ST/S00615X/1].This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.R.G. would like to thank S. Linden for sharing their work on the dynamical modeling of the M101-NGC 5474 interaction.We also thank the anonymous referee for a detailed report that helped improve this paper.2013,2018), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), cigale (Noll et al. 2009;Boquien et al. 2019), SciPy (Virtanen et al. 2020), emcee (Foreman-Mackey et al. 2013) , piecewise-regression (Pilgrim 2021)

Figure 1 .
Figure 1.A color-magnitude diagram for sources in the 2MASS All-Sky Catalog of Point Sources (Skrutskie et al. 2006) within a 30 ×30 cutout centered on M101.Sources are colored according to their Hα EW.The box indicated by the solid black lines indicates those sources that are consistent with stars that we subsequently mask in our segmentation routine.

Figure 2 .
Figure 2. The segmentation maps for M101 (left), NGC 5477 (top right), and NGC 5474 (bottom right).Each shaded region is a detected H II region after star removal and deblending as described in the text.The M101 map measures 30 × 30 , the NGC 5477 map measures 5 × 5 , and the NGC 5474 map measures 10 × 10 .In all images, north is up and east is to the left.make the net correction based on a weighted combination of these two ratios, weighting by the off-band flux ratio of each H II region and its surrounding environment in a similar fashion to what was done with the Hβ correction.Finally, we must make a correction for the fact that our Hα on-band filter measures Hα+[N II] emission while the Hα off-band filter measures [S II] emission.The relative [N II]/Hα and [S II]/Hα line ratios vary with metallicity, electron temperature, ionization state, etc., and not necessarily in the same way(Burbidge & Burbidge 1962;Baldwin et al. 1981;Osterbrock 1989).Our previous work(Garner et al. 2021) corrected the Hα emission using a flux ratio of [N II]/Hα = 0.33(Kennicutt 1983;Jansen et al. 2000).However,James et al. (2005) demonstrated that this ratio works best for [N II]selected H II regions, and overestimates the flux ratios of other regions and galaxies as a whole.Instead, we utilize published spectroscopic line ratios for H II regions in M101 to determine the relation between [N II] and [S II].In particular, we use the data of:Kennicutt & Garnett (1996) totaling 41 regions; Li et al. (2013) totaling 28 regions; and the CHemical Abundances of Spirals group (Croxall et al. 2016, hereafter CHAOS) totaling 77 regions.We investigated this data for a simple linear relation between [S II] and [N II], but there was considerable scatter arguing that additional tween [N II] and [S II].In particular, we use the data of: Kennicutt & Garnett (1996) totaling 41 regions; Li et al. (2013) totaling 28 regions; and the CHemical Abundances of Spirals group (Croxall et al. 2016, hereafter CHAOS) totaling 77 regions.We investigated this data for a simple linear relation between [S II] and [N II], but there was considerable scatter arguing that additional parameters may be in play.CHAOS found a strong radial gradient in the N/O ratio for M101 with a flattening at R/R 25 1.35, while there was no radial gradient in the S/O ratio.This suggests that there should be a radial gradient in the [N II]/[S II] ratio, which we show for the combined spectroscopic dataset in Figure 4. Notably, this radial gradient flattens at 1.25R 25 , consistent with the radius noted by CHAOS above, where the ratio assumes a constant value of [N II]/[S II]≈ 10 −0.25 .

Figure 3 .
Figure 3.The narrowband filters used to image the M101 Group overlaid on top of a synthetic H II region spectrum from cigale.Top: the [O II] filters.Middle: the Hβ and [O III] filters.Bottom: the Hα filters.In all cases, the onband filter sits on top of their respective emission lines, while the off-band filter sits redward of the emission line except for the Hβ-off filter which sits blueward.Particular emission and absorption lines are marked by vertical lines.

Figure 4 .
Figure 4.The [N II]/[S II] residuals as a function of radius.Colored points indicate the survey: Croxall et al. (2016, green),Li et al. (2013, orange), andKennicutt & Garnett  (1996, blue).The solid black line indicates the best fit to the data within 1.25R25, with the adjoining solid gray lines showing the 1σ uncertainty.

Figure 5 .
Figure 5.The Hα, Hβ, [O III], and [O II] fluxes plotted against radius from the center of M101.H II regions in M101 are marked by circles and H II regions in NGC 5474 are marked by plus signs; both are colored by local density of points in the plot.H II regions in NGC 5477 are marked by blue stars.The orange lines are the median values for M101 in bins of 0.1R25, the standard deviation on each of the median values is ∼0.6, and the gray dashed lines mark the inner-outer disk of M101.In all plots, the corrections described in Section 3.2 have been applied.
For the [O II], [O III], and Hβ flux comparison, we calculate the line fluxes as outlined in Section 3.2.However, for Hα, where the CHAOS data explicitly reports Hα, [N II], and [S II] fluxes individually, instead of correcting our Hα flux for [N II] and [S II], we compare our uncorrected Hα fluxes to the combined CHAOS fluxes for the lines, i.e., Hα + [N ii] − [S ii].Finally, since the CHAOS fluxes are corrected for internal reddening, we correct our fluxes for reddening following the procedure outlined in Section 3.2.

Figure 6 .
Figure 6.From left to right, top to bottom: the Hα fluxes, Hβ fluxes, [O III]/Hβ ratio, R23 ratio, [O III] fluxes, [O II] fluxes, [O II]/Hβ ratio, and O32 ratio for the 77 H II regions in M101 selected by the CHAOS group.The fluxes are measured in erg s −1 cm −2 .All quantities shown are logarithmic.The points are colored by distance from the center of M101 in R25.The dashed lines are 1 : 1 lines, and the values reported in the legend are the scatter from the 1 : 1 line, σ.The CHAOS spectroscopic flux is measured in 1 slits, while our narrowband fluxes are measured in 3 apertures, leading to the systematic offset seen in the line fluxes.The CHAOS Hα fluxes simulate our Hα flux measurements by adding the [N II] fluxes and subtracting the [S II] fluxes.

Figure 7 .
Figure7.Top: log(R23) as a function of radius using each galaxy's own R25.Fitted lines are given to guide the eye to radial trends.Bottom: log(O32) as a function of radius.In both plots, blue points are regions in M101 in our data, purple circles are regions in M101 from CHAOS, green plus signs are regions in NGC 5474, and orange stars are regions in NGC 5477.Black points with error bars show characteristic errors for both strong-line ratios as a function of radius and ratio.

Figure 8 .
Figure8.Our measured R23 data (black points) for M101 overlaid on the R23 probability density contours for the simple exponential model (left) and broken exponential (right) generated using the best-fit Bayesian parameters described in the text.For both plots, the KK04 calibration was used to generate theoretical R23 values.

Figure 9 .
Figure 9. KK04 metallicities and fitted radial gradients for the two satellites, NGC 5477 (top) and NGC 5474 (bottom).The solid purple and brown lines are the best-fit solutions, and the faded orange and green lines are 100 random walk fits that reflect the posterior distributions.Abundances are calculated via the upper branch of the R23 relation (see text).

Figure 12 .
Figure 12.The spatial distribution of KK04 abundances for M101 (left), NGC 5477 (top right), and NGC 5474 (bottom right).Points are colored by KK04 abundance (color bar to the left), and angular and physical units are given.

6. 1 .Figure 13 .
Figure13.The gradient of the outer regions of M101 in a quadrant rotated by 5°.Only those H II regions that satisfy 0.5 < R/R25 < 1.8 were used in calculating gradients.The calibrations are as shown in the legend.

Figure 14 .
Figure 14.The M101-NGC 5474 interaction from the Linden model.Time progresses from left to right, top to bottom.The orbit of NGC 5474 is traced in black and the stars in both galaxies are colored by local density of points.The red circle is of radius 2 kpc centered on the center of M101.Negative time indicates time before closest approach (∼14 kpc) and the last panel (bottom right) shows the current day arrangement of the galaxies.

Table 1 .
Properties of M101 and Nearby Companions

Table 2 .
Narrowband Imaging Datasets −18 8.0 × 10 −19 (Martins et al. 2005)t al. 2021).While the faintest H II regions in the satellite galaxies can be powered by approximately a single O7V star, the faintest region in M101 can be powered by approximately a single O9V star(Martins et al. 2005).Compare these values to those obtained by CHAOS: 77 H II regions targeted by spectroscopy, with a minimum Hα luminosity of 2.3 × 10 36 erg s −1 , an SFR of 1.2 × 10 −5 M yr −1 , potentially powered by a single O8V star.

Table 3 .
Bayesian Priors and Fits

Table 4 .
Simple exponential Fits

Table 5 .
M101 Broken Exponential Fits The fitted data was limited to only those regions for which this calibration strictly applies, i.e., P ≥ 0.4.Note-The inner/outer distinction refers to the data fitted inside and outside the breakpoint listed in the last column.