The Milky Way’s plane of satellites is consistent with ΛCDM

The Milky Way is surrounded by 11 ‘classical’ satellite galaxies in a remarkable configuration: a thin plane that is possibly rotationally supported. Such a structure is thought to be highly unlikely to arise in the standard (ΛCDM) cosmological model (Λ cold dark matter model, where Λ is the cosmological constant). While other apparent discrepancies between predictions and observations of Milky Way satellite galaxies may be explained either through baryonic effects or by invoking alternative forms of dark matter particles, there is no known mechanism for making rotating satellite planes within the dispersion-supported dark matter haloes predicted to surround galaxies such as the Milky Way. This is the so-called ‘plane of satellites problem’, which challenges not only the ΛCDM model but the entire concept of dark matter. Here we show that the reportedly exceptional anisotropy of the Milky Way satellites is explained, in large part, by their lopsided radial distribution combined with the temporary conjunction of the two most distant satellites, Leo I and Leo II. Using Gaia proper motions, we show that the orbital pole alignment is much more common than previously reported, and reveal the plane of satellites to be transient rather than rotationally supported. Comparing with new simulations, where such short-lived planes are common, we find the Milky Way satellites to be compatible with standard model expectations. The ‘plane of satellite galaxies’ surrounding our Milky Way seemed to defy dark matter theory for 40 years. Observations now suggest that the alignment is transient, while new simulations form similar structures far more often than previously thought.

The Milky Way is surrounded by 11 'classical' satellite galaxies in a remarkable configuration: a thin plane that is possibly rotationally supported. Such a structure is thought to be highly unlikely to arise in the standard (ΛCDM) cosmological model (Λ cold dark matter model, where Λ is the cosmological constant). While other apparent discrepancies between predictions and observations of Milky Way satellite galaxies may be explained either through baryonic effects or by invoking alternative forms of dark matter particles, there is no known mechanism for making rotating satellite planes within the dispersion-supported dark matter haloes predicted to surround galaxies such as the Milky Way. This is the so-called 'plane of satellites problem', which challenges not only the ΛCDM model but the entire concept of dark matter. Here we show that the reportedly exceptional anisotropy of the Milky Way satellites is explained, in large part, by their lopsided radial distribution combined with the temporary conjunction of the two most distant satellites, Leo I and Leo II. Using Gaia proper motions, we show that the orbital pole alignment is much more common than previously reported, and reveal the plane of satellites to be transient rather than rotationally supported. Comparing with new simulations, where such short-lived planes are common, we find the Milky Way satellites to be compatible with standard model expectations.
The original discovery of the Milky Way's 'plane of satellites' (then just five galaxies) 1 preceded the advent of the Λ cold dark matter (ΛCDM) model, where Λ is the cosmological constant, as the paradigm for galaxy and structure formation 2 . A key ΛCDM prediction is that galaxies such as the Milky Way (MW) are surrounded by a dispersion-supported dark matter halo and by satellite galaxies formed within its substructures. However, while several ΛCDM predictions, including the discovery of dozens of additional MW satellites, have since been borne out 3 , the 'plane of satellites problem' 4,5 has emerged as its most persistent challenge [6][7][8][9][10][11] .
That the 'plane of satellites' problem has so far eluded resolution is not for lack of trying. Planes of satellites form with the same (low) frequency in collisionless and hydrodynamic cosmological simulations [12][13][14] and in MW analogues in isolation or in pairs 15,16 , with no significant correlation with other properties of the host halo 16 . There is evidence that filamentary accretion 17,18 , a compact satellite system 19 or the presence of massive satellites 20 can generate some anisotropy, but systems as thin as the Milky Way's are still very rare 12 . Moreover, any planes that do form in ΛCDM are transient, chance alignments of substructures 14,[20][21][22] rather than long-lived, rotationally supported disks. Article https://doi.org/10.1038/s41550-022-01856-z angular component, independently of the radial distribution. For both c/a and (c/a) red , smaller values imply greater anisotropy; for example, a sphere has c/a = 1 while a perfectly thin disk has c/a = 0. Note that for small N, the expectation values of c/a and (c/a) red decrease, regardless of the underlying anisotropy 25 . Figure 1 shows the present positions and estimated orbits of the 11 brightest MW satellites projected along the principal axes. We measure c/a = 0.183 ± 0.004 and (c/a) red = 0.3676 ± 0.0004, mean ± s.d.

