Warps and Bars from the External Tidal Torques of Tumbling Dark Halos
Abstract
The dark matter halos in CDM cosmological simulations are triaxial and highly flattened. In many cases, these triaxial equilibria are also tumbling slowly, typically about their short axes, with periods of order a Hubble time. Halos may therefore exert a slowly changing external torque on spiral galaxies that can affect their dynamical evolution in interesting ways. We examine the effect of the external torques exerted by a tumbling quadrupolar tidal field on the evolution of spiral galaxies using Nbody simulations with realistic, disk galaxy models. We measure the amplitude of the external quadrupole moments of dark halos in cosmological simulations and use these to force disk galaxy models in a series of Nbody experiments for a range of pattern speeds. We find that the torques are strong enough to induce long lived transient warps in disks similar to those observed in real spirals and also induce the bar instability at later times in some galaxy models that are otherwise stable for long periods of time in isolation. We also observe forced spiral structure near the edge of the disk where normally selfgravity is too weak to be responsible for such structure. This overlooked influence of dark halos may well be responsible for many of the peculiar aspects of disk galaxy dynamics.
Subject headings:
methods: nbody simulations – methods: numerical –cosmology: dark matter halo1. Introduction
The nearly flat rotation curves of spiral galaxies are direct evidence for the existence of massive dark halos around the galaxies, within the paradigm of NewtonianEinsteinian gravity (Roberts & Whitehurst, 1975; Rubin et al., 1979, 1982, 1985; van Albada et al., 1985; Salucci & Burkert, 2000; Sofue & Rubin, 2001). The dark halos forming in cold dark matter cosmological models have been identified with the inferred dark halos of galaxies; a substantial amount of work has gone into trying to test for the consistency between reality and simulation (e.g., Barnes & Efstathiou, 1987; Frenk et al., 1988; Dubinski & Carlberg, 1991; Navarro et al., 1996; Bullock et al., 2001; Neto et al., 2007). Most work has focused on the sphericallyaveraged density profile and the implications of this profile for galactic rotation curves. Although debate still continues on the value of the slope of the inner cusp (Gentile et al., 2004; Ferreras et al., 2007) and whether or not the inferred rotation curves are consistent with real ones (Moore, 1994; de Blok & McGaugh, 1998; de Blok et al., 2003), it seems that the consistency between the real and simulated dark halos is quite strong and that the cold dark matter (CDM) paradigm is in reasonable shape in this sense (Primack, 2007).
Other work has focused on secondary dynamical effects of dark halos that arise from their triaxial nature and their implication for other properties in disks, namely, oval distortions, warps and perhaps bars (Binney, 1978; Sparke, 1984; Franx & de Zeeuw, 1992; Sackett et al., 1994; Kuijken & Tremaine, 1994; Weinberg, 1998; Sellwood, 2003; Berentzen et al., 2006). The flattened potential of triaxial halos has been recognized in some work to induce oval distortions to otherwise circular disks, as well as vary the velocity along the orbit (Hayashi et al., 2007). The orbits of stars in the disk become flattened in a direction orthogonal to the halo density. However, the detection of such oval distortions and noncircular velocities is marginal at best, suggesting that halos may be nearly axisymmetric in the plane of an ordinary, high surfacebrightness spiral galaxy. On the other hand, the effects of noncircularity in dwarfs and lowsurface brightness galaxies may be stronger (Valenzuela et al., 2007; Hayashi & Navarro, 2006).
Another source of complex dynamical effects is the probable misalignment between a disk and a triaxial dark halo. Diskhalo misalignment has been argued to be the origin for warps seen in most edgeon spiral galaxies. A misalignment implies that the disk will feel a torque from the halo and the nearly circular orbits will therefore precess. If the disk has selfgravity and the halo potential is static, warped modes can arise (e.g., Toomre, 1983; Sparke & Casertano, 1988; Kuijken, 1991). A fatal flaw in this hypothesis is that gravitational interactions between the disk and the halo are strong, with the precessing disk experiencing dynamical friction from the halo (e.g., Nelson & Tremaine, 1995; Binney et al., 1998). One therefore expects the disk and halo to become aligned with one another within a few dynamical times and this is indeed borne out in experiments (Dubinski & Kuijken, 1995). Thus it has been suggested that it is the outer halo that is misaligned with the disk, while a tight coupling exists between the inner halo and the disk (Bailin et al., 2005; Binney, 2007).
The origin of galactic warps then remains a perplexing issue but alternative scenarios suggest a transient origin. Other ideas suggest that the cosmic infall of gas and dark matter alters the relative orientation the disk and dark halo due to addition of angular momentum (e.g., Ostriker & Binney, 1989; Debattista & Sellwood, 1999; Jiang & Binney, 1999) and the resulting torques on the disk lead to transient warping. van der Kruit (2007) also considers the possibility of the onset of a late warped outer disk, due to gas infall, in a configuration that is independent of the relatively younger inner truncated stellar disk. The infall picture has also been recently used by Shen & Sellwood (2006) to reproduce the warping of a simulated disk galaxy; in fact, their simulations corroborate the result that the lineofnode of the warp forms a leading spiral.
Tidal fields from satellite galaxies are known to excite dramatic behavior in spirals; Kalirai et al. (2006) invoke the possibility of a warp or overdensity induced by satellite interaction to explain the origin of the secondary cold population that they identify in M31 while the grand design spiral in M51 is inferred to be the result of it’s interaction with its companion NGC 5194 (Toomre & Toomre, 1972; Salo & Laurikainen, 2000). The Magellanic origin of the Galactic warp (Weinberg, 1998; GarcíaRuiz et al., 2002) is revisited by Weinberg & Blitz (2006) in which they suggest that the warp is formed as a joint handiwork of the tidal field of the Magellanic Clouds as well as the effect of the distortions that such has on the Galactic halo. The overall picture that emerges from these studies is that the source of the warping torque lies with asymmetries in the halo, at radii beyond the gravitational influence of the disk.
Dark halos formed in simulations performed within the CDM scenario suggest that not only do they often settle into triaxial figures of equilibrium but that they can also tumble, much like a rigid body (Dubinski, 1992; Bureau et al., 1999; Bailin & Steinmetz, 2004). Static and tumbling equilibrium are both perfectly valid in analogy to the two solutions for the homogeneous Jacobi ellipsoids (Chandrasekhar, 1969) which include a tumbling and static solution with internal circulation. Bailin & Steinmetz (2004) find a log normal distribution of halo pattern speeds with a mean value of about , corresponding to a tumbling period of 44 Gyr. A typical halo will then tumble through 120 degrees over a Hubble time. The tumbling periods of dark halos are rather long but nevertheless, a picture emerges in which disk galaxies are embedded within slowly tumbling, triaxial halos with rotation axes that are probably slightly misaligned with the disk rotation axes. It is likely that the inner halo and disk are aligned with each other but the outer halo (say beyond 100 kpc) is tumbling slowly and thereby affecting the disk with its tidal field. A pertinent question to ask then is how strong are the torques that a tumbling dark halo exerts on the disk at the center and whether or not these torques can have a significant influence on the evolution of a spiral galaxy over a Hubble time?
In this paper, we carry out a series of experiments to study the effect of a slowly tumbling, external quadrupolar potential, on the evolution of a disk, with the primary goal of inducing galactic warps. We first measure the expected strength of the quadrupolar tidal field directly from dark matter halos within CDM cosmological simulations and subsequently estimate the torques expected on exponential disks. We find that the typical torques from cosmological dark halos are generally an order of magnitude larger than the effect of the LMC, rendering them more effective in trigerring something dynamically interesting.
Guided by such experimentally motivated numbers, we set up simulations of ideal galactic models embedded in nearly spherical halos, forced slowly by an external quadrupolar tidal field. We discover in these experiments that it is fairly simple to make warps of amplitudes similar to the observed ones.
We also find that the gradual tidal forcing of tumbling dark halos has the effect of pushing an otherwise stable disk into a region where it becomes subject to the bar instability at late times. This is a new bar triggering mechanism and may be an additional factor that affects the evolution of the fraction of barred galaxies over cosmic history (e.g., Curir et al., 2008; Sheth et al., 2008).
The plan of the paper is as follows. In §2, we calculate the amplitude of the torque on a galactic disk, expected from the external tidal field of the triaxial dark halo, as extracted from cosmological body simulations. In §3, we present body models of M31 and the MilkyWay that are forced by the external quadrupolar field with the same strength expected from cosmological simulations. In §4, we present the results of these simulations while in §5, present analytical models of rigid disks forced by external tidal fields, to quantify these effects. In §6 we conclude with a discussion of the implications of these results for warps, bars, and spirals in real galaxies.
2. Tidal Torques on Disks from Cosmological Dark Halos
The dark halos in CDM cosmological models are highly flattened and triaxial with typical axis ratios and (e.g., Barnes & Efstathiou, 1987; Frenk et al., 1988; Dubinski & Carlberg, 1991; Warren et al., 1992; Jing & Suto, 2002). Analysis of the halo spin parameter shows a tendency for alignment with the minor axis, suggesting that a disk that forms within a dark halo is likely to have its own spin angular momentum closely aligned with the halo minor axis (e.g., Dubinski, 1992; Warren et al., 1992; Dubinski & Kuijken, 1995; Binney, 2007). The angular momentum within halos is distributed between internal streaming and tumbling motion. Early studies showed that halos are slowly tumbling through space with long periods. More recent quantitative analyses in standard CDM cosmological models show that the distribution of tumbling frequencies is roughly lognormal (Bailin & Steinmetz, 2004) and peaked at about with . The dissipative infall of the baryonic matter during galaxy formation can also modify the halo, both increasing its central concentration (e.g., Blumenthal et al., 1984) as well as rounding out the shape by increasing within the region of influence of the forming disk while leaving about the same (Dubinski, 1994; Kazantzidis et al., 2004). The axis ratios of the volumetric density contours at radii beyond the influence of the disk are less affected and remain about the same as the initial values.
Disks and halos are also found to be tightly coupled through dynamical friction in the inner regions of galaxies (Dubinski & Kuijken 1995) so one would expect the disk to lie in the principal plane of the halo, probably with the spin vectors of the disk and halo aligned (Bailin et al., 2005; Libeskind et al., 2007). However, misalignment may persist at larger radii and the slow tumbling of a halo could lead to an external torque on the disk. Although current selfconsistent simulations of galaxy formation have not fully addressed this point, we make the hypothesis that the disks are aligned with the principal planes of their inner halos through dynamical friction while the outer halos remain misaligned and may slowly be changing orientation with the tumbling frequencies measured in cosmological dark matter simulations.
We estimate this torque through the analysis of more than 2000 dark halos extracted from a cosmological dark matter simulation with near standard parameters , , and (Spergel et al., 2007). The simulation contains particles in a cube of dimension Mpc and was run using the GOTPM cosmological Nbody code (Dubinski et al., 2004). The simulation was run from using 2800 equal time steps. The comoving particle softening length was set to 3.5kpc. Halos are extracted in a two stage process; we use a friendsoffriends algorithm to determine the dense centers of candidate halos. These regions could either be associated with an isolated halo or merging pair or group. The particles in a spherical region surrounding these dense centers are extracted out to roughly the distance where the density drops to 200 times the critical density and then analyzed separately. We fit Navarro, Frenk & White (1996, NFW) profiles to all these candidates but retain only those which have fit profiles that lie within a threshold . In this way, the original sample of about 2700 halos is reduced to a sample of about 2200 halos that have welldefined profiles, undisturbed by major mergers or significant asymmetries. We then take a careful look at this sample to learn something about halo shapes and torques.
Many studies have determined the distribution of ellipsoidal shape of halos (Frenk et al., 1988; Dubinski & Carlberg, 1991; Warren et al., 1992; Cole & Lacey, 1996; Jing & Suto, 2002; Bailin & Steinmetz, 2005). We take this analysis one step further and calculate the distribution in expected strengths of tidal torques on a disk that may be misaligned with the outer halo. An explicit assumption is that the disk and halo are aligned and will remain so because of dynamical friction in the inner regions while misalignments between the outer halo will persist and change dynamically as the halo tumbles through space. We find that the relative orientations of the inner () and outer halo () is greater than 20 for half of the halos in our cosmological simulation (Fig. 1). We therefore expect significant misalignment that can lead to torques on the disk from the outer halo. We proceed by analysing the potential of the outer halo using a standard multipole expansion method.
We are interested in the possible torques at the centre of a disk tilted with respect to the outer halo. Previous work has shown that disks and inner halos are aligned within a radius roughly equal to the extent of the disk, though such alignment may extend further (Dubinski & Kuijken, 1995; Binney et al., 1998). Any torque acting on the disk from the halo probably then arises from misalignments beyond this radius. It is useful to quantify these effects more concretely. According to the NFW model, dark matter halos extend to the virial radius where is the concentration, is the scale radius of the density profile and the mean enclosed density is 200 times the critical density . At , halos have characteristic mass and circular velocity . Typical concentrations for galactic scale halos are (e.g., Bullock et al., 2001) though they can vary considerably with a scatter . A galaxy halo is usually characterized by the peak value of the rotation curve rather than . Analysis of the NFW rotation curve shows that the maximum value of the rotation curve occurs at and is related to :
(1) 
with
(2) 
Given the parameters and we find and then can determine (or ) and through the usual identities (Navarro et al., 1996):
(3) 
(4) 
As a specific example, consider a Milky Way sized NFW dark halo with with a and concentrations with the range . Assuming , these models will have virial masses in the range and scale radii from kpc. The scale radius is roughly the size of the disk and since the disk and halo are tightly coupled within this region we expect alignment. We can expect misalignment and tidal torques from the halo mass distribution for halo mass beyond the edge of the disk and so we need to calculate the external component of the potential beyond to estimate the tidal torques.
For a halo particle distribution, the external potential from matter with at a point with is given by the expression (e.g., Binney and Tremaine 1987),
(5) 
where is the usual spherical harmonic function and the coefficients are evaluated from the matter beyond through:
(6) 
with indexing a list of the particles that have .
For a dark halo, it is natural to compute these coefficients for the particles with and , as a reasonable measure of the external torquing potential of a dark halo.
Taking the lead from CMB analysis, a useful way of quantifying the strength of the halo tidal potential is through the parameter defined as:
(7) 
The monopole term will lead to zero torque and the dipole terms are probably unimportant since dark halos tend to be ellipsoidal and do not show a significant lopsided mode. The most important terms for torquing are the quadrupole terms ; one expects successive even terms to have some effect but with increasingly lower strengths. We therefore focus on the quadrupole terms that have a strength given by the parameter as the main source of torque on a misaligned disk. We now calculate this value for our sample of halos and determine its importance for disk torques.
Since the value of depends on the mass of the halo, we can rescale all halos to some reference mass. We are interested mainly in the distribution of the relative strength of the external tidal potential in typical spiral galaxies; so we first renormalize the halos to have the same scale radius and scale mass. The following formalism is then undertaken with these halos.

