Dear Magellanic Clouds, welcome back!

We model both galaxies as pure $N$-body systems, neglecting any effect of the gas. This is justified by the expectation that the orbit is determined by the total gravitational potential, which is dominated by dark matter, with a secondary contribution of stars, and an even smaller contribution of gas. Tepper-García et al. (2019) find that the addition of a hot gas corona in a hydrodynamical simulation does not appreciably change the orbit of the LMC compared to the purely collisionless case (although, of course, it affects the properties of the Magellanic gas stream, which we cannot study in our setup).

Both galaxies are set up as equilibrium models initially, using various tools from the stellar-dynamical toolbox Agama (Vasiliev, 2019). This is trivial for the LMC, which is represented by a spherical density profile with the distribution function computed from the Eddington inversion formula. We do not explicitly include the stellar disc or halo of the LMC, because these are sub-dominant in the total mass distribution (van der Marel et al. 2002 estimated the stellar mass to be $3\times 10^{9}\,M_{\odot}$) and located sufficiently deep in the potential well to avoid being stripped. The density profile of the LMC is designed to match the circular-velocity curve in the inner few kpc determined from stellar kinematics (Vasiliev, 2018) and the enclosed mass constraint of $1.7\times 10^{10}\,M_{\odot}$ at 8.7 kpc from the LMC centre (van der Marel & Kallivayalil, 2014), so the inner region of our LMC model effectively represents both stars and dark matter. The density follows the Navarro–Frenk–White (NFW) profile with a smooth but rather sharp cutoff:

The Milky Way models consist of a spherical dark matter halo with the same tapered NFW density profile, and two stellar components: a spherical bulge with density profile

and a total mass $M_{\mathrm{bulge}}=1.2\times 10^{10}\,M_{\odot}$, and a single exponential disc with density profile

and a total mass $M_{\mathrm{disc}}=5\times 10^{10}\,M_{\odot}$. The circular velocity at the Solar radius is $\sim 230$ km s^-1, of which stars contribute $\sim 180$ km s^-1. The bulge and the halo have isotropic velocity distributions by default, but we also created a variant of the Milky Way model with a radially anisotropic halo having $\beta\equiv 1-\sigma_{\mathrm{tan}}^{2}/(2\sigma_{\mathrm{rad}}^{2})=0.5$. The two-component equilibrium model for the Milky Way is constructed with the Schwarzschild orbit-superposition method, as implemented in Agama.

We consider two variants of models for each galaxy, with the halo parameters listed in Table 1. The Milky Way’s mass profile, shown in Figure 1, matches contemporary observational constraints (see e.g. Wang et al. 2020 for a compilation of recent estimates), with a virial mass being either $1.0\times 10^{12}\,M_{\odot}$ for the lighter model $\mathcal{M}$10 or $1.1\times 10^{12}\,M_{\odot}$ for the heavier model $\mathcal{M}$11. The initial mass of the LMC model is higher than the current estimates (see Figure 1 in V23 for a summary), but as it is tidally stripped on the first pericentre passage, its mass drops down to the acceptable range.

The simulations are conducted using the fast-multipole code gyrfalcON (Dehnen, 2000), augmented with a custom plugin for on-the-fly tracing of the centres of both galaxies. We use up to $10^{7}$ particles in the final highest-resolution runs, but the early stages are run with a reduced number of particles (by factors of 4 or 20). We set the softening length to 0.5 kpc for halo components and 0.2 kpc for the Milky Way stars (the equivalent Plummer softening length is $\sim 1.5\times$ smaller), and use a block timestep scheme with the smallest and largest timesteps of 0.5 Myr and 8 Myr, respectively. Since the public version of the code is not parallelized, a single full-resolution simulation takes $\sim 40$ hours.

We adopt the LMC centre position $\alpha=81.28^{\circ}$, $\delta=-69.78^{\circ}$ and PM $\mu_{\alpha,\delta}=\{1.858,0.385\}$ mas yr^-1 from Gaia Collaboration (2021), distance of 49.6 kpc from Pietrzyński et al. (2019), and heliocentric line-of-sight velocity of 262.5 km s^-1 from van der Marel et al. (2002). Using a right-handed Galactocentric reference system with the Solar position at {$-8.12$, 0, 0.02} kpc and velocity {12.9, 245.6, 7.8} km s^-1 (Astropy collaboration, 2022; Drimmel & Poggio, 2018), this translates into the LMC position {$-0.6$, $-41.0$, $-26.8$} kpc and velocity {$-69.9$, $-221.7$, $214.2$} km s^-1. Of course, these quantities have non-negligible uncertainties, which have a significant impact on the past orbit reconstruction (see Figure 2 in V23 and associated discussion), but to facilitate the comparison between different simulations, we aim at reproducing the adopted position and velocity as precisely as possible.

The first guess for the initial phase-space coordinates of the LMC $\boldsymbol{w}_{\mathrm{init}}\equiv\{\boldsymbol{x}_{\mathrm{init}},\boldsymbol{v}_{\mathrm{init}}\}$ may come from integrating its test-particle orbit back in time in the static Milky Way potential, but this is very inaccurate for the massive LMC. A better guess could be obtained by integrating the coupled equations of motion of both galaxies as if they were rigid (non-deforming) but moving points sourcing their respective potentials, adding a drag force responsible for the dynamical friction (Gómez et al., 2015; Jethwa et al., 2016; Erkal et al., 2019). This approach is adopted in many studies (e.g., Laporte et al., 2018b; Garavito-Camargo et al., 2019; Petersen & Peñarrubia, 2020), but since the galaxies in live $N$-body simulations are distorted in the interaction, their orbits deviate from the adopted ones. Previous studies attempted to minimize the deviation by tuning the amplitude of dynamical friction, or by iteratively refining the initial conditions using a grid search in the velocity space and picking up the closest point (e.g., Donaldson et al., 2022; Lilleengen et al., 2023)¹¹1Instead of a grid search, Guglielmo et al. (2014) used a genetic algorithm to tune the initial conditions, however, their simulations were not fully self-consistent (the Milky Way was modelled as a static potential).. However, the present-day position and velocity of the LMC $\boldsymbol{w}_{\mathrm{final}}\equiv\{\boldsymbol{x}_{\mathrm{final}},\boldsymbol{v}_{\mathrm{final}}\}$ are typically still off by a few kpc and few tens of km s^-1– an error too large for an accurate reconstruction of the orbit.