The effect of the radial distribution
Earlier work 26 has found that only 0.7% of ΛCDM simulations produce systems as anisotropic as the MW. However, we find this likely to be a severe underestimate caused by the artificial disruption of satellites in numerical simulations, which can cause artificially extended radial profiles [27][28][29][30] . To identify analogues of the classical MW satellites in simulations, it is common to select at z = 0 the 11 satellites with the highest value of v peak , the peak value of the maximum circular velocity, v max , across a halo's history 17 . In our own analysis of 202 MW analogues resimulated in the SIBELIUS constrained simulations intended to reproduce the structures observed in the Local Universe 31 , including only surviving subhaloes, we find only 2 (1%) systems with c/a as low as the MW, and none that reproduce the radial distribution. However, this is an 'incomplete' sample. Accounting for artificial disruption 32 and selecting the 11 satellites among a 'complete' sample containing both the We examine here the contention that the MW contains an exceptional plane of satellites, explain the origin of the observed anisotropy and study its time evolution in light of proper motion measurements by the Gaia space telescope 23 .

The present MW plane of satellites
The Milky Way's 'plane of satellites' canonically consists of the 11 'classical' satellites, the brightest within a radius of r = 300 kpc of the Galactic Centre, believed to constitute a complete sample. To characterize the spatial anisotropy of a system of N satellites, it is customary to consider the inertia tensor, defined as where x n are the coordinates of the nth satellite relative to the centre of positions, and i and j index the three spatial dimensions. We label the square roots of its eigenvalues as a, b and c, corresponding to the dispersions in position along the unit eigenvectors x a , x b and x c . A related metric is the 'reduced' inertia tensor 24 , defined after projection of the positions onto a unit sphere. We label the square roots of its eigenvalues as a red , b red and c red . Both c/a and (c/a) red ≡ c red /a red parametrize the spatial anisotropy but differ in the weight attached to each galaxy. c/a measures the full spatial anisotropy, whereas (c/a) red considers only its x a (kpc) x a (kpc) x c (kpc)  show the orbits of the 4 satellites with Galactocentric distances beyond 100 kpc and the right panels show the remaining 7 orbits in an inset of ±100 kpc around the Galactic Centre. With the Gaia EDR3 measurements, the proper motions are very well constrained, with the exception of the LMC and SMC. It can also be seen that the MW satellites are highly concentrated, with 7 out of 11 within 100 kpc and only 2, Leo I and Leo II, at r > 200 kpc. Several galaxies, including Leo I and II, are presently crossing the 'plane' (indicated by the grey horizontal lines in the bottom panels), which soon disperses as a result.    surviving and the recovered, artificially disrupted satellites, we find radial distributions resembling the MW's, as shown in the left panel of Fig. 2. This reproduces the results of very-high-resolution simulations 30,33 . As each satellite contributes to the inertia (equation (1)) in proportion to r 2 i , c/a is sensitive to the radial profile. To quantify this relationship, we introduce the Gini coefficient formalism. The central panel of Fig. 2 shows the summed weights of the closest i satellites from the centre, , normalized by the total weight of all 11 satellites, The area between each curve and the diagonal measures the inequality of the satellites' contributions to the inertia, or the sample Gini coefficient of inertia, Compared to a more equal distribution, the Milky Way's centrally concentrated radial profile is equivalent to sampling a system with fewer points, lowering the expectation value of c/a.
The right panel of Fig. 2 shows the relationship between G and c/a. Systems with higher central concentration (higher G) tend to be more anisotropic (lower c/a). Accounting for artificial disruption (filled circles), 58% of ΛCDM systems have G above the MW and 11 (~5%) have c/a < 0.183. Neglecting this effect (faint crosses) produces no systems with G as high the MW and only two (1%) with as low c/a.
The Milky Way's anisotropy additionally results from the fact that its two outermost satellites, Leo I and Leo II, which contribute two thirds of the total inertia, are currently in close proximity. However, as is already apparent from Fig. 1, this constellation is short-lived.

