The Stability of Some Galaxy Disks is Still Perplexing

The problem of how some disk galaxies avoid forming bars remains unsolved. Many galaxy models having reasonable properties continue to manifest vigorous instabilities that rapidly form strong bars and no widely accepted idea has yet been advanced to account for how some disk galaxies manage to avoid this instability. It is encouraging that not all galaxies formed in recent cosmological simulations possess bars, but the dynamical explanation for this result is unclear. The unstable mode that creates a bar is understood as a standing wave in a cavity that reflects off the disk center and the corotation radius, with amplification at corotation. Here we use simulations to address one further idea that may inhibit the feedback loop and therefore contribute to stability, which is to make the disk center dynamically hot and/or to taper away mass from the inner disk, which could be masked by a bulge. Unfortunately, we find that neither strategy makes much difference to the global stability of the disk in the models we have tried. While deep density cutouts do indeed prevent feedback through the center, they still reflect incoming waves and thereby provoke a slightly different instability that again leads to a strong bar.


INTRODUCTION
Bars are a common feature of disk galaxies: visual classifications (e.g.Willett et al. 2013) suggest 25%-30% of disk galaxies are barred, consistent with older estimates of the strongly barred fraction (reviewed by Sellwood & Wilkinson 1993).Remarkably, the estimated barred fraction rises to near 70% in galaxy images taken in the near-IR (Buta et al. 2015;Erwin 2018).However, these surveys all agree that a substantial minority of well-observed disk galaxies lack a bar, which presents the problem we address here.For example, Sellwood et al. (2019) were unable to account for the absence of a bar in M33, a particularly well studied galaxy having a gently rising inner rotation curve, a heavy disk, and two prominent spiral arms but no large bar.These properties of M33 are by no means unique; the SPARC sample (Lelli et al. 2016) includes several sellwood@arizona.eduraymond.carlberg@utoronto.caunbarred galaxies such as NGC3877, Holmberg IV, and UGC11557 that also have both heavy disks and gently rising rotation curves.

Disk instabiilities
The bar instability of massive disks has been known for many years (Hohl 1971;Ostriker & Peebles 1973).These early simulations were buttressed by Kalnajs (1978), who presented a global mode analysis of the full-mass isochrone disk, and his predicted mode was confirmed in simulations by Earn & Sellwood (1995).Sellwood & Athanassoula (1986) used the dominant mode of the Kuzmin-Toomre disk as a test case for their mode-fitting procedure, and Toomre (1981) provided examples of the modes of the Gaussian disk.In all these cases, the massive disks of the unperturbed models had gently rising inner rotation curves.Toomre (1981, see also Binney & Tremaine 2008) convincingly accounted for the instability as a cavity mode, with feedback through the center and amplification at corotation, which we describe more fully in §3.2.He also predicted that cutting the feedback loop, e.g. by inserting a dense bulge, which would cause waves to be damped at an inner Lindblad resonance (ILR), could stabilize the entire disk, as Zang (1976) had apparently found for linear instabilities having m ≥ 2 in the Mestel disk.Efstathiou et al. (1982) studied how bar-formation in an exponential disk having constant Q was affected by changes to the rotation curve.They reported that barstability was dependent on the amplitude of the rotation curve, reaching the much-cited conclusion that the halo in bar-stable models should be "the dominant contributor to the total mass."Embedding the disk in a dense halo suppresses the bar instability because m = 2 waves are no longer amplified.However, this is not a satisfactory explanation for the absence of bars in galaxies because the disk should then manifest multi-arm spiral patterns; two-armed spirals, which are the most common patterns in galaxies (Davis et al. 2012;Hart et al. 2016;Yu & Ho 2018), would be suppressed by the dominant halo for the same reason.Efstathiou et al. (1982) also reported that disk stability was independent of the halo core radius, although that second finding was privately challenged by Toomre because it was at variance with his suggestion that a dense bulge should cause most reasonable patterns to be damped at an ILR.Sellwood (1989) indeed found that the small-N simulations by Efstathiou et al. (1982) were affected by non-linear amplified shot noise that overwhelmed the ILR and, in more careful experiments, he confirmed Toomre's contention that a dense bulgelike mass can stabilize the disk.However, the insertion of an ILR is not a panacea for the bar instability, partly because few galaxies can be sufficiently smooth and quiescent over a long period of time for all swingamplified disturbances to be damped at an ILR and partly because Sellwood (2012) later reported that even initially very smooth simulations of the Mestel disk suffered from secular growth of non-axisymmetric stuctures due to impedance changes in the bulk of the disk, caused by weak resonant scattering, which eventually led to a strong bar.
Furthermore, Athanassoula (2002Athanassoula ( , 2008)), Saha & Naab (2013), and Berrier & Sellwood (2016) found that the bar instability is yet more vigorous in simulations that employ a halo composed of mobile particles, rather than a rigid mass distribution.This is because the global mode in the disk is able to elicit a supporting response from the halo that varies in the expected manner with the anisotropy of the halo velocity distribution (Sellwood 2015).
Disks also support other types of mode for which the mechanism is not a standing wave.The most no-table are edge (Toomre 1981;Papaloizou & Lin 1989) and groove (Sellwood & Kahn 1991) modes that are driven from the corotation resonance, but neither creates a bar.Bars can also be formed by modes related to the radial orbit instability (e.g.Lynden-Bell 1979;Polyachenko & Polyachenko 1994) that may operate in globally stable disks, but these ideas do not help to account for the absence of bars in some disks.Bertin et al. (1989) presented a global stability analysis of a large family of disk-halo models, finding that those having cool, low-mass disks and dense bulges supported slowly growing spiral modes.Those authors interpreted the instabilities as cavity modes also, with the more mild being of the type proposed by Mark (1977).These WASER modes invoke travelling waves refracting off a "Q-barrier" from the short-to the long-wave branches of the WKB dispersion relation (see Binney & Tremaine 2008;Sellwood & Masters 2022).Though both trailing waves, the short-and longwaves propagate radially in opposite directions, allowing a feed-back loop that was closed by mild amplification at corotation.As far as we are aware, the only direct test of one of their cases was presented by Sellwood (2011) who reported a long-lived wave having a constant pattern speed in a simulation in which disturbance forces were restricted to m = 2.However, Sellwood (2011) also reported that more vigorous instabilities rapidly emerged in the same model when disturbance forces from higher sectoral harmonics were included.

