Bubbles and OB associations

1 Introduction

Stars are born in groups and clusters with a mass function d N / d M ∼ M − α , α = 1.25 – 2.25 in mass limits 300 and 3 × 1 0 7 M ⊙ (e.g. Krumholz et al. 2019). Clusters with moderate and large mass M ≳ 1 0 3 M ⊙ are more or less easily identified, whereas the lower mass limit of M min ∼ 300 M ⊙ is poorly-defined, especially in the inner part of the Galactic disc at the galactocentric distances r ≲ 3 – 5 kpc with a rather crowd stellar environment. The problem exacerbates when opaqueness due to dust extinction is accounted for – it prevents observations of OB-associations towards the Galactic centre beyond 3 kpc (Bronfman et al. 2000). However, the lower limit for stellar clusters is of great importance for understanding star formation process. In particular, it determines the smallest spatial scale at which an interrelation between the star formation rate and the gas mass does not operate (Kruijssen and Longmore 2014).

Wood and Churchwell (1989a,b) first suggested to associate bright far infrared point-like sources with ultracompact HII regions tracing sites of OB star formation. Later on, Wouterloot and Brand (1996, and references therein) and Bronfman et al. (2000) have used molecular (CO and CS) lines to measure the kinematic distances to respective star formation regimes. In both cases, these observations are biased towards dense ( ∼ 1 0 3 – 1 0 5 cm − 3 ) and relatively young stellar clusters, and thus more diffuse ones are missed. From this point of view, complementary approaches to identify sites of star formation would be desirable. In this article, we present arguments in favour of SNe-driven superbubbles and supershells around compact stellar clusters with age t OB ≥ 10  Myr as signatures of an OB association.

A small fraction of stars in clusters are massive and end their life as supernovae (SNe) explosions. As a result, outflows of kpc-sizes are formed around massive stellar clusters (see, e.g. Bagetakos et al. 2011), while such explosions in less populated stellar clusters may form chimneys and worms (Heiles 1984, Heiles et al. 1996, Norman and Ikeuchi 1989, Normandeau et al. 1996). The Galactic interstellar disc can be broken out by even only a few SNe explosions if they are members of stellar clusters located higher than a half of the disc scale height (Bronfman et al. 2000, Zari et al. 2018, Damiani et al. 2019). The signatures of such breakouts can be found in H α surveys (see, e.g. Haffner et al. 2003). Such out-of-plane gas structures can serve to reveal underlying stellar clusters. It is reasonable to argue that their morphology and emission properties are sensitive to characteristics of ambient gas, and as such might be used for rough estimates of the clusters’ mass.

2 Model

We carry out 3D hydrodynamic simulations (Cartesian geometry) of SNe explosions inside a cluster located in the galactic disc. We set up a gaseous disc to be initially in the hydrostatic equilibrium in a gravitational potential (as in many previous papers, see, e.g. de Avillez 2000, Hill et al. 2012, Walch et al. 2015, Li et al. 2017), which consists of two components: a dark matter (DM) halo and a baryonic disc. In cylindrical coordinates with z -axis perpendicular to the disc, the z -component of the gravitational acceleration due to the DM halo is calculated from a Navarro–Frenk–White profile: g DM ( z ) = G M DM ( r ) z / r 3 , where M DM is the DM mass enclosed within radius r , and we adopt virial radius of the halo equal to 200 kpc and concentration parameter c = 12 as for the Milky Way case (Navarro et al. 1997).

