CENTRAL REGIONS OF BARRED GALAXIES: TWO-DIMENSIONAL NON-SELF-GRAVITATING HYDRODYNAMIC SIMULATIONS

The real interstellar gas is multiphase and turbulent, with temperatures differing by a few orders of magnitude (e.g., Field et al. 1969; McKee & Ostriker 1977, 2007). For simplicity, we model this highly inhomogeneous gas using an isothermal equation of state with an effective sound speed c_s that includes a contribution due to turbulent motions. We have calculated 15 different models in which we vary both c_s and the initial mass M_BH(0) of the central BH as parameters. Table 1 lists the properties of each calculation. The sound speed is chosen to vary between 5 and 20 km s⁻¹. Models with a postfix bh0 do not initially possess a central BH, and are similar to those in Piner et al. (1995). Models with the postfix bh7 or bh8 have a BH with mass M_BH(0) = 4 × 10⁷ M_☉ and 4 × 10⁸ M_☉, respectively; they are analogous to models in Maciejewski (2004b) and Ann & Thakur (2005). For most models, we fix the BH mass to its initial value, but we also consider three additional models (cs05bh0t, cs20bh0t, and cs20bh7t) in which the BH mass is varied with time according to $M_{\rm BH}(t) = M_{\rm BH}(0) + \int _0^t \dot{M}(t^\prime)dt^\prime$, where $\dot{M}$ is the mass inflow rate across the inner boundary (see below). These time-varying M_BH models allow us to study the effect of BH growth due to the gas accretion on bar substructures.

Figure 1 plots the net circular rotation curves together with a contribution from each component when M_BH = 4 × 10⁸ M_☉. The solid and dashed lines are along the bar major and minor axes, respectively. The circular velocity is almost flat at ∼200 km s⁻¹ in the outer parts. Without the BH, the rotation curve v_c would rise linearly with R close to the center, but the presence of the BH results in v_c∝R^−1/2 for R ≲ 0.1 kpc. This rapid increase of v_c will render the gaseous orbits in the very central regions highly resistant to pressure perturbations, resulting in smaller mass inflow rates than the cases without it, as we show below.

Zoom In Zoom Out Reset image size

Figure 2 shows the characteristic angular frequencies, Ω − κ/2, Ω, and Ω + κ/2 along the bar major and minor axes as solid and dashed lines, respectively.⁸ Here, Ω² ≡ R⁻¹dΦ_ext/dR and κ² ≡ R⁻³d(R⁴Ω²)/dR denote the angular and epicyclic frequencies, respectively. The horizontal dotted line in each panel represents the bar pattern speed of Ω_b = 33 km s⁻¹ kpc⁻¹, with the corotation resonance (CR) located at R_CR = 6 kpc for all models. For bh0 models with no BH, the Ω − κ/2 curve peaks at R_max = 0.53 kpc and is equal to Ω_b at the two ILRs with radii of R_IILR = 0.19 kpc and R_OILR ≈ 2 kpc. Because the rotation curve rises steeply toward the center, bh7 and bh8 models have only a single ILR at R_ILR ≈ 2 kpc. The Ω − κ/2 curve in bh7 models attains its local maximum and minimum at R_max = 0.53 kpc and R_min = 0.19 kpc, respectively. In bh8 models, on the other hand, it increases monotonically with decreasing R since the BH dominates the gravitational potential. We will show in Section 4 that the shape of nuclear spirals depends critically on the sign of d(Ω − κ/2)/dR (e.g., Buta & Combes 1996).

Zoom In Zoom Out Reset image size

To solve Equations (1) and (2), we use the two-dimensional grid-based code CMHOG in cylindrical geometry (Piner et al. 1995). CMHOG implements the piecewise parabolic method in its Lagrangian remap formulation (Colella & Woodward 1984), which is third-order accurate in space and has very little numerical diffusion (viscosity). All the runs are carried out in a frame corotating with a bar whose major axis is aligned along the y-axis (i.e., ϕ = ±π/2), so that the bar potential remains stationary in the simulation domain. By assuming a reflection symmetry with respect to the galaxy center, the simulations were performed on a half-plane with −π/2 ⩽ ϕ ⩽ π/2 constructed by making a cut along the bar major axis.

As mentioned in Section 1, the original version of CMHOG used by Piner et al. (1995) contained a serious bug in the way the gravitational forces were added to the hydrodynamic equations. The CMHOG code places a bar potential Φ_bar with the major axis aligned along the x-axis, calculates the bar forces f_x = −∂Φ_bar/∂x and f_y = −∂Φ_bar/∂y, and then transforms them into cylindrical coordinates. In the original version of the code, the incorrect transform relations f_R = f_xcos ϕ + f_ysin ϕ and f_ϕ = f_xsin ϕ − f_ycos ϕ were used. In fact, the correct transformation rule for the azimuthal force should be f_ϕ = −f_xsin ϕ + f_ycos ϕ. With the sign of the azimuthal force reversed (but the radial force correct), the flows in models computed using the original CMHOG code behave as if the bar potential were aligned parallel to the y-axis, but with forces quite different from the intended ones. Other than these force transformations, the complex hydrodynamic algorithms in CMHOG are all implemented correctly, and were well tested in the original paper by Piner et al. (1995). Therefore, previous numerical studies based on CMHOG should remain valid as long as they do not adopt the incorrect transformations of the bar forces inherited from Piner et al. (1995). Unfortunately, the results of Piner et al. (1995) were compromised by a trivial sign error in the coordinate transform relations for the bar forces.

Figure 3 compares the results for a typical simulation run with the original and corrected version of CMHOG. For this test, the grid resolution is taken identical to that in Piner et al. (1995) with 251 and 154 zones in the radial and azimuthal directions, respectively. The figure plots the logarithm of the surface density at t = 300 Myr. Several differences are apparent. For example, the gas around the CR at R = 6 kpc in the left panel is largely evacuated and has corrugated streams linked to the ends of the bar major axis, whereas the CR region is relatively featureless in the right panel. The nuclear ring in the left panel is fairly smooth, while it is quite clumpy in the right panel. Most importantly, the very central region inside the ring is almost unperturbed in the left panel, while the right panel shows spiral structures in the central region. These differences suggest that the original CMHOG code is unable to properly model the flow in the central regions, especially weak nuclear spirals. We have also run the same model using other grid-based codes adopting Cartesian coordinates, such as TVD (Kim et al. 1999) and Antares (Yuan & Yen 2005) as well as the particle-based GADGET code (Springel et al. 2001), in all of which the bar forces are calculated by taking finite differences of Φ_bar rather than using f_x and f_y directly. The new version of the CMHOG code used in this work gives results which are much more similar to the results of these other codes, which gives us further confidence that the gravitational forces due to the bar are now being treated correctly.

Zoom In Zoom Out Reset image size

To resolve the central regions accurately, we set up a non-uniform, logarithmically spaced cylindrical grid with 1024 radial zones extending from 0.02 kpc to 16 kpc and 480 azimuthal zones covering the half-plane. This makes the zones approximately square-shaped throughout the grid (i.e., ΔR = RΔϕ). The resulting grid spacing is ΔR = 0.13, 6, and 100 pc at the inner radial boundary, R = 1 kpc where nuclear rings typically form, and the outer radial boundary, respectively. This increases the resolution in the inner regions by over an order of magnitude, in comparison to the models presented in Piner et al. (1995). This level of grid resolution is crucial to resolve nuclear spiral structures within R = 1 kpc. We use outflow and continuous boundary conditions at the inner and outer radial boundaries, respectively, while the azimuthal boundaries are periodic. The gas crossing the inner boundary is considered lost from the simulation domain. We keep track of the total mass crossing the inner boundary in order to study the mass inflow rates into the galactic nucleus.

