ACCRETION OF THE MAGELLANIC SYSTEM ONTO THE GALAXY

The orbit of the LMC and SMC will be affected by their current positions, velocities, and by their masses. The positions and velocities used here for the LMC and SMC are given in Table 1. Particular note should be taken of most recent measurements of the proper motions of the LMC and SMC by Kallivayalil et al. (2006), which give a relative velocity between the clouds at the current epoch of 105 ± 42 km s⁻¹. These values imply that the LMC tangential velocity is approximately 100 km s⁻¹ higher than previously thought (Kallivayalil et al. 2006).

Given the low uncertainty in position, we assume that the LMC and SMC are at fixed locations given by the best estimates of observations. The velocity, with higher uncertainty, is randomly altered for each run assuming independent Gaussian uncertainties in the measurements, the line-of-sight velocity and the proper motions are calculated and then transformed to Cartesian velocities around the Galaxy (for transformations to the Cartesian coordinates see van der Marel et al. 2002)

The masses of the LMC and SMC have larger uncertainties than that of the position and velocity. van der Marel et al. (2002) obtained a mass for the LMC of (8.7 ± 4.3) × 10⁹ M_☉ within 8.9 kpc using an analysis of carbon stars. This mass is less than half that estimated by Schommer et al. (1992) who derive a mass of 2 × 10¹⁰ M_☉ using radial velocities for several of the oldest star clusters in the LMC that lie well beyond 6 kpc of its center. We match the circular velocities of van der Marel et al. (2002) with that of an Einasto halo at 8.9 kpc to get a virial mass of 2 × 10¹¹ M_☉ for a halo virialized today. Using a simplified model of tidal radius r³_tid/R_peri³ = m_tid/M_Galaxy we find that the LMC will have a mass of 5.6 × 10¹⁰ M_☉ within a tidal radius of 30 kpc, only slightly more massive than that of Schommer et al. (1992).

For the SMC, Hardy et al. (1989) obtained a lower mass limit of 1.0 × 10⁹ M_☉ from observations of carbon stars, while Dopita et al. (1985) obtained 0.9 × 10⁹ from similar observations of planetary nebulae. More recently, Harris & Zaritsky (2006) determined a mass of (2.7–5.1) × 10⁹ M_☉. We use the velocity dispersion from Harris & Zaritsky (2006) to calculate an Einasto halo virial mass of 4.5 × 10¹⁰ M_☉, with a mass 5.7 × 10⁹ M_☉ within a tidal radius of 19 kpc.

These masses, somewhat above the calculated masses at smaller radii, also approximate the mass of a theorized common halo between the LMC and SMC (see, e.g., Bekki 2008), without the problem of a large subhalo within a subhalo (in common halo models, the LMC halo will have a mass of ∼20% of the common halo). These masses are also similar to those calculated from the stellar mass content and seen inside N-body simulations (Boylan-Kolchin et al. 2011).

For the model, we ignore any mass loss from the clouds (in particular that due to tidal stripping), although mass loss is known to occur as evidenced by the amount of matter in the Magellanic Stream; (∼2 × 10⁸ M_☉; Putman et al. 2003). Dynamical friction between the LMC, SMC, and the Galaxy is also calculated using the tidal masses of the LMC and SMC. The effect of dynamical friction between the LMC and SMC is ignored, although LMC/SMC binary systems are still found in N-body simulations that include dynamical friction (Boylan-Kolchin et al. 2011) suggesting that the survival time with dynamical friction is sufficiently large for the Magellanic system to reach the present day with these masses.

The LMC/SMC and dwarf galaxy orbits occur within the potential of the Galaxy and are altered due to interactions between the Magellanic clouds themselves and the dwarfs.