The baryonic disc is assumed to be self-gravitating with an isothermal velocity dispersion. The acceleration perpendicular to the disc is g ∗ ( z ) = 2 π G Σ ∗ tanh ( z / z ∗ ) , where Σ ∗ and z ∗ are the stellar surface density and scale height of the stellar disc, respectively. Contribution from a gaseous disc is included following Li et al. (2017), by dividing g ∗ ( z ) by f ∗ = Σ ∗ / ( Σ ∗ + Σ gas ) , which implicitly suggests the gaseous disc to be non-self-gravitating. We assume the cylindrical radial distance R = 3  kpc, for which the value of z ∗ is assumed to be equal to 0.3 kpc as it has been estimated in the solar neighbourhood (Gilmore and Reid 1983), the stellar surface density equals to Σ ∗ = 180 M ⊙ pc − 2 , and gas surface density Σ g = 11.5 M ⊙ pc − 2 . The latter is very close to the observed density in the Milky Way disc at the Solar circle (Kalberla 2003, Kalberla and Kerp 2009), and the one considered in Section 4.3 of Li et al. (2017); the mid-plane density is n 0 = 3 cm − 3 . Our model is consistent with that of gravitational potential built by Kalberla (2003) at the cylindrical radial distance 3 kpc from the Milky Way centre (Figure 7 in Kalberla 2003). Our choice is motivated by the following two reasons: (a) a sufficient number of unresolved stellar clusters are expected at this galactocentric radius (see, e.g. Bronfman et al. 2000), and (b) the relative scale heights of the gaseous and stellar discs are most favourable for SNe explosions in young stellar clusters to produce vertical gaseous structures extending above the disc.

The gas density profile is found from the hydrostatic equilibrium. At large heights, the number density is kept uniform at 1 0 − 2 cm − 3 . Initially, the gas temperature is 9 × 1 0 3 K ; the gas metallicity is kept constant equal to the solar value within the entire computational domain. Standard simulations are performed with a physical cell size of 4 pc in the computational domain with 128 × 128 × 192 cells, which corresponds to 512 × 512 × 768 pc 3 . In several runs, the domain has been expanded to 768 × 768 × 1024 pc 3 and 102 4 3 pc 3 . This resolution is close to the spatial limit sufficient for the Sedov–Taylor phase to capture consistently (Kim and Ostriker 2015): the ratio of cooling radius to injection is equal to or greater than a factor of 3. In our simulations, the cooling radius for solar-metallicity gas behind shock waves with T ∼ 1 0 7 K is l c ∼ 12 / n 0 pc for cooling rate adopted following Vasiliev (2013). We have conducted simulations for one of the models presented in this article with twice better resolution, and we did not find any significant changes in physical quantities of interest.

SNe are distributed randomly in a cluster of radius 10 pc. For the spatial resolution adopted here, a cluster volume consists of around 100 cells. Within this volume, 10–20 SNe explode in ∼ 20 Myr with a rate of ∼ 0.5 – 1 SN per 1 Myr. A typical example of a sequence of SN positions is shown in Figure 1. Even though the position of subsequent SNe explosions is random the overall dynamics of a growing superbubble is stable against a particular realization. This circumstance can be readily understood because the radius of a single remnant from an SN can grow in 1–2 Myr up to R ∼ 30 – 50 pc , i.e. much larger than the cluster radius. From this point of view, the randomization is unnecessary. However, it provides a practical advantage for numerical procedure. Random distribution allows us to easily avoid the numerical strong-shock instability, also referred to as the odd-even instability described by Quirk (1994). We consider three models: in the first, the centre of the cluster is located at the midplane ( z = 0 ), and in the others, it is shifted to some height z : 20 and 60 pc. Both values of the shifted clusters are lower than 1/2 of the disc scale height for the galactocentric distance equal to 3 kpc. The mass and energy of each SN are injected in one cell for a standard spatial resolution of 4 pc or a region of radius 4 pc in case of a better resolution; SN energy is assumed to be 1 0 51 erg. Masses of massive stars – SNe progenitors – in a cluster are distributed randomly within 8–40 M ⊙ according to the Salpeter initial mass function. We begin our run when the most massive SN explodes. The intervals between following SNe correspond on average to the lifetime of massive stars, which relates stellar mass to log ( t l , year ) = 10.04 – 3.8054 log M + 1.0646 log 2 M (Iben 2012). Because we consider a small number of massive members in a cluster (especially for a cluster of 10 SNe), the shape of a cumulative bubble slightly depends on the initial spatial SNe distribution; however, the cluster radius is quite small and the differences in SNe spatial positions within the cluster are not crucial for our purposes so that the global dynamics of the bubble is determined by the disc properties and the number of SNe only.

