A Study of the Properties and Dynamics of the Disk of Satellites in a Milky-Way-like Galaxy System

The dynamics of the satellite systems of Milky-Way-like galaxies offer a useful means by which to study the galaxy formation process in the cosmological context. It has been suggested that the currently observed anisotropic distribution of the satellites in such galaxy systems is inconsistent with the concordance $\Lambda CDM$ cosmology model on the galactic scale if the observed satellites are random samples of the dark matter (DM) sub-halos that are nearly isotropically distributed around the central galaxy. In this study, we present original high-resolution zoom-in studies of central galaxies and satellite systems based upon initial conditions for the DM distribution from the Aquarius simulations but with substantial high-resolution baryon physics added. We find that the galaxy most like the Milky Way in this study does indeed contain a disk of satellites (DOS). Although one galaxy DOS system does not answer the question of how common such disks are, it does allow the opportunity to explore the properties and dynamics of the DOS system. Our investigation centers on the spatial arrangement (distances, angles, etc.) of satellites in this Milky-Way-like galaxy system with a specific emphasis on identifying and analyzing the disk-like structure along with its dynamical and morphological properties. Among the conclusions from this study, we find that the satellites and DM sub-halos in the galaxy simulations are anisotropically distributed. The dynamical properties of the satellites, however, indicate that the direction of the angular momentum vector of the whole satellite system is different from the normal direction of the fitted DOS and from the normal direction of the velocity dispersion of the system. Hence, the fitted DOS appears to be comprised of infalling sub-halos and is not a rotationally supported system.


