Nuclear Spiral Shocks and Induced Gas Inflows in Weak Oval Potentials

TL;DR

This study uses high-resolution hydrodynamic simulations to analyze how weak oval potentials in galaxy centers can generate nuclear spirals and shocks, leading to gas inflows that may fuel black hole activity.

Contribution

It demonstrates that even weak bar-like potentials can produce significant nuclear spirals and shocks, influencing gas inflow rates in galaxy centers.

Findings

Weak potentials can induce shocks at small radii.

Low shear results in loosely wound spirals and higher inflow.

High shear produces tightly wound spirals and circumnuclear disks.

Abstract

Nuclear spirals are ubiquitous in galaxy centers. They exist not only in strong barred galaxies but also in galaxies without noticeable bars. We use high-resolution hydrodynamic simulations to study the properties of nuclear gas spirals driven by weak bar-like and oval potentials. The amplitude of the spirals increases toward the center by a geometric effect, readily developing into shocks at small radii even for very weak potentials. The shape of the spirals and shocks depends rather sensitively on the background shear. When shear is low, the nuclear spirals are loosely wound and the shocks are almost straight, resulting in large mass inflows toward the center. When shear is high, on the other hand, the spirals are tightly wound and the shocks are oblique, forming a circumnuclear disk through which gas flows inward at a relatively lower rate. The induced mass inflow rates are enough to power black hole accretion in various types of Seyfert galaxies as well as to drive supersonic turbulence at small radii.

Figures

Equations

$$V=65+95\tanh \left(\frac{R-0.07}{0.06}\right)-50\log R+1.5(\log R+3{)}^3,$$

$$V^2=\frac{G{M}_{\mathrm{BH}}}{R}+{\left(\frac{V_{0}R}{R_0+R}\right)}^2,$$

$$\begin{array}{r}\lambda =\frac{L_{\mathrm{bol}}}{L_{\mathrm{Edd}}}=1.5\times 10^{-2}f\left(\frac{ϵ}{0.1}\right)\left(\frac{{\Phi }_a{V}_0^{-2}}{2\times 10^{-4}}\right)\\ \left(\frac{{\Sigma }_0}{10^2{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}}\right){\left(\frac{M_{\mathrm{BH}}}{3\times 10^6{\mathrm{M}}_{\odot }}\right)}^{-1},\end{array}$$

$$\sigma \sim \frac{(\dot{M}/\Sigma )V^2}{G{\pi }^2\delta R^2\Sigma }.$$

$$\sigma \sim 230\frac{V_{100}^2}{\delta R_{10}^2{\Sigma }_{100}}\mathrm{km}{\mathrm{s}}^{-1},$$

Full text

††slugcomment: Accepted for the publication in the ApJL

1. Introduction

Nuclear dust and gas spirals are observed on scales of several hundred parsecs in the central parts of the Milky Way (Roberts & Goss, 1993; Sofue, 1995; Sawada et al., 2004; Hsieh et al., 2015; Ridley et al., 2017) and other disk galaxies (Buta & Crocker, 1993; Knapen et al., 1995). Sometimes these spirals have irregular shapes (Elmegreen et al., 1998; Martini et al., 2003a, b) reminiscent of turbulence (Elmegreen et al., 2002). Figure 1 shows two examples of irregular spirals. The more regular cases have been attributed to gas flows in bars or nuclear bars (Simkin et al., 1980; Shlosman et al., 1989b; Martini et al., 2003b), while the irregular spirals could be from gravitational (Shlosman et al., 1989a), magnetic (Balbus & Hawley, 1991), or hydrodynamic (Montenegro et al., 1999; Krumholz & Kruijssen, 2015) instabilities. Since many nuclear disks have Toomre stability parameters much larger than unity (Elmegreen et al., 1998; Martini & Pogge, 1999; Lin et al., 2016), the spirals in these cases are unlikely to form by two-dimensional (2D) gravitational processes. This leaves the question of how nuclear turbulence might be generated and how it can become nonlinear to the point where shocks form, giving the appearance of sharply delineated, irregular, dust lanes.