Cosmological simulations
Cosmological simulations with hydrodynamics are developing apace and create objects that bear some resemblance to galaxies (for a review, see Crain & van de Voort 2023).Algorry et al. (2017) and others have examined the frequency and properties of bars in the simulated galaxies, finding the bar fraction to be somewhat lower than that observed.However, this could simply be due to inadequate resolution, as Zhou et al. (2020) find a higher bar frequency in the Illustris TNG100 models, especially among the higher mass galaxies (Zhao et al. 2020).
Furthermore, Algorry et al. (2017), Roshan et al. (2021), and others have noted that bars in these simulations generally have too low a pattern speed, in the sense that the corotation radius is much larger than the bar semi-major axis, which differs from the properties of observed bars in the nearby universe (Aguerri et al. 2015, see also Buttitta et al. 2023).Once again, Frankel et al. (2022) find the discrepancies of bar properties between simulations and observations are lessened, but not as yet eliminated, as the resolution of the simulations is improved.
Slow bars in the simulations are probably a consequence of dynamical friction from too high a density of dark matter near the centers of galaxies (Weinberg 1985;Debattista & Sellwood 2000, and much subsequent work).Marasco et al. (2020) highlighted halo domination as a shortcoming of the simulations on other grounds.Furthermore, Navarro et al. (2018) studied a multi-arm spiral disk from their simulation, which as already noted, is a symptom of an overly dominant halo; the relatively weak heating and radial migration in the model they studied is likely also a consequence of the mild, multi-arm spirals in that sub-maximal disk.
This very incomplete summary of the relevant literature indicates that the frequency of barred "galaxies" in cosmological simulations is perhaps lower than among galaxies in the local universe, and the bars have properties that differ from those of observed galaxies.However, trends with improving resolution suggest some of these differences may ultimately go away.But the reasons that a substantial fraction are unbarred include (1) stabilization by overly dense halos (e.g.Reddish et al. 2022), ( 2) that infalling sub-halos may destroy bars, though that would both heat and thicken the disk of the host galaxy.As far as we are aware, this is a topic that has yet to be thoroughly addressed (but see e.g.Ghosh et al. 2021).A final possibility (3) is that there is some as yet unknown stabilizing factor at work in both the simulations and in the real universe that prevents some galaxies from forming bars.
But in relation to the question we address in this paper, we do not find that cosmological simulations have yielded any clear understanding of why a galaxy has, or lacks, a bar.Indeed Zhou et al. (2020) concede in their summary that the morphologies of "individual galaxies are subject to the combined effects of environment and internal baryonic physics and are often not predictable."It is important to keep improving the resolution, and perhaps also tweaking the feedback recipe etc., in order to reproduce the observed frequency and properties of bars, but understanding from such complicated simulations the dynamical reason for the eventual match of the models with the observed facts will be extremely challenging.We therefore pursue a parallel investigation using idealized models in which we have some hope of developing deeper insight into this complicated question of disk dynamics.

This paper
We note that the extensive set of models presented by Bertin et al. (1989) explored regimes that have not otherwise been carefully examined.While employing the same type of galaxy model as did Efstathiou et al. (1982), they considered disks that were dynamically hot in the center and/or in which the inner disk surface density had been tapered away, implying a rigid bulgelike component to maintain the adopted rotation curve.They reported that these properties had profound effects on the shapes and growth rates of the dominant linear instabilities.
Our purpose here is to test whether bar-formation can be averted in a moderately heavy stellar disk,1 either by increasing the velocity spread of the stars in the inner disk and/or by reducing the surface density of the inner disk.We therefore present a stability study of models that resemble, but do not exactly match, a small subset of those explored by Bertin et al. (1989).We find fairly vigorous bi-symmetric global modes in most cases with pattern speeds high enough to avoid ILRs.The linear bar-mode can be tightly wrapped in the inner part, while feedback clearly includes a reflection from trailing to leading.Reflection off the center is inhibited in strongly cutout disks, but a new type of instability replaces the classic bar-mode.In both cases, the dominant mode again leads to a large bar.

MASS MODELS
Following Bertin et al. (1989), we consider a family of idealized disk-halo galaxy models in which we determine the dominant instability.In our case, we use 2D quietstart (Sellwood 1983) N -body simulations to follow the evolution of the initial equilibrium model, and fit the dominant mode to the simulation data using the technique described by Sellwood & Athanassoula (1986).

Baseline model
Our baseline model is of the type originally proposed by Fall & Efstathiou (1980).The surface density of a flat, axisymmetric exponential disk has the radial profile where R d is the disk scale length and M d is the nominal mass of the infinite disk.We limit its radial extent using a cubic function to taper the surface density from Σ(5R d ) to zero at R = 6R d .The rotation curve is that of a cored isothermal sphere with r c being the core radius, although we will be interested exclusively in quantities in the disk plane where r = R.The implied halo density is whatever is required, when combined with the disk attraction, to achieve this rotation curve in the disk plane.We relate the rotation curve to the disk properties by setting V 0 = 0.9(GM d /R d ) 1/2 , and generally choose r c = 0.5R d .
As the maximum circular speed arising from a razor-thin exponential disk is 0.622(GM d /R d ) 1/2 , the disk, though quite heavy, has less than the required mass to account for the central attraction at any radius.
Since we are here interested in bisymmetric linear instabilities of the disk only, we generally do not compute the axisymmetric central attraction of the disk.We compute only the m = 2 disturbance forces from the disk particles, and then add the central acceleration −V 2 (r)/r to every particle at every step.Note that this strategy implies that we represent the halo component as a rigid mass distribution.
We set the radial velocity dispersion of the disk particles using the Toomre (1964) (3) and κ is the local epicyclic frequency (Binney & Tremaine 2008).Note that we generally use a radially dependent Q function, although in all our models we set Q OD = 1.2 in the outer disk.
We adopt the distribution function for the disk component using the form proposed by Shu (1969) ) Here E is the excess energy of a particle above E c , which is that of a circular orbit at the guiding center radius R g (L z ).Although this DF assumes no retrograde stars, we later reverse the angular momentum of some low-L z particles in order to smooth the discontinuity in f (E, L z ) at L z = 0, which does not affect the equilibrium.The function (4) clearly assumes a Gaussian velocity distribution at all radii, which Shu (1969) argues is the appropriate form for a (partially) relaxed disk.We select particles from this DF using the method described in the appendix of Debattista & Sellwood (2000).
The function F (L z ) has to be determined numerically and the procedure we adopt is described in the on-line manual (Sellwood 2014).As there are many possible functions F that fit the adopted disk surface density, we impose two extra requirements.Not only are rapid fluctuations of F with L z physically unreasonable, but we have also found that even mild "ripples" in the function F (L z ) can seed additional disk instabilities related to groove modes (Sellwood & Kahn 1991).We therefore penalize the fit to the surface density also to minimize (5) Note that the numerical search for the optimum F seeks a balance between fitting the disk surface density while also minimizing T , and finding the optimum balance is something of an art.It is also required that F (L z ) ≥ 0 for all L z , although we find that this requirement is generally satisfied for a smooth F without imposing an additional constraint.
In order to create a quiet start, we place three copies of each particle almost regularly around a half-circle.By restricting disturbance forces to m = 2 only, these three particles mimic an initially smooth, circular wire of uniform mass per unit length that osillates radially at the epicyclic frequency and distorts in response to particle dynamics in the global gravitational potential.
Here, and throughout the paper, we use units such that . For those who prefer physical units, a possible scaling is to set R d = 2 kpc, and τ dyn = 10 Myr, which implies V 0 ≃ 176 km s −1 and M d ≃ 1.78 × 10 10 M ⊙ .

Variants of the baseline model
Following Bertin et al. (1989), we consider models with prescribed Q-profiles and disks having central mass cutouts, and in two cases we also increase the halo core radius r c .
The Q-profile has the functional form where R Q is the radial extent of the higher Q values, Q jump is the difference between Q(0) and Q OD , the argument x = 1 − 2R/R Q varies from 1 ≥ x ≥ −1 and the Fermi-like taper function rises smoothly from zero at x = −1, to unity at x = +1.The central surface density is cut away in the following manner where the central depth D = Σ c (0)/Σ 0 , and in this case x = 2R/R cut − 1.Since the rotation curve is unaffected, deeper cutouts imply the galaxy model has a greater bulge mass.Note that σ R will also be reduced as more mass is cut away (see eq. 3), provided that the prescribed function Q(R) is unchanged.The functions eqs.( 6) and (8), illustrated in Figure 1, introduce four new parameters: R Q and R cut which are the radial extents of the enhanced Q values and of the surface density cutout, and Q jump and D. In principle, R Q and R cut could differ, but we have found it convenient to keep R Q = R cut , which reduces the number of extra parameters by one.

