© ESO 2020

1 Introduction

The mechanisms leading to the formation and growth of SMBHs in the centers of their host galaxies is an open unsolved problem in astrophysics. The mass of the SMBH is correlated with the mass of the bulge (Kormendy & Ho 2013) and some extraordinary starburst activities of active galactic nuclei (AGN) may be connected to accreting SMBHs (Yesuf & Ho 2020). However, this starburst activity usually occurs in regions enshrouded by dust, hiding completely their emission in optical wave bands. The energy is radiated away usually at infrared wavelengths and galaxies with abundant accretion onto a central SMBH appear as ultra luminous infrared galaxies (ULIRG; Sanders & Mirabel 1996). Interstellar matter can be funneled into central regions of galaxies by mergers and interactions of gas-rich galaxies that trigger the starburst activity. The feedback driven by young and massive stars leads to outflows (Davies et al. 2014), thereby reducing the interstellar matter (ISM) supply to the central SMBH, limited also by the conservation of angular momentum.

Nuclear star clusters (NSCs) are found nestled in the cores of galaxies with SMBHs at the bottom of the galaxy’s potential wells. Their typical half-mass radii are (2–5) pc, with the full extent going out to (5–50) pc and estimated masses are in the range (10⁵–10⁹) M_⊙ (Spengler et al. 2017). NSC masses appear to be related to masses of their host galaxies sharing a relationship similar to that between SMBH masses and galactic bulge masses (Kormendy & Ho 2013). The formation of a NSC follows two possible scenarios. The first assumes that globular clusters formed in the galaxy outside of the galactic center and then spiralled down into the galaxy’s core because of dynamical friction before merging there to form the NSC (Tremaine et al. 1975). However, it appears that there are not enough RR Lyrae stars in the NSC of our Galaxy to be consistent with this first scenario: if the NSC were formed from globular clusters, the expected number of RR Lyrae stars would be almost an order of magnitude higher compared to numbers identified by the Hubble space telescope observations of the center of the Milky Way (Dong et al. 2017). The alternative scenario proposes in situ formation of the NSC developing from gas migrating into the galactic center due to internal processes in the galactic disc or as a result of interactions with other galaxies (Schödel et al. 2008).

We investigate here the evolution of shells expanding into the interstellar medium surrounding the SMBH in the centers of galaxies inside of NSCs. The expanding shells are the result of energy and momentum inserted by young and massive stars of the NSC in the form of winds, radiation, and supernova explosions. The evolution of shells expanding into the interstellar medium resulting from either supernovae or stellar winds has been described in many papers and books, for example, by Ostriker & McKee (1988) and Bisnovatyi-Kogan & Silich (1995). Here we summarize very briefly how the blast wave resulting from a supernova explosion evolves. The first stage is the phase of free expansion followed by the so-called Sedov-Taylor phase characterized by constant thermal energy inside the blast wave driving the expansion. When the radiative losses in the shell of the swept-up gas become important and the heat energy created by the compression of interstellar gas is radiated away, the dense wall of the shell behind the shock shrinks to a thin layer and the structure enters the pressure-driven, thin-shell phase. Later, when the pressure inside the shell drops, the wave is no longer driven by the interior pressure and it only keeps its momentum. This is the so-called snowplow phase. Even later, when the expansion velocity decreases below the local sound speed, the mass accumulation ceases and the thin shell dissolves.

The main intention of this paper is to investigate whether supernovae (SNe) occurring in NSCs can deliver mass into the central region inside of the circumnuclear disk surrounding the SMBHs in the centers of galaxies. We explore the importance of the SN position relative to the SMBH and we analyze how the inserted energy, mass, and the density of the interstellar medium influence the mass delivery to the SMBH. Here we assume that the gasinside the expanding SN remnant is non-relativistic (gas velocity is below a few thousand km s⁻¹), and that shells are able to survive encounters with random density fluctuations and stay coherent when they are deformed by the tidal forces of the SMBH and of the NSC. To find suitable energies and positions of SNe for mass delivery near the central SMBH, in this paper we perform simulations with the simplified hydrodynamical code RING.

In the 3D code RING we assume that the wall of the shell is infinitesimally thin, meaning that its thickness is much smaller than its diameter. In this and future papers, simulations with the fast code RING will map the expanding shell properties for different SNe energies, positions relative to the SMBH and NSC, and densities of the ISM. This approach was developed by Kompaneets (1960) and Bisnovatyj-Kogan & Blinnikov (1982), and has been used by Tenorio-Tagle & Palous (1987), Ehlerová & Palous (1996), Silich et al. (1996), Ehlerová et al. (1997), Efremov et al. (1999), Elmegreen et al. (2002), and others.

The paper is organized as follows: in Sects. 2 and 3 we describe the gravitational potential of the SMBH and of the NSC, the density and motion of the ISM, and the position of the supernova relative to the SMBH. The formation of the thin shell, as simulated in 1D with the hydrodynamical code FLASH, is presented in Sect. 4. In Sect. 5 we introduce the 3D code RING, which is used in Sect. 6, where we show the time evolution of the thin shell. Results are given in Sect. 7, and discussed in Sect. 8. Conclusions are shown in Sect. 9.

2 Gravitational potential of the SMBH and NSC

The gravitational potential of the SMBH, Ψ_SMBH, is given as $${\Psi }_{\text{SMBH}}=-\frac{G{M}_{\text{SMBH}}}{(r_{\text{GC}}+\varepsilon )},$$(1)

where G is the gravitational constant, M_SMBH is the SMBH mass, and r_GC is the distance from the galactic center. Here we assume that the SMBH is located exactly at the Galactic center and that the NSC is similarly centered there. A small constant ɛ has a value of 10⁻⁶ pc preventing infinite values of the potential and of its derivatives. The SMBH mass, M_SMBH, is a free parameter with a value in the interval (10⁶ –10⁹) M_⊙. In this paper, we keep M_SMBH constant, since we assume the mass delivered to it during the interaction is much less than the mass M_SMBH itself. In the case of the SMBH in the center of the Milky Way, we take M_SMBH = 4 × 10⁶ M_⊙.