Figure 1

Projected positions of SNe explosions into ( x , y ) , ( y , z ) , and ( x , z ) planes (left to right), dashed circles show the radius of a cluster, colour bar encodes the time when a next SN in the sequence explodes, and squares show projections of the numerical cells. The exact positions are generated randomly and do not coincide with the cell centre. However, the energy of an exploded SN is assumed to be injected into the entire cubic cell.

The code is based on the unsplit total variation diminishing approach that provides high-resolution capturing of shocks and prevents unphysical oscillations. We have implemented the Monotonic Upstream-Centred Scheme for Conservation Laws (MUSCL-Hancock scheme) and the Harten–Lax–van Leer contact (HLLC) method (see, e.g. Toro 1999) as an approximate Riemann solver. This code has successfully passed the whole set of tests proposed in Klingenberg et al. (2007). For instance, the Sedov blast test for a single SN remnant is reproduced well (see Figure 1 in Vasiliev et al. 2017). Growth of an SN remnant radius at ejecta-dominated phase and its transition to adiabatic expansion afterwards (see Appendix A in Vasiliev et al. 2017) is in a good agreement with the well-known analytic and numerical solutions (e.g. Truelove and McKee 1999). The hydrodynamic driver of our code is identical to that described in Khoperskov et al. (2014), where the results of the other several tests can be found. The boundary conditions are periodic for the x - and y -directions, and outflowing for z . The time step is chosen as a minimum value from the Courant–Friedrichs–Lewy condition and cooling time over the whole grid.

Simulations are run with radiative cooling processes. A tabulated non-equilibrium cooling function is taken from Vasiliev (2011, 2013). The fitted function is obtained for gas cooling isochorically from 1 0 8 K down to 10 K. The non-equilibrium calculation (Vasiliev 2011, 2013) includes kinetics of all ionization states of H, He, C, N, O, Ne, Mg, Si, Fe, as well as kinetics of molecular hydrogen at T < 1 0 4 K . We apply a diffuse heating term representing the photoelectric heating of dust grains (Bakes and Tielens 1994), which is thought to be the dominant heating mechanism in the interstellar medium, namely in the warm neutral gas considered here (Wolfire et al. 1995). In our simulations, the heating rate is assumed to be time-independent and exponentially decreasing in the vertical direction with the scale height of the ISM disc. Such an assumption is sufficient to stabilize radiative cooling of ambient gas at T = 9 × 1 0 3 K . Any deviation of the heating rate in the unperturbed gas violates the balance between cooling and heating and stimulates thermal instability, resulting in redistribution of gas mass in the interstellar disc (see, e.g. in de Avillez 2000, Hill et al. 2012). In order to avoid contamination of such effects, we follow Li et al. (2017) to assume the heating rate exponentially decreasing upwards across the entire computational domain. Thus, in the absence of strong perturbations from SNe explosions the initial density and temperature distributions remain in the hydrostatic equilibrium within more than 30 Myr with a high accuracy – deviations do not grow above a few percents. Thus, the initial hydrostatic equilibrium of the interstellar gas layer under the external gravitational potential from stellar disc and DM is highly stable, and possible small deviations do not contaminate overall dynamical scenario.

3 Results

3.1 Dynamics of bubbles

Vertical slices in the ( x , z )-plane in Figures 2 and 3 present the density and temperature evolution for two clusters consisting of 20 and 10 SNe, respectively, with the centres located in the midplane ( z = 0 ).^[1] Figure 4 presents an example of a 3D appearance of density and temperature structure in the bubble corresponding to the model in Figure 2. At the age of ∼ 10 Myr , the bubble around the cluster of 20 SNe expands to the radius comparable to the disc scale height. Afterwards, it continues expanding predominantly in the vertical direction. At the time t = 12 Myr , the bubble already shows evidence of forming an outflow, which occurs at t ≃ 14 Myr . Furthermore, the major part of the energy injected by next SNe is transferred directly into the halo (see plates for t = 16 and 18 Myr), because of a very low gas density in the hole created at previous stages. One expects that after t ∼ 24 Myr , when all SNe in the cluster are exhausted, the hole is filled by gas from the disc within the next ∼ 1 Myr , while hot gas above the disc keeps its high temperature during the next several Myr. The fragments from the destroyed shell will fall on to the disc within a much longer period (Fraternali 2017). Note that the gas unperturbed by SNe explosions keeps the initial hydrostatic distribution even at late times t = 18 Myr , as seen, e.g. in the distributions at all heights close to the vertical borders of the box in Figures 2–6.