We first fit the scale parameters and to the NFW density profile in space with the potentialdensity pair of the spherical NFW profile given by:
(8) (9) 
Next we determine as the radius containing the mean density and find and so compute , and from the potential.

We rescale all halos to a putative Milky Way model for our determination and call the rescaled tidal parameter . We use the parameters , a concentration of which imply kpc and from the equations above.
To correlate this parameter to a specific example, we compute the expected value of the quadrupole tidal parameter for a satellite of the type of the Large Magellanic Cloud, (mass M at a distance of kpc). We refer to this tidal parameter as . The LMC has a minor tidal effect on the dynamics of the Galaxy at its current distance (e.g.,. Hunter & Toomre, 1969; GarcíaRuiz et al., 2002) but simulations of disk galaxy satellite encounters indicate that interactions in which the orbit of the LMC intersects the disk during a close encounter in the past may excite a strong tidally generated spiral (Weinberg & Blitz, 2006).
A close encounter with the LMC that brings it to half its current distance is expected to correspond to a quadrupolar tidal field that is about 8 times stronger since . Thus, the ratio is a good indicator of the active strength of the halo tidal field on Galactic dynamical evolution. If then the tidal field is comparable to the effect of the LMC on the Galaxy and so is quite weak while if , we might expect significant tidal torques on the Galaxy that may cause interesting evolution in the form of spirals, bars and warps.
Figure 2 shows the distribution of values of the parameter measured for the sample of halos using two different outer radii, and , plotted against the axis ratio . The halo axis ratios are determined using the normalized inertia tensor of halo particles within using an iterative method that gives a measure of the best fit ellipsoid of the distribution (Dubinski & Carlberg, 1991). There appears to be no strong correlation between the value of and halo shape since the scatter is quite large. The mean value of is about 5 for compared to about 25 for a cutoff at . The implication then is that a Milky Way disk that is misaligned with the outer halo will experience a tidal field that is perhaps 5 to 25 times the strength of the tidal field of an LMCtype satellite, depending on the radial location where the misalignment begins. For halos with a larger concentration , both and are smaller for the MW model but the values of vary within a factor of a few. This result is somewhat surprising for it implies that the external tidal field acting on a disk galaxy from a triaxial halo is comparable in strength to the tidal field of a large satellite having a close encounter with a galactic disk.
Finally, the disk will only be torqued if there is a significant misalignment with the halo independent of the strength of . One might imagine that more flattened halos with the largest are more closely aligned and thus have a weaker potential for torquing. We therefore looked for a correlation between and the relative alignment of the inner and outer halo . We see no such effect (Fig. 3) and conclude that significant perturbations can arise from external torquing due to this misalignment.
We now go on to demonstrate some of these effects in controlled Nbody simulations.
3. Simulations
We perform a series of body simulations using a modified version of a parallelized treecode (Dubinski, 1996). The code has been changed to include an external quadrupolar potential field that is tumbling at a fixed pattern speed . Forces on particles are the sum of the body force and the external potential field. If particles stray beyond a radius then the field is shut off. For the simulations described here, the radius is kpc. A model of a galactic disk plus surrounding dark halo is simulated within this tidal field. We also make sure the model remains centered on the quadrupole field. This prevents the outer dark matter halo in the Nbody model from being significantly distorted.
The model galaxies in question represent the Milky Way and Andromeda (M31), as given by Widrow & Dubinski (2005) (WD herein). These models both have exponential disks with Hernquisttype bulges embedded within a cuspy dark halo similar to an NFW profile. We use the most stable versions of these models called MWb and M31a in WD but we refer to them simply as MW and M31 in this paper. The M31 model has been tested at a variety of resolutions and is known to be stable against bar formation after run times of 10 Gyr (Gauthier et al., 2006). The MW model on the other hand is mildly unstable and develops a bar later on after Gyr as we shall see below. The particles are distributed with the following numbers in both the M 31 and MW runs: =1,500,000, =500,000 and =2,000,000. The model parameters are provided in Table 1.
Model  