The stellar density distribution in the NSC, ρ_NSC, is given as a combination of two γ-models proposed byDehnen (1993) in the form suggested by Chatzopoulos et al. (2015): $${\rho }_{\text{NSC}}={\rho }_{{\text{NSC}}_{\text{1}}}+{\rho }_{{\text{NSC}}_{\text{2}}},$$(2) $${\rho }_{{\text{NSC}}_{\text{i}}}(r_{\text{GC}})=\frac{3-{\gamma }_i}{4\pi }\frac{M_{{\text{NSC}}_{\text{i}}}{a}_i}{r_{\text{GC}}^{{\gamma }_i}{(r_{\text{GC}}+a_i)}^{4-{\gamma }_i}},$$(3)

where i = 1, 2 and M_NSCi is the mass of the two individual nuclear stellar subclusters.

For simulations of the shell expanding in the NSC of the Milky Way, following Chatzopoulos et al. (2015) we adopt $$\begin{array}{l}{M}_{{\text{NSC}}_{\text{1}}}=2.7\times {10}^7{M}_{\odot },{\gamma }_1=0.51,a_1=3.9\text{pc};\\ M_{{\text{NSC}}_{\text{2}}}=2.8\times {10}^9{M}_{\odot },{\gamma }_2=0.07,a_2=94.4\text{pc}.\end{array}$$

The total NSC mass is M_NSC = M_NSC1 + M_NSC2. The mass NSC₂ appears to coincide with the one described by Serabyn & Morris (1996). In future simulations, the division between the two stellar subclusters will be varied together with the total value of M_NSC: (10⁵ < M_NSC < 10⁹) M_⊙.

The corresponding gravitational potential of the NSC is $${\Psi }_{\text{NSC}}(r_{\text{GC}})=-{\sum }_{i=1}^2\frac{G{M}_{{\text{NSC}}_{\text{i}}}}{a_i}\frac{1}{(2-{\gamma }_i)}\left(1-{\left(\frac{r_{\text{GC}}}{r_{\text{GC}}+a_i}\right)}^{2-{\gamma }_i}\right).$$(6)

In this paper we neglect the contribution of the ISM to the overall gravitational field, and the total gravitational potential, Ψ(r_GC), is given as the combination of the SMBH and NSC only: $$\Psi (r_{\text{GC}})={\Psi }_{\text{SMBH}}+{\Psi }_{\text{NSC}}.$$(7)

3 Density and motion of the ISM and positions of the supernova

In this initial exploration of how a supernova will affect gas dynamics near a supermassive black hole, we simplify the physical conditions by assuming that the ambient ISM is distributed in a homogeneous way with a constant density, n_out, which is a free parameter. We explore three values: n_out = 10³, 10⁴, and 10⁵ cm⁻³ (Antonini et al. 2015; Neumayer et al. 2020). In subsequent communications we shall also analyze shells expanding in lower or higher ISM densities.

The ambient ISM rotates in the gravitational field of the SMBH and NSC with the plane of rotation parallel to the symmetry plane of the hosting galaxy. The circular rotation velocity, V_rot (R, z), is given as $$V_{\text{rot}}(R,z)={\left(\frac{\text{d}\Psi (r_{\text{GC}})}{\text{d}{r}_{\text{GC}}}{r}_{\text{GC}}\right)}^{1/2}\frac{R}{r_{\text{GC}}},$$(8)

where $r_{\text{GC}}=\sqrt{R^2+z^2}$ is the galactocentric distance and (R, z) are the galactocentric cylindrical coordinates with the rotational axis perpendicular to the z = 0 plane. We assume that in the z = 0 plane the rotation of the unperturbed gas completely balances the gravitational attraction. At the distance z from the plane, rotation slows down since it only balances the component of gravitational force parallel to the plane ofrotation; however, the ISM density is kept homogeneous and independent of the distance z. The rotational velocity as a function of cylindrical coordinates (R, z) is given in Fig. 1.

We analyze the expansion of SNe shells from different initial positions within the NSC close to the SMBH. The initial values (R₀, z₀) are varied inside the volume $$1<r_{\text{GC}}<150\text{pc}.$$(9)

In subsequent papers we shall also consider the SNe motion relative to the ISM at the instant of the SNe explosion.

Fig. 1 Interstellar matter (ISM) rotation velocity, Vrot, is shown as afunction of R for different distances z from the symmetry plane.

4 Formation of the thin shell

Over a few seconds, a supernova injects energy E_tot into the surrounding ISM creating an expanding shock wave. Here we adopt the canonical value E_tot = 10⁵¹ erg; other values of E_tot will be explored in the future. To set the initial conditions for the three-dimensional code RING we first performed one-dimensional simulations with the grid-based hydrodynamical code FLASH v4.3 (Fryxell et al. 2000) using the Piecewise Parabolic Method (Colella & Woodward 1984) with the time-step controlled by the Courant-Friedrichs-Lewy criterion and with equilibrium radiative cooling using the prescription of Schure et al. (2009).

The simulation starts from a small spherical volume into which we inserted SNe ejecta of mass M_ej with total energy E_tot. The initial conditions of the free-expansion phase are given by Tang & Chevalier (2017) with the density ρ_ej and velocity v_ej profiles specified as $${\rho }_{\text{ej}}=\frac{1}{4\pi }\frac{M_{\text{ej}}}{r_{\text{ej}}^3}{\left(\frac{r_{\text{ej}}}{r}\right)}^2$$(10) $$v_{\text{ej}}={\left(6\frac{E_{\text{SN}}}{M_{\text{ej}}}\right)}^{1/2}\left(\frac{r}{r_{\text{ej}}}\right),$$(11)

where r is the distance from the SN explosion, and r_ej is the radius of the spherical volume where the ejecta are distributed. Here we take M_ej = 10 M_⊙ and r_ej = 0.24 pc. The shock wave expands into a warm (T = 10⁴ K) ambient medium of density ρ_out = μ n_out with μ = 0.609 m_H, where m_H is the atomic mass unit, and n_out = 10³, 10⁴, or 10⁵ cm⁻³.