Figure 2

2D slices showing gas density (log( n , cm − 3 ), upper panels) and temperature (log( T , K ), lower panels) for a cluster of 20 SNe located at the midplane ( z = 0 ) at 12, 14, 16, and 18 Myr (from left to right) since the first SN has exploded. The cluster centre is located at ( x , z ) = ( 0 , 0 ) , the coordinates are in kpc.

Figure 3

The same as in Figure 2, but for a cluster consisting of 10 SNe. The cluster centre is located at ( x , z ) = ( 0 , 0 ) , the coordinates are in kpc.

Figure 4

3D image of the growing bubble shown in Figure 2 at the age 18 Myr – left and right panels depict the density and the temperature structure, correspondingly; for the sake of clarity, the density is limited in a range n = 1 0 − 4 − 1 cm − 3 , the temperature is within T = 1 0 3 – 1 0 5 K . In order to eliminate rather strong contaminations from the gaseous disc its contribution has been subtracted in the left panel; the thin elliptically shaped artefact remained from the subtraction procedure. We use the YT package (Turk et al. 2011) to produce this figure.

For a less populated cluster, with 10 SNe, the time delay between subsequent explosions is close to merging between SN remnant and the ISM.^[2] After each merging phase, a weak perturbation of density/temperature appears. For instance, such a relic from previous explosions is found in the left panel of Figure 3, the largest one located in the disc has a radius around 100 pc. The gas density inside this relic is close to the unperturbed value in the disc so that a subsequent SN forms a shell in a gas perturbed by previous SNe – the inner circle with an enhanced density shell. The bubble consists of several hot regions at t ≲ 12 Myr . Afterwards, these regions merge at t ∼ 14 Myr and form a narrow hot channel perpendicular to the disc. Subsequent SNe transfer their energy through this channel to higher heights. Although the hot gas reaches distances more than the scale height, the delay between SNe is longer than the time needed for filling this channel by gas from the disc ( t ∼ 16 Myr ), and next SN explodes in gas with a density close to typical in the galactic midplane ( t ∼ 18 Myr ). Two features distinguish this case from the one with higher SN rate: (i) even though the gas evacuated by the initial episode of SNe explosions falls back and fills the cavity the disc height does not relax to the initial height, but becomes thinner around the position of the SNe cluster, (ii) small hot regions located at ∼ 3 scale heights above and below the cluster in the disc. Both structures are short-lived and practically disappear in the next several Myr. Thus, low-populated clusters at radial distances of ∼ 3  kpc from the Galactic centre hardly leave observable relicts/remnants in the gas distribution.

As mentioned above, stellar clusters are spread in heights above/below the midplane. SNe in such clusters may explode at heights that are considerably greater than the cluster radius. Young ( ≤ 5 Myr) stellar clusters have been found up to ∼ 100 pc in the solar neighbourhood (Zari et al. 2018). Such height is roughly half of the gaseous scale height at the solar galactocentric distance. As the stellar disc becomes thicker closer to the Galactic centre, one can expect that young stellar clusters may be located at similar or even higher positions. Thus, it might be interesting to consider SNe explosions in stellar clusters located close to the midplane and slightly larger than a half of the scale height.

