Three supernova shells around a young star cluster in M33

Using a specialized technique sensitive to the presence of expanding ionized gas we have detected a set of three concentric expanding shells in an HII region in the nearby spiral galaxy M33. After mapping the kinematics in H{\alpha} with Fabry-Perot spectroscopy we used slit spectra to measure the intensities of the [SII] doublet at {\lambda}{\lambda} 671.9, 673.1 nm and the [NII] doublet at {\lambda}{\lambda} 645.8, 658.3 nm to corroborate the kinematics and apply diagnostic tests using line ratios. These showed that the expanding shells are shock dominated as would be the case if they had originated with supernova explosions. Estimating their kinetic energies we find fairly low values, indicating a fairly advanced stage of evolution. We obtain density, mass and parent star mass estimates, which, along with the kinetic energies, are inconsistent with the simplest models of shock-interstellar medium interaction. We propose that the presence and properties of an inhomogeneous medium offer a scenario which can account for these observations, and discuss the implications. Comparing our results with data from the literature supports the combined presence of an HII region and supernova remnant material at the observed position.


INTRODUCTION
It is recognized that feedback from stellar winds and supernovae has had a decisive effect on galaxy formation and evolution (Dekel & Silk 1986). Galaxy formation models without strong stellar feedback convert their gas into stars too rapidly by two orders of magnitude (e.g. Hopkins, Quataert, & Murray (2011)) compared to observations (Kennicutt 1998), form too many stars by factors up to 10 3 , (e.g. Faucher-Giguère, Kereš, & Ma (2011)) and cannot explain the distribution of heavy elements in the intergalactic medium (e.g. Wiersma et al. (2010)). Approximations to the effects of supernovae on the interstellar medium are now a feature of simulation models (e.g. Marinacci, Pakmor, & Springel (2014))which, as well as making ad hoc approximations to ensure that the feedback is sufficiently strong to reproduce the properties of observed galaxies, have traditionally assumed as simplifying limitations that the surrounding ISM is homogeneous, and supernovae occur singly (Chevalier 1974;Thornton et al. 1998). These limitations have been gradually relaxed, as computing power has increased.
It is known that the real ISM is highly inhomogeneous ( al. 2004) and this has been used in a recent set of models by Martizzi, Faucher-Giguère, & Quataert (2015). The effect of multiple supernovae occurring within a stellar cluster was modelled by Kim, Kim, & Ostriker (2011) and more recently by Hennebelle & Iffrig (2014). Martizzi, Faucher-Giguère, & Quataert (2015) produced models with both of these features. The aim of all the models is to reproduce the balance between star formation enhancement and quenching due to the effects of the supernovae on their surroundings, and thereby to predict the macroscopic properties of galaxies. As is frequent in all issues related to star formation, particularly massive star formation, a key difficulty for modellers is the lack of observational constraints; in this article we report observations which take a step forward in improving them.

OBSERVATIONS
We observed a field on the southern arm of M33 centered on the HII region studied in this article during on 31st October 2014 using the Galaxy Hα Fabry-Perot System (GHαFaS, see Hernandez et al. (2008); Fathi et al. (2008)), and followed this up using the Intermediate dispersion Spectrograph and Imaging System (ISIS) on 25th December 2014. The region in question is listed as BCLMP17-d, the central part of the BCLMP17 complex, in the 1974 Boulesteix HII region catalogue (Boulesteix et al. 1974;Hodge, Skelton, & Ashizawa 2002).
GHαFaS is an integral field spectrometer mounted at the Nasmyth focus of the 4.2m William Herschel Telescope (WHT) at the Observatorio del Roque de los Muchachos (ORM, La Palma), producing a 3.4x3.4 arcmin data cube with seeing-limited angular resolution (spaxel size ∼ 0.2 arcsec) and a spectrum at each point on the object, with 48 channels covering a spectral range close to ∼ 400 km/s, yielding a nominal velocity resolution of ∼8 km/s for Hα. The data reduction procedures for GHαFaS have been well described in the literature. Here we give a very brief summary. GHαFaS does not use an optical derotator so a correction is applied by software, as described in detail in Blasco-Herrera et al. (2010). Velocity calibration, phase correction, sky removal and adaptive binning are performed as described in Daigle et al. (2006) producing the final cube used here.
ISIS is a high-efficiency, double-armed, medium-resolution (8-120Å/mm) spectrograph,at the Cassegrain focus of the WHT, capable of providing up to 4' slit length and 22" slit width. We observed the region with a 10" wide slit using the red arm R1200R grating centred at 6800Å. This allowed us to observe simultaneously the Hα, SII and NII lines with a spectral resolution of ∼12 km/s. We used a slit width of 10" to encompass in a single observation the central part of the region where we detected the bubbles with GHαFaS.
We reduced the spectrum using IRAF. As we use only line flux ratios no flux calibration was necessary.