Introduction
The alignment of the 11 originally discovered most luminous, or "classical", satellites of the Milky Way relative to the Milky Way stellar disk [1][2][3] has been known yet unexplained for decades.For example, Kroupa et al. [3] fit the spatial distribution of the then available satellite samples into a thin disk that is highly inclined relative to the Milky Way stellar disk.It was then claimed [3] that this observation is inconsistent with the hierarchical structure formation in ΛCDM cosmology on the galactic scale at a high statistical significance level (99.5 percent).This is true as long as the observed satellites only trace a distribution of the dark matter (DM) sub-halos that reflects the shape of their host DM halo (see Kroupa et al. [4] and Kroupa [5] for extended reviews).However, this might not be true if, for example, the satellite distribution reflects the asymmetry of filaments in the cosmic web.
It has also been suggested that at least some of the satellites have a coherent motion in phase space [3,12,[17][18][19][20][21].However, it is not altogether clear whether this indicates that they may have formed from the same origin or that they arrived individually.
Many efforts have been made to give a natural explanation of the DOS phenomenon.Early studies [32,33] suggested that satellites in disk galaxy systems have a tendency to avoid the equatorial region relative to the disk with a spatial distribution along a certain special direction.However, data from the 2dF Galaxy Redshift Survey and the Sloan Digital Sky Survey (SDSS) have suggested that this may only be true for host galaxies with low present-day star formation rates (spectral type η < −1.4) and for satellites with low radial velocity (<160 km s −1 ) relative to the primary [20,[34][35][36][37][38][39][40].
Developments in numerical simulation techniques for galaxy formation and evolution, especially for the baryonic physical processes, have made it possible to study the DOS phenomenon numerically in detail.Early simulations were based on using pure dark matter [41] or semi-analytic galaxy formation [42,43] simulations in a ΛCDM cosmology.Those studies found that certain satellite systems formed in the simulations are indeed distributed along a special direction that resembles the spatial distribution of the satellites of the Milky Way.They also suggested that the DOS phenomenon may naturally exist in a universe that is described by a ΛCDM cosmology model.
Subsequently, however, it was argued that, in contrast to the assumption made by [3], the distribution of the satellites does not trace the overall distribution of the DM sub-halos.Rather, it is close to that of the sub-halos with the most massive progenitors in the accretion process [42].The direction of the alignment of the satellites has been found to be close to the major axis of the host DM halo in the simulations of [20,[42][43][44][45][46][47][48].With a large sample size of 436 and 431 parent haloes, respectively, Libeskind et al. [46] and Deason et al. [47] find that statistically the particular type of the spatial alignment of the satellites of the Milky Way is not uncommon (20%) among similar systems.However, Wang et al. [48] found that the possibility is much lower (∼1%).
The possible coherent motion of the satellites in phase space has been investigated by [21,[45][46][47]49,50] with simulations.It was found in their simulations that in some galaxy systems, the direction of the net angular momentum of the satellite system can be close to that of the main galaxy.However, it is rare (<10%) [46,47,51] for satellite systems in the simulations to have an orbiting pole alignment like the one found for the satellites of the Milky Way [18].
An alternative explanation to the coherent motion of the satellites in the phase space is that the observed satellites are the tidal dwarf galaxies (TDGs) that formed in a close encounter between two galaxies in the Local Group [1,3,18,[52][53][54].The formation of TDGs in galaxy interactions has been shown to be quite common (e.g., [55,56]).Using a large sample of general galaxy interaction simulations, Refs.[29,57] have shown that the resulting TDGs can easily have a similar spatial distribution and orbiting pole to those of the satellites of the Milky Way.Furthermore, using constrained simulations [58][59][60] of the interaction between the Milky Way and the tidal tail of the Andromeda galaxy, Fouquet et al. [61] found that the resulting tidal tail matches the spatial distribution and phase space properties of the Milky Way satellites.However, there are still several problems with this theory: First, TDGs are generally believed to contain only baryonic matter (gas) [56,62].So, if a majority of the satellites of the Milky Way are indeed TDGs, they cannot have a very high mass-to-light ratio contrary to current observations that show a ratio of up to a few thousand (e.g., [63,64]) except in some particular cases [53,65,66].Second, one needs realistic simulations to study the detailed physical properties of these gas-rich TDGs in the proposed scenarios beyond the ones in [57,61].And finally, the properties of the proposed galaxy candidates in the galactic interactions simulated in [57,61] are still in need of close examination as to whether they match the observations.
On the basis of previous and new (Gaia DR2) proper motion data of the 11 classical satellites of the Milky Way, Pawlowski and Kroupa [50] showed that there is indeed an alignment of the orbital poles of 7 out of 11 classical satellites.Using the TNG100 simulation data, they argued that such an alignment is extremely rare (0.1%) in ΛCDM cosmological simulations.
However, Samuel et al. [67] used FIRE-2 simulations [68] to show that spatially thin and/or kinematically coherent satellite planes in Milky-Way-like galaxies are not very rare at the level of a few percent.They also found that an LMC-like satellite will greatly improve the odds of finding such systems.The conundrum of satellites that are aligned with the DOS but do not appear to be satellites of the LMC has been explored in the context of a single LMC pericentric passage (e.g., [69]), and a second LMC pericentric passage [70].The authors find that the latter scenario leads to a larger debris field and could result in an additional 4-6 dwarf galaxies lost from the Magellanic system.
A recent study [21] of the new proper motion data for M31's satellite, NGC 147 and NGC 185 has led to the suggestion of the Great Plane of Andromeda (GPoA).It was shown that these two satellites indeed have an orbital plane aligned with the GPoA.Additionally, they showed that while it is not very rare (more than 4%) to find two co-orbiting satellites around a Milky-Way-like galaxy in cosmological simulations (e.g., TNG50, TNG100, ELVIS, and PhatELVIS), it is extremely rare (1 in 1000) for these two co-orbiting satellites to also reside within a thin satellite plane around the main galaxy.
In another recent work, Sawala et al.
[31] made use of a collection of zoom-in constrained simulations.That paper claimed that a thin plane of satellites like that of the Milky Way is relatively common (5.5%) but short-lived if numerically disrupted satellites are included.
Also, using semi-analytical galaxy formation models, Gu et al. [28] found that 4.7% of galaxy clusters and 13.1% of MW-mass systems have planes of satellites that are thinner than that of the current Milky Way in their simulations.Although the satellite accretion is highly anisotropic at higher redshift, whether it can account for the narrow aspect ratio of the currently observed satellite planes remains inconclusive.
More recently, using Magneticum Pathfinder simulations and a new plane of satellites finding algorithm, Förster et al. [71] found that thin planes of satellites exist in almost all simulated galaxy systems.
By studying the accretion history of the Local Group in the HESTIA simulation suite, Refs.[72,73] also found that satellites in the Local Group accreted along the axis of the slowest collapse of the local filament.
Most recently, Pham et al.
[74] used zoom-in high-resolution dark-matter-only simulations to show that the satellite system of the Milky Way is unusual but still consistent with the ΛCDM model at the 2-3σ level, although it also depends upon the radial distribution of the satellites.
Finally, Xu et al. [75] did find in TNG50 simulations the presence of one satellite system around a Milky-Way-like galaxy that closely resembles the Milky Way both spatially and dynamically.Although the majority of satellites seem to accrete along the local filaments, the satellite plane is transient in nature.
Here, we similarly examine a DOS system within the context of higher resolution simulations of a Milky-Way-like galaxy system.The original simulations presented in this paper follow the methodology established by the Aquila Comparison Project [76,77].These adopted the initial conditions from the pure dark matter Aquarius simulations [78].Our investigation focuses on the spatial arrangement of satellites in the Milky-Way-like galaxy system that most resembles the Milky Way, that is the C halo.Interestingly, in our simulations, this galaxy contains a DOS system.This allows the opportunity to conduct a dynamical and morphological analysis in this higher resolution context.Our specific goal is to shed light on the nature of the anisotropy, alignment and the potential implications it holds for our understanding of the formation and evolution of the DOS system around a Milky-Way-like galaxy within the framework of ΛCDM cosmology.
This paper is organized as follows: we describe the details of our simulations in Section 2, the results are presented in Section 3, followed by a summary and conclusions in Section 4.

The Numerical Simulations
In order to study the formation and evolution of Milky-Way-like galaxy disk systems, we utilized original Smoothed-Particle hydrodynamics (SPH) simulations that were performed using the Research Computing and Cyber-Infrastrucure unit at Penn State University.These simulations are independent of but guided by results from the Aquila Comparison Project [76,77].As in the Aquila project, the original simulations reported here utilized initial conditions from the pure dark matter Aquarius study [78].A number of Milky-Way like DM halos with similar mass and evolution history to those of the Milky Way were identified in the Aquarius simulations.In the work reported here we used the initial condition of the "C" halo.This is because the final main galaxy formed in the Aquila Comparison Project with this initial condition has the best defined stellar disk [76,77].The simulations reported here utilized sufficient resolution of baryon particles to enable identification of the surrounding dwarf galaxies.Interestingly, in our analysis, we find that a DOS system exists.Although the existence of a galaxy with a DOS does not answer the question as to how frequent such systems are, the existence of this system allows an opportunity to study the properties and evolution of the DOS system.That is the purpose of the present paper.
The simulations focused on a zoom-in high-resolution region surrounding the selected halos.Unlike the Aquarius simulations, which only used dark matter particles, both dark matter and gas particles were used in the SPH simulations.We adopted cosmological parameters from the Aquarius simulations, with Ω dm = 0.21 and Ω b = 0.04 in the highresolution region.Assuming a ΛCDM cosmology, we set the cosmological parameters to Ω m = 0.25, Ω Λ = 0.75, H 0 = 73 km s −1 Mpc −1 and σ 8 = 0.9.Periodic boundary conditions were applied with a box size of 100 Mpc h −1 (137 Mpc for h = 0.73) in comoving coordinates.The simulations began at z = 127 with the chosen initial condition s and were evolved to the current epoch at z = 0.
The simulation code used in this study is the Tree-PM SPH code GADGET-3 [78,79], which is an improved version of the widely used simulation code GADGET-2 [79].The GADGET-3 code incorporates most of the important baryonic physical processes such as metallicity-dependent cooling, star formation, supernova feedback, black hole physics and a multi-phase gas model.Adopting the notation from the Aquila project, the "C−4" and "C−5" refer to the baryonic mass resolution level (cf.[76,80,81]).As noted in Table 1, in the zoom-in region of the "C−4" simulation (higher resolution one), a mass resolution of 2.57 × 10 4 M ⊙ was achieved for gas and star particles, and 2.70 × 10 5 M ⊙ for dark matter particles.For the "C−5" simulation (lower resolution one), a mass resolution of 2.06 × 10 5 M ⊙ was achieved for gas and star particles, and 2.16 × 10 6 M ⊙ for dark matter particles.Although we present results for the "C−5" simulations in Figures 1 and 2 as well as Table 1 below, this is only to demonstrate that the general characteristics do not change when going to higher resolution.However, since our goal is to consider satellites as small as 10 5 M ⊙ , the "C−5" resolution would contain fewer than one gas particle, making it inadequate for this study.Hence, we focused mainly on using the "C−4" simulation for the investigations in the present work.This resolution sufficiently supports a detailed study of low-luminosity satellites.By sufficient resolution, we mean that with a star particle mass of ∼3 × 10 4 M ⊙ , halos with 10 6 M ⊙ contain at least 30 star particles.To identify structures at various scales (main halo, sub-halos, central galaxy and satellites) in the simulations, a group finding algorithm was initially applied to the particle data.For this step, a friends-of-friends (FOF) algorithm was used to link the dark matter particles based on their mean particle spacing.Subsequently, gas and star particles were linked to the dark matter particle groups using the same algorithm.The most massive group identified was considered the main galaxy.For detecting gravitationally bound substructures within the galaxy system, we employed the SUBFIND algorithm [82].This method first identifies the local over-densities via the SPH kernel interpolation, followed by the iterative removal of gravitationally unbound particles.Any remaining groups with over 20 particles were classified as physically bound substructures.The characteristic parameters of the simulations and the main galaxies at redshift z = 0 in this study are listed in Table 1. a Mass of the gas, dark matter and star particles at high-resolution regions in the simulations.b Gravitational softening length.c Virial radius of the main galaxy defined by the radius of a sphere within which the mean matter density is 200 times the critical density.d Mass enclosed within the sphere of the virial radius.e Mass of the gas, dark matter and stellar components in the main galaxy excluding the satellites.f Number of the luminous satellites around the main galaxy within its virial radius or 1 Mpc radius.
Figure 1 shows the main galaxy formed in the simulations.Here, we render the edge-on view of the two-dimensional density projection map of the stellar, dark matter and gas components in the main galaxy system within its virial radius in the "C−4" simulation at redshift z = 0.In this figure and throughout the paper, we define the galactic coordinate system of the main galaxy in the simulations by choosing the direction of the angular momentum of the stellar component within three times the stellar half-mass radius as the normal direction of the stellar disk.Thus, for the edge-on view of the galaxy, the projection is perpendicular to the direction of the stellar component's angular momentum.It is worth noting that the angular momentum vector is not necessarily perpendicular to the stellar disk, leading to the slightly tilted appearance of the edge-on views in Figure 1.In order to show the structures clearly, we only include the particles in the physically bound substructures identified by SUBFIND.
As one can see from this figure, by comparing the distributions of the stellar and DM components in the system, only some of the DM sub-halos have baryonic luminous counterparts, and thus could be observed as satellites of the main galaxy.Hence, pure dark matter simulations are not sufficient to fully study the DOS phenomenon.The implementation of baryonic physics is essential for such an investigation.This aspect is related to the "missing satellite problem" [83][84][85][86] (and references therein).However, here we focus on the spatial distribution, not the population of the satellites.Figure 1 also shows that the gas component only exists in the main galaxy and that very few satellites appear at redshift z = 0. Most satellites are gas free.
Figure 2 shows the face-on and edge-on views of the central part of the main galaxy's stellar component within one-tenth of its virial radius in the "C−5" and "C−4" simulations at redshift z = 0.For a clear illustration of the structure of the main galaxy, we only include the particles in the main galaxy but not the ones in the satellites.It can be seen from this figure that substantial stellar disk structures are formed in the main galaxies in both of the "C−4" and "C−5" simulations.The galaxy in the "C−5" simulation has a more distinct disk, similar to the one observed in the previous "C−5" simulations [76,77], while the galaxy in the "C−4" simulation exhibits a less distinct disk due to an ongoing minor merger event at z = 0 in this simulation.This is evidenced, for example, by the small clump at (x, y) = (10, −5) in the "C−4" face-on view in Figure 2. We present both the "C−4" and "C−5" simulations here to illustrate that the fundamental features remain consistent while going to higher resolution.Therefore, for the remainder of our analyses, we focus on the higher resolution "C−4" simulations, where the satellite features are more distinctively defined.

Results
In this section, we present our results in three parts.First, we examine the spatial distribution of the satellites surrounding the main galaxy.Next, we investigate the dynamical properties of the main galaxy system, analyzing the relative motion of satellites to the main galaxy and the angular momentum across different components in the whole system.For these analyses, our focus remains on the current system at z = 0. subsequently, we explore the time evolution of the entire galaxy system to find possible explanations for the findings from the first two parts.

Spatial and Mass Distributions
We first investigate the effects of the satellite luminosity (proxied by mass in the simulation) on the observed anisotropic distribution of the satellites because the early discovered classical satellites of the Milky Way that were included in the previous investigations [1][2][3]8,9,17] (and references therein) are also the more luminous ones.We investigate various satellite samples defined by different lower mass thresholds (10 8 M ⊙ , 10 7 M ⊙ , 10 6 M ⊙ and 10 5 M ⊙ ) to determine whether the disk of satellite phenomenon is consistently observed across all samples, indicating its presence as a generic property for the entire satellite population and not solely a result of selection bias from observations.For reference, Figure 3 shows a histogram of the satellite mass function.There is a declining number of satellites with increasing mass, aligning with the expectations of ΛCDM models [87].Of course, one should consider the question of the detectability of satellites over this mass range within the radius of the simulations in this study.Considering the lowest mass (and therefore lowest luminosity) bin, if we adopt one solar luminosity for one solar mass, then a 10 5 M ⊙ satellite would have an absolute visual magnitude of −7.7.Even for a limiting G band magnitude of 21 as in the GAIA survey [88] (and ignoring the small G-band bolometric correction), these lightest satellites should be detectable out to more than 5 Mpc, which is well beyond the 1 Mpc size of the largest satellite distribution considered below.
In Figure 4, we plot the edge-on views of the main galaxy's stellar disk and the dashed line indicates the orientation of the DOS formed by the satellite samples with minimum stellar masses of 10 8 M ⊙ , 10 7 M ⊙ , 10 6 M ⊙ and 10 5 M ⊙ within the main galaxy's virial radius in the "C−4" simulation at redshift z = 0.
The disk plane of the satellites was fit with the Hesse normal form, where ⃗ n DOS is the unit normal vector of the plane, ⃗ x is the position vector of a point and d is the distance from the plane to the coordinate origin.The fitting uses a least square method on the positions of the satellites relative to the galaxy center as described in [3].With this method, the fitted planes in this study are not forced to contain the galaxy center as suggested by [12].The galaxy system in Figure 4 is oriented in such a way that the edge-on view is perpendicular to the angular momentum vector as in Figure 1, but the projection direction is now rotated so that the fitted DOS plane is always seen edge on.To offer a parallel perspective in the dark matter environment, in Figure 5, we also plot the edge-on view of the main galaxy's DM halo (defined as perpendicular to the DM halo angular momentum vector) and the DM DOS formed by all DM sub-halos within the virial radius.As in Figure 4, the dashed line marks the orientation of the DM DOS with the projection direction rotated so that the fitted DM DOS plane is seen edge-on.As noted below, it is interesting that the orientation of the DM DOS plane is similar to stellar DOS planes shown in Figure 4.
Figure 6 depicts the angle formed between the normal direction of the DOS and the angular momentum direction of both the main galaxy's stellar disk (red line) and the host dark matter halo (black dotted line), respectively.This angle is plotted against the minimum stellar mass criterion of the satellite sample.These measurements provide a quantitative assessment of the relative orientation shown in Figure 4. To determine the orientation of the stellar disk, we first chose the disk stellar particles with similar but simpler criteria than those used in [76].A stellar particle is categorized as a disk particle if: (1) it is located inside the main galaxy's disk radius (one-tenth of the main galaxy's virial radius) and outside the galactic center 1 (one-tenth of the main galaxy's disk radius); (2) the angle θ between the angular momentum of the stellar particle and that of the stellar disk satisfies, cos θ ≥ 0.7; and (3) the difference between the stellar particle's circular speed and the circular speed of the disk at the same location is no larger than 30% of that of the disk.While a more complicated method [76] is available, our current choice is sufficient.Our primary goal is to determine the disk's orientation within a few degrees, which is not much affected even if a few particles are misidentified.We then diagonalized the second-moment tensor of the disk stellar particles, to find the orientation of the disk using a similar algorithm to the one that we used to find the orientation of the DOS.
Figure 7 provides a schematic representation of the various angles depicted in Figures 6, 10, 12, and 14.Here, we take ⃗ n main ≈ ⃗ L main ( ⃗ L main illustrates the angular momen- tum vector of the main disk, while ⃗ n main is the direction perpendicular to the stellar disk of the main host).The quantity θ DOS approximately represents the angle between ⃗ n main and the normal direction to the DOS, ⃗ n DOS .The quantity ⃗ σ DOS = ∑ i ⃗ σ i is the sum of the dispersion vectors (minimum bulk motion) ⃗ σ i of the satellite galaxies in the sample.This direction is determined by summing the velocities of the particles in the satellite.The direction maximum sum indicates a bulk flow.Conversely, the direction of the minimum sum is only showing the random dispersion and is perpendicular to the direction of any bulk flow.Hence, one can use the sum of the dispersion to identify a bulk flow of the DOS.As one can see from Figure 4 and quantitatively from Figure 6, the orientations of the DOS formed by the satellite samples with different luminosity thresholds are close to each other.In the luminosity range studied, the DOSs always have inclination angles 2 above about 50 • relative to the stellar disk of the main galaxy.For some particular luminosities (e.g., M * > 10 7 M ⊙ ), it can be nearly perpendicular to ⃗ n main .The angles θ DOS between the normal direction of the DOS and the normal to the host dark matter halo are much smaller in the same luminosity range than the same angles for the stellar component.This suggests that the DOS exhibits a stronger alignment with the main DM halo compared to the stellar disk of the main galaxy.This alignment seems to imply a tendency for the DOS to be preferentially accreted along the major axis of the DM halo.We speculate that this may point toward the accretion of satellites along filaments.However, this would require a study out to much larger distances in the cosmic web.
For the DM sub-halos, we use the same fitting algorithm to find a DM DOS around the main galaxy.In this case all DM sub-halos with at least 20 particles (i.e., M DM ∼ 10 6 M ⊙ ) are counted.As shown in Figure 5, this DM DOS has a similar orientation to those of the luminous satellites in Figure 4.In particular, for the most direct comparison (the plot with M * > 10 6 M ⊙ ), the orientations of the DOS indicated by the dashed lines in Figures 4 and 5 are almost identical.Indeed, this suggests that the orientation of the DOS formed by the luminous satellites may be close to the preferred direction in the spatial distribution of the dark matter around the main galaxy.These two interesting phenomena prompted us to further study the distribution of the satellites and dark matter sub-halos on a larger scale simply because the DM extends to larger scales.
We next investigated the spatial distribution of the satellites 3 within a 1 Mpc radius around the main galaxy.We chose this scale because it is comparable to the scale of the Local Group.It is important to acknowledge that, as we extend our analysis to a larger radius, the issue of how distance affects the detectability of satellites becomes more pronounced for observational data.As noted above, the satellites discussed here should all be detectable by current surveys out to much more than one Mpc.Figures 8-10 are the identical plots on a 1 Mpc scale to that of the virial-radius scale of Figures 4-6, respectively.As on the smaller scale of the main galaxy's virial radius, it is clearly shown in Figures 8 and 10 that the orientations of the DOS formed by the satellites with different luminosity thresholds are close to each other on a 1 Mpc scale.The orientations of the DOS formed by the luminous satellites and the DM sub-halos are also very close to each other on a 1 Mpc scale as indicated in Figures 8 and 9.
This is further illustrated in Figure 10 by the trend of the angle θ DOS between the normal direction of the DOS and the angular momentum vector of the host dark-matter halo.The angle decreases monotonically from about 50 • for M ⋆ = 10 8 M ⊙ to about 10 • for M ⋆ = 10 4 M ⊙ .This result suggests that if we fit a complete sample of the satellites with the lowest possible luminosity using the same fitting algorithm, the fitted DOS plane should be very close to the long axis of the host dark matter halo.To compare the spatial distribution of the satellites in our simulations to observations, we plot the Hammer-Aitoff projection of the mock full sky map in Figure 11.This figure shows the main galaxy's stellar disk and the satellites within the virial radius (top plot) and those within 1 Mpc radius (bottom plot) in galactic coordinates using our simulation data.For this sky projection map, the galactic coordinates are defined by positioning the Sun on a line that passes through the galactic center.This line runs parallel to where the DOS plane intersects the main galaxy's stellar disk.We assume the Sun has a distance of 8.0 kpc from the galactic center.One easily recognizable feature of Figure 11 is the almost identical locations of the common satellite points between the virial radius map and the 1 Mpc radius map.This suggests that the DOS formed by the satellites on these two scales are very close to each other.Using the fact that the satellites with low galactic latitudes are often blocked by the Milky Way's own stellar disk and that the observed satellites of the Milky Way [4] in all have galactic latitudes higher than 20 • .We define the "Zone of Avoidance" of the main galaxy in this study as a region with galactic latitudes lower than 20 • as marked by dashed lines on Figure 11.
For a quantitative comparison, we list some of the key characteristic parameters of the DOS fitted in our study and those from some observational investigations in Table 2. Their definitions are explained in the footnote.For the galactic latitude b of the normal direction of the fitted DOS plane, we notice that the "C-4 Vir ZOA" data set (satellites within the virial radius but excluding those within the zone of avoidance) has a b value that is comparable to the values from the fitted DOS of the Milky Way and the "C-4 Vir" data set has a b value that is close to the one from the Andromeda galaxy (for further details on the specific characteristics of these data sets, please see the footnote below Table 2).a Names of the data sets that are defined as follows: "C-4 Vir" and "C-4 1Mp": satellite data of the main galaxy with minimum stellar mass of 10 5 M ⊙ within the virial radius and 1 Mpc radius in the "C−4" simulation presented in this study; "C-4 Vir ZOA" and "C-4 1Mp ZOA": same as "C-4 Vir" and "C-4 1Mp" except that the satellites in the "Zone of Avoidance" of the main galaxy are excluded in the fitting; "KTB(05) MW Vir" and "KTB(05) MW 1Mp": satellite data of the Milky Way from Table 1 of Kroupa et al. [3] up to a distance of 254 kpc and 956 kpc; "MKJ(07) MW": satellite data of the Milky Way from Table 2 of Metz et al. [12] with most satellites and the algebraic least-squares(ALS) fitting method; "MKJ(07) And.MI" and "MKJ(07) And.KG": satellite data of the Andromeda galaxy from Table 2 of Metz et al. [12] that are originally from Koch and Grebel [10] and McConnachie and Irwin [11] with the ALS fitting method; "MKJ(09) MW": satellite data of the Milky Way from Table 1 of Metz et al. [6] plus those in Metz et al. [12]; "MKJ(09) And.": satellite data of the Andromeda galaxy from Table 2 of Metz et al. [6] plus those in Metz et al. [12]; "Kroupa(10)": satellite data of the Milky Way from Table 2 of Kroupa et al. [4].b The furthest satellite's distance to the center of the main galaxy in the fitting.c Number of the satellites in the fitting.d Galactic longitude of the normal direction of the fitted DOS plane as defined in the data set.e Galactic latitude of the normal direction of the fitted DOS plane as defined in the data set.f Root mean square height of the fitted DOS plane.g Disk aspect ratio that describes the thinness of the DOS.h Distance from the fitted plane to the center of the main galaxy.
The reason that we use different data sets for the Milky Way and Andromeda galaxies is that the Andromeda galaxy does not have an observational "Zone of Avoidance" as the Milky Way does.The comparability of values obtained from the "C-4 Vir ZOA" and "C-4 Vir" data sets with those from the Milky Way and Andromeda, respectively, is intriguing because we only have one realization of Milky-Way-like galaxy systems from the "C" halo initial condition and yet we get values that are similar to those from the observations.At the very least, it shows that the high inclination angle of the fitted DOS plane in such a system is not unlikely because we only choose the "C" halo to have a similar mass and evolution history to the Milky Way with other factors chosen randomly.
The distance from the fitted plane to the center of the main galaxy is much less than the root mean square height of the fitted DOS plane in all of our data sets.Thus, the fitted DOS planes in our study meet the consistency check [12] that constituent satellites orbit within the host potential such that a disk made up of a virialized satellite population passes near the origin (or host galaxy).Disks with distances greater than one disk height from the host are considered unphysical.We measure "∆/R cut ", the flattening parameter defined in [3], where ∆ is the rms height of the disk and R cut is the farthest distance to the Galactic Centre in the satellite sample.This disk aspect ratio describes the thinness of the DOSs in our study.It is generally larger (∼10-30%) than the ones from observation (∼10-20%).This means that the fitted DOSs are thicker using the data from our simulations.However, when considering a scale of 1 Mpc, the fitted DOS appears thinner.This is because, at this larger scale, the satellite distributions align more closely with the shape of the local dark matter distribution.

Dynamical Properties
From observational data on the proper motion of the satellites, it has been argued that some of the Milky Way's satellites may have coherent motion in phase space.This coherence implies that they may share a common origin [3,12,17,18].While some numerical studies have shown [45][46][47]49] that the direction of the satellites' angular momentum can align closely with that of the main galaxy, it remains challenging [46,47] to identify a satellite system in the simulations that matches the orbiting pole alignment observed in [18].
The bottom black solid line on Figure 12 represents the angular separation between the normal direction of the DOS and the angular momentum vector direction of the DOS, plotted as a function of the minimum stellar mass of the satellites within the virial radius in the C−4 simulation at redshift z = 0.If the satellites within the fitted DOS are rotationally supported, these two vectors should closely, if not perfectly, align.However, the angular separation between these two vectors ranges between 20 • and 70 • showing no discernible trend with satellite mass, and the angle never approaches zero.This suggests that, at least within the context of the DOS fitted with satellites in this study's simulation, the satellites do not exhibit coherent rotation within the DOS plane.
The middle line in Figure 12 further supports this observation.This shows the angle between the direction of the angular momentum with respect to the center of the main galaxy and the normal direction of the velocity dispersion ⃗ σ DOS of the DOS.This is plotted as a function of minimum stellar mass in the satellites while keeping other parameters the same.Given that the normal direction of the velocity dispersion represents the direction of the least bulk velocity flow of the satellites, the very large angle of about 90 • in Figure 12 indicates a substantial bulk velocity flow along a direction parallel to the angular momentum vector L Main of the main galaxy.This suggests that the bulk motion of the fitted DOS in this simulation is characteristic of an infalling system, rather than an orbiting one.For an orbiting motion, the bulk motion would point parallel to the plane of the DOS, not perpendicular to it.Thus, there is no indication that the fitted DOS plane has a common orbital plane with the satellites within it.
To examine the relation between the directions of the angular momenta of the DOS and the main galaxy, in the top solid line on Figure 12, we plot the angle (solid line) between the direction of the angular momentum vectors of the DOS and that of the stellar disk of the main galaxy L main as a function of minimum stellar mass in the satellites.The dashed line shows the angle with respect to the host dark matter halo.There are several distinct features in this figure: First, the angle between the directions of the angular momentum of the DOS and the plane of the stellar disk of the main galaxy or the host dark matter halo is independent of the stellar mass in the satellites.Second, the main galaxy and the host dark matter halo have almost the same angular momentum direction.Most importantly, one can see in Figure 12 that the angle between the directions of the angular momentum of the DOS, the stellar disk of the main galaxy and the host dark matter halo is always about 100 • .This indicates that the direction of the angular momentum of the DOS with respect to the center of the main galaxy is almost perpendicular to that of the main galaxy and the DM sub-halos.This shows a very distinct difference in the dynamics between the satellites and the main galaxy.The dynamics of the satellites reflect the dynamics of infall toward the main disk.This differs from a rotationally supported system orbiting around the plane of the DOS.

Time Evolution
Until now, we have been concerned with the properties of the satellite systems during the current epoch at redshift z = 0.This is because most of the observational data of the satellites are from the current satellite systems of the Milky Way and Andromeda galaxy.More recent observations of extragalactic disk galaxy systems in data from the 2dF Galaxy Redshift Survey and SDSS have extended this analysis to higher redshift s to probe deviations of the satellite spatial distribution and orientation from those of the Milky Way and Andromeda systems [89,90].Therefore, it is also essential to study the time evolution of the satellite systems around the main galaxy if one wishes to understand the large-scale dynamics of such systems.
In Figure 13, we render the edge-on views of the main galaxy's stellar disk and the DOS formed by the luminous satellites within the virial radius (left column) and a 1 Mpc radius (middle column) and the DOS formed by all dark matter sub-halos within a 1 Mpc radius (right column) from redshift z = 5 to z = 0. From this figure one can clearly see that the DOS formed by the satellites on different scales and the DOS formed by the DM sub-halos have similar orientations relative to the stellar disk of the main galaxy at each redshift.This is especially true at lower redshifts (z < 3) as is manifested in Figure 14 in which we plot the angle between the normal directions of the DOS formed by the luminous satellites and that by all dark matter sub-halos within a 1 Mpc radius as a function of redshift z.As shown in Figure 14, the angle gradually decreases from about 50 • at z = 5 to almost zero at z = 0.This strongly suggests that the spatial distribution of the luminous satellites traces that of the dark matter on a 1 Mpc scale at low redshifts (z < 3).Before z = 0, the system has had sufficient time to evolve and settle into a more stable configuration.The abrupt change in the angle at z = 2 reflects the fact that a major accretion event happened around that time.It appears that after a merger event starts, the angle increases rapidly due to the complex dynamics induced by ongoing mergers, causing disruption and deviations from the previously formed alignment.As time progresses, the system undergoes a transitional period during which the combined influence of ongoing disruptions and the emergence of a more coherent structure might contribute to the observed variability in the angle.Eventually, around z = 0.7, as the system settles into a more stabilized configuration, the dynamics become dominated by the dark matter component, and the angle between the normal directions declines rapidly.This decline indicates a restoration of alignment between the satellite distribution and the underlying dark matter distribution.This trend of alignment near z = 0 suggests a compelling association between the dynamical evolution of satellite galaxies and the underlying DM sub-halo structure.

Summary and Conclusions
In summary, we utilized original high-resolution hydrodynamical simulations to examine the satellite system surrounding a Milky-Way-like galaxy.The incorporated baryonic physics in these original simulations enabled a comprehensive study of the spatial and dynamical properties of luminous satellites in such a system.The high mass resolution (∼3 × 10 4 M ⊙ for star particles) in the simulations makes the capture of the formation of low luminosity satellites with M > 10 5 M ⊙ possible in our study.Confining our study to the galaxy with the most Milky-Way like well-formed disk, we found, intriguingly, the presence of a DOS.While the existence of a DOS in a single galaxy does not address its likelihood of occurrence, it allowed an opportunity to investigate the properties and dynamics of the simulated DOS.
Our conclusions from this study are as follows: 1. We find that the properties of the fitted DOS planes are almost independent of the lower mass thresholds (proxy for luminosity) applied to the satellite samples, from 10 5 M ⊙ to 10 8 M ⊙ .

2.
The DOS formed by the luminous satellites and by the dark matter sub-halos have similar orientations.That is, the orientation of the DM DOS of Figure 5 and the most similar stellar DOS of Figure 4 are almost identical.

3.
Comparing our results with some previous observational studies, we find that the fitted DOS plane in our simulations can exhibit a significant inclination angle relative to the galactic disk which is slightly larger than the one observed in the Milky Way.However, the fitted DOS disk in the simulations is also thicker than that of the Milky Way.That is, the DOS aspect ratios given in Table 2 for the simulations range from 0.1 to 0.3, whereas the observed aspect ratio for the Milky Way DOS is ∼0.1-0.2.

4.
We also find that the "Zone of Avoidance" of the Milky Way may have some impacts on the fitted parameters.Indeed, when including the ZOA in the simulations, the properties (angles, distances, etc.) from the simulations more closely resemble those from observations.5.
Perhaps the most important conclusion from this study is the pronounced difference in the direction of the angular momentum of satellites in the DOS plane compared to the main galaxy.In the simulations, they are nearly perpendicular to each other suggesting that the satellite motion is consistent with infall from a larger scale.This is consistent with the observational data, indicating that, at least in our current study, the DOS is not rotationally supported.

6.
Our analysis shows that the anisotropic distribution of the satellites in the simulations can have characteristic properties such as angles and distances similar to the characteristics of the Milky Way's fitted DOS from observations.
Based upon these conclusions we speculate that the DOS derived from the positions of the nearby satellites within the main galaxy's virial radius may reflect the distribution of the satellites and DM sub-halos on a larger scale.Also, we suggest that the fitted DOS plane has been greatly affected by major accretion processes that happened in the evolution history of the main galaxy.The newly formed more luminous satellites accreted in these processes may have dominant roles in the formation of the observed DOS.
What is needed are much larger samples of higher-resolution realistic simulations of Milky-Way-like galaxies to substantiate that the observed DOS can be naturally formed in a standard ΛCDM cosmology [31].We will address this in a future work (cf.[91]).Further study of the DOS will also benefit from recent and forthcoming observational efforts including the potential discovery of new satellites (e.g., LSST), the measurements of the proper motions of satellites (e.g., Gaia [88]) and their line-of-sight velocities (e.g., DESI Milky Way Survey, [92]).The latter effort will also lead to the observation of extended features around sixteen MW satellites in the survey footprint [92].

Figure 1 .
Figure 1.The edge-on view of the log-scaled density projection map of the stellar, dark matter and gas components in the main galaxy system within its virial radius in the "C−4" simulation at redshift z = 0.

Figure 2 .
Figure 2. The face-on and edge-on views of the density projection map of the stellar component in the central part of the main galaxy in the "C−5" and "C−4" simulations at redshift z = 0.

Figure 3 .
Figure 3. Histogram of the number of satellites in various mass bins from the "C−4" simulation.

Figure 4 .
Figure 4.The edge-on views of the main galaxy's stellar disk within the virial radius in the "C−4" simulation at redshift z = 0.The dashed lines indicate the orientation of the plane of the DOS.The different DOS data sets are obtained by setting the minimum stellar mass of the satellites to be 10 8 M ⊙ , 10 7 M ⊙ , 10 6 M ⊙ and 10 5 M ⊙ as labeled.

Figure 5 .
Figure 5.The edge-on view of the main galaxy's dark matter halo and the DOS formed by all dark matter sub-halos within the virial radius in the "C−4" simulation at redshift z = 0.The dashed line indicates the orientation of the plane of the DOS.

Figure 6 .
Figure 6.The angle θ DOS between the normal direction of the DOS ⃗ n DOS and the angular momentum vector ⃗ L main of the stellar disk of the main galaxy (red line), and the host dark matter halo (black dotted line), respectively.This angle's variation is presented with respect to the minimum stellar mass criterion of the satellite sample within the virial radius in the "C−4" simulation at redshift z = 0.

Figure 8 .
Figure 8.The edge-on views of the main galaxy's stellar disk and the DOS formed by the satellites with minimum stellar mass of 10 8 M ⊙ , 10 7 M ⊙ , 10 6 M ⊙ and 10 5 M ⊙ within 1 Mpc radius in the "C−4" simulation at redshift z = 0.The dashed line indicates the orientation of the plane of the DOS.

Figure 9 .
Figure 9.The edge-on view of the main galaxy's dark matter halo and the DOS formed by all dark matter sub-halos within 1 Mpc radius in the "C−4" simulation at redshift z = 0.The dashed line indicates the orientation of the plane of the DOS.

Figure 10 .
Figure 10.The angle between the normal direction ⃗ n DOS of the DOS and the angular momentum vector ⃗ L main of the stellar disk of the main galaxy (red line) and to the host dark matter halo (black dotted line) as a function of minimum stellar mass in the satellites within 1 Mpc radius in the "C−4" simulation at redshift z = 0.

Figure 11 .
Figure 11.The Hammer-Aitoff projection of the full sky map of the main galaxy's stellar disk and the satellites with a minimum stellar mass of 10 5 M ⊙ within the virial radius (top) and 1 Mpc radius (bottom) (in RA(deg) and Dec(deg)) in the "C−4" simulation at redshift z = 0.The satellites in the "Zone of Avoidance" of the main galaxy are indicated by the dashed contour.

Figure 12 .
Figure 12.Angular relationships relative to the direction of the angular momentum of the DOS as a function of minimum stellar mass in the satellites within the virial radius in the "C−4" simulation at redshift z = 0.The line in solid blue represents the angle between the directions of the angular momenta of the DOS and the stellar disk of the main galaxy.The line in dashed blue represents the angle between the directions of the angular momenta of the DOS and the host dark matter halo.The red line represents the angle between the angular momentum vector of the DOS and the normal direction of the velocity dispersion vector, denoted as ⃗ σ DOS , of the DOS.The black line represents the angle between the normal direction of the DOS, represented by ⃗ n DOS , and the direction of the angular momentum vector of the DOS.

Figure 13 .
Figure 13.The edge-on views of the main galaxy's stellar disk and the DOS formed by the satellites with a minimum stellar mass of 10 5 M ⊙ within the virial radius (left column) and 1 Mpc radius (middle column) and all dark matter sub-halos within 1 Mpc radius (right column) in the "C−4" simulation at various redshifts.

Figure 14 .
Figure 14.The angle between the normal directions of the DOS ⃗ n DOS formed by the satellites with minimum stellar mass of 10 5 M ⊙ and all dark matter sub-halos within 1 Mpc radius as a function of redshift z in the "C−4" simulation.

Table 1 .
List of the characteristic parameters of the simulations and the main galaxies at redshift z = 0.

Table 2 .
List of the characteristic parameters of the data sets and the fitted DOS at redshift z = 0.