Simplifying assumptions
Our simulations embody several approximations to reduce the computational cost.As this is a preliminary search for stabilizing influences, any models that may turn out to be promising would need to be resimulated without these approximations.Specifically: • Particles are restricted to motion in a plane as a first approximation.Note large-scale nonaxisymmetric modes are only slightly weaker in a moderately thickened disk than in a 2D disk, and gravity softening provides some allowance for disk thickness.
• The restriction of disturbance forces to a single sectoral harmonic does not matter, because each sectoral hamonic behaves independently until the amplitude becomes large, and here we wish to measure the mode frequency in the linear regime.
We confirm in §4.2 that including extra harmonics does not alter the linear mode and causes differences at large amplitude only.
• We employ a rigid halo since we are searching for stabilizing influences.If a disk is unstable with a rigid halo, then it has no chance of being stable in a live one, but any apparently stable cases could be re-simulated later using a live halo.
• We have also neglected the role of gas.Recall Sellwood et al. (2019) reported that the global stability of the disk of M33 was very little affected when the observed gas component was treated either as collisionless particles or as gas in a number of different ways.Note that gas comprises 30% of the total disk mass in that galaxy, though the atomic component is more spread out than are the stars.
We critically re-examine our adopted model and simplifiying assumptions in §6.
is available for download.Table 1 gives the values of the numerical parameters adopted for most simulations presented in this paper.Fourier analysis of the mass distribution on each grid ring, which separates the solution for the field into different sectoral hamonics, makes it easy to restrict the disturbance forces to those arising from m = 2 distortions of the particle distribution.We report checks in which we vary these parameters in §5 below, and a test with a 2D Cartesian grid in §3.4.
As usual, we measure non-axisymmetric distortions of the distribution of the N particles using an expansion in logarithmic spirals: where (R j , φ j ) are the polar coordinates of the jth particle at time t, m is the sectoral harmonic, and γ is the (radially constant) angle of the spiral component to the radius vector, which is the complement to the spiral pitch angle.

Mode fitting
Recall that a normal mode is a self-sustaining, sinusoidal disturbance of fixed frequency and constant shape.The perturbed surface density of a mode in a galaxy disk is the real part of where the frequency, ω, is complex when the mode grows or decays.The complex function A m (R), which is independent of time, describes the radial variation of amplitude and phase of the mode.Here, ω = mΩ p + iβ, with Ω p being the pattern speed and β the growth rate.
The results in Table 2 were obtained by fitting the function (eq.10) to the data from the simulation using the least-squares procedure described by Sellwood & Athanassoula (1986).The input data to the fit are either the logarithmic spiral transforms (eq.9) or the amplitude and phase of the m = 2 component of the mass assigned to the points on each grid ring.A quiet start reduces the seed amplitude of the mode to well below the level that would be expected from particle shot noise, and is essential to allow a long enough period of linear growth to obtain a credible measurement of the growth rate.In a few cases, a fit of two superposed modes seems to be preferred, but we list only that with the highest growth rate.
We quote uncertainties in the measured frequencies in Table 2 that span the entire spread of values from fits to the data of both types and over slightly different time ranges in a single simulation.However, they are probably severe underestimates because different real-Figure 2. The black lines in the upper panel present the amplitude of the logarithmic spiral transform of the particle positions (eq.9) at intervals of 10τ dyn in the baseline model.The red, dotted curves show the fit of a single gowing mode to the data at the same times.The contours in the lower panel are of the best fit mode taken from a fit over the same time interval as above, but to the density transforms on the grid rings.The full-drawn and dotted circles mark respectively the radii of the CR and the OLR -the pattern speed of this mode is high enough to be well clear of possible ILRs.
izations of the same model turn out to have unstable mode frequencies that differ by some ∼ 10%.Choosing different parameter values, such as the weight given to the smoothing term T in eq. ( 5), results in a different realization that is not exactly the same model.Even though the fitted surface density and Q-profiles of two or more realizations are barely distinguishable, we found that minor local differences in the density of particles as a function of L z affect the frequency of the fundamental mode to a surprising extent.

WKB waves
We interpret the properties of the fitted modes in terms of a local dispersion relation for density waves in disks, which connects the frequency of a steady wave ω to its radial wavenumber |k|.The formula is independent of the sign of k, which is conventionally taken as positive for trailing waves and negative for leading.It embodies the WKB approximation that computes the self-gravity of the spiral as that of a plane wave in a thin sheet (Binney & Tremaine 2008).The version given by Lin & Shu (1966) which states that the self-gravity term decreases the wave frequency ω = m(Ω p − Ω), here assumed to be purely real, below the natural frequency of radial oscillation κ.The "reduction factor" F ≤ 1 (given by Binney & Tremaine 2008, their Appendix K) depends upon Q, k, and ω, and quantifies the extent to which the self-gravity term is weakened by random motion.Eq. ( 11) has severe limitations (Sellwood & Masters 2022).It applies to waves that are sufficiently tightlywrapped that the spiral pitch angle can be neglected, allowing the radial frequency ω, to be related the Dopplershifted frequency at which stars encounter an m-fold symmetric spiral ω = m(Ω p − Ω) (Lin & Shu 1966).It also applies equally to leading and trailing waves because it also omits any hint of swing-amplification (Toomre 1981) near corotation.Despite these limitations, the short wave branch (Binney & Tremaine 2008) does give some qualitative indication of the behaviour spiral waves away from corotation.Toomre (1969) pointed out that a spiral wave packet propagates radially across the disk at a group velocity v g = ∂ω/∂k, which may be computed from eq. ( 11).He found that short waves travel away from or toward the corotation resonance (CR) when they are respectively trailing or leading and that, other things being equal, v g should be decreased somewhat when Q is raised.His numerical calculations showed that eq. ( 11) breaks down over a broad region around the CR, but it makes reason-ably accurate predictions elsewhere.The wave carries angular momentum at the group velocity.

THE MODE OF OUR BASELINE MODEL
Our baseline model is globally unstable, consistent with the earlier study by Efstathiou et al. (1982).The disk is massive enough and the rotation curve rises slowly enough that a global cavity mode of the type described by Toomre (1981) must be expected.

Linear mode
Figure 2 illustrates the fitting procedure for the simulation of our baseline model.The solid black curves in the upper panel show the amplitude of the logarithmic spiral transforms (eq.9) of the particle positions at intervals of 10τ dyn during the period of linear growth i.e. excluding early times that were noise-dominated, and later times when the mode saturates.The logarithmic amplitude scale reveals approximate equal spacing of these curves over time, indicating steady exponential growth.The red, dotted curves mark the amplitude of the fitted function (eq.10) at the same times, and are in reasonable agreement with the data.Note that the data used in the fit are from transforms taken five times more frequently than those illustrated.A bias towards higher amplitude of trailing waves (tan γ > 0) over the leading components (tan γ < 0) is evident, but the mode has significant amplitude over the range −5 tan γ 5.
We also fitted a mode to the density transforms on the grid rings over the same time interval, and draw the fitted function A m (R) in the lower panel; note that the contours are of positive overdensity only, the underdense part is not contoured.We have drawn a solid circle to mark the radius of the corotation resonance (CR), and a dotted circle at the radius of the outer Lindblad resonance (OLR); the pattern speed of this, and all other modes in this paper (Table 2), is high enough that there are no ILRs.
The mode has a tightly wrapped spiral appearance in the inner parts (lower panel of Fig. 2), that differs from the open spiral of the domininant mode of fully self-gravitating disks that has been reported by others.We will show that the different apppearance of the linear mode is largely due to the adopted rotation curve, which raises the orbital frequencies in the inner disk above those that arise in most self-gravitating disk models.The mode in Fig. 2 is more open in the outer parts, however.