The clustering of orbital poles
Supporting the notion that the satellite plane constitutes a spinning disk 34 , the orbital poles of 7 of the 11 classical satellites (the Large Magellanic Cloud (LMC), the Small Magellanic Cloud (SMC), Fornax, Leo II, Carina, Draco and Ursa Minor) are reported to cluster with a standard deviation in direction of only 16.0° (ref. 26 ), found in only 0.04% of ΛCDM systems. Using the more precise proper motions from Gaia Early Data Release 3 (EDR3) for the same satellites, we find that this angle increases to 23.2 ∘ +3.5 −2.8 . We also repeated the analysis, and find that a different subset (that includes Leo I instead of Leo II) yields a smaller dispersion of 18.9 ∘ +1.9 −1.4 . Among our sample of 202 simulated systems, adopting either 18.9° or 23.2° and accounting for a minimum look-elsewhere effect, we find three or five systems with subsets of satellites with a smaller probability to arise from isotropic distributions. That is, we find that ~2% of ΛCDM systems contain satellites whose orbital poles are even more anisotropic than the most clustered subset of the Milky Way, an ~50-fold increase over previous results. The orbital clustering of a subset of the Milky Way satellites is unusual in ΛCDM, but not astronomically improbable. Notably, one of our systems has both a c/a value below the MW value (0.166 compared to 0.183 for MW) and an orbital pole clustering more unlikely to arise out of isotropy than the MW's (probability of 0.005 compared to 0.009 for the MW). With no strong correlation between orbital pole clustering and anisotropy, and estimated probabilities below 10 −3 and 10 −2 , respectively 26 , the combination of the two had previously been considered extremely unlikely.
Importantly, while the 'plane of satellites' includes all 11 classical satellites, the orbital anisotropy only concerns a subset, which is in fact more spatially isotropic than the system as a whole. The orbital pole clustering does not drive the spatial anisotropy.

The plane of satellites is transient
Another defining feature of a rotationally supported disk would be a significantly higher velocity dispersion parallel to the plane (σ v∥ ) than perpendicular to it (σ v⊥ ). However, for the MW's classical satellites, we measure σ v∥ = 165.1 ± 1.2 km s −1 and σ v⊥ = 121.6 ± 0.4 km s −1 . The ratio, σ v∥ /σ v⊥ = 1.36, is indistinguishable from the purely geometrical factor of √2. By this basic measure, the plane is not rotation supported.
The longevity of the plane can also be tested directly via orbital integration. This method was first, but inconclusively, applied using pre-Gaia data 35 . In this work, we benefit from significantly more precise observations, including Gaia EDR3 proper motions and more accurate distances.
In Fig. 1, we saw that several satellites are presently crossing the plane, while Leo I and II, which dominate the inertia, are moving apart.
To elucidate the impact of such 'fortuitous alignments' on the anisotropy, we show in Fig. 3 the effect on c/a and (c/a) red of randomizing the position of each galaxy on the sky at the observed radius, keeping all other galaxies at their observed positions. For Sagittarius, the satellite with the smallest distance, the c/a distribution is extremely narrow: Sagittarius contributes less than 1% to the inertia tensor. However, randomizing the angular position of either of the two outer satellites, Leo I or Leo II, raises the median value of c/a to 0.28 and 0.31, respectively, with maxima of 0.53 and 0.63, far above the ΛCDM median.
While placing satellites at random sky coordinates highlights the sensitivity of c/a to individual galaxies, it is not a physical process. However, we also show at the top of Fig. 4 the anisotropy distributions when each satellite simply moves along its orbit. The time-averaged anisotropy of the system is then calculated over one full orbital period centred on the present time. Depending on the orbital phase of Leo II alone, c/a could be as high as 0.39, more isotropic than most ΛCDM systems. In other words, almost the entire anisotropy of the classical MW satellites is due to the orbital phases of Leo I and Leo II. Simply omitting both Leo I and Leo II from the analysis would yield c/a = 0.279, also more isotropic than 37% of ΛCDM systems of 9 satellites.
The bottom panels of Fig. 4 show time-averaged probability densities of c/a and (c/a) red when all satellites evolve simultaneously. The current value of c/a is significantly lower than in the recent past: over the past 0.5 and 1 Gyr, the time-averaged medians of c/a are 0.23 and 0.27 respectively, greater than 13% and 23% of ΛCDM systems. (c/a) red has varied widely. Neither metric is an invariant of the satellite system, both are sensitive to the orbital phases.
The four panels of Fig. 5 show the evolution of c/a, (c/a) red and of the orientations of the planes defined by the full and reduced inertia tensors, which we parametrize by the angles between the vectors normal to the planes, x c and x c,red , and their present day equivalents, We further see that the orientation of the Milky Way's plane of satellites is not stable, but has tilted by ~17° over the past 0.5 Gyr (and ~40° for the reduced definition). This is readily understood considering that the plane is largely defined by only two outlying satellites. This is demonstrated clearly in Fig. 6, where we show the orbits and positions of the 11 classical satellites projected edge-on, in the frame defined by the minor and major axes of the inertia tensor at z = 0. The central panel, which shows the positions of the satellites at z = 0, shows the plane aligned with the major axis, analogous to the bottom left panel of Fig. 1. The other four panels show the evolution up to 1 Gyr in the past and the future. The orientation of the plane of satellites is changing, so that at each moment, it points toward the current position of the outermost satellites. Rather than satellites orbiting inside a stable plane, the plane tilts as it tracks the positions of its most distant members.