EXPANSION MAPS
To search for expansive components we analysed the GHαFaS data cube with BUBBLY (Camps-Fariña et al. 2015), a program designed to find and fit multi-component spectral line data cubes. The analysis is automatic, using a mathematical property of Gaussians to detect multiple components, which are then fitted (if they are above a threshold signal-to-noise) and checked to look for pairs of components equidistant in velocity (within uncertainty limits) from the central peak representing the bulk of the HII region. These detections are then mapped in expansion velocity giving "expansion maps", which we use to find the bubbles.
We used the default configuration in the BUBBLY parameter file with a minimum S:N ratio of 5 to detect the expansion components. Because of the relatively low signal-to-noise ratio threshold we added visual inspection of the expansion maps to ensure that detection was spatially coherent. In Fig. 1 we show the line profile of Hα emission which indicates the presence of the three shells, with a central main peak and three pairs of secondary components. This profile comes from only a single resolved element of the velocity map. Figure 2 was produced by mapping all the pixels where a given pair of peaks was detected, which displays the projected extent of the shell; the process was repeated for each of the three shells. We can see that the smallest shell has the most rapid expansion velocity, and the largest shell has the slowest, as expected for shells originating in the centre and doing work on the surrounding interstellar medium during their expansion.

PROPERTIES OF THE BUBBLES
We can measure the radius and expansion velocity directly on the expansion maps and, using the Hα flux of the shells, can calculate their density (Relaño & Beckman 2005): We cannot independently calculate the shell width ∆R, so we took a lower limit arising from the evolution of a strong shock disregarding instabilities (Lozinskaya 1992). A lower limit in shell thickness means a lower limit in the total mass of the shell. The fact that we measure a lower limit in mass will be important later in the discussion.
The flux is measured calculating the fraction of flux the pair of secondary components emits and using a calibrated Hα image of M33 (Hoopes, Walterbos, & Bothun 2001). We chose this method because we are more confident in the calibration of the photometry of this image rather than calibrating GHαFaS.
We can also calculate mass and energy and using these three properties: The age is estimated using the relation of t = 2/7 * (R/v), from the analytic expression for the evolution of a supernova remnant in the snowplough phase (McKee & Ostriker 1977). Using this expression is justified in the following sections, where we also discuss the validity of the values calculated.
Lastly, we calculate the mean pre-shock density for each bubble assuming the shells swept out all of the gas they encountered, dividing the total mass of the shell by the total spherical bubble volume. We disregard the masses of the SN ejecta (a few solar masses), as the shell masses have much larger values.
The ISIS spectrum was used to compute emission line ratios, to probe the origin of the shells. To obtain results for each shell separately we applied a modified version of BUBBLY to each spectrum over the region in the spatial direction for each of these lines:   Notes. See section 4 for details on how the properties were determined. ne is the electron density inside the shell, while n 0 is the average ambient particle density it encountered. The properties marked with * (Age and ambient density) are only valid for a shock evolving in a homogeneous medium. The calculations use a lower limit for shell thickness, this implies that mass and kinetic energy are lower limits while density is an upper limit. Even then, the dependence on shell thickness is √ ∆R, so a factor 10 increase in the thickness would increase mass and energy a factor ∼ 3.
To improve the statistical significance of the ratios we repeated the detection using binning of 2, 3, 5, 6, 15 and 30 elements in the spatial direction, and used the mean value with the standard deviation as the error after eliminating clearly discrepant values.
We could in principle also use the [SII]6716/[SII]6731 ratio as an indicator of density, but unfortunately the values measured are not within the range in which the ratio is sufficiently sensitive to density. Table 1 lists all parameters measured in the bubbles, both from GHαFaS and ISIS data.

