ROTATION OF THE MILKY WAY AND THE FORMATION OF THE MAGELLANIC STREAM

The number of parameters required to describe the interaction involving the Milky Way, the LMC, and the SMC depends on the physical model adopted. We have employed a restricted N-body numerical code assuming the gravitational interaction only. Even such a simplified view of the interaction involves ≈25 independent parameters including the initial conditions of the LMC and the SMC motion, their total masses, parameters of mass distributions, particle disk radii, and orientation angles, and also the parameters defining the gravitational potential of the Milky Way, i.e., the flattening parameter of the dark matter (DM) halo, the LSR circular velocity, and the LSR angular rotation rate (see Section 5). Some of the parameters were constrained by theoretical studies (including scale radii of the LMC/SMC halos, the Coulomb logarithm Λ for the dynamical friction in the Milky Way halo, and the halo flattening parameter q for the model of the gravitational potential of the Galaxy). Their mean values and searched errors were discussed in Růžička et al. (2007). This section introduces the observationally estimated parameters of the LMC–SMC–Galaxy interaction (see also Table 1).

Table 1. Free Parameters of the LMC–SMC–Milky Way Interaction

Parameter	Value	Comment
(αlmc, δlmc)	(8190 ± 098, − 6987 ± 041)	Equatorial coordinates
(αsmc, δsmc)	(132 ± 03, − 725 ± 03)
(m − M)lmc	18.5 ± 0.1	Distance moduli
(m − M)smc	18.85 ± 0.10
μlmcW(mas yr-1)	〈 − 2.11, − 1.92〉	Proper motion components
μlmcN(mas yr-1)	〈 + 0.39, + 0.49〉
μsmcW(mas yr-1)	〈 − 1.34, − 0.69〉
μsmcN(mas yr-1)	〈 − 1.35, − 0.99〉
υlmcrad(km s-1)	262.2 ± 3.4	Line-of-sight velocities
υsmcrad(km s-1)	146.0 ± 0.60
rlmcdisk(kpc)	〈9.0, 12.0〉	Initial radii of the particle disks
rsmcdisk(kpc)	〈5.0, 8.0〉
(ilmc, plmc)	(347 ± 62, 1299 ± 60)	Inclination and position angles of the particle disks
(ismc, psmc)	(60° ± 20°, p = 45° ± 20°)
Θ0(km s-1)	〈210, 260〉	LSR circular velocity
q	〈0.71, 1.30〉	Flattening of the DM halo of the Milky Way

Download table as:  ASCIITypeset image

It has been shown by Růžička et al. (2009) that the current proper motions of the Magellanic Clouds are the most critical parameters of their interaction with the Milky Way. Figure 1 shows the proper motions of the LMC and the SMC as estimated by the latest analysis by Kallivayalil et al. (2006a, 2006b) and Piatek et al. (2008) based on the measurements by the Hubble Space Telescope (HST). For convenience, the LMC/SMC proper motion vectors are decomposed into two components related to the local directions of west (μ_W) and north (μ_N). The reason for involving two results derived from a single data set is the discrepancy between the studies by Kallivayalil et al. (2006a) and Piatek et al. (2008) regarding the SMC western proper motion μ_W. The resulting proper motion space for the Magellanic Clouds explored in this study is a union of the sets defined by the 1σ errors of the data by Kallivayalil et al. (2006a, 2006b) and Piatek et al. (2008) for the LMC and the SMC, respectively (Figure 1):

Detailed analysis by Růžička et al. (2009) has quantified the sensitivity of the dynamical evolution of the Magellanic System to the variation in different free parameters. Except for the LMC/SMC current spatial velocities and positions, the remaining parameters had a rather weak impact on the interaction. As long as the parameters are independent adding other parameters to the original set does not change the conclusion by Růžička et al. (2009) regarding their importance. However, once parameter independence cannot be guaranteed, it is a good practice to keep the dimensionality of the parameter space as high as possible. Therefore, we will introduce the LSR circular velocity as a new dimension of the Magellanic parameter space established in a way similar to Růžička et al. (2007, 2009).

Unlike the proper motion of the Clouds, their LSR radial velocities could be measured with high accuracy. Following van der Marel et al. (2002), we set υ^lmc_rad = 262.2 ±  3.4 km s⁻¹. The SMC radial velocity error was estimated by Harris & Zaritsky (2006) as υ^smc_rad = 146.0 ±  0.60 km s⁻¹.

The heliocentric position vector of the LMC was adopted from van der Marel et al. (2002), i.e., the equatorial coordinates are (α_lmc, δ_lmc) = (8190 ±  098, − 6987 ±  041), its distance modulus is (m − M)_lmc = 18.5 ±  0.1. The equatorial coordinates of the SMC were set to the ranges (α_smc, δ_smc) = (132 ±  03, − 725 ±  03) (see Stanimirović et al. 2004, and references therein). A range of distance determinations for the SMC was provided by van den Bergh (2000) and we used his resulting distance modulus (m − M)_smc = 18.85 ±  0.10.

Several observational determinations of the LMC disk plane orientation have been published so far (see, e.g., Lin et al. 1995; van der Marel et al. 2002). In our parameter study, the LMC inclination i and position angle p and their errors agree with van der Marel et al. (2002), i.e., i = 347 ±  62 and p = 1299 ±  60. As the SMC misses a well-defined disk, the orientation and the position angle usually refer to the SMC "bar" defined by Gardiner & Noguchi (1996). Based on the estimates by van den Bergh (2000) or Stanimirović et al. (2004), we adopted the error ranges i = 60° ±  20° and p = 45° ±  20° for the SMC initial disk inclination and position angle, respectively.

Gardiner et al. (1994) analyzed the H i surface contour map of the Clouds to estimate the initial LMC and SMC disk radii r^lmc_disk and r^smc_disk, respectively. With the use of their work and of the study by Brüns et al. (2005), we varied the LMC/SMC disk radii within the ranges of 10.5 ±  1.5 kpc (LMC) and 6.5 ±  1.5 kpc (see Table 1). Regarding the absence of a clearly defined disk of the SMC and possible significant mass redistribution in the Clouds during their evolution, the results require careful treatment.

Current total masses m_lmc and m_smc follow the estimates by van den Bergh (2000). The masses of the Clouds are functions of time and evolve due to the LMC–SMC exchange of matter and as a consequence of the interaction between the Clouds and the MW. Our test-particle model does not allow for a reasonable treatment of a time-dependent mass loss. Therefore, the masses of the Clouds are considered constant in time and their initial values at the starting epoch of simulations are approximated by the current LMC and SMC masses.

The gravitational potential of the Milky Way is modeled by the superposition of three static components including the axially symmetric logarithmic potential (Binney & Tremaine 1987) of the DM halo, the Miyamoto–Nagai potential of the Galactic disk, and the Hernquist bulge. Since the logarithmic halo is the dominant component of our model, its flattening parameter was treated as a free parameter. The flattening q was varied within the range of 〈0.71, 1.30〉. The lower limit is based on the condition required to avoid negative mass densities in the logarithmic halo. The upper limit was set up to introduce a convenient symmetry with respect to the spherical shape (q = 1). The scaling velocity factor υ₀ of the logarithmic model (more details in Binney & Tremaine 1987) is a function of the LSR circular velocity and the solar galactocentric distance.

The freedom in the values of the parameters q and υ₀ introduced a spread in the total mass of the Galaxy m_MW = 〈0.59, 5.90〉 × 10¹²M_☉ within the radius of 250 kpc. The estimates of the total mass of the Milky Way exceeding ≈3 × 10¹²M_☉ have neither observational not theoretical support (see, e.g., Binney & Tremaine 1987; Li & White 2008), but this fact will be taken into account when the results are interpreted.

Unlike the studies by Růžička et al. (2007, 2009), we have treated the LSR circular velocity Θ₀ as a free parameter. This decision was motivated by recent findings by Reid & Brunthaler (2004). Their measurements of the proper motion of Sag A* yielded the estimate for the LSR circular velocity Θ₀ = 236 ±  15 km s⁻¹. At the same time, they give a revisited value of the LSR angular rotation rate Ω₀ of 29.45 km s⁻¹ kpc⁻¹. These two quantities then define the variations of the solar galactocentric distance, as R₀ = Θ₀/Ω₀.

Reid & Brunthaler (2004) performed very precise position measurements of the famous radio source Sag A* with respect to two background extragalactic radio sources. The data collected over a period of eight years yielded values of the apparent motion of Sag A*. This quantity is a composition of the reflected solar galactocentric motion and of the peculiar motion of Sag A* which is quite small. The fraction of the apparent motion corresponding to solar rotation around the Galactic Center may be treated as a combination of the LSR circular velocity Θ₀ and the deviation of the solar motion from the circular orbit (Dehnen & Binney 1998). Finally, Reid & Brunthaler (2004) were able to determine the latter which yielded an estimate for the LSR circular velocity Θ₀ = 236 ±  15 km s⁻¹ assuming the solar distance to the Galactic center (R₀) of 8.0 ±  0.5 kpc. Notably, the uncertainty of Θ₀ originating from the measurement error was only 1 km s⁻¹. The rest appears due to the uncertainty of R₀ (Reid & Brunthaler 2004). Calculating the LSR angular rotation rate Ω₀ ≡ Θ₀/R₀ yields a value with a very small uncertainty: Ω₀ = 29.45 ±  0.15 km s⁻¹ kpc⁻¹. The corresponding relative error is only 0.15/29.45 = 5 × 10⁻³ which is by 1 order of magnitude below the relative error of any of the parameters examined in our study. Therefore, the LSR angular rotation rate Ω₀ will be treated as a constant in this paper, i.e., Ω₀ = 29.45 km s⁻¹ kpc⁻¹.

Although we are already familiar with the impact of most of the discussed parameters on the behavior of the restricted N-body model of the LMC–SMC–Milky Way interaction (see Růžička et al. 2009), this does not apply on the LSR circular velocity. Even though it is a natural choice, taking the 1σ error box of Θ₀ by Reid & Brunthaler (2004) for the range explored by the parameter study is arguably not a good practice. We want to avoid the case of important models being localized at either end of the studied parameter range. Such a coincidence would then require an extension of the studied range of the LSR circular velocity followed by another time-consuming run of the automated parameter search. The described procedure is necessary to uncover and understand the behavior of the interacting system in such a parameter region of apparently special features. Therefore, the range of the new input parameter Θ₀ was obtained by extending the 1σ uncertainty of the LSR circular velocity as published in Reid & Brunthaler (2004) by 10 km s⁻¹, i.e., Θ₀ = 〈210, 260〉 km s⁻¹, corresponding to the range of the solar galactocentric distance R₀ of 〈7.13, 8.83〉 kpc.

The model itself is an advanced version of the scheme by Gardiner et al. (1994): it is a restricted N-body (i.e., test particle) code describing the gravitational interaction between the Galaxy and its dwarf companions. The potential of the Milky Way is dominated by the flattened DM halo and dynamical friction is exerted on the Magellanic Clouds as they move through the halo (Binney 1977). The LMC and the SMC are represented by Plummer spheres, initially surrounded by test-particle disks. For further details, see Růžička et al. (2007, 2009).

The exploration of the parameter space for the interaction involving the Galaxy and the Magellanic Clouds was performed by a genetic algorithm. Genetic algorithms belong to the class of evolutionary optimizers that mimic the selection strategy of natural evolution. Holland (1975) was the first one who proposed the application of genetic algorithms on optimization problems in mathematics. Recently, the performance of genetic algorithms was studied for galaxies in interaction (see, e.g., Theis 1999). Theis & Kohle (2001) analyzed the parameter space of two observed interacting galaxies—NGC 4449 and DDO 125. Genetic algorithms turned out to be very robust tools for such a task if the routine comparing the observational and modeled data is appropriately defined. The approach by Theis & Kohle (2001) was later adopted and improved in order to explore the interaction of the Magellanic Clouds and the Galaxy (Růžička et al. 2007, 2009).

The comparison between the model and observations became more efficient by involving an explicit search for the shapes in the data (Růžička et al. 2009). Also, the significant system-specific features (such as a special geometry and kinematics) were taken into account, further improving the performance of the genetic algorithm for exploration of the LMC–SMC–Galaxy interaction. More detailed information is to be found in Section 4.3 of this paper.

The automated search of the parameter space is driven by a routine comparing the modeled and observed distributions of H i associated with the Magellanic System (Brüns et al. 2005). The match is measured by the fitness function (f) which is a function of all input parameters, as every parameter set determines the resulting simulated H i data cube. Our fitness function returns a floating-point number between 0.0 (complete disagreement) and 1.0 (perfect match) and consists of four different comparisons, including search for structures and analysis of local kinematics. For details on the fitness function used for this study, see the Appendix.

The fitness function f is the only part of the genetic algorithm that reflects the nature of the studied problem. Although the automated evolutionary optimization is a robust method applicable on a remarkable variety of systems, its actual performance and efficiency are critically dependent on the proper choice for the fitness function. If the function is defined in a sensible way, the number of parameter combinations examined by the genetic algorithm is minimized. However, the performance of the genetic algorithm also depends on its convergence rate expressed as , where i = 1, 2, ... stands for the number of the actual generation and is the mean value of the fitness function in the ith generation. The convergence rate is generally proportional to the ratio of the generation size N_gen and of the dimension n of the studied parameter space over which the fitness function is defined. Unfortunately, the value of N_gen/n is very low if genetic algorithm-based optimizer is applied on the parameter spaces of interacting galaxies. The number of parameters involved is always high (see Section 3). On the other hand, the maximum generation size has to be limited in order to keep the total computational time requested by the parameter search reasonable.

In order to overcome such a difficulty, we searched the parameter space repeatedly in a fixed number of optimization steps, i.e., generations of models, and localized ≈100 high-fidelity models. Such a procedure is not likely to reveal the global maximum of the fitness function, but it results in a map of distribution of good models over the entire parameter space. In principle, every region of the parameter space allowing for the satisfactory reproduction of the large-scale structures associated with the Magellanic Clouds may be identified.

The reader might ask why there is such emphasis placed on the fitness function itself if it, in fact, does not seem to provide any physical information about the interacting system of the Galaxy and the Magellanic Clouds. Indeed, the function f serves primarily as a driver to the genetic algorithm. However, the search for good models of the observed Magellanic System is efficient only if relevant astrophysical data are supplied as the input to f. As already mentioned, our study deals with detailed morphological and kinematic information from the 21 cm survey by Brüns et al. (2005) and with the corresponding modeled data. The fitness function then makes a link between the observable data and the initial state of the Magellanic System.

The rotation curve of a galaxy

introduces a fundamental constraint on the models of its overall potential Φ_tot(R, z). In our study, the following three component model of the potential of the Milky Way has been used:

where

is the logarithmic model of the Galactic DM halo,

is the Miyamoto–Nagai potential of the Galactic disk with the total mass M_d, and

is the Hernquist formula for the spherically symmetric potential of the Milky Way bulge with the mass M_b.

At the galactocentric distance where the Magellanic Clouds have resided, i.e., D ⪆ 50 kpc, the gravitational field of the Galaxy is dominated by its DM halo component. The role of the DM halo becomes more prominent with the increasing distance to the galactic center, because ∇Φ_L ∼ 1/r while ∇Φ_MN ∼ 1/r² and also ∇Φ_B ∼ 1/r². The necessary adjustment of the Milky Way potential Φ_tot required by Equation (2), if the LSR circular velocity Θ₀ is changed, was achieved by varying the velocity constant υ₀ in Equation (4).

The choice for the remaining parameters in Equations (4)–(6) follows Fellhauer et al. (2006) and Růžička et al. (2007, 2009) and so we have set a = 0.7 kpc, b = 6.5 kpc, c = 0.26 kpc, M_b = 3.4 10¹⁰M_☉, M_d = 10¹¹M_☉, and R_c = 12.0 kpc.

The value of the characteristic radius R_c determining the concentration of the Galactic halo is five times smaller than the present galactocentric distances of the Magellanic Clouds, and by the order of 10 smaller compared to the typical LMC/SMC galactocentric distances over the last several Gyr. Hence, its impact on the density profile of the region of the DM halo of the Milky Way, where the Clouds have resided, is weak. Regarding the total mass of the Galaxy enclosed within a given radius, it is determined by the velocity parameter υ₀, as

where the density of the logarithmic halo of the Galaxy ρ_L(R, z) ∼ υ²₀, while ρ_L(R, z) ∼ 1/R²_c. We did not consider the scale radius R_c as a free parameter because its influence on the mass distribution and the total mass of the DM halo of the Milky Way is significantly lower compared to the velocity parameter υ₀.

The flattening q of the DM halo potential entered our simulations of the Magellanic System as a free parameter (see Section 3). The velocity parameter υ₀ of the logarithmic halo was considered a function of the LSR circular velocity Θ₀. The corresponding formula comes out of the Equations (2) and (3) after substituting for Φ_L, Φ_MN, and Φ_B:

where R₀ = Θ₀/Ω₀.

The analysis of the Equation (8) shows that if

and

then

It can be easily seen that the conditions (9) and (10) are always satisfied for the values of the parameters a, b, c, and R_c assumed in our study. As we have Equation (7) for the logarithmic halo (Binney & Tremaine 1987), the total mass of the Milky Way enclosed within an arbitrary radius r must grow if the LSR circular velocity Θ₀ is increased.

A convenient transformation of the position and velocity vectors of the Magellanic Clouds from the heliocentric spherical coordinates to the galactocentric Cartesian frame (Figure 2) was derived by van der Marel et al. (2002). Let (l, b) and (α, δ) be the corresponding pair of the galactic and equatorial coordinates of the center of mass of either of the Magellanic Clouds, respectively. If the Cloud's present heliocentric distance is D, its galactocentric Cartesian coordinates are

where R_☉ is the solar galactocentric position vector and u₀ is the heliocentric unit position vector of the Cloud.

The heliocentric spatial velocity of an object in space is usually expressed in terms of its proper motion and the line-of-sight systemic velocity υ_sys. The proper motions in the directions of west and north are defined as (van der Marel et al. 2002)

The transformation from the heliocentric velocity coordinates μ_W, μ_N, and υ_sys to the galactocentric velocity components may be then expressed as

where

and υ_☉ is the solar galactocentric velocity vector.

Taking Equations (11) and (13) into account together with the fact that R_☉ = (−Θ₀/Ω₀, 0, 0) and υ_☉ = (0, Θ₀, 0) (see Figure 2), it is obvious that the current galactocentric positions and velocities of the Magellanic Clouds depend on the values of the LSR circular velocity Θ₀ and of the LSR angular rotation rate Ω₀. If either of these two parameters is varied, a change to the galactocentric position vectors of the Clouds occurs due to the change of their x-components:

Similarly, the galactocentric spatial velocities of the LMC and the SMC would be influenced, because their y-components involve the explicit dependence on Θ₀:

Note that the phase-space configuration of the system involving the Sun, the Magellanic Clouds, and the Galactic Center is such that υ¹ ∼ Θ₀, while |υ¹| ∼ −Θ₀. As a consequence, the magnitude of the galactocentric velocity

of the LMC decreases if the LSR circular velocity is increased. Figure 3 reveals clearly that the same conclusion holds for the SMC.

Zoom In Zoom Out Reset image size

Variation of the LSR circular velocity Θ₀ and of the LSR angular rotation rate Ω₀ modifies both the magnitudes and the directions of the position and velocity vectors of the Clouds. Specifically, increase in the LSR circular velocity causes the reduction of the magnitude of the LMC/SMC galactocentric velocities.

With the use of the genetic algorithm, ∼10⁶ parameter combinations, i.e., individual restricted N-body simulations, were examined in total, and ∼100 sets providing the highest fitness models were collected. In this section, the features of these models will be discussed with respect to the parameters influencing the motion of the Magellanic Clouds.

We have treated the LSR circular velocity as a free parameter in order to analyze its impact on the tidal interaction in the Magellanic System. The automated exploration of the modified parameter space revealed a significant qualitative change regarding the preferred proper motions of the Clouds, compared to the work by Růžička et al. (2009).

They have done a detailed analysis of how the models of the LMC–SMC–Galaxy interaction depend on different parameters. It was clearly shown that the interacting system is very sensitive to the parameters defining the past positions of the Clouds' orbits in phase space. In the work by Růžička et al. (2009), the set of the critical parameters involved the current heliocentric positions and velocities of the Magellanic Clouds and, to some extent, the flattening parameter of the Milky Way halo. However, introducing the new value of the LSR angular rotation rate by Reid & Brunthaler (2004) and treating the LSR circular velocity as a free parameter has changed the old picture.

Figure 4 shows the distribution of the local peaks of the fitness function f (i.e., the fitness of every model identified by the genetic algorithm) over the studied ranges of the western and northern components of the proper motions for both Magellanic Clouds. In contrast to Růžička et al. (2009), there is no apparent preference for specific values of the proper motions for either of the Clouds. Models of a very similar quality have been located over the entire studied proper motion ranges of the LMC and the SMC. Notably, the high-fitness models fall into either of two groups clearly separated by the fitness values of their members. The majority of the models have the fitness exceeding the value of 0.55 while most of the rest are described by the relation f < 0.50. Only ∼5% of all the models are located within the fitness range of 0.50 < f < 0.55. Later in Section 6.7, the lower limit for the fitness of the satisfactory models will be precisely established and discussed. For the moment, we will just mention that the almost empty belt between the fitness values of 0.50 and 0.55 separates the acceptable models of the Magellanic System from those reproducing the H i observational data by Brüns et al. (2005) insufficiently.

Zoom In Zoom Out Reset image size

It can be easily seen from Equations (4) and (8) that the LSR circular velocity Θ₀ and the flattening q of the Galactic DM halo determine the total mass of the Milky Way. We have focused on the distribution of all the high-fidelity models of the LMC–SMC–Galaxy interaction over the q–Θ₀ plane of the studied parameter space. Figure 5 visualizes the result. The actual value of the total enclosed mass of the Milky Way does not introduce any limitation regarding the ability of our tidal model to reproduce the observed H i large-scale Magellanic structures. Such behavior exists due to the fact that local peaks of the fitness function (i.e., good models) of a similar quality have been localized over the entire range of the LSR circular velocity.

Zoom In Zoom Out Reset image size

The distribution of the genetic algorithm fits with respect to the flattening parameter q is quite different. A strong preference for oblate (q < 1) configurations of the gravitational potential of the Milky Way halo exists. Such a result agrees with the conclusions by Růžička et al. (2007). In the next paragraphs, the role of the LSR circular velocity and the flattening of the Galactic halo will be addressed regarding the impact of these parameters on the orbital history of the Magellanic Clouds.

We have revealed an interesting anti-correlation in the output of the automated parameter search by the genetic algorithm. Figure 6 depicts the relation between the magnitude of the current galactocentric velocity of the LMC and the total mass of the Milky Way for all the high-fitness models. The general trend of the dependence, i.e., the galactocentric velocity of the LMC decreasing as the Galactic mass increases, is driven by the LSR circular velocity Θ₀ due to the way it links the Equations (7), (8), and (17), and together with the fact that the mass of the Milky Way is dominated by the DM halo. Equation (7) may be expressed as

and for the galactocentric velocity we obtain

where i = lmc. Equations (18) and (19) define the parametric representation of a curve in the m_MW–|υ_lmc| plane. If all parameters except the LSR circular velocity Θ₀ are assumed constant and Θ₀ is varied within the range of 〈210, 260〉 km s⁻¹, the parametric equations yield the curves depicted in Figure 6. The curves are associated with randomly selected high-fitness models by taking their parameter values for the input to Equations (18) and (19). Apparently, the curves are close to linear. Their slope is controlled by the LSR circular velocity and they are well aligned with the distribution of models in the m_MW–|υ_lmc| plane which provides the key to Figure 6. The distribution of models in the plot is primarily controlled by the LSR circular velocity which sets the linearly decreasing trend of the velocity–mass dependence. The vertical (velocity) and horizontal (mass) placement of the curves depends on the heliocentric positions and velocities of the LMC and on the flattening of the DM halo of the Milky Way, respectively.

Zoom In Zoom Out Reset image size

Figure 7 will help us to understand how the variations of the flattening parameter q and of the LSR circular velocity are reflected in the past orbits of the Clouds. One of the parameters was varied within its entire range considered in this study, while the second one was fixed to the value corresponding to a high-fitness model of the interaction.

Zoom In Zoom Out Reset image size

Figure 7 reveals a qualitative difference in the way the specified parameters affect the motion of the Magellanic Clouds. The influence of the LSR circular velocity is significantly stronger than the role of the flattening of the Galactic potential. This is not surprising. While the flattening parameter only has an impact on the distribution of mass of the Galaxy, variations in the LSR circular velocity also lead to changes of the LMC/SMC galactocentric velocities. Moreover, these effects amplify each other regarding their impact on the orbits of the Clouds (see Section 5). If Θ₀ is increased, the total mass of the Milky Way increases implying a higher gravitational potential energy of the Clouds with respect to the Galaxy. Regarding Equation (17) and the related discussion, the increasing LSR circular velocity reduces the galactocentric velocities of the Clouds, i.e., their kinetic energy. Thus, the total energy of the LMC–Galaxy and the SMC–Galaxy pairs changes with Θ₀ rapidly.

The lower row of Figure 7 shows that the orbits of the Magellanic Clouds projected to the plane of sky are always offset from the position of the Magellanic Stream, regardless of the choice for the [q, Θ₀] pair. A closer look at Figure 15 reveals that such a deviation is controlled by the northern component μ^lmc_N of the LMC proper motion and cannot be removed unless the results by Kallivayalil et al. (2006a, 2006b) and Piatek et al. (2008) are revisited. The particles forming the leading and trailing streams in our restricted N-body model (and generally in every tidal model) essentially follow the orbits of their progenitors. This purely tidal approach is thus limited. As long as the LMC/SMC projected orbits are offset from the observed Magellanic Stream, it is difficult to model the full extent of the large-scale structure of the Stream in the position–position–LSR radial velocity space (see Section 6.8).

Choi et al. (2007) have demonstrated that tidal tails originating in massive satellite galaxies become offset from their progenitors' orbits to some level. However, unless additional physical processes introduce the dissipation of energy and angular momentum of the particles in the tidal tails, such a phase-space deviation between the stream and the orbit of its progenitor is limited. Figure 8 shows how this is reflected by the kinematics of the Magellanic Clouds and the Stream. The high-fitness models are split into the groups of f > 0.60 and f < 0.50, respectively. It is a common feature of the best models (f > 0.60) that the LSR radial velocity along the LMC orbital track agrees well with the famous linear velocity profile of the Magellanic Stream (Brüns et al. 2005).

Zoom In Zoom Out Reset image size

Despite the previously discussed behavior of tidal models, the apparent offset of the recent (−300 Myr ≲ t ⩽ 0) LMC/SMC orbits from the Magellanic Stream was reproduced partially by our simulations. It did not occur due to a large deviation between the tidal particle tail and the orbits of the Clouds. The modeled projection of the Magellanic Stream to the plane of sky is still well aligned with the previous revolutions of the LMC/SMC orbits (see Figure 15), but the preceding orbital cycles are projected to different positions in the plane of sky reflecting the geometry and physical properties of the LMC–SMC–Milky Way system.

The subsequent revolutions of the orbits of the Magellanic Clouds are shifted with respect to each other in our good models of the interaction with the Milky Way. This is to be attributed to the projection effect of the configuration of the system involving the observer (Sun), the Galactic Center, and the Magellanic Clouds (Figure 2). The flattening of the Galactic potential must be taken into account, as it prevents the LMC from having an orbit confined within a two-dimensional plane. Figure 9 shows the intensity of the radial and the axial components of the gravitational field of the logarithmic halo. The ratio of the components depends strongly on the position [R, z] with respect to the Galactic Center. That leads to the formation of a three-dimensional orbit of the LMC. We did not mention the SMC, because its low mass makes it dependent not only on the gravitational field of the Galaxy, but also on the attraction by the LMC.

Zoom In Zoom Out Reset image size

The previous studies of the dynamical evolution of the Magellanic Clouds assumed a good alignment of the LCM/SMC orbits with the Magellanic Stream in the position–position–LSR radial velocity space (see, e.g., Gardiner & Noguchi 1996; Mastropietro et al. 2005; Connors et al. 2006). However, such an assumption is at odds with the HST proper motion measurements for the Clouds (Besla et al. 2007). We have seen that nothing changes even for the revisited values of Θ₀ and Ω₀ by Reid & Brunthaler (2004). What it represents for modeling the large-scale features of the Magellanic Clouds was briefly discussed in this section and will be further addressed later in Sections 6.8 and 6.9.

The Magellanic Clouds share a common low-density gaseous envelope (see Brüns et al. 2005) which is usually considered one of the signs indicating that the Clouds have formed a gravitationally bound couple even for several Gyr. Gravitational binding of the Clouds lasting several Gyr turned out to be extremely rare in our simulations. However, the proper motion data by Kallivayalil et al. (2006a, 2006b) and by Piatek et al. (2008) imply that the distance between the LMC and the SMC is very likely to have been shorter than today for several hundred Myr. This may be sufficient to explain the presence of the H i envelope, as our simulations have shown that an envelope surrounding the Magellanic Clouds may be formed rather quickly on the specified timescale.

The information about the LMC–SMC distance in the recent past comes directly from the combination of Figures 10 and 11. We have calculated the scalar product of the present galactocentric velocity vectors of the Magellanic Clouds for all the high-fitness models. Figure 10 shows the fitness of these models as a function of the actual value of the scalar product. The current LMC and SMC galactocentric velocity vectors are nearly parallel to each other in all cases and their angular deviation does not exceed 20° with the mean value of only ≈10°. It is natural to ask whether the spatial separation of the Clouds is increasing or decreasing at present. Figure 11 offers the answer. If the time dependence of the distance between the Clouds is plotted over the last 700 Myr for the models of f > 0.55, one may clearly see that the LMC and the SMC have been receding from one another for ≈100 Myr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To quantify the rate at which the LMC–SMC separation D_l-s changes we have calculated the function dD_l-s/dt at the present time t = t₀ = 0 Gyr. If the function dD_l-s/dt is evaluated for the models plotted in Figure 11, it returns the values between 50 kpc Gyr⁻¹ and 120 kpc Gyr⁻¹.

The preceding paragraphs lead to the conclusion that the orbital history of the Magellanic Clouds in the high-fitness models always involves a close encounter between the LMC and the SMC at the time of −250 Myr ≲ t ≲ −80 Myr. It is quite interesting to see whether such behavior is outstanding or typical within the LMC/SMC proper motion ranges established by the recent HST measurements.

Figure 12 illustrates how the rate at which the current spatial separation of the Magellanic Clouds changes depends on the proper motion of the SMC. We have examined two cases. In one case, the proper motion components of the LMC were fixed to the values of a selected very good model (f > 0.60). The model was chosen in such a way that its LMC proper motion components are close to the midpoint values of the intervals established by Equations (1). The SMC proper motion was varied within the entire ranges for μ^smc_W and μ^smc_N involved in our study. In all cases, the Magellanic Clouds were found to be receding from one another (Figure 12, left plot), i.e.,

In the other case, we have kept the values of all the parameters equal to those ones of the mentioned high-fitness model, but the LMC proper motion was modified to reach the minimum galactocentric velocity. Calculating the rate dD_l-s/dt then yielded the result depicted in the right-hand plot of Figure 12. Obviously, no change occurred and Equation (20) still holds.

Zoom In Zoom Out Reset image size

Thus, we have seen that strong support exists for a remarkable feature of the kinematics of the Magellanic Clouds. The recent estimates of the LMC/SMC proper motions by Kallivayalil et al. (2006a, 2006b) and by Piatek et al. (2008) are very likely to introduce a close (D_l-s ≈ 10 kpc) past encounter between the Magellanic Clouds at the time of ≈−100 Myr for an arbitrary choice for the LMC/SMC proper motions.

It is reasonable to expect a long-term gravitational binding between the Clouds of the order of 10⁹ yr to be a natural condition for repeated encounters in the LMC–SMC system. The corresponding alignment of the past orbital tracks of the Clouds in the phase space requires a very special setup of the parameters for the model of the interacting system, making such a case very unlikely. However, Figure 13 reveals a surprising picture of the LMC/SMC orbital motion in the deep past. If the separation between the Clouds for the best models (f > 0.60) is plotted as a function of time, similar behavior is found in all cases. After the previously discussed close approach a couple of Myr ago, the mutual distance of the Magellanic Clouds reaches 150 kpc, exceeding the value of 200 kpc for some configurations. The Clouds become unbound shortly after their most recent encounter. Nevertheless, every model contains a second close approach at the time t < −2.5 Gyr with the separation of the Magellanic Clouds dropping below ≈15 kpc. As we will show later, this encounter is the event triggering the formation of an extended trailing tidal tail (the Magellanic Stream) and of its leading counterpart (the Leading Arm).

Zoom In Zoom Out Reset image size

We have stated in Section 6.1 that only a subset of the models identified by the genetic algorithm in the parameter space may be considered satisfactory reproductions of the observational data (Brüns et al. 2005). In general, an acceptable model of the interaction between the Magellanic Clouds and the Milky Way is supposed to reproduce the extent of the observed Magellanic Stream in the three-dimensional space composed of the LSR radial velocity and of the position in the plane of sky. The second large-scale H i structure associated with the Clouds—the Leading Arm—is not used to separate good models from the unsuccessful ones.

The Leading Arm appears in tidal models of the LMC–SMC–Galaxy interaction as a natural counterpart of the Magellanic Stream (see Gardiner & Noguchi 1996; Connors et al. 2006), but its acceptable reproduction by the means of numerical simulations remains a challenge. Our parameter study has not changed that picture—the modeled distribution and kinematics of H i in the Leading Arm region and around the main LMC and SMC bodies remains similarly unsatisfactory over the parameter space, especially failing to reproduce the observed morphology of the Leading Arm. In terms of the fitness function, the value of f never exceeds ≈0.35 if calculated only for the Leading Arm.

Reliable modeling of the evolution of the central regions of the Magellanic Clouds currently eludes the level of sophistication of our numerical code. The inner parts of the Clouds contain high-density baryonic mass and in order to study their evolution a detailed model involving self-gravity and hydrodynamical processes which account for the mass exchange cycle between stars and gas is required. Our restricted N-body code overestimates the amount of gas tidally stripped from the LMC/SMC central regions because the tides are not balanced by self-gravity and the dissipative behavior of gas. As a consequence, the modeled column density is lower than observed values in the inner 5 kpc of the Clouds. This effect is stronger in the case of the SMC as it is the dominant source of matter for the Magellanic Stream and the Leading Arm (see Figure 14 or 16). The particle distribution in the LMC is affected by tidal heating by the Milky Way and the original particle disk is turned into a three-dimensional spherical structure.

Zoom In Zoom Out Reset image size

The Magellanic Stream turned out to be the appropriate indicator of the quality of our models. Its formation is very sensitive to the choice of the model parameters and critically influences the resulting fitness. Figure 14 illustrates these facts. While the model of f = 0.61 is able to fit the basic features of the Magellanic Stream both in the projected H i distribution and the LSR radial velocity profile, a typical simulation representing the model group of f < 0.50 places the Magellanic Stream to the position–position–LSR radial velocity space incorrectly. The Stream is extended insufficiently both in position and the LSR radial velocity (upper row of Figure 14), and the slope of the simulated profile of the LSR radial velocity along the Stream exceeds the already observed one at the Magellanic longitude of ≈−30° (lower right plot of Figure 14). Generally speaking, the described behavior of the modeled trailing stream is responsible for the resulting fitness of a given model, and was used to define the desired threshold level of the fitness value.

The obvious lack of models within the fitness range of (0.50, 0.55), clearly visible in Figures 4 and 10, was combined with the specified features of the high-fitness models. This yielded the desired threshold level of the fitness function quite naturally. Figure 4 shows that the fitness of the model establishing the upper limit to the mentioned gap in the distribution of the fitness values lies slightly below 0.55 (it is 0.546 actually). However, the transition between the good and unsatisfactory models is always gradual, as we have seen, and thus the value f = 0.55 can be taken for the threshold value of the fitness function with no loss of generality.

In order to discuss our results regarding the reproduction of the large-scale distribution of H i associated with the Magellanic Clouds, we compared a selected high-fidelity model with the low-resolution data cube of the H i observations by Brüns et al. (2005). The simulated three-dimensional data cube depicted in Figure 15 (already shown in Figure 14) was produced by a model with the following values of the LSR circular velocity and the flattening of the Galactic halo: Θ = 232 km s⁻¹, q = 0.81. These values yield the total enclosed mass of the Milky Way within the radius of 250 kpc of 2.20 × 10¹² M_☉.

Zoom In Zoom Out Reset image size

Assuming the recent measurements of the LMC and the SMC proper motions, our tidal model was able to produce a trailing tail similar to the observed Magellanic Stream regarding both the LSR radial velocity profile and the extent in the projected position. While the LSR radial velocities measured along the Magellanic Stream are reproduced very well, the far tip of the simulated stream at the position of b ≳ −60° and l ≈ 70° is displaced compared to the observations (Figure 15). Such a displacement appears in every high-fitness model identified by the genetic algorithm in the studied parameter space.

The column density of H i in the Magellanic Stream decreases rather smoothly towards its far tip. This feature also exists in our high-fidelity simulations of the tidally induced formation of the large-scale Magellanic structures. The middle row of Figure 15 shows how the mean H i column density on the line perpendicular to the great circle of the Magellanic longitude (defined in Wannier & Wrixon 1972) changes with the longitude. The simulated profile of the column density (right-hand plot) has a steeper slope than the observed Stream and its gradual decrease ends at the Magellanic longitude of ≈−70° compared to ≈−90° for the observations.

The plots in the lower row of Figure 15 also offer the view of the area in the plane of sky where the H i structure called the Leading Arm is located. Brüns et al. (2005) introduced the division of the structure into three systems of gaseous clumps named LA I, LA II, and LA III. Their approximate positions in the Galactic coordinates [b, l] are [ − 15°, 300°] (LA I), [15°, 290°] (LA II), and [15°, 270°] (LA III). The simulation shown in Figure 15 was able to reproduce the part LA III and a fraction of the clump LA II. Some of our high-fitness models offer quite a good match of the overall kinematics and spatial extent for both observed structures LA II and LA III. Nevertheless, no model reproducing both the system LA I and the Magellanic Stream at the same time was found.

It is interesting to take a look at the composition of the simulated Magellanic Stream, and at the process of its buildup over the course of time. The lower left plot of Figure 14 reveals that most of the particles forming the Magellanic Stream in the selected high-fitness model are of SMC origin. A number of the LMC particles are spread over the full extent of the Stream as well. This is a typical result and it is based on the fact that the LMC particles are more tightly bound due to the significantly higher mass of the LMC compared to the SMC—the LMC:SMC mass ratio is 10:1 in our simulations.

It has already been discussed by Connors et al. (2006) that the Magellanic Stream might consist of several filaments projected to the same region in the plane of sky. Such filaments would be the tidal debris spread along the past orbits of the Clouds. Our study has confirmed such a scenario. Figure 15 shows the projected orbits of the Magellanic Clouds from t = −3 Gyr to the present, plotted over the contours of the integrated column density of H i (lower right plot). The SMC has crossed the area currently occupied by the Magellanic Stream three times, while two passages have occurred in case of the LMC. If we disentangle the projected orbital paths, we will find that the previous orbital revolutions appear more diagonal in our plot. They are placed at lower Galactic longitudes. The resulting trailing tail is thus a mixture of particles released from the LMC/SMC due to the tidal stripping during several epochs covering the time range from present to ≈−3.3 Gyr when the redistribution matter was triggered by a close LMC–SMC encounter.

The complex structure of the modeled Stream can be seen in the lower plot of Figure 16. The LMC/SMC particles are color-coded according to the epoch when they were stripped from the Clouds. The Stream and the Magellanic Bridge are dominated by the SMC particles, while most of the LMC particles remain in the vicinity of the main LMC body and have been gradually tidally heated due to the gravitational field of the Milky Way.

Zoom In Zoom Out Reset image size

The past orbital history of the Magellanic Clouds for the discussed high-fitness model can be seen in the lower plot of Figure 16. For the last 4 Gyr, the Clouds have moved within the Galactocentric radius of 200 kpc, undergoing two close (D_l-s ≈ 10 kpc) encounters at the times of −3.3 Gyr and −0.15 Gyr, respectively. The impact of the encounters on the distribution of the LMC/SMC particles was different.

To quantify the influence of the encounters in the LMC–SMC–Galaxy system, we have counted the number of particles shifted from the original orbits of the radius r_i to the distance of 2r_i from the center of their progenitor. If this function is evaluated for the time bins of 500 Myr from t = −4 Gyr, we obtain the upper plot of Figure 16. Apparently, the number of the disturbed particles was already increased in the time bin of 〈 − 3.5, − 3.0〉 Gyr as a consequence of the first LMC–SMC close approach. This level was later only mildly decreasing, as more of the disturbed particles were stretched in the tidal field of the Milky Way, satisfying the above stated definition of strongly affected particles.

The second close LMC–SMC encounter was significantly more dramatic, disturbing more than 10% of all particles. This event is the origin of the common envelope surrounding both the Clouds, and of the turbulent gaseous filament connecting the main bodies of the Clouds, named the Magellanic Bridge (Brüns et al. 2005). Substantial fraction of the particles contributed to the Magellanic Stream as well.

The results summarized in the preceding paragraphs and depicted in Figure 16 have a remarkable consequence regarding the formation history of the Magellanic Stream. Over the last ≈3 Gyr, it was supported at an almost constant rate by the gas tidally stripped from the Clouds. However, this rate was increased by a factor ≈20 due to a recent encounter between the Clouds at t = −0.15 Gyr. We would like to mention that on the qualitative level these results are an appropriate description of all the high-fitness models found in the parameter space for the LMC–SMC–Milky Way interaction.

We have found that the high-fitness simulations prefer the oblate (q < 1) shape of the logarithmic halo of the Galaxy over the prolate configuration (see Figure 5), but several reasonable models of q as high as 1.10 exist. As long as good models of the interacting system exist for the prolate halo, this shape of the Galactic gravitational potential cannot be ruled out. However, we have shown that such a possibility is substantially less likely than the case of q < 1. Moreover, only two acceptable models with significantly prolate (q ≳ 1.10) DM halos of the Milky Way have been identified for a reasonable total mass of the halo, i.e., below the limit of ≈3 × 12 M_☉ (see Section 3.1).

The impact of the flattening parameter on the phase-space positioning of the LMC/SMC orbits has been addressed in Section 6.4 of this paper. The shape of the Galactic halo also influences the redistribution of the matter associated with the Magellanic Clouds due to the resulting tidal field of the Milky Way. This effect is the key to the preference for the oblate configuration of the DM halo of our Galaxy. We have calculated the tidal force exerted by the Milky Way on two points separated by Δr = 10 kpc and moving around the Galactic Center on circular polar orbits of the radii 50 kpc, 100 kpc, and 150 kpc. Figure 17 shows the result. In general, the oblate configuration offers the most efficient tidal stripping, improving the probability that a sufficient redistribution of the LMC/SMC particles occurs within the available timescale of several Gyr.

Zoom In Zoom Out Reset image size