Each model starts from a uniform disk with surface density Σ₀ = 10 M_☉ pc⁻² that is rotating in force balance with an axisymmetric gravitational potential without a bar. In order to avoid strong transients in the fluid flow caused by a sudden introduction of the bar, we slowly introduce the bar potential over one bar revolution time of 2π/Ω_b = 186 Myr. This is accomplished by increasing the bar central density ρ_bar linearly with time and decreasing the bulge central density ρ_bul, while keeping the net central density ρ_bar + ρ_bul fixed. This ensures that the shape of the total gravitational potential Φ_ext, when averaged along the azimuthal direction, is unchanged with time. All the models are run until 500 Myr. This corresponds to 1.2 × 10⁴ and 10 orbits at the inner and outer radial boundaries, respectively, for bh8 models with M_BH = 4 × 10⁸ M_☉.

Figure 4 plots snapshots of the logarithm of the gas density at a few selected epochs in the inner regions of Model cs05bh0. The bar is oriented vertically along the y-axis, and the gas is rotating in the counterclockwise direction relative to the bar. The images extend to 6 kpc on either side of the center, corresponding to the CR radius, outside of which the gas remains almost unperturbed.⁹ Piner et al. (1995) found that the CR regions exhibit time-dependent flow structures, as reproduced in Figure 3(a). On the other hand, Figure 4 shows that the CR regions in our simulations are quite stable and exhibit only at late time low-amplitude wavelike features entrained by the dense gas located at the bar ends, similar to the results of SPH simulations (e.g., Englmaier & Gerhard 1997; Patsis & Athanassoula 2000). This indicates that the complicated structures seen in Piner et al. (1995) were likely an artifact of the errors in their force transformation.

Zoom In Zoom Out Reset image size

A striking feature of each of the snapshots shown in Figure 4 at times greater than 120 Myr is large amplitude oscillations in the density in the ring and dust lanes. These features are also seen in the simulations of Wada & Koda (2004) for spiral shocks, and are attributed to the "wiggle instability." As we will discuss below, this shock instability appears to be caused by vorticity generation in curved shocks.

An introduction of the non-axisymmetric bar potential induces perturbations on the gas orbits, causing them to deviate from circular trajectories. The gas density increases (decreases) in regions where neighboring orbits come close together (diverge). When t = 60 Myr, the overdense regions are preferentially located downstream from the bar major axis (Figure 4(a)). At this time, the overdensity produced by orbit crowding is largest at R ∼ 3 kpc. Since the perturbing force by the bar is weak inside R ∼ 1 kpc, the overdensity there is correspondingly small and not readily discernible. Over time, the overdense regions become narrower and sharper as the bar amplitude grows and eventually develop into shock fronts at around t = 100 Myr. In what follows, we term these narrow shocks off-axis shocks.

A nuclear ring is beginning to shape at this time, as well. To illustrate the formation of nuclear rings in our models, we plot as solid lines in Figure 5 instantaneous streamlines of the gas that starts from Point A marked at (x, y) = (1.5, −2.7) kpc in Model cs05bh0 and from (x, y) = (1.4, −2.5) kpc in Model cs20bh0 with c_s = 20 km s⁻¹ and no BH at t = 100 Myr. The two dotted circles in Figure 5(a) indicate the inner and outer ILRs at R_IILR = 0.19 kpc and R_OILR = 2 kpc. Note that the thick lines representing the overdense ridges in both models directly cross the outer ILR. The changes of the azimuthal and radial velocities along the streamlines are shown in Figures 5(b) and (c), where the dotted line indicates the equilibrium rotation curve of the model galaxy with no BH. On emerging from the overdense region (Point A), the gas moves radially inward on its epicycle orbit and increases (decreases) its azimuthal (radial) velocity due to the Coriolis force. It reaches Point B closest to the center when it attains v_R = 0 and largest v_ϕ. After this point, it moves radially outward, decreasing v_ϕ until it hits the off-axis shock at Point C. The gas loses a significant amount of angular momentum at the shock and begins to fall in. In addition, the shocked gas is swept by other shocked gas flowing from the bar-end regions along the shocks. Note that the shape of the off-axis shocks shown in Figure 5(a) coincides with the gas streamline from Points C to D, indicating that all the gas after crossing the shocks moves radially in along the shock fronts in the developing stage of the nuclear rings.

Zoom In Zoom Out Reset image size

As the shocked gas moves along the shock fronts from Point C, it gradually rotates faster again. When the azimuthal velocity of the gas is increased to the level comparable to the equilibrium circular velocity at some radius R, it begins to follow a closed orbit (Point D), forming a nuclear ring at that radius. In other words, the centrifugal barrier inhibits the inflowing gas from moving further in. Regardless of the BH mass, this happens at R ∼ 0.8–1.2  kpc in models with c_s = 5 km s⁻¹ and R ∼ 0.4–0.6 kpc in models with c_s = 20 km s⁻¹ (Figures 4(b) and (c)). The facts that the off-axis shocks penetrate the ILR in bh7/bh8 models, the outer ILR in bh0 models, and that the ring positions are almost independent of M_BH when the rings are beginning to form suggest that the ring formation is unlikely to be governed by resonances. For models with low sound speed, the shape of nuclear rings is similar to an x₂-orbit. Clearly, the presence of the nuclear ring prevents the shocked gas from infalling directly to the nucleus.

Since the off-axis shocks are curved, Crocco's theorem ensures that vorticity can be generated at the shock fronts. Figure 6 plots snapshots of the potential vorticity $\xi \equiv |\nabla \times \mathbf {u} + 2\mathbf {\Omega }|/\Sigma$ relative to the initial value ξ₀ near the shocks at t = 90, 100, 110, and 120 Myr of the standard model. At t = 90 Myr, ξ/ξ₀ is largest along the shocks. Vorticity produced at the shocks is advected with the background flows and enters the shock fronts at the opposite side after a half revolution. Vorticity grows secularly with time by successive passages across the shocks. When vorticity achieves substantial amplitudes, it causes the shock fronts to wiggle and fragment into small clumps with high vorticity (see Figure 4(c)). The process of clump formation along the shocks in our models bears remarkable resemblance to the wiggle instability of spiral shocks found by Wada & Koda (2004) (see also Kim & Ostriker 2006). These clumps are carried radially inward and add to the nuclear ring, making the latter fairly inhomogeneous (Figures 4(d) and (e)).

Zoom In Zoom Out Reset image size

The off-axis shocks shown in Figure 5 are not stationary largely because the bar potential is not yet fully turned on. As the strength of the bar potential keeps increasing, they become stronger and move slowly toward the bar major axis. Gas that is added to the ring from the off-axis shocks has increasingly lower angular momentum with time, causing the ring to shrink in radius with time. After the bar potential is fully turned on, the off-axis shocks become gradually weaker as the amount of gas in the mid-bar regions lost to the ring increases with time. At the same time, orbital phase mixing and frequent clump collisions in the ring make the latter rounder and align its semimajor axis in the direction perpendicular to the bar major axis. Note that the rings are always attached to the inner end of the off-axis shocks.