Bar mode mechanism
As noted in the introduction, Toomre (1981) elucidated the mechanism of a bar mode as that of a cavity mode, or standing wave, between the disk center and corotation.More specifically, a trailing wavetrain propagates inward at the group velocity and reflects off the center into an outwardly propagating leading wavetrain, but the second reflection off the CR causes the wavetrain to be strongly swing-amplified into an amplified trailing wave, providing positive feedback that leads to an unstable run away -but only until the failure of the smallamplitude approximation that underlies the feedback loop.The growing wave in the cavity between CR and the center removes angular momentum from the inner disk (Lynden-Bell & Kalnajs 1972;Binney & Tremaine 2008;Sellwood & Masters 2022), and since the net angular momentum of an isolated disk cannot change, the outer disk must accept that removed from the inner disk by the disturbance.An exclusively trailing wave carrying positive angular momentum propagates outward where it wraps more tightly as it is finally absorbed by wave-particle interactions at the OLR (Lynden-Bell & Kalnajs 1972).
It seems plausible that this mechanism operates in the mode reported in Fig. 2, since the disturbance reaches the very center of the disk, where it reflects.The bias of the leading components of the transform in the upper panel being weaker than the trailing components, as expected from swing amplification at the CR, gives rise to the overall trailing appearance of the disturbance.But the fundamental modes of disks having gently rising rotation curves do not manifest the multiple sub-peaks along the mode ridge line that we find here.These subpeaks are regularly spaced 90 • apart, suggesting interference between the leading and trailing waves, similar to that Toomre (1981) reported for the overtone modes of the Gaussian disk.The mode transform in the upper panel also has more structure than is typically seen for open bar modes, doubtless as a consequence of interference between the leading and trailing waves.
The frequency of the mode is determined by the phase closure constraint, which means that the complete feedback cycle must encompass an integral number of wave periods, else the travelling waves will selfinterfere.Normally, in simple full-mass disk models such as the isochrone, the dominant mode is the fundamental with a single antinode in the cavity.Overtones having two, three, and more antinodes were reported by Toomre (1981) for the Gaussian disk, but they all had lower growth rates than the fundamental.In our case, we see several antinodes already in this, the dominant mode of our baseline model, presumably because the time re- quired for the wavetrain to travel from CR to the center and back is far longer than the oscillation period of the wave.But it unclear to us why this particular mode, of all overtones that are possible in this model, should have the highest growth rate, and therefore stand out in our simulation.
The mechanism proposed by Toomre (1981) suggests that the mode growth rate is determined by the amplification factor at the CR divided by the time required for the wavetrain to travel to the inner reflection radius and back to the CR.However, it would be hard to predict what this should be because the group velocity of the wavetrain (Toomre 1969;Binney & Tremaine 2008;Sellwood & Masters 2022) cannot even be estimated from WKB theory for part of the required radial range because eq. ( 11) has no solutions in the "forbidden region" when Q > 1.While it is known that growing waves do propagate in this region, we do not have a dispersion relation for them from which we could estimate their v g .It is clear that the model forms a strong bar by t = 500, but it is interesting that the first feature to become visible at t = 400 is the small inner spiral which later settles into a short inner bar that rotates with the large bar.

A further numerical check
It seemed possible that the choice of overtone referred to above, and other aspects of the mode, could have been preferred by our polar grid.We therefore reran the same model using a 2D Cartesian grid, with two-fold rotational symmetry imposed to suppress odd sectoral harmonics and a strategy to kill forces that may arise from m = 4 density variations, as described by Sellwood (2020).In order to maintain the same central attraction as on the polar grid, we assigned the mass of a smooth, tapered exponential disk to a separate copy of the grid and then subtracted this unchanging smooth disk from the masses of the moving particles that were assigned to the working grid at each step, before solving for the gravitational attraction of the residual density variations (Sellwood 2020).Thus the accelerations acting on the particles are those arising from the disturbance density distribution together with the fixed attraction implied by the rotation curve, as in our baseline simulation.The six-fold rotational symmetry of the quiet start, does not affect the dynamics on the polar grid, where forces were restricted to m = 2, but is a nuisance on the Cartesian grid.Therefore, from the same basic file of particles, we placed 20 instead of three on each half-ring to make the unperturbed mass distribution 40-fold symmetric, which suppressed all likely non-axisymmetric forces except for those arising from bisymmetric disturbances.
We employed a square grid with 1024 cells on a side, set the disk scalelength R d = 80 mesh spaces and the Plummer gravity softening length ǫ = R d /20 = 4 mesh spaces to be consistent with disturbance forces on the polar grid.Our numerical concern was allayed when we found that the dominant disturbance in this simulation had the same pattern speed, to within 0.2%, growth rate to within 1%, and mode transform that was indistinguishable from that in our baseline run on the polar grid (Fig. 2, top panel).

Reducing the inner slope of the rotation curve
We can understand the tight winding of the inner part of the mode in Fig 2 through the WKB dispersion rela-tion.Rearranging eq. ( 11), we find Note that this is not an explicit expression for k, because the factor F in the denominator also depends on k, but it strongly suggests that one way to prefer tightly-wrapped waves, or large |k|, is to increase κ.Thus in the inner disk of our model, the numerator in eq. ( 12) is increased as r c is decreased, which is the main reason that the mode in Fig. 2 is more tightly wrapped than modes in disks having more gently rising rotation curves.
In order to confirm that the tight winding of the inner part of the mode was due to the smallish core radius of the global potential, we present in Fig. 4 the fitted linear modes of models Rc0.75 (left) and Rc1.0 (right), in which we increased the core radius from the baseline value (r c = 0.5) to respectively r c = 0.75 and r c = 1.0, but we preserved a flat Q = 1.2 at all radii and did not include a cutout.The inner slope of the rotation curve in model Rc1.0 is very slightly below that expected from the disk mass alone, and so implies a hollow halo at radii r 0.3, but as the axisymmetric self-gravity of the disk is ignored, the particles simply move in the global  0) is increased; in all cases, RQ = 1.5 and there is no cutout.There was a hint of a second mode in the data for the LH panel but the fit presented is for a single mode.We have fitted two modes to the data in the upper middle panel, but present only the dominant mode in the lower middle panel.The increasing size of the mode shapes in the lower panels reflects their decreasing pattern speeds, which are given in Table 2.The growth rates of the modes in HC1.0 and HC1.5 are about the same as that in the baseline model, but the mode HC2.0 (right) grows less rapidly.
potential and the disk instability is unaffected by this undesirable property.
Notice from the top right hand panel of Fig. 4 that the logarithmic spiral transform of the mode differs qualitatively from those from models with smaller r c , although the mode still has a substantial leading component.The mode shape in the lower right panel is more open, having just a single extra antinode near the center, confirming our suggestion that the steeper inner rise to the rotation curve is indeed the cause of the tightly-wrapped inner spiral.Again, the coherent disturbance continues into the very centers of these two simulations, indicating that these modes also reflect off the disk center.