Archival properties of the cluster
We have supporting information from the literature obtained by querying the position with VizieR; the presence of a supernova remnant in this HII region detected using radio observations at lower angular resolution is reported by Gordon et al. (Gordon et al. 1999) (Tables 1,2 #90). Also relevant is a survey of HII regions and molecular gas by Miura et al. (2012) which used archival Hubble Legacy Archive images to date the accompanying star cluster by fitting isochrones to a color-magnitude diagram extracted using aperture photometry (see section 3.2.2). This gives an age range of 6-12 Myr for the cluster, which is also reported to have 51 O stars (Table 3 #27). This survey also shows a 2.2 · 10 6 M molecular cloud adjacent to the cluster (Table 2 #24).

Identification as supernova remnants
The identification as supernova remnants was made using forbidden line ratios [SII]/Hα and [NII]/Hα to place our bubbles within diagnostic plots which distinguish between supernova remnants, HII regions and planetary nebulae (Sabbadin, Minello, & Bianchini 1977;Garcia Lario et al. 1991). The classic diagnostic of [SII]/Hα > 0.4 is included in these plots as the region for SNR is delimited by this value on the X axis (see Fig. 3). All three fall well within the defined regions for SNR, though for bubbles 1 and 2 the error reaches the area for photoionized gas. Significantly, the same analysis for the central peak places it, as expected, within the area for HII regions.
Though the ratios indicate otherwise, an alternative origin for the bubbles might be stellar winds, but this is ruled out in practice; as winds are continuously injecting energy since the stars' birth they cannot produce the segregated concentric structure observed, as the individual bubbles would blend and merge into one structure before growing to their current size, which already encompasses most of the cluster. The only possibility for winds to produce bubbles afterwards would be a red supergiant wind accelerated by a subsequent Wolf-Rayet phase (see Freyer, Hensler, & Yorke (2006)), but the stars capable of producing this effect have masses over 25 M and could no longer be present, barring multiple star formation bursts. . Diagnostic diagram used to distinguish between supernova remnants (SNR), HII regions (HII) and planetary nebulae (PNe) based on the relation between the SII to Hα and NII to Hα line ratios. We plot the remnants with their identifier numbers as well as the HII region.

Physical properties
The evolution of supernova remnants is divided into three main phases, an initial free expansion while the mass of the ejecta dominates, an adiabatic expansion once the swept-up material mass is similar to that of the ejecta and a radiative phase, also called snowplough phase, in which the losses by radiation become significant and the gas cools. Given that the initial energy of a supernova is around 10 51 erg our remnants currently have lost some 99% of their energy, which means that they must have been radiative for some time. We can check this using the expressions for the adiabatic phase (Hnatyk, Petruk, & Telezhyns'Kyi 2007): where the left expression is the evolution of the radius and on the right we have the transition time to radiative phase. For our measured ambient density we find a transition time of 54 kyr which already does not fit the data, two of our remnants should not have transitioned yet. Furthermore, even with the faster expansion of the adiabatic phase we find predicted sizes roughly double the measured radii. We will address this apparent inconsistency again shortly.
The mass range for the progenitor stars, calculated from the age of the cluster, is from ∼15 to 19 M , which is not unreasonable if we find three similar mass stars in a cluster of this size. Given the relatively short intervals between the supernova explosions, their progenitors must have had very similar masses.
The swept-up masses of the shells are very interesting. Their values are not consistent with the simplest models but they give us a clue to solve the inconsistency in the physical parameters. In a homogeneous medium the expanding shock should sweep out almost all the surrounding gas into the shell leaving a very rarefied medium, so that a second supernova soon after the first should sweep up hardly any material, lengthening the free-expansion phase. In spite of this, we find that each consecutive remnant has accumulated tens of solar masses of material. This is so even though our value for the thickness of each shell is a lower limit, meaning the masses are also lower limits as they scale with the square root of the thickness. A mechanism is needed to replenish gas inside a remnant for the next shell to encounter sufficient material; the most plausible is mass-loading from cold dense clouds. It has long been suspected that the interior of an HII region is not homogeneous nor fully ionized, and that small molecular clouds can survive for a long time inside it as they slowly evaporate (Bok & Reilly 1947).
Hydrodynamic ablation by shocks or supersonic winds can also disrupt such clouds, so each supernova shell feeds gas to the next one. The effects of the interaction on the shell are also important. Simulations show (Chevalier 1974) that upon encountering a cold, dense cloud, the shell will tend towards radiative evolution, rapidly slowing down and losing energy, greatly advancing its stage of evolution. In the cited paper a simple model with a single density jump at a fixed radius is used, while we expect clusters of cloudlets instead. Either way, the effect of a density jump on the shell is what we have observed in our remnants, so adding the condition of inhomogeneity to our model explains the apparent inconsistencies. It also means that the ages estimated assuming a homogeneous medium are not correct. An inhomogeneous medium slows down the remnants, so our values are upper limits, but they should be within a small factor of the real expansion ages.

DISCUSSION
We have presented results which seem to indicate the presence of three currently expanding supernova remnants within a single young HII region. Our key claim is that the expansion detected must arise from three different objects. Initially, we considered other geometrical configurations, such as typical complex structures in HII regions with several arches and cavities, but we detect very clearly an approaching and receding component for all three along the same lines of sight, so they are either nested or aligned on the line of sight. Given the similar size of the bubbles relative to the cluster, and their concentricity, we believe that nested shells offer the simpler explanation.
We found that the physical properties appeared inconsistent given their masses, kinetic energies and ages. This is solved considering an inhomogeneous medium with embedded dense clumps of molecular gas, which would rob the remnants of energy as well as refill the ambient gas to be swept up by the subsequent supernova explosions. Because of this, our interpretation requires such clumps to survive for a long time inside the HII region.
The survival of clumps inside HII regions has long since been documented, (Bok & Reilly 1947), and they have been found even in extreme examples such as 30 Doradus (Indebetouw et al. 2013), where a complex of molecular clouds is found associated with the brightest star cluster.
The effect of stellar winds and radiation on their parent molecular clouds has been studied in simulations by Rogers & Pittard (2013) and Dale et al. (2014). In the former, three massive ∼30 M stars affect a surrounding 4 pc massive molecular cloud. The main finding is that the winds carve low density channels through which the gas escapes, making the denser regions surprisingly resistant to ablation. Another important result is that dense clumps are largely impervious to supernova shocks. They find that the dense material is either destroyed or pushed outside the original 4 pc radius in just over 6 Myr, which coincides with the lower limit of the age for our HII region. Dale et al. (2014) perform cluster-scale simulations and find results compatible to those of Rogers & Pittard (2013) with a very structured dense medium. They use molecular clouds with a variety of physical parameters finding that the denser, more massive clouds are much harder to disrupt. For lower mass clouds the stars carve irregular cavities on scales of 10 pc. The molecular cloud associated with our cluster is similar to the densest and most massive of Dale et al. (2014)'s. Both studies show that the presence of embedded clumps is expected in clusters born from massive molecular clouds, and that they can easily survive until the stars explode as supernovae. We also need to estimate the rate of cloud evaporation required to refill the gas and reproduce the masses in the observed shell. From the masses and ages of the remnants we require an injection rate of about 25 M per 10 kyr. Dale et al. (2014) find that except at the onset of star formation, winds hardly affect the disruption of the molecular cloud; with photoevaporation dominating this process, so here we need consider only cloud destruction by photoevaporation.
Using the relation from Pittard (2007) for photoevaporation and taking a sample cloud with radius 1 pc at a distance of 10 pc from the centre of the cluster, we estimate the flux of ionizing photons from the cluster using its Hα luminosity and obtain an injection rate of 0.8 M per 10 kyr. We can therefore easily reproduce the required mass injection without an overly large number of clouds, even accounting for the fact that our masses are lower limits.
Feedback from stellar winds and supernovae plays an important role in many processes relevant to galaxy evolution. As well as those mentioned above there is the possible dissipation of nuclear dark matter cusps (Pontzen & Governato 2012), the dissipation of metals both within galaxies (Spitoni et al. 2009) and from galaxies into the intergalactic medium (Heckman, Armus, & Miley 1990), and the enhancement of the infall rate of low metallicity intergalactic gas to the disc, by its interaction with supernova ejectae in the galactic halo (Marasco, Fraternali, & Binney 2012). To incorporate any of these effects into models of the formation and evolution of galaxies, it is clearly important to have not only model inputs to the dynamics of the interaction of the winds and supernovae with the surrounding ISM, but measurements which can help quantify the process. As the stars in massive star clusters form virtually simultaneously, cases of multiple supernova explosions must be common enough to be incorporated into the relevant models.But much more generally to make realistic models mass loading of the expansion via the dissipation of dense cool clumps must be taken into consideration. The triple supernova has allowed us to illustrate its importance, and make quantitative estimates of relevant parameters for much wider application. Sabbadin F., Minello S., Bianchini A., 1977, A&A, 60, 147 Spitoni E., Matteucci F., Recchi S., Cescutti G., Pipino A., 2009, A&A, 504, 87 Thornton K., Gaudlitz M., Janka H.-T., Steinmetz M., 1998 500, 95 Wiersma R. P. C., Schaye J., Dalla Vecchia C., Booth C. M., Theuns T., Aguirre A., 2010, MNRAS, 409, 132