Early-growing Supermassive Black Holes Strengthen Bars and Boxy/Peanut Bulges

We use the grid-based N-body simulation package GALAXY ⁶ (Sellwood 2014) to simulate the growth of point masses representing an SMBH at various stages in the formation of the bar. All models discussed in this paper were evolved from initial conditions used in previous works (Debattista et al. 2017, 2020; Anderson et al. 2022) and were generated using GALACTICS (Kuijken & Dubinski 1995; Widrow & Dubinski 2005; Widrow et al. 2008). These models consist of exponential disks within a modified Navarro–Frenk–White (NFW; Navarro et al. 1996) live dark matter halo. For all initial conditions the spherical halo density distribution ρ(r) is described by

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{2}^{2-\gamma }{\sigma }_{h}^{2}}{4\pi {a}_{h}^{2}}\displaystyle \frac{C(r)}{{\left(r/{a}_{h}\right)}^{\gamma }{\left(1+r/{a}_{h}\right)}^{3-\gamma }},\end{eqnarray} \tag{ 1 }$$

where σ_h characterizes the halo velocity dispersion, and a_h characterizes the halo scale radius. C(r) is a function that smoothly truncates the model at a finite radius (Widrow et al. 2008):

$$\begin{eqnarray}&&C(r)=\displaystyle \frac{1}{2}\mathrm{erfc}\left(\displaystyle \frac{r-{r}_{h}}{\sqrt{2}\delta {r}_{h}}\right).\end{eqnarray} \tag{ 2 }$$

When γ = 1 and r_h → ∞, Equation (1) is exactly the NFW distribution. For all models considered in this paper, the halo parameters are σ_h = 400 km s⁻¹, a_h = 16.7 kpc, γ = 0.873, r_h = 100 kpc, and δ r_h = 25 kpc (Debattista et al. 2020). The live dark matter halo consists of 4 × 10⁶ particles of mass ≃1.7 × 10⁵ M_⊙ each.

The disk has an exponential radial density profile and an isothermal vertical density profile described by

$$\begin{eqnarray}&&{\rm{\Sigma }}(R,z)={{\rm{\Sigma }}}_{0}\exp (-R/{R}_{d})\,\ {{\rm{sech}} }^{2}(z/{z}_{d})\end{eqnarray} \tag{ 3 }$$

where R_d (z_d ) is the disk scale length (height), set to 2.4 kpc (0.3 kpc; Debattista et al. 2020). The initial kinematics of the disk are such that the radial velocity dispersion decreases exponentially as

$$\begin{eqnarray}&&{\sigma }^{2}(R)={\sigma }_{R0}^{2}\exp (-R/{R}_{\sigma }),\end{eqnarray} \tag{ 4 }$$

in which R_σ is fixed at 2.5 kpc and σ_R0 is the disk's central radial velocity dispersion (Debattista et al. 2020). The disk consists of 6 × 10⁶ equal mass particles contributing to a total disk ⁷ mass of ≃5.37 × 10¹⁰ M_⊙.

We assign the disk particles to a cylindrical polar grid nested within a larger spherical grid to which the halo particles are assigned. These grids share an origin that is relocated to the disk particle centroid at regular intervals to ensure the greatest spatial resolution at the central region of greatest density. This model does not contain a classical bulge component, and only forms a boxy or peanut-shaped bulge as a result of secular evolution.

The results in the main body of this work are evolved from the initial conditions for Model 2 from Debattista et al. (2020), for which σ_R0 = 128 km s⁻¹, and which is our principal control model (Model C). In order to evaluate the robustness of our results against stochastic effects, in Section 3.3 we present results from three azimuthally scrambled versions of Model C, and in Appendix A we also briefly consider results from two additional models from Debattista et al. (2020)—their Model 1 and Model 3 (also known, respectively, as D5 and D2 in Debattista et al. 2017). These models differ by having a value of σ_R0 = 90 km s⁻¹ and σ_R0 = 165 km s⁻¹, respectively, with all other parameters identical to the control model.

The SMBH is represented as a Plummer potential using a built-in function in GALAXY and is initialized with a small nonzero SMBH mass. The mass of the SMBH ( M_BH) grows as a function of time to a final mass M_fin according to the equation:

$$\begin{eqnarray}\,{M}_{\mathrm{BH}}(t)=\left\{\begin{array}{ll}\,{M}_{\mathrm{fin}}(0.02+0.98{\sin }^{2}(\pi \tau /2) & 0\leqslant \tau \leqslant 1\\ \,{M}_{\mathrm{fin}} & \tau \gt 1\end{array}\right.\end{eqnarray} \tag{ 5 }$$

where τ = t/ t_grow, and t_grow is the timescale over which the SMBH grows.

Since this is a pure N-body simulation, accretion by the SMBH is not modeled. Rather, our SMBH grows strictly due to an artificial deepening of the Plummer potential as given by Equation (5). The final mass of the SMBH in all cases is 0.0014 M_disk, which is comparable to the black hole mass fraction in M31 (Bender et al. 2005; Tamm et al. 2012). Since the SMBH potential is free to move due to accelerations from other particles, we introduce it with an initial mass of 0.02 · M_fin. Through testing we have found this mass (equal to 168 times the disk particle mass) is the minimum necessary to reduce early random motion of the SMBH and is required to ensure the SMBH does not accelerate away from the galaxy center while at ∼0 mass.

To introduce the SMBH into a model at a specific time, we use a snapshot of Model 2 evolved with no SMBH until the desired time as the initial conditions for a new model. Our SMBH potential is then added and evolved as in Equation (5). In this way we may grow an SMBH at a known stage in the evolution of a model by branching a new model off from the case where no SMBH is present.

In all models we keep the final mass (M_fin), the initial mass, the growth period (t_grow), and Plummer potential softening length fixed (see values listed in Table 1) and only vary the time of introduction of the SMBH. Previous works have shown that long-term effects on measurable bar quantities are fairly independent of t_grow, but depend strongly on M_fin and softening length (Shen & Sellwood 2004). Since our tests showed similar results, we consider a single growth period of t_grow = 50 dynamical times in simulation units, equivalent to 378 Myr.

We use an SMBH softening length ε = 33.33 pc, or two-thirds of the global ε = 50 pc (for both disk and halo particles). ⁸ This is more diffuse than the most compact CMCs meant to represent an SMBH with ε ∼ few parsecs, such as in Shen & Sellwood (2004), but still much more compact than simulated gas concentrations and star clusters (Shen & Sellwood 2004; Athanassoula et al. 2005; Sellwood & Gerhard 2020). This choice was motivated primarily by computational considerations; reducing ε for the given mass greatly increases the time resolution requirements in the vicinity of the SMBH. We choose a base time resolution of 0.03784 Myr (0.005 dynamic times) for all models. At larger radii, particle motion is calculated every 2^N time steps for N zones beginning at 2.4, 7.2, 12, and 19.2 kpc (multiples of scale radii). The criterion of Shen & Sellwood (2004) demonstrates that our chosen softening length is large enough that a circular orbit on the scale of ε will be sufficiently resolved, i.e., an orbital period will take ≳100 time steps. Increasing base time resolution to the extent required to reduce the SMBH softening length by a factor of 5–10 was prohibitively expensive. ^{9 , 10}

Simulation parameters required to evolve these models using GALAXY version 15.4, including full details of grid structure, are included in Appendix B. Simulation snapshots used for analysis were saved every 800 time steps; therefore, the time resolution for results presented in plots is 30.3 Myr unless otherwise noted.

The results presented in most of this paper are based on a small handful of models that we refer to as the principal set, detailed in Table 2. Each represents a scenario in which an SMBH is introduced into a snapshot of the control model, Model C, which is evolved from the initial conditions of Debattista et al. (2020) Model 2 with no SMBH. All other models in Table 2 are grown from snapshots of Model C with the SMBH introduced at the indicated time, and advanced forward to reach the same duration of total evolution, T_final = 7.568 Gyr.

In Model BF₀, the SMBH is introduced at t = 0, before bar formation i.e., before the m = 2 amplitude begins to increase from zero (SMBH growth is completed by t = 0.387 Gyr). The three models in which the SMBH is introduced before Model C experiences bar buckling are labeled Model BB₁, BB₂, and BB₃. The buckling time is considered to be the time of peak buckling amplitude (defined below).

These three models branch from Model C at t = 0.575, 1.150, and 1.877 Gyr, respectively (corresponding to convenient times in internal simulation units). Finally, to compare with previous work (e.g., Hasan et al. 1993; Norman et al. 1996; Shen & Sellwood 2004; Athanassoula et al. 2005; Hozumi & Hernquist 2005; Du et al. 2017) in which the SMBH was grown after the bar had reached a quasi-equilibrium state, we also ran a model with the SMBH introduced in Model C at t = 3.784 Gyr, after buckling has occurred (Model C buckles at 2.85 Gyr), which we call Model AB₁. Section 3.3 and Appendix A describe additional models that we explored to assess the generality of our results.

We first examine the effect of SMBH growth on the bar strength (quantified by bar amplitude A_{m = 2}) both prior to buckling (when the bar is forming and growing) and at late times. Figure 1 shows A_{m = 2} (top) and buckling amplitude A_buck (bottom) as a function of time. Each curve corresponds to a different model from the principal set in Table 2. Model C (blue) is the control model without an SMBH. The other curves show models with a growing SMBH, and the vertical bands show the period during which the SMBH grows in the model with the curve of the same color. In all models, bar buckling is characterized by a sharp drop in A_{m = 2} in the top panel and a corresponding spike in A_buck in the lower panel.

Zoom In Zoom Out Reset image size

The early rise and drop in A_{m = 2} prior to ∼1.1 Gyr is associated with formation then dissolution of a strong m = 2 spiral pattern. Although a weak bar first forms during this time, its contribution to the measured m = 2 mode is subordinate to the spiral. During this time no buckling occurs (despite the drop in A_{m = 2}), as is evident from the fact that A_buck ∼ 0. After the drop at t ∼ 1.1 Gyr, the bar becomes the dominant m = 2 structure, and remains dominant for the remainder of the time.

Model C (no SMBH) buckles the earliest with the lowest bar amplitude prior to buckling. The lower panel also shows that this model has the greatest buckling amplitude. After buckling, its bar ceases growth and quickly reaches an approximately steady state in bar amplitude that lasts throughout the simulation period. In contrast, Model BF₀ (orange), where the SMBH starts to grow at t = 0, buckles later than all other models with a much lower spike in A_buck than Model C. Model BF₀ has the largest A_{m = 2} at late times. This difference in bar amplitude is also apparent shortly after buckling, before any further bar growth occurs.

Figure 1 (top) shows similar trends for the other models with SMBHs introduced before buckling: SMBHs that start growing later (but before bar buckling) have smaller late-time bar amplitudes than Model BF₀, but all have higher A_{m = 2} than Model C. Consistent with previous work (e.g., Shen & Sellwood 2004; Athanassoula et al. 2005; Hozumi & Hernquist 2005), Model AB₁ (brown curve) in which the SMBH is grown after bar buckling, when bar growth has ceased, is the only model in which A_{m = 2} at late times is lower than in Model C.

Figure 1 (bottom) shows that in cases in which the SMBH is introduced prior to bar buckling, the buckling event is partially suppressed. We define "suppression" as either the reduction of buckling amplitude (A_buck), a delay of the buckling event (location of the peak in A_buck), or both, relative to Model C. In Models BF₀, BB₁, and BB₃, buckling is both weakened and delayed. Model BB₂ produces a strong buckling event, almost as strong as the buckling experienced by Model C; however, it is still significantly delayed compared to Model C such that the final bar amplitude at late times is still significantly larger than in Model C. We assert that since buckling weakens the bar, the delay in buckling allows the models to develop a stronger bar prior to buckling. Simultaneously a generally weaker buckling event is associated with a smaller decrease in bar amplitude during buckling. Therefore, both effects contribute independently to a greater final bar strength than in Model C. The reasons for this suppression of buckling are explored in depth in Section 4.

Finally we draw attention to the long-term A_buck behavior of Model BF₀, which after buckling remains elevated relative to other models and gradually increases with time for the duration of the simulation (Figure 1 lower panel). This is indicative of a persistent m = 2 asymmetry about the midplane (in both the bar and disk). In other words, the bar and disk remain bent throughout the duration of the simulation rather than only bending out of plane for a short period during buckling. This is similar to behavior encountered, but typically not elaborated upon, in multiple other works (e.g., Gardner et al. 2014; Li & Shen 2015; Debattista et al. 2017; Smirnov & Sotnikova 2018), and is poorly understood. Cuomo et al. (2023) carried out a detailed study of such long-term asymmetry about the midplane in both simulations and observations. We do not investigate this asymmetry further in this paper, although we note that similar long-term bending was also observed in a few other simulations with other initial conditions that we simulated (see, e.g., Figures 2 and 3). This long-term bending is seen in simulations with and without SMBH and has also been found by other authors (all without SMBH).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Sellwood & Debattista (2009) have shown that large N-body systems are inherently chaotic and thus susceptible to stochastic effects, where even small differences in the initial conditions, even for the same model parameters can build to significant differences late in the simulation (affecting quantities such as bar growth rates and susceptibility to repeated buckling). These authors assert that all outcomes are, however, equally valid representations of the evolution of the system. In the previous subsections, Model C was evolved from the same set of initial conditions as Model 2 of Debattista et al. (2020), which are evolved for a longer period by Anderson et al. (2022). An observant reader will notice that while the bar amplitude of Model C does not increase after the buckling event at t = 2.85 Gyr, the bar amplitude in Model 2 of Anderson et al. (2022; which first buckles at a similar time t = 2.9 Gyr) continues to grow and eventually buckles a second time. We draw attention to this difference as a possible example of stochasticity, where the two models diverge in behavior gradually and relatively late in their evolution (in this case, a small difference in how hot the disk is at the time after buckling). Sellwood & Debattista (2009) demonstrated that such effects can grow from differences as minute as a difference in memory allocation between two processors or a change in the order in which particle coordinates are listed in the initial conditions, with all other factors being equal. Therefore, a late time divergence of model behaviors (even beginning from exactly the same initial conditions) evolved using entirely different simulation codes, such as GALAXY by us and PKDGRAV (Stadel 2001) by Anderson et al. (2022) is virtually unavoidable. This was already demonstrated by Sellwood & Debattista (2009) in a particularly relevant direct comparison of initial conditions evolved using GALAXY and PKDGRAV.

As a further test for sensitivity to stochastic effects with the same code (GALAXY) and same parameters, we ran two additional tests starting with the initial conditions for Model C presented in the previous two subsections.

• 1.  
  We use a feature provided by the GALAXY code that allows various sectoral harmonics of the Fourier expansion of the potential to be turned on/off, thereby turning on/off particular forces in calculation and allowing for enforcement of certain symmetries. We disabled all modes other than m = 0, and evolved Model C for 48 dynamical times (363 Myr). In this case, we did not grow an SMBH during the time the higher Fourier modes were turned off. This model effectively resulted in a new set of axisymmetric equilibrium initial conditions with overall model parameters identical to Model C. We note that the suppression of all modes except m = 0 results in a disk with somewhat lower radial velocity dispersion than in Model C, but with otherwise identical properties. This simulation was treated as a separate set of initial conditions, and evolved further without an SMBH (Model $C^{\prime}$) and with an SMBH inserted at t = 0 analogous to BF₀ (Model ${{BF}}_{0}^{\prime}$). Since Model $C^{\prime}$ has smaller radial velocity dispersion than Model C, the bar strength initially increases more slowly than in Model C and buckles significantly later than in Model C (see orange curve in Figure 2). As can be seen in Figure 2, the bar in Model $C^{\prime}$ continues to grow in amplitude after buckling (orange curve) unlike Model C (blue curve). In Model ${{BF}}_{0}^{\prime}$ (not shown), the final bar strength was stronger than the bar in Model $C^{\prime}$ at the final time step and is also stronger than for Model BF₀.
• 2.  
  We azimuthally scrambled the initial conditions of Model C, creating two more axisymmetric models by generating random azimuthal angles for all particles (and appropriately rotating their Cartesian planar velocities). The two models dubbed ${\widetilde{C}}_{i}$ and ${\widetilde{C}}_{{ii}}$ were created by the same process, with the only difference being the starting value of the random seed used to scramble (see green and red curves in Figure 2). As can be seen in this figure both these scrambled initial conditions buckle slightly later than Model C (but before Model $C^{\prime}$). In particular Model ${\widetilde{C}}_{{ii}}$ buckles only weakly and ends up with a higher bar amplitude than any of the other Model C analogs.

We also studied the effect of both early-growing black holes (BF₀ analogs) and black holes introduced after bar buckling (AB₁ analogs) on the final bar strength. Since the scrambled control models (${\widetilde{C}}_{i},{\widetilde{C}}_{{ii}}$) do not reach steady state (unlike Model C), for the AB₁ analogs we introduce the SMBH the same number of dynamical time units after bar buckling as we do for model AB₁. Figure 3 shows Model C analogs (${\widetilde{C}}_{i},{\widetilde{C}}_{{ii}}$), Model BF₀ analogs (${\widetilde{{BF}}}_{0i},{\widetilde{{BF}}}_{0{ii}}$), and Model AB₁ analogs (${\widetilde{{AB}}}_{1i},{\widetilde{{AB}}}_{1{ii}}$). To avoid crowding, we do not show models $C^{\prime}$, ${{BF}}_{0}^{\prime}$, and ${{AB}}_{1}^{\prime}$ but the latter two models show a similar degree of stochasticity to the models shown in Figure 3.

Both ${\widetilde{{BF}}}_{0i}$ and ${\widetilde{{BF}}}_{0{ii}}$ buckle later and more weakly than their corresponding control models. As in BF₀, the bar amplitudes in ${\widetilde{{BF}}}_{0i}$ and ${\widetilde{{BF}}}_{0{ii}}$ continue to grow after buckling. Since Model ${\widetilde{{BF}}}_{0i}$ buckles later than Model BF₀, it does not attain as high an amplitude as BF₀ but it still has a stronger bar than its counterpart without an SMBH (Model ${\widetilde{C}}_{i}$). In contrast, both Model ${\widetilde{C}}_{{ii}}$(red) and Model ${\widetilde{{BF}}}_{0{ii}}$ (purple) buckle so weakly that they have similar bar strengths at the end of the simulation. However, even in this case, the buckling amplitude is lower for the model with the SMBH, and buckling is delayed relative to the control, confirming our findings in Section 3.2.

The Model AB₁ analogs also show interesting results. In the case of ${\widetilde{{AB}}}_{1i}$ (green), the bar is slightly weakened relative to the control model (${\widetilde{C}}_{i}$) (blue). However, introduction of the SMBH does not stop bar growth; therefore, it does not weaken or destabilize the bar in absolute terms. In the case of ${\widetilde{{AB}}}_{1{ii}}$ (brown), the final bar strength is identical to that of the control model ${\widetilde{C}}_{{ii}}$ (red), which, as we saw in Figure 2, has the strongest bar amplitude of all of the control models.

There is clearly a great deal of stochasticity in the models we have presented in this section. However, we can confidently state that an SMBH of ≲0.2% of disk mass that grows prior to bar formation never destabilizes or weakens a bar. The theoretical analysis presented in Sections 4 and 5 sheds light on the empirical results found in this section. Although an SMBH introduced after buckling might slightly weaken a bar relative to a case with no SMBH, it does not always do so, especially if the bar is still growing (as we see for Model ${\widetilde{{AB}}}_{1{ii}}$).

Since the previous generation of simulations all introduced an SMBH after the bar had attained a steady state (e.g., Hasan et al. 1993; Norman et al. 1996; Shen & Sellwood 2004; Athanassoula et al. 2005; Hozumi & Hernquist 2005; Du et al. 2017) and since our principal control Model C is the only one that attains a steady state after buckling, in the remainder of this paper we will primarily focus on Model C and the models derived from it that are listed in Table 2.

We now examine how the growth of an SMBH affects the strength of the BP/X-shaped bulges that form in our simulations. Figure 4 shows x–z projections of the surface density of final snapshots of four models. The BP/X-shaped bulge is evident, and the models with stronger bars (Models BF₀, BB₁) clearly show stronger BP/X shapes than models with weaker bars (Models C, AB₁).

Zoom In Zoom Out Reset image size

We parameterize the shape of the BP/X bulges using a new bar-deprojection method developed by Dattathri et al. (2023), which we briefly describe here. We construct parametric 3D bar models that are added as additional modules to the surface-brightness fitting routine IMFIT (Erwin 2015). IMFIT is able to take any input 3D distribution and project it at an arbitrary orientation to derive the best-fit parameters of the 3D model required to fit the image. In the new BP/X bar fitting module the bar is assumed to have a sech² profile in a scaled radius given by

$$\begin{eqnarray}&&\rho ={\rho }_{0}\,{{\rm{sech}} }^{2}(-{R}_{s}),\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{R}_{s}={\left({\left[{\left(\displaystyle \frac{x}{{X}_{\mathrm{bar}}}\right)}^{{c}_{\perp }}+{\left(\displaystyle \frac{y}{{Y}_{\mathrm{bar}}}\right)}^{{c}_{\perp }}\right]}^{{c}_{\parallel }/{c}_{\perp }}+{\left(\displaystyle \frac{z}{{Z}_{\mathrm{bar}}}\right)}^{{c}_{\parallel }}\right)}^{1/{c}_{\parallel }}.\end{eqnarray} \tag{ 9 }$$

The origin of the coordinate system is the galaxy center; X_bar, Z_bar, and Y_bar are the semiaxis lengths of the bar along the long-axis (x), the axis perpendicular to the disk plane (z), and y is the other axis in the disk plane (these describe an ellipsoidal bar). c_∥ and c_⊥ control the diskiness/boxiness of the bar (the 3D analog of the parameters proposed by Athanassoula et al. 1990). While Dattathri et al. (2023) explored a range of values for c_∥ and c_⊥, here we set c_∥ = c_⊥ = 2 (corresponding to an ellipsoidal bar). The BP/X shape is determined by a scale height perpendicular to the disk that depends on the position in the x − y plane represented by a double Gaussian centered at the galactic center:

$$\begin{eqnarray}\begin{array}{rcl}{Z}_{\mathrm{bar}}(x,y) & = & {A}_{\mathrm{pea}}\exp \left(-\displaystyle \frac{{\left(x-{R}_{\mathrm{pea}}\right)}^{2}}{2{\sigma }_{\mathrm{pea}}^{2}}-\displaystyle \frac{{y}^{2}}{2{\sigma }_{\mathrm{pea}}^{2}}\right)\\ & & +\,{A}_{\mathrm{pea}}\exp \left(-\displaystyle \frac{{\left(x+{R}_{\mathrm{pea}}\right)}^{2}}{2{\sigma }_{\mathrm{pea}}^{2}}-\displaystyle \frac{{y}^{2}}{2{\sigma }_{\mathrm{pea}}^{2}}\right)+{z}_{0},\end{array}\end{eqnarray} \tag{ 10 }$$

where R_pea and σ_pea are the distance of the center of the peanuts from the galaxy center and the width of each peanut, respectively. A_pea measures the vertical extent of the peanut above the ellipsoidal bar z₀ (not the distance from the midplane), since z₀ is the base vertical scale height of the bar. The maximum distance of the peanut from the midplane is therefore given by

$$\begin{eqnarray}&&{h}_{\mathrm{pea}}\approx {z}_{0}+{A}_{\mathrm{pea}},\end{eqnarray} \tag{ 11 }$$

which we refer to as the "peanut height," which is the maximum value of Z_bar. We note that in our model the peanut is entirely associated with the bar, and there is no separate bulge component.

Equation (10) is similar to the "peanut height function" proposed previously by Fragkoudi et al. (2015), with the additional assumption that both halves of the peanut are symmetric about the origin, and the peanuts are aligned with the major axis of the bar.

The shape of the bar is controlled primarily by three parameters: R_pea, h_pea, and σ_pea. Dattathri et al. (2023) demonstrated that these three parameters offer great versatility for modeling a variety of shapes for the BP/X feature. They also show, by fitting mock 2D surface-brightness data generated from some of the N-body models in this work, that the 3D density distribution and potential arising from this parameterization of the bar and BP/X bulge closely match the N-body density and potentials. Dattathri et al. (2023) also showed that the BP/X shape parameters can also be obtained for an N-body model by fitting Equation (10) directly to the stellar particle distribution in order to obtain the "true" values of R_pea, h_pea, and σ_pea instead of going through the deprojection algorithm. Dattathri et al. (2023) found that this 3D fit to the N-body particle distribution gives a better representation of the gravitational potential.

Here we use this direct 3D fitting to the N-body stellar particle distribution to quantify the strength of the BP/X bulge at the final time step of 12 different models. The key parameters that determine the BP/X shape: R_pea, h_pea, and σ_pea, are recovered for each model at the final time step.

Figure 5 shows the values of the parameters R_pea and h_pea obtained via 3D fitting as a function of the bar amplitude of the model (σ_pea shows no correlation with bar amplitude and hence is not shown). It is clear that the strength of the peanuts as parameterized by h_pea and R_pea are strongly correlated with bar amplitude. This strong correlation between the bar amplitude and the parameters (R_pea and h_pea) characterizing the BP/X shape does not depend on whether or not a model includes a black hole (the points with black edges are control models while the others have an SMBH). The cause of this correlation is investigated in Section 6. While further simulations are needed to confirm this correlation in hydrodynamical simulations, the existence of these correlations has an important implication. For edge-on disks in which the m = 2 bar amplitude is impossible to measure by traditional means, the use of the BP/X fitting function in Equation (10) could enable a determination of bar strength (if the bar major axis is 25°–90° to the line of sight). This will be investigated further in the future.

Zoom In Zoom Out Reset image size

We use the AGAMA package (Vasiliev 2019) to estimate the frequencies assuming a triaxial potential for the star particles and an axisymmetric potential for the halo. Then, at each radius we compute the average between the frequencies along the major and minor axes of the bar.

In order to estimate the approximate location of the vILR and ILR, we compute the pattern speed of the bar at various times. This is done by fitting a straight line to the time evolution of the phase angle of the bar (assumed to increase monotonically with time) for 10 snapshots around the snapshot of interest. The obtained pattern speeds are shown as horizontal dotted lines in each panel of Figure 7 (with lengths corresponding to the bar lengths in that snapshot). We compute the approximate location of the vILR, ILR, and corotation resonances as the radii where the curves for Ω(R) (thin solid), Ω(R) − κ(R)/2 (thick solid), and Ω(R) − ν(R)/2 (dashed–dotted) intersect the average bar pattern speed Ω_p around that snapshot—see Figure 7.

Zoom In Zoom Out Reset image size

At t = 0 all of the models are identical, so only Model C is shown. The next four panels show the three simulations at additional snapshots. The bottom-right panel shows the approximate location of the ILR and vILR as a function of time in each of the three models.

At t = 1.51 Gyr, in Model BF₀ (orange), Ω rises sharply at small radii (a consequence of the early-growing SMBH and an increase in the stellar density around it) and both Ω(R) − κ(R)/2 and Ω(R) − ν(R)/2 have already begun to increase relative to Model C. At t = 2.27 Gyr, both the ILR and vILR are present at small radii in the model with an early SMBH but not in Model C. This is because of the combined effect of the increase in Ω and because the SMBH has already started heating the disk vertically, increasing σ_z , and therefore decreasing ν enough to cause Ω − ν/2 to increase. These two factors cause the appearance of both the ILR and vILR at small radii, nearly 0.5 Gyr prior to the appearance of these two resonances in Model C (see bottom-right panel). In model BF₀, the ILR and vILR appear at R ≲ 1 kpc by around t ≃ 2 Gyr, and their radial location moves rapidly outward between first appearance and the final snapshots (where it is beyond 5 kpc). In Model C (blue), these resonances appear later (t = 2.5–3.0 Gyr), and these resonances stay within R ∼ 3.5 kpc. Model AB₁ is only shown after the SMBH finishes growing (after t = 4.232 Gyr). In this model, both ILR and vILR have R < 3.5 kpc throughout the evolution, similar to Model C. The resonance locations in both models C and AB₁ show small oscillations around a constant value after ∼3.5 Gyr, which we verified as being due to small oscillations in the estimated pattern speeds.

We note that at t = 1.51 Gyr, the vILR is not present in any model. By t = 4.54 Gyr, all models have buckled and the vILR is present in all models thereafter, but it is clear that the SMBH, despite having a small nominal sphere of influence (r_BH ≲ 0.175 kpc), is able to cause vertical heating in the bar for stars that travel to much larger radii due to the high radial anisotropy of the bar. The bottom-right panel suggests that the vILR does not appear until the bar buckles in Model C, but it appears far prior to bar buckling in Model BF₀. Although not shown, we note that in Model BB₁ too, the vILR appears before the bar buckles, although it appears slightly later than in Model BF₀, as would be expected for BB models in relation to BF₀ from our previous results.

These results suggest that, since the vILR and ILR in Model BF₀ sweep over large radial range, they can resonantly capture a significantly larger fraction of disk stars than they do in Models C and AB₁, where they sweep a much smaller radial range.

We now briefly describe the orbital frequency analysis first introduced by Binney & Spergel (1982, 1984). The method was significantly improved and refined by Laskar (1990, 1993) who referred to their algorithm as "Numerical Analysis of Fundamental Frequencies (NAFF)."

In Hamiltonian dynamics, the angle variables and their canonically conjugate actions J_i uniquely define a regular orbit (Binney & Tremaine 2008). In a 3D potential, the J_i , i = 1,..,3 are conserved, and time derivatives of the angle variables ${{\rm{\Omega }}}_{i}=\dot{{\theta }_{i}}(t),i\,=\,1,..,3$ are also constants of motion. Since regular orbits in galaxies are quasiperiodic, their space and velocity coordinates can be represented by time series of the form:

$$\begin{eqnarray}&&x(t)=\displaystyle \sum _{k=1}^{{k}_{\max }}{A}_{k}{e}^{i{\omega }_{k}t}\end{eqnarray} \tag{ 12 }$$

with similar expressions for y(t), z(t) and velocity components, V_x (t), V_y (t), and V_z (t). The amplitudes A_k of the N (typically N = 10–15) largest peaks in the spectrum and their corresponding frequencies ω_k are obtained by taking a Fourier transform of a complex time series, e.g., f_x (t) = x(t) + iV_x (t) (and similarly for f_y , f_z ) constructed from the spatial and velocity coordinates of an orbit. This is followed by Gram–Schmidt orthogonalization to properly extract additional, lower-amplitude frequencies in the spectrum. Eventually, we obtain the basis set of three linearly independent frequencies Ω_i . All of the other frequencies ω_k in the three frequency spectra (one each for f_x , f_y , and f_z ) can be written as integer linear combinations of these three frequencies; therefore, the Ω_i , i = 1,..,3 are referred to as "fundamental frequencies."

Here we use our implementation ¹² of the NAFF algorithm to recover orbital fundamental frequencies (Valluri & Merritt 1998; Valluri et al. 2010) and to classify bar orbits (as described in Valluri et al. 2016). With this code, the frequency components in the spectrum can be recovered with high accuracy (1 part in 10⁵ or better) in ∼20–30 orbital periods.

Previous work has shown (Valluri et al. 2010, 2013, 2016) that when frequency analysis is applied to a large representative sample of orbits drawn from a self-consistent distribution function, the analysis of the frequency maps (plots showing ratios of orbital fundamental frequencies plotted against each other) and automated orbit classification enables a quantitative assessment of the relative importance of different resonances and different types of orbits to the phase-space structure of the galaxy. It also enables one to understand how the orbital distribution function is altered by evolution of the galactic potential.

Frequency analysis is also a powerful means to assess whether orbits are regular or chaotic. Since regular orbits conserve actions, their frequencies also remain constant in a static potential. By computing each of the three fundamental frequencies Ω_i over two separate time intervals T₁, T₂ one can compute

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{{\rm{\Omega }}}_{i}}=\left|\displaystyle \frac{{{\rm{\Omega }}}_{i}({T}_{1})-{{\rm{\Omega }}}_{i}({T}_{2})}{{{\rm{\Omega }}}_{i}({T}_{1})}\right|\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{\Omega }}}=\max ({{\rm{\Delta }}}_{{{\rm{\Omega }}}_{i}}),i=1,\ldots ,3\end{eqnarray} \tag{ 14 }$$