Here we consider the possibility that turbulent wavelets moving inward will strengthen by geometric convergence, as predicted by the $1/R$ terms for radius $R$ in the nonlinear equations for spiral density waves (Bertin et al., 1989). Montenegro et al. (1999) called this a curvature instability although it is more of an amplification of pre-existing waves than an exponential growth from noise. We simulate 2D nuclear spirals driven by oval distortions of various strengths and rates of shear in the background rotation curve. Even with a weak distortion, the resulting gas structures readily form shocks at small radii, as predicted by the analytic theory. These shocks then torque the gas and cause it to accrete at a fairly high rate. We suggest that in three-dimensional models with turbulent energy dissipation, the gravitational binding energy released by this accretion will sustain the turbulence.

There have been many previous studies of bar-driven nuclear spirals starting with Sanders & Huntley (1976), Roberts et al. (1979), Athanassoula (1992), Piner et al. (1995), and others. Simulations by Maciejewski (2004), Thakur et al. (2009) and others traced bar-driven spiral shocks all the way to the vicinity of a nuclear black hole. These shocks were driven by strong bars and formed the usual long dust-lanes that are observed in barred galaxies. Kim et al. (2012a, b) and Li et al. (2015) also had a strong bar to study the formation of nuclear rings. Strong forcings like these overwhelm the geometric effects of shock steepening during inward wave propagation. Other works also did not consider our low-shear case, which brings out the curvature steepening even more as the resulting spiral pitch angles are large.

2. Models

We consider bar- and oval-driven spirals in two models with different background shear: HS and LS models. The HS model is to represent the central regions of galaxies with high shear like the Milky Way, whose mass distribution is dominated by a nuclear bulge and a supermassive black hole (SMBH) (Launhardt et al., 2002). The resulting rotation velocity (e.g., Krumholz & Kruijssen 2015) can be well fitted by

[TABLE]

for $\;0.01\leq R\leq 10$, where $V$ is in units of km s*-1* and $R$ in kpc. On the other hand, the LS model is designed to simulate disk galaxies with low shear, like NGC 3041 (Erroz-Ferrer et al., 2016), for which we take

[TABLE]

where $V_{0}=220\rm\;km\;s^{-1}$, $R_{0}=0.3{\rm\,kpc}$, and a black hole mass of $M_{\rm BH}=3\times 10^{6}\;M_{\odot}$.

Figure 2 plots the radial distributions of the rotational velocity, shear parameter $q\equiv-d\ln\Omega/d\ln R$ with angular frequency $\Omega=V/R$, and various frequencies of the HS and LS models. The gas is taken to be infinitesimally thin, initially uniform with surface density $\Sigma_{0}$, and isothermal with a constant speed of sound ${c_{s}}$. We take ${c_{s}}=10{\rm\,km\,s^{-1}}$ as our fiducial value, but also consider the cases with ${c_{s}}=20{\rm\,km\,s^{-1}}$. For the perturbing gravitational potential, we adopt a simple bi-symmetric sinusoidal form $\Phi_{\rm ptb}(R,\phi,t)=\Phi_{a}\cos(2\phi-2\Omega_{p}t)$ with constant amplitude $\Phi_{a}$. We consider two cases with $\Phi_{a}/V_{0}^{2}=2\times 10^{-5}$ for an oval distortion in an unbarred galaxy and $\Phi_{a}/V_{0}^{2}=2\times 10^{-4}$ for a weak bar in a barred galaxy. These may also represent a weak secondary bar or a more-or-less isotropic central portion of a strong bar in more strongly barred galaxies. By considering a weak bar, our models cannot capture interactions of nuclear spirals with a nuclear ring that forms by a strong bar.

The potential rotates rigidly with pattern speed $\Omega_{p}=60$ and $30{{\rm\,km\,s^{-1}}{\rm\,kpc}^{-1}}$ for the HS and LS models, respectively. We use Athena++, a newly developed grid-based code utilizing a higher-order Godunov scheme (Stone et al. 2017, in preparation), to evolve the gas in 2D logarithmic cylindrical coordinates (Kim et al., 2012b). Our simulation domain extends from the inner radial boundary, $R_{\rm in}=10{\rm\,pc}$, to the corotation radius (CR), and covers $\phi=0-2\pi$. The number of zones in the radial and azimuthal directions is $1024\times 516$ and $1024\times 548$ in the HS and LS models, respectively, making all the zones approximately square-shaped.