Discussion
The high reported anisotropy of the MW satellite system can largely be attributed to its high central concentration, not previously reproduced in simulations, combined with the close but fleeting contiguity of its two most distant members. Accounting for the radial distribution reveals the MW satellites to be consistent with ΛCDM expectations. Compared with previous work, we also find a much higher likelihood of subsets whose orbital poles are as clustered as the MW. Although the Milky Way contains such a subset, the plane of satellites does not constitute a rotationally supported disk. Instead, it evolves on  Our orbital integration assumes a static MW potential with satellites as massless test particles. We Monte Carlo sample all sources of observational error and also vary the components of the potential within the uncertainties (Methods). We find our results to be robust, and while the real potential is more complex (for example, due to the presence of the LMC), these simplifications are valid within a dynamical time (~2 Gyr) of the halo, particularly for the important outer satellites 36,37 . Furthermore, a more complex potential would only accelerate the dissolution of a planar configuration 38 . Our simulations do not include the potential disruption of satellites by the disk of the MW. This might slightly extend the radial distributions, but not enough to change our conclusions. Importantly, among systems that have radial distributions similar to the MW, flattened 'planes of satellites' are quite common.
This work only directly addresses the archetypal 'plane of satellites' around the MW. Anisotropic satellite distributions have also been reported around several other galaxies 14,39 with different criteria, which can exacerbate the look-elsewhere effect. Assessing their significance requires a careful statistical analysis 12,40 . While not all criteria are equally sensitive to the radial distribution, we also expect that the much higher anisotropy we report here for simulated MW analogues will apply to ΛCDM analogues of other, similarly defined systems.
After centuries of observations, the Milky Way and its satellites are the best studied system of galaxies in the Universe. Viewed in sufficient detail, every system inevitably reveals some remarkable features. However, based on the best currently available data, there is no evidence for a plane of satellites incompatible with, or even particularly remarkable in, ΛCDM. On the contrary, as judged by the spatial anisotropy of the brightest satellite galaxies, we appear to live in a fairly regular ΛCDM system.