with Δ_Ω defined as the frequency drift. Therefore, frequency drift provides a measure of how chaotic (irregular) an orbit is.

For the analysis of the N-body simulations studied in this paper, we randomly select 10⁵ star particles from Model C in the initial snapshot (t = 0). We integrate the orbits of this same set of particles in all of the models at t = 1.51, 2.27, 3.03, 3.78, 4.54, and 7.57 Gyr (snapshots 200, 300, 400, 500, 600, and 1000, where the integers denote the number of dynamical times elapsed since the start of the simulation) in the respective potentials and in the presence of the rotating bars. The potentials are estimated with AGAMA, assuming a triaxial distribution for star particles, an axisymmetric distribution for DM particles and a Plummer potential for the SMBH (when one is present). While not all of these star particles eventually end up in the bar in the final snapshot of any specific model, we consider this selection to be sufficiently random and uniform to enable us to compare the results of orbital frequency analysis of the various models in an unbiased manner.

We first compute the orbital frequencies in cylindrical coordinates in the inertial frame of the galaxy, which allows us to properly identify the corotation resonance (Beraldo e Silva et al. 2023). We use the standard cylindrical polar coordinates R, v_R , ϕ, v_ϕ , z, v_z . R, v_R , and z, v_z are the canonically conjugate radial and vertical coordinates and momenta, respectively, and hence can be used to define a complex time series f_R (t) = R(t) + iv_R (t) and f_z (t) = z(t) + iv_z (t). However, since ϕ is an angular coordinate and v_ϕ is an angular velocity (not linear coordinate and linear momentum), this pair does not yield the correct tangential frequencies when used to construct the time series for frequency analysis. Instead we use the angular momentum L_z = xv_y − yv_x and the Poincaré symplectic polar variables $\sqrt{2| {L}_{z}| }\cos \phi$ and $\sqrt{2| {L}_{z}| }\sin \phi$ to define the complex time series ${f}_{\phi }(t)\,=\sqrt{2| {L}_{z}| }(\cos \phi +i\sin \phi )$ (Papaphilippou & Laskar 1996; Valluri et al. 2012). Beraldo e Silva et al. (2023) further discussed the advantages of these complex combinations.