The imposed gravitational potential perturbs an otherwise uniform gas disk. The perturbations organize into spiral waves that grow as thermal pressure tends to align the apocenters of perturbed gas orbits inside the inner Lindblad resonance (ILR), reinforcing the perturbations (Montenegro et al., 1999). This is opposite to the case of spirals outside the ILR where gravity aligns the perturbed orbits and pressure resists them.

3. Results

Figure 3 plots the snapshots of the gas surface density $\Sigma$ at selected times for the HS (left) and LS (right) models with $\Phi_{a}/V_{0}^{2}=2\times 10^{-4}$. Figure 4 displays the distributions of $\Sigma$ at $t=1.0{\rm\,Gyr}$ in the $x$–$y$ plane, with successive zoom-in views. The nuclear spirals are limited to inside the CR, with very weak density variations outside the CR. They are tightly wound in the HS model due to strong background shear, and loosely wound in the LS model. The spirals are piecewise logarithmic, with a pitch angle of $i_{p}\sim 10^{\circ}$ at $R>0.15{\rm\,kpc}$ in the HS model and $i_{p}\sim 35^{\circ}$ at $R>0.3{\rm\,kpc}$ in the LS model.

Initially, the shape of the nuclear spirals is controlled by the radial gradient of $\Theta\equiv\Omega-\kappa/2$ such that they are leading in the regions where $d\Theta/dR>0$ and trailing where $d\Theta/dR<0$ (Buta & Combes, 1996; Kim et al., 2012b), a characteristics of kinematic density waves (Wada, 1994; Sormani et al., 2015). Favored by shear, trailing spirals usually grow faster and more strongly than leading spirals (Kim et al., 2012b). In addition, the regions (between $R_{1}$ and $R_{2}$ in Fig. 3) for leading spirals are so narrow that the density inside the ILR becomes soon dominated by trailing spirals. Due to a geometric effect (Montenegro et al., 1999), trailing spirals increase in amplitude as they converge to the origin, analogous to the amplification of sound waves in sonoluminescence (Kondic et al., 1995). Thus, the nuclear spirals become nonlinear earlier at smaller $R$, readily developing into shock waves. Strong background shear makes the shock fronts inclined to the gas streamlines in the HS model, resulting in relatively low perpendicular Mach numbers of $\mathcal{M}_{\perp}\sim 1.5$. In the LS model, the shock fronts are less inclined and become stronger as they unwind over time, resulting in $\mathcal{M}_{\perp}\sim 2.5$. This unwinding is presumably caused by a nonlinear increase in the angular momentum flux (Lee & Goodman, 1999). We find that nuclear spirals in models with ${c_{s}}=20{\rm\,km\,s^{-1}}$ evolve qualitatively similarly to those in the lower ${c_{s}}$ counterparts, although they form earlier in the former due to stronger thermal pressure.

The formation of shocks has two effects: the resulting torques and orbital energy dissipation cause gas accretion to the center, while back reactions to the torques modify the trailing spirals. Because of the first effect, gas encountering the shocks loses its angular momentum and moves inward, eventually passing through the inner radial boundary of our simulation. At the same time, the regions with shocks grow in size radially as the spirals become stronger, reaching $R\sim 70{\rm\,pc}$ at $t=0.5{\rm\,Gyr}$ in the HS model and $R\sim 80{\rm\,pc}$ at $t=0.8{\rm\,Gyr}$ when the spirals saturate. The boundary of the shocked regions at $t=1{\rm\,Gyr}$ is indicated by the dotted circles in Figures 4c and 4f. In the HS model, the spiral shocks are oblique and in a logarithmic shape. The inflowing gas follows more-or-less circular orbits with small-amplitude epicycle motions, as indicated by instantaneous streamlines in Figure 3, and piles up near the center due to the geometric convergence effect, forming a circumnuclear disk with radius of $\sim 20-30{\rm\,pc}$. In the LS model, on the other hand, the shocks are almost perpendicular and gas after the shocks makes quite a radial orbit especially at $R\lesssim 25{\rm\,pc}$, plunging almost directly to the inner radial boundary, without forming a circumnuclear disk.

