{'fraction_non_alphanumeric': 0.04046362175809043, 'fraction_numerical': 0.02740417718613725, 'mean_word_length': 4.129620908876854, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Feedback from Multiple Supernova Explosions inside a Wind-Blown Bubble

We study the evolution of multiple supernova (SN) explosions inside a pre-exiting cavity blown by winds from massive progenitor stars. Hydrodynamic simulations in one-dimensional spherical geometry, including radiative cooling and thermal conduction, are carried out to follow first the development of the wind-blown bubble during the main sequence and then the evolution of the SN-driven bubble. We find the size and mass of the SN-driven bubble shell depend on the structure of the pre-existing wind bubble as well as the SN explosion energy E_{SN} (= N_{SN} 10^{51} ergs). The hot cavity inside the bubble is 2-3 times bigger in volume and hotter than that of a bubble created by SNe exploded in a uniform interstellar medium (ISM). For an association with 10 massive stars in the average ISM, the SN-driven shell has an outer radius of R_{ss} ~ (85 pc) N_{SN}^{0.1} and a mass of M_{ss} ~ (10^{4.8} Msun) N_{SN}^{0.3}at 10^6 years after the explosion. By that time most of the explosion energy is lost via radiative cooling, while ~10% remains as kinetic energy and ~10% as thermal energy. We also calculate the total integrated spectrum of diffuse radiation emitted by the shock-heated gas of the SN bubble. For the models with 0.1 solar metalicity, the radiative energy loss is smaller and the fraction of non-ionizing photons is larger, compared to those with solar metalicity. We conclude the photoionization/heating by diffuse radiation is the most dominant form of feedback from SN explosions into the surrounding medium.


Introduction
Massive stars inject a significant amount of energy into the surrounding medium in various forms: photoionization/heating by stellar radiation and mechanical energy through stellar winds and supernova (SN) explosions, some of which later transforms into gas thermal energy and radiation. For example, OB stars deposit stellar radiation energy into the interstellar medium (ISM) with the luminosity of 10 4 < ∼ L * /L ⊙ < ∼ 10 6 and the H Lymann-limit photon luminosity (hν > 13.6eV) of 10 47 s −1 < ∼ Q H < ∼ 10 50 s −1 (Smith, 2006). So the radiative feedback of an early type star during the main sequence (τ MS ∼ 10 6 yrs) may amount up to ∼ 10 51−52 ergs, which is comparable to typical SN explosion energy.
These stars also inject mechanical energy into the ISM via stellar winds during the main sequence (MS) and the red supergiant (RSG) or luminous blue variable (LBV) phase, and also Wolf-Rayet stage in the case of very massive stars with M * > 25 M ⊙ (Chu et al., 2004). The MS winds have a typical wind speed of v w ∼ 1000−2500 km s −1 and a mass loss rate ofṀ w ∼ 10 −7 −10 −6 M ⊙ yr −1 , while the RSG winds have much lower speed, v w ∼ 10 − 50 km s −1 , but higher mass loss rate,Ṁ ∼ 10 −5 − 10 −4 M ⊙ yr −1 (Chu et al., 2004;Dwarkadas, 2005;Smith, 2006). Owing to much longer duration the MS wind generates a bigger bubble, while the RSG wind produces a much smaller but denser shell inside the MS bubble created earlier. The mechanical energy luminosity of the MS winds ranges 100 < ∼ L wind /L ⊙ < ∼ 3000 (Smith, 2006).
Massive stars with mass M > ∼ 8M ⊙ end their lives as core-collapse SNe and deposit typically E SN ∼ 10 51 ergs in the form of kinetic energy of the ejecta. Initially the SN ejecta expands freely into the wind-blown cavity and later interacts with the shell produced previously by the winds, generating various shocks, e.g., transmitted and reflected shocks. Most of the explosion energy is thermalized via such shocks. After the shock-heated gas cools adiabatically down to 10 6 K or after the forward ejecta shock catches up with the wind bubble shell, the supernova remnant (SNR) begins to lose a significant fraction of thermal energy via radiative cooling. The final energy budget deposited into the surrounding medium in the forms of thermal, kinetic, and radiative energy depends on the detailed history of radiative cooling and heat transport by thermal conduction, which in turn depends on the size and structure of the preexisting wind bubble as well as on the ambient ISM. For example, if a SN explodes in a dense environment, the explosion energy would be lost mostly by radiative cooling and the impact to the ISM should be minimal. According to Tenorio et al. (1990b), in the case of SN explosions inside wind-blown bubbles in the typical ISM with n H,0 = 1 cm −3 , the remaining kinetic energy fraction is E k ∼ 0.3E SN at the end of SNR evolution, while the remaining thermal energy fraction, E th ∼ 0.2 − 0.4E SN .
Feedback from massive stars should play a significant role in the formation and evolution of galaxies. Especially, feedback from the first-generation SNe affects the formation and evolution of proto-galaxies and the large scale structure in the early Universe (e.g., Kay et al., 2002;Scannapieco et al., 2006, and references therein). Although SN feedback is an important ingredient for galaxy formation, it occurs typically on length and time scales much smaller than grid spacings and time steps of numerical simulations. So such subgrid physics is often treated through phenomenological prescriptions in numerical simulations for galaxy formation and large scale structure formation (Cen & Ostriker, 2006). Either 'thermal' feedback (e.g., Katz, 1992) or 'kinetic' feedback (e.g., Navarro & White, 1993) with a feedback efficiency parameter 0 < ǫ < 1, which is the fraction of SN explosion energy deposited in the form thermal or kinetic energy, were tried with varying degrees of success. For a thermal feedback model, the radiative cooling must be turned off artificially for a period of time in order to prevent the dense gas from re-radiating the deposited energy immediately. Kay et al. (2002) showed that both thermal and kinetic feedback models can alleviate the so-called "catastrophic cooling" problem of galaxy formation. Thus in order to construct a physical subgrid prescription for SN feedback in numerical studies of galaxy formation it is necessary to estimate quantitatively the feedback efficiency parameters, ǫ th , ǫ kin , and ǫ rad for thermal, kinetic and radiation energy, respectively.
Massive stars often form in an association, so multiple explosions of corecollapse SNe can create a large superbubble (Tomisaka & Ikeuchi, 1986;Tenorio et al., 1990a;Smith & Cox, 2001). Multiple SN explosions in the proto-galactic environment, in particular, may have played a significant role in the formation of halo globular clusters. In a self-enrichment model of metal-poor globular clusters, proto-globular cluster clouds are assumed to be enriched by multiple SN explosions from massive first-generation stars contained within the clouds and later second-generation stars form in the cooled shell driven by the SN explosions (e.g., Parmetier et al., 1999, and references therein). However, some have argued that the binding energy of a proto-globular cluster cloud may not be large enough to survive a few SNe, and that the homogeneity of metalicity within a globular cluster indicates that pre-enrichment must have occurred within larger systems (Searle & Zinn, 1978). According to numerical studies of early structure formation, first objects might form at a redshift of ∼ 30 with a typical mass of 10 6 M ⊙ and first-generation stars enrich the protogalactic clouds before they merge into larger clouds of 10 12 M ⊙ via gravitational clustering (Bromm & Larson, 2004). This favors the model that combines the preenrichment by first stars before the formation of protogalactic halos and the self-enrichment by next generation SNe within proto-globular cluster clouds. In this scenario the stability of the cooled shell and further star formation critically depends on how efficiently the SN energy is radiated away during the SNR evolution.
Thus in the present work we attempt to estimate feedback effects of corecollapse SNe, including multiple detonation, through numerical hydrodynamic simulations. We first calculate the creation of a bubble by stellar winds from multiple massive stars (N star = 1 − 50) during the precedent MS stage. The end results of the MS wind stage are adopted as the initial states of the SN explosion simulations. We then make quantitative estimates of the fractions of explosion energy transferred into the ISM in different forms, i.e., thermal, kinetic, and radiation energy. Diffuse radiation emitted by the shock-heated gas, especially, ionizing radiation for hydrogen and helium atoms are calculated.
In §2, the basic features of our numerical simulations, including radiative cooling and thermal conduction, are described. We then present the numerical results of wind bubble models in §3, and multiple supernova models in §4. A brief summary is given in §5.

Hydrodynamic Equations
We solve the hydrodynamic conservation equations of the ideal gas in onedimensional spherical coordinate system, including radiative cooling and thermal conduction : where u is the radial velocity, e = P/[ρ(γ − 1)] + u 2 /2 is the total energy of the gas per unit mass, Λ is the energy loss rate per unit volume, and S is the radial heat flux. A ratio of specific heats is assumed to be γ = 5/3. The gas temperature is determined by T = µm H P/(ρk B ), where k B is Boltzmann's constant, m H is the mass of hydrogen nucleus and µ is the mean molecular weight.
We use a grid-based, Eulerian hydrodynamics code in the spherical symmetric geometry, based on the Total Variation Diminishing (TVD) scheme (Ryu et al., 1993). The TVD scheme solves a hyperbolic system of gasdynamical conservation equations with a second-order accuracy. It has been widely used for a variety of astrophysical problems (e.g., Baek et al., 2003).

Radiative Cooling and Thermal Conduction
We define the cooling rate as where n H is the number density of hydrogen nuclei. For the cooling rate coefficient, L(T ), we adopt the non-equilibrium radiative cooling rate for an optically thin gas that is calculated by following the non-equilibrium collisional ionization of the gas cooling from 10 8.5 K to 10 4 K under the isobaric condition (Sutherland & Dopita, 1993). We consider cases with the solar metalicity ( Z ⊙ ) and with 1/10 of the solar metalicity (0.1 Z ⊙ ).
Although we calculate diffuse radiation emitted by the shock-heated gas, we do not explicitly follow the radiative transfer of diffuse radiation and photoionization of the gas. Instead we assume L(T ) = 0 for T ≤ 10 4 K and set the minimum temperature at T min = 10 4 K. This is in effect equivalent to assuming thermal equilibrium at 10 4 K between radiative cooling and photoionization heating.
After completing the hydrodynamic part, updating hydrodynamical quantities, radiative cooling is applied to the thermal energy per unit volume, E, as a separate step (Baek et al., 2003): where t cool = E/Λ is the cooling time scale and ∆t is the integration time step. This is designed to follow approximately the thermal history of the gas that cools radiatively and approaches the specified minimum temperature in a time scale much shorter than the dynamical time scale (LeVeque, 1997).
In some astrophysical plasmas such as solar corona and stellar winds thermal conduction can be important in heat transport. We include the so-called classical Spitzer conduction for a fully ionized gas (Spitzer, 1962): where ln Λ c = 29.7 + ln n −1/2 e (T /10 6 K) and n e is the electron number density. This classical formula is based on the assumption that the mean free path of the electrons is short compared to temperature scale height. However, when temperature scale height is comparable or shorter than the electron mean free path, thermal conduction reaches the saturation limit and heat flux is no longer equal to k∇ r T . In the saturation limit, heat flux is given approximately by where m e is the electron mass (Cowie & McKee, 1977). So we set S = min[S cl , S sat ] as local heat flux.
With the classical formula the thermal conduction term becomes a diffusion term, − → ∇ r · (k − → ∇ r T ). This type of partial differential equation can be solved implicitly by Cranck-Nicholson Scheme (Press et al., 1992). For the saturated formula, however, this cannot be applied. So we solve the thermal conduction term explicitly and the integration time step is chosen as where c s = (γP/ρ) 1/2 is the sound speed.
In grid-based Eulerian simulations a shock transition is typically resolved over several grid zones with a strong temperature gradient, although the physical shock thickness may be much smaller than the size of one grid zone. So unphysically rapid radiative cooling and heat transport can occur across the numerical shock transition zones. To avoid such numerical artifacts, we identify zones that undergo shocking process and turn off radiative cooling and heat transport by thermal conduction inside those zones.
We neglect magnetic fields and dust in our calculations. If there exist dynamically significant, regular magnetic fields in the ambient medium, the assumed spherical symmetry would brake down. Also coherent tangential magnetic fields would substantially reduce the radial thermal conduction flux. Dust may play some important roles through radiation pressure on dust and dust cooling in the cooled bubble shell in a more realistic situation. Proper treatments of such processes require radiative transfer, chemistry and dust formation/destruction, which are beyond the scope of this study. However, if the bubble shell stays photoionized at 10 4 K by diffuse radiation, as we assume here, then the effects of dust are likely to be minimal. Smith & Cox (2001), for example, included dust cooling in their simulations of multiple SNe and found that the dust cooling effects were only minor.

Wind Models
We first consider the development of a wind bubble to study the pre-supernova circumstellar environment. Massive early-type stars emit a steady stellar wind with constant terminal velocity V w and mass-loss rateṀ w = dM w /dt. We assume that the progenitor star has an initial mass 15M ⊙ and spends ∆t w = 4 × 10 6 yr in the main sequence phase. The MS wind is characterized with a mass-loss rate,Ṁ w = 2.5×10 −7 M ⊙ yr −1 , and wind velocity, V w = 2000 km s −1 (Smith, 2006). Here we consider only the MS wind, since the RSG wind produces a small shell with the radius of a few pc inside the bubble created by the MS wind. The RSG wind shell would be swept up quickly by the SN ejecta during the early free expansion stage, so interactions of the ejecta with this shell can be ignored.
We calculate the evolution of a bubble blown by collective stellar winds from multiple massive stars in an association. The luminosity of the collective wind is the key physical parameter that determines the properties and evolution of the bubble and it is given by We select the number of stars, N star = 1, 5, 10, 30 and 50, because a typical OB association contains 10-100 early-type stars (see Table 1). In a realistic model, one needs to specify the initial mass function of the association to obtain a total wind luminosity. In our simplified model, however, L w , is parameterized with N star . So a larger value of L w can be regarded as a higher wind luminosity either from a more massive star or from multiple stars.
The wind is realized numerically by injecting the gas flow with V w at the inner boundary (r w ) of the spherical simulation grid. The wind density at the inner boundary is given by ρ w = N star · [Ṁ w /(4πV w r 2 w )]. We use N r = 3000 uniform spherical grid zones in 150 pc radius with a spatial grid spacing of ∆r = 0.05 pc. Numerical convergence is often an important issue for hydrodynamic simulations with radiative cooling, because the cooling rate has a strong density dependence (i.e., Λ ∝ ρ 2 ). So in order to test numerical convergence of the simulations WF01 model is calculated with several different grid spacings. Fig.  1 shows the structure of the bubble in the simulations with 3 different grid spacings. The contact discontinuity (CD) between the shocked wind and the shocked ISM is around 23 pc, while the location of the conduction front (CF) is different in the three simulations. The size of the hot cavity (∼ 20pc) is smaller than the so-called Field length λ F = [kT /(n 2 H L)] 1/2 ≈ 56 pc, where T ∼ 10 6 K and n H ∼ 0.01 cm −3 are used (Begelman & McKee, 1990). So heat transport by thermal conduction is dominant over radiative cooling inside the cavity. Heat transport from the hot shocked wind gas to the warm shocked ISM reduces the temperature of the shocked wind gas, triggering radiative cooling and compression. These processes are not followed accurately enough to be converged in the simulations, so the density and temperature distributions are significantly different inside the cavity among these simulations. However, the pressure distribution is very similar in the three calculations and so the dynamical evolution of the wind bubble shell seems to be converged. Radial positions of the CD and outer compression wave seem to be converged within a few percent, while the total energy inside the simulation volume seems to be converged within ∼ 10 %. This fast convergence is achieved, because the density compression factor at the shell is not very large in the simulations (ρ shell /ρ 0 < ∼ 10). We also find that the differences in the density and temperature structures inside the cavity among the simulations with different grid spacings do not affect significantly the evolution of SN explosions that will be described in §4.2.

Evolution of Wind Bubbles
A typical bubble has four distinct regions from the central star outward: 1) freely expanding wind with the density profile, ρ ∝ r −2 , 2) a wind termination shock facing inward and the shocked wind region, 3) a contact discontinuity (CD) between the shocked wind region and the shocked ISM region, 4) a forward shock propagating into the ambient ISM (Weaver et al., 1977). While the termination shock is strong, a compression wave may form, instead of the outer forward shock, depending on the ratio of wind ram pressure to the ambient pressure. Fig. 2 shows the gas density and temperature profiles at t = 10 6 years in the models listed in Table 1. Although the outermost structure is initially developed as a weak shock for all models, it becomes a compression wave with the propagation speed, v < ∼ c s in some models, as the bubble expands out into the ambient medium. For WF01 model, for example, the pressure increases only by a factor of 1.5 across the compression wave. Even for WF50 model the outer shock speeds is only v s ∼ 20 km s −1 . Also the shock transition is much wider than that typically found in adiabatic simulations using the same TVD code. In WF50 model the inner termination shock is at 10 pc, the conduction front at ∼ 25 pc, the CD at 41 pc, and the outer shock wave at R bs ∼ 45 pc. Without thermal conduction, the hot cavity inside the CD would be almost isothermal with T ≈ 10 7.6−7.8 K. But heat transport by thermal conduction enhances radiative cooling near the conduction front, leading to a strong temperature gradient inside the cavity.
Near the CD, the gas density changes by more than two orders of magnitude, so the gas cools quickly down to T min = 10 4 K and forms a shell. The thickness of the shell is rather thick (∆R bs ∼ 0.3R bs ). The continuous injection of wind energy seems to be radiated away efficiently near the conduction front. For a spherical shell with inner radius, R bs,i = αR bs , outer radius, R bs , and uniform density, the shell density can be derived from the condition that the shell mass is equal to the swept-up mass: So the density enhancement factor is small for thick shells found in our simulations. For WF01 model, for example, α ≈ 0.7 and n shell /n H,0 ≈ 1.5.
From Fig. 2 we can see that the size of bubble increases with N star and decreases with n H,0 . Weaver et al. (1977) showed that the bubble would evolve self-similarly and the outer radius of the shell would grow as R bs ∝ L 1/5 w ρ −1/5 0 t 3/5 , if radiative cooling and the ambient pressure were negligible. They also derived a slightly modified time dependence of R bs ∝ t 0.58 by numerical integrations of simplified conservation equations where an approximate cooling rate and an ambient pressure with T ISM =8000K were adopted. In our simulations the ambient pressure is not negligible and the outward moving front is not a shock, so the self-similarity seems to be broken. Fig. 3 shows the evolution of R bs and swept-up shell mass, M bs = (4π/3)R 3 bs ρ 0 , for all the models in Table 1. For our models with L w = N star (3.17×10 35 ergs −1 ) and n H,0 = 1 cm −3 , the radius and mass of wind-blown bubble shell can be approximated by and M bs ≈ (1.6 × 10 3 M ⊙ )N 0.54 star ( t 10 6 yrs ) 1.2 .
These relations are slightly different from what Weaver et al. (1977) derived through simplified numerical integrations described above.
Right panels of Fig. 3 show how R bs and M bs depend on n H,0 or the metalicity Z. These quantities do not scale with the ambient density in a simple way, because the radiative energy loss depends not only on the gas density but also on the detailed history of how heat is transported from the hot cavity to the warm shell via thermal conduction. It seems that for n H,0 < 1 cm −3 the swept-up mass is nearly independent of the ISM density, while R bs ∝ n −1/3 H,0 . For high density medium (n H,0 ≫ 1 cm −3 ), the cooling time scale becomes much shorter than the dynamical time, and so the numerical simulation becomes increasingly difficult. For WH01 model with n H,0 = 3 cm −3 , M bs is about 3 times that of WF01 models. However, this trend cannot be generalized to higher density cases because of nonlinear behavior of heat transport and radiative cooling. For WZ01 with lower metalicity (0.1 Z ⊙ ) the cooling is less efficient and so the results are similar to those of WL01 model with lower ambient density.
Bottom panels of Fig. 3 show the fraction of energy lost due to radiative cooling. The radiation energy is calculated as the energy lost from the simulation volume, E rad (t) = L w · t + E tot,i − E tot (t), where E tot = E k + E th is the total energy inside the simulation volume and E tot,i is its initial value. Here the thermal energy E th excludes the initial thermal energy of the ambient gas. The energy loss due to radiative cooling is significant: about 60 % for WL01 and WLZ01 models and about 60-80 % for WH01 and WF01-WF50 models. For WF models the fraction of E rad is reduced somewhat for higher L w (larger N star ), because higher power wind creates hotter cavities.