At t = 300 Myr, Model cs05bh0 reaches a quasi-steady state in the sense that temporal changes in the overall flow pattern are very slow, and the locations of the shocks and rings do not vary much with time. Figure 4(f) overplots some of the x₁-orbits (dotted curves aligned vertically) and x₂-orbits (solid curves aligned horizontally), showing that the shape of the nuclear ring matches well with an x₂-orbit, while the off-axis shocks closely follow one of the x₁-orbits over the whole length of the shocks. This is because when c_s = 5 km s⁻¹ the impact of thermal pressure on the gas orbits is much smaller than that of the gravitational and centrifugal forces, so that pure orbit theory (neglecting pressure forces) is a good description. When c_s ≳ 15 km s⁻¹, however, thermal pressure gradients strongly affect gas orbits, and thus the morphology of substructures in the central regions is modified, as we will discuss below.

Even if the gravitational potential is the same, the flow morphology and velocity fields differ considerably depending on the sound speed. Figure 7 shows instantaneous streamlines plotted over the logarithm of the density distribution in Models cs05bh0 and cs20bh0 at t = 300 Myr. Red curves denote the streamlines that go through the off-axis shocks, while those enveloping the off-axis shocks are represented by green curves. In all models, the off-axis shocks are almost parallel to x₁-orbits. They start from the bar major axis at the outer ends, offset toward downstream in the mid-bar regions, and connect to the nuclear rings at the inner ends. The mean offset of off-axis shocks from the bar major axis is larger for models with smaller c_s.

Zoom In Zoom Out Reset image size

The outer end regions of the off-axis shocks have complicated density structures including the "4/1-spiral shocks" marked with a red arrow in Figure 7(a) (Englmaier & Gerhard 1997) and the enhanced density ridges (a white arrow) termed "smudges" by Patsis & Athanassoula (2000). As discussed by Englmaier & Gerhard (1997), the 4/1-spiral shocks are produced by collisions of gas moving on x₁-orbits with that on the 4/1-resonant family (e.g., Contopoulos & Grosbøl 1989). The gas loses angular momentum at the shocks and subsequently switches to lower orbits. As the streamlines in green display, in models with c_s = 5 km s⁻¹, the 4/1-spiral shocks are quite strong and spatially extended, so that the orbits after the shocks become relatively radial and converge at the opposite side of the bar, building a smudge after about a half revolution. Collisions of streams off the 4/1-spiral shock and the smudge on the same side of the bar funnel the gas at the intersections to an x₁-orbit, which are the starting points of the off-axis shocks. When c_s = 20 km s⁻¹, on the other hand, the 4/1-spiral shocks are short and weak, and the streamlines off the shocks diverge, so that a dense ridge does not form. Since the gas becomes less compressible with increasing c_s, steady off-axis shocks in models with large c_s can be supported only in inner regions where the bar perturbations are sufficiently strong. This explains why the mean offset of the off-axis shocks from the bar major axis becomes smaller as c_s increases (e.g., Englmaier & Gerhard 1997). With weak 4/1-spiral shocks and no smudge, the gas in the bar-end regions in a model with c_s = 20 km s⁻¹ is comparatively unsteady, sometimes generating small dense blobs that move inward along the off-axis shocks.

To quantify the shock properties, we place a slit in each model, indicated by the short white lines in Figure 7. The slit starts from (x, y) = (0, −1.5 kpc) and runs perpendicular to the local segment of the off-axis shocks, roughly at R ∼ 1.5–1.8  kpc. Figure 8 plots the profiles along the slit of surface density, velocities, and the compression factor,

$$\begin{equation} \alpha \equiv -(\nabla \cdot \mathbf {v})\Delta R/c_s, \end{equation} \tag{ 3 }$$

the last of which can be used as an effective measure of the shock strength (e.g., Maciejewski 2004b; Thakur et al. 2009).¹⁰ The gas is flowing from left to right in the increasing direction of s, where s measures the distance along the slit from the starting point. Table 2 gives the shock properties for all models at t = 300 Myr: s_sh is the position of the off-axis shocks along the slit, Σ_max is the peak density after s_sh, ${\mathcal M}_\bot \equiv v_\bot /c_s$ and ${\mathcal M}_\Vert \equiv v_\Vert /c_s$ are the Mach numbers of the incident flows in the directions perpendicular and parallel to the shocks, respectively, $i_{\rm sh}\equiv \tan ^{-1}({\mathcal M}_\bot /{\mathcal M}_\Vert)$ is the inclination angle of the pre-shock velocity relative to the shock front, $dv_\Vert /ds$ quantifies the velocity shear in the post-shock region, and α_max is the maximum value of the compression factor occurring at the shock front. It is apparent that the off-axis shocks tend to move toward the bar major axis with increasing c_s, while there is no clear dependence of s_sh on the BH mass. The compression factor at the shock is insensitive to M_BH and scales roughly with c_s as α_max ∼ 7.7(c_s/5 km s⁻¹)^0.92.

Zoom In Zoom Out Reset image size

Table 2. Properties of Off-axis Shocks

Model	ssh	Σmax/Σ0	${\mathcal M}_\bot$	${\mathcal M}_\Vert$	ish	$dv_\Vert /ds$	αmax
	(kpc)				(deg)	(103 km s−1 kpc−1)
cs05bh0	0.98	5.7	12.1	10.9	47.7	2.7	7.2
cs10bh0	0.72	3.6	6.0	4.9	50.8	1.8	3.5
cs15bh0	0.55	5.7	5.7	1.0	80.4	2.3	2.6
cs20bh0	0.45	5.0	3.5	1.2	71.3	1.8	1.6
cs05bh7	0.93	5.8	15.4	6.3	67.9	1.8	7.4
cs10bh7	0.63	9.8	8.1	−1.3	−81.4	1.6	4.7
cs15bh7	0.47	9.0	6.6	−1.2	−79.8	1.7	3.0
cs20bh7	0.31	5.7	2.3	1.3	59.4	1.0	1.9
cs05bh8	0.96	6.1	16.7	7.2	66.7	2.0	7.5
cs10bh8	0.74	7.4	9.6	−0.1	−89.5	2.3	6.3
cs15bh8	0.43	6.8	7.2	−0.8	−83.3	1.4	3.0
cs20bh8	0.32	21.6	5.9	−2.6	−66.6	1.4	2.9

Notes. s_sh is the shock position along the slit; Σ_max is the maximum density attained immediately after s_sh; ${\mathcal M}_\bot$ and ${\mathcal M}_\Vert$ are the Mach numbers of the incident flow perpendicular and parallel to the shock, respectively; i_sh is the inclination angle of the incident flow relative to the shock; $dv_\Vert /ds$ is the velocity shear in the post-shock region; α_max is the maximum value (occurring at the shock front) of the compression factor $\alpha \equiv -(\nabla \cdot \mathbf {v}) \Delta R/c_s$.

Download table as:  ASCIITypeset image