OTHER RESULTS
Table 2 lists the parameter values for most of the simulations we report here.It also gives the estimated eigenfrequency, ω, of the dominant linear mode.In all simulations other than the two just presented, we kept the core radius of the potential, r c = 0.5, as in our baseline model.
The purpose of these simulations is to explore models having higher Q and/or surface density cutouts in the inner parts of the disk.It has proved informative to compare the frequencies, mode transforms, and mode shapes in sequences of runs from Table 2 in which all but one parameter were held fixed.

Disks with hot centers
We first present a sequence in which we varied the central value of Q, keeping R Q = 1.5 and with no disk cutout.We chose Q(0) = 1.0 (HC1.0)(lower than the baseline model), Q(0) = 1.5 (HC1.5) and Q(0) = 2.0 (HC2.0).The linear modes of these three cases are presented in Fig. 5, and the baseline model fits into this sequence with Q(0) = 1.2.
The principal consequence of increasing Q(0) is to slow the pattern speed, which moves the resonances to large radii and increases the spatial scale of the mode (lower panels).Somewhat surprisingly, the growth rate was little affected as Q(0) increased from 1 to 1.5, but it decreased when Q(0) = 2 (Table 2).Furthermore, the shapes of the mode transforms in the upper panels all resemble that of our baseline model, with strong leading components and, despite the higher Q in the inner Figure 6.Fits to the linear evolution of three simulations having increasingly deep central cutouts.In all cases Rcut = 1.5, while Q(R) = 1.2 everywhere, and the baseline model (Fig. 2) having no cutout is the first of this sequence.
disk, the modes all appear to reflect off the center (lower panels).

Cutting away disk mass
Fig. 6 presents results from a sequence of simulations in which we tapered away increasing quantities of the disk surface density.We varied D = 0.8 to D = 0.2, keeping R cut = 1.5 and maintaining a flat Q = 1.2 at all radii.Again, the baseline model begins this sequence with D = 1.0.
We have very successfully fitted two modes to the logarithmic spiral transforms in the upper left panel, but draw only the dominant mode in the lower left.The character of the mode transform (upper panels) differs in the middle and right from that in the left panel, and this is again reflected most clearly in the mode shape in the bottom right, which manifestly does not extend into the center -we have vainly searched for any coherent waves in the region R < 0.5 in 1DC0.2.The inner edge of the mode in 1DC0.5 (middle panel) is not as decisively outside the center, but we could find little evidence for a coherent wave inside R < 0.4.Also we do not find a clear trend in the pattern speeds as the cutout is deepened, but the vertical spacing of the curves in the upper panel gives a clear visual indication that the growth rates of the dominant mode rises from left to right, as confirmed by the numerical values in Table 2.
One consequence of cutting away mass from the inner disk is that the denominator of eq. ( 12) contains the factor Σ, suggesting that k should increase and spiral modes would become still more tightly wrapped.We do not see this happen in Fig. 6, but instead the disk center, the mode would be most tightly wrapped, does not appear to support a disturbance at all.We address this result in §5.
Fig. 8 displays the last part of the evolution of a rerun of simulation 2DC0.2 to show how a deep inner cutout affects the non-linear formation of a bar.In addition to the usual m = 2 disturbance forces, this simulation included the self-consistent axisymmetric term (m = 0) of the disk particles as well as the m = 4 term, in order to capture their influence on the non-linear evolution.Note, the inclusion of the m = 4 term required us to increase the number of particles on each half-ring from Figure 7.A parallel sequence to that in Fig. 6 but with the central value of Q also raised to Q(0) = 2.Note Rcut = RQ = 1.5 in all three cases.The first member of this sequence, is HC2.0 in the right hand panel of Fig. 5 and the final member of the sequence, 2DC0.2, is presented in Fig. 11.
3 to 5 in order to maintain a quiet start.The extra force terms have no effect on the m = 2 mode frequency or shape, as expected from theory, and we were unable to detect an instability at m = 4.The limiting amplitude of the bar is increased over that in the m = 2 only case and the resulting bar, which has a density minimum at its center, is a little thinner.In reality, the density minimum would probably be obsured by a dense bulge.

Varying the radial extent of the cutout
Fig. 9 shows the consequence of varying the radial extent of the cutout.In all three cases, Q(0) = 3 and D = 0.2, but R cut = 1.0 in RQ1.0 (left), R cut = 1.5 in RQ1.5 (middle) and R cut = 2.0 in RQ2.0 (right).
We found from the modes in Fig 6 that a deep central cutout out erases the multiple peaks in the mode transform for the baseline model (top panel of Fig 2) and similar cases, as it appears to have done here, but here we have also raised the central Q(0).Despite the simple appearance of the mode transform, the modes in Fig. 9 are still quite tightly wrapped, though only in the cool (Q = 1.2) outer disk beyond R Q = R cut .The growth rate of the dominant mode is little affected by increasing R cut , but the pattern speed is reduced and, while the spatial scale of the mode rises as a consequence, it appears to do so roughly homologously.Note that none of these modes extends into the center, and we find that the inner edge of the mode moves to larger radius as the radial extent of the cutout is increased.

MODES IN CUTOUT DISKS
We reported in §3 that the mode in our baseline model reflected off the center.Fig. 10 presents additional evidence in support this statement because the mode amplitude indeed remains significant to small radii and declines steeply to R = 0.
However, in some of the cases reported in §4, it seemed that the mode amplitude dropped to near zero inside a finite radius, and we here show another particularly clear example in Fig. 11.This model, 2DC0.2, continues the sequence shown in Fig. 7 to the case for which D = 0.2, and has no coherent wave at radii R 0.3.We have added the bottom panel confirming that the amplitude inside this radius is ≪ 1% of the peak and seems consistent with noise, which implies that the disturbance is confined to radii R 0.3.Note that the mode still possesses a strong leading component (top panel of Fig. 11).

Numerical checks
We have checked carefully that the inner edge of the mode it is not a numerical artifact.It seemed possible that representation of a very tightly wrapped trailing wave on a polar grid may alias as a similarly tightly wrapped leading wave, causing the disturbance to reflect off the grid.In order to check that this is not the case, we reran the same initial file of particles in two additional simulations in which we first doubled and then halved the number of grid points in each dimension.These changes will have, respectively, halved and doubled the spacing of the grid rings, but neither the inner edge of the mode, as shown in the left panel of Fig. 12, nor the mode frequency changed thereby ruling out reflection due to grid aliasing.
Another possibility is that gravity softening may be playing a role, since it increasingly attenuates the selfgravity of density waves by the factor e −|k|ǫ , and may prevent them from propagating at all when the wavelength (= 2π/|k|) decreases below some possible minimum value.We have therefore tried both halving and increasing the softening length by factors of two and four.Note that a change to the softening length affects the entire mode, changing both parts of the eigenfrequency, the positions of the resonances, and the radial wavelength of the mode.However, the mode in each case again had tiny amplitude at radii R < 0.2, though in two cases it was small but greater than zero between 0.2 < R < 0.3.We have also verified that halving the time step made no difference to the mode shape or frequency.These successful checks indicate that the mode not extending into the center is physically real.