Multiple Supernova Models
Next we calculate the evolution of multiple supernovae exploded inside the preexisting wind bubble created by the MS winds. So we take the final structures of the wind bubble simulations at t = 4 × 10 6 yrs, described in the previous section, as the states of the circumstellar medium prior to SN explosions.
Each SN explosion is characterized by the explosion energy, E ej = 10 51 erg, and the ejecta mass, M ej = 10 M ⊙ . Then the ejecta can be approximated by a uniform density core (ρ 1 ) and an outer envelope with a steep power-law density profile (ρ 2 ∝ r −n ), expanding freely with the radial velocity, v = r/t (Chevalier & Liang, 1989). For a multiple explosion where N SN supernovae detonate simultaneously, i.e., total explosion energy E SN = N SN E ej , the density profile of the ejecta is defined by where and v t = 10 3 We adopt n = 10, v t = 3.16×10 3 km s −1 and F = N SN ·(8.14×10 7 g cm −3 s 3 ). We note that overall results of our calculations do not depend sensitively on the specific values of M ej and E ej adopted. As shown below, the evolution of a superbubble is controlled mainly by the total explosion energy E SN rather than E ej . Also a SNR should lose its memory of ejecta mass after the swept-up mass becomes much greater than M ej . Table 2 shows the multiple SN explosion models considered here. The number of SNe is given by N SN 1 and varies from 1 to N star , since N star is the total number of possible SNe inside a given bubble. As in the wind bubble stage, we consider three values of the ISM density: n H,0 = 1 cm −3 as a fiducial value (SNF models), 0.3 cm −3 (SNL01) and 3 cm −3 (SNH01). Since SN explosions can occur sequentially with a time interval, we also consider SNF02a model in which a second SN detonates 10 5 years after the first SN. This will be compared with SNF02b model in which two SNe explode simultaneously. Again SNZ01 and SNZ10 models have a lower metalicity, 0.1 Z ⊙ .
We use the same grid setup as the wind simulations: i.e., N r = 3000 uniform grid zones in 150 pc radius. Again two additional simulations are carried out for SNF01 model with smaller grid spacings (N r = 6000 and 12000) to test numerical convergence. Radial position of the outer shock is converged within 0.5 %, while the total energy is converged within 5 %. The simulations are terminated at t f = 10 6 years. Fig. 4 shows the radial profile of the remnant for the single SN case (SNF01 model). Initially, the SN profile described by Eqs. (14)- (17) is superimposed onto the final output of the corresponding wind bubble model (WF01): the wind termination shock is at 3 pc and the hot cavity extends out to 25 pc (solid lines in Fig. 4). The compression wave has weakened and the shell thickness is large (α ≈ 0.43), so the density of swept up shell, n shell = 1.1 cm −3 , is slightly higher than the ambient density (see Eq.

Evolution of Supernovae Remnants
[11]). The SN blast shock wave passes through the pressure barrier produced by the wind termination shock and travels quickly through the hot cavity. The forward shock first interacts with the bubble shell at ∼ 10 4 years, pushing out the inner shell boundary and generating two secondary waves, i.e., a transmitted shock propagating outward into the bubble shell and a weak reflected compression wave moving inward into the ejecta. A reverse shock forms in the SN ejecta around 5 × 10 3 years and is reflected at the center around 7.5 × 10 3 years. Afterward it becomes a traveling compression wave, because the shock heated gas inside the cavity is hot. When this wave runs into the shell at ∼ 10 5 years, the shell is heated for the second time.
By the termination time (10 6 years) the shell is pushed out by the SNR to R ss =67 pc with inner radius R ss,i =54 pc and the swept up shell gas has T = 10 4 K and n shell ≈ 4 cm −3 , not much higher than the ambient medium.
On the other hand, the cavity is filled with the hot gas of T ≈ 10 6 K and n H ≈ 0.01 cm −3 . The outer shock weakens with expansion and slows down to ∼ c s . Fig. 5 compares the two different models of double explosions, i.e., sequential (SNF02a) and simultaneous (SNF02b) explosion models. In SNF02a model, the second SN explosion is loaded 10 5 years after the first SN explosion. Three snapshots are shown: 1) At 6.0 × 10 4 years the outer shock in SNF02b model has expanded further out than the shock in SNF02a model. 2) At 1.2 × 10 5 years, a little after loading of the second SN in SNF02a model, the second SN shock expands out to 35 pc, still well inside the shell at 42-46 pc generated by the first SN. 3) At 2.0 × 10 5 years the SN-driven shell has expanded slightly further out (∼50 pc) in SNF02b model than in SNF02a model. However, the difference is small enough, so we say the bubble size is similar for the two models.
The left panels of Fig. 6 also show the comparison of the two double explosion models. The outer shock radius and shell mass of the two models are very similar. Here E th and E k are total thermal and kinetic energy inside the simulation volume, respectively, excluding the initial thermal energy of the ambient gas. The time evolution of E th and E k shows that the initial kinetic energy of the SN ejecta is transferred to the thermal energy gradually in the beginning as the blast wave expands, but rather abruptly at ∼ 10 4 years when the blast wave runs into the preexisting wind-blown shell. The thermalization process happens again at later times in SNF02a model, when the second SN detonates at 10 5 years and when its blast wave arrives at the first SN shell at ∼ 10 5.2 years. Interestingly, the latter event occurs at the similar time when the aforementioned compression wave hits the shell in SNF02b model, causing the second period of active radiative cooling.
Most of thermal energy is lost through radiative cooling. Here the radiation energy is calculated by E rad (t) = E tot,i − E tot (t). We note that the initial total energy, E tot,i , includes both the SN energy and the energy of the preexisting wind bubble, so (E tot +E rad )/E SN is greater than one. At the termination time of our simulations, the final energy budget in the two models are similar: 6% in thermal energy, 7% in kinetic energy, and 90% lost through radiation energy. Figs. 5-6 imply that, for given total explosion energy, the different fashion of multiple explosion, i.e., either sequential or simultaneous detonation, does not make significant differences in the size and mass of the final bubble shell.
The right panels of Fig. 6 compare the models with different ambient density and different metalicity. At lower density or with lower metal abundance, the bubble size is bigger and E rad is smaller, compared to higher density models, but the shell mass is similar for all the models considered here.
Simultaneous explosions with 1, 3, 5, and 10 SNe inside a wind bubble generated by 10 massive stars are considered as SNF10a-d models. Initial conditions are taken from WF10 model simulation: the cavity is bounded by the CD at 56 pc and the outer compression wave is at 82 pc. The upper two panels of Fig. 7 show the bubble structure at 10 6 years for SNF10a-d models. More energetic explosion pushes the shell further out and produces a faster shock, leading to a larger and denser shell. The lower panels of Fig. 7 show the profiles of SNRs in a uniform ISM (n H,0 = 1 cm −3 ) for comparison. It shows that SN-driven bubbles are larger and hotter, if SNe explode inside a preexisting bubble. For example, the size of hot bubble for SNF10a-d models is 33-55 % larger than that for SNe inside a uniform ISM, resulting in 2.4 -3.7 times bigger volume of hot gas. Also the cavity temperature is 2-3 times higher in SNF10a-d models. This implies that the volume filling factor and the temperature of the hot bubbles generated by core-collapse SNe would be underestimated by a factor of 2-4, if the preexisting bubble is not properly accounted for.
The left panels of Fig. 8 show the evolution of the bubble properties for SNF10a-d models, while the right panels compare the SNZ10a-d models with lower metalicity. The radial position of the outermost shock increases as R ss ∝ t 0.6−0.8 before the arrival of the blast wave at the shell, and afterward it increases very slowly, R ss ∝ t 0.1−0.15 . The size of SN-driven bubbles depends on both the SN explosion energy and the size of pre-existing wind bubbles. According to the numerical results shown in Fig. 8, R ss and M ss depends rather weakly on E SN , i.e., R ss ≈ (85pc)N 0.1 SN and M ss ≈ (10 4.8 M ⊙ )N 0.3 SN .
Again the evolution of E th , E k , and E rad indicates the complex thermalization history via shocks. The first major episode of radiative cooling occurs when the SN blast wave first hits the wind bubble at t ∼ 10 4.5 years. When the compression wave generated by the reverse shock runs into the shell at t ∼ 10 5.4 years, the shell is heated and compressed, leading to the second episode of rapid cooling. In SNF10d model this behavior is most distinct. For SNF10a-d models, about 80-90 % of E SN is lost via radiative cooling within 10 6 years, while only 10 % remains as thermal energy and < ∼ 10 % as kinetic energy. For SNZ10a-d models with lower metal abundance the radiative cooling is less efficient compared to SNF10a-d models, so about 70 % is lost, 20-40 % remains as thermal energy and 10 % as kinetic energy. Since the shell has a significant kinetic energy at that time, it would continue to expand slowly. Exploration of later evolutionary stage requires more realistic treatments of radiative transfer, photoionization/heating, and non-equilibrium radiative cooling below 10 4 K, which is beyond the scope of this study.
In the self-enrichment model of globular cluster formation, the shell of gas swept up by the blast wave from multiple SNe must cool and be confined gravitationally within the parent proto-globular cloud in order for the enrichment process to be effective (Morgan & Lake, 1989;Parmetier et al., 1999;Parmentier, 2004). Thus the success of this model lies on how efficiently the SN explosion energy is radiated away to form gravitationally bound system. According to Fig. 8 (SNZ10 models), in an average ISM with n H,0 = 1 cm −3 and P ISM = 2.3 × 10 −12 erg cm −3 , massive stars in an association with N star ≈ 50 would generate a bubble shell with R ss ≈ 120 pc and M ss ≈ 10 5.5 M ⊙ . Incidentally this mass is similar to the characteristic mass of globular clusters. If 10 % of explosion energy is left as the kinetic energy of the shell, then E k ≈ 5 × 10 51 ergs. With the gas temperature of T = 10 4 K, the total thermal energy of the photoionized shell is E th ≈ 10 51 ergs. But the gravitational energy, |E G | ∼ GM 2 ss /R ss ≈ 10 50 ergs, is too small to confine the shell gravitationally. Thus the self-enrichment model would not work inside such low density clouds, but may operate in much higher density environment. If we blindly take the simple scaling relations shown in Figs. 6 and 8, i.e., R ss ∝ n −1/3 H,0 , but M ss ≈ constant, then the gravitational energy of the shell would increase with the ambient density as |E G | ∝ n 1/3 H,0 . In a much denser cloud higher ambient pressure would also help the shell to be confined, possibly resulting in a gravitationally bound system.