Naively, one would expect that the off-axis shocks become weaker as the sound speed increases, since the density jump in planer isothermal shocks is proportional to ${\mathcal M}_\bot ^2$. However, this is not the case, as Figure 8 and Table 2 demonstrate. There are several reasons for this. First, it is the velocity component normal to the shock front v_⊥ that determines the shock jump conditions, and because the inclination angle of the streamlines which enter the shock varies with location and with c_s, v_⊥ varies in a complicated fashion. For example, Figure 8 shows that for off-axis shocks formed at R ∼ 1.5–1.8  kpc, Model cs15bh0 with c_s = 15 km s⁻¹ has the largest peak density as well as the largest v_⊥ = 85 km s⁻¹ and i_sh = 80°. On the other hand, Model cs10bh0 has the smallest v_⊥ = 60 km s⁻¹ (with i_sh = 51°) and thus the lowest density enhancement. Since the sound speed is lower, Model cs05bh0 with v_⊥ = 60 km s⁻¹ produces Σ_max comparable to that in Model cs15bh0. Second, we note that the Rankine–Hugoniot jump conditions for stationary planar shocks are not applicable to the curved and two-dimensional off-axis shocks formed in our simulations. As Figure 7 displays, the flows are fully two dimensional in the sense that streamlines diverge before the shocks and converge after the shock with the radial inflow coming from the regions near the end of the bar. The fact that the compression factor α measured at the shock front is smaller than ${\mathcal M}_\bot -{\mathcal M}_\bot ^{-1}$ expected from planer isothermal shocks also indicates that the shocks are two dimensional. Finally, the density and velocity fluctuations generated by the vortex-generating instability are important around the off-axis shocks especially for models with small c_s, so that the flows are not strictly stationary.¹¹ Note that the shocked gas has strong velocity shear, amounting to $dv_\Vert /ds\sim (1\hbox{--}3)\times 10^3\; {\rm km}\;{\rm s}^{-1}\;{\rm kpc}^{-1}$, which is about 10 times larger than the velocity shear arising from galaxy rotation in the solar neighborhood. Such strong shear can stabilize the high-density, off-axis shocks against self-gravity.

Gas that loses angular momentum at the off-axis shocks flows radially inward and forms a nuclear ring in the central regions. Figure 9 shows diverse morphological features produced in the regions with |x|, |y| ⩽ 1 kpc for all models at t = 300 Myr. Figure 10 overplots instantaneous streamlines for a few selected models. Some models (with low c_s) have a nuclear ring together with inner spiral structures, some models (with high c_s and small M_BH) have a ring with no apparent spirals, while others (with high c_s and large M_BH) possess only nuclear spirals without an appreciable ring.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When c_s = 5 km s⁻¹, the nuclear rings are quite narrow and clearly decoupled from the nuclear spirals. Even though the ring has a large density, the low sound speed makes the effect of thermal pressure on the gas orbits insignificant. The gas around the ring in Model cs05bh0 thus follows x₂-orbits fairly well and the shape of the ring does not deviate considerably from x₂-orbits (Figure 10(a)). When c_s = 20 km s⁻¹, on the other hand, the pressure force in the central regions becomes important and affects the shape of the gas streamlines. Even the inflowing gas that arrives at the contact points between the off-axis shocks and the nuclear ring takes very different orbits depending on its location. Some gas at the outer parts of the contact points is pushed out by the pressure gradient and follows trajectories that are much rounder than x₂-orbits, while the gas in the inner parts is forced to take inner highly eccentric orbits (Figure 10(b)). Consequently, the gas in the central regions in Models cs20bh0 spreads out spatially and forms a ring that is more circular and broader than in Model cs05bh0. Since the presence of a central BH increases the initial angular momentum of the gas in the central regions, the pressure effect becomes less important as M_BH increases. Figure 10 shows that the pressure distortion of x₂-orbits is still significant for M_BH = 4 × 10⁷ M_☉, while the gas orbits in the very central parts at R ≲ 0.2 kpc remain almost intact when M_BH = 4 × 10⁸ M_☉. In Model cs20bh8, some gas at R ∼ 0.5–0.8 kpc temporarily moves radially outward due to the radial pressure gradient built up by the background gas and is subsequently swept inward by other gas flowing in along the off-axis shocks.

Figure 11 plots the radial distribution of the averaged gas surface density 〈Σ〉, averaged both azimuthally and temporally over t = 300–500 Myr for all models. The locations of the ILRs as well as Rmax and Rmin corresponding to the local maximum and minimum of the Ω − κ/2 curve are indicated by arrows on the abscissa. Table 3 gives the inner and outer radii, Rin and Rout, of the ring, the mass-weighted ring radius Rring = ∫RoutRin〈Σ〉RdR/∫RoutRin〈Σ〉dR, the peak density 〈Σ〉max, and the mean density Σring = ∫RoutRin〈Σ〉dR/(Rout − Rin) of the ring in each model. Here, Rin and Rout are defined as the radii where 〈Σ〉 = 〈Σ〉max/5. Note that Models cs15bh8 and cs20bh8 do not harbor a well-defined nuclear ring, as will be discussed in the next section.

Zoom In Zoom Out Reset image size

All rings that form are located within R_OILR if there are two ILRs or R_ILR if there is a single ILR, indicating that the formation of a nuclear ring does not require the presence of two ILRs. However, there is in general no direct connection between the ring positions and R_OILR or R_ILR. When c_s = 5 km s⁻¹, the rings are all located at R ∼ 0.6–1.2 kpc, independent of M_BH. The mass-weighted radius is R_ring ∼ 0.90–0.92  kpc, indicating that the ring position is not governed by the shape of the Ω − κ/2 curve. When c_s = 10 km s⁻¹, the ring radius decreases to R_ring ∼ 0.67–0.68  kpc, again insensitive to the BH mass, consistent with the tendency for the off-axis shocks to move closer to the bar major axis as c_s increases. Rings with c_s = 10 km s⁻¹ have a larger surface density than those with c_s = 5 km s⁻¹, corresponding approximately to a constant ring mass (i.e., Σ_ring∝R⁻¹_ring). As c_s increases further, high thermal pressure provides strong perturbations for x₂-orbits and tends to spread out the gas in the central parts, resulting in a broad distribution of 〈Σ〉 at R ≲ 0.5 kpc. When the BH is not massive enough, these perturbations wipe out coherent, weak spiral structures that formed earlier in the nuclear regions.

Because the presence of a central BH dominates the potential only in the central region, the allowance for the growth of the BH due to mass accretion does not make a significant difference in the regions outside the ring. Even inside the ring, Figure 11 and Table 3 show that the changes in the ring size R_ring and the ring density Σ_ring caused by the temporal change of M_BH are less than 4% and 10%, respectively. We will show below that BH growth in our simulations also does not significantly affect the properties of nuclear spirals and mass inflow rates.

To study the spiral features, it is convenient to show the logarithm of the gas surface density in logarithmic polar coordinates. Figure 13 plots snapshots of gas surface density of Models cs05bh0 and cs10bh0 on the ϕ − log R plane. Any coherent features with a positive (negative) slope on this plane are leading (trailing) waves. Only the regions with R ⩽ 1 kpc are shown. Two horizontal lines mark R_IILR and R_max where the Ω − κ/2 curve is a locally maximum. At early time, the non-axisymmetric bar potential induces weak m = 2 perturbations in the central regions. Perturbed gas elements in the galactic plane follow slightly elliptical orbits, which are closed in a frame rotating at Ω − κ/2 and thus precess at a rate Ω − κ/2 near the ILRs when seen in a stationary bar frame. As succinctly depicted in Buta & Combes (1996), due to collisional dissipation of gas kinetic energy occurring on converging orbits, the gas forms spiral structures whose shape depends critically on the sign of d(Ω − κ/2)/dR such that spirals are leading (trailing) where d(Ω − κ/2)/dR is positive (negative). Figure 13 indeed shows that when t = 50 Myr the perturbed density is leading at R < R_max and trailing at R > R_max, although the trailing features are soon overwhelmed by the nuclear ring. Located away from the nuclear ring, however, the inner leading waves are able to grow with time and eventually develop into shock waves at t ∼ 200 Myr. These nuclear spirals are short, extending over R ∼ 0.13–0.50 kpc, quite open with a pitch angle of i_p = −30°, and almost completely detached from the nuclear ring.