There has been much work on Galaxy potential models since the first basic schemes devised by Murai & Fujimoto (1980) and Fujimoto & Sofue (1976). In contrast to other recent approaches that take simple spherical Galactic potentials, we choose to model the Galaxy using a multicomponent model of the Galactic potential developed by Flynn et al. (1996). This model matches known Galactic parameters such as the rotation curve, local disk density, and disk scale length to high accuracy. The lack of noticeable warp in the Magellanic Stream either in position on the sky or radial velocity—see Lin et al. (1995)—suggests that the disk and bulge likely have small or negligible effects on the orbits of the clouds, allowing us to have greater faith in our derived orbits. In this model the Galactic potential Φ(R, z) is given in cylindrical coordinates, where R is the planar Galactocentric radius and z is the distance above the plane of the disk. The total potential Φ is modeled by the sum of the three different potentials: the dark halo Φ_H, a central component Φ_C, and a disk Φ_D; thus,

$$\begin{equation} \Phi =\Phi _H+\Phi _C+\Phi _D. \end{equation} \tag{ 1 }$$

The potential of the dark halo Φ_H is assumed to be spherical and of the form

$$\begin{equation} \Phi _H=\frac{1}{2}V_H^2\ln \big(r^2+r_0^2\big), \end{equation} \tag{ 2 }$$

where r is the Galactocentric radius (r² = R² + z²). The potential has a core radius r₀, and V_H is the circular velocity at large r. The potential of the central component Φ_C is modeled by two spherical components, representing the bulge/stellar-halo and inner core components:

$$\begin{equation} \Phi _C=\frac{GM_{C_1}}{\sqrt{r^2+r_{C_1}^2}}-\frac{GM_{C_2}}{\sqrt{r^2+r_{C_2}^2}}, \end{equation} \tag{ 3 }$$

where G is the gravitational constant, $M_{C_1}$ and $r_{C_1}$ are the mass and core radius of the bulge/stellar-halo term, and $M_{C_2}$ and $r_{C_2}$ are the mass and core radius of the inner core. The disk potential Φ_D is modeled using an analytical form that is a combination of three Miyamoto–Nagai potentials (Miyamoto & Nagai 1975):

$$\begin{equation} \Phi _D=\Phi _{D_1}+\Phi _{D_2}+\Phi _{D_3}, \end{equation} \tag{ 4 }$$

where

$$\begin{equation} \Phi _D=\frac{-GM_{D_n}}{\sqrt{(R^2+[a_n+\sqrt{(z^2+b^2)}]^2)}}\quad n=1,2,3, \ldots \end{equation} \tag{ 5 }$$

and the parameter b is related to the disk scale height, a_n to the disk scale length, and $M_{D_n}$ the masses of the three disk components. Flynn et al. (1996) discuss in detail the justification of this model and its close fit to observational data. We assume hereafter that it reasonably accurately models the Galactic potential with the adopted parameters listed in Table 2.

Table 2. Adopted Parameters for the Galactic Potential

Component	Parameter	Value
Dark halo	r0	8.5 kpc
	VH	220 km s−1
Bulge/stellar-halo	$r_{C_1}$	2.7 kpc
	$M_{C_1}$	3.0 × 109 M☉
Central component	$r_{C_2}$	0.42 kpc
	$M_{C_2}$	1.6 × 1010 M☉
Disk	b	0.3 kpc
	$M_{D_1}$	6.6 × 1010 M☉
	a1	5.81 kpc
	$M_{D_2}$	−2.9 × 1010 M☉
	a2	17.43 kpc
	$M_{D_3}$	3.3 × 109 M☉
	a3	34.86 kpc

Download table as:  ASCIITypeset image

Solving Poisson's equation with the Galactic potential and integrating to 100 kpc gives a total mass of M_MW = 1.43 × 10¹² M_☉. This value agrees well with the Galactic mass derived in previous work for radii greater than 20 kpc (see Table 3).