Summary of the evidence
We report in Fig. 13 the radial amplitudes of fitted eigenmodes as the physical parameters D and Q(0) are varied.The sequence in the top panel indicates that the modes do not reach the center when D ≥ 0.5, a finding that seems to hold in all simulations reported in Table 2, irrespective of either R cut or Q(0).
In fact, the evidence in the lower panel of Fig. 13, which shows a sequence with increasing Q(0), indicates that the inner radius of the mode is little affected by increasing random motion.It is true that the inner radius of mode does move inwards (slightly) as Q rises -the cyan line is from the model with the hottest center, the red from the coolest.However, the important density gradient is probably that of the guiding centers, which is approximately unchanged in this sequence of models, and the sharp inner edge of the mode in the cooler disks is increasingly blurred by larger epicycles in the hotter disks.Note also that we reported in the right panel of  Fig. 12 that the inner radius of the mode is also little affected by the value of ǫ.

Mode mechanism
Since it was provoked by a steep inner density gradient, it seemed natural to suppose that the inner mode was related to the better-known outer edge modes of disks (Toomre 1981;Papaloizou & Lin 1989).Outer edge modes occur where mild non-axisymmetric dis-tortions in a steep surface density gradient create coorbiting over-densities.Each over-density induces a swing amplified (Julian & Toomre 1966;Binney 2020) supporting wake in the interior disk that gives angular momentum to the edge overdensity, causing it to rise outwards, and therefore the edge distortion grows.
Thus, a possible inner edge mode could operate by a similar mechanism.Co-orbiting overdensities could be created by mild non-axisymmetric distortions at the radius of steepest density gradient that would induce a strong trailing response in the higher surface density region.The attraction from this wake would cause the overdensities to be pulled back and therefore to sink towards the rotation center, creating a similar instability.However, it cannot be the mechanism of the modes we report from our simulations, because then the corotation radius of the modes would be close to the radius of the steepest gradient, which it manifestly is not in Figs. 6,7,9,and 11.Similar inner modes were reported by Zang (1976, see also Evans & Read 1998;Sellwood & Evans 2001), who found instabilities in the Mestel disk when the center was cut out too sharply.In his case also, the slow pattern speed of the m = 2 mode placed corotation at more than twice the radius of the steepest surface density gradient, as is the case for the modes we find.Zang noted in his thesis "Although we have some ideas about the details of this process, we are presently unable to offer a clear (sic) physical mechanism" for the instability provoked by the sharp inner cutout.Evans & Read (1998) remark "If,  2. All curves were rescaled to span 0 ≤ |A2(R)| ≤ 1. Left panel: the red and blue lines are from cases where the resolution was respectively doubled and halved from our standard grid, which is reproduced as the green line.Right panel: the red, blue and cyan lines are from cases in which ǫ was respectively halved, and increased by factors of two and four from our standard value, which is reproduced as the green line.
however, the cut-out is sufficiently sharp (as for larger cut-out indices), it presents a barrier which reflects the incoming trailing waves."Thus they clearly view the mode as a slightly different type of cavity mode, but they do not offer any further explanation for the nature of the barrier or the reflection.
The strength of the leading wave component in every one of the modes in our simulations is clear evidence that swing amplification plays a major role, even when the mode does not reach the center.As explained above, an outer edge mode does not require a feedback loop, because the instability is excited at corotation, but it is hard to see how the slow modes, found by Zang (1976), Evans & Read (1998), and here in this paper, could grow without feedback.An inner reflection off the center is excluded for m ≥ 2 waves in cusped potentials, and we observe that the modes we find in cutout disks do not reach the center either.So we concur with Evans & Read (1998) that an inner reflection must occur at a finite radius in order to turn inwardly propagating trailing waves into outwardly propagating leading waves.

Possible reflection mechanism
We offer the following speculative idea for how an inner reflection could happen.WKB wave theory (eq.12) suggests that spiral waves should become more tightly wrapped in models having deeper surface density cutouts.Thus inwardly travelling waves will wrap ever more tightly as they enter the cutout region.As this happens, the ability of a disk with random motion to carry the wave must dwindle because stars whose epicycle radii are greater than the diminishing wavelength of the wave cannot provide much of a supporting response.A wave supported by an ever decreasing fraction of stars in the center of the velocity distribution cannot sustain its amplitude as |k| rises.The collisionless nature of a stellar disk excludes possible dissipation and, if there is no nearby resonance where wave-particle interactions can occur (Lynden-Bell & Kalnajs 1972), wave action must be conserved and the only possible outcome is that the wave must bounce.
To follow up this idea, we have attempted to use the measured pattern speeds from Table 2 to solve eq. ( 11) for k(R), including the full functional form for F as well as a factor e −|k|ǫ that further attenuates the gravity term.Unfortunately, setting Q(R) from eq. ( 6), together with the gravity softening term, prevented us from finding solutions for |k| over almost the entire radial range inside the CR, because the forbidden region, an artifact of eq. ( 11) that stems from assuming a steady wave, is so broad.We were therefore unable to make even a very rough estimate of radius at which the wave would bounce, which anyway would be inconsistent with the infinite plane wave approximation that underlies eq. ( 11).
Not only is the idea that the wave might reflect off the decreasing density gradient highly speculative, but it is also hard to reconcile with the evidence in the bottom panel of Fig. 13, which indicates that the reflection radius is little changed as the central Q is increased from 1.2 to 3, with everything else held fixed; the slight change could be consistent with epicyclic blurring, as we note above.It seems unlikely that the reflection radius would be independent of the degree of random motion, since |k| for steady waves is strongly affected by Q, but we have no other idea to offer.

Gentler inner tapers
Since the instability of disks in cusped potentials disappears when the taper is made more gentle, it is natural to ask whether the same happens in our case.Our inner taper function (eq.8) differs from that used by Zang (1976) and Evans & Read (1998), but we can make it more gentle by increasing R cut .We therefore present two more simulations, listed as the last two entries in Table 2.These two models differ in just one respect from 2DC0.4, reported in the right hand panels of Fig. 7: R cut = 2.5 in model MR2.5 and R cut = 3.5 in model MR3.5.
The upper panel of Fig. 14 illustrates the slower and milder instability, in comparison with that in 2DC0.4,we obtained from simulation MR2.5.Notice from the lower panel that the reflection radius is a little closer to the center than that in model 2DC0.4 even though R cut has been increased.We could not detect a growing disturbance in MR3.5.Thus we do find that gentler inner tapers are stabilizing, as did Zang (1976) and Evans & Read (1998), but we believe it is for a different reason, as follows.
The m = 2 instability of our baseline model is a standing wave between the center and corotation where it amplifies, as in the usual bar-mode picture.When we apply a deep cut out, we find that reflection can occur away from the center, but we still have a cavity mode with amplification at corotation.Zang (1976) and Evans & Read (1998), on the other hand, calculated disk modes in a cusped potential, which precludes feedback through the center (except for m = 1), finding that models having gentle cutouts had no global instabilities for m ≥ 2. They also found that sharp inner cutouts provoked instabilities, and Evans & Read (1998) asserted that swing amplification farther out in the disk could combine with reflection off a sharp cutout to allow a cavity mode.So the modes in sharply cutout models, in both their case and in ours, apparently have the same mechanism.
However, when we make a gentler inner cutout by increasing R cut at a fixed central depth, we also cut away a lot more mass from the disk, which eventually turns off m = 2 swing amplification because X > 3 over the whole disk, as shown in Fig. 15.The different taper function in their otherwise smooth, self-similar disks leaves the swing amplifier still eager to operate in the outer disk, but it is denied feedback because the ILR absorbs infinitesimal incoming waves when the inner taper does not make a barrier.Thus the explanation for stability in cusped models is quite different: the stability of MR3.5 results from no other physical mechanism than that argued by Efstathiou et al. (1982), who found that increased halo crushes m = 2 swing amplified waves!