Diffuse Radiation Emitted by Supernova Driven Bubbles
As shown in the previous section, most of mechanical SN energy is transformed into thermal energy of the bubble shell and then lost via diffuse radiation. So radiative feedback is a significant part of SN feedback into the surrounding ISM. In order to study the photoionization/heating feedback of SN explosions we calculate the diffuse radiation emitted by the bubble.
First, we generate both continuum and line emissivities, f ν = 4πj ν (T )/n 2 H (ergs cm 3 s −1 eV −1 ) and f line = 4πj line (T )/n 2 H (ergs cm 3 s −1 ), for n H = 1 cm −3 as a function of gas temperature for photon energy E ph (eV) = hν = [E min , E max ] = [1.36×10 −4 , 10 8 ], using CLOUDY90 code (Ferland et al. 1998). Collisional ionization equilibrium is assumed and no background radiation field is adopted. The emissivity tables are generated for two values of metalicity, Z ⊙ and 0.1 Z ⊙ . Volume integrated spectrum of the diffuse radiation emitted by the bubble during the integration time t is calculated by where F ν is given in units of ergs eV −1 and line emissivities divided by the width of a given photon energy bin are added to the corresponding continuum emissivity. In practice, only V (T, t) ≡ 4π dt n 2 H (T )r 2 dr is calculated during simulations and then F ν (t) = V (T, t) · [f ν (T ) + f line /∆E ph ] is calculated in a post-processing step.
The upper panel of Fig. 9 shows the evolution of F ν (t) for a single SN model. Most of hard X-ray photons with E ph > 1 keV are emitted early before the blast wave runs into the wind bubble shell (t < 8 × 10 3 years). Afterward, the fluxes of lower energy photons increase with time, as the SN energy is deposited as heat into the bubble shell and radiated away mainly in UV and optical ranges. The lower panel of Fig. 9 compares F ν for SNF01 and SNZ01 models integrated to the final simulation time. In SNZ01 model metal lines are weaker, but H and He recombination continuum radiation is stronger, because the temperature is higher due to lower cooling rate.
We also calculate the integrated energies and photon numbers of diffuse radiation that can ionize hydrogen and helium atoms as follows: where E i = 13.6, 24.6, and 54.4 eV is the ionization energy for H, He, and He + , respectively. These quantities are summarized in Table 3. For SNF01 model, the total number of H ionizing photons is Φ 13.6 ≈ 10 61 and the total energy carried away with those photons is ε 13.6 ≈ 0.55E SN . The total number of H ionizing photons scales roughly with the SN explosion energy as Φ 13.6 ∼ 10 61 N SN .
Diffuse radiation is enhanced overall in the models where a denser and hotter shell is produced. Comparing the results of SNL01, SNF01 and SNH01 models in Table 3, we can see that more diffuse photons are produced at denser environments. A SN with higher explosion energy generates a faster blast wave and produces a denser and hotter shell, emitting harder diffuse radiation.
The results of SNF10a-d models show such trends. On the other hand, the comparison of SNF02a and SNF02b models indicate that the simultaneous double explosion generates more ionizing photons than the sequential double explosion, although the total radiation energy (E rad ) is similar in the two models. SNZ models with lower metalicity produce less ionizing radiation and more non-ionizing radiation, compared to the models with solar metalicity.