Figures 5 and 6 show the density evolution for two clusters consisting of 10 SNe located at 20 and 60 pc above the midplane, respectively. These are around 1/5 and 3/5 of the scale height at a radial distance of ∼ 3 kpc from the Galactic centre. One can see that bubbles formed by these SNe clusters are merged to the ISM at around t ∼ 16 – 18 Myr , similar to those in the midplane (Figure 3). However, for clusters located above/below the midplane new remarkable structures form above/below the disc at ∼ 1 – 3 scale heights: they contain warm diffuse gas with n ∼ 0.1 cm − 3 and exist for a longer characteristic time ( t ∼ 18 Myr ) than those formed at similar heights for midplane clusters (see the right panel in Figure 3). Moreover, the asymmetry in the positioning of these structures – either above or below the midplane may provide additional information about the height at which the driving energy source – an SNe cluster – is located. One can note that for the cluster located at 60 pc the gas perturbed by SNe explosions extends close to the computational domain boundaries at t ∼ 18 Myr (see right panel in Figure 6). However, it is not important for further evolution as soon as the gas is subsonic, and its mass flux is negligible.

Figure 5

The same as in Figure 2, but for a cluster consisting of 10 SNe located at height z = 20 pc above the midplane, i.e. the cluster centre is located at ( x , z ) = ( 0 , 0.02 ) , the coordinates are in kpc.

Figure 6

The same as in Figure 2, but for a cluster consisting of 10 SNe located at height z = 60 pc above the midplane, i.e. the cluster centre is at ( x , z ) = ( 0 , 0.06 ) , the coordinates are in kpc.

3.2 Emission from gas

The structures described above may be detected in various lines and wavelength ranges. To calculate emission features the hydrogen radiative coefficients and Gaunt factors are taken from Draine (2011), we use the ionization gas composition corresponding to the cooling rate for a given temperature (Vasiliev 2013). With these assumptions we simulate emission and velocity distributions along each line-of-sight and obtain emission and velocity maps (see also Vasiliev and Nath 2015, Vasiliev and Moiseev 2015, Vasiliev et al. 2017). Figures 7–10 present the H α intensity, ionized gas velocity dispersion, X-ray intensity in the 0.7–1.2 keV band, and HI column density for the models shown in Figures 2–6 at time t = 16 Myr (right middle panel). Note that Figures 7 and 8 show the maps for the whole disc, whereas Figures 9 and 10 present the emission features only above the midplane. That is because the gas below the midplane is only slightly perturbed by SNe clusters located above the midplane. In each model, one can note the region close to the midplane with slightly enhanced emissivity. These features are produced by the density perturbations originated from the merging phases of early SN explosions (see the density slices in Figures 2–6).

Figure 7

The H α intensity lg( I ( H α ) , erg cm − 2 s − 1 arcsec − 2 ), ionized gas velocity dispersion ( σ , km/s), X-ray intensity in the 0.7–1.2 keV band lg(I(0.7–1.2 keV), keV cm − 2 s − 1 str − 1 ), and HI column density (lg(N(HI), cm − 2 ), panels from left to right) for a cluster of 20 SNe located at the midplane ( z = 0 ) at 16 Myr since the first SN has exploded.

Figure 8

The same as in Figure 7, but for a cluster of 10 SNe located at the midplane ( z = 0 ) at 16 Myr since the first SN has exploded.

Figure 9

The H α intensity, ionized gas velocity dispersion, X-ray intensity in the 0.7–1.2 keV band and HI column density (from left to right) for a cluster of 10 SNe located at height z = 20 pc above the midplane at 16 Myr since the first SN has exploded. Only the maps above the midplane are shown.

Figure 10

The same as in Figure 9, but for a cluster of 10 SNe located at height z = 60 pc above the midplane at 16 Myr since the first SN has exploded.