Zoom In Zoom Out Reset image size

Figure 14 plots the azimuthal distributions of surface density and velocities of the nuclear spirals at R = 0.25 kpc in Models cs05bh0 (solid lines) and cs05bh0t (dotted lines) when t = 500 Myr. Over the course of the orbits, the changes of the radial and azimuthal velocities associated with the spirals amount to ∼100 km s⁻¹, which is indeed large enough to induce shocks. The peak densities occurring at ϕ ∼ 100°, 280° for Model cs05bh0 correspond to shock fronts with a compression factor of α ∼ 3.4. The density bumps at ϕ ∼ 135°, 315° are produced by waves launched from the inner boundary. In Model cs05bh0t, the BH mass is increased to M_BH ∼ 10⁵ M_☉ due to mass inflow, which supports slightly stronger, more leading spiral shocks than in Model cs05bh0 with no BH. Despite continual perturbations by traveling trailing waves propagating from the inner boundary, the nuclear spirals in these models last until the end of the run. That leading nuclear spirals are persistent when the sound speed is small is consistent with the results of SPH simulations reported by Ann & Thakur (2005).

Zoom In Zoom Out Reset image size

The usual WKB dispersion relation for tightly wound, linear-amplitude waves in a non-self-gravitating medium reads

$$\begin{equation} (\omega -m\Omega)^2 = k^2c_s^2 + \kappa ^2, \end{equation} \tag{ 4 }$$

where ω is the wave frequency and k and m are the radial and azimuthal wavenumbers (e.g., Goldreich & Tremaine 1978). For m = 2 waves corotating with a bar (i.e., ω = 2Ω_b), Equation (4) becomes

$$\begin{equation} k = \pm \frac{2}{c_s} \sqrt{ (\Omega + \kappa /2 - \Omega _b)(\Omega -\kappa /2-\Omega _b)}, \end{equation} \tag{ 5 }$$

(e.g., Englmaier & Shlosman 2000; Maciejewski 2004a), indicating that nuclear spirals corotating with the bar can exist only in the regions between R_IILR and R_OILR when there are two ILRs.¹² However, the top row in Figure 13 reveals that the nuclear spirals in Model cs05bh0 extend slightly inward of R_IILR. This seemingly contradicting result is due to nonlinear effects which are not captured by the WKB theory. The velocity perturbations associated with these spirals are so large that fluid elements just outside R_IILR can move across the inner ILR over the course of their epicycle orbits, providing perturbations for the gas at R < R_IILR that responds passively. In Model cs05bh0, the radial velocity perturbation is Δv_R = 52 km s⁻¹. Since the epicycle frequency is κ = 1100 km s⁻¹ kpc⁻¹ at R = R_IILR, the corresponding radial amplitude of the epicycle orbits is estimated to be ΔR = Δv_R/κ = 0.05 kpc, which is in good agreement with the radial extent of the nuclear spirals inward of R_IILR.

Models without a BH and with c_s ⩾ 10 km s⁻¹ do form nuclear spirals at early time, but they are all transient, lasting less than 200 Myr. The bottom row of Figure 9 shows how nuclear spirals are destroyed in Model cs10bh0. The nuclear ring in this model is not only located inside R_max but also has large thermal pressure, continuously generating sonic perturbations that propagate radially inward. Because of the background shear, the perturbations are preferentially in the form of trailing waves which interact destructively with the leading spirals that formed at R < R_max, destroying the latter. The destruction of nuclear spirals happens at t = 185 Myr for Model cs10bh0. This occurs earlier as c_s increases, since disturbing pressure perturbations are correspondingly stronger.

Since the presence of a central BH greatly changes the Ω − κ/2 curve in the central regions, it is interesting to explore how the morphologies of nuclear spirals depend on the BH mass. In bh7 and bh8 models with a single ILR, stationary waves in the bar frame can exist inside R_ILR all the way down to the center. Figure 15 plots snapshots of gas surface density for Models cs05bh7 and cs10bh7 on the logarithmic polar plane. The positions of R_max and R_min are indicated by the horizontal lines. As expected, the overdense perturbations produced by orbit crowding at early time (t = 50 Myr) have leading configurations at R_min < R < R_max and trailing configurations at R < R_min or R > R_max. The overdense regions at R > R_max are subsequently wiped out as the nuclear ring forms, while those at R < R_max grow into trailing nuclear spirals. In Model cs05bh7, the leading spirals at R_min < R < R_max also grow slightly until t ∼ 150 Myr to temporarily form "inner-trailing and outer-leading" structures represented by the double cross-hatching in Figure 12. As trailing perturbations from both the inner trailing spirals and the outer nuclear ring propagate and interfere with the leading spirals, it becomes increasingly difficult to identify coherent spiral structures at R_min < R < R_max. On the other hand, the trailing spirals at R < R_min keep growing until t ∼ 230 Myr after which their amplitude of Σ_peak/Σ_avg ≈ 2.0 remains more or less constant. They are approximately logarithmic in shape with a pitch angle of i_p = 85. Unlike in Model cs05bh0 where leading spirals are actually shocks, the trailing nuclear spirals in Models cs05bh7 are relatively weak (with the maximum compression factor of α_max ∼ 0.28 at t = 500 Myr) and never develop into shocks.

Zoom In Zoom Out Reset image size

In Model cs10bh7 (bottom row of Figure 15), the nuclear spirals have larger k and thus are more open than those in the c_s = 5 km s⁻¹ counterparts (see, e.g., Equation (5)). Since the nuclear ring in this model forms at R_ring ≈ 0.5 kpc, it can directly destroy the leading spiral at R_min < R < R_max, and feeds the inner trailing spirals by supplying trailing perturbations. Thus, the inner spirals grow both in strength and spatial extent to make contact with the nuclear ring at t = 210 Myr. At this time, all parts of the nuclear spirals turn to shocks with the maximum density of Σ_max/Σ₀ = 3.8 and the corresponding compression factor of α = 0.8 at R = 0.1 kpc. Gas passing through the spiral shocks loses angular momentum, increasing the mass inflow rate at the inner boundary. As the amount of gas lost in the nuclear regions increases, the density of the nuclear spirals decreases with time, but they remain as shocks with the compression factor of α ∼ 1.3 until the end of the run.

In Models cs15bh7 and cs20bh7, nuclear spirals start out with inner-trailing and outer-leading shapes, and evolve into trailing-only configurations, as in the other bh7 models. However, the highly eccentric orbits of the gas affected by thermal pressure dismantle the inner spiral structures almost completely in these models: thermal perturbations are so strong that a central BH with M_BH = 4 × 10⁷ M_☉ cannot enforce circular orbits in the very nuclear regions (see Figure 10(c)).

Since the Ω − κ/2 curve decreases monotonically with R < 1 kpc in bh8 models, nuclear spirals, if they exist, are all trailing as evident in Figure 16. For Model cs05bh8 (top row of Figure 16), the spirals evolve almost independently of, and are well separated from, the nuclear ring. At t = 500 Myr, they have a very small pitch angle i_p = 35 corresponding to $kR=2\cot i_p=33$ owing to the strong background shear. They are also very weak, with an amplitude of Σ_peak/Σ_avg = 2.1 at R ≈ 0.1 kpc. When c_s = 10 km s⁻¹, the pitch angle is increased to i_p = 58 and the nuclear spirals extend outward all the way to the nuclear ring. Except near the contact points, the spirals are still decoupled from the ring. With quite a large value of kR = 20, the component of the rotational velocity perpendicular to the spirals is not large enough to induce shocks.