The second effect of the shocks changes the density and velocity fields significantly – not only in the shocked regions but also the surrounding regions. In particular, the velocity induced by the shocks is able to boost the epicycle amplitudes of gas elements just outside the shocked regions. The perturbations propagate outward in the form of pressure-modified inertial waves and interfere with the trailing spirals there. In the HS model, the perturbations are so weak that their effect is moderate, creating only mild modulation of the spirals. On the other hand, the perturbations in the LS model are quite strong and make the trailing spirals branch out, resulting in V-shaped structures, similar to those seen in NGC 2207 (Fig. 1). After the shocks saturate, gas just outside the shocked regions is pushed inward due to pressure gradients and experiences more shocks successively in orbital motions, moving toward the center. The infall radial velocity is typically $\sim 30{\rm\,km\,s^{-1}}$ and $\sim 50{\rm\,km\,s^{-1}}$ in the HS and LS models, respectively.

Figure 5 plots the mass inflow rate $\dot{M}(t)=-\displaystyle\int\Sigma v_{R}Rd\phi$ measured at $R=R_{\rm in}$ divided by $\Sigma_{0}$, where $v_{R}$ is the radial velocity, for all models with both ${c_{s}}=10$ and $20{\rm\,km\,s^{-1}}$. In the HS models, the nuclear spirals grow fast and initiate gas inflows when the shocks first form near the inner boundary. Since the spirals are tightly wound, however, the shocks are relatively weak and the induced mass inflow rate remains small in these models. In the LS models, on the other hand, spirals are more loosely wound and grow more slowly, resulting in stronger shocks and larger gas inflows. Note that the mass inflow rate in our models with a weak bar/oval is insensitive to ${c_{s}}$. This is in sharp contrast to the cases with a strong bar in which thermal pressure spreads out the gas in a nuclear ring to enhance the gas density at the center, resulting in larger $\dot{M}$ for larger ${c_{s}}$ (Ann & Thakur, 2005; Kim et al., 2012b). In all models, the mass inflow rate is almost proportional to the amplitude of the imposed external potential, and persists for a long period of time.

The time-averaged mass inflow rates in our models can be summarized as $\langle\dot{M}\rangle=10^{-5}f\Sigma_{0}(\Phi_{a}V_{0}^{-2}/2\times 10^{-4})$ $M_{\odot}$ yr*-1*, where $f$ is a factor, varying less than an order of magnitude, arising from the difference in the background shear, and $\Sigma_{0}$ is in ${{\rm\,M_{\odot}}{\rm\,pc}^{-2}}$. If all the inflowing mass is accreted to an SMBH with mass $M_{\rm BH}$, this corresponds to an Eddington ratio of

[TABLE]

where $L_{\rm bol}=\epsilon\langle\dot{M}\rangle c^{2}$ is the bolometric luminosity of an AGN with efficiency $\epsilon$, and $L_{\rm Edd}$ is the Eddington luminosity, $3.2\times 10^{4}\,(M_{\rm BH}/M_{\odot})\,L_{\odot}$. Observations show that $\lambda\lesssim 0.1$ for classical Seyfert 1 galaxies with broad emission lines (Meyer-Hofmeister & Meyer, 2011), $\lambda\sim 10^{-2}$ for Seyfert 2 galaxies with narrow lines (Bian & Gu, 2007), and $\lambda\sim 10^{-3}$ for low-luminosity Seyfert 1 galaxies (Ho, 2008). These values suggest that gas inflows driven by shock fronts in nonlinear nuclear spirals can account for observed levels of AGN activity in Seyfert galaxies, depending on the amount of gas in the disk and the amplitude of the external potential.

4. Turbulence Generation