Figure 8 shows histograms of (Ω_ϕ − Ω_P )/Ω_R (left) and (Ω_ϕ − Ω_P )/Ω_z (right) at different times, color coded by the mean ${z}_{\max }$ in each bin (with width 0.0025)—note the different y-axis scales. The two left columns refer to Model C. The main resonances are easily recognized from the sharp peaks: the corotation ((Ω_ϕ − Ω_P ) = 0), the ILR ((Ω_ϕ − Ω_P )/Ω_R = 0.5; see Athanassoula 2003), and the vILR ((Ω_ϕ − Ω_P )/Ω_z = 0.5).

Zoom In Zoom Out Reset image size

The ILR is known as the main resonance supporting bars (e.g., Contopoulos & Papayannopoulos 1980; Athanassoula 2003). The early significant development of this resonance in Model C (upper-left panel of Figure 8) agrees with the early development of the bar observed in Figure 1. It is interesting to note that this early development of the ILR in Model C was not detected in the simpler quasiaxisymmetric approximation—see Figure 7.

In the second column, we see that in Model C the vILR is almost unpopulated at 2.27 Gyr, but a significant number of orbits have crossed this resonance by 4.54 Gyr, i.e., have (Ω_ϕ − Ω_P )/Ω_z ≳ 0.5. Interestingly, the mean ${z}_{\max }$ increases abruptly for these orbits, in agreement with the theoretical expectation of the excitement of vertical motion by this resonance (Binney 1981; Pfenniger & Friedli 1991)—see also Beraldo e Silva et al. (2023).

