,
(4-1)
To investigate the effect of the orbiting satellite, I employ an N-body particle simulation, which represents the galaxies as collections of discrete particles. Although this sounds reasonable (galaxies are collections of stars), in practice the number of particles used to represent a galaxy (~105) is a tiny fraction of the actual number of stars (~1011). Therefore, instead of representing stars, particles serve to sample phase-space, providing an average value for a region of space. This is acceptable because the motion of stars in a galaxy are primarily determined by the mean potential, and not by two-body interactions; a collisionless system. Regardless, this averaging causes a slow accumulation of error, limiting the accuracy of the simulation. While increasing the number of particles mitigates this error, a greater number of particles causes the simulation to run more slowly.
All particle models are realized as binary files containing each particle's mass, position, and velocity. My simulations use a total of 152,000 particles: 100,000 in the halo; 10,000 in each the bulge and satellite; and 32,000 in the disk. Running on a timeshared DEC Alpha workstation, a 3 Gyr run requires roughly 24 hours of computation.
The algorithm employed is a hybrid self-consistent field (SCF) and tree-code written by John Dubinski (Dubinski & Nelson 1996). While the algorithm's specifics are beyond the scope of this paper, the SCF code (Hernquist & Ostriker 1992) calculates the potential of the spherical components, the satellite, halo, and bulge, using a low order expansion of a basis set which well represents R1/4 density distributions. The tree-code (Barnes & Hut 1986) calculates the potential of the disk, running in O(N lg N) time. The SCF method is very efficient, running in O(N) time, but its use is limited to situations in which the components retain a uniform density distribution. Due to the majority of particles in spherical components, the hybrid algorithm runs much faster than a pure tree-code.
Unlike standard N-body simulations, which use discrete particles that are free in all dimensions, Dubinski's (1996) model represents the disk with a collection of rigid concentric rings centered on a common origin. Only the origin may move in space, individual rings are free to rotate about the origin. Forces acting on each ring correspond to orbit-averaged forces. The advantages of a ring model are both that the ring-disk is sensitive to slight perturbations and that the characteristics of the ring-disk are easily analyzed. Instead of fitting density contours to a collection of particles, inclination and phase are immediately obtained from each ring's rotation matrix. The disk space required to represent a ring model is also much smaller than a corresponding particle model.
The ring system evolves in the following manner. At each timestep, the simulation code decomposes each ring into a number of particles (32 in these simulations), spaced evenly about the ring's perimeter. The SCF code uses the spherical component's density to find their potentials. The treecode computes the potential of the disk. The code then finds the net force on all particles in the simulation (Equation C-5). With normal particles, the code integrates the force using leapfrog integration to find the new position and velocity (Equations C-6, C-7). With ring particles, the code sums the force vectors and finds the net torque on the ring. Integration of this torque provides the motion along the ring's axes, which in turn modifies the Euler coordinates of the ring.
Does a rigid-ring model reproduce the dynamics of a full N-body model? May and James (1984) compare the two models, concluding that ring models evolve similarly to N-body models, with the exception of radial phase mixing. Unlike real galaxies and N-body models, ring models do not permit mass from one orbit to mix with adjacent orbits. Therefore, we expect that warps will endure longer in a ring model, because of the relative isolation of each ring. Rings will also display characteristics of gyroscopic motion. Once misaligned from the gravitational potential, rings should precess and nutate.
I use a galactic model designed by Kuijken and Dubinski (1995) in their paper regarding the construction of nearly self-consistent galaxy models. Their technique creates a particle distribution by randomly sampling the distribution function of the galaxy's components. This sort of model is particularly useful in the study of warps because it is initially close to equilibrium. In general, galactic models take several simulation timesteps to settle, producing an initial model with a particle distribution different from that desired. The properties of the galactic model are listed in Table 4. This is a fully self-consistent model, with a hybrid ring-disk system as noted above. This model generally reproduces the features of the Milky Way.
All calculations are performed with a dimensionless system of units in which the mass of the disk is 0.82 M, and the gravitational constant G is 1. In this system, the natural units for mass, velocity, and length are respectively: M = 4.5 x 1010 Mu, V = 220 km s-1, and Rd = 4 kpc. One unit of time is roughly 1.78 x 107 years, 17.8 million years. . Galaxy Model Parameters
| Disk | Bulge | Halo | |||||
| M | r,0 | Re/Rd | M | Re/Rd | M | Re/Rd | |
| Model | (1) | (2) | (3) | (4) | (5) | (6) | (7) |
| MW-B+ | 0.82 | 0.47 | 6.0 | 0.43 | 1.0 | 10.2 | 30.1 |
(1) disk mass, (2) disk central radial velocity dispersion, (3) disk radial extent, (4) bulge mass, (5) bulge radial extent, (6) halo mass, (7) halo radial extent
The bulge is realized as a King model (see §4.2.4 below), and represented with 10,000 particles. It extends to roughly 1 Rd with 1/2 of its mass within 0.4 Rd. Its total mass is 0.43 M, individual particles weigh 4.3 x 10-5 M.
The ring system is composed of 1000 concentric rings with a combined mass of 0.81744 M. The rings span radii from ~0 to 8 Rd; the rings from 6 to 8 Rd (rings 750-1000) are essentially massless, intended to represent gas. The inter-ring separation is 0.008 Rd, or 32 pc, constant across all but the first few rings. The mass per ring peaks with 0.00252 M at ~1 Rd, ring 127; the half mass radius is ~1.6 Rd, ring 200 (see Figure 7). The median ring mass is 9.2 x 10-4 M. Orbital time at M1/2 is 10.4 simulation times, or ~1.8 x 108 yr. For comparison to the Milky Way, the orbital time at the earth's radius is ~2.29 x 108 yr., roughly 230 million years. As seen in Figure 8, the tangential velocity of the disk is roughly constant.
The halo is represented by 100,000 particles. It extends to roughly 30 Rd, with 99% of the mass within 26 Rd (Figure 9). The total halo mass is 10.21 M, each particle weighs 1.02 x 10-4 M, or 4.59 x 106 Mu. The half mass radius is 6.86 Rd. The total halo mass is 4.59 x 1011 Mu within the cutoff of 120 kpc; this is slightly lighter than the 6 x 1011 Mu halo used by Weinberg (1995). As discussed in §2.2.3, the mass and extent of the halo may be underestimated in this model. The halo is only slightly flattened. In a system with a spherical halo there are no fundamental planes, and therefore no restoring force. As discussed in §3.2, a realignment force exists between a flattened halo and a tilted disk. The disk warps as the components align from the core outwards. In this model, we hope to observe a disk tilt.
Halo particles of order 106 Mu are effectively orbiting black holes (Lacey & J. Ostriker 1985). As they pass through the disk, their large mass heats the disk, increasing the scale height. I expect that halo particle bombardment will produce most of the random noise in the disk model. While black holes of this mass are theoretically capable of exciting a warp (Nelson & Tremaine 1995), I normalize the runs to remove the background noise and examine only the effect from the satellite.
The self-consistent satellite galaxy follows a spherical King density profile (King 1966). The distribution function of a King model is that of an isothermal model, but truncated at finite radius. King models allow alteration of the concentration, the variance of which produces brightness profiles that fit globular clusters and some elliptical galaxies well. The parameters that describe a King model are the King radius r0 , the tidal radius rt , and the concentration c. The King radius r0 is defined as
,
(4-1)
where 0 is the central density and the velocity dispersion parameter. The tidal radius rt specifies an approximate size of the satellite, beyond which material becomes unbound. Finally the concentration c is defined as
c log10(rt/r0). (4-2)
The satellite models are each composed of 10,000 particles. The two models weigh 0.108 M and 0.323 M; their properties are summarized in Table 5. These values approximate the predicted range for the LMC's mass as measured by Meatheringham et al. (1988) and Schommer at al (1992): 4.5 x 109 Mu and 1.44 x 1010 Mu respectively. . Satellite Model Properties
| Mass | c | Rtidal | R1/2 | |
| Run1 | 0.108 | 0.6 | 1.25 | 0.25 |
| Run2 | 0.323 | 0.6 | 2.16 | 0.44 |
(1) mass, (2) concentration, (3) tidal radius, (4) half-mass radius
To determine the initial position and velocity for the satellite models, I first take the present-day position and velocity measured by Jones et al. (1994). The inclination of the satellite orbit is -76°, retrograde with respect to the disk's rotation. To run the satellite 'backwards' in time, I reverse the present-day velocity and integrate a single particle orbit within the fixed axisymmetric galactic potential used to generate the halo model (see §4.2). I estimate the effects of dynamical friction using the halo density at the particle's position to apply Chandrasekhar's formula (see Appendix D) in the opposite sense.
The reverse orbit produces perigalacticon at 44 kpc about 168 million years ago, similar to the predictions of Gardiner et al. (1994). I run the reverse orbit back through two orbital times; the end-point of the reverse orbit becomes the starting point for the particle model. Finally, I recenter the satellite particle model to its initial position and velocity opposite that of the reverse orbit.
Although the reverse orbit starts at the present position of the LMC, the simulations do not identically reproduce the forward orbit. This occurs because the reverse orbit is calculated in a static potential, and is unable to estimate density variances produced during the simulation. These differences may be observed by comparing Figures 10 and 19.
To generate the initial conditions for Run1 and Run2, I merge the appropriate satellite model with the galactic model (Table 6). I then recenter the entire particle model such that the center of mass is at the origin. The simulation code then operates on the resulting binary file.
The initial position of the galactic model's combined elements are shown in Figure 10, including the estimated satellite orbit. While Figure 10 omits the halo for clarity, half of the halo mass resides within the radius of the disk, and over 90% within the area of the plot (see Figure 9).
Table 6. Particle Model Parameters
| N | Max R | R1/2 | M | |
| (1) | (2) | (3) | (4) | |
| Disk | 1,000 x 32 | 8 Rd | 1.6 Rd | 0.81744 M |
| Bulge | 10,000 | 1 Rd | 0.247 Rd | 0.43 M |
| Halo | 100,000 | 30 Rd | 6.86 Rd | 10.2 M |
(1) particles, (2) radial extent, (3) half-mass radius,
(4) total mass
[an error occurred while processing this directive]