A rotating satellite plane around Milky Way-like galaxy from the TNG50 simulation

We study the Satellite Plane Problem of the Milky Way\ (MW) by using the recently published simulation data of TNG50-1. Here, we only consider the satellite plane consisting of the brightest 14 MW satellites \ (11 classical satellites plus Canes Venatici I\ (CVn I), Crater II and Antlia II). One halo\ (haloID=395, at z=0, hereafter halo395 ) of 231 MW like candidates, possesses a satellite plane as spatially thin and kinematically coherent as the observed one has been found. Halo395 resembles the MW in a number of intriguing ways: it hosts a spiral central galaxy and its satellite plane is almost ($\sim 87^{\circ}$)perpendicular to the central stellar disk. In addition, halo395 is embedded in a sheet plane, with a void on the top and bottom, similar to the local environment of MW. More interestingly, we found that the 11 of 14 of the satellites on the plane of halo395, arise precisely from the peculiar geometry of its large-scale environment\ (e.g. sheet and voids). The remaining three members appeared at the right place with the right velocity by chance at z=0. Our results support previous studies wherein the Satellite Plane Problem is not seen as a serious challenge to the $\Lambda$CDM model and its formation is ascribed to the peculiarities of our environment.


INTRODUCTION
The accepted cosmological paradigm, known as the ΛCDM model, has been immensely successful at describing a myriad of observations.As a theory rooted in empirical observation and based on first principles, it can account for the formation of cosmic structures on a variety of scales: from galaxies like the Milky Way, to clusters and ultimately to the large scale structure of the Universe.Its success however has been challenged by observations of the smallest cosmological objects: dwarf galaxies.The careful observations of these objects in the vicinity of the Milky Way, when compared against theoretical models for how we expect them to form, have thrown up a number of problems or challenges (e.g.reviews Bullock & Boylan-Kolchin (2017); Sales et al. (2022) ).The Satellite Plane Problem may be the most stubborn one (Kroupa et al. 2005;Kroupa 2012;Bullock & Boylan-Kolchin 2017;Pawlowski 2018Pawlowski , 2021;;Boylan-Corresponding author: Yingzhong Xu, Xi Kang xuyingzhong@zju.edu.cn,kangxi@zju.edu.cnKolchin 2021; Sales et al. 2022;Sawala et al. 2022;Kanehisa et al. 2023).Although it was proposed decades ago (Kunkel & Demers 1976;Lynden-Bell 1976), its solutions are still being debated.
Despite the fact that such configurations exist rarely, many mechanisms have been proposed to understand the formation of the thin rotating planes.For example, Zentner et al. (2005); Libeskind et al. (2005); Shao et al. (2018) argue that satellites accreting along the filaments will naturally lead to planar configuration, while others (e.g.Pawlowski et al. (2012)) disagree.Li & Helmi (2008); Santos-Santos et al. (2020) suggest that accretion in a group may explain the kinematic coherence of the plane.Furthermore, recent studies (Samuel et al. 2021) have shown that the presence of Large Magellanic Cloud(LMC)-like satellites will increase the chances of forming planar structures, since many smaller satellites will fall together with LMC-like objects.Besides that, some studies (Bílek et al. 2018(Bílek et al. , 2021;;Banik et al. 2022) have proposed a more radical solution: satellite planes are generated from the tidal tails of interacting galaxies in the context of modified Newtonian dynamics instead of ΛCDM.
Previous studies on the small-scale challenges that use hydro-simulations, always have to make a choice: either a large cosmological box or zoom-in simulations.On one hand larger simulated volumes lack sufficient resolution to study the dwarf galaxies at the heart of the small scale challenges, while on the other hand smaller volumes necessarily have smaller sample sizes thus limiting the scope of such studies.It is also worth noting that the recent N-body simulations can guarantee both high resolution and large volume (e.g.Ishiyama et al. (2021)).However, we choose to use hydro-simulations in this work.Because they include baryonic processes, we can try to find a more realistic MW-like system: a thin rotating satellite plane around and perpendicular to the central spiral stellar disk.
The recently released hydro-simulation TNG50-1 (Nelson et al. 2019a,b;Pillepich et al. 2019) (Section 2) balances resolution and sample size.The IllustrisTNG project is a suite of state-of-the-art cosmological hydro-dynamical simulations with well calibrated baryonic sub-grid physics implemented to reproduce the properties of local galaxy population, such as galaxy morphology and the stellar mass function (Springel et al. 2018).In this work, we aim to study the prevalence, origin and time evolution of the MW's Satellite Plane by using TNG50-1.
The structure of this article is organized as below: The observation data, numerical data and all kinds of measures we used in this work are enumerated in Section 2. Section 3 presents the criteria used to select the model halo.Detailed analyses of the selected model halos are described in Section 4. Section 5 gives a discussion.Finally, Section 6 summaries the main findings.

Observation data.
The 14 brightest MW satellites (11 classical satellites plus Canes Venatici I (CVn I), Crater II and Antlia II) are considered in this work.Although many more satellites are observed around the Milky Way, the census at lower stellar mass is not observationally complete and the simulation can't resolve these masses in any case(see Section 2.2).Many previous studies (i.e.Koposov et al. (2007); Tollerud et al. (2008); Walsh et al. (2009); Carlsten et al. (2020Carlsten et al. ( , 2021) ) )have assumed the census of satellites with M ⋆ ≥ 10 5 M ⊙ is complete throughout the virial volume, while some disagree (Yniguez et al. 2014;Samuel et al. 2020).Here we follow the accepted assumption that the satellite sample with M ⋆ ≥ 10 5 M ⊙ or M V ≲ −8 is complete.Thus the brightest 14 satellites serve as the point of convergence between simulation resolution and observational completeness.
For the 11 classical satellites and CVn I, sky coordinates, heliocentric distance and line-of-sight heliocentric velocities are taken from McConnachie (2012), while for Crater II and Antlia II, these data are taken from Torrealba et al. (2016Torrealba et al. ( , 2019)).Proper motion data for all the 14 satellites are taken from McConnachie & Venn (2020) which uses the Gaia Early Data Release-3 (EDR3).In order to estimate the observed uncertainties of some quantities (e.g.c/a, ∆(k).see the following section), the same methods as Samuel et al. (2021) are used here.For each satellite, the heliocentric distance, line-of-sight heliocentric velocities and proper motion are sampled 1000 times by assuming a Gaussian error distribution, then converted to a 6-D (positions and velocities) Cartesian Galactocentric coordinates by using Astropy (Astropy Collaboration et al. 2013, 2018).

Numerical Simulation.
Throughout this work, we compare observations with the simulation data from TNG50-1.This has a box size ∼ 51.7 Mpc, dark matter resolution ∼ 4.5 × 10 5 M ⊙ and stellar particle resolution ∼ 5 × 10 4 M ⊙ .It adopts the best fitting cosmological parameters in Planck Collaboration et al. (2016): dark energy density Ω Λ = 0.6911, matter density Ω m = 0.3089, baryon density Ω b = 0.0486, Hubble constant H 0 = 67.74kms−1 Mpc −1 , normalization σ 8 = 0.8159 and the spectral index n s = 0.9667.The compromise between sample size and simulation resolution of TNG50-1 make it ideal for analyzing the satellite galaxies in MW-mass halos.

Measures of disc thickness.
We use two methods to measure the thickness of the satellite plane: minor to major axis ratio (c/a) and rootmean-square (RMS) height h.
Ratio c/a is calculated as the ratio of square root of the minimal and maximum eigenvalues of the inertia matrix M: where r i is the position vector of i th satellite relative to the halo center and r 0 is the geometric center of the system of satellites: r 0 = 1 Nsat i r i .The normal of the plane nc is defined as the direction of the minor axis (i.e. the eigenvector corresponding to the smallest eigenvalue c).
As for the RMS height h, we use a similar definition as Samuel et al. (2021).We randomly generate 10 4 planes centered on the geometric center of satellites r 0 , then iteratively calculate the RMS height: where n⊥ is the unit normal vector of the random plane and ⃗ x i is the position vector of i th satellite with the origin is on the geometric center of those satellites.The minimum value of h is the RMS height of satellites and the corresponding plane is defined as the midplane.The smaller the values of h and c/a, the thinner the plane.

Measures of kinematic coherence.
We use the k − ∆ k relation (Pawlowski & Kroupa 2020) to characterize the kinematic coherence.For a given group of N sat satellites, we flip the direction of angular momentum of retrograde satellites at first (i.e.(−1) • n retrograde ), as we treat retrograde and prograde satellites equally here.Specifically, we first randomly generate 1000 direction vectors.For each random direction n random, j , (j = 1, 2, • • • , 1000) , we calculate the RMS angle distance D j between the angular momentum direction of those N sat satellites and the random direction n random, j based on the formula: where n i represents the angular momentum direction of i th satellite and we note that there is a absolute sign in the formula.Then the random direction with the minimal angle distance (i.e.min {D j }) is defined as the best direction n best and the angular momentum of a satellite would be flipped if the angle between the direction of its angular momentum and the best direction n best (i.e.arccos (n best • n i )) is larger than 90 • .After flipping the angular momentum of retrograde satellites, the k − ∆ k relation is obtained through the following steps: 1.For a positive integer k with k ≤ N sat , we can get Nsat k different combinations of k satellites among the N sat samples.
2. We calculate ∆(k) for each combination according to the following formula: Where n i is the angular momentum direction of i th satellite and ⟨n⟩ k represents the average direction of the k satellites.
3. Then ∆ k is defined as the minimal ∆(k Repeating those processes for a series of k values, we will get curve k − ∆ k . Based on this curve, only those model galaxies which meet ∆ k,model ≤ 84th percentiles of ∆ k,obs (i.e.≤ median + 1σ), ∀k ∈ [4, 5, • • • , 14] are considered to be kinematically comparable to MW.
In addition to the angular dispersion ∆ k , we also identify whether the satellites are orbiting in the plane they formed.This is important as it determines if this planar configuration is stable for a long time.We calculate angles θ k between the normal of the plane nc and the average angular momentum direction ⟨n⟩ k of the k satellites which have minima ∆(k) among all the Nsat k combinations (i.e.min{∆(k)}): (5)

SELECTING MODEL HALOS
The first step in attempting to find systems to compare with the observations is to identify MW "like" halos in the simulation at z = 0.These are (generously) selected with M 200 ∈ [0.3, 3] × 10 12 M ⊙ and M ⋆ ∈ [1, 10] × 10 10 M ⊙ , where M 200 is the total mass of the halo enclosed in a sphere with an average density 200 times the critical density of the universe and M ⋆ is the total stellar mass of the central galaxy within twice the stellar half mass radius.An additional criterion is applied to ensure that the MW halos that have been identified are located in environments similar to ours.An isolation requirement is imposed that ensured no halo more massive than 0.5 × 10 12 M ⊙ exists within 780 kpc, the distance between the MW and M31.Considering both the observational completeness limits of the census of the satellite galaxy population and the resolution of TNG50-1, we select the top 14 brightest satellites (as done by (Samuel et al. 2021)), which have M ⋆ ≥ 10 5 M ⊙ and are between 15 kpc and 300 kpc from the center of the host galaxy.Any halo which doesn't meet these criteria (i.e.doesn't have at least 14 satellites in this volume) is omitted from consideration.For the remaining halos, the 14 satellites with the most massive stellar mass are selected for further analysis.For reference: 548 halos in the simulation fall into the stated halo and stellar mass range; of these 404 are isolated and have no massive halo closer than the distance to Andromeda; of these 231 have at least 14 satellites that are above the stated stellar mass limit, and sit outside the host galaxy disc but within 300kpc.
Only one single halo at z = 0 survives all these constraints and matches the thinness and orbital coherence of the MW system.Namely it is as massive as the MW in both DM and stellar mass, is isolated according to the criteria stated earlier and has a similarly thin and co-rotating satellite structure.The ID of the halo is 395 (i.e.haloID=395, at z = 0).In the following we present detailed analysis of this halo395 and investigate how its satellite plane forms in the simulation.

Basic Characteristics
Halo395 has a dark matter mass M 200 = 4.9 × 10 11 M ⊙ , which is at the lower limit of observational constraints on the MW halo mass (Wang et al. 2020), though some recent results favor a lower MW halo mass (e.g.Bird et al. (2022)).Its stellar mass is M ⋆ = 2 × 10 10 M ⊙ , which is similar to the accepted values for the MW (e.g.M ⋆ = 6 × 10 10 M ⊙ (Licquia & Newman 2015)).The halo hosts a spiral disk galaxy as well (see panel(a) of Figure 5).The stellar mass (within two half-light radii) of the most luminous satellite is ∼ 3 × 10 7 M ⊙ , thus this model galaxy doesn't possess satellites as massive as the Large Magellanic Cloud (LMC), whose mass is estimated at M ⋆ ∼ 1.1 × 10 9 M ⊙ or the Small Magellanic Cloud (SMC), whose mass is estimated at

The Satellite Plane Problem
Figure 1 shows an edge-on view of the satellite plane (both halo395 and MW), with satellites colored red or blue depending on whether their tangential velocity is approaching or receding in this projection.From this plot the satellites of halo395 under consideration are dis-tributed in a thin plane, and most satellites are orbiting in this plane in a manner that is consistent with rotation.The ratio c/a and RMS height h of this halo's satellite population is 0.2 and 27 kpc respectively while, in the case of observed MW, the values are 0.25 and 27 kpc respectively .As for orbital coherence, the curve k − ∆ k of the model halo is below that of observed one when considering a 1σ uncertainty (see Figure 2 and Section 2). Figure 3 shows the k − θ k relation of the model galaxy and of observation with 1σ uncertainty.Obviously, the model galaxy behaves better than the MW at all values of k.Furthermore, the angle between the average orbit pole of the satellites and the norm of the satellite plane is also ∼ 21 • which is smaller than the observed value of ∼ 31 • (see Figure 3).Besides that, in Figure 4, the satellites radial distributions of halo395 and MW are compared.It is clear that the satellites radial distribution of halo395 is roughly consistent with that of the MW.
The satellite plane of halo395 is almost perpendicular (∼ 87 • ) to the central stellar disk (see the left panel of Figure 1).In the case of the observed MW, the angle between the satellite plane and the MW disk is ∼ 81 • (see the right panel of Figure 1).Thus not only does this halo have a thin rotating plane, but it is oriented in the same sense as the observations namely roughly perpendicular the MW's disk.In sum: halo395 which resembles the MW in terms of its halo mass, its stellar mass, its morphology, its isolation and its satellites radial distribution, also bears a thin satellite plane that's spinning in a coherent manner and is oriented exactly in the same way with respect to the Galaxy itself.

Environment of the Halo
More surprising results are seen in Figure 5   close to the data (770 ± 40 kpc (Holland 1998;Joshi et al. 2003;Ribas et al. 2005;McConnachie & Irwin 2006) and 3.8 ± 0.1 Mpc (Harris et al. 2010) respectively).Again, we note that our model 'M31' also has a slightly lower mass than the real one (∼ 5 × 10 11 M ⊙ vs ∼ [6−20]×10 11 M ⊙ ) but the mass ratio between the simulated and observed LGs is nearly identical.On large scales, this 'local group' is embedded in a sheet with a void above and below, as shown in panel (e).In short, the model galaxy and its environment indeed resembles  the cosmography of the MW and the local group (Tully et al. 2019).

Possible Origin of the Satellite Plane
Previous studies (Garrison-Kimmel et al. 2019b;Brooks & Zolotov 2014;Dutton et al. 2016;Buck et al. 2019;Sales et al. 2022) have shown that other smallscale challenges can be solved with well-designed subgrid models, while the Satellite Plane Problem has often been an obstacle and solutions are in debate.It was suggested (Libeskind et al. 2015) that the plane of satellites originated from the geometry of the cosmic web and the same result was also obtained here.The upper panel of Figure 6 shows that 11 of 14 satellites (marked with a dotted box) fell into the model galaxy along the plane of a large cosmic sheet.The dispersion of the satellite's angular momentum decreased as the sheet, they were formed in, become thinner due to the gravitational collapse of the sheet (and expansion of the voids above and beneath the sheet, see Figure 7 and Figure 6).After most of (10/11) the 11 satellites enter and evolve in halo395, their angular momentum remain unchanged, as the potential of the halo is almost isotropic (see the middle panel of Figure 7), which leads to two important results: (1) The dispersion of each satellite's orbital angular momentum keeps a small value throughout the history of this halo and (2) the ratio c/a of the satellites oscillates around a fairly low value.The reason for point (2) is that angular momentum conservation keeps the satellites moving in their respective orbits while the relative positions of the satellites change periodically thus leading to an oscillation of the ratio c/a.Because the angular dispersion is small (as point ( 1) shows), the oscillation of the ratio c/a is around a small value of ∼ 0.2 and has a fairly small amplitude (see Figure 7).In this case, a rotating thin plane is more likely to be 'observed' at present time (i.e. at z = 0), although the ratio c/a is oscillating and today's value is by no means its lowest (i.e. at a lookback time of ∼ 6 Gyrs, the value of c/a ≈ 0.1).The reader will note that the other three satellites just happen to appear at the right positions with the right velocities at z = 0 (see Figure 6 and the orange curve in the third row of Figure 7), culminating in a MW-like rotating thin plane (with 14 members) at that time.We conclude that, for halo395, the rotating thin plane is transient(as the value of c/a is oscillating and the appearance of the other three satellites is by chance) and that the accretion along the local super-galactic sheets can indeed lead to a rotating satellite plane.It is worth noting that Sawala et al. (2022) proves that the satellite plane of MW is also transient by using EDR3 proper motion data.

DISCUSSION
Unlike many previous studies(e.g.Pawlowski (2018)), we don't account for the effect of obscuration of the host stellar disk here.In these studies, they often use a simple correction when selecting satellites: exclude everything that lies within |b| ⩽ 12 • where b is the galactocentric latitude.We don't do the correction for two reasons: 1.As mensioned in Section 2.1, we assumed the satellites sample with M ⋆ ≥ 10 5 M ⊙ is complete, like many other works.correction will contaminate the analysis of the origin of satellite planes, however, which is one of the main goals of this work.
The halo mass M 200 = 4.9 × 10 11 M ⊙ of halo395 may be a concern, as it is at the lower end of observed data.However as some previous studies (Pawlowski et al. 2019b;Samuel et al. 2021) found, there are only weak or no correlations between the satellite plane metrics (h, c/a and orbital dispersion) and host mass.We conclude that this isn't a particularly critical problem.Furthermore if the halo mass is low, its only off by a factor of < 2.
As shown by Sawala et al. (2022), the metric c/a correlates with the satellites radial distribution, since the metric c/a is sensitive to the radius of satellites.The more centrally concentrated the satellites(higher Gini coefficient, see Sawala et al. (2022)), the lower the c/a value.An alternative radius-independent metric is the reduced inertial tensor method (Bailin & Steinmetz 2005)  .The evolution of c/a, the angular momentum of 11 satellites and their dispersion.We consider 11 satellites which are the red dots surrounded by a dotted box in Figure 6, and we obtain the evolution of their individual angular momentum, θ, and the dispersion of their angular momentum (σ θ ).The upper panel shows the evolution of σ θ , where prograde and retrograde satellites are treated equally.The middle panel shows the evolution of angular momentum θ for each individual satellite.Here θ is measured between the angular momentum of a satellite and an arbitrary fixed direction vector.The lower panel shows the evolution of c/a of the satellite plane composed of 11(blue line) and full 14(orange line) satellites.The orange line starts at Lookback Time =9.5 Gyr, because one of the 14 satellites didn't form until then, see the caption of Figur 6.Note that only here we consider all 14 satellites.The four vertical dotted lines (with black, green, red and purple color) indicate the epochs when the 1st, 5th, 10th and the 11th (the last) satellite are accreted by the model galaxy respectively.The accretion time is defined as the time when the satellite first crossed the halo virial radius.Actually, the 11th satellite never crossed the virial radius, but it was within 300 kpc at z = 0.It is also worth noting that, in the middle panel, the satellites represented by the three solid lines (with θ ∼ 150 • ) at the top are retrograde in the plane, and others are prograde.It is seen that for most of the lookback time, the angular momentum of satellites are almost kept fixed since their accretion.
in the same way as the inertia tensor, except it uses unit position vectors.In this work, we use the common metric c/a, as we want to find satellite systems that are the same as MW's in terms of full spatial anisotropy (by using c/a), not just in terms of angular anisotropy (by using (c/a) Red ).Interestingly, halo395, selected by using c/a, possess a similar satellites radial distribution as that of MW (see Figure 4).Interestingly, it is found that the (c/a) Red of MW and halo395, at z=0, are 0.342 and 0.319 respectively.This means that either c/a of halo395 indicates that the satellites are indeed in a thin plane.
It is noted that we use a slightly different metric to measure kinematic coherence (i.e.k − ∆ k see Section 2) when compared with some other works (e.g.Samuel et al. (2021)).In those works, they use the same for-mula as us, but they only consider a specific k.Since there is no better motivation for any particular value of k over any other, we consider all integers ranging from 4 to 14 (i.e. the full satellite sample size).This means the criteria we used is both more strict and more general, and the result we obtain is more robust.
We have shown that halo395 is embedded in a sheet plane (Figure 5 and Figure 6) and the formation of satellite plane of halo395 is related to its large scale environment (i.e.voids and sheet).Moreover, the norm of the satellite plane of halo395 points to the void rightly above the sheet plane (can be directly seen by eye from Figure 5 and Figure 6).Interestingly, Libeskind et al. (2015) showed that MW is embedded in Local Sheet and the norm of its satellite plane also points towards the local void above and below the Local Sheet, which suggest that the scenario of the formation of the satellite plane found here may also apply to the real MW.
It is also worth noting that as angular momentum is a vector, one may want to distinguish the prograde or retrograde orbits differently, as some previous studies (Pawlowski & Kroupa 2020) have done.Here we also check the k − ∆ k relation using Eq. ( 4) without flipping the angular momentum of retrograde satellites and the result is shown in Figure 8.It shows that the dispersion ∆ k of the model galaxy is larger than that of MW when k ≥ 9.This is because some model satellites are prograde and the rest are retrograde, but all of them are orbiting in a thin plane as shown in the Figure 1 and Figure 2. Although, halo395 is a bit different from the MW in this context, the mechanism we proposed (see Section 4.4) here may still apply to the satellite plane of MW.However, some recent studies (e.g.Li et al. (2021)) found that most fainter members of the satellite plane of MW are co-orbiting.In this case, it may be hard to explain it with our scenario alone, but Santos-Santos et al. (2023) pointed that the late capture of LMC-like object could increase the fraction of co-orbiting satellites (notably, there are no LMC analogs in halo395, see Section 4.1).So, the late capture of massive satellites like LMC may be also needed to explain the asymmetry in the direction of orbits.
Of course, we can not draw as a conclusion that the Satellite Plane Problem is solved, by only considering the plane around the MW, as some other satel-lite planes were found around other galaxies like M31 and CenA ; there are hints of satellite planes around more distant galaxies well beyond the Local Group e.g.M81 (Chiboucas et al. 2013), M83 (Müller et al. 2018), M101 (Müller et al. 2017).So, to truly understand the Satellite Plane Problem requires 1. large statistically significant observational samples which include more distant galaxies (e.g.Ibata et al. (2014); Mao et al. (2021)).
2. high resolution cosmological simulations with large volumes.

SUMMARY
In this work, we studied the Satellite Plane Problem of the MW, by using numerical data of the newly published simulation TNG50-1 which has high resolution and large volume at the same time.We searched the simulation to see if any halo had a a thin, rotating satellite plane like the MW, and if such halos do exist, how do the planes formed and how long can it last.Below are our major findings: 1.One halo (haloID= 395,z=0) is similar to the MW in the sense of halo Mass, stellar mass, morphology, environment and the satellites radial distribution and does not suffer from the Satellite Plane Problem (out of 231 candidates) in TNG50-1.Thus the MW-like satellite plane can be found in TNG50-1 simulations based on ΛCDM model.
2. The spatially thin and kinematically coherent satellite plane found at z=0 is transient (see Figure 7 and Section 4.4).This is in agreement with the findings of Sawala et al. (2022).
3. The major part (11/14) of the satellite plane of halo395 accreted along the sheets.Thus, for our special case (halo395), the accretion of satellites along the sheets leads the formation of a thin and rotating plane (see Figure 5, Figure 6 and Section 4.4).However, large samples of MW-like halos in different environments are needed to prove whether a special geometry (i.e.sheet and voids ) can boost the formation of thin and rotating plane, which is beyond the scope of our current work.
We thank the anonymous referee for very useful comments.We thank Marcel S. Pawlowski for his helpful suggestions and comments.We are also grateful to Yaoxin Chen for nice discussions.

Figure 1 .
Figure 1.Edge-on view of the satellite plane of halo395 and MW.This figure is for Satellite Plane problem, which shows the positions of the top 14 brightest model/MW satellites (blue/red stars), and central stellar disk (thick black solid line).The view is oriented to show the edge-on of both the central stellar disk and the satellite plane (long black solid line).Color blue/red indicate that satellites are moving out of/into the screen.The long black dashed lines indicate the root-mean-square height of the plane (Section 2).The left/right panel depict halo395/MW.

Figure 2 .
Figure 2. the k − ∆ k curve.Yellow data points with error bars show the observational results of the MW with 1σ uncertainty.The magenta triangles are for the model galaxy, halo395.

Figure 3
Figure 3. the k − θ k relation.Yellow data points with error bar show the result of the MW with 1σ uncertainty, and magenta triangles for the model galaxy.Here θ k is the angle between the average of the satellite angular momentum and the normal of the satellite plane.This quantity measures how well is the orbital plane aligned with the spatial plane.

Figure 4 .
Figure 4. the satellites radial distribution of MW and halo395.The black dashed line with gray shaded region shows the observational result of MW with 1σ uncertainty, and the orange solid line is that of halo395.Apparently, the satellites radial distribution of halo395 broadly match the observation.

Figure 5 .
Figure 5. the zoom-in view of the model galaxy and its surroundings at z=0.Here, the views of all panels are oriented to show the edge-on view of the satellite plane.From panel (a) to panel (e), the upper right panel shows the image of the central spiral galaxy of the model halo (halo395).Panel (b) shows the distribution of all luminous satellites ( M⋆ > 10 5 M⊙), and where the top 14 brightest satellites are marked by white solid circles.The central cyan solid circle indicates the central stellar disk, which is almost perpendicular to the satellite plane.Panel (c) shows a 'M31' analogue (M200 = 5.2 × 10 11 M⊙).The 'Centaurus A/M83 group' ( with halo mass M200 = 1.4 × 10 13 M⊙ ) is shown as big circle in panel (d) where the points represents all subhalos with M⋆ > 10 8.5 M⊙.Panel (e) shows the large-scale environment around the model galaxy.The dotted circles in each panel are scaled by the virial radius of the halos in that panel.The whiteness of points in each panel are scaled to the stellar mass of the satellites.The thickness of slice (e) and (d) are 3 cMpc and 8 cMpc respectively, where cMpc is the co-moving radius.

Figure 6 .
Figure6. the large-scale structure around the model galaxy and its evolution.The view of this figure is rotated so that the normal of the satellite plane is aligned with the Z-axis.Each column shows the projection to the X-Y, X-Z and Y-Z planes.Each row is for different epoch(redshift: 5.85, 3.28, 2.32, 0 or Cosmic Age: 0.96, 1.94, 2.84, 13.8  Gyr) of the evolution.The green cross and the red dots represent the central galaxy and satellites (top 14 brightest) respectively.In all panels, only subhalos with total mass ≥ 4.5 × 10 8 M⊙ are shown.The blue, orange and green circle indicate the virial radii of 'Centaurus A/M83 group', 'M31' and the model galaxy.Red dots surrounded by a dotted rectangle in the upper left panel fell on the model galaxy along a sheet plane and are analyzed in Figure7.It is worth noting that only 13 red dots were shown on the first three rows, as one of the 14 satellites formed very recently until z = 1.41 or Cosmic Age = 4.5 Gyr.The thickness of all slices (panels) are 3 cMpc.

Figure 8
Figure 8. the k − ∆ k curve, computed using Eq.(4) without flipping the angular momentum of retrograde satellites.All data points are similar to Figure 2.