Accretion at a rate ${\dot{M}}$ corresponds to a rate of energy input to the gas equal to about $\text@frac{1}{2}{\dot{M}}V^{2}$ for orbital speed $V$ (see Krumholz & Kruijssen, 2015, for a more detailed description). The energy dissipation rate from a uniform disk with surface density $\Sigma$, radius $R$, mass $M=\pi R^{2}\Sigma$, velocity dispersion $\sigma$, and scale height $H=\sigma^{2}/(\pi G\Sigma)$ is assumed to equal a factor $\delta<1$ times the turbulent energy divided by the crossing time over a disk thickness: $\text@frac{1}{2}\delta M\sigma^{3}/H=\text@frac{1}{2}G\pi^{2}\delta R^{2}\Sigma^{2}\sigma$. If we set the energy input rate equal to the energy dissipation rate and rearrange, we get the velocity dispersion in the inner disk,

[TABLE]

Figure 5 shows that ${\dot{M}}/\Sigma\sim 10^{-5}{\rm\,pc}^{2}{\rm\,yr}^{-1}$ for $\Phi_{a}/V_{0}^{2}=2\times 10^{-4}$. Inserting standard values, we obtain

[TABLE]

where $V_{100}=V/(100{\rm\,km\,s^{-1}})$, $R_{10}=R/(10{\rm\,pc})$, and $\Sigma_{100}=\Sigma/(100{{\rm\,M_{\odot}}{\rm\,pc}^{-2}})$. For radii up to tens of parsecs, $\sigma$ exceeds $\sim 20{\rm\,km\,s^{-1}}$. For realistic $\delta<1$, $\sigma$ is even higher.

This high value of $\sigma$ proposed for the three-dimensional case suggests that the accretion rate in our 2D simulations is high enough to pump supersonic turbulence at small radii. Such turbulence would presumably give a Kolmogorov power spectrum like that observed for nuclear dust structure in NGC 4450 and NGC 4736 (Elmegreen et al., 2002). Our simulations do not show it directly because the dissipation rate in the simulations is much higher than in the derivation above, as the model gas is forced to retain the same velocity dispersion (i.e., isothermal) with disk thickness $H=0$.

5. Discussion

We have presented the results of 2D hydrodynamic simulations for nuclear spirals driven by an oval distortion or a weak bar potential. Due to the geometric effect, nuclear spirals turn to shocks at small radii even if the potential perturbations are very weak. The background shear affects the morphologies of the nuclear spirals and shocks significantly, such that spirals are tightly (loosely) wound and the shocks are oblique (perpendicular) when shear is high (low). Our high-shear HS model forms a circumnuclear disk embedded with loosely-wound spiral shocks by accumulating inflowing gas. Circumnuclear disks with spirals seen in external galaxies, for example, in NGC 1097 (Hsieh et al., 2011; Onishi et al., 2015), could form by a similar process. We note however that NGC 1097 and Milky Way are strongly barred galaxies, so that nuclear rings produced by the strong bars might have affected the formation of circumnuclear disks.

Most previous theoretical work that studied gas accretion onto SMBHs considered giant elliptical galaxies embedded in the hot gaseous halos of their galaxy clusters, finding that gas accretion is chaotic and intermittent, occurring as the cold clouds formed by cooling and thermal instability in the hot gas rain down (Pizzolato & Soker, 2005; Sharma et al., 2012; Gaspari et al., 2013; Voit et al., 2015). Indeed, a cold, clumpy accretion flow has recently been observed in the nucleus of the brightest cluster galaxy in Abell 2597 (Tremblay et al., 2016). Our results suggest that nuclear spiral shocks at the centers of gas-rich disk galaxies, rather than thermal processes, can feed their SMBHs. These shocks are limited to the very central regions where driven pressure waves converge and steepen to become nonlinear. The radial inflow velocity of $\sim 50{\rm\,km\,s^{-1}}$ observed in the center of NGC 1097 (Fathi et al., 2006) is consistent with our model.

We are grateful to D. Elmegreen for help in making Figure 1 and to J. Stone for permission to use the Athea++ code at the developing stage of the code. This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 3348-20160021). The computation of this work was performed on the Linux cluster at KASI (Korea Astronomy and Space Science Institute).