The two right columns of Figure 8 show the histograms for the model BF₀. The ILR (left panels) is similarly populated over time, but it is more populated than Model C already at 1.51 Gyr, and the peak reaches significantly larger values at t = 4.54 Gyr. At the final snapshot (bottom row), we identify 54,934 orbits within Δ∣(Ω_ϕ − Ω_P )/Ω_R ∣ ≤ 0.01 from the ILR, in comparison to 28,482 identified in the same range for Model C, which agrees with the higher bar amplitude observed for model BF₀—see Figure 1.

The right panels show how the vILR is populated over time, as it is already detected in Model BF₀ at 2.27 Gyr, and that the mean ${z}_{\max }$ increases abruptly for orbits crossing this resonance (i.e., moving rightward in the histogram). In particular, we note that this resonance is very strongly populated at t = 4.54 Gyr and at t = 7.57 Gyr (much stronger than in Model C) and that at the final snapshot the mean ${z}_{\max }$ for orbits close to the resonance is substantially larger than in Model C, although the number of stars strictly at the resonance is smaller. The number of orbits that crossed the vILR is also significantly larger than in Model C. Taking into account that the BP/X bulge is stronger in this simulation, this suggests that after stars are released from (or cross) the vILR, they continue to support the BP/X. This seems to agree with the theoretical scenario described by Quillen et al. (2014), where stars do not stay trapped by the vILR for a long time, but crossing this resonance is enough to promote stars to orbits with high ${z}_{\max }$ supporting the BP/X shape—see also Sellwood & Gerhard (2020).