For the cluster consisting of 20 SNe low-intensity H α structures can be found up to ∼ 3 scales of height, or ∣ z ∣ ∼ 0.3 kpc (left panel in Figure 7), and the velocity dispersion of these structures is higher than ∼ 50 km s − 1 (left middle panel), then, such gas has a specific locus on the “H α intensity – velocity dispersion” diagram – an area that corresponds to both low H α intensities and high values of velocity dispersion (Martínez-Delgado et al. 2007, Moiseev and Lozinskaya 2012, Egorov et al. 2014, Vasiliev and Moiseev 2015, Egorov et al. 2017). Hot gas outflowing from the disc has sufficiently high temperature to emit in X-rays. The HI column density in the disc is significant for attenuating X-ray intensity below 0.3 keV, so that here we consider more energetic range, 0.7–1.2 keV (right middle panel). One can find X-ray spots located symmetrically above and below the midplane. Note that the X-ray spot in the midplane disappears within less than 1 Myr after exhaustion of SNe explosions, because of efficient gas merging in the midplane, whereas the X-ray spots above and below the disc remain for next several Myr. Their shapes are very similar to X-ray structures observed in the Galactic centre (see, e.g. Figures 6 and 1 in Extended data by Ponti et al. 2015, Ponti et al. 2019, respectively). The X-ray structures formed by less SNe number (Figures 8–10) look like narrow chimneys and similar to small filaments seen in Figure 6 by Ponti et al. (2015). The gas in these structures is highly ionized, such that only a few lumps of neutral gas can be found in HI maps (right panel).

There are no distinguishable emission structures after SNe explosions in a 10 SNe cluster located at the disc midplane (Figure 8). During the evolution, one may only find bubbles formed by individual SN explosions, because the time delay between subsequent explosions is longer than the duration of the merging phase. Figure 8 presents the distributions for time moment after 4 Myr when the last SN has exploded (for the evolution see Figure 3), so that these emission distributions present the latest period of the merging phase.

For a cluster displaced out of the midplane several remarkable gaseous structures are appeared (Figures 5 and 6), and they might be found in H α and HI emission (see left and right panels in Figures 9 and 10). They are found at ∼ 3 scale heights ( ∼ 0.3 kpc ), though even more remarkable for clusters at moderate heights ∼ 1 / 5 − 1 / 2 of scale height, or ∼ 20 – 60 pc . Clusters at smaller ( ≲ 20 pc) or larger ( ≳ 60 pc) heights do not give any structures, because at smaller heights SNe cannot form collective outflow (see the description above and Figure 3), whereas at larger ones the swept-up mass is too small. Note that for clusters at larger heights, the energy injected by SNe is easily channelled into the vertical direction, and the volume of high-velocity gas is greater than for clusters at lower heights. As a result, a high velocity dispersion region appears at larger heights (see left middle panel in Figure 10).

3.3 Dust emission in hot gas

Grains in hot gas are heated both by collisions behind shocks and by ultraviolet (UV) starlight produced by stellar population in the disc. To calculate the total infrared (IR) emission of dust grains we assume a steady-state heating by UV photons and collisions. We model the dust population consisting of graphites and silicates of five sizes: 30, 100, 300, 1,000, and 3,000 Å distributed according to Mathis et al. (1977). The dust-to-gas ratio is taken equal to 1/120. The local UV energy density within one scale height of the disc is set the value u ∗ = 1.5 × 1 0 − 13 erg cm − 3 , it is three times higher than the field near the solar neighbourhood (Habing 1968), because the stellar surface density in our simulations is taken about three times higher than that near the Sun as well. The dust absorption and emission coefficients are taken from Draine (2011) and Laor and Draine (1993). The heating rate due to collisions with electrons is adopted from Dwek (1987). Small dust grains can be destructed in hot plasma, we take it into account in the manner described by Drozdov and Vasiliev (2021).

Figure 11 presents the total IR ( λ ∼ 3 –3,000 μ m ) surface brightness (left) and average temperature (right panel) of dust grains in hot gas of the bubble formed by a cluster of 20 SNe located at the midplane at 16 Myr since the first SN has exploded. Note that here we do not depict the dust emission of the disc.

Figure 11

Total IR surface brightness (log( L d , erg s − 1 cm − 3 ), left panel) and temperature ( T d K , right panel) of dust grains in hot gas of a cluster of 20 SNe located at the midplane ( z = 0 ) at 16 Myr since the first SN has exploded.