MW  3.3  2.81  0.44  8.8  345  435  0.88  30  1  1.3 at 2.5 
M31  7.7  5.58  0.60  12.9  337  461  1.83  30  1  1.25 at 2 

Col. 1: Model for galaxy, Col. 2: Disk mass (10 M), Col. 3: Disk scale length (kpc), Col. 4: Disk scale height (kpc), Col. 5: NFWhalo scale radius (kpc), Col. 6: NFWhalo scale velocity (km s), defined in Section 2, Col. 7: characteristic bulge velocity (km s), Col. 8: bulge scale length, Col. 9: disk truncation radius (kpc), Col. 10: truncation width (kpc), Col. 11: Toomre parameter.
We also populate the disks with particles out to large radii to better investigate edge effects. The MW model extends to 10 radial scale lengths while the M31 model goes out to 7 scale lengths. We note that the scale radii of the MW and M31 halos are kpc and somewhat smaller than our putative MW halo in the discussion above. The parameters are fitted values for observed rotation curves and inferred brightness profiles of the two galaxies. The real inner halos of M31 and MW are likely contracted compared to the pure dark matter models. The mass model fits to the real data have scale radii that are smaller by a factor of two in accordance with expectations of the contraction of halos by dissipative processes during galaxy formation.
We have found it convenient to select a set of simulation units that set the scale lengths and rotation velocities to values near unity. We scale the simulation units to the physical units in the following way:

Gravitational constant =1,

1 simulation length unit = 4 kpc

1 simulation velocity unit = 220 km s,

1 simulation mass unit = 4.5 M,

1 simulation time unit = 17.7 Myr
We express the general external tidal field as:
(10) 
The coefficients and describe the strength of the field and are derived from Equation 6. If we assume that the quadrupole is aligned with the principle axes and is due to a planesymmetric distribution then the terms are zero. For these simulations we are only interested in the quadrupole terms. The pure external quadrupole field then becomes:
(11) 
We estimate a suitable tidal field using the following procedure. We first take an NFW profile that has been flattened into a perfect ellipsoidal shape with axis ratios and . We then determine the spherically averaged density profile of this distribution and fit a spherical NFW profile to determine the scale radius, and mass as is done in cosmological simulations. We are interested in the tidal field due to material beyond a radius, , so we compute the coefficients of the quadrupolar tidal field for the material beyond this radius. As an example, we find that for a flattened NFW halo with and with a scale length kpc and that the quadrupolar tidal field coefficients in simulation units are:

= 0.00050,

= 0.0,

= 0.00015.

The quadrupole amplitude is .
(If we used a more concentrated halo with then kpc and leading to a field strength that would be about 3 times larger for .) The Nbody galaxy models have a nearly spherical halo at large radii so we model the effect of the triaxial halo by this quadrupolar field. We assume that the inner halo has aligned with the disk and is essentially axisymmetric and remains aligned for the course of the simulation. As a comparison, the quadrupole coefficients for a LMClike satellite of M, at a distance of 50 kpc placed on the axis, is in simulation units so that close to expectation value for the cosmological distribution of halos for material with ..
To introduce a misalignment, the disk is tilted by 30 with respect to the field and the field is made to rotate about the axis with different pattern speeds. This choice of the tilt angle is arbitrary but is approximately the expected angle for at least 30% of dark halos, according to the distribution shown in Figure 1. The 4 pattern speeds that we choose are and 0.88 corresponding to tumbling periods of and 7 Gyr respectively. We note that these are very slow compared to the pattern speeds of barred galaxies () and so the corotation radii are mainly larger than (The circular frequency at in the NFW model is .) These periods are approximately 4.0, 2.0, 1.0 and 0.5 times the age of the universe. We recall that Bailin & Steinmetz (2004) find a lognormal distribution peaking at , corresponding to a tumbling period of 44 Gyr. Thus, our simulations probe the fast tumbling tail of this distribution – the half of the distribution that is most likely to yield interesting dynamical effects. The quadrupole coefficients and pattern speeds for the MW and M31 runs are summarized in Table 2.
Model  

MWcon  0.0  0.0  0.0  
MW0  5.0  1.5  0.11  56 
MW1  5.0  1.5  0.22  28 
MW2  5.0  1.5  0.44  14 
MW3  5.0  1.5  0.88  7 
MW2a  2.5  0.75  0.44  28 
MW2b  1.3  0.38  0.44  28 
MW2c  0.63  0.19  0.44  28 
M31con  0.0  0.0  0.0  
M310  5.0  1.5  0.11  56 
M311  5.0  1.5  0.22  28 
M312  5.0  1.5  0.44  14 
M313  5.0  1.5  0.88  7 