Finally, Figure 9 shows the histograms for Model AB₁. It is clear that at the final snapshot, both the ILR and the vILR are less populated than in Model C. The figure also suggests that the number of orbits that crossed the vILR and the mean ${z}_{\max }$ around the vILR are smaller than in the control Model C, in agreement with the weakness of the BP/X bulge in this simulation.

Zoom In Zoom Out Reset image size

In summary, these figures confirm the qualitative inferences drawn from the quasiaxisymmetric approximation in Figure 7, particularly the early development of the ILR and the vILR in Model BF₀, and that the effect of the ILR and vILR sweeping outward in the disk is to capture a significantly greater fraction of orbits into these resonances, strengthening the bar and causing greater levitation of orbits to high ∣z∣ where they support the BP/X bulge.

To further understand why an early-growing SMBH (e.g., in Model BF₀) can strengthen both the bar and the BP/X bulge, while the SMBH grown after bar buckling weakens a steady-state bar, and the specific resonant orbit families that support the BP/X shape, we now examine frequency maps in Cartesian coordinates.

Valluri et al. (2016) showed that in frequency maps in Cartesian coordinates, for orbits integrated in the rotating frame of the bar, bar-supporting orbits are largely clustered in a cloud with 0.45 < Ω_x /Ω_z < 0.75 and 0.65 < Ω_y /Ω_z < 0.95, approximately the same region occupied by box orbits in stationary (nonrotating) triaxial potentials. These orbits (in the frame rotating with the bar) primarily originate from bifurcations of the linear long-axis orbit, which is the parent orbit of the "box" orbit family in stationary triaxial potentials and are referred to as the x1-orbit family in bars (Valluri et al. 2016). The BP/X-shaped bulge is associated with several families of resonant orbits, although it is not populated exclusively by resonant orbits (Abbott et al. 2017).

The introduction of a central point mass (representing an SMBH) in a (static) triaxial potential causes an increase in the number of resonances populated by box-like orbits (Valluri & Merritt 1998). As the mass of the SMBH increases, some of the resonances on the frequency map become thicker (due to increased resonant capture) while others are broken up or "fractured" due to the increased overlap of resonances. It has been well known since the early work of Chirikov (1979) that increasing the strength of a perturbation in a potential increases the number of resonances and that resonant overlap is an important factor driving the increase in the fraction of chaotic orbits. These chaotic orbits undergo mixing that drives the potential to a new dynamical equilibrium with fewer chaotic orbits (Merritt & Valluri 1998).