Although the potential well of the Galaxy is by far the deepest, the LMC, SMC, and dwarfs will interact with each other. We assume that the LMC and SMC (and each dwarf) consist of an Einasto halo, with a sharp boundary at the tidal radius. The acceleration experienced by one of these objects at x₁ by another at x₂, separated by a distance r = ||x₁ − x₂|| is then (Nichols & Bland-Hawthorn 2009)

$$\begin{equation} \frac{{d}\Phi }{d{\bf x}} =\left\lbrace \begin{array}{@{}ll}3v_s^{2}r_s2^{-2-3/\alpha }\exp (2/\alpha)\alpha ^{-1+3/\alpha }\\ \times \quad \gamma (3/\alpha,2[r/r_s]^{\alpha }/\alpha)\frac{({\bf x}_1 - {\bf x}_2)}{r^3} &{\rm if \ } r < r_{{\rm tidal},2}\\ \\ \frac{GM_{2}({\bf x}_1 - {\bf x}_2)}{r^{3/2}} &{\rm if \ } r \ge {}r_{{\rm tidal},2}, \end{array}\right. \end{equation} \tag{ 6 }$$

where v_s is the halo scale velocity, r_s is the halo scale radius, and α is the Einasto scale factor (set here at α = 0.18).

The relatively narrow spread in dwarf galaxy positions suggests that if dwarf galaxies did fall in with the Magellanic system they would be more tightly bound than most dwarf associations seen today. We hence assume that the dwarfs were in a bound orbit of the most massive object in the Magellanic system, the LMC, at a previous apogalacticon of the LMC.

Each dwarf is randomly assigned a periapsis and circularity (and hence eccentricity) according to the distributions of Wetzel (2011) extended down to the lower mass of the LMC. These distributions lead to a cumulative distribution function (CDF) for periapsis of

$$\begin{eqnarray} F\left(\frac{r_{\rm peri}}{{\rm kpc}}\right) &=& 48.7\lbrace \gamma ^{-1}\big (3/17,\gamma _{\rm min}\nonumber \\ &&+[r_{\rm peri}/r_{\rm vir, LMC}]^{0.85}[\gamma _{\rm max} - \gamma _{\rm min}]\big)\rbrace ^{1/0.85}, \quad \end{eqnarray} \tag{ 7 }$$

where γ⁻¹ is the inverse lower incomplete gamma function, and

$$\begin{equation} \gamma _{\rm min/max} = \gamma (3/17,[0.32r_{\rm min/max}/r_{\rm vir, LMC}]^{0.85}), \end{equation} \tag{ 8 }$$

where r_min/max is the minimum and maximum allowed periapsis of the dwarfs. We choose an r_max of the tidal radius, as any dwarf that finds itself outside the tidal radius at closest approach will experience a greater attraction to the Milky Way than the LMC and will hence be definitely lost. A minimum periapsis is chosen to be three times the exponential scale length of the LMC disk (disk scale length ∼1.4 kpc; van der Marel et al. 2002), and orbits close to the disk are likely to produce warping that is not observed within the LMC disk.

The circularity CDF is given by

$$\begin{equation} F(\eta) = 2.25\eta ^{2.05} \phantom{.}_2F_1(2.05,-0.66;3.05;\eta), \end{equation} \tag{ 9 }$$

where η is the circularity of the orbit and $\phantom{.}_2F_1(a,b;c;z)$ is the Gaussian hypogeometric function.

Any orbits which would extend past the calculated tidal radius are excluded and a new one calculated.

The position of the dwarf within the orbit is determined by approximating the rosette orbit that arises from the spherical potential as a Keplerian orbit (the distributions also approximate it as Keplerian, but an extended Navarro–Frenk–White halo sees no significant change outside of low-periapsis orbits; see Wetzel 2011) and is found by solving

$$\begin{eqnarray} P\left(r < xr_{\rm peri}\right) &=& \frac{1}{\pi }\left(\frac{\pi }{2} - \beta - \arctan \left[\frac{1-x(1-\epsilon)}{\beta }\right]\right),\nonumber\\ \beta &=& \sqrt{(x-1)(1-\epsilon)(1+\epsilon -x[1-\epsilon ])}, \end{eqnarray} \tag{ 10 }$$