Col 1: model name, Col 2: component of the external quadrupolar field (10 Ms), Col 3: component of the external quadrupolar field (10 Ms), Col 4: halo pattern speed (km s kpc), Col 5: tumbling period (Gyrs)
We simulate each run for 8000 equal timesteps with units, corresponding to about 7.1 Gyr in physical units. We soften gravity with a Plummer softening kernel with radius of 40 pc for the stars and 100 pc for the dark matter. Energy is conserved typically to within 0.5% and total angular momentum to within 1%.
4. Results
Before exploring the results for externally torqued disks, we first examine the evolution of control models. Figure 4 shows the final state for the MW and M 31 models in the absence of any external torquing potential. The disks remain coplanar and there are minimal signs of disk thickening and warping resulting from the amplification of the Poisson noise in the Nbody disk. There are some remnants of spiral structure in the disks which arise from swing amplification of the noise in the disk. This transient spiral structure heats the disk azimuthally and slowly fades away as the simulations progress.
When we turn on the external quadrupole potential with the expected amplitudes from cosmological dark halos, we clearly see the various effects predicted for disk torquing. Figures 5 and 6 show the evolution of representative models of the MW and M 31 runs from 3 perpendicular views. Associated videos^{1}^{1}1Animations are available at the website http://www.cita.utoronto.ca/dubinski/warpmovies also show how the external torque causes the disk to precess and nutate like a tilted gyroscope on a table top. The evolution of the direction of the spin vector of the disks described by the normalized and components of the angular momentum quantifies this behavior (Fig. 7) and the nutational frequency depends on the pattern speed of the forcing quadrupole potential. All runs indicate that the model disk precesses through an angle of 3060 over 7 Gyr as it responds to the external quadrupolar field implying precession periods of more than Gyr or precession rates less than 0.15 rad Gyr. We see below that these precession rates are in accord with analytical estimates from a rigid disk model forced by an external potential. These precession rates are about an order of magnitude smaller than those found in models with a disk misaligned within a triaxial potential (e.g., Kuijken, 1991; Jeon et al., 2009) used in models of warping behavior that postulate that the main source of torque on a disk is due to he potential of the inner part of the triaxial halo. The phenomenology that we are exploring is therefore different and the resulting effects are slower and more subtle. So while the precession periods are long compared to the Hubble time, a disk can slew through and a significant angle over its lifetime and this dynamical evolution leads to a more gradual transient warping at the disk edge in a mechanism similar to those invoking the accretion of angular momentum (Ostriker & Binney, 1989). We conclude that this must be an important process in disk dynamical evolution over the age of the universe. For real galaxies, we expect a range of halo tumbling pattern speeds, tidal field strengths and initial tilt angles so there should be also be a range of warping behavior in the population at large.
4.1. Warps
While the the disk stars within approximately 5 radial scalelengths stay within the plane and act like a rigid body because of their selfgravity, stars near the edge of the disk begin to precess at different frequencies and begin to warp away from the disk. Our simulations indicate that the model disks develop strong warping early on in the runs and that such warping is sustained through the length of the run (about 7 Gyrs) (Fig. 8).
The disk demonstrates strong flocculent spiral pattern, to even the diskedge where selfgravity is considered too weak to excite such structure. This is advanced as the handiwork of the external quadrupole since we do not notice prominent outer spiral patterns in the control runs. This is particularly noticeable in the MW runs that have a more extended disk.
Following Levine et al. (2006), we quantify the disk warps that develop in the simulations by analyzing the vertical deviation of the disk from a plane at different radii through the function,
(12) 
where is the particle azimuth.
The inner, nearly coplanar disk is considered to be in the plane. The particles inside the ring  of radius and width  are assigned the height from the midplane with defined as the corresponding particle azimuth. We first rotate into the frame of the inner disk before the analysis. We determine the midplane of the disk by computing the moment of inertia of the disk particles within 3 scale lengths where the disk is nearly planar and then diagonalizing the inertia tensor to give the orientation of the disk. We determine the coefficients using a Fourier analysis of the particle distribution through,
(13) 
The amplitude of the deviation is given by and the phase angle through . We can then construct functions of versus to examine disk warping. The function corresponds to bowlshaped deviations, corresponds to the usual integralsign shaped warps while are second order “scalloped” warps.
Figure 9 shows the results of this analysis for the 4 models of both the MW and M31 system. The dominant vertical deviation is the integralsigned warp. For the MW models, the warp turns upwards at kpc and continues to kpc reaching about 5 kpc above the plane. The M 31 models shows a qualitatively similar behavior though the warping begins at the larger radius of kpc but only reaches about 3 kpc above the plane at a radius of kpc.
The M31 model is more extended and massive than the MW model while being forced by the same external tidal field and so shows a smaller degree of warping behavior.
We also examined the alignment of the warps by measuring the line of nodes (LON) from the phase angle determined from the vertical Fourier analysis. We present the tipLON plots for our models following Briggs (1990) (Fig. 10). The warps begin with a relatively straight line of nodes when they begin to move out of the plane at the disk edge but generally show a leading spiral curving out to following the observations from Brigg’s analysis. At larger radii, the warping becomes less coherent resulting mainly by the onset of differential precession of the stellar disk orbits and the modulated forcing of the tumbling external quadrupole potential. While we are only modeling a pure Nbody system one expects that the mechanism will generate coherent warping in the gas as well at least out to 6 scale lengths.
4.2. Bars
The tidal distortion due to a galaxy interaction in a flyby can trigger the bar instability under some conditions (e.g. Noguchi, 1987; Gerin et al., 1990). The effect is due to a transient external quadrupole that disturbs the stellar orbits in the center of the galaxy. Under our hypothesis of a tumbling misaligned halo, there is also a changing external quadrupole potential that can perturb the disk and potentially trigger the bar. The physical situation is somewhat different than a impulsive galactic flyby however, in that rate of change of the orientation is slow and the amplitude of the tidal potential is constant. Nevertheless, the Milky Way disk is found to develop a fast bar at about 3 Gyrs. Once formed, the bar is sustained and its pattern speed declines from an initial value of 50 to 35 by the end of the run in accord with typical expectations for models like this (Fig 11) (e.g., Debattista & Sellwood, 1998; O’Neill & Dubinski, 2003; Dubinski et al., 2009).
Figure 12 shows the evolution of the strength of the bar that is triggered in the Milky Way for different models by the torquing provided by the external halo. We find a welldeveloped bar, of nearly similar strengths, forming by about 3 Gyrs in all the runs, indicating that the rate of tumbling of the halo is immaterial to the bar strength. However, of the quadrupole appears to have a weak effect on the time of onset of the bar instability  this is hastened with increasing tumbling period , i.e. falling pattern speed of the external halo, though this is not marked; the fractional difference in bar onset times between the runs done with the slowest and fastest halo in the Milky Way model ( and , respectively) is about 16. This small difference in the bar onset time is also noticed in the snapshots from these two runs, shown in Figure 5.
Instead, we find that it is the strength of the quadrupole of the external halo that is the crucial factor in controlling the bar instability. In fact, a control run performed with a null external quadrupole leads to a much later ( 5.8 Gyrs) onset of a very weak bar in the Milky Way disk (maximum bar strength is less than 0.4 times the bar strength reached in the other runs); see Figure 4. Thus, we assure ourselves that the inherent bar unstable nature of the Milky Way model is very weak, albeit nonzero. Additionally, the control runs for both models result in significantly weaker spirality, especially at the edges of the disk and no warp.
The M 31 disk is found to be unaffected by the halo torquing, in regard to bar formation. This is apparent from the pictures of the M 31 disk at the end of the simulations (Figure 6); the result is found independent of the halo pattern speed used in the run.
5. Analytic Rigid Disk Model
It is interesting to gauge the effects of the external quadrupolar term in the potential of a tumbling triaxial dark halo in an analytical model using a rigid, thin exponential disk for comparison to the phenomenology in the simulations. We emphasize again that the inner disk and halo are tightly coupled through dynamical friction and so remain aligned. The only possible source of torquing on the disk is a misaligned and possibly tumbling outer halo. The assumption of a rigid disk is a reasonable assumption since the simulations show that the disk maintains coherence because of selfgravity within even while it tips in response to the external tidal torque. We use the classical Euler equations of motions for a rigid spinning disk. The dynamic variables are the standard Euler angle () where is the inclination of the disk, is the longitude of ascending node and the angle about the disk’s symmetry axis. In the presence of a time dependent torque, the disk will precess and nutate and so and will evolve in time. In these calculations, we set up the disk with an initial inclination of and set the halo quadrupole potential tumbling about the axis with for direct comparison to the results of the simulations..
5.1. Euler Equations of Motion
After Goldstein et al. (2002) (see p.210 and equation 5.52), the Lagrangian of the disk is
(14) 
where is the moment of inertia about the symmetry axis are the moments measured about an axis in the plane of the disk. The potential energy is determined from the interaction between a rigid thin exponential disk and a tumbling external quadrupole potential. The equations of motion can then be written as:
(15)  
(16) 
where (the disk angular momentum) is a constant of motion and is the potential energy between the disk and the tumbling halo. A tilted disk will naturally precess. For example, with an axisymmetric potential , the expected precession frequency can be derived from the first Euler equation by assuming that is constant. For small we expect
(17) 
(cf. Goldstein et al., 2002). If the potential also depends on and is rotating about the axis, we can expect a more complex evolution combining both precession and nutation.
5.1.1 Direct Solution of the Euler Equations for a Rigid Exponential Disk
We first consider the interaction of a rigid exponential disk with an external quadrupolar tidal field for comparison to the results of the Nbody simulations. We point out that this process is quite different from the disk precession examined in work that has treated a triaxial halo as a rigid background potential (Kuijken, 1991; Jeon et al., 2009). In those treatments, the precession rate is an order of magnitude larger due to the much larger torques felt by the inner halo. However, the large torque in these simulations is artificially high since we expect dynamical friction to align the disk and inner halo within a few disk dynamical times (Dubinski & Kuijken, 1995). In this study, we assume that only the outer halo can exert any torque and subsequently the strength of the torque is considerably smaller, typically by an order of magnitude.
In Appendix A, we compute the time dependent interaction potential between a rigid exponential disk and an external quadrupole potential to plug into the Euler equations. The potential is:
(18)  
where and are the mass and exponential scale radius of the disk, and are the quadrupole terms for the external terms of a triaxial potential in equation 11, is the tumbling frequency of the quadrupole potential assumed to be about the axis with and we assume that at . We can compute the moments of inertia and and the spin directly from the Nbody model and so fully specify the Euler equations for our case.
For the M31 model, we solve the Euler equations of motion numerically using the same initial state as the Nbody model along with four different pattern speeds for the quadrupole potential. We use the values of and expected for an exponential disk as well as the spin (disk angular momentum) measured from the Nbody models. The specific parameters for the M31 model in dimensionless units are , so that , and . Using equation 17, we can estimate the precession period assuming that timeaverages to zero for sufficiently large . We find a value of which is comparable in magnitude to our chosen pattern speeds but opposite in sign. This halo corresponds to our choice of of about 6.3. We then integrate the Euler equations for approximately 53 Gyr corresponding to a little more than a typical precession period to get a complete picture of the evolution of the disk for comparison to the Nbody simulations. Figure 13 shows the evolution of the spin axis of the disk for comparison to the similar plots for the Nbody simulation (Fig. 7). As in Fig. 7, Figure 13 presents this evolution by tracing the normalized, component of the disk angular momentum against the component of the same.
The agreement in behavior for the four pattern speeds is excellent suggesting that the Nbody disk is behaving in a similar manner to a rigid body over this time period. The nutation induced by the triaxiality is readily apparent again in this plot. We notice that the nutation period is shorter for a higher pattern speed. Thus, for larger pattern speeds, as the nutation period is small, we might expect a spiral galaxy to wobble several times over its lifetime if it is embedded in a rapidly tumbling dark halo. For the slower pattern speeds, the nutation period becomes comparable to the precession period and the resulting effect is simply a steady change of the disk orientation. Stars and gas on the outer edge of the disk that have weaker selfgravity are prone to precess differentially and this may be an additional origin of galactic disk warps.
5.2. Resonant Interactions
For the flattened halo discussed here the precession rate of disk is negative while the tumbling rate is positive and in the same sense as the disk. While this case is the most reasonable, it is possible that the sign of the precession rate and halo tumbling rate might be the same and so resonances can occur when . In the first investigations of orbits in a rotating triaxial potential, Binney (1978, 1981) discovered that under certain conditions retrograde orbits in the rotating frame can be unstable leading to large motions due to resonant interactions with the potential and suggested that this might be a mechanism to produce galactic warps. Heisler et al. (1982) demonstrated the existence of Binney’s resonance in orbital studies using an analytic triaxial potential. Tremaine & Yu (2000) have even suggested that this resonance may provide a mechanism for creating polar ring galaxies. Accreted gas and forming stars counterrotating with respect to main disk may be levitated out of the plane if there is an underlying rotating triaxial potential.
So far the study in this paper has elucidated how selfconsistent Nbody disks evolve in tumbling triaxial potential rotating in the same direction as the disk. In the other case of a counterrotating disk, it is interesting to see if Binney’s resonance manifests itself in disk evolution in our models.
Let us first consider the simple rigid disk model. Consider the Euler equations again for the case where and corresponding to steady precession at a fixed angle at a constant rate. For this to be the case in general, we need to find the condition such that at all times. Differentiating equation 18 with respect to and setting it to zero implies that . If then this implies a resonance when . With this resonant condition, we can use the first Euler equation to find from the quadratic equation
(19) 
where is known from equation A25 (see Appendix 1). For a flat rotation curve with speed the disk angular momentum is given by . If we assume that is small then to a good approximation it is given by
(20) 
For our model of M31, with our estimate of and for the coefficients of an external quadrupole potential and an initial disk inclination of , the resonance occurs when . Thus we see that the size of this typical resonant precession frequency is comparable to the expected tumbling frequencies of dark halos.
We can examine the behavior of the disk precessional evolution near this resonance. In Fig. 14, we represent the disk spin axis evolution to bring out the unstable behavior of rigid models with set to and . The radically different nearresonance disk behavior noted in this figure, marks the position of a region of instability. Here the disk inclination can change by large amount and the precession and nutational behavior are quite irregular.
We also consider an additional M31 Nbody model for comparison to the rigid disk model. We evolve it in the same way as previous models but here examine a counterrotating external quadrupole potential with to see if the irregular behavior that is seen in the simple rigid disk model is reproduced in a selfconsistent Nbody model. Figure 15 shows a comparison of the time evolution of the polar angles representing the direction of the disk angular momentum vector for the rigid disk and Nbody disk models. (Here is offset by 90 degrees from the usual definition as the angle of the line of nodes in the Euler equation.) For the rigid disk, we adjust the initial inclination slightly to and find excellent agreement with the Nbody disk evolution. We have integrated for the unusually long time of 25 Gyr to see the full tilt evolution. However, the amplitude of the quadrupole potential could easily be twice as large based on the observed variance of in cosmological halos and so this behavior would occur in half the time since from equation 20 we see that the precession rate is directly proportional to the external quadrupole amplitude. It is remarkable that in this counterrotating case, the disk inclination grows by 90 before reversing course.
While the Binney resonance is derived in the context of closed collisionless orbits in rotating triaxial potentials, a version of it also appears to be operating collectively in a selfgravitating disk in the same context. Our idealized model of a Nbody disk forced by a rotating external quadrupole potential tips through 90 when counterrotating tumbling frequencies are near the resonant value. The case of counterrotation is probably rare but not impossible given the existence of some disks with counterrotating components (e.g., Rubin et al., 1992). It may also help explain the phenomenon of polar rings as suggested by Tremaine & Yu (2000). Accreted gas that is counterrotating with respect to the tumbling halo might rise out of the halo principal plane within a Hubble time given a large enough quadrupole amplitude and suitable tumbling frequency. One caveat to be aware of is that given the long timescales of this process, dynamical friction may begin to align even the outer parts of the halo so the resulting torque may diminish with time and so weaken the effect.
6. Conclusions
In this paper we discuss the effect of the tidal field of a tumbling and triaxial external halo on the dynamical evolution of a galactic disk. We have shown that the disk precesses in response to its misalignment with the external halo and the gradual realignment of the disk with the external tumbling halo causes it to warp. The disk is also seen to develop strong spirality to its edge in the potential of the external halo that it sits in. It is also suggested that if a disk feels a significant external quadrupole that is essentially static or at most slowly tumbling, it might trigger a bar instability, late in the disk.
The conclusions achieved in this work are based on analytical calculations that are backed up by Nbody simulations of the Milky Way and M 31 models, subject to the quadrupolar tug of an external tumbling and triaxial halo. The quadrupolar strength is judged via cosmological simulations, while a range of halo tumbling periods are scanned through in the simulations. We advance this dynamical mechanism as an important cause of nonaxisymmetric structures to set in galactic disks. An important result of our work is that the prodding of the disk by a external quadrupole may not always produce bars in the disks. Thus, the M 31 disk was found to be robust to the bar instability, though the Milky Way disk was found to be much more susceptible to the bar being triggered by this mechanism. However, both disks were found to warp. The dispersion in cosmological halo properties imply that the external quadrupole is rather weak in some systems and so may not be effective in driving disk evolution.
Our work elucidates the importance of the quadrupolar term in the potential of the external halo, particularly, in terms of the trigerring of the bar and warp instabilities. Given the viability of such tidal prodding of the disk by the external halo, detailed investigation of this mechanism in galaxy formation simulations is suggested.
Acknowledgments
We acknowledge help from the CITA visitor programme which aided in this collaboration. We thank an anonymous referee for detailed comments that greatly improved this paper. We acknowledge NSERC for financial assistance and SHARCNET for providing supercomputer time. DC acknowledges the support of a Royal Society Dorothy Hodgkin Research Fellowship.
Appendix A
We calculate the interaction potential between an exponential disk and a timedependent external quadrupole potential rotating with pattern speed . If we first assume the potential is static and aligned with the body axes then we can write the potential from equation 11 in Cartesian coordinates as:
(A1) 
with
(A2)  
(A3)  
(A4) 
where one can show easily that the coefficients , and are defined in terms of the harmonic expansion in equation 11 through:
(A5)  
(A6)  
(A7) 
If the potential is rotating counterclockwise about the axis with pattern speed then we can introduce the timedependence of the potential into the coordinates using inertial frame coordinates where
(A8)  
(A9)  
(A10) 
so that potential becomes
(A11) 
with the coefficients
(A12)  
(A13)  
(A14)  
(A15) 
We can thus derive the timedependent interaction potential between a tilted exponential disk and this tumbling quadrupole potential. The parametric equation of a ring of radius with inclination angle and line of nodes is given by:
(A16)  
(A17)  
(A18) 
where is the angle measured along the ring. The mass element of the ring is given by
(A19) 
so we can derive the interaction potential by integrating over and through
(A20) 
If we assume an exponential disk with
(A21) 
and substitute the expressions for in equations A18 and do the integrals we arrive at the intermediate form of the interaction potential:
(A22) 
This can be reduced to the following simpler form in terms of the original quadrupole coefficients and :
(A23) 
The derivatives of the interaction potential used in the Euler equations are then:
(A24)  
(A25) 
References
 Bailin et al. (2005) Bailin, J., Kawata, D., Gibson, B. K., Steinmetz, M., Navarro, J. F., Brook, C. B., Gill, S. P. D., Ibata, R. A., Knebe, A., Lewis, G. F., & Okamoto, T. 2005, ApJLett., 627, 17
 Bailin & Steinmetz (2004) Bailin, J., & Steinmetz, M. 2004, ApJ, 616, 27
 Bailin & Steinmetz (2005) —. 2005, ApJ, 627, 647
 Barnes & Efstathiou (1987) Barnes, J., & Efstathiou, G. 1987, ApJ, 319, 575
 Berentzen et al. (2006) Berentzen, I., Shlosman, I., & Jogee, S. 2006, ApJ, 637, 582
 Binney (1978) Binney, J. 1978, MNRAS, 183, 779
 Binney (1981) —. 1981, MNRAS, 196, 455
 Binney (2007) —. 2007, Dynamics of Disks (Island Universes  Structure and Evolution of Disk Galaxies), 67–+
 Binney et al. (1998) Binney, J., Jiang, I.G., & Dutta, S. 1998, MNRAS, 297, 1237
 Blumenthal et al. (1984) Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Nature, 311, 517
 Briggs (1990) Briggs, F. H. 1990, ApJ, 352, 15
 Bullock et al. (2001) Bullock, J. S., Kolatt, T. S., Sigad, Y., Somerville, R. S., Kravtsov, A. V., Klypin, A. A., Primack, J. R., & Dekel, A. 2001, MNRAS, 321, 559
 Bureau et al. (1999) Bureau, M., Freeman, K. C., Pfitzner, D. W., & Meurer, G. R. 1999, AJ, 118, 2158
 Chandrasekhar (1969) Chandrasekhar, S. 1969, Ellipsoidal figures of equilibrium (The Silliman Foundation Lectures, New Haven: Yale University Press, 1969)
 Cole & Lacey (1996) Cole, S., & Lacey, C. 1996, MNRAS, 281, 716
 Curir et al. (2008) Curir, A., Mazzei, P., & Murante, G. 2008, A&A, 481, 651
 de Blok et al. (2003) de Blok, W. J. G., Bosma, A., & McGaugh, S. 2003, MNRAS, 340, 657
 de Blok & McGaugh (1998) de Blok, W. J. G., & McGaugh, S. S. 1998, ApJ, 508, 132
 Debattista & Sellwood (1998) Debattista, V. P., & Sellwood, J. A. 1998, ApJ, 493, L5+
 Debattista & Sellwood (1999) —. 1999, ApJ, 513, L107
 Dubinski (1992) Dubinski, J. 1992, ApJ, 401, 441
 Dubinski (1994) —. 1994, ApJ, 431, 617
 Dubinski (1996) —. 1996, New Astronomy, 1, 133
 Dubinski et al. (2009) Dubinski, J., Berentzen, I., & Shlosman, I. 2009, ApJ, 697, 293
 Dubinski & Carlberg (1991) Dubinski, J., & Carlberg, R. G. 1991, ApJ, 378, 496
 Dubinski et al. (2004) Dubinski, J., Kim, J., Park, C., & Humble, R. 2004, New Astronomy, 9, 111
 Dubinski & Kuijken (1995) Dubinski, J., & Kuijken, K. 1995, ApJ, 442, 492
 Ferreras et al. (2007) Ferreras, I., Saha, P., Williams, L. L. R., & Burles, S. 2007, ArXiv eprints, 708
 Franx & de Zeeuw (1992) Franx, M., & de Zeeuw, T. 1992, ApJ, 392, L47
 Frenk et al. (1988) Frenk, C. S., White, S. D. M., Davis, M., & Efstathiou, G. 1988, ApJ, 327, 507
 GarcíaRuiz et al. (2002) GarcíaRuiz, I., Kuijken, K., & Dubinski, J. 2002, MNRAS, 337, 459
 Gauthier et al. (2006) Gauthier, J.R., Dubinski, J., & Widrow, L. M. 2006, ApJ, 653, 1180
 Gentile et al. (2004) Gentile, G., Salucci, P., Klein, U., Vergani, D., & Kalberla, P. 2004, MNRAS, 351, 903
 Gerin et al. (1990) Gerin, M., Combes, F., & Athanassoula, E. 1990, A&A, 230, 37
 Goldstein et al. (2002) Goldstein, H., Poole, C., & Safko, J. 2002, Classical mechanics (AddisonWesley)
 Hayashi & Navarro (2006) Hayashi, E., & Navarro, J. F. 2006, MNRAS, 373, 1117
 Hayashi et al. (2007) Hayashi, E., Navarro, J. F., & Springel, V. 2007, MNRAS, 377, 50
 Heisler et al. (1982) Heisler, J., Merritt, D., & Schwarzschild, M. 1982, ApJ, 258, 490
 Hunter & Toomre (1969) Hunter, C., & Toomre, A. 1969, ApJ, 155, 747
 Jeon et al. (2009) Jeon, M., Kim, S. S., & Ann, H. B. 2009, ApJ, 696, 1899
 Jiang & Binney (1999) Jiang, I.G., & Binney, J. 1999, MNRAS, 303, L7
 Jing & Suto (2002) Jing, Y. P., & Suto, Y. 2002, ApJ, 574, 538
 Kalirai et al. (2006) Kalirai, J. S., Guhathakurta, P., Gilbert, K. M., Reitzel, D. B., Majewski, S. R., Rich, R. M., & Cooper, M. C. 2006, ApJ, 641, 268
 Kazantzidis et al. (2004) Kazantzidis, S., Kravtsov, A. V., Zentner, A. R., Allgood, B., Nagai, D., & Moore, B. 2004, ApJ, 611, L73
 Kuijken (1991) Kuijken, K. 1991, ApJ, 376, 467
 Kuijken & Tremaine (1994) Kuijken, K., & Tremaine, S. 1994, ApJ, 421, 178
 Levine et al. (2006) Levine, E. S., Blitz, L., & Heiles, C. 2006, ApJ, 643, 881
 Libeskind et al. (2007) Libeskind, N. I., Cole, S., Frenk, C. S., Okamoto, T., & Jenkins, A. 2007, MNRAS, 374, 16
 Moore (1994) Moore, B. 1994, Nature, 370, 629
 Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
 Nelson & Tremaine (1995) Nelson, R. W., & Tremaine, S. 1995, MNRAS, 275, 897
 Neto et al. (2007) Neto, A. F., Gao, L., Bett, P., Cole, S., Navarro, J. F., Frenk, C. S., White, S. D. M., Springel, V., & Jenkins, A. 2007, MNRAS, 381, 1450
 Noguchi (1987) Noguchi, M. 1987, MNRAS, 228, 635
 O’Neill & Dubinski (2003) O’Neill, J. K., & Dubinski, J. 2003, MNRAS, 346, 251
 Ostriker & Binney (1989) Ostriker, E. C., & Binney, J. J. 1989, MNRAS, 237, 785
 Primack (2007) Primack, J. R. 2007, Nuclear Physics B Proceedings Supplements, 173, 1
 Roberts & Whitehurst (1975) Roberts, M. S., & Whitehurst, R. N. 1975, ApJ, 201, 327
 Rubin et al. (1985) Rubin, V. C., Burstein, D., Ford, Jr., W. K., & Thonnard, N. 1985, ApJ, 289, 81
 Rubin et al. (1982) Rubin, V. C., Ford, Jr., W. K., Thonnard, N., & Burstein, D. 1982, ApJ, 261, 439
 Rubin et al. (1992) Rubin, V. C., Graham, J. A., & Kenney, J. D. P. 1992, ApJ, 394, L9
 Rubin et al. (1979) Rubin, V. C., Roberts, M. S., & Ford, Jr., W. K. 1979, ApJ, 230, 35
 Sackett et al. (1994) Sackett, P. D., Rix, H.W., Jarvis, B. J., & Freeman, K. C. 1994, ApJ, 436, 629
 Salo & Laurikainen (2000) Salo, H., & Laurikainen, E. 2000, MNRAS, 319, 393
 Salucci & Burkert (2000) Salucci, P., & Burkert, A. 2000, ApJLett., 537, L9
 Sellwood (2003) Sellwood, J. A. 2003, ApJ, 587, 638
 Shen & Sellwood (2006) Shen, J., & Sellwood, J. A. 2006, MNRAS, 370, 2
 Sheth et al. (2008) Sheth, K., Elmegreen, D. M., Elmegreen, B. G., Capak, P., Abraham, R. G., Athanassoula, E., Ellis, R. S., Mobasher, B., Salvato, M., Schinnerer, E., Scoville, N. Z., Spalsbury, L., Strubbe, L., Carollo, M., Rich, M., & West, A. A. 2008, ApJ, 675, 1141
 Sofue & Rubin (2001) Sofue, Y., & Rubin, V. 2001, ARA&A, 39, 137
 Sparke (1984) Sparke, L. S. 1984, MNRAS, 211, 911
 Sparke & Casertano (1988) Sparke, L. S., & Casertano, S. 1988, MNRAS, 234, 873
 Spergel et al. (2007) Spergel, D. N., Bean, R., Doré, O., Nolta, M. R., Bennett, C. L., Dunkley, J., Hinshaw, G., Jarosik, N., Komatsu, E., Page, L., Peiris, H. V., Verde, L., Halpern, M., Hill, R. S., Kogut, A., Limon, M., Meyer, S. S., Odegard, N., Tucker, G. S., Weiland, J. L., Wollack, E., & Wright, E. L. 2007, ApJS, 170, 377
 Toomre (1983) Toomre, A. 1983, in IAU Symposium, Vol. 100, Internal Kinematics and Dynamics of Galaxies, ed. E. Athanassoula, 177–185
 Toomre & Toomre (1972) Toomre, A., & Toomre, J. 1972, ApJ, 178, 623
 Tremaine & Yu (2000) Tremaine, S., & Yu, Q. 2000, MNRAS, 319, 1
 Valenzuela et al. (2007) Valenzuela, O., Rhee, G., Klypin, A., Governato, F., Stinson, G., Quinn, T., & Wadsley, J. 2007, ApJ, 657, 773
 van Albada et al. (1985) van Albada, T. S., Bahcall, J. N., Begeman, K., & Sancisi, R. 1985, ApJ, 295, 305
 van der Kruit (2007) van der Kruit, P. C. 2007, A&A, 466, 883
 Warren et al. (1992) Warren, M. S., Quinn, P. J., Salmon, J. K., & Zurek, W. H. 1992, ApJ, 399, 405
 Weinberg (1998) Weinberg, M. D. 1998, MNRAS, 297, 101
 Weinberg & Blitz (2006) Weinberg, M. D., & Blitz, L. 2006, ApJLett., 641, L33
 Widrow & Dubinski (2005) Widrow, L. M., & Dubinski, J. 2005, ApJ, 631, 838