DISCUSSION
The simulations presented in this paper capture the important gravitational dynamics of the disk, but are otherwise highly idealized, with several simplifying assumptions summarized in §2.3, whereas the real universe is far more complicated.In the light of this consideration, and at the suggestion of the referee, we here re-examine our assumptions and findings from isolated simulations.

The case of M33
The local group galaxy, M33, is perhaps the poster child of the bar stability problem.Sellwood et al. (2019) closely examined the stability of M33 in a set of models that matched all available data, but were unable to account for the absence of a bar in that well-studied galaxy.Their fully self-consistent 3D simulations with live halos modeled the gas in a variety of ways, none of which prevented the rapid formation of a bar.The stellar and gaseous density profiles of the disk were constrained by observations, which exclude possible cut outs and, although the stellar velocity dispersion is known only in the center, that is the most important constraint.The regular gas kinematics of the inner disk (Corbeli et al. 2014;Kam et al. 2017) indicates that it has not been disturbed by infalling substructure in the recent past and the pronounced warp is sufficiently far out to have no plausible effect on the stability of the inner disk.While all their models had smooth halos, some rotated both with, or counter to, the rotation sense of the disk.Reducing the stellar massto-light ratio to some quite unreasonably low value did lead to milder instabilities, but the disk then supported multi-arm spirals, whereas near IR images of the galaxy (e.g.Kam et al. 2015) indicate two major spiral arms with perhaps a weaker third arm.Their attempts to induce a open bi-symmetric spiral in a lower mass disk through a tidal interaction were not successful.Smercina et al. (2023) posted a study of star counts in M33 which indicated a bar-like structure in the inner 1kpc of the disk.We note that the strong bar created by global instability in the simulations by Sellwood et al. (2019) had a semi-major axis of ∼ 3kpc, and therefore the finding of a short, weak bar in M33 does not alter the conclusion that the apparent global stability of M33 presents a dynamical challenge.
Thus the study by Sellwood et al. (2019) did not make many of the simplifying assumptions listed in §2.3, yet all their plausible models still formed a strong bar within 1 Gyr.A possible weakness of their study is that the halos lacked the substructure expected in a hierarchical universe.However, massive subhalos would disrupt the disk, which clearly has not happened recently in M33, while low-mass subhalos would have little effect on disk stability and anyway probably would not survive at all in the inner halo (e.g.Sawala et al. 2017).

The present study
Our purpose in this present study was to revert to highly idealized simulations in which the dynamics is firmly under the control of the experimenter, to determine whether the properties of the disk center could be changed in such a way as to inhibit the bar instability.Although we made a number of simplifications, listed in §2.3, in order to obtain reliable results at low computional cost, those approximations have been examined before and most are known to have little or no effect on the stability of the entire disk; the exception being the assumption of a rigid halo that, when relaxed, is known to enhance the growth rate of the bar mode, as noted in the introduction.
Unfortunately, we have found that neither deep cut outs nor a high central Q have much effect on the global stability of the models we have tried.These models can be thought of as more closely corresponding to galaxies, such as the Milky Way, that have higher mass than does M33, though we made no attempt to model any particular galaxy, neither do they bear a close resemblance to those created in cosmological simulations.It is therefore reasonable to ask whether more realistic models would yield a different result?Our baseline model adopts an exponential disk and a cored isothermal potential.
It seems unlikely to us that a different disk mass profile would make a qualitative difference to global stability.We base this expectation on Toomre's mechanism for the bar mode, which requires only swing amplification at corotation and an inner reflection of inwardly traveling trailing waves into outwardly traveling leading waves.There seems little reason to expect that the modes of any other disk mass profile would differ fundamentally, as long the waves are able to propagate inside corotation and the disk is heavy enough to amplify m = 2 disturbances.Indeed, as reviewed in the introduction, we have evidence in the literature that the isochrone, Kuzmin-Toomre, and Gaussian disks all possess bar modes of the type we are exploring.Note also that our paper already explored disks of other mass profiles, because we changed the baseline exponential profile by cut outs of various depths and extents.
The cored isothermal potential is perhaps more questionable, since we know from the studies of Zang (1976) and of Evans & Read (1998) that disks in cusped potentials can be linearly stable to all modes of m ≥ 2, as long as there are no sharp edges or grooves in the disk.The global stability of such disks is precarious, however, because even quite mild non-linear effects can lead eventually to a strong bar (e.g.Sellwood 2012).We chose a cored model so as not to repeat these earlier studies.Cosmological simulations could also motivate the choice of a cusped halo mass profile, but numerous studies of galaxies have questioned whether real halos are cusped, with the balance of the evidence favoring cores (e.g.Weinberg et al. 2015;Li et al. 2020).Our cored isothermal potential also implies a halo density profile that drops as r −2 at large radii, which is slower than the r −3 power expected from cosmologically simulated halos, but this cannot affect disk stability because spherically distributed matter well outside the disk exerts no forces on the interior.Thus we doubt that our conclusion that deep disk cut outs and/or hot disk centers have little effect on global stability is dependent on our adopted potential.
Of course, this simple analytic potential is perfectly spherical, smooth, and unresponsive -properties that are not expected for real halos formed in a hierarchical universe.However, we might expect that a lumpy halo would promote disk instability, since non-axisymmetric disturbances in the disk induced by passing halo inhomogeneities would be swing-amplified and thereby contribute to the growth of bars.Also, as noted above, a halo composed of mobile particles is able to enhance bar growth, so strong instabilities in our rigid halo would grow yet more rapidly in a live one, which does not help the problem at hand.Note also that halo rotation does affect bar growth (Saha & Naab 2013;Sellwood 2015;Collier et al. 2019), but has not so far been found to stabilize the disk.These arguments do not, however, exclude an, as yet unknown, factor that may allow unbarred disks to survive in cosmological simulations.