where $\epsilon =\sqrt{1-\eta ^2}$ is the eccentricity of the orbit.

By randomly determining the angular position of the orbit and the radial direction of travel, the position and velocity can then be determined in the Cartesian coordinates.

The mass of the dwarfs was set at either 1.1 × 10⁷ M_☉ within a tidal radius of 300 pc, consistent with the common mass of dwarf halos from Strigari et al. (2008), or a mass of 2.36 × 10⁸ M_☉ within the same tidal radius which represents a dwarf halo 10 times more massive at the virial radius. The tidal radius of the dwarfs was set arbitrarily to enclose this mass, although any dependence on mass is small (Section 3.1).

The amount of dwarfs will determine the amount of dwarf–dwarf interactions that occur, and within the comparatively small halo of the LMC these dwarfs could gravitationally disturb each other. The frequency of interactions between small dwarfs (with a mass M_dwarf/M_LMC ∼ 10⁻³) bound to the LMC will be a factor of their orbital distribution—if they themselves are in a disk-of-satellites around the LMC, not modeled here, interactions will be more frequent—and the amount of dwarfs present.

The frequency and effect of these dwarf–dwarf interactions are examined by comparing the final energies of the dwarfs around the LMC in models A and C and models B and D. We consider a dwarf to be unbound at a point in time if it has enough energy such that it would escape the LMC considering only a two-body interaction. A two-sample Kolmogorov–Smirnov (K-S) test could not rule out the null hypothesis that the distributions in final energies from the LMC of models A and C or models B and D are consistent today (with probabilities that the distributions are consistent of p_AC > 0.3, p_BD > 0.4) or at +500 Myr (p_AC > 0.3, p_BD > 0.9). This null result suggests that dwarf–dwarf interactions are rare with only a few dwarfs bound to the LMC, or that such interactions cause only slight perturbations in other bound dwarfs orbits.

The effect of increasing dwarf tidal mass by a factor of ∼20 is measured by the comparison of the energy of the dwarfs in models C and E. Applying a two-sample K-S test to the energy of the dwarfs from models C and E suggests the null hypothesis cannot be ruled out with probabilities of the same distribution being p_CE ∼ 0.03 today and p_CE > 0.2 at +500 Myr. The increase in p value at +500 Myr for all model comparisons may be a consequence of the strongest LMC–Galaxy–dwarf interactions occurring near perigalacticon, approximately today, with any variance along the LMC–SMC orbit having a bigger effect now than in +500 Myr.

Due to the large differences in the percentage of dwarf galaxies that become unbound between models that began at the last or second last apogalacticon or that had no apogalacticon, see Table 4, it would be expected that there are significant variations in the energy distributions of the dwarfs. Using a two-sample K-S test shows significant differences between the final energy of dwarfs in models A and B today with the probability of the null hypothesis, that is the final energies of dwarfs in models A and B coming from the same distribution, being p_AB < 10⁻³ with a K-S statistic of D = 0.70 and between models C and D today p_CD < 10⁻³ with a K-S statistic of D = 0.69. These energy distributions are also significant between models C and F (p_CF < 10⁻³, D = 0.20) and models D and F (p_DF < 10⁻³, D = 0.81). Similarly the differences are significant at +500 Myr (p < 10⁻³ between A and B, C and D, C and F, and D and F). Even within bound dwarfs, the distributions in energy are significantly different (p < 10⁻³; A and B, C and D, C and F, and D and F). As a run lost a higher percentage of dwarfs, the remaining bound dwarfs were of lower energy; that is, bound dwarfs in model F were of higher energy than dwarfs in model C which were of higher energy than model D. This loss of high-energy dwarfs occurs as the LMC/SMC approaches perigalacticon. Runs in model F, which had orbits most altered by dwarfs due to the large timescales, were more likely to lose dwarfs the closer the perigalacticon passage (occurring approximately today).