Zoom In Zoom Out Reset image size

We have seen earlier that the large thermal pressure in the rings of models with c_s ⩾ 15 km s⁻¹ provides strong perturbations that make the gas orbits near the galaxy center highly eccentric, destroying nuclear spirals in bh0 and bh7 models. In bh8 models, however, the situation is quite different since a central BH with M_BH = 4 × 10⁸ M_☉ dominates the potential, keeping the orbits almost circular in the very central regions. Even with a large thermal pressure, the gas orbits there cannot be very eccentric, so that the nuclear spirals are protected from disruptive pressure perturbations. On the other hand, the ring material is quite distributed because of the large pressure gradients, feeding a trailing nuclear spiral that grows strongly in Models cs15bh8 and cs20bh8.

At t = 120 Myr, the spirals in Models cs15bh8 and cs20bh8 turn into shocks and touch the densest parts of the ring located at R ∼ 0.4 kpc. Because of the larger pitch angles, the shocks are stronger in the outer parts; the portions at R ≲ 0.1 kpc are weak shocks with density jumps of only ∼2. Similarly, the shocks in Model cs20bh8 are stronger than those in Model cs15bh8 since the former has more open spirals and a denser ring. In both models, the shocks near the ring are so strong that even the gas constituting the ring suffers from a significant loss of angular momentum at the intersecting points. At t = 150 Myr, the ring material is essentially dissected by the trailing spiral shocks and gradually moves toward the center. As the ring material continues to flow in, the spirals appear as a direct continuation of the off-axis shocks (t = 350 Myr), consistent with the results of Maciejewski (2004b).

Figure 17 plots the radial distributions of the maximum and mean densities as well as the maximum compression factor in the inner 1 kpc regions of Models cs15bh8 and cs20bh8 at t = 500 Myr. The inflow of the ring material in Model cs20bh8 is sufficiently strong that the gas is collected in the nuclear regions with R ∼ 0.02–0.1 kpc; the amount of the gas that goes out through the inner boundary is much smaller than that which comes in. With weaker shocks and thus less angular momentum loss, on the other hand, the destroyed ring gas in Model cs15bh8 is still mostly at R > 0.1 kpc at this time. Note that the density enhancement at R = 0.05 kpc resulting from the dissolution of the ring is ∼8Σ₀ and ∼170Σ₀ in Models cs15bh8 and cs20bh8, respectively: the mean contrast of the spirals that are weak shocks with the compression factor α ∼ 0.3 is Σ_peak/Σ_avg ∼ 1.6–1.8 at R = 0.05 kpc in both models. This increase of the gas surface density in the central regions is the primary reason for enhanced mass inflow rates at late time in bh8 models with large c_s.

Zoom In Zoom Out Reset image size

We have presented detailed numerical models that explore the formation of substructures produced by the gas flow in barred galaxies. Previous models based on particle simulations (e.g., Englmaier & Gerhard 1997; Ann & Thakur 2005; Thakur et al. 2009) did not have sufficient resolution to resolve nuclear spirals. On the other hand, studies that used the grid-based code CMHOG (e.g., Piner et al. 1995; Maciejewski 2004b) unknowingly made mistakes in the force evaluation for the bar potential, so that the results needed to be recomputed.

In this paper, we have corrected the errors in the original CMHOG code and run high-resolution hydrodynamical simulations. To resolve the nuclear regions, we employed a logarithmically spaced cylindrical grid, with a zone size of ΔR ⩽ 6 pc at R ⩽ 1 kpc where nuclear rings and spirals form. We have included the potential from a central BH and studied the flow properties as the mass of the BH M_BH and the sound speed c_s in the gas are varied. For simplicity, the effects of gaseous self-gravity and magnetic fields are not included. The main results of the present work are summarized as follows.

1. Off-axis shocks. The imposed non-axisymmetric bar potential provides gravitational torques to the otherwise circular-rotating gas, perturbing its orbit. The perturbed orbits crowd at the downstream sides of the bar major axis and produce overdense ridges that eventually develop into off-axis shocks. At a quasi-steady state, the off-axis shocks are overall almost parallel to x₁-orbits: they start from the bar major axis at the outer ends, are gradually displaced downstream as they move inward, and connect to the nuclear ring at the inner ends. While the positions of the off-axis shocks are almost independent of M_BH, since the effect of a BH is negligible at large radii, they depend on c_s in such a way that the shocks are, on average, located closer to the bar major axis as c_s increases. This is primarily because gas with larger c_s should be more strongly perturbed to induce shocks, which occurs deeper in the potential well and thus results in shocks on lower x₁-orbits (Englmaier & Gerhard 1997).

The off-axis shocks are in general curved. Flow streamlines are complicated near the shocks in that they diverge before the shocks and are promptly swept inward by inflowing gas right after the shocks. Therefore, the usual Rankine–Hugoniot jump conditions for planar, one-dimensional shocks are not applicable to the off-axis shocks. In fact, the shock strength, as measured by the peak density Σ_max at R ∼ 1.5–1.8  kpc from the center, is Σ_max/Σ₀ ∼ 3–6 for models with no BH and does not sensitively depend on c_s. The compression factor of the off-axis shocks is insensitive to the BH mass and depends on c_s roughly as α_max ∼ 7.7(c_s/5 km s⁻¹)^0.92. The off-axis shocks have very strong velocity shear amounting to ∼(1–3) × 10³ km s⁻¹ kpc⁻¹. This strong shear may suppress the growth of gravitational instability in the high-density off-axis shocks when self-gravity is included.

2. Nuclear rings. When gas passes through the off-axis shocks, it loses angular momentum, flows inward, and forms a nuclear ring where the centrifugal force balances the external gravity. The ring is attached to the inner ends of the off-axis shocks and thus becomes smaller in size as c_s increases. The rings that form in our models are all located inside the (outer) ILR, but this does not imply that the ring formation is related to the ILRs. When c_s is small, the pressure perturbations on the gas orbits are so weak that the rings are quite narrow and their shape is well described by x₂-orbits. The mean radius is R_ring ∼ 0.8–0.9  kpc when c_s = 5 km s⁻¹ and R_ring ∼ 0.5–0.6  kpc when c_s = 10 km s⁻¹, independent of the BH mass. This suggests that the ring position is not determined by the Ω − κ/2 curve, hence by the ILRs. When c_s ⩾ 15 km s⁻¹, on the other hand, large thermal pressure strongly affects the gas orbits in the nuclear rings. For example, some gas near the contact points between the ring and the off-axis shocks is forced out to follow relatively round orbits, while other gas is pulled in radially to make very eccentric orbits. These diverse gas orbits near the contact points tend to spread out the ring material, making the rings much broader than in models with smaller c_s.