Observations
We adopt the sky positions and radial velocities from the McConnachie 41 catalogue, and combine these, where available, with McConnachie and Venn 23 proper motion measurements based on Gaia EDR3 (ref. 42 ). The systemic proper motions were measured within a Bayesian framework that combines information from stars with full astrometric data with information from stars with only photometric and proper motion data. The method is a mixture model that associates a probability for each candidate star to be associated with the target galaxy taking into account foreground and background contaminants. It constitutes the best technique currently available in the literature to determine proper motions and provides the most precise measurements so far 37 .  Table 2.
We also repeated our analysis using the Gaia EDR3 proper motions of Battaglia et al. 37 and the Gaia DR2 proper motions described in Riley et al. 43 . A comparison of the evolution of c/a and (c/a) red , analogous to Fig. 5, is shown in Extended Data Fig. 1. Unsurprisingly, the evolution based on the two EDR3 datasets are in excellent agreement. The main difference when using the DR2 data is the larger uncertainty (the astrometry errors of EDR3 are reduced by approximately a factor of two compared with DR2), but the evolution of both c/a and (c/a) red is essentially the same in all three datasets, and the results are consistent within the respective errors.

Monte Carlo samples.
We account for measurement errors by generating Monte Carlo samples of the satellites in the space of observed quantities: distance modulus, radial velocity and proper motions, as well as the position of the Sun relative to the Galactic Centre. We model each observable as a Gaussian distribution with the mean value and standard deviation given by the measurements and their quoted errors. For the Sun's distance from the Galactic Centre we assume R ⊙ = (8.178 ± 0.022) kpc (ref. 51 ), for the circular velocity at the Sun's position V circ = (234.7 ± 1.7) km s −1 (ref. 52 ) and for Sun's motion with respect to the local standard of rest (U, V, W) = (11.10 ± 0.72, 12.24 ± 0.47, 7.25 ± 0.37) km s −1 (ref. 53 ), where U is defined positive toward the Galactic center, V is positive in the direction of Galactic rotation and W is positive toward the North Galactic Pole.
The clustering of orbital poles. To characterize the clustering of orbital poles, we adopt the orbital pole dispersion for a subset of N s satellites, Δ std , defined by Pawlowski and Kroupa 26 as: where θ i is the angle between the orbital pole of the ith satellite and the mean orbital pole of the satellites in the subset. To compute the clustering of an observed system relative to expectation, the same analysis is performed on the observations and simulations.  x c,0 (kpc) x a,0 (kpc) x a,0 (kpc) x a,0 (kpc) x a,0 (kpc) x a,0 (kpc) Article https://doi.org/10.1038/s41550-022-01856-z N s = 3...11, and discovered that N s = 7 yielded the most unusual configuration. However, there is no a priori reason to specifically consider N s = 7. When considering only a proper subset of satellites, the interpretation of Δ std (N s ) as evidence for unusual clustering is subject to the 'look-elsewhere effect'. To account for this, we follow here the method of ref. 12 , which involves performing the same analysis for the simulated systems. As in the observations, we consider all subsets of size N s = 3...11 in each simulated system, and identify the most unlikely to arise by chance from an isotropic distribution, which we calculate based on 10 5 isotropic distributions of N = 11 points, and the probability distributions of Δ std (N s ) for all N s = 3...11 possible subsets. In Extended Data Fig. 2, we show a Hammer equal area projection of the most likely orbital poles and their Monte Carlo sampled uncertainty. The 3 black circles identify the clustering of 7 of the 11 satellites. The dotted line corresponds to Δ std (7) = 16.0°, the dispersion calculated by Pawlowski and Kroupa using pre-Gaia EDR3 data. The solid line shows our results, Δ std (7) = 23.2°, for the same satellites using Gaia EDR3 data. The dashed line shows our result for the most clustered subset Δ std (7) = 18.9°, exchanging Leo II for Leo I.