The time evolution of the position of the forward shock, r_sh, and of its expansion velocity, v_sh, and the time evolution of the thermal energy inside of the shock, E_th, and of the kinetic energy of the shock, E_kin, are shown in Fig. 2. After some time a thin shell forms behind the leading shock. At this thin shell formation time, t_sh, the thin shell has radius r_sh and mass m_sh. The thermal energy inside the cavity of theshell, E_th, and the kinetic energy of the shell, E_kin, at time t_sh are given in Table 1. The rest of the initial injected energy, E_tot, was radiated away. At this stage of SN expansion, when the thin shell is formed, we switch to the code RING using the infinitesimally thin shell approximation. In Fig. 3 we show profiles of density, temperature, velocity, speed of sound, and pressure for n_out = 10⁴ cm³ at the time t_sh = 502 yr.

The 3D simulation with the code RING starts from the radius r_sh at the time t_sh given in Table 1 when, according to the 1D FLASH simulation, the thin shell begins to form. Formation of the thin shell is connected to the rapid loss of thermal energy due to the increased density that leads to enhanced cooling. The time of the shell formation was set on the basis of FLASH simulations and the choice was tested by a comparison of FLASH and RING results. Into this initial volume we inserted at the time t_sh the initial thermal energy E_th. The shell of total mass m_sh (see Table 1) gets the initial expansion velocity corresponding to kinetic energy E_kin defined by FLASH simulations. The initial velocities of individual elements have a constant value in all directions from the SN position, meaning that in all simulations of expanding shells we disregard the rotation of the ISM in setting the initial velocities of the shell elements. This is a valid assumption because the shell size at t_sh is still small compared to the scale over which the local velocity field changes.

5 Expansion of the shell: code RING