3. Nuclear spirals. Since even the non-axisymmetric bar potential is nearly axisymmetric in the central parts, the gaseous responses inside a nuclear ring are not as dramatic as in off-axis shocks. Nevertheless, non-axisymmetric m = 2 perturbations are able to grow inside a ring and develop into nuclear spirals that persist for a long period of time, provided that c_s is small or M_BH is large. Although all models have weak spiral structures at early time, in models with large c_s they are soon destroyed by eccentric gas orbits as well as perturbations induced by the large pressure in the rings unless the BH mass is very large. When M_BH = 4 × 10⁸ M_☉, the disruptive pressure perturbations from the rings cannot penetrate the very central parts where the gas has extremely large initial angular momentum. In this case, the nuclear spirals are protected from the surrounding and thus are long-lived.

The shape of nuclear spirals is determined by the sign of d(Ω − κ/2)/dR such that spirals that form in the regions where d(Ω − κ/2)/dR is positive (negative) are leading (trailing), confirming the theoretical expectations of Buta & Combes (1996). With no BH, only the model with c_s = 5 km s⁻¹ has persistent leading spirals in the regions where the Ω − κ/2 curve is an increasing function of R. The leading spirals in this model are quite strong and develop into shocks with the peak density Σ_peak/Σ_avg ∼ 7.9 and the compression factor α ∼ 3.4 at R = 0.25 kpc at the end of the run. Models with M_BH = 4 × 10⁷ M_☉ initially have hybrid features comprising trailing spirals at R < R_min and leading spirals at R_min < R < R_max, where R_min and R_max refer to the radii where the Ω − κ/2 curve attains a local minimum and maximum, respectively. When c_s ⩽ 10 km s⁻¹, however, the leading parts in the hybrid spirals are destroyed by the trailing pressure waves launched by the ring, leaving only the weak trailing spirals behind. When c_s ⩾ 15 km s⁻¹, both leading and trailing spirals are destructed completely by the pressure perturbations. In models with M_BH = 4 × 10⁸ M_☉, d(Ω − κ/2)/dR < 0 in the whole nuclear regions, so that nuclear rings are always trailing. When c_s ⩽ 10 km s⁻¹, the spirals well separated from the ring are weak and tightly wound. When c_s ⩾ 15 km s⁻¹, on the other hand, the trailing spirals are fed by the gas in the nuclear ring and grow into shocks, the outer ends of which join the inner ends of the off-axis shocks smoothly. Although the nuclear spirals in these models have only modest density contrasts of Σ_peak/Σ_avg ∼ 1.7 at R = 0.05 kpc, the background density resulting from the inflows of the ring material is greatly enhanced (by more than two orders of magnitude when c_s = 20 km s⁻¹).

4. Mass inflow rates. Although gas experiences a significant loss in angular momentum when it meets off-axis shocks during galaxy rotation, this does not necessarily translate into the mass inflow all the way to the galaxy center. In models with c_s ⩽ 10 km s⁻¹, a narrow nuclear ring formed by gas with x₂-orbits inhibits further inflows of the gas, resulting in the mass inflow rate $\dot{M}$ through the inner boundary less than 0.01 M_☉ yr⁻¹, regardless of the BH mass. The mass inflow rates are greatly enhanced to $\dot{M}> 0.01\; {M}_{\odot }\; {\rm yr}^{-1}$ in models with M_BH ⩽ 4 × 10⁷ M_☉ and c_s ⩾ 15 km s⁻¹ or M_BH = 4 × 10⁸ M_☉ and c_s = 20 km s⁻¹ for different reasons depending on the BH mass. When M_BH ⩽ 4 × 10⁷ M_☉, the gas orbits are affected by thermal pressure and some gas in the ring can take highly eccentric orbits, directly falling into the inner boundary. In models with M_BH = 4 × 10⁸ M_☉, on the other hand, the gas orbits near the center are more or less circular, but the density in the nuclear regions is greatly enhanced because of strong nuclear spirals, increasing $\dot{M}$.

Nuclear rings play an important role in the evolution of barred galaxies by providing sites of active star formation near the centers (e.g., Buta 1986; Garcia-Barreto et al. 1991; Barth et al. 1995; Maoz et al. 2001; Mazzuca et al. 2008). Rings certainly consist of gas that migrates from outer parts inward by losing angular momentum at off-axis shocks, but what stops further migration to form a ring remains a matter of debate. As a trapping mechanism of the ring material, Combes (1996) proposed the non-axisymmetric bar torque that forces gas to accumulate between two ILRs or at a single ILR depending on the shape of the gravitational potential (see also Buta & Combes 1996), while Regan & Teuben (2003, 2004) favored the gas transitions from x₁- to x₂-orbits rather than orbital resonances. However, our numerical results show that a ring is formed at early time because of the centrifugal barrier that the migrating material feels. Later on, nuclear rings slowly shrink in size as gas with lower angular momentum gas is continuously added. Although nuclear rings are located in between two ILRs in models with no BH, this is a coincidence. Models with a central BH show that the specific ring positions are insensitive to the shape of the Ω − κ/2 curves, suggesting that nuclear ring formation is not a consequence of the orbital resonances. In fact, the gas flows that produce the ring are not in force balance and have a large radial velocity, so that the concept of resonances and ILRs is not applicable to nuclear rings (e.g., Regan & Teuben 2003). In addition, the notion of x₂-orbits as the gas trapping locations is meaningful only for small c_s.¹³ When the sound speed is large, thermal pressure at the contact points between the off-axis shocks and nuclear ring causes the gas orbits to deviate from x₂-orbits considerably.

The results of our simulations suggest that not all barred galaxies possess nuclear spirals at their centers. Long-lasting nuclear spirals exist only when either the sound speed is small or the BH mass is large: they do not survive in models with both large c_s and small M_BH. Two common views regarding the nature of nuclear spirals are low-amplitude density waves and strong gaseous shocks (see, e.g., Englmaier & Shlosman 2000; Maciejewski et al. 2002; Maciejewski 2004a, 2004b; Ann & Thakur 2005; Thakur et al. 2009). And our simulations indeed show that they are either tightly wound density waves or shocks when the pitch angle is relatively large (i_p > 6°). One may speculate that nuclear spirals are more likely to be shocks rather than density waves when c_s is small. In contrast to this prediction, however, nuclear spirals in models with M_BH = 4 × 10⁸ M_☉ are density waves when c_s is small and shocks when c_s is large. This is of course because as c_s increases, waves in nuclear regions tend to be more open (with smaller |k|) in the beginning, and they are subsequently supplied with more gas from the rings as they grow. This is entirely consistent with the results of Ann & Thakur (2005) who used SPH simulations to show that nuclear spirals are supported by shocks when c_s ≳ 15 km s⁻¹ in models with a massive BH. We note however that weak trailing spirals seen in our Model cs10bh8 are absent in Model M2 (c_s = 10 km s⁻¹ and M_BH = 4 × 10⁸ M_☉) of Ann & Thakur (2005), which is presumably due to an insufficient number of particles to resolve nuclear spirals in their SPH simulations.

Of 12 models with differing c_s and M_BH(0) that we have considered, only 1 model possesses leading spirals, suggesting that galaxies with leading nuclear spirals would be very uncommon in nature. To our knowledge, only two galaxies, NGC 1241 and NGC 6902, are known to possess leading features in the nuclear regions (Díaz et al. 2003; Grosbøl 2003).¹⁴ Based on our simulations, the existence of leading spirals at centers requires two stringent conditions: (1) the gas should be dynamically cold enough to protect nuclear spirals from nuclear rings and (2) there should be a wide range of radii with d(Ω − κ/2)/dR < 0 in the central parts, which can easily be accomplished when there are two ILRs (or without a strong central mass concentration). The second condition is consistent with the linear theory that predicts short leading waves propagating outward from the inner ILR (Maciejewski 2004a). The facts that NGC 6902 is a barred-spiral galaxy (Laurikainen et al. 2004) and does not show significant X-ray emissions indicative of AGN activities (Desroches & Ho 2009) are not inconsistent with the second requirement for the existence of leading nuclear spirals. Since NGC 1241 is a Seyfert 2 galaxy with an estimated BH mass of log (M_BH/ M_☉) = 7.46 (Bian & Gu 2007), however, the nuclear star-forming regions in this galaxy are unlikely to be associated with gaseous nuclear spirals studied in this work.