In Figure 10 we show frequency maps for Models C, AB₁, and BF₀ at t = 4.54 Gyr and t = 7.57 Gyr constructed from the integration of the same sets of 10⁵ randomly selected star particles considered in Figure 8. The points on the frequency map are colored by the logarithm (base 10) of the pericenter radius of each orbit (over the integration time of 100 orbital periods).

Zoom In Zoom Out Reset image size

The top row of Figure 10 shows that the number and strength of most resonances in the control Model C do not change significantly between t = 4.54 Gyr and t = 7.57 Gyr (during which time the bar strength remains almost constant). In Model AB₁ (middle row) the frequency map at t = 4.54 Gyr (only 0.75 Gyr after the SMBH was introduced) already shows some differences from Model C, with several new weak resonance lines and increased clustering/scattering around the intersections of resonances. In addition, some strong resonances (e.g., (1, −2, 1) and (3, −5, 2)) appear "broken up" or "fractured" at the places where they intersect other resonances. By t = 7.57 Gyr it is clear from the colors of the points that all of the orbits in Model AB₁ have significantly higher pericenter radii than they do in Model C or even in the previous snapshot of Model AB₁ at t = 4.54 Gyr. Furthermore, as can be seen in Figure 12 at t = 7.57 Gyr, certain resonances, e.g., (1, −3, 2), (1, −2, 1), and (2, 0, −1) in Model AB₁ are depopulated relative to their appearance at t = 4.54 Gyr.

The bottom row shows frequency maps for Model BF₀ (SMBH introduced at t = 0). Recall that in this model the bar formed and grew around the fully grown preexisting SMBH. The frequency maps show that some resonances, e.g., (1, −3, −2), (2,0, −1) are more strongly populated than in Model C and some additional resonances (e.g., (3, 0, −5), (0, −1, 1)) have appeared. However, some strong resonances in the top two rows (e.g., (1, −2, 1) and (3, −5, 2)) are depopulated.

The overall structure of the frequency map in Model BF₀ is also quite different from that in the top two rows, reflecting the different orbit populations in this model. Since the bar forms and grows around a fully grown SMBH, the growing bar captures stars that then populate orbit families that are not destabilized by the presence of the central SMBH. The frequency map for Model BF₀ at t = 7.57 suggests that the bar-supporting resonant orbit families are populated in a different manner than in Model C. Additionally in Model BF₀ at t = 7.57 Gyr the diagonal resonance line ((1, −1, 0); mostly occupied by disk stars) is much thinner than in the other panels—further evidence that the bar, having captured a significantly larger fraction of disk orbits, is much stronger in this model than in any other.

Although we do not show it, the Cartesian frequency map for Model BB₁ at t = 7.57 Gyr appears intermediate between the maps for Models C and BF₀. Although the bar has already begun to form prior to SMBH introduction, its continued growth allows the stars to populate stable orbits.

Figure 11 shows histograms of the spherical pericenter radius $\mathrm{log}({r}_{\mathrm{per}})$ for the 10⁵ orbits shown in the frequency maps in Figure 10. At t = 4.54 (left), Models C and AB₁ differ very little in the distribution of $\mathrm{log}({r}_{\mathrm{per}})$, and all three models show a distinct peak of orbits at very small radii $\mathrm{log}({r}_{\mathrm{per}}/\mathrm{kpc})\lt -1$. Note that this range of radii is comparable to the nominal r_BH ∼ 0.1–0.175 kpc of the SMBH in Model AB₁ and Model BF₀. It is important to point out that Model BF₀ (with the SMBH grown at t = 0) has a slightly higher fraction of orbits at the smallest r_per than even Model C (at t = 4.54), which is evidence that a growing bar has absolutely no difficulty in finding ways to populate stable orbit families that pass quite close to the SMBH. Furthermore, the peak at the smallest r_per lies well inside the SMBH sphere of influence in Model BF₀ of 0.175 kpc.

Zoom In Zoom Out Reset image size

By t = 7.57 Gyr, the peak in r_per/kpc < −0.75 has decreased slightly, consistent with the expectation that these orbits were scattered by the SMBH to large radii. Despite the fact that Model BF₀ has an SMBH of exactly the same mass as AB₁, the peak at the smallest radii is more highly populated at t = 7.57 Gyr, and has significantly more stars in this innermost region than Model C (which has no SMBH).

The picture that emerges from the preceding analysis is that in Model AB₁, the SMBH that was introduced after the bar formed and reached equilibrium, scatters orbits with log(r_per/kpc) < −0.75 to larger pericenter radii. Presumably since many of these were bar-supporting, their scattering weakened the bar and, because the bar potential is no longer growing at this point, stars cannot be captured onto alternative stable bar-supporting orbits. However, when the bar forms and/or grows around the SMBH (e.g., in Model BF₀), scattered orbits as well as newly captured orbits succeed in populating alternative centrophobic resonant orbit families that are bar-supporting.

The quantitative differences in the occupancy of resonant bar-supporting orbits in Models C, AB₁, and BF₀ can be seen in Figure 12, which shows the number of orbits with frequency ratios within 5 × 10⁻³ of the major resonances that appeared in the frequency maps in Figure 10. Already at t = 4.54 Gyr, we see that the occupancy of some resonances (e.g., (1, −3, 2), (2, 0, −1)) is significantly greater in Model BF₀ compared to Models C, AB₁. Other resonances such as (1, −2, 1), (3, −5, 2), and (3, 0, −2) are depopulated. At t = 4.54 Gyr, the differences between Model C and AB₁ are fairly small, but we begin to see changes in the population of resonant orbits. By the end of the simulation (right panel), we see that in Model AB₁ the number of orbits associated with several resonances ((1, −3, 2), (1, −2, 1), (3, 0, −5), and (2, 0, −1)) is reduced relative to the numbers in Model C (although one resonance is more strongly populated (3, −5, 2), which interestingly is the resonance least populated in Model BF₀). Figure 11 indicates that introduction of the SMBH in Model AB₁ weakens the bar, primarily by decreasing the population of many bar-supporting resonances with the smallest pericenter radius. Because the bar is no longer growing after the introduction of the SMBH in Model AB₁ (and the ILR and vILR do not move in radius), there is no active mechanism by which the bar can adapt by capturing stars onto orbits that are stable in the presence of the SMBH. In contrast we saw in Figure 3 that if the bar amplitude continues to grow after buckling, then introducing the SMBH later may have little or no effect on the final bar strength.

Zoom In Zoom Out Reset image size

We note in passing that we also analyzed orbital frequency drift Δ_Ω (Equation (14)) of all 10⁵ orbits in several models at t = 7.57 Gyr. Interestingly we find no difference in the distributions of this measure of chaoticity, and all models have a similar (small) fraction of chaotic orbits. This implies that although stars in Model AB₁ are scattered by the SMBH, by t = 7.57 Gyr they have successfully populated stable regular orbits (primarily at larger distances from the center).