Large ( a ≳ 300 Å) grains in the hot gas at low heights are predominantly heated by the UV starlight, whereas smaller ones ( a ≲ 100  Å) are kept hotter by collisions in high-temperature interior of the bubble. Because dust emissivity is mainly determined by smaller grains, the IR brightness and dust temperature in hot plasma of the bubble are significantly contributed by collisions. One can note that bright “caps” reach heights of ∣ z ∣ ∼ 0.4 kpc . They are more pronounced in comparison with the total H α emission (see left panel in Figure 7). The dust temperature is distributed similar to the X-ray emission in the range of 0.7–1.2 keV (see third panel in Figure 7). One can see that a significant volume filling factor of hot gas in bubbles can be found in the bubble formed by a cluster of 20 SNe. Therefore, dust emission can be used in identifying clusters consisting of more than 20 SNe or more massive than ∼ 3 × 1 0 3 M ⊙ (assuming one SN per 150 M ⊙ ).

One can note that the dust temperature at low and middle galactic latitudes is usually below 20 K (see middle panel of Figure 9 in Planck Collaboration et al. 2013). These temperature values are thought to be determined by a local UV background (e.g. Section 9 in Draine 2003). At the same time in the all-sky dust temperature map (see middle panel of Figure 9 in Planck Collaboration et al. 2013), there are many small patches with higher dust temperature values: 23–27 K. In Figure 11 one can see that the dust temperature in the “caps” reaches ∼ 25 K and even higher, that is close to the values in such observed patches with typical angular size about several degrees or physical size of 0.2–0.4 kpc at distance of several kpc. A detailed comparison of the simulated “caps” and observed patches will be done elsewhere. However, one can expect that the emission of dusty “caps” can be used in searching SNe explosions in low-mass OB clusters.

4 Discussion

Throughout the article we implicitly assumed that the radial density gradient has negligible effects on the dynamics of growing bubbles and their emission characteristics. Such an assumption is motivated by the fact that in all cases we considered the characteristic scales of SN remnants to merge and form a local vertical outflow are much shorter (from ∼ 1 pc within a cluster to ∼ 100 pc within the ISM disc) as compared to spatial scales in the radial direction (a few kpc for spiral galaxies). On the contrary, the effects of radial density gradient in the disc might have effects in dwarf galaxies where the ISM discs are amorphous and fluffy with the vertical to the radial scale ratios ( ≳ 0.1 ) being larger than in disc galaxies ( ≲ 0.03 ) (see, e.g. Puche et al. 1992). We also assumed the ISM disc implicitly to be in the hydrostatic equilibrium and neglected local turbulent motions. The presence of turbulence leads to a faster washing out of the density perturbations formed due to the merging phase (see density slices in Figures 2–6); however, it is unlikely to change the global dynamics of SN bubble (see Fielding et al. 2018, Vasiliev et al. 2019 for comparison).

5 Conclusion

We have studied the evolution of collective remnants formed by several (up to 20) SNe in the galactic disc. Such a number of SNe can be attributed to small OB clusters. The morphology of a remnant depends on both disc and cluster properties, e.g. disc thickness, midplane density, and height at which the cluster centre is located above the disc midplane, cluster radius and mass function of massive stars, and so on. We have found the following

1. the energy injected by a cluster of ∼ 20 SNe at the midplane is sufficient to break the disc out at radial distance ∼ 3 kpc in the Galaxy; the outflow can reach 2–3 scale heights and emit in H α and X-ray (0.7–1.2 keV) during several Myr;

2. a cluster of ∼ 10 SNe at the midplane is not able to drive gas outflow from the disc;

3. ∼ 10 SNe clusters located at ∼ 0.2 – 0.5 scale heights above the midplane produce H α structures extending ∼ 0.5 – 2 scale heights, these structures can be identified by velocity dispersion of ionized gas;

4. A cluster of ≳ 20 SNe at the midplane can produce the IR-bright dusty “caps” at heights ∼ 0.2 – 0.4 kpc , the dust temperature in the “caps” is ≳ 20 K , which is mainly controlled by small grains heated presumably by collisions in hot plasma of SNe bubble.