Finally, we discuss the mass inflow rates derived in our models in regard to powering AGNs in Seyfert galaxies. The mass accretion rate is often measured by the Eddington ratio defined by $\lambda \equiv L_{\rm bol}/L_{\rm Edd} = 4.5\times 10^{-2} (\epsilon /0.1) (\dot{M}/10^{-2}\; {M}_{\odot }\; {\rm yr}^{-1})(M_{\rm BH}/10^7\; {M}_{\odot })^{-1}$, where L_bol and L_Edd denote the bolometric and Eddington luminosities of an AGN and is the mass-to-energy conversion efficiency of the accreted material. Observations indicate that λ ≲ 0.1 for classical Seyfert 1 galaxies with broad iron emission lines (e.g., Meyer-Hofmeister & Meyer 2011; see also Peterson 1997), while λ ∼ 10⁻³ for low-luminosity Seyfert 1 AGNs (e.g., Ho 2008). In our numerical models, the mass inflow rates are larger for models with smaller M_BH and larger c_s. Taking ≈ 0.1 (e.g., Yu & Tremaine 2002), the mass inflow rates amount to λ ≲ 10⁻⁴ for M_BH = 4 × 10⁸ M_☉ and c_s ≲ 15 km s⁻¹, λ ∼ 10⁻³ for M_BH = 4 × 10⁸ M_☉ and c_s = 20 km s⁻¹ or for M_BH = 4 × 10⁷ M_☉ and c_s = 10 km s⁻¹, and λ ∼ 0.02–0.4 for M_BH = 4 × 10⁷ M_☉ and c_s ≳ 15 km s⁻¹. For classical Seyfert galaxies with strong bars, this suggests that the masses of central BHs are likely to be less than 10⁸ M_☉, which appears to be consistent with the measured values from the relation between the BH masses and the velocity dispersions of stellar bulges (e.g., Watabe et al. 2008) and reverberation mapping techniques (e.g., Gültekin et al. 2009; Denney et al. 2010). Of course, this result may depend on many factors such as the axis ratio and strength of the bar, presence of self-gravity and magnetic fields, gas cooling and heating, turbulence, etc., all of which would affect gas dynamic significantly. Extending the present work to include these physical ingredients would be an important direction of future research.

In this appendix, we describe the gravitational potential of the model galaxy that maintains the rotation of the (non-self-gravitating) gas disk. The potential is comprised of four components: a stellar disk, a bulge, a bar, and a BH. For the disk, we take a Kuzmin–Toomre model with surface density

$$\begin{equation} \sigma (R) = \frac{v_0^2}{2\pi G R_0} \left(1 + \frac{R^2}{R_0^2} \right)^{-3/2}, \end{equation} \tag{ A1 }$$

where v₀ and R₀ are constants (Kuzmin 1956; Toomre 1963). The corresponding gravitational potential at the disk midplane is

$$\begin{equation} \Phi _{\rm disk} = - \frac{v_0^2R_0}{\sqrt{R^2 + R_0^2}}, \end{equation} \tag{ A2 }$$

(Binney & Tremaine 2008). We take v₀ = 260 km s⁻¹ and R₀ = 14.1 kpc in our simulations.¹⁵ The total disk mass is M_disk = v²₀R₀/G = 2.2 × 10¹¹ M_☉.

For the bulge, we use a modified Hubble profile with volume density

$$\begin{equation} \rho (R) = \rho _{\rm bul}\left(1 + \frac{R^2}{R_b^2}\right)^{-3/2}, \end{equation} \tag{ A3 }$$

and the potential

$$\begin{equation} \Phi _{\rm bul} = - \frac{4\pi G \rho _{\rm bul}R_b^3}{R} \ln \left(\frac{R}{R_b} + \sqrt{1 + \frac{R^2}{R_b^2}}\right), \end{equation} \tag{ A4 }$$

where ρ_bul and R_b are the central density and the characteristic size of the bulge, respectively. We take ρ_bul = 2.4 × 10¹⁰ M_☉ kpc⁻³ and R_b = 0.33 kpc, with the corresponding bulge mass of M_bul = 2.8 × 10¹⁰ M_☉ within R = 6 kpc.

The bar is modeled by a Ferrers (1887) ellipsoid with volume density

$$\begin{eqnarray} \rho = \left\lbrace \begin{array}{@{}ll}\rho _{\rm bar}(1-g^2)^n & {\rm for}\ g<1, \\\\ 0 & {\rm elsewhere}, \end{array} \right. \end{eqnarray} \tag{ A5 }$$

where ρ_bar is the bar central density, g² = y²/a² + (x² + z²)/b² with (x, y, z) being the Cartesian coordinates, and a and b are the semimajor and semiminor axes of the bar, respectively. The exponent n measures the central concentration of the bar density distribution. In this work, we fix the bar parameters to n = 1, a = 5 kpc, b = 2 kpc, and ρ_bar = 4.5 × 10⁸ M_☉ kpc⁻³ for all models. The total mass of the bar is then M_bar = 2^{2n + 3}πab²ρ_barΓ(n  +  1)Γ(n  +  2)/Γ(2n  +  4) = 1.5 × 10¹⁰ M_☉ and the bar quadrupole moment is Q_m = M_bara²[1 − (b/a)²]/(5  +  2n) = 4.5 × 10¹⁰ M_☉ kpc², corresponding to a strong bar. When n = 1, the bar potential is given explicitly as

$$\begin{eqnarray} \Phi _{\rm bar}(x,y,z) &=& - \frac{\pi G ab^2 \rho _{\rm bar}}{2}\nonumber\\ &&\times\, [ W_{000} + x^2(x^2W_{200} + 2y^2W_{110}-2W_{100})\nonumber \\ && +\, y^2(y^2W_{020} + 2z^2W_{011}-2W_{010})\nonumber \\ && +\, z^2(z^2W_{002} + 2x^2W_{101}-2W_{001})], \end{eqnarray} \tag{ A6 }$$

where the coefficients W_ijk's are as defined in Pfenniger (1984). As the bar pattern speed, we choose Ω_b = 33 km s⁻¹ kpc⁻¹ which places the CR at R_CR = 6 kpc.

Finally, for a central BH with mass M_BH, we use a Plummer potential

$$\begin{equation} \Phi _{\rm BH} = - \frac{GM_{\rm BH}}{\left(R^2 + R_s^2\right)^{1/2}}, \end{equation} \tag{ A7 }$$

with the softening radius R_s = 1 pc. We take M_BH = 0, 4 × 10⁷ M_☉, and 4 × 10⁸ M_☉ to study the effects of the central mass concentration on the bar substructures. With M_BH ≪ M_disk, M_bul, M_bar, the BH affects the rotation curve only in the very central regions (with R ≲ 0.5 kpc). Figures 1 and 2 plot the rotational velocity for models with M_BH = 4 × 10⁸ M_☉ and the angular frequency curves for all models, respectively.