Time evolution and time integration
To infer the time evolution of the Milky Way satellite system, the orbits of the satellites are integrated numerically as massless test particles in a static potential using the Gala package 54 . The potential consists of a disk, stellar nucleus and bulge, and a dark matter halo. The disk is modelled as an axisymmetric Miyamoto-Nagai disk 55 , which, for our default model, has disk mass 5.17 × 10 10 M ⊙ , a = 3 kpc and b = 0.028 kpc (ref. 56 ). The nucleus and stellar bulge are both modelled as spherically symmetric Hernquist profiles 57 . For the nucleus we assume a mass of 1.71 × 10 9 M ⊙ and a scale radius a = 0.07 kpc, and for the bulge we assume a mass of 5.0 × 10 9 M ⊙ and a = 1.0 kpc. For the dark matter halo, we assume a spherically symmetric Navarro, Frenk and White (NFW) 58 potential.
Until recently, the Milky Way halo mass may have been a prohibitive source of uncertainty for calculating the orbital evolution of the satellites, as its value was known only to within a factor of two 59 . However, the Galactic halo mass has now been estimated with an uncertainty of only about 20% using Gaia data. Multiple dynamical probes, such as the stellar rotation curve, the escape velocity, the motions of halo stars, globular clusters and satellite galaxies 60-64 , consistently imply a dark matter halo mass for the MW of M 200 = (1.0 ± 0.2) × 10 12 M ⊙ and NFW concentration, c 200 = 11 ± 2.
Based on these results, we adopt a reference MW halo of mass 1.0 × 10 12 M ⊙ and a concentration parameter c 200 = 11, corresponding to an NFW scale radius of r s = 19.2 kpc. The positions and velocities relative to the plane of satellites, and the orbital periods and apocentre distances for the default potential, are listed in Supplementary Table 1, where the quoted uncertainties reflect 68% confidence intervals for all quantities based on Monte Carlo sampling. Varying the MW potential within the observational uncertainties does not significantly affect the conclusions of our study, as we show in Extended Data Fig. 3 where we vary the halo mass between 0.8 and 1.2 × 10 12 M ⊙ . The inferred anisotropies, as measured by c/a and (c/a) red , begin to diverge for lookback times beyond ~500 Myr, but the behaviour is qualitatively similar across the halo mass range. As shown in Extended Data Fig. 4, we also tested the impact of varying the other three mass components and the concentration parameter within the observational uncertainties, and found no significant effects. The true Milky Way potential evolves with time, but the dynamical time of the halo (~2 Gyr at z = 0) is much longer than the timescale for the reported dissolution of the plane of satellites (several hundred Myr). A further possible source of uncertainty may be the mass of the LMC and its perturbation of the potential. However, at a distance of 49 kpc from the centre, even a massive LMC does not significantly perturb the acceleration field at distances above 150 kpc (ref. 36 ) and would not significantly affect the orbits of the two most distant satellites, Leo I and Leo II.

Numerical simulations
Initial conditions. The simulations used in this work are cosmological zoom-in-constrained simulations, based on initial conditions created for the SIBELIUS project 31 and designed to reproduce Local Group (LG) analogues within the observed large-scale structure. The simulations assume a ΛCDM cosmology with density parameters Ω 0 = 0.307 and Ω Λ = 0.693 for matter and dark energy, respectively, an r.m.s. density variation on a scale of 8h −1 Mpc of σ 8 = 0.8288 and a Hubble parameter of h = 0.6777. We use physical units throughout this work. Building on an octree representation of the phase information 65 , we used the methods described in ref. 66 to supplement the observationally constrained scales in the initial density field 67,68 with independent random information below 3.2 comoving Mpc (cMpc). In total, we generated 60,000 simulations, resulting in several thousand loosely defined Local Group analogues. From these, we selected 112 for the high-resolution resimulations used in this work. All initial conditions refine a Lagrangian region extending to at least r = 3 Mpc around the centre of the LG at z = 0, with the remainder of a 1,000 3 cMpc 3 volume populated with progressively more massive particles. The simulations start at z = 127. More details about the initial conditions may be found in ref. 31 .
Simulations. All simulations used in this work were performed with the public version of the Gadget-4 code 69 , on the Cosmology Machine at the University of Durham and at the Finnish IT Center for Science. Extended Data Fig. 5 shows a comparison of one LG analogue at two different resolutions. The left panel corresponds to the resolution of the set of 60,000 simulations of the SIBELIUS project (particle mass 2.0 × 10 9 M ⊙ ). At this resolution, a MW analogue halo contains approximately 500 particles, and only the largest substructures are identifiable. The right panel shows a simulation with a particle mass of 1.0 × 10 6 M ⊙ . At this resolution, a MW analogue halo contains approximately 10 6 particles, and an average of ~200 subhaloes down to 2 × 10 7 M ⊙ can be identified within 300 kpc from the centre. All results presented in this paper are based on simulations at this resolution.