Summary
We study multiple supernova (SN) explosions inside a preexisting cavity blown up by stellar winds from massive progenitor stars through numerical hydrodynamic simulations. Calculations are performed with a grid-based Eulerian code in one dimensional spherical symmetry, including radiative cooling and thermal conduction. Heat transport by thermal conduction is important inside the hot cavity, while radiative cooling is dominant in the bubble shell. We consider models with a wide range of physical parameters: the ambient density, 0.3 ≤ n H,0 ≤ 3 cm −3 , the metalicity, Z ⊙ and 0.1 Z ⊙ , the number of massive stars, N star = 1 − 50, and the number of supernovae, N SN = 1 − 10.
In order to study the radiative feedback effects, we also calculate the ionizing diffuse radiation emitted by the SN-driven bubble.
We first calculate the development of a bubble driven by collective winds from N star massive stars with L w = N star · (3.17 × 10 35 ergs −1 ). Characteristic structures of wind bubbles are developed initially: an inward facing termination shock, a conduction front near the contact discontinuity, and a forward expanding shock. The outermost shock is not strong and becomes a compression wave, because the wind ram pressure decreases as the bubble expands. At the end of the main sequence (∼ 4 × 10 6 years), the outer boundary of the bubble shell (i.e., radius of the compression wave) is R bs ≈ (60pc)N 0.18 star , and the inner boundary is R bs,i ≈ (0.7 − 0.9)R bs for n H,0 = 1 cm −3 . So the shell is rather thick and the enhancement factor of the shell density is small, i.e., n shell /n H,0 = 1.5 − 3.7. The mass of the bubble shell is M bs ≈ (2 × 10 4 M ⊙ )N 0.54 star .
The evolution of a supernova remnant inside a wind bubble is much more complex than that in a uniform medium. It depends on the structure and size of the preexisting bubble as well as the total explosion energy, E SN = N SN ·10 51 ergs. The forward shock travels quickly through the cavity and runs into the bubble shell, resulting in a transmitted forward shock. On the other hand, the reverse shock in the ejecta is reflected at the center and becomes a compression wave that travels outward. When this wave runs into the shell, the shell is heated for the second time, followed by the second period of active radiative cooling. The wave bounces off the shell and moves inward with a much reduced amplitude.
At the termination time of the simulations (t f = 10 6 years), the outer shock radius is R ss ≈ 67 pc for a single SN inside a bubble blown by one massive star (SNF01 model), while R ss ≈ 85 pc for a single SN inside a bubble blown by 10 massive stars (SNF10a model). In comparison, the shock radius for a single SN in a uniform ISM of the same density at the same time is about 60 pc. The radii of the hot cavity and the bubble shell can be bigger up to 50 % and the cavity temperature is higher by 2-3 times for SN explosions in a wind bubble considered here, compared to the SN explosions inside a uniform ISM. For multiple SN explosions inside a bubble blown by 10 massive stars, the size and mass of the shell increase rather weakly with the total explosion energy as R ss ≈ (85pc)N 0.1 SN and a mass of M ss ≈ (10 4.8 M ⊙ )N 0.3 SN . For N SN = 50, for example, the SN-driven shell would have typically R ss ∼ 125 pc and M ss ∼ 10 5.3 M ⊙ for the average ISM.
Kinetic energy of the SN is transfered to thermal energy via various shocks. About 20-60 % of the explosion energy is radiated away after the SN shock hits the wind bubble shell, and then 20-40 % is lost after the compression wave crushes into the shell later. Although the final energy budget at the termination time of our simulations depends on both the wind bubble structure and E SN , in most cases with solar metalicity, E rad /E SN ≈ 0.8−0.9, E th /E SN ∼ 0.1 and E k /E SN < ∼ 0.1. Ionizing photons of diffuse radiation account for up to ε 13.6 /E SN ∼ 0.55, while the total number of ionizing photons scales as Φ 13.6 ∼ 10 61 N SN . For the cases with 0.1 Z ⊙ , the radiative loss is reduced a little bit, i.e., E rad /E SN ∼ 0.7, but the fraction of ionizing radiation is reduced a factor of 2-3, compared to the models with solar metalicity. So a larger fraction of energy is lost via non-ionizing radiation at metal-poor environments.
We suggest physically correct prescriptions for SN feedback in numerical simulations of galaxy formation should reflect the natures described in this study.  Fig. 1. Gas density, temperature, and pressure profiles at t = 2.5 × 10 6 years in WF01 model calculated with N r = 3000 (solid lines), 6000 (dotted), and 12000 (dashed) radial grid zones. The density is given in units of ρ o = 2.34 × 10 −24 g cm −3 , and the pressure in units of P o = 9.36 × 10 −8 erg cm −3 .  . Three profiles correspond to the remnants at t = 6 × 10 4 yr, 1.2 × 10 5 yr, and 2 × 10 5 yr, from inside to outside, respectively. The normalization constants are the same as Fig. 3. Fig. 6. Left panels: Radial position of the outer shock, swept-up mass of the SNdriven shell, thermal energy, kinetic energy, and radiation energy (from top to bottom, respectively) are shown for SNF02a (solid lines) and SNF02b (dashed lines) models. Right panels: Same as the left panels except that the results for are shown for SNF01 (solid lines), SNL01 (dotted lines), and SNH01 (dashed), and SNZ01 (long dashed) models.  Radial position of the outer shock wave, swept-up mass of the SN-driven shell, thermal energy, kinetic energy, and the energy lost due to radiative cooling (from top to bottom, respectively) for SNF10a (solid lines), SNF10b (dotted), SNF10c (dashed), and SNF10d (long dashed) are shown. Right panels: Same as the left panels except that Z = 0.1Z ⊙ is assumed. Fig. 9. Upper panel: Integrated spectrum of diffuse radiation emitted by the bubble for SNF01 model at t = 10 4 years (red), 5 × 10 4 years (blue), and 10 6 years (green). Lower panel: Integrated spectrum of diffuse radiation emitted by the bubble in SNF01 model with Z = Z ⊙ (red) and SNZ01 with Z = 0.1Z ⊙ (blue) models at t = 10 6 years.