In the present study, we aim to radically improve the precision of matching the final phase-space coordinates by iteratively adjusting the initial conditions $\boldsymbol{w}_{\mathrm{init}}$. In each iteration, we run a master simulation with a given $\boldsymbol{w}_{\mathrm{init}}$ and a suite of companion simulations with slightly different initial conditions, and determine their final phase-space coordinates. Then we measure the Jacobian of the mapping from $\boldsymbol{w}_{\mathrm{init}}$ to $\boldsymbol{w}_{\mathrm{final}}$ and use it to assign the next guess for $\boldsymbol{w}_{\mathrm{init}}$ (a new value not chosen from the current suite). The entire process is then repeated several times until the mismatch in $\boldsymbol{w}_{\mathrm{final}}$ drops to an acceptable level ($\lesssim 1$ kpc and 1 km s^-1). The details of this procedure and associated challenges are described in the Appendix.

The LMC potential represents all particles that originated from the LMC, but many of them are no longer bound to it. To determine the bound mass, we use the following procedure: starting with all LMC particles, we construct their total potential and compute the particle energies $E$ in the reference frame associated with the LMC (i.e., subtracting its centre velocity from particle velocities when computing the kinetic energy). Then particles with $E>0$ are discarded, and the potential is recomputed from the remaining particles, followed by another round of eliminating particles with $E>0$ and repeating the cycle until convergence. The concept of bound mass is ill-defined in the case of an eccentric orbit, because some particles classified as unbound near the pericentre may be subsequently recaptured, but the above procedure gives a reasonable qualitative approximation. The lower set of curves in Figure 1 illustrate the enclosed mass profiles for two models, $\mathcal{L}$2 and $\mathcal{L}$3, in three variants: initial snapshot, current distribution of all particles originally belonging to the LMC system, and only the currently bound ones. These curves largely coincide within 20 kpc from the LMC centre, but start to deviate further out, with the bound mass essentially reaching a limit around 50 kpc. A comparison with the gray-shaded region representing dynamical constraints on the enclosed mass from the Orphan–Chenab stream modelling (Koposov et al., 2023) demonstrates that the two sets of LMC models are on the opposite ends of this range, and likely bracket the actual profile.

For subsequent analysis, we construct smooth potential approximations for each snapshot in the simulation (stored every 64 Myr), separately for particles belonging to the LMC and to the Milky Way halo. It is convenient to work in the reference frame associated with the Milky Way centre. Because this frame is non-inertial, the total potential contains a spatially uniform but time-dependent UniformAcceleration term, which is obtained from the second derivative of the Milky Way trajectory (it is essential to smooth the trajectory beforehand, as otherwise the second derivative is far too noisy, see section 3.1.2 in Sanders et al. 2020). The Galactic stellar disc is approximated by the initial potential at all times, the Galactic halo is represented by a sequence of Multipole potentials pinned to origin, and the LMC is represented by a sequence of Multipole potentials moving along the smooth trajectory extracted from the simulation. We also created alternative BasisSet potential expansions, which use the Hernquist & Ostriker (1992) basis set and are compatible with other dynamical modelling packages such as Gala (Price-Whelan, 2017), though these are less efficient to evaluate than the Multipole expansions.

A natural question is whether such a time-dependent potential approximation provides a good description of the actual orbits of particles. Lowing et al. (2011) and Sanders et al. (2020) performed such accuracy analysis for the case of a single host halo from a cosmological simulation, in which there are no LMC-like massive satellites flying around, but numerous subhaloes create small-scale structures in the potential not captured by the low-order basis-set expansion. Nevertheless, it was found that most particle trajectories are well reproduced over the duration of the simulation. In our case, there are no small-scale substructures (apart from the dynamical friction-induced wake in the Milky Way halo), and the LMC is granted its own potential expansion, which moves along a smooth trajectory, so we may expect the fidelity of orbit reconstruction to be even better. Figure 2 shows a few examples of orbits from the actual simulation and their reconstructed counterparts, demonstrating that they stay fairly close (to within a few kpc) for the entire duration of simulation in most cases. Even if individual orbits may not be perfectly reconstructed, the ensemble of nearby orbits corresponding to a given satellite (Section 4.2) is very well reconstructed in the statistical sense. Re-running the entire simulation as an ensemble of non-interacting particles closely reproduces the current distribution and kinematics of Magellanic debris and the corresponding perturbations in the Milky Way halo.

After experimenting with different orders of expansion, we find that $\ell_{\mathrm{max}}=6$ is more than sufficient, and we use 25 log-spaced nodes between 0.1 and 1000 kpc in the radial grids. Sanders et al. (2020) also compared the accuracy of orbit approximation in the case when the potential was taken from the nearest snapshot or linearly interpolated between snapshots, and found that the latter makes little difference but is twice more expensive (unless the expansion coefficients themselves are linearly interpolated before computing the potential, which is not available in Agama). However, we noticed that since the orbit integration in Agama uses a high-order method (dop853), the extra cost of a double potential evaluation is more than compensated by its increased smoothness, allowing the integrator to take larger timesteps. Even so, the computational cost (number of timesteps) scales with the requested accuracy $\epsilon$ as $\epsilon^{-1/n}$ with $n\simeq 4$ for the nearest-timestamp potential and $n\simeq 6$ for the linearly interpolated potential, faster than the $n=8$ scaling expected for this integrator in an infinitely smooth potential. Nevertheless, a value of $\epsilon=10^{-4}$, while much coarser than the default $\epsilon=10^{-8}$ used for analytic potentials, is sufficient for the relatively short duration of orbit integration (at most a few tens dynamical times).

We first examine the past orbits of the LMC for different combinations of its mass and Galactic potential. Figure 3, top panel, shows the Galactocentric distance of the LMC in four simulations, together with test-particle orbits in the two Milky Way potentials $\mathcal{M}$10, $\mathcal{M}$11, shown by dot-dashed curves. It confirms the trend reported in V23: more massive LMC have smaller apocentre distances and shorter periods. This results from the interplay of two opposite effects. On the one hand, a more massive LMC experiences stronger dynamical friction, thus it must have started with a higher orbital energy to arrive to the same present-day phase-space point. On the other hand, though, a more massive LMC induces a larger reflex motion in the Milky Way, thus to achieve the same measured relative velocity, it needs to have a lower velocity in the centre-of-mass frame, meaning lower total energy and shorter period. Equivalently, the orbital period of the binary system depends on its total mass and is shorter for higher-mass LMC models. It turns out that the second factor outweighs the first one, though the difference in orbital period is rather minor between the two LMC models.

By contrast, the dependence of the orbital period on the Milky Way potential is quite dramatic, with $\gtrsim 30\%$ difference between $\mathcal{M}$10 and $\mathcal{M}$11 (which differ by 10% in virial mass). As discussed in Section 3.2 of V23, this is due to the orbit being only marginally bound and the LMC currently located close to its pericentre, therefore its kinetic energy is almost equal in magnitude to its potential energy, and even a small variation in the latter leads to dramatic changes in the total energy. It is certainly possible to make the period comparable or longer than the Hubble time (i.e., put the LMC onto a first-approach trajectory), but the period could also be as short as 5–6 Gyr for a reasonable choice of Galactic potential. In these idealised simulations, we neglected the cosmological evolution of the Milky Way mass profile, which would have lengthened the orbital period (see Figure 11 in Kallivayalil et al. 2013), but also neglected other factors that could shorten it – the perturbation of the LMC velocity by the SMC or a favourable orientation of a non-spherical Milky Way potential (both discussed in V23). Overall, it is plausible that the previous passage of the LMC could have occurred between 5 and 10 Gyr ago at a distance just over 100 kpc, although the details may differ from our setup.

The middle and lower panels show the evolution of angular momentum and bound mass in these four simulations. Both quantities experience stepwise drops around each pericentre passage, but interestingly, the decrease of angular momentum occurs mostly after the pericentre, whereas the decrease of bound mass mostly precedes it. The bound mass during the past orbital period stays in the range (1.2–2)${}\times 10^{11}\,M_{\odot}$, compatible with the current estimates (see Figure 1 in V23 for a compilation of literature values).

A popular approach for reconstructing the past orbit of the LMC, taking into account dynamical friction from the Galactic halo and reflex motion that the LMC induces on the Milky Way, is the solution of a coupled ODE system for two extended masses (Equation 1 in V23), using the gradient of gravitational potential of each galaxy at the current distance $r(t)$ to compute the accelerations of both galaxies. It has been followed in many previous studies, e.g., Gómez et al. (2015), Jethwa et al. (2016), Erkal et al. (2019), Patel et al. (2020), Correa Magnus & Vasiliev (2022). As the Coulomb logarithm is generally treated as a tunable parameter, Hashimoto et al. (2003) found that a distance-dependent value $\ln\Lambda=\ln(r/r_{\mathrm{min}})$ provides a better match the results of actual $N$-body simulations, with $r_{\mathrm{min}}$ comparable to the scale radius $r_{\mathrm{s}}$ of the satellite. The complication usually ignored by these studies is that the LMC mass also decreases with time, especially in the second-passage scenario. Figure 4 compares the actual orbits in our $N$-body simulations (solid blue) with the solutions of the approximate ODE system (dashed red), in which the LMC mass has been set to half its initial value and $r_{\mathrm{min}}$ set to $0.8\,r_{\mathrm{s}}$. There are noticeable deviations in the reconstructed orbital periods (top row), and the LMC-induced acceleration of the Milky Way is overestimated by $\sim$20% at its peak (bottom row), because this prescription ignores the deformation of both galaxies. Although with two free parameters, we managed to match the grid of 4 simulations reasonably well, there is no guarantee that this recipe would be adequate for all reasonable choices of Milky Way and LMC mass profiles. It is clear that a more accurate approximation is highly desirable for a fast and reliable quantitative analysis, short of running the actual simulations.

As established in several recent studies (Garavito-Camargo et al., 2019, 2021a; Cunningham et al., 2020; Erkal et al., 2020, 2021; Petersen & Peñarrubia, 2020, 2021), a massive LMC induces significant perturbations in the Galactic halo at distances $\gtrsim 20$–30 kpc (see Section 4.2 in V23 for an extended discussion). A natural question is, therefore, whether the second-passage scenario produces different signatures from the first-passage scenario considered in the above papers.

Figure 5 shows the contours of Milky Way halo density projected onto the LMC orbital plane (top panel) and heliocentric sky maps of density and line-of-sight velocity perturbations (middle and bottom rows) in the $\mathcal{L}$3–$\mathcal{M}$10 simulation. At a glance, these plots look very similar to the ones from the first-passage simulations (e.g., Figure 5 in V23). For a more detailed analysis, we constructed another simulation of a first passage, using the same Galactic model $\mathcal{M}$10 and a lower-mass LMC model $\mathcal{L}$2, started near the apocentre at time 4 Gyr ago. Of course, such a trajectory extrapolated further into the past would have had another pericentre passage, but we ignore this inconsistency for the sake of experiment. The density and kinematic perturbation maps in this scenario look so similar to Figure 5 that it is not even worth showing them separately. In other words, any possible signatures of the first passage either dissipate or are totally obscured by the most recent passage, which occurred at a twice smaller distance.

Models with the same LMC mass but different Milky Way potentials again look quite similar, and the variation of LMC mass only changes the amplitude of perturbations, but not their qualitative appearance. On the other hand, there are interesting differences between the models with isotropic Milky Way halo discussed so far, and a model with a radially anisotropic halo ($\beta=0.5$, typical for the outer region of NFW haloes, e.g., Hansen & Moore 2006), which is shown in Figure 6. The velocity perturbations remain largely the same as in the previous figure, but the dynamical friction-induced overdensity along the past orbit of the LMC is much more prominent, and even more intriguingly, there is another overdensity along the future trajectory (upper left corner of the sky map). The large-scale dipole asymmetry is not significantly affected, but the local wake is strongly enhanced, in complete agreement with the study of Rozier et al. (2022), who used a very different method (linear response theory instead of conventional $N$-body simulations). It is unclear which of the two cases better matches the current observations: a density asymmetry is much more difficult to measure reliably than a kinematic perturbation, and the only study which attempted this (Conroy et al., 2021) does not match either plot in detail (although they do note that the measured density contrast in the local wake is significantly higher than even in their simulation with the highest LMC mass of $2.5\times 10^{11}\,M_{\odot}$). The previous orbital period of the LMC in the radially anisotropic halo is similar to its isotropic counterpart, but the initial angular momentum is higher and the distance of the first pericentre is correspondingly larger by $\sim 20\%$; in other words, radial anisotropy of the host galaxy enhances the radialization of satellite orbits, as already mentioned in Vasiliev et al. (2022).

It is also worth noting that the observationally determined direction of the dipole perturbation in line-of-sight velocities (Petersen & Peñarrubia, 2021) differs by more than 45^∘ from the roughly north–south direction found in all simulations, including those in the present paper; the reasons for this discrepancy are not clear and may point to some deficiencies of current models, first or second passage alike.

The LMC also imparts a noticeable perturbation on the Milky Way disc, likely contributing to the warp observed in its outer regions both in stars and gas. One should keep in mind that the very definition of the disc plane is not unambiguous. Here we use the particles within 10 kpc from the Galactic centre to construct a plane that minimizes the scatter of their $z^{\prime}$ coordinates relative to that plane. Although the initial models had $z^{\prime}$ aligned with the $z$ axis, by $t=0$ this best-fitting plane is already tilted by 1–1.5^∘, which has an unintended consequence that the LMC coordinates in this rotated system are slightly different from the input values. Fortunately, the tilt direction is such that this offset roughly corresponds to a shift along the same trajectory by $\sim$3–4 Myr, so we ignore this minor deviation instead of trying to fit its present-day position in the rotated coordinate system.

Figure 7 shows the mean $z$ coordinate ($z_{0}$) of stars across the Galactic disc in the $\mathcal{L}$3–$\mathcal{M}$11 simulation, which is qualitatively similar to other simulations. The noticeable warp in the top left corner is nearly opposite to the current location of the LMC ($x\approx-1,y\approx-41$ kpc), and is caused by the same effect as the distortions in the Galactic halo, namely, that the inner part of the disc has shifted downward relative to the regions further away from the LMC. This displacement map qualitatively resembles Figure 6 in Laporte et al. (2018a), who conducted a similar set of simulations, though without matching the present-day LMC coordinates as precisely as in the present study. We stress that the perturbations rapidly evolve, e.g., the direction of the warp changes by $\sim\!90^{\circ}$ over 0.1 Gyr, so for a quantitative analysis and comparison with observations, it is important to have a precisely fitted model. Here we do not attempt a detailed comparison, only noting that the warp in the model is similar to the one observed in stars (e.g., Figure 7 in Romero-Gómez et al. 2019, or Figure 2 in Lemasle et al. 2022), but somewhat smaller in amplitude. Laporte et al. (2018b) found that a combination of perturbations from the LMC and Sagittarius may create a stronger warp signature than a simple superposition, thus we do not expect our LMC-only models to match the data perfectly.

Given that our simulations produce an LMC model not only matching its observed present-day position and velocity, but also dynamically consistent with its prior evolution (i.e., with appropriate degree of tidal stripping), it makes sense to explore the future fate of the LMC and its effect on the Milky Way. We run the simulations for another 5 Gyr, and depending on the LMC mass, it completes another 3–4 orbits before being entirely disrupted (although the lowest-mass combination, $\mathcal{L}$2–$\mathcal{M}$10, is not fully merged by the end of this interval). Figure 3 shows the future evolution of its mass and Galactocentric distance to the right of the vertical gray line.

Although the current magnitude of perturbations in the Galactic halo will not be exceeded in future pericentre passages due to the rapidly decreasing LMC mass, its effect on the Galactic disc will be more prominent as the pericentre radius drops to $\sim$25 kpc on the next passage, and then further to 10–15 kpc. Figure 8, top row, shows the maps of the disc warp similar to Figure 7, but at 5 Gyr in the future for the most massive combination $\mathcal{L}$3–$\mathcal{M}$11, in which the LMC is fully disrupted by that time. The perturbations in the disc are much stronger than now, and it has become significantly hotter (lower right panel shows that the thickness has doubled across all radii) and tilted by $\gtrsim 10^{\circ}$ relative to its initial orientation, but not destroyed.

Lower left panel of that figure shows the density profile of Magellanic debris, which has two characteristic break radii at 10 kpc (the radius of the last pericentre) and 40–50 kpc (the current pericentre radius). Although the initial LMC models did not have a separate stellar component, one could get a crude but sensible approximation by tracking 1% of most bound particles (the mass $3\times 10^{9}\,M_{\odot}$ roughly matches the estimate of total stellar mass from van der Marel et al. 2002). The density of these particles is shown by red dashed line, and stays well below the density of Milky Way disc (dashed green) in the inner Galaxy, but may be comparable or even exceed the density of the pre-existing stellar halo at large distances. The total mass of present-day Milky Way stellar halo is estimated to be $\sim 10^{9}\,M_{\odot}$ (Deason et al., 2019), so the LMC will quadruple it when fully disrupted. These findings agree well with the predictions of Cautun et al. (2019) based on Milky Way–LMC analogues in the EAGLE simulation, although in our experiments, it takes somewhat longer (4–5 Gyr) to fully destroy the LMC.

Lynden-Bell (1976) noticed that the dwarf galaxies Draco, Ursa Minor and Sculptor lie close to the plane of the Magellanic gas stream, itself discovered shortly before (Mathewson et al., 1974), and proposed that they are tidally stripped satellites of the postulated Great Magellanic Galaxy – the precursor of the LMC itself. Lynden-Bell (1982a) added Carina dSph to the list of galaxies possibly associated with the Great Magellanic Galaxy, while Lynden-Bell (1982b) introduced another dSph grouping – Fornax, Leo I, Leo II and Sculptor (FLS), but treated them separately from the Magellanic group, and Majewski (1994) added Sextans and Phoenix dSph to the FLS group. However, the planes of two groups differ by only $30^{\circ}$, so Kroupa et al. (2005) proposed that they form a single group, later dubbed “Vast Polar Structure” by Pawlowski et al. (2012). Although a spatial alignment could happen by chance, it gradually became clear that most of these satellites also have similar directions of angular momentum (orbital poles), making this structure dynamically coherent (Metz et al., 2008; Pawlowski & Kroupa, 2013, 2020; Fritz et al., 2018).

D’Onghia & Lake (2008) proposed a Magellanic group consisting of 7 galaxies (two Clouds, Draco, Leo II, Sagittarius, Sextans and Ursa Minor), although they did not comment on why the orbital plane of Sagittarius is nearly orthogonal to that of the LMC. Metz et al. (2008) presented the now canonical list of clustered orbital poles of classical dwarfs: the Clouds, Carina, Draco, Fornax, Ursa Minor, and noted that if the Clouds were on their first passage and unrelated to the other four galaxies that have been orbiting the Milky Way for a while, it would be a rather unlikely coincidence for them to share the same orbital plane just by chance. Although the idea that a relatively thin and kinematically coherent disc of satellites could result from an accretion of a massive group of galaxies early in the Milky Way history is quite natural (e.g., Li & Helmi, 2008; Smith et al., 2016), somehow it did not get much traction, with Metz et al. (2009) arguing against it on the grounds that the observed plane of satellites in the Milky Way is too compact compared to the expectations from group accretion. Instead, Pawlowski et al. (2011) proposed a scenario in which the VPOS consists of “tidal dwarf galaxies” formed out of gas expelled during a tidal interaction between the Milky Way and another galaxy early in its evolution, somewhat reminiscent of the now abandoned Jeans–Jeffreys tidal theory of the Solar system formation (Jeffreys, 1952). They considered the possibility that the proto-LMC could be this intruder galaxy, though in their scenario the pericentre passage leading to the tidal stripping would have occurred at a much smaller distance ($\sim 10$ kpc) than the most recent pericentre distance of the LMC (45–50 kpc).

Until about 10 years ago, the census of Milky Way satellites was limited to relatively bright objects, whose association with the Magellanic system is difficult to establish with confidence (apart from the aforementioned similarity of orbital pole directions). Table 2 summarizes the proposed claims of Magellanic origin for classical dwarfs. However, since then many new ultrafaint galaxies were discovered in deep photometric surveys such as DES, and quite a few of them are very likely to be LMC satellites (see Table 1 in V23 for a compilation of results).

Numerous studies that considered the association of dSph and ultrafaint galaxies with the LMC system can be broadly divided into several categories:

1. 1.

   One may take the present-day positions and velocities of these objects as well as those of the LMC itself, then reconstruct the past orbit of the LMC in the assumed Galactic potential (ideally, not as a test particle moving in a fixed potential, but taking into account the reflex motion it induces on the Milky Way), and then follow the orbits of dwarf galaxies in the combined potential of the Milky Way and the LMC. An object can be declared an LMC satellite if it is currently bound to it, or a former satellite if it was bound at some point in the past (i.e., its relative velocity w.r.t. the LMC was below the escape velocity of the latter). Examples of this approach include Erkal & Belokurov (2020), Patel et al. (2020), Battaglia et al. (2022), Correa Magnus & Vasiliev (2022), Pace et al. (2022). The limitation of this approach is that it needs high-accuracy 6d phase-space coordinates of all objects for a reliable reconstruction of past orbits, which are rarely available, and even then the reconstruction becomes increasingly less reliable further into the past (e.g., D’Souza & Bell, 2022). This creates a bias against the association with the LMC system, since it is much more likely that the measurement errors will scatter an orbit away from the LMC than closer to it.

2. 2.

   To combat the latter problem, one may instead perform forward orbit modelling, namely, consider a large sample of particles initially bound to the LMC system, evolve their orbits forward up to the present time, and then identify which of the observed satellites are close to any of these LMC-associated particles within uncertainties. This alone is not enough to establish the likelihood of association – one needs to compare it with the alternative hypothesis that an object comes from the Milky Way system (or was accreted onto it independently from the LMC). This could be achieved by sampling an equally large number of particles from the Milky Way halo and counting the fraction of particles in the spatial and kinematic neighbourhood of each galaxy coming from the two alternative sources. The first part of this approach was followed in Nichols et al. (2011), Yozin & Bekki (2015), Jethwa et al. (2016), but these studies only examined whether a given dSph lies in the region occupied by the LMC debris (not even comparing the velocities due to a lack of measurements) without considering the alternative (non-LMC) origin.

3. 3.

   One can take this idea further and perform a full $N$-body simulation of the two interacting galaxies, comparing the positions and velocities of observed dSph with the spatial and kinematic distribution of particles originating in both galaxies. Several groups have conducted simulations of this kind (Garavito-Camargo et al., 2019; Petersen & Peñarrubia, 2020; Vasiliev et al., 2021; Pawlowski et al., 2022; Makarov et al., 2023), but none so far have used them to assess the satellite population of the Magellanic system, although Vasiliev et al. (2021) employed this approach to determine the association of globular clusters with the Sagittarius galaxy (which was also simulated as a live $N$-body system in that study). The challenge is to match the current position and velocity of the LMC with sufficient precision for this analysis to be meaningful, and only the latter study took great effort to meet it.

4. 4.

   Finally, one can resort to the analysis of cosmological simulations, either in large volume with thousands of Milky Way analogues, such as Millenium-II (Boylan-Kolchin et al., 2011), or more tailored Local Group-scale simulations, such as Aquarius (Sales et al., 2011, 2017; Kallivayalil et al., 2018), ELVIS (Deason et al., 2015), FIRE (Jahn et al., 2019), Auriga (Pardy et al., 2020) or Apostle (Santos-Santos et al., 2021). In this case the main challenge is to find suitable analogues of the Milky Way–LMC system, abandoning any hope for a detailed match, but only qualitatively reproducing its main properties. Earlier studies, such as Sales et al. (2011), did not find conclusive evidence of Magellanic association for any of the then-known dSph galaxies (except the SMC, of course), but with the more recent discoveries of numerous ultrafaint dwarfs in the vicinity of the LMC, several of them were proposed to be its satellites based on the similarity of their location (and sometimes velocity) to particles from LMC analogues in the simulations. Jahn et al. (2019) and Pardy et al. (2020) also proposed Carina and Fornax as possible former satellites of the LMC, based on their updated PM from Gaia DR2, but Santos-Santos et al. (2021) found this possibility rather tenuous.

In this study, we take the third approach, analyzing the possible origin of each dSph by comparing its phase-space coordinates with particles in the simulation. We also use elements of the second approach, rewinding dSph orbits in the pre-recorded potential of the simulation approximated by smooth potential expansions and comparing the probability of arriving with the Magellanic system against originating in the Milky Way. We take the coordinates, distances and line-of-sight velocities from the catalogue of McConnachie (2012, updated in 2021)²²2https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/community/nearby/, and Gaia eDR3 PM from Battaglia et al. (2022); the SMC PM is taken from Gaia Collaboration (2021).

Let $\boldsymbol{X}$ be the 6d vector of observable phase-space coordinates of a given object (longitude $l$, latitude $b$, distance $D$, two PM components $\mu_{l}$, $\mu_{b}$, and line-of-sight velocity $V$)³³3Here longitude and latitude are not in the Galactic coordinate system, but rather in a system rotated such that the target object is at origin, minimizing distortions of the PM field that could appear for objects near celestial poles., and $\mathsf{E}$ be the covariance matrix of measurement uncertainties (in case of uncorrelated errors, it is diagonal with squared standard deviations $\epsilon_{c}$ in each coordinate $c$). Although the position is measured essentially perfectly, we need to introduce an artificial uncertainty in celestial coordinates in order to have a sufficient number of particles to match, setting $\epsilon_{l,b}=4^{\circ}$. We also put a lower limit on the uncertainty in distance ($\epsilon_{D}\geq 0.1\,D$), line-of-sight velocity ($\epsilon_{V}\geq 10\,$km s^-1) and PM ($\epsilon_{\mu}\geq 0.05\,$mas yr^-1), both to increase the chance of matching and to account for various systematic deviations between the model and the real Milky Way. The $N$-body snapshot is then converted into the same coordinates $\boldsymbol{X}^{(i)}$, $i=1..N_{\mathrm{body}}$. To determine the probability of the object to belong to either the LMC or the Milky Way, we use three different approaches.

In the first method (particle matching), we compute the uncertainty-weighted Mahalanobis distance between the object and each $N$-body particle: $d^{(i)}\equiv\sqrt{\big{(}\boldsymbol{X}-\boldsymbol{X}^{(i)}\big{)}^{T}\;\mathsf{E}^{-1}\;\big{(}\boldsymbol{X}-\boldsymbol{X}^{(i)}\big{)}}$. The matching likelihood for $i$-th particle is $p^{(i)}\propto m^{(i)}\,\exp\big{[}-\frac{1}{2}(d^{(i)})^{2}\big{]}$, where $m^{(i)}$ is the mass of this particle, and the overall coefficient of proportionality determined from the condition that $\sum_{i=1}^{N_{\mathrm{body}}}p^{(i)}=1$. Then the probability of association with the LMC is simply the sum of $p^{(i)}$ over particles that originally belonged to the LMC. One can obtain a representative sample of past orbits by selecting particles in proportion to their likelihood of matching and retrieving their trajectories from the original simulation. This method is the most direct, but may suffer from low-number statistics even in the highest-resolution snapshots, if the number of particles with sufficiently high matching probabilities is too small. In these cases, we artificially increase the uncertainties on kinematic measurements ($\mu_{l,b}$ and $V$) until the highest matching probability is below 0.1 (i.e., there are at least 10 different matches).

In the second method, we limit the matching to the spatial dimensions ($l$, $b$, $D$), obtaining a much larger sample (typically a few thousand particles) that are close to the target object in space. We select particles with sufficiently small Mahalanobis distance from the target object in the three spatial dimensions ($d_{\mathrm{l,b,D}}^{(i)}<4$), and their contribution to the distribution of each population $g$ (LMC or Milky Way) is weighted by the probability of spatial matching $p_{\mathrm{l,b,D}}^{(i)}\equiv\exp\big{[}-\frac{1}{2}(d_{l,b,D}^{(i)})^{2}\big{]}$, so that the total mass of particles in each population is $M_{g}=\sum_{i\in g}m^{(i)}\,p_{\mathrm{l,b,D}}^{(i)}$. Then we construct two Gaussian mixture models, one for the Milky Way particles and the other for the particles originating in the Magellanic system. Each population is described by $N_{\mathrm{comp}}=2$ Gaussian components in the entire 6d space of observables, with amplitudes $A_{g,c}$, means $\overline{\boldsymbol{X}}_{g,c}$ and covariance matrices $\mathsf{S}_{g,c}$, where the sum of amplitudes for each population $g$ is $M_{g}$. Since particle weights are pre-multiplied by the probability of spatial matching, the distributions in these three dimensions are nearly Gaussian by construction; nevertheless, this procedure still allows one to recover correlations between spatial and kinematic dimensions, which are sometimes important for the LMC debris. The probability distribution of the given population $g$ (LMC or Milky Way) in the entire 6d space is given by

To account for observational uncertainties, this probability distribution should be convolved with the Gaussian error distribution described by the error covariance matrix $\mathsf{E}$, which for a Gaussian mixture model amounts to adding $\mathsf{E}$ to $\mathsf{S}_{g,c}$ in the above equation. Here the uncertainties in the spatial dimensions should be to a very large number (formally, infinity), since we have already preselected particles according to their spatial matching probability, so no further convolution is needed (this effectively amounts to marginalizing over these dimensions). Finally, the probability of the object to belong to $g$-th population is simply $p_{g}/\sum_{g}p_{g}$. The advantage of this method over the first one is that the smooth probability distributions are constructed from a much larger number of particles, but evaluated only in a possibly narrow range of kinematic coordinates, thus effectively interpolating them in a region where the actual number of matching particles could be rather small. The method is illustrated in Figure 9 for two galaxies, Fornax and Leo I.

The third method does not use the $N$-body snapshot directly, but instead relies on orbit integrations in the evolving gravitational potential described by a series of multipole expansions, as detailed in Section 2.3. Namely, we sample $N_{\mathrm{samp}}=\mathcal{O}(10^{3})$ points from the uncertainty distribution described by the mean coordinates $\boldsymbol{X}$ and the error covariance matrix $\mathsf{E}$. Then we integrate the orbits back in time from these initial conditions to the starting point of the simulation (10–11 Gyr ago), when the LMC was still approaching the Milky Way from a distance of a few hundred kpc. At this point, one may be tempted to simply count the fraction of orbits that end up being bound to the LMC at that time, as was done in most previous studies, but this will not give a correct answer, because not all points in our initial conditions are equally probable a posteriori. To illustrate the fallacy, consider the case of very large PM uncertainties (e.g., Pisces II) and/or missing line-of-sight velocity $V$ (e.g., Delve 2). Most of the rewound orbits would miss the LMC by a large margin, but they would also be unbound from the Milky Way or have apocentres well beyond the virial radius, making these orbits unlikely to appear in our simulation. Even if the uncertainties are small, the unweighted fraction of LMC-bound orbits tends to be underestimated, sometimes quite severely. Instead, one needs to weigh each orbit with the probability of finding it among the initial conditions of our simulations, which consist of two equilibrium galaxy models. So we again use a mixture model approach, but instead of Gaussians in the present-day observable coordinates, we use the distribution functions of initial models $f_{g}(\boldsymbol{x},\boldsymbol{v})$ normalized such that the integral of $f_{g}$ over the entire 6d phase space is the total mass of each galaxy (or just its dark halo, in the case of Milky Way). Of course, in this evaluation we need to consider relative positions and velocities of particles with respect to each galaxy’s centre $\{\boldsymbol{x},\boldsymbol{v}\}_{g}^{\mathrm{(cen)}}$, and set $f_{g}$ to zero if the relative velocity exceeds the escape velocity. The un-normalized likelihood of association with either population is $p_{g}=\sum_{s=1}^{N_{\mathrm{samp}}}f_{g}\big{(}\boldsymbol{x}^{(s)}-\boldsymbol{x}_{g}^{\mathrm{cen}},\boldsymbol{v}^{(s)}-\boldsymbol{v}_{g}^{\mathrm{(cen)}}\big{)}$, and the normalized probability of association is again $p_{g}/\sum_{g}p_{g}$. This approach is similar to the one used in Correa Magnus & Vasiliev (2022), except that here we have a mixture of two populations bound to each galaxy, instead of just the Milky Way halo plus the “unmixed” population as in that paper. Since relatively few objects (mostly distant dSph) were found to be unmixed, and none of them are likely to be associated with the LMC, we ignore this component in the present study.

Reassuringly, all three methods agree remarkably well (to within 10% difference in probability) for the vast majority of objects, so in the summary Table 3 we provide a single value (average between three methods) for each galaxy in each of the four simulations, skipping objects with less than 10% probability of Magellanic association. These objects can be grouped into two classes. Eleven of them are likely to be Magellanic satellites at present time: Carina II, Carina III, Delve 2, Eridanus III, Horologium I, Horologium II, Hydrus I, Phoenix II, Pictor II, Reticulum II and SMC. The probabilities exceed 90% (except Reticulum II, which is only $\sim$80%), despite three of these galaxies lacking line-of-sight velocity measurements (Delve 2, Eridanus III and Pictor II). This classification broadly agrees with most recent studies (e.g., Fritz et al. 2019, Erkal & Belokurov 2020, Patel et al. 2020, Battaglia et al. 2022, Correa Magnus & Vasiliev 2022, Pace et al. 2022, see Table 1 in V23 for a compilation). Figure 10 shows that these galaxies, shown by blue colour, are indeed located within the region occupied by particles still bound to the LMC (brown dashed contour) and move in the same direction.

Perhaps more interestingly, many of the remaining galaxies in Table 3, shown by shades of red in the above figure, have a significant probability of association with the LMC, even though they are not currently bound to it. In particular, four of the classical dSph – Carina, Draco, Fornax and Ursa Minor – have a roughly equal chance to either be accreted with the Magellanic system or originate from elsewhere (in case of Carina, the chances are even higher, $\sim 2:1$). A distant galaxy Leo II has a probability of Magellanic association between 0.3 and 0.5, and even the peculiar Leo I, the least bound object in this sample, has a non-negligible chance of association in the $\mathcal{L}$3 series of models⁴⁴4We remind that regardless of whether Leo I is associated with the LMC or not, its current velocity relative to the Milky Way centre has been increased by a few tens km s^-1 due to the reflex motion of the Milky Way caused by the recent LMC flyby (Erkal et al., 2020; Correa Magnus & Vasiliev, 2022).. Even some galaxies with missing velocity components might be associated with the LMC under certain conditions; Table 4 lists the predicted PM and line-of-sight velocity values assuming such an association. Summing up the probabilities of galaxies in this group, we may conjecture that the original Magellanic system might have brought in some 4–6 galaxies that became Milky Way satellites, half of them in the mass range of classical dSph. They would populate the unusual 12 magnitudes luminosity gap between the SMC and other ultrafaint LMC satellites highlighted by Dooley et al. (2017) and Pardy et al. (2020), and indeed the latter study suggested that Carina and/or Fornax might fill this gap.

Figure 11 shows the distribution of Magellanic and Milky Way particles in different kinematic spaces: orbital pole orientation (top panel), energy–angular momentum (middle panel) and Galactocentric radius–velocity (bottom panel), and the locations of observed galaxies in these spaces. Particles still bound to the LMC and galaxies from the first group are well localized in each of these spaces, but possible former LMC satellites are considerably scattered and mixed with native Milky Way satellites in any particular projection, justifying the need for a full 6d analysis as in the present study. The distribution of former satellites in orbital phase is quite asymmetric, with all but one object (Pisces II) populating the leading arm of the tidal debris stream (the trailing arm is visible as the red band in the lower part of the bottom panel). A similar situation occurs with the globular clusters stripped from the Sagittarius galaxy, and the reason for this asymmetry is currently unknown. As the census of Milky Way satellites is very likely still incomplete (Koposov et al., 2009; Jethwa et al., 2016), one may use these “finding charts” to assess the possibility of Magellanic association of any object discovered in the future, though only if we have kinematic information. Figure 12 instead shows just the observable sky-plane distribution of LMC debris in different distance bins, both for first-passage (left) and second-passage (right) simulations. Unfortunately, the difference between these scenarios is not very dramatic, and mostly seen in the outermost distance bin and in the Northern hemisphere, which cannot be reached by debris from a single passage.

Figure 13 shows possible past orbits of selected galaxies in one of the simulations ($\mathcal{L}$3–$\mathcal{M}$11), illustrating that the galaxies from the second group have been stripped after the first pericentre passage, and their subsequent orbits do not resemble that of the LMC, but still stay roughly in the same orbital plane. This kinematic peculiarity, of course, has been known for a long time, as discussed in Section 4.1, but here we provide an explicit demonstration that such an arrangement is a natural outcome of the tidal stripping of a massive Magellanic system.

Pawlowski et al. (2011) proposed that the disc of satellites consists of tidal dwarf galaxies formed out of gas that was tidally stripped from the progenitor galaxy (possibly the proto-Magellanic system) during its pericentre passage. In our view, this scenario is conceptually similar to the one considered in the present paper, except that we assume that these galaxies were already formed before the first pericentre passage of the Magellanic system, but the stripping mechanism is the same, so we remain agnostic as to whether these galaxies are expected to be “primordial” (formed in the standard cosmological paradigm) or “second generation” (formed out of tidally stripped gas without necessarily be embedded in dark matter haloes). Full hydrodynamical simulations may be needed to resolve this question.

Although thin satellite planes appear to be quite rare around generic Milky Way-mass host galaxies in cosmological simulations (see Section 2.2.3 in Pawlowski 2021 for a review of relevant literature), Samuel et al. (2021) found that the fraction of satellite planes significantly increases if one selects host galaxies that have an LMC analogue near its pericentre. Garavito-Camargo et al. (2021b) suggested that the dynamical perturbation from the LMC may enhance the apparent clustering of orbital poles due to direct deflection of orbits by the LMC flyby and the global displacement of the Milky Way-centred reference frame relative to the outer halo. On the other hand, Pawlowski et al. (2022) and Correa Magnus & Vasiliev (2022) found that this effect alone is insufficient to account for the observed VPOS thinness. After excluding the current LMC satellites and rewinding the orbits of remaining objects in the combined Milky Way–LMC potential back in time to undo the LMC perturbation, the clustering of orbital poles was not significantly affected. To complement this analysis, Figure 14 shows the asymmetries in the spatial distribution (axis ratio of the moment of inertia tensor) and in kinematics (mean angle $\theta$ between the angular momentum of particles relative to that of the LMC) of particles in our simulations, separately for the Milky Way (black) and Magellanic particles (red). It is clear that the perturbations in the native Milky Way halo population are minimal – the orbital poles of particles stay very close to isotropic ($\theta\approx 90^{\circ}$) within 100 kpc and only show a slight enhancement in the direction of the LMC orbital pole further out, confirming the results of Pawlowski et al. (2022). By contrast, particles belonging to the Magellanic debris have a prominent concentration of their orbital poles ($\theta\lesssim 40^{\circ}$) and are more flattened (axis ratio 1:2), although the observed satellite plane is still thinner.

However, these idealised simulations may not be representative of the evolution of Galactic satellite distribution in a realistic cosmological context. Following upon the work of Samuel et al. (2021), Garavito-Camargo et al. (2023) considered Milky Way–LMC analogues in the Latte cosmological simulations, and found that the enhancement of orbital pole concentration is more prominent than in isolated simulations of Milky Way–LMC encounter, although the displacement of the inner Galaxy w.r.t. the outer halo remains the driving mechanism. However, this enhancement persist only for a short time around the pericentre passage of the LMC analogue: using the IllustrisTNG cosmological simulations, Kanehisa et al. (2023) found that such encounters do not produce long-lived satellite planes. This result is at odds with our conclusion that the objects stripped on the first passage remain close to the LMC orbital plane; the discrepancy might be caused by the idealised (non-cosmological) nature of our simulations.

Figure 15 shows the possible origins of galaxies that were part of the Magellanic system, expressed as the percentile of their initial binding energy among all LMC particles. Naturally, the eleven galaxies from the first group (current LMC satellites) must have started rather deep inside its potential in order to avoid being stripped on the first passage, but not too tightly bound (with the exception of Pictor II) to avoid tidal disruption by the LMC itself. Galaxies from the second group have generally broader distributions of possible initial binding energies, but still avoiding least bound orbits, which may explain their stronger spatial and dynamical coherence compared to the entire distribution of Magellanic debris shown on the previous figure.

We also examined the probability of Magellanic association of dSphs for the simulation of a single passage ($\mathcal{L}$2$\mathcal{M}$10 run for 4 Gyr). In this case, all eleven galaxies from the first group are still highly likely to be current LMC satellites, but among the second group, only Carina retains a significant (40–50%) probability of association, Grus II, Hydra II, Pisces II and Tucana II stay at the level 20–40%, and remaining galaxies are not compatible with the distribution of Magellanic debris. Thus the ability to explain the existence of the satellite plane is an important argument in favour of the second-approach model in the standard cosmological paradigm. The alternative scenario advocated by Pawlowski & Kroupa (2013) and subsequent studies involves modified Newtonian gravity theories (e.g., MOND, Milgrom 1983), in which case it is the early encounter with M31 rather than LMC that could give rise to the dynamically coherent second-generation tidal dwarf galaxy population (Bílek et al., 2018; Banik et al., 2022).