The expansion of the infinitesimally thin shell in 3D is described with the code RING. The supersonically expanding shell accumulatesthe ambient medium and slows down. The code RING was originally written by Jan Palouš (2D version: Tenorio-Tagle & Palous 1987; 3D version: Palouš 1990) and later extended by Soňa Ehlerová, most recently used in Ehlerová & Palouš (2018). The code RING is currently not publicly available, and a more detailed description will be given on the web page (http:// galaxy.asu.cas.cz/pages/ring). There have been other codes that use the Kompaneets approximation (Silich et al. 1996; Mashchenko et al. 1999); the main difference between these and RING lies in the treatment of conditions inside the bubble.

The infinitesimally thin shell is divided into N_elem elements. They are originally distributed equally across the initial sphere; we use the HEALPIX (Górski et al. 2005) distribution with 192 elements. The subsequent evolution may distort the shell so that the distribution of elements does not remain homogeneous and their sizes (= areas) vary.

For each of the elements we solve the momentum and mass conservation equations. The momentum equation is $$\frac{\text{d}}{\text{d}t}(m_{\text{i}}{v}_{\text{i}})+{\rho }_{\text{out}}[S_{\text{i}}\cdot (v_{\text{i}}-v_{\text{out}})]v_{\text{out}}=(P_{\text{int}}-P_{\text{out}})S_{\text{i}}+m_{\text{i}}g,$$(12)

where m_i is the mass of the ith element, S_i is its surface area (i.e., the surface area multiplied by the normal vector), v_i is its velocity of expansion, v_out is the velocity of the medium outside of it, P_out and ρ_out are the pressure and density of the medium outside of the shell, g is the gravitational acceleration of the NSC and SMBH, and t is the time.

The shell sweeps up the ambient medium as it expands: for each surface element, when the velocity component perpendicular to the element, $v_{\perp }={(v_{\text{i}}-v_{\text{out}})}_{\perp }$, exceeds the local speed of sound, the mass accumulation is given as $$\frac{\text{d}{m}_{\text{i}}}{\text{d}t}={\rho }_{\text{out}}[(v_{\text{i}}-v_{\text{out}})\cdot S_{\text{i}}].$$(13)

When v_⊥ drops below the sound speed in the ISM, the accumulation of mass is stopped and the element continues its expansion in the gravitational potential in a ballistic way without accumulating mass further.

The infinitesimally thin shell approximation carries with it the assumption that the shell has such a high density that it is able to radiate away all the thermal energy produced by compression of the ambient medium by the shock. The pressure inside the cavity of the shell P_int follows the equation $$P_{\text{int}}=\frac{2}{3}\frac{E_{\text{th}}}{V_{\text{int}}},$$(14)

where V_int is the total volume inside the shell.

We disregard the radiative cooling of the hot and low density medium inside the cavity of the shell. Its thermal energy, E_th, follows the adiabatic energy balance equation $$\frac{\text{d}{E}_{\text{th}}}{\text{d}t}=-\frac{\text{d}{V}_{\text{int}}}{\text{d}t}{P}_{\text{int}}.$$(15)

Integrations are done using the Runge-Kutta method (a fourth order, RK4) with an adaptive step size. The surface areas of expanding elements are computed at every time step by interpolation between the positions of neighbors.

To check the importance of radiative cooling from the interior cavity of the shell, we follow the shell expansion with the 1D code FLASH, and we record the time t_cool when the medium inside the shell, originally heated by the reverse shock to high temperatures, cools down to 10⁴ K (see Table 1). For t > t_cool the internal pressure P_int drops sufficiently to be ignored, and the shell enters the momentum-driven snowplow phase. We examine the results of simulations with the code RING when the internal pressure is switched off at t_cool. In this case the RING simulations give the total mass collected in the shell as being only ~ 1% lower compared to the cases when the internal pressure acts all the time. This demonstrates that the importance of cooling from the interior of the shell is small: the shell expansion is mainly governed by momentum with little contribution from the internal pressure.

As an additional test of the validity of the transition from FLASH to RING, in the static medium without the external gravity of the NSC and SMBH, we computed the shell expansion during the first 3000 yr. Figure 4 compares the radius of the shock as a function of time from FLASH and RING. As is visible in Fig. 4, there is a close correspondence in shell radii: the difference at 3000 yr is still below 2% demonstrating a very good agreement achieved between these two very different simulations.

Fig. 2 One-dimensional FLASH simulations. Left panel: time evolution of the position of the forward shock, rsh (r3, r4, r5), and of its expansion velocity, vsh (v3, v4, v5). Right panel: time evolution of the thermal energy inside the shock, Eth (Et3, Et4, Et5), and of the kinetic energy of the shock, Ekin (Ek3, Ek4, Ek5). Individual quantities are labeled with the value of the power index of the ISM density nout = 103, 104, and 105 cm−3.

Fig. 3 One-dimensional FLASH simulations: profiles of the density n, the temperature T, the velocity v, the sound speed cs and the pressure p∕kB for r < 1 pc at time = 502 yr in a medium of density nout = 104 cm−3.

Fig. 4 Shell radius, rsh, as a function of time during the first 3000 yr of expansion. The 1D FLASH is joined by the 3D RING at tsh : absolute values (left panel), relative values of the difference (rsh,F − rsh,R)∕rsh,F (right panel), where rsh,F is the radius of the shell according to FLASH and rsh,R is the radiusof the shell according to RING simulations. The number in the line label shows the power index of the value nout = 103, 104, and 105 cm−3.

6 Later evolution: simulations with the code RING

The evolution of the expanding thin shell after its formation at t_sh is followed with the code RING. To explore expansions from different positions relative to the SMBH, we changed the initial galactocentric distance, r_GC, and the elevation angle, θ_GC, which is the altitude of the SN explosion above the galactic plane as seen from the galactic center. We computed the total mass collected in the shell, m_{sh, tot}, at the end of mass accumulation, when all its elements slow down and expand with subsonic velocity.

As an example, in Fig. 5 we show the projected column density of the 3D shell along three viewing axes in galactocentric Cartesian coordinates (x, y, z) at a time of 58 kyr after the formation of the thin shell. The size of the box is 20 × 20 × 20 pc and the thin shell expands into the homogeneous ambient medium of density n_out = 10⁴ cm⁻³ from SN explosions at r_GC = 5 pc, θ_GC = 0°, 70°, 90°.

As shown in Fig. 5 (top panel), at θ_GC = 0°, where the galactocentric gravitational attraction is balanced with the centrifugal force, the expanding shell is deformed by differential rotation, creating elongated shapes when projected onto the galactic plane, wrapping around the Galactic center. In this case, at 200 kyr, when all the shell elements become subsonic, m_{sh, tot} = 6300 M_⊙.

Another deformation is due to the gravitational attraction, F_z, parallel to the rotational axis, which is an important part of the total gravitational force not balanced by the rotation, stretching the shell toward the symmetry plane of the galaxy. For θ_GC = 70° (Fig. 5, middle panel), the shell elements follow helical orbits with parts of the shell spiraling toward the center of gravity. In this case, its total mass is m_{sh, tot} = 15 000 M_⊙.

In the case in which θ_GC = 90° (Fig. 5, bottom panel), when the SN explodes on the rotational axis, there is almost no rotation and the expansion is influenced by the galactocentric gravitational force only. Its total mass is m_{sh, tot} = 10 500 M_⊙.

The mass of the shell corresponds to the volume originally occupied by the ISM that was collected in the shell elements. It is influenced by the length of the path on which an element has moved supersonically relative to the motion of the ISM. The longest supersonic path is at middle elevations around θ_GC ~ 70°, where the rotation combines with the motion towards the galactic plane (see the next section for more detailed analysis).

7 Results

7.1 Angular momentum of the swept up gas

We explore the distribution of the angular momentum, $L{(R,z)}_i$, of shell elements at the time when their expansion velocity drops below the speed of sound. At this time, the mass accumulation is stopped and the mass m_i of an element and its angular momentum L_i are fixed. The angular momentum, $L{(R,z)}_i$, of an individual shell element defines the “angular momentum radius”, R_L, which is the distance from the rotational axis where a body at the circular orbit at z = 0 has the same specific angular momentum (angular momentum per unit mass): $$R_L{V}_{\text{rot}}(R_L,z=0)=L{(R,z)}_i/m_i.$$(16)

The distribution of mass as a function of angular momentum radius R_L represents the mass distribution in a flat disk, where the two components of velocity, the component away from the rotational axis and the component parallel to the rotational axis, are suppressed. The only remaining component is the circular rotational component. Actually, R_L reflects the centrifugal barrier of the matter collected in expanding SN shells.

We compare the distribution of mass as a function of R_L : in Fig. 6 we plot the fraction of mass Δm distributed in ΔR_L bins as a function of R_L in the following cases:

• 1.

  ISM: it gives the mass of a homogeneous ISM with rotation given by Eq. (8) at galactocentric distances r_GC < 50 pc.

• 2.

  SN-ISM: it shows the total mass of elements of expanding shells at the time of crossing the sound speed from SNe distributed randomly within the computational region following a uniform density distribution. The total mass and angular momentum depend on the total number of SNe included in the computation. Here, the total number of SNe is adjusted so that the total mass in all expanding SNe elements (at the time when they stop mass accumulation) is the same as the mass of the homogeneous medium inside the computational region that was used in the previous case, taking into account the density of the ISM.

• 3.

  SN-NSC: it shows the total mass of elements of expanding shells at the time of crossing the sound speed from SNe distributed with a radial profile corresponding to the density distribution in the young NSC₁ as described by Eq. (3). The total number of SNe is the same as in the case SN-ISM.

These three cases are compared in Fig. 6 where we show the mass, Δm, in ΔR_L bins as a function of R_L. The mass distribution as a function of R_L for the homogeneous ISM shows a narrow central peak composed of the matter residing along the rotational axis with L_i = 0, and some more ISM elements with small L_i, since the L_i of elements of the rotating ISM at a given R decreases with z proportional to $R/\sqrt{R^2+z^2}$. Consequently the region of small L_i < R_LV_rot(R_L, z = 0) has a conical shape around the rotational axis. This region, where the mass with a small angular momentum resides, is separated from the regions with a mass of higher angular momentum by a separation line at $$\frac{V_{\text{rot}}(R,z)R}{V_{\text{rot}}(R_L,z=0)}=R_L.$$(17)

The homogeneous ISM shows an increasing profile of $\frac{\Delta m}{\Delta R_L}(R_L)$ (see Fig. 6); its slope reflects the shape of the rotation curve as the rotation slows down with distance z off the symmetry plane. The distribution has a maximum near ~35 pc and it decreases for R > 40 pc, which is due to the finite size of the computational region (in our case 50 pc from the galactic center). The three densities of the ambient medium are compared in Fig. 6. In the case of the homogeneous ISM distribution, they show identical profiles (providing they are normalized to the same total mass).

In the second case, with the homogeneous distribution of SNe within the computational region, for n_out = 10⁴ cm⁻³ there is an increase up to the maximum at ~10 pc followed by a decrease for R_L > 10 pc (see Fig. 6). With n_out = 10⁵ cm⁻³, the maximum is at a distance of ~7 pc and it has a higher value, and with n_out = 10³ cm⁻³ the profile is much flatter, with a maximum at ~17 pc at the lower maximum value.

Compared to the previous case of a homogeneous ISM distribution, some of the matter with R_L > 25 pc is shifted to lower R_L values forming the maximum at lower values of R_L. This is due to the fact that the motion of those elements of the expanding SNe shells at high z distances from the galactic plane with low angular momentum is accelerated by the part of the galactocentric gravitational attraction of the SMBH and NSC not balanced by rotation. This force, which is perpendicular the galactic plane, is the reason why those elements moving toward the galactic plane collect more mass than they would collect without this perpendicular force. The motion of elements at low z distances with high angular momentum is not influenced by this perpendicular attraction (or less compared to the high z cases) and the mass and momentum accumulation is stopped earlier.

In the third case, the total number of SNe is the same as in the case of the homogeneous distribution. As one may see in Fig. 6, there is less mass particularly in the outer parts at R_L > 5 pc, and the maxima shift below R_L ~ 5 pc, close to the size of the core of the density profile of the young NSC₁, where most of the SNe explode.

Fig. 5 Three-dimensional simulations with the code RING. We show the projected shell column density at 60 kyr after the formation of the shell. The SN explodes at rGC = 5 pc and θGC = 0° (top), 70° (middle), 90° (bottom) in an ambient medium of density nout = 104 cm−3, with MSMBH = 4 × 106 M⊙ and the NSC as given in Sect. 2. The color bar gives the projected density in units of cm−2. The meshed region shows the central region with a radius of 1 pc. The dashed line shows the reference circle with a radius of 5 pc and the cross gives the original position of the SN explosion.

Fig. 6 Distribution of mass as a function of RL. We show a fraction of the total mass Δm in the ΔRL bins as a function of RL for three values of nout, with MSMBH = 4 × 106 M⊙ and the NSC as given in Sect. 2. The number in the line label shows the power index of the value nout = 103, 104, and 105 cm−3. The vertical scale holds for nout = 104 cm−3; for a lower or higher ISM density it should be changed to ten times lower or higher values, respectively.

Fig. 7 Border lines for feeding the central 1 pc in the (R, z) plane. Black continuous line: RL = 1 pc line separating the low and high angular momentum places. Colors represent different values of nout in particles per cm3: blue – 103; black – 104; red – 105, with MSMBH = 4 × 106M⊙ and the NSC as given in Sect. 2. The open symbols show feeding of the central 1 pc by SNe shell elements, solid symbols show feeding by supersonic elements only. The thin dashed lines give places of constant rGC (in 5 pc steps) and θGC (in 10° steps).

7.2 Feeding the central region

We derive how much of the shell mass is delivered into the central region, defined as a sphere with radius, r_feed, which is a free parameter in this model. The physical motivation behind it is that SMBHs are often surrounded by gaseous circumnuclear disks, and any gas moving through the volume occupied by such a disk will probably impact it and merge with it in dissipative shocks and then eventually will be accreted onto the SMBH. In the Milky Way, the radius of the circumnuclear disk is on the order of several parsecs (Tsuboi et al. 2018). Here we stay conservative and set r_feed = 1 pc.

The mass feeding depends on the position of the SN explosion. Figure 7 shows the line in the (R, z) plane where $$\frac{V_{\text{rot}}(R,z)R}{V_{\text{rot}}(1,0)}=1\text{pc}$$(18)

(black continuous line) separating the places of low angular momentum (to the left of it) from high angular momentum places (to the right from it). It is compared to lines based on RING simulations separating the places of the SNe explosions from where the shell elements reach the central 1 pc (to the left of the zigzag curve), and places from where the shell elements do not reach the central 1 pc (to the right of the zigzag curve). We show the results for different values of n_out and we distinguish the mass delivered by supersonic elements (solid lines) from mass that is delivered to the central 1 pc by all elements (dashed lines).

The zigzag lines in Fig. 7 generally follow the solid line separating the low and high angular momentum places, but some feeding of the central 1 pc can also be achieved even from high angular momentum places. This is more likely close to the galactic center at r_GC < 20 pc for expansions in the media having lower values of n_out. It can be explained as being due to the larger size of shells expanding in media with a lower density compared to shells expanding in media with a higher density. At larger distances from the galactic center at r_GC > 20 pc, the size of the shells is small compared to the galactocentric distance for any of the tested densities and the delivery to the central parsec follows the line separating low and high angular momentum places.

There is another result visible in Fig. 7. The mass feeding into the central 1 pc with supersonically expanding shell elements is restricted by the initial galactocentric distance, r_GC: it happens from SNe inside of a maximum r_GC,feed,max, which depends on the density of the ambient medium, n_out. It is 20, 50, or 140 pc for n_out = 10⁵, 10⁴, or 10³ cm⁻³; for larger distances no supersonic elements reach the central 1 pc.

We examine the mass m_{sh, tot}, which is the total mass collected in a single expanding SN shell at the end of mass accumulation when all elements move with subsonic velocity. In Fig. 8, we show m_{sh, tot} at different (r_GC,  θ_GC) positions together with the fraction that is delivered by the supersonically moving elements of a shell into the central 1 pc. We show the results for three different values of n_out. The largest values of m_{sh, tot} are achieved at certain values of θ_GC, which grow with decreasing n_out from 70° to 90°. The maximum value of m_{sh, tot} increases from 10⁵ M_⊙ for n_out = 10⁵ cm⁻³ to 10⁶ M_⊙ for n_out = 10³ cm⁻³, and its position shifts away from the galactic center for lower values of n_out.

As can be seen in Fig. 8, there is a delivery region from where the supersonic elements reach the central 1 pc. We compute the average mass fed per SN, $\overline{m_{\text{feed}}^{100}}$, which is the average mass transported by one SN explosion taking place in the region r_GC ≤ 100 pc to the central region r_GC ≤ r_feed = 1 pc. SNe are positioned randomly in the regionr_GC ≤ 100 pc with a specified radial distribution (homogeneous, NSC₁). The product of the SN rate within the central 100 pc and the average mass fed per supernova, $\overline{m_{\text{feed}}^{100}}$, gives the mass flux ${\dot{M}}_{\text{feed}}^{100}$, which is the amount of mass transported per unit of time from the region within 100 pc into the central 1 pc via expanding SN shells.

The values of $\overline{m_{\text{feed}}^{100}}$ are given in Table 2. There we also give the fraction of the volume inside 100 pc occupied by an average expanding SN shell in supersonic expansion, $\overline{f_{\text{SN}}^{100}}$, the fraction of volume inside 100 pc occupied by the delivery region, $f_{\text{delivery}}^{100}$, and an average time interval during which the shell is at least partly supersonic, $\overline{t_{\text{sonic}}^{100}}$. The last quantity multiplied by the SN rate within the central 100 pc gives the average number of SNe that should be visible within the central 100 pc.

Fig. 8 The total mass collected in a single expanding shell at the end of mass accumulation, msh, tot, as a functionof the galactocentric distance rGC and the elevation angle θGC for nout = 103, 104, and 105 cm−3, with MSMBH = 4 × 106M⊙ and the NSC as given in Sect. 2. The mass delivery to the central parsec is shown with hatching: the areas without hatching do not deliver any mass, single hatched areas deliver up to 10% of msh, tot, the double-hatched areas deliver more than 10% of msh, tot.

8 Discussion

Observations of Seyfert galaxies and AGNs (Diamond-Stanic & Rieke 2012) show a link between the nuclear star formation rate, SFR, and the black hole accretion rate, Ṁ_SMBH, that is approximated as $\text{SFR}\propto {\dot{M}}_{\text{SMBH}}^{0.8}$. AGNs at z < 0.35 exhibit an SFR correlated with X-ray luminosity, demonstrating the presence of AGN activity (Zhuang & Ho 2020). The correlation between SFR and Ṁ_SMBH was discussed by Zhuang et al. (2020, and references therein): both quantities depend on the molecular gas mass within the central region, raising the question of how significant the above correlation is. However, the AGN activity evolves on much shorter timescales, < 100 Myr, than the star formation, >100 Myr, and rapid AGN variability occurs on even shorter timescales (Hickox et al. 2014), thus presenting a challenge for models connecting Ṁ_SMBH to the SFR.

Here, we propose a connection between SFR and Ṁ_SMBH based on expanding SN shells. The SN rate is certainly connected to the SFR: with a Salpeter initial mass function, or with a top-heavy IMF, the SN rate is in any case proportional to the SFR. Within the central region, we do not separate the clumpy star formation inisolated molecular clouds from more distributed star formation inside the entire region: we assume that the star formation happens in the central region with a homogeneous or centrally condensed distribution similar to the scale of the NSC₁.

Observations of the Galactic center show the SFR = 3 × 10⁻³ M_⊙ yr⁻¹ (Figer et al. 2004; Rockefeller et al. 2004; Zubovas et al. 2013), which corresponds to a SN rate of 3 × 10⁻⁵ SN yr⁻¹. Similarly, a burst forming 10⁵ M_⊙ of stars within 100 pc of the galactic center will produce ~1000 massive stars that will explode during the following 30 Myrs as SNe, which also gives an average SFR of 3 × 10⁻³ M_⊙ yr⁻¹, and a corresponding SN rate within the central 100 pc of 3 × 10⁻⁵ yr⁻¹ during those 30 Myr. For homogeneously distributed star formation in an ISM of density 10³ −10⁵ cm⁻³, the resulting supernovae would feed a total mass into the central 1 pc of (9 × 10³−10⁵) M_⊙ with an inflow rate of $\stackrel{.}{\text{M}}=(3\times {10}^{-4}-4\times {10}^{-3})$ M_⊙ yr⁻¹. If star formation were distributed in the same way as the NSC₁ density distribution, the total mass fed to the central parsec would be (2 × 10⁵−10⁶) M_⊙, with an inflow rate of $\stackrel{.}{\text{M}}=(7\times {10}^{-3}-4\times {10}^{-2})$ M_⊙ yr⁻¹.

The SN rate within the central 100 pc – 3 × 10⁻⁵ yr⁻¹ may be multiplied by the average duration of a SN shell as given in Table 2, to give the number of expanding shells we would expect to see after a burst of star formation. This number varies between 6 and 35 for different ISM densities and for alternative choices of NSC₁ or homogeneous stellar distributions. The expanding SNe shells only fill a small fraction of the volume.

We show that an expanding SN shell can deliver considerable mass into the galactic center during epochs of star formation. The bursting AGN activity may be connected to mass delivery into the circumnuclear disk by individual expanding SN shells. Probably the SN shells feed the circumnuclear disk, which provides the mass to the central accretion disk. However, this accretion is limited by stellar winds and outflows due to feedback from the SMBH and young stars in the central star cluster. Thus, part of the inflowing mass may be pushed away, never reaching the nuclear SMBH.

However, AGNs show much higher values (Zhuang & Ho 2020): SFR = 1−100 M_⊙ yr⁻¹. With such a large nuclear SFR, the mass inflow via expanding shells may increase to 0.1−1000 M_⊙ yr⁻¹. This is just a crude estimate, since the shell versus shell interactions and supersonic shell collisions will influence the results and their importance should be evaluated.

The mass delivery by expanding SN shells also depends on the size of the circumnuclear disk. It can be only a fraction of a parsec as in NGC 4258 (Menezes et al. 2018) or it can be very large, ~ 100 pc, as in gas-rich AGNs (Combes et al. 2019; Combes 2020). We explored the dependence of the mass delivery rate on the size of the delivery region for M_SMBH = 4 × 10⁶M_⊙: the total mass delivered scales with square of the radius of the circumnuclear disk, which probably reflects the geometry of the “capture region”. The mass delivery by expanding shells as a function of the size of the sphere of influence of the SMBH will be analyzed in a separate paper.

A studyof SN remnants in the vicinity of a SMBH was also described by Rimoldi et al. (2015), with an emphasis on their X-ray properties. They adopt a similar method using the Kompaneets (1960) approximation that assumes the wall of the expanding SN shell is thin. This assumption is valid when the thin expanding shell forms, which is connected to the rapid loss of thermal energy due to the increased density leading to enhanced cooling. Rimoldi et al. (2015) use the thin shell approximation to describe also the earlier Sedov-Taylor phases when the expansion is still adiabatic, where the Kompaneetsapproximation may not be accurate. However, there are several differences compared with our treatment. In Rimoldi et al. (2015): (1) the ISM density distribution has a steep radial gradient, in particular outside 1 pc from the galactic center; (2) the gravity includes the SMBH only, the NSC is not considered; (3) the distribution of massive stars and corresponding SNe is more concentrated in the galactic center, even when compared to the NSC₁ distribution adopted here. Also, they do not account for the angular momentum of the ISM, which is critical to our treatment. Rimoldi et al. (2015) show SN remnant morphologies for different masses of the SMBH and for explosions at different distances from the SMBH. Here, we perform a more statistical study aimed at elucidating the inflow rates as a function of ISM density, SNe radial distribution, and SNe location above the galactic plane.

Supernovae expanding shells with a low angular momentum that occur in a conical region near a rotational axis of a galaxy deliver mass to the central region, where it may settle into a circumnuclear disk. When shell material passes within 1 pc of the central black hole, its orbital time,10⁵ yr $\times ({\frac{{10}^6 M_{\odot }}{M_{\text{BH}}})}^{1/2}$, approaches the expansion time of the supernova shell, and the specific angular momentum of the shell material will cause it to cease its inward motion at the centrifugal barrier. At that point, the shearing shell will soon circularize, creating a circumnuclear disk.

The SMBH is fed out of the disk, however its viscous time may be much greater compared to the disk feeding time by expanding SNe shells. Bursts of star formation invoke the scenario where the circumnuclear disk is formed quickly within a few tens of millions of years followed by longer periods of dissolution and continuous SMBH feeding out of the circumnuclear disk.

An alternative idealized physical situation could be one in which a pressure gradient supports the vertical stratification of the gas. However, there is little evidence that the gas layer in galactic nuclei follows such a simple description, except on relatively small scales (several parsecs) where circumnuclear disks can be present. A more realistic description of the interstellar medium would be one in which there are pronounced density inhomogeneities and in which the vertical structure of the gas on ~ 10−20 pc scales is determined by a stochastic velocity field caused by transient flows and ejections caused by supernovae and winds from star formation complexes, by accretion events, and by hydromagnetic activity. Our simplified treatment of the interstellar medium is adopted in order to identify the important dynamical implications of a supernova occurring in the strongly anisotropic gravitational potential of a galactic nucleus. In later papers, we will investigate how specific gas configurations such as circumnuclear disks or a large-scale vertically stratified medium would affect the dynamics of a supernova remnant. Ultimately, the consequences of a relatively strong, ordered magnetic field will also be explored.

If the central parsec contains massive stars with strong winds, as our Galactic center does, or if accretion onto the SMBH leads to an outflowing wind, then the interaction of such winds with an incoming supernova shell will alter the dynamics, depending on the strength of the wind, but it would not be likely to significantly change the amount of ultimately accreted matter. This is because the momentum of such winds is predominantly radial, so it does not affect the specific angular momentum of the material with which it interacts, except possibly to reduce it (Blank et al. 2016). We have not included such complications in our models, but note that in specific cases such as the Galactic center, it can embellish the evolutionary picture (see Rockefeller et al. 2005).

Another complication that might be presented by a galactic nuclear region is that, during the expansion time of the models that we have considered, there is a fairly significant chance that other SNe will occur in the overall volume that we are considering (for example, a glance at the radio image published by Heywood et al. 2019 reveals a relatively high projected surface density of SN remnants that can be seen toward the Galactic center). So while we only consider the dynamics induced by a single SN, the more general case of multiple SNe might be occurring a significant fraction of the time, and could have consequences such as strong shocks colliding with the shell, leading to some redistribution of angular momentum that might favor mass flow into the central parsec, or reduce the mass delivery predicted by the isolated SN model.

If a supernova takes place inside the central parsec, it might not have time to create much of a shell, depending on the ambient density (cf., Fig. 3), in which case the blast wave would pass over the SMBH, perhaps diverted to some extent by the SMBH magnetosphere or by any wind produced by the SMBH, but relatively little mass would be deposited there to be accreted. Indeed, a SN occurring within the central parsec could have the effect of partially or fully clearing that region, depending on how much interstellar gas is present there. The relevance of SN in the central parsec is highlighted by the fact that a massive young cluster with ~100 O stars or massive post-main-sequence stars occupies the central half-parsec of our Galaxy, and furthermore, the overall stellar density is maximized in the central parsec, so that even Type Ia SNe would be concentrated there if their progenitors are present. In the case of the Galactic center, the existing circumnuclear disk is too massive to be expelled by a SN, but in any case, its inner edge has a ~1 parsec radius, so recent SNe might be partially responsible for the relatively evacuated central cavity.

The very central part at r_GC < 0.01 pc will be strongly influenced by the activity of the SMBH. The timescale for further dynamical evolution of the material inside r_GC = 1 pc will then be the relatively long viscous time of the disk, and the fate of much of the disk will be to migrate inward and eventually lead to direct accretion onto the SMBH plus a concomitant outflow of material energized by the accretion.

9 Conclusions

We performed simulations of SN shells expanding inside the NSC near the SMBH with the code RING using the infinitesimally thin shell approximation in 3D. In our model, the homogeneously distributed ISM rotates following the rotation curve (Eq. (8)), which partially balances the gravitational attraction of the NSC and SMBH. We tested three different ISM densities, n_out = 10³, 10⁴, and 10⁵ cm⁻³, and identified the delivery region from where the SN expanding shells deliver mass into the central 1 pc. It is a region of small angular momentum with a low value of R_L centered around the rotational axis (see Eq. (16)).

The simulations show that SNe occurring within a conical region around the rotational axis of the galaxy can feed the central accretion disk surrounding the SMBH. A burst forming 10⁵ M_⊙ of stars within the central 100 pc will produce ~ 1000 SNe during the subsequent 30 Myr, which will deliver with their expanding shells a total mass of (1− 2) × 10⁵ M_⊙ to the central parsec, for ISM density n_out = 10³ cm⁻³, or (9 × 10³−10⁶) M_⊙ for ISM density n_out = 10⁵ cm⁻³. The corresponding continuous SFR = 0.003 M_⊙ yr⁻¹ leads to an average mass flux into the central parsec of (3− 7) × 10⁻³ or (0.3−40) × 10⁻³ M_⊙ yr⁻¹. A high SFR ~ 30 M_⊙ yr⁻¹ inside the central 100 pc of the galaxy will increase this mass inflow rate to ≤ 3−400 M_⊙ yr⁻¹ depending on the ISM density and on the radial distribution of SNe.

Acknowledgements

We are grateful to the anonymous referee whose constructive comments improved the manuscript. This study has been supported by Czech Science Foundation Grant 19-15480S and by the project RVO:67985815.

References

All Tables

All Figures

	Fig. 1 Interstellar matter (ISM) rotation velocity, Vrot, is shown as afunction of R for different distances z from the symmetry plane.
In the text

Fig. 2 One-dimensional FLASH simulations. Left panel: time evolution of the position of the forward shock, rsh (r3, r4, r5), and of its expansion velocity, vsh (v3, v4, v5). Right panel: time evolution of the thermal energy inside the shock, Eth (Et3, Et4, Et5), and of the kinetic energy of the shock, Ekin (Ek3, Ek4, Ek5). Individual quantities are labeled with the value of the power index of the ISM density nout = 103, 104, and 105 cm−3.

In the text

	Fig. 3 One-dimensional FLASH simulations: profiles of the density n, the temperature T, the velocity v, the sound speed cs and the pressure p∕kB for r < 1 pc at time = 502 yr in a medium of density nout = 104 cm−3.
In the text

Fig. 4 Shell radius, rsh, as a function of time during the first 3000 yr of expansion. The 1D FLASH is joined by the 3D RING at tsh : absolute values (left panel), relative values of the difference (rsh,F − rsh,R)∕rsh,F (right panel), where rsh,F is the radius of the shell according to FLASH and rsh,R is the radiusof the shell according to RING simulations. The number in the line label shows the power index of the value nout = 103, 104, and 105 cm−3.

In the text

Fig. 5 Three-dimensional simulations with the code RING. We show the projected shell column density at 60 kyr after the formation of the shell. The SN explodes at rGC = 5 pc and θGC = 0° (top), 70° (middle), 90° (bottom) in an ambient medium of density nout = 104 cm−3, with MSMBH = 4 × 106 M⊙ and the NSC as given in Sect. 2. The color bar gives the projected density in units of cm−2. The meshed region shows the central region with a radius of 1 pc. The dashed line shows the reference circle with a radius of 5 pc and the cross gives the original position of the SN explosion.

In the text

Fig. 6 Distribution of mass as a function of RL. We show a fraction of the total mass Δm in the ΔRL bins as a function of RL for three values of nout, with MSMBH = 4 × 106 M⊙ and the NSC as given in Sect. 2. The number in the line label shows the power index of the value nout = 103, 104, and 105 cm−3. The vertical scale holds for nout = 104 cm−3; for a lower or higher ISM density it should be changed to ten times lower or higher values, respectively.

In the text

Fig. 7 Border lines for feeding the central 1 pc in the (R, z) plane. Black continuous line: RL = 1 pc line separating the low and high angular momentum places. Colors represent different values of nout in particles per cm3: blue – 103; black – 104; red – 105, with MSMBH = 4 × 106M⊙ and the NSC as given in Sect. 2. The open symbols show feeding of the central 1 pc by SNe shell elements, solid symbols show feeding by supersonic elements only. The thin dashed lines give places of constant rGC (in 5 pc steps) and θGC (in 10° steps).

In the text

Fig. 8 The total mass collected in a single expanding shell at the end of mass accumulation, msh, tot, as a functionof the galactocentric distance rGC and the elevation angle θGC for nout = 103, 104, and 105 cm−3, with MSMBH = 4 × 106M⊙ and the NSC as given in Sect. 2. The mass delivery to the central parsec is shown with hatching: the areas without hatching do not deliver any mass, single hatched areas deliver up to 10% of msh, tot, the double-hatched areas deliver more than 10% of msh, tot.

In the text