The limited suite of simulations presented (as well as others we briefly discuss in Appendix A) show that if a bar reaches a steady state after buckling, then an SMBH that grows before buckling results in a greater m = 2 amplitude as well as a stronger BP/X shape (characterized by larger values of R_pea and h_pea) than the control simulation without an SMBH or the simulation in which the SMBH grows after the bar reaches a steady state. However, we also saw from our scrambled control models that if the bar continues to grow after buckling, the late introduction of an SMBH may have little or no effect. We now ask if the early versus late growth of an SMBH can be inferred from observations of a galaxy at a single time, for instance from line-of-sight kinematics in external galaxies.

Sellwood & Gerhard (2020) argued that it is possible to distinguish between bars that buckle and cases where a BP/X bulge formed purely via resonant capture (no buckling) by observing differences in the h₄ component of the Gauss–Hermite expansion of the line-of-sight velocity distribution in the two cases. Restricting their measurement to observations of a binned vertical region ∣z∣ < 200 pc, with the disk viewed face-on, they find that as a result of buckling, h₄ becomes strongly negative as previously noted (Debattista et al. 2005). Sellwood & Gerhard (2020) found for "the bar that did not buckle," that h₄ remained positive or became more strongly positive. However, they found that the negative h₄ appeared to trend to 0 with time, and any systematic differences between cases became lost to noise by ∣z∣ ≳ 500 pc.

With the goal of searching for kinematic evidence for early-/late-growing SMBHs, we analyzed 2D kinematic maps (mimicking the Voronoi binned kinematics typically obtained for external galaxies) generated from our N-body snapshots both at early and late times in the evolution. Examples of kinematics maps are presented for four galaxies at the final snapshot (t = 7.57 Gyr) for both the face-on (Figure 13) and edge-on (Figure 14) orientations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Unlike the case in Sellwood & Gerhard (2020), all of our models undergo buckling, although the magnitude of buckling varies and is dependent on when the SMBH is introduced. Therefore, we do not see a clear kinematic difference between the systems in which resonant trapping into the vILR played a stronger role (e.g., Models BF₀, BB₁), and models that buckled more strongly (Models C, AB₁). As previously noted (Debattista et al. 2005; Sellwood & Gerhard 2020), the face-on kinematic maps show strongly negative h₄ values in the central region, changing to positive h₄ values at the ends of the bar. In Figure 13 (face-on) we see that since the bar in Model BF₀ is much longer than in Model C, the distance from the center at which h₄ becomes positive is at or near the edge of the nominal field of view in our kinematic maps, but the behavior is otherwise similar. In Figure 14 (disk edge-on, bar side-on), we see in Model C a positive h₄ in a central ring with two positive h₄ extensions in the disk plane. In Models BF₀ and BB₁, the ring-like h₄ region is broken up and now shows a quadrupolar structure. Although we do not show Models BF₀ and BB₁ at earlier times in their evolution when the BP/X shape was weaker, both of these models showed the same type of complete central ring-like structure with positive h₄, which appeared shortly after buckling (as in Model C), but changed to the quadrupolar structure as the BP/X shape grew stronger. We therefore believe that these differences only reflect the fact that the models with the quadrupolar h₄ have significantly stronger BP/X bulges than the models with ring-like h₄ structures. We also see that the h₃ bimodality along the bar is stronger in Model C than in the models with early-growing SMBHs; but that too appears to be a consequence of the differences in the strength of the BP/X shape. It is outside the scope of this paper to determine exactly which features of orbital kinematics give rise to the different structures in h₃, h₄, and if signatures of evolutionary history (early versus late SMBH growth) may be extracted from a more detailed kinematic analysis.

As discussed in Section 1, it has been known for a long time that the tight M_BH–σ relation observed in elliptical galaxies and classical bulges does not appear to hold in the so-called pseudobulges (with low Sérsic index n or exponential density profiles) or in BP/X-shaped bulges, both of which are thought to have formed by secular evolution. While the central velocity dispersion σ in pseudobulges shows little correlation with M_BH (e.g., Kormendy et al. 2011), σ may be systematically higher in barred spiral galaxies than in spirals with classical bulges (e.g., Graham 2008; Hartmann et al. 2014). The reason for the increase in scatter in pseudobulges and possible offset in barred spirals is still hotly debated.

The simulations presented in this paper have a fixed M_BH for all of the models, but the strength, and therefore mass, of the BP/X structure parameters R_pea and h_pea are correlated with m = 2 bar amplitude. Also, as shown in Figure 6, both the radial and vertical central stellar velocity dispersions (σ_R , σ_z ) are higher in models BF₀ and BB₁ than in Model AB₁, despite the latter having an SMBH of exactly the same mass. Consequently, our simulations have resulted in galaxies with the same M_BH but with BP/X bulges of different masses and different central velocity dispersions. We find that the average 3D velocity dispersion within r = 1.5 kpc in Model AB₁ at late times is 135 km s⁻¹ while in Model BF₀ it is 144 km s⁻¹, which is ∼7% larger for the same M_BH. Furthermore, the average 3D velocity dispersion of Model C at this time is 134 km s⁻¹, nearly identical to Model AB₁, despite not having an SMBH.

This supports previous arguments (e.g., Graham 2008; Brown et al. 2013; Hartmann et al. 2014) that the scatter in the M_BH–σ relation for barred galaxies is due to differences in bar strengths. Furthermore, quantities such as the half-light radius (within which σ is generally computed) are harder to define when dealing with a BP/X-shaped structure rather than a classical bulge or elliptical galaxy since it is not well described by an ellipsoidal light distribution.

The pure N-body simulations studied in this paper showed that the growth of an SMBH early in the life of a disk galaxy dramatically alters the formation and evolution of the bar. The SMBHs in these simulations were grown in an ad hoc manner (similar to previous N-body simulations of SMBH growth in bars), and therefore the simulations did not conserve mass, as one would expect in the case of a real galaxy. Future hydrodynamical simulations including star formation and more realistic prescriptions for SMBH accretion should be carried out to assess how general the results of our work are. In particular, since gas "cools" the disk by producing young stars on nearly circular orbits, disks with gas tend to have lower velocity dispersion. Therefore, simulations with gas and star formation may show less suppression of bar buckling. On the other hand, gas being collisional and dissipative is more likely to experience loss of angular momentum due to shocks in gas arising from their noncircular motions around the bar. Angular momentum transport by the bar could increase the growth rate of the SMBH and somewhat alter the evolution we have observed here. The most important point raised by our experiments is that, contrary to popular belief, SMBHs do not always weaken/destroy bars, and in fact bars have no difficultly forming around fully grown SMBHs or growing co-evally with SMBHs. In general we find that the presence of the SMBH in a growing bar, far from weakening the bar, may actually strengthen it by causing vertical and radial heating that strengthen both the ILR and vILR.