CONCLUSIONS
All models but the last (MR3.5)presented here have turned out to be moderately to strongly unstable, and we have found little evidence from this study that surface density cut outs or hot centers contribute to disk stability.This finding is independent of changes to grid resolution, number of particles, time step, or grid geometry.2 However, the unstable mode is not always due to a standing wave in a cavity between the center and corotation, with amplification at corotation in the manner described by Toomre (1981).The modes in deeply cutout disks reflect outside the center, although these slightly different instabilities also saturate as strong bars.In both cases, the mode transforms include substantial leading components that evolve to amplified trailing waves, causing the net overall trailing spiral appearance for the mode.We therefore conclude that both types of mode rely on swing amplification, which is vigorous for bisymmetric waves in the somewhat heavy disks we adopt.
Our adopted rotation curve rises sufficiently gently that ILRs can be avoided enabling the cavity mode cycle described by Toomre (1981) in disks that do not have deep density cutouts.But those modes generally had a more tightly-wrapped spiral shape in the inner disk than is usual for bar instabilities, and we demonstrated this was because the unresponsive mass in a bulge or inner halo, required by the adopted rotation curve, raises κ and decreases the preferred wavelength in the inner disk.
Our objective in this study was to test whether the bar instability could be quelled by changing the properties of the inner disk to impede the transmission of waves through the center.However, we have found that deep inner cutouts in the disk are not stabilizing because the cavity mode in this case reflects off the sharp cutout, and we have shown that this once again leads to a strong bar.We tried more gentle cutouts, which did eventually inhibit the bar-instability, but only because the disk mass was so drastically reduced by the extensive cutout that swing-amplification of m = 2 vaves was inhibited at all radii.It is, however, rather surprising that disks having hot centers are also quite unstable.
Although our models do not exactly match those studied by Bertin et al. (1989), we do not reproduce their predicted more mild instabilities in a disk having a hot center and/or disk cutout.Neither have we seen the putative inner "refraction" from the trailing short-to long-wave branches of the dispersion relation proposed by Mark (1977); as already noted, every mode we have found had a strong leading-wave component.
The bar instability of model galaxies, for which theoretical understanding is steadily improving, persists as an unsolved problem for real galaxies.The preference in galaxies for bisymmetric spiral patterns (see Sellwood & Masters 2022, for a review) is interpreted by theorists as evidence for heavy disks (Sellwood & Carlberg 1984;Athanassoula et al. 1987), but how those same galaxies can avoid forming a bar still has no satisfactory explanation, as was highlighted by Sellwood et al. (2019) for the case of M33.Embedding the disk in a sufficiently dense halo that swing amplification of m = 2 waves is curtailed does indeed inhibit bar formation (Ostriker & Peebles 1973;Efstathiou et al. 1982), but would also suppress two armed spirals, which are the most common patterns in real galaxies (Davis et al. 2012;Hart et al. 2016;Yu & Ho 2018).The suggestion by Toomre (1981) that cutting the feedback loop of the cavity mode would stabilize the disk seemed promising, but the damping of a mode at an ILR can be overwhelmed by noise (Sellwood 1989(Sellwood , 2011)).Furthermore, this paper has, if anything, exacerbated the puzzle by showing that global stability is little affected by a dynamically hot center to the disk, or a deep cutout of responsive mass from the inner disk.
Thus the puzzle of why some galaxies lack bars remains unsolved.Although real galaxies are clearly affected by non-gravitational physics, the global stability of galaxy disks must be dominated by gravitational dynamics.We therefore believe that the absence of bars in some galaxies formed in the cosmological simulations (see §1.2) has an as yet unidentified dynamical explanation.One way forward may then be to try to identify possible stabilizing factors in those complicated simulations, and to test them one-by-one in more controlled experiments.

Figure 1 .
Figure1.The upper panel shows the function Q(R) (eq.6) and the lower panel the ratio Σc(R)/Σ(R) (eq.8).The lines in both panels use most adopted values of Qjump and D.

Figure 3 .
Figure 3.The evolution of the baseline model from 400 ≤ t ≤ 500, which forms a strong bar.The color scale shows the logarithm of the mass surface density, which initially drops to zero at R = 5R d .Notice the inner spiral and bar features, and the tightly wrapped trailing spiral near the OLR in the last two panels.

Fig. 3
Fig.3presents snapshots from the later part of the baseline simulation that show the non-linear evolution, in which perturbing forces are still restricted to m =

Figure 4 .
Figure 4. Fits to the linear evolution of two simulations in which rc = 0.75R d (left) and rc = R d (right).As in the baseline model, both have flat Q = 1.2 and no cutouts.The two logarithmic spiral transforms of the mode (upper panels) differ qualitatively, as do the mode shapes in the inner parts (lower panels).

Figure 5 .
Figure5.Fits to the linear evolution of a sequence of simulations in which the central value of Q(0) is increased; in all cases, RQ = 1.5 and there is no cutout.There was a hint of a second mode in the data for the LH panel but the fit presented is for a single mode.We have fitted two modes to the data in the upper middle panel, but present only the dominant mode in the lower middle panel.The increasing size of the mode shapes in the lower panels reflects their decreasing pattern speeds, which are given in Table2.The growth rates of the modes in HC1.0 and HC1.5 are about the same as that in the baseline model, but the mode HC2.0 (right) grows less rapidly.

Figure 8 .
Figure8.The evolution of the disk component in a rerun of 2DC0.2 in which self-consistent forces from m = 0, 2 and 4 all contributed.A large bar formed over the interval 200 ≤ t ≤ 300.Notice the inner density minimum, in which a tightly wrapped spiral developed, that persists to the last moment shown.However, the model requires a dense bulge (not shown) to support the inner rotation curve, which would make inner minimum and disk features hard to observe without sophisticated imaging and kinematic data to separate disk from bulge stars.

Figure 9 .
Figure9.Fits to the linear evolution of simulations in which the outer radius of the inner density cutout and enhanced Q is varied, while D = 0.2 and Q(0) = 3 are held fixed in all three cases.Note Rcut = RQ.

Figure 10 .
Figure 10.The radial variation of the fitted mode amplitude in our baseline model.The solid curve shows ℜ{Am(R)}, the dashed curve ℑ{Am(R)}, while the dotted curves mark the mode envelope ±|Am(R)|.The vertical lines mark the radii of the resonances: CR solid line and OLR dotted.

Figure 11 .
Figure 11.The linear mode of run 2DC0.2,which is the continuation of the sequence shown in Fig. 7 to the case for which D = 0.2.As in the rest of the sequence, Q(0) = 2 and Rcut = RQ = 1.5.The bottom panel shows the radial variation of the fitted mode amplitude, as described for Fig 10.Note that the mode amplitude is flat and near zero for R 0.3.

Figure 12 .
Figure 12.The radial variation of the inner part of |A2(R)| of the fitted mode to tests of the numerical parameters in simulations with the same initial file of particles used for run 2DC0.2(Fig 11).These simulation tests are not listed in Table2.All curves were rescaled to span 0 ≤ |A2(R)| ≤ 1. Left panel: the red and blue lines are from cases where the resolution was respectively doubled and halved from our standard grid, which is reproduced as the green line.Right panel: the red, blue and cyan lines are from cases in which ǫ was respectively halved, and increased by factors of two and four from our standard value, which is reproduced as the green line.

Figure 13 .
Figure13.The radial variation of |A2(R)| of the fitted mode as the physical parameters were changed.Again all curves were rescaled to span 0 ≤ |A2(R)| ≤ 1.The lines in the upper panel are from the indicated runs in which 1 ≥ D ≥ 0.2 in steps of 0.2 while Rcut = 1.5 and Q(0) = 2 are fixed.Notice that the mode extends into the center when D ≥ 0.6, but appears not to when D = 0.4 (cyan) and D = 0.2 (magenta).The lines in the lower panel are from runs in which Rcut = 1.5 and D = 0.2 while Q(0) is increased from Q(0) = 1.2 to Q(0) = 3.It is remarkable that the central value of Q has such a small effect on the inner edge of the mode.

Figure 14 .
Figure14.Above: the fitted mode of model MR2.5.Below: comparison between the mode amplitude profiles as Rcut is increased.Note that the inner reflection is closer to the center in MR2.5 than in 2DC0.4,for which Rcut = 1.5.

Figure 15 .
Figure15.The function Xm = 2πR/λcrit for the indicated three models.Since m = 2, we see that the swingamplification parameter X > 3 at all radii for the stable model MR3.5.

Table 1 .
Default numerical parameters