Due to the varying initial positions and velocities of the LMC and SMC, the orbits differ substantially between runs, with times for apogalacticons varying by hundreds of millions of years. Within model C, the SMC is bound to the LMC at the last apogalacticon in 13,840 (55%) runs, within model D the SMC is bound to the LMC two apogalacticons ago in 2691 (10%) runs, and within model F the SMC is bound to the LMC 2500 Myr ago in 559 (22%) runs. A bound orbit for model D is shown in Figure 2 from t ∼ −2500 Myr, with the current velocity for both the LMC and SMC shown as a vector emanating from their current position. An unbound orbit for model D is shown in Figure 3 from t ∼ −3000 Myr, with the initial velocity again as a vector. Due to the large number of orbits in which the SMC may be unbound, this orbit is not necessarily typical of all runs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If the SMC was unbound at the last apogalacticon while the dwarfs were bound, the dwarfs were more likely to be ejected (p < 10⁻³), possibly as a result of the SMC and LMC becoming bound. However, this effect is small, with on average 0.08 more dwarfs ejected per run in model C when the SMC was unbound at the last apogalacticon, with no impact if the SMC was bound at the second to last apogalacticon (model D, p > 0.05). This suggests that the SMC becoming bound to the LMC may only remove dwarfs that would in any case become unbound at a future point.

Despite being gravitationally bound at a previous apogalacticon, a large fraction of dwarfs become unbound from the LMC, often ending up a great distance away from the Magellanic system at the present day.

The final positions today (shown in Figure 4 for model C, Figure 5 for model D, and Figure 6 for model F) show clear differences due to the starting apogalacticon.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For dwarfs that were bound at the last apogalacticon (model C), ∼30% of these will be unbound today; however they remain relatively close to the LMC and are spread out along the orbital path of the Magellanic system. These unbound dwarfs form a bimodal distribution leading and trailing the LMC. This bimodality illustrates the difficulty for a dwarf to remain unbound near the LMC, whose gravity will recapture most dwarfs close by.

Of the dwarfs that were bound two apogalacticons ago (model D), a much larger fraction (∼90%) have become unbound by the present day and have begun their own orbits around the Galaxy. These unbound dwarfs occupy three main types of orbits, a relatively small (∼100 kpc) low-eccentricity orbit around the Galaxy, a much larger (∼300 kpc) highly eccentric orbit, with others remaining close to the LMC/SMC and traveling along similar orbits. As in model C, dwarfs that are close to the LMC and unbound have a bimodal distribution leading and trailing, again showing that is difficult for a dwarf to remain unbound near the LMC's gravitational influence.

Very few of the dwarfs that accompanied the LMC in an originally unbound orbit (model F) have become unbound in the model today or at +500 Myr. These dwarfs are much more likely to have a significant effect on the orbit of the LMC and SMC and correspondingly show a larger region of bound dwarfs as the model LMC's orbit is altered from the initial backward integration. The distribution of unbound dwarfs is similar, however, to those dwarfs bound at the LMC at the last apogalacticon with a leading group and a slight overdensity trailing the LMC.

Despite these dwarfs being ejected from the system, they remain in the same plane as the orbit of the LMC and SMC and within the disk-of-satellites seen around the Galaxy today; these unbound dwarfs also end up retaining their original direction of travel despite potentially being ejected at high velocity in a complicated three-body interaction. All but one dwarf from models C and D are moving clockwise in the Y − Z plane, with the one dwarf traveling counterclockwise having only a small velocity in the plane (v_YZ ∼ 8 km s⁻¹).

Of note is that there are no dwarfs in model C, D, or F which comes close to the current position of Leo I (−122 kpc, −117 kpc, 189 kpc) (Kroupa et al. 2005). In addition, with Leo I likely being unbound to the Galaxy, if it did infall bound to the Magellanic system an extremely rare multi-body reaction must have occurred to boost it away from the Galaxy. Similar ejections have been seen in large N-body simulations (Sales et al. 2007) between satellites of the host system. Leo I could therefore still have been associated with the Magellanic system, but unbound to the LMC when entering the Galactic halo.

The large fractions of simulated dwarfs that were bound in the past and unbound today mean that the lack of dwarfs observed to be orbiting the LMC today does not imply there were no dwarfs bound to the Magellanic system in the past.

With a large fraction of dwarfs becoming unbound, many of which leave the Magellanic system, the properties that allow dwarfs to become unbound can be examined: in particular, the contribution of the orbits of the LMC and SMC and the initial orbit of the dwarfs themselves around the LMC.

The effect of variations in the Magellanic system's orbit on dwarf loss is examined by looking at how many dwarfs were lost in each run. If the LMC/SMC orbits have significantly different loss rates, then future refinement of the orbit of the clouds may give evidence against dwarfs falling in with the clouds. If no physical property of the LMC/SMC orbits correlates significantly with dwarf ejection probabilities, then we can model ejections as a Poissonian distribution. Hence if ∼30% of all dwarfs were lost, as in model C, and with no contribution from the orbits of the LMC and SMC, the number of runs which lose 0–4 dwarfs would be a Poissonian distribution with ratios of (6748, 10530, 6162, 1602, 156). In fact, the runs from model C lose these dwarfs in the proportions of (6866, 10379, 6137, 1637, 181), a difference that is statistically insignificant (a χ²-test gives p > 0.05). The effect of the orbits is also not seen after the second to last apogalacticon in model D; again modeling as a Poisson process the number of dwarfs lost per run would be expected to be in the proportions (6, 181, 1895, 8820, 15397), while the results from model D give (6, 225, 1891, 8697, 15481) (a χ²-test gives p > 0.01). This suggests that the positions of the dwarfs at apogalacticon are much more important than the orbits of the LMC and SMC (excluding whether the SMC is bound initially, see Section 4.1).

The impact of different orbits of dwarfs around the LMC at apogalacticon is also examined. Dwarfs which originally possess a low energy—that is, are deep within the potential well of the LMC—are less likely to be stripped either when the SMC passes close by or at perigalacticon, when the gravitational attraction of the Galaxy is largest. A K-S test comparing the initial energy of the dwarf orbits that remain bound and those that are unbound indicates that the distributions of initial energy are different for model C today and at +500 Myr (the null hypothesis that the distributions are the same has probabilities of p_today < 10⁻³ and p_+500 Myr < 10⁻³) and for model D today and at +500 Myr (p_today < 10⁻³, p_+500 Myr < 10⁻³). The fraction of dwarfs that remain bound as a function of initial energy is shown in Figure 7. Dwarfs that were bound one apogalacticon ago (model C) are still undergoing removal from the Magellanic system as it passes close to the Galactic center. Dwarfs that were bound two apogalacticons ago (model D) are no longer undergoing removal, with dwarfs that are bound now remaining bound over the next 500 Myr, with a handful being recaptured as they pass close by the LMC during their own orbit of the Galaxy. No matter when the dwarfs were originally bound, dwarfs that were deep inside the potential well originally are much more likely to be bound today and remain bound over the next 500 Myr, while dwarfs of energy above about −1.8 × 10¹⁴ erg are all equally likely to be removed.

Zoom In Zoom Out Reset image size

The large tail of dwarfs with highly eccentric orbits seen in model D is also reflected in the energy–angular-momentum diagram (Figure 8), with these dwarfs occupying high angular momentum, high energy positions on this diagram. Relatively few unbound dwarfs occupy the same position on this diagram as dwarfs that end up bound, and any detected dwarfs that occupy this position will be less likely to have fallen in with the LMC on a previous orbit.

Zoom In Zoom Out Reset image size