Structure finding.
Structures and self-bound substructures were identified using the Friends-of-Friends and Subfind algorithms implemented in Gadget-4 at 60 snapshots equally spaced in time, from z = 4 until a lookback time of 1 Gyr, and a further 40 snapshots equally spaced over the final 1 Gyr up to z = 0. Given their mass and separation, the two most massive self-bound substructures of the LG analogues can either belong to the same or to separate Friends-of-Friends structures. Throughout this work, we refer to the two principal self-bound substructures of each LG analogue at z = 0 simply as 'haloes' and to the lower mass substructures within 300 kpc of the centre of potential of each halo as 'satellites'. We select Local Group analogues as pairs of haloes with individual masses in the range 0.6-2.2 × 10 12 M ⊙ , separated by 500-1,000 kpc, with radial velocity −150 < v r <−50 km s −1 and transverse velocity v t < 70 km s −1 . In total, our set of high-resolution simulations contains 101 LG analogues, and, for the purposes of this work, we consider both haloes as a MW analogue.
We use Gadget's on-the-fly merger tree construction to find the progenitors of these subhaloes at previous snapshots. We cut the chain of links when a subhalo's progenitor is no longer found, or when a clear discontinuity in mass and position indicates that a satellite's progenitor has been erroneously identified as the main halo. At each snapshot, we record the maximum circular velocity of each subhalo, x c,red and define v peak as the highest value of v max of a subhalo and its progenitors over time. Following ref. 17 , we use the standard procedure to rank satellites by v peak , and identify the top 11 within 300 kpc of each MW analogue at z = 0 as analogues to the classical MW satellites.
Obtaining a complete satellite sample. As noted above, the radial distribution of satellites is important for the anisotropy. Numerical simulations suffer from the artificial disruption of substructures, that Article https://doi.org/10.1038/s41550-022-01856-z can affect subhaloes far beyond the particle number limit at which they can theoretically be identified 28,70 . According to van den Bosch and Ogiya 71 , the main cause of this artificial tidal disruption is inaccurate force softening, which can cause force errors that, for particles within substructures, do not cancel out: once a particle is lost to the main halo, it cannot be recovered. Additionally, the amplification of discrete noise in the presence of a strong tidal field near the centre of a halo can lead to a runaway instability that can lead to the complete, but purely numerical, disruption of a subhalo.
These effects can, however, be mitigated using semi-analytical models (which populate merger trees constructed from simulated dark matter subhaloes with galaxies). These models include so-called 'orphan' galaxies, that is, galaxies whose dark matter subhalo has been numerically disrupted. After the subhalo is disrupted numerically, its subsequent evolution is followed by tracing the positions of its most bound particle 32 . Our 'complete' sample includes these 'orphan' subhaloes.
One important result of this work is that the 'incomplete' and 'complete' samples of satellite haloes have different radial distributions. Even though our high-resolution simulations resolve, on average, 200 surviving satellite haloes inside 300 kpc of each MW analogue at z = 0 and we rank the satellites by v peak (v max being more strongly affected by tidal stripping), we find that the radial distribution of the top 11 surviving satellites in the 'incomplete' samples are systematically and significantly less centrally concentrated than the MW's brightest satellites. Figure 2 in the main text showed consistency between the radial distributions of the complete sample and the Milky Way satellites. It is reproduced in the top row of Extended Data Fig. 6. By contrast, the bottom row of Extended Data Fig. 6 shows how our comparison between simulations and observations would have looked had we considered only subhaloes from the incomplete samples. In the bottom left, it can be seen that their radial distribution is systematically offset from the MW data (shown as thick black line). For example, the innermost nine MW satellites are found within a distance of 140 kpc from the Galactic Centre, but none of the 202 incomplete samples has 9 satellites within this radius. The bottom-centre panel shows that the more uniform radial distributions of the incomplete samples lead to more equal contributions to the moment of inertia than in the case of the MW satellites. In fact, as shown on the bottom-right panel, none of the incomplete satellite systems have Gini coefficients as high as the MW's and only two have c/a as low or lower. This is in line with previous studies (for example, refs. 22,26 ), which, presumably using incomplete subhalo populations, have found c/a values as low as the MW's to be very rare.
To demonstrate the dependence of the anisotropy on the Gini coefficient and the radial concentration independently of our simulations, we also repeat the analysis on synthetic random data. The distributions shown in the top row of Extended Data Fig. 7 use a parent radial distribution that is uniform in r 1/2 , while those in the bottom row use one that is uniform in r; both parent angular distributions are isotropic. The relation between c/a and G is independent of the radial distribution, but the more centrally concentrated r 1/2 distribution attains larger values of G and higher anisotropies than the more extended one.
Baryonic effects. In hydrodynamic simulations, Milky Way mass haloes tend to be slightly more spherical than their counterparts in collisionless simulations 72 , which might reduce the anisotropy of subhaloes. Additionally, the adiabatic contraction induced by baryons can increase the concentration of the dark matter halo 73 , which can lead to a more central concentration of subhaloes but, conversely, could also enhance the disruption of subhaloes near the centre.
The presence of a stellar disk could also lead to additional disruption of subhaloes close to the centre of the galaxy [74][75][76][77] . To estimate its possible impact, we have calculated the fraction of subhaloes and orphans that pass through a cylinder whose radius, r = 4.3 kpc, is twice the scale-length and whose height, h = 0.5 kpc, is twice the scale-height of the MW disk during the past 1.3 Gyr. We find that in approximately half of the systems, none of the 11 satellites that we consider analogues of the classical satellites have passed through this disk, and at most one satellite has in ~83% of systems. Since satellites that passed through the disk are more likely to be found near the centre, this disruption could make the radial profiles slightly more extended. We find that this does not have a significant effect on our results: removing all satellites identified as having passed through the disk and replacing them with the next highest by v peak only reduces the median Gini coefficient from 0.63 to 0.61, while it increases the median value of c/a from 0.35 to 0.36.
Some authors have argued that, in aggregate, the various baryonic effects that could individually lead to higher or lower anisotropy can increase the anisotropy of satellite systems in ΛCDM 13,35,78 . On the other hand, refs. 26 and 20 concluded that baryons only have a negligible effect. A fully realistic simulation of the Milky Way satellite system would clearly include baryons, but we believe that the inclusion of baryons would, at the very least, not decrease the anisotropy significantly. In other words, including baryons is unlikely to make the 'plane of satellite problem' worse.

Data availability
The observational data that we use are listed with references in Supplementary Table 2. Both the observational data and the simulation data necessary to recreate all the figures and tables in the paper are available at https://github.com/TillSawala/plane-of-satellites. The entirety of the raw simulation data, comprising 23 TB, are archived at the DiRAC Data Centric system at Durham University. Access may be provided by reasonable and specific request to the corresponding author.

Code availability
The analysis in this paper was performed in Python 3 and makes extensive use of open-source libraries, including Matplotlib 3.4.2, NumPy 1.21.1 (ref. 79 ), SciPy 1.7.0 (ref. 80 ), GalPy 1.7.0 (ref. 81 ), Py-SPHViewer 82 , TensorFlow 83 and Gala 1.4.1 (ref. 54 ). The complete analysis code is available at https://github.com/TillSawala/plane-of-satellites. indicates the dispersion that we find for the same set based on Gaia EDR3, the dashed black curve indicates the minimum dispersion that we find for seven satellites exchanging Leo I and Leo II. The orbital poles of the MW satellites are significantly clustered, but several of our simulated ΛCDM systems contain equally or more strongly clustered satellite systems.