Planes of Satellites around Simulated Disk Galaxies. I. Finding High-quality Planar Configurations from Positional Information and Their Comparison to MW/M31 Data

The so-called "small-scale problems in ΛCDM" refer to the discrepancies between the predictions for dwarf galaxies in the standard cosmological model as first revealed by dark matter (DM)-only cosmological simulations, and the actual observed properties dwarfs present (see Bullock & Boylan-Kolchin 2017, for a review). Among them, the planar configurations satellite galaxies show around their hosts ("Planes of satellites problem," see Pawlowski 2018 for a review), observed in the local universe, have long been considered as one of the most challenging.

The high degree of anisotropy of Milky Way (MW) satellite positions, which appear forming a common plane approximately perpendicular to the Galactic disk, was noted several decades ago (Kunkel & Demers 1976; Lynden-Bell 1976) and first quantified by Kroupa et al. (2005) with the then-known 11 "classical"¹¹ satellites. With the addition of globular clusters, streams, and newly discovered fainter satellites (especially with SDSS; York et al. 2000), the anisotropy increased even more so, as these objects contributed further to a "vast polar structure (VPOS)" around the MW (Pawlowski et al. 2012). Anisotropy among Andromeda's (M31) satellites was first noted by Koch & Grebel (2006) and Metz et al. (2007), and later confirmed with a larger sample of satellites including those recently discovered with the Pan-Andromeda Archaeological Survey (PAndAS survey; Ibata et al. 2013). It was found that a majority of satellites are lopsided toward the MW's side (McConnachie & Irwin 2006). In addition, M31 satellites do not define only one main planar structure, like in the MW; instead a thin plane of satellites including approximately half of the satellite sample was singled out (Conn et al. 2013; Ibata et al. 2013; Pawlowski et al. 2013, hereafter "Ibata-Conn-14" plane). Finally, new star clusters and dwarf galaxy candidates have been recently found in other nearby galactic systems in the local universe like CenA or the M101 group of galaxies. Studies suggest as well an anisotropical 3D-spatial distribution (Tully et al. 2015; Müller et al. 2017, 2018).

Additionally, proper motion data has revealed that a high fraction of MW satellites present orbital angular momentum vectors mostly perpendicular to the Galactic disk axis (Metz et al. 2008; Pawlowski & Kroupa 2013; Fritz et al. 2018; Gaia Collaboration et al. 2018). In particular, Fritz et al. (2018) used recent proper motion data from GAIA DR2 to calculate the orbital poles of objects orbiting within 420 kpc around the MW. According to their results, approximately $\sim 40 \%$ (lower limit) of the confirmed MW satellites present orbital poles within an area of 10% of the sphere around the normal direction to the "VPOS," which they define as a co-orbitation criteria. Also, claims for co-rotation of satellites in the M31 "Ibata-Conn-14" plane have been made based on the direction of radial (i.e., line of sight) velocities (Metz et al. 2007; Ibata et al. 2013) as no proper motion data is yet available. It has been shown, however, that these are generally not a representative measure of the true 3D-velocities (Buck et al. 2016).

Flattened spatial configurations of satellites in the MW and M31 have been well studied and quantified using information from only the three-dimensional positions of the satellites. In particular, these positional analyses have used sampling techniques like "bootstrapping" (Metz et al. 2007) or the "four-galaxy-normal density plot" method (hereafter 4GND plot; Pawlowski et al. 2013) to statistically show the existence of predominant planar configurations of satellite positions in the MW and M31. These planes are then accurately characterized by their normal vectors, axis ratios, and root-mean-square heights (Δrms), computed from the eigenvalues of the Tensor of Inertia (ToI) plane-fitting technique (Metz et al. 2007; Pawlowski et al. 2013, 2015). In particular, the 4GND plot method consists in fitting planes to subsamples of four different satellites, and projecting the normal vectors on the sphere, creating a density map. The over-density regions that appear as a consequence of the accumulation of normal vectors broadly point in the normal direction to a predominant planar configuration of satellites. Following this method, Pawlowski et al. (2013) detected the specific subsamples of satellites that mostly contribute to planar configurations in the MW and M31, defining the "VPOS-3" plane of satellites in the MW and the "GPoA" plane in M31.

In a recent paper, Santos-Santos et al. (2019, hereafter Paper I), have extended the 4GND plot method to allow for an identification, systematic cataloging, and more detailed quality study of the planar configurations of satellites in the MW and M31 systems. Rather than deriving a unique plane of satellites per over-density in the 4GND plot found with the previous method, the extension yields a collection of planes of satellites, each with an increasing number of members. In this way, it is possible to identify the highest-quality planes in terms of the ToI parameters and the number of satellites considered. New to previous findings, in Paper I, it was shown that two distinct planes of satellites are present in M31: the "GPoA" and another plane (labeled "M31-2-18" in Paper I). The two planar structures present very similar characteristics and are, interestingly, oriented perpendicularly to each other.

Since the very advent of these discoveries, theoretical studies have tried to assess the frequency of planar satellite configurations like those observed in cosmological simulations within the ΛCDM paradigm. These studies have mostly made use of large-volume DM-only simulations and pay attention to determining the significance¹² of the observed planes. In particular, Libeskind et al. (2009), Wang et al. (2013), Bahl & Baumgardt (2014), and Cautun et al. (2015) have analyzed different versions of the large-volume Millenium DM-only simulation (Springel et al. 2005; Boylan-Kolchin et al. 2009), populating subhalos with galaxies following semi-analytic models. Using different methods for plane-identification, they find that planes of subhalos as thin and even thinner than the ones observed in the MW and in M31 are expected in ΛCDM. However, they acknowledge that these are not the mean case found, but only consistent with the tail of the predicted flattening distribution (see Pawlowski et al. 2014, for a different view). On the other hand, Cautun et al. (2015) show that not accounting for the "look-elsewhere effect" results in an important overestimation of the significance of planes in the MW and M31 systems of factors of 30 and 100, respectively. Indeed, according to that work, ∼10% of MW-like mass halos in ΛCDM simulations have planes of satellites that are more prominent than those observed in the MW or M31 systems, presenting a large diversity when characterized by their thickness and number of satellites. On the other hand, a different approach has been followed by Buck et al. (2015). Instead of a large volume, they use several DMO zoom-in simulations to show that a thin plane, as the one around M31 with 15 satellites (Ibata et al. 2013), is not a challenge for the ΛCDM paradigm. However, neither the VPOS-3 (with 24 satellites) or VPOSall (27) planes in the MW are recovered in their analysis.

While some insight has been gained from collisionless N-body simulations, these experiments do not allow the formation of galactic disks. Having well-behaved massive MW-like simulated disks, however, may be critical to the planes of satellites issue (see Ahmed et al. 2017; and for a different perspective, Pawlowski et al. 2019). Indeed, the dynamical effect of a live disk potential on satellite planes could change the frequencies alluded to above, due, for example, to the torques that satellites suffer from galaxy disks—except when they are on planar or polar orbits, or far away from the disk plane (Danovich et al. 2015; Welker et al. 2018). Another effect is that galaxy disks preferentially destroy satellites on radial orbits that pass near them (Garrison-Kimmel et al. 2017; Sawala et al. 2017). Riley et al. (2019) show that massive disk potential effects cause the velocity anisotropy parameter, β, to decrease as compared to less massive disks. Including these relevant effects in planes of satellites studies is thus necessary for a fair comparison with results from the MW/M31 disk galaxy systems. Another relevant point is that there are $\gtrsim 30$ confirmed satellites in the MW and M31 (McConnachie 2012): a proper comparison demands as well the analysis of simulated disks surrounded by roughly as many resolved satellites.

Meeting all of the previous requirements is currently a situation not found in large-volume hydrodynamical simulations. For example, in their analysis of the EAGLE simulation, Shao et al. (2019) analyze planes of 11 satellites around central galaxies of any morphology and compare them to the MW "classical" plane—thus, the motivation to analyze zoom-in high-resolution hydrodynamical simulations in order to study planes of satellites.

A few studies exist using zoom-in hydro-simulations. In Gillet et al. 2015 and Ahmed et al. 2017, the method used for optimal plane searching has consisted in fitting a vast number of pre-defined planes with a given constant thickness to the satellite sample. These planes are forced to pass through the center of the main galaxy (see also, Buck et al. 2015). In Libeskind et al. 2007, Maji et al. 2017, and Garaldi et al. 2018, the ToI method has been used, mainly on the subsample of the 10-11 most massive satellites. In this way, planes of satellites have been found. The best ones, however, are not thin and populated enough so as to reproduce the VPOS-3 and VPOSall planes in the MW (see, e.g., Maji et al. 2017) or even the Ibata-Conn-14 plane in M31 (see, e.g., Gillet et al. 2015).

In this project, we further contribute to the study of zoom-in simulations by introducing a novel systematic and detailed plane searching technique to study the properties and quality of planes of satellites. This methodology is suitable to analyze both observational data and simulations. We aim to gain further insight on the properties of planes of satellites one can find in well-behaved massive disk galaxies (like the MW and M31), thus formed in zoom-in high-resolution cosmological hydrodynamical simulations, where a high enough number of resolved satellites can be identified. This is an important step, whose outcomes are needed prior to attempting any physical interpretation of the origin and/or evolution of planes.

Specifically, in this paper, we develop a detailed analysis of planar configurations of satellites from positional information by applying, for the first time, the 4GND plot method (Pawlowski et al. 2013) and its extension (introduced in Paper I) to simulations. For the reasons explained above, we focus on a set of zoom-in hydrodynamical simulations where well-behaved MW-like disk galaxies and their satellite systems form. Moreover, as the effects of galaxy companions either in binary systems or groups add complexity to this analysis, we consider only isolated galaxies. As mentioned previously, this method allows us to identify the constituent satellites forming planar spatial configurations and analyze their quality in terms of their population ${N}_{\mathrm{sat}}$ (or, equivalently, the fraction of satellites ${f}_{\mathrm{sat}}$) and the ellipsoid of concentration axes ($a,b,c$, with $a\gt b\gt c$). Planes of high quality are those with a high ${f}_{\mathrm{sat}}$ and a low c/a, meaning they are populated and thin. These analyses are done over the entire galaxy's evolution after halo virialization. This allows us to glean important insights into the different kinds of planar structures one can find in simulations of disk galaxies, and how they compare to the observed planes of satellites.

We note that in the "planes of satellites problem," the persistence issue (i.e., is there a same group of satellites that is spatially distributed in a planar-like configuration across time?) is closely related to the kinematical character of planes. Therefore, in a forthcoming paper (I. Santos-Santos et al. 2020, in preparation, Paper III), the full six-dimensional phase-space information on satellite motions will be used to carry out kinematically based analyses as an optimal methodology to address satellite plane persistence.

The paper is organized as follows. In Section 2, we introduce the simulations analyzed, while in Section 3, their corresponding satellite samples and some of their properties are presented. Section 4 describes the method used for positional plane searching and plane quality analysis. Results obtained are reported in Sections 5 and 6. In Section 7, we assess the possible co-orbitation of satellites within the planes they form. In Section 8, we discuss the implications of our results. Finally, Section 9 summarizes the results and exposes the conclusions reached.

We have chosen to analyze planes of satellites orbiting around isolated, simulated galaxies that resemble the MW system. In particular, we demand that the simulation meets the following requirements:

• a)  
  to contain a central galaxy with a thin gaseous and stellar disk at redshift $z\sim 0$, with a large radial extent ($R\,=15\mbox{--}25$ kpc).

We note that very thin disks are not that common in hydro-simulations yet.

Moreover, as merger events could destabilize the galaxy+satellites system, complicating clean plane detections as well as the possibility of reaching conclusions concerning the origin of these planes, we as well require:

• (b)  
  an overall quiet assembly history, i.e., free of major merger events after virialization.

This is in line with the current understanding of the MW's disk formation and accretion history (Belokurov et al. 2018; Helmi et al. 2018).

To allow for a proper statistical comparison of the results obtained with those coming from observations of the MW and M31 at z = 0, which harbor (at least) around 30 satellites each,¹³ the system must also:

• (c)  
  host a numerous (∼30) satellite population around the central galaxy.

Finally, in order to accurately compute the center of mass and the orbital angular momentum of the baryonic component of a simulated satellite, and rely on it as physical,

• (d)  
  we demand satellite objects must include more than 50 baryonic particles.

We have pre-analyzed a set of different zoom-in cosmological hydro-simulations, finding among them two that reach the previous prerequisites: Aquarius-C${}^{\alpha }$ and PDEVA-5004. Both simulations follow the "zoom-in" technique but make use of very different initial conditions, codes, and physics prescriptions. This fact will allow us to reach conclusions that are independent of simulation modeling.

2.1. Codes and Host Galaxies

3.1. Identification

The samples of satellites in each simulation have been selected following the same criteria, that is, to take all objects with stars (${M}_{\star }\gt 0$) within a distance of 350 kpc from their host that are bound to the host galaxy, with a resolution limit of presenting at least 50 baryonic particles (${M}_{\mathrm{bar}}\approx 1\times {10}^{7}{M}_{\odot }$). This selection has been made at redshift $z\sim 0.5$, to include satellites that may end up accreted by the disk at z = 0. To prove if a given object is indeed a satellite (i.e., is bound to its host), we have computed its orbit. This fixed sample of satellites has then been followed back and forth in time.¹⁴

For a proper comparison with MW results, we take into account Galactic obscuration—which prevents us from observing satellites orbiting in the plane of the disk of the MW, by applying an observational bias to the simulated satellite sample at each timestep. Following Pawlowski (2016), we have chosen it to hide objects with projected positions on the sphere at latitudes $| b| \lt 12^\circ$ (angular distance as measured from the plane of the galaxy disk). This is a first approximation to try to mimic the obscuration effects of the MW's disk; however, we acknowledge that a more precise model to account for Galactic obscuration should depend on satellite distances and luminosities.

The selection of satellites in Aq-C${}^{\alpha }$ has been done using a Friends-of-Friends algorithm to identify structures and then the SubFind halo finder (Springel et al. 2001) to construct subhalo catalogs at each timestep. Particle IDs have been used to trace back in time the selected satellites. The tool used for the selection of satellites in PDEVA-5004 has been IRHYS (by H. Artal, under development). This visualization and analysis tool permits the selection of objects (i.e., satellites) as sets of particles and enables us to trace them back and forth in time.

A total number of ${N}_{\mathrm{tot}}$ = 30 (35) satellites have been detected in Aq-C${}^{\alpha }$ (PDEVA-5004) at selection time ($z\sim 0.5$), of which 25 (27) survive until the last analyzed timestep, respectively. In Figure 1, such numbers are plotted as a function of the universe age ${T}_{\mathrm{uni}}$. ${N}_{\mathrm{tot}}$ changes because satellites disappear as they are accreted by the central disk galaxy. Also, satellites are not considered during periods where they orbit at distances $\gt 450\,\mathrm{kpc}$ (this happens with a couple of backsplash galaxy cases). In the case when obscuration in the plane of the disk is considered ("bias"), ${N}_{\mathrm{tot}}$ varies additionally because satellites go into and out of the avoidance volume.

Zoom In Zoom Out Reset image size

3.2. Satellite Radial Distances and Their Distributions

The suite of Aq-C${}^{\alpha }$ and PDEVA-5004 satellites presents different properties and evolutionary histories that reflect in a variety of orbits. We find some satellites that progressively lose angular momentum and are eventually accreted by the disk, some that follow orbits where successive apocenter and pericenter distances do not show important variations, and some backsplash galaxies that have just recently been captured by the halo and orbit at long distances occasionally even outside the virial radius. Interestingly, when all radial distance histories are plotted together, a coincidence of pericenters is clear at certain moments. In this way, as the simulations evolve, there are moments of maximum spatial expansion and moments of a maximum compression of the satellite systems.

These effects can be observed when analyzing the evolution of the radial distribution of satellites with cosmic time. Figure 2 (top panel: Aq-C${}^{\alpha }$, bottom panel: PDEVA-5004) shows the fraction of the total number of satellites within a certain distance from the center of the main galaxy, compared to the MW and M31 distributions at z = 0.¹⁵ Two colored lines are shown per panel, which represent the results obtained using all of the satellites, or taking into account the observational obscuration bias, at each timestep. The distributions change with time, showing periods where there is a higher concentration of satellites at shorter distances (corresponding precisely to the moments of maximum collective approach) and others where there is a higher expansion of the system. As an example, the PDEVA-5004 system (in the case where all satellites are considered, i.e., "no bias") is more compact (higher ${f}_{\mathrm{sat}}$ within 100 kpc) at ${T}_{\mathrm{uni}}=8.9\,\mathrm{Gyr}$, and more expanded at ${T}_{\mathrm{uni}}=11\,\mathrm{Gyr}$. On the other hand, in the Aq-C${}^{\alpha }$ system, it is not until ${T}_{\mathrm{uni}}\sim 9\,\mathrm{Gyr}$ that the complete sample of satellites is within a distance of ∼350 kpc. A moment of maximum compactness is ${T}_{\mathrm{uni}}\sim 10.3\,\mathrm{Gyr}$. Curiously, PDEVA-5004's satellite radial distribution resembles very well that of the MW for a long period of time, while Aq-C${}^{\alpha }$'s is similar to that of M31 at ${T}_{\mathrm{uni}}\sim 10\,\mathrm{Gyr}$. These resemblances are kept when using ${N}_{\mathrm{sat}}$ instead of ${f}_{\mathrm{sat}}$ in the case that total satellite numbers of simulations and observations are matched (see Sections 8.1 and 8.2). We explore the effect of radial compactness on the quality of planes of satellites in Section 8.3.

Zoom In Zoom Out Reset image size

3.3. Mass Distribution of Satellites

To search for planar positional configurations of satellites in our simulations we have followed the 4-galaxy-normal density plot (4GND plot) method presented in Pawlowski et al. (2013), extended as explained below (see also Paper I). This method allows us to determine if there is a subsample out of a given sample of ${N}_{\mathrm{tot}}$ satellites that defines a planar arrangement in terms of the outputs of a fitting technique based on the Tensor of Inertia (ToI; Metz et al. 2007; Pawlowski et al. 2013). Satellite planes are searched for through a regression method that minimizes orthogonal distances from the points to the optimal plane solution. Apart from the plane (or equivalently, its normal vector ${\boldsymbol{n}}$), the outputs of the regression can be characterized by the following parameters in terms of the corresponding ellipsoids of concentration (Cramér 1999):

• (i)  
  ${N}_{\mathrm{sat}}$: the number of satellites in the subset (or, the fraction of satellites it involves ${f}_{\mathrm{sat}}\equiv {N}_{\mathrm{sat}}/{N}_{\mathrm{tot}}$);
• (ii)  
  c/a: the ellipsoid short-to-long axis ratio;
• (iii)  
  b/a: the ellipsoid intermediate-to-long axis ratio;
• (iv)  
  Δ rms: the root-mean-square thickness perpendicular to the best-fitting plane;
• (v)  
  ${D}_{\mathrm{CG}}$: the distance from the center-of-mass of the central galaxy to the plane.

These outputs are used to quantify the quality of the best-fitting structures to a subsystem of ${N}_{\mathrm{sat}}$ satellites. First of all, assuming $c/a\lt 1$, b/a indicates whether the distribution is planar ($b/a\sim 1$) or filament-like ($b/a\ll 1$). High-quality planes are those that involve many satellites and are thin, therefore demanding high ${f}_{\mathrm{sat}}$ and low c/a (or equivalently low Δrms, a quantity that most often is correlated with c/a once the system acquires its stable size). Low ${D}_{\mathrm{CG}}$ planes pass near the disk center, a requirement asked of a physically consistent satellite system when its gravitational potential minimum lies approximately at this center. Finally, ${\boldsymbol{n}}$ allows us to visualize the plane orientation with respect to a given reference frame, for example, the host galaxy disk plane. At the end of this section, the quantification of plane quality, as well as how to compare the qualities of two or more planes, is described in more detail (see also Paper I).

4.1. Method: Four-galaxy-normal Density Plots

We have applied the 4GND plot method to each timestep of the two simulations. A thorough description follows (see also Section 2.4 in Pawlowski et al. 2013, and Paper I).

• 1.  
  A plane is fitted (using the ToI method) to the positions of every combination of four different satellites taken from the total sample of ${N}_{\mathrm{tot}}$ satellites. As three points always define a plane, four is the lowest possible amount to take into consideration under the condition of making the number of combinations high enough to get a good outcome signal.¹⁹
• 2.  
  The axes sizes $a,b,c$ and normal vector directions (i.e., four-galaxy-normals) of the planes fitted to each combination of four satellites are stored. Then, all of the four-galaxy-normals are plotted in a galactocentric coordinate system such that the central disk's spin points toward the south pole, and a density map (i.e., 2D-histogram) is drawn from their projections on a regularly binned sphere with ${N}_{\mathrm{bin}}$ bins. Spherical projections are shown with Aitoff diagrams in Galactic (longitude l, latitude b) coordinates in a $l=[-90^\circ ,+90^\circ ]$ projection because opposite normal vectors are equivalent. As in Pawlowski et al. (2013), each normal is weighted by $w=\mathrm{log}\left(\tfrac{a+b}{c}\right)$ to emphasize planar arrangements of satellites over filament-like or spherical-like ones. In these plots, over-densities (or density peaks, i.e., areas of four-galaxy-normal vector accumulation) are signaling groups of four satellites contributing to the same dominant planar space-configuration. As illustrated in Figures 3 and 4, we show examples of 4GND plots for the Aq-C${}^{\alpha }$ and PDEVA-5004 simulations, respectively. In some cases, these over-dense areas are more extended, and in others, more concentrated. (Note that the expectation from a random distribution of satellites is a density map with equal density in each bin.)
• 3.  
  Density peaks are differentiated and isolated. A set of ${N}_{\mathrm{peaks}}$ density peaks are selected around the highest-density bins of the 4GND plots, with the requirement that they are separated more than 15° from the center of all of the (N_peaks − 1) over-densities. The specific peak location in ($l,b$) is given by the center of the corresponding high-density bin.²⁰
• 4.  
  We determine how much a given satellite s contributes to the pth specific peak (i.e., its respective contribution-number, ${C}_{p,s}$, with p = 1, ..., ${N}_{\mathrm{peak}}$ and s = 1, ..., ${N}_{\mathrm{tot}}$). To this end, we select all of the four-galaxy-normals that fall within an aperture angle of 15°²¹ around the pth peak location. Each of the four satellites contributing to these four-galaxy-normals is counted as contributing the four-galaxy-normal's weight to the peak. This has been normalized using ${C}_{N,\mathrm{all}}$, the total weighted number of four-galaxy-normals, including those that are not within 15° of some peak center, such that the sum ${\sum }_{p,s}\,{C}_{p,s}=1$. Such normalization is necessary for a meaningful comparison of results at different timesteps (where ${N}_{\mathrm{tot}}$ varies) and, also, with observational data. At fixed p, ${C}_{p,s}$ is high when satellite s contributes to many of the four-galaxy-normals laying within 15° of the peak center. We note that over-densities that are located close to each other on the sphere generally share many of the satellites that contribute most to four-galaxy-normals within 15° of the peak. However, the peak isolation criterion used in the previous step cures our peaks of such redundancies.
• 5.  
  For a given peak p, we order all satellites by decreasing ${C}_{p,s}$ to it. This is done for all of the isolated peaks found in a given 4GND plot, and for each of them, we obtain an ordered list of contributing satellites. Examples of such lists are shown in Figure 5 for two different peaks defined at a given timestep in PDEVA-5004. The x-axis shows satellite IDs in decreasing ${C}_{p,s}$ order (see y-axis values).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. An Extension to the Method

5.1. Four-galaxy-normal Density (4GND) Plots

5.2. Satellite Contribution-numbers to Peaks

5.3. Peak Strength Analysis

The ${C}_{p}$ peak strengths of Peak 1 and Peak 2 (i.e., C₁ and C₂) for the MW and M31 are given in Table 1. Errors are 1σ deviations over 100 random realizations of their radial distance uncertainties, as explained in Paper I. These are especially large in the case of M31 satellites; as a result, peaks in its 4GND plot are more blurred as compared to the MW ones. This contributes to lower C₁ and C₂ values in M31.

Table 1.  Peak Strengths C_p for the Main Two Peaks Found in the MW and M31 4GND Plots (see Paper I)

	${C}_{1}\pm \sigma$ (%)	${C}_{2}\pm \sigma$ (%)
MW	22.92 ± 0.26	14.31 ± 0.20
M31	10.53 ± 0.62	10.52 ± 1.62

Note. Peak strength is computed as ${C}_{p}\equiv {{\rm{\Sigma }}}_{s}{C}_{s,p}$, where ${C}_{s,p}$ is normalized to the total weighted number of four-galaxy-normals (see Section 5.3). Results shown are means and 1σ standard deviations calculated from 100 random realizations using the radial distance uncertainties. Peaks #1 are the strongest ones and Peaks #2 follow in strength.

Download table as:  ASCIITypeset image

In the upper panels of Figure 6, we present the value of C₁ at each timestep for Aq-C${}^{\alpha }$ and PDEVA-5004. C₁ fluctuates, reaching values that can be even higher than those of the MW or M31 at z = 0. In general, the application of the observational obscuration bias enhances the strength value of the main peak.

Zoom In Zoom Out Reset image size

The universe ages ${T}_{\mathrm{uni}}$ where the respective values of C₁ (in the case where all satellites are considered, "no bias") have maxima (minima) are marked by green (magenta) vertical lines. These time intervals of local maxima and minima have an average duration of 0.5–1 Gyr (consistent with the values Shao et al. (2019) find in their analyses of the EAGLE simulation). This has been estimated from their FWHM, where we take the mean C₁ as the floor value. These periods will be related with plane quality in the next sections.

Another interesting possibility that the peak strength ${C}_{p}$ allows is to determine the number of peaks with strengths above given thresholds or within given intervals, at different universe ages. This is shown in the lower panels of Figure 6, with respect to the strengths of the peaks labeled as #1 in the MW and M31, i.e., ${C}_{1,\mathrm{MW}}$ = 22.9% and ${C}_{1,{\rm{M}}31}$ = 10.5%. At given times, there are a few peaks encompassing a high % of four-galaxy-normals (high C₁) that then break into several different peaks with lower strengths. These later collimate into high ${C}_{p}$ peaks again. That is, the number of peaks with ${C}_{p}$ within given peak strength intervals fluctuates with time. We recall that the background is also accounted for to normalize the peak strengths at fixed timesteps.

Summing up, as measured with C₁, the peak strengths in observations and simulations are consistent within given time intervals (∼3 periods with peaks stronger than the MW in the 'bias' case, and 1 such period in the 'no-bias' case.). Regarding the number of peaks, we see that when C₁ takes high values, the number of weak peaks decreases and that of stronger peaks increases. This happens especially at the last timesteps analyzed. In particular, the number of strong peaks (${C}_{1}\gt {C}_{1,\mathrm{MW}}$) reaches 1–3 at most, (occurring at late times), consistent with observations. Again, we can see that the obscuration bias favors the emergence of strong peaks in both simulations.

6.1. Comparing to the MW and M31 Satellite Systems

6.2. Quality of Simulated Planes in Terms of the Satellite Fraction Involved

6.2.1. b/a versus f_sat

Concerning the application of the method to the simulation data, we first address the planar ($b/a\sim 1$) or filamentary ($b/a\ll 1$) character of the best-fitting structures found with the ToI analysis, where b/a is the intermediate-to-long axis ratio in the ToI scheme. As an illustrative example of our results, in Figure 7, we plot b/a versus ${f}_{\mathrm{sat}}$ for the main density peaks found in Aq-C${}^{\alpha }$ at different timesteps. In this figure, and in the following ones that compare the changes of a ToI output with ${f}_{\mathrm{sat}}$, each panel corresponds to a given timestep. Lines of different colors stand for the collections of planes of satellites associated with the respective peaks numbered with the same color coding in Figure 3. Based on Figure 6 and the 4GND plots shown in Figures 3 and 4, the consideration of a number of ${N}_{\mathrm{peak}}$ = 3 peaks for Aq-C${}^{\alpha }$ and ${N}_{\mathrm{peak}}$ = 5 peaks for PDEVA-5004 seems a reasonable choice ensuring the exploration of all possibly relevant planar configurations.

Zoom In Zoom Out Reset image size

Observational data results are shown as gray lines and points. Points show the specific values for MW/M31 observed planes of satellites at z = 0 mentioned in the literature (i.e., MW: classical, VPOS-3, VPOSall; M31: Ibata-Conn-14, GPoA; see values in Table 1 of Paper I). When comparing to the MW, we show simulated results where the obscuration bias is applied; when comparing to M31, we show results considering all satellites.²⁵

We find no filamentary structures; in fact, at all timesteps, a planar structure exists whose b/a is larger than the observational case, at all ${f}_{\mathrm{sat}}$. The general behavior of b/a, both for Aq-C${}^{\alpha }$ and PDEVA-5004, and in the "bias" and "no bias" cases, is that b/a changes only slightly when new satellites are added to the fit, giving rise to wide ${f}_{\mathrm{sat}}$ intervals where b/a is almost constant. This means that the planar character of the spatial configuration of satellites does not depend very much on the number of satellites involved. This behavior is also found in the MW and M31 (see gray lines in Figure 7).

6.2.2. c/a versus f_sat

Having confirmed that the structures found in our simulations are indeed planar, we can proceed with the study of the quality of such planes through c/a. As explained previously, quality is assessed by a two-parameter notion (c/a, ${f}_{\mathrm{sat}}$) such that at a given ${f}_{\mathrm{sat}}$, the plane with the lowest c/a presents the highest quality. In particular, we define that a strong consistency exists between plane collections from simulations and observations when there is one colored line from simulations with similar or lower c/a values than that of the MW/M31 gray line at all ${f}_{\mathrm{sat}}$. A weaker condition refers to consistency between an observed plane and one detected in simulations with the same particular ${f}_{\mathrm{sat}}$. In this case, the peak assuring consistency between data and simulations can vary from ${f}_{\mathrm{sat}}$ to ${f}_{\mathrm{sat}}$.

We have carried out the analysis of plane population and thickness (c/a versus ${f}_{\mathrm{sat}}$) for the main peaks found in Aq-C${}^{\alpha }$ and PDEVA-5004, comparing to both the MW and M31. As examples of the results obtained, in Figure 8, we present the results for Aq-C${}^{\alpha }$ versus MW (top panel; the observational obscuration bias is applied) and PDEVA-5004 versus M31 (bottom figure; all satellites are considered).

Zoom In Zoom Out Reset image size

Independently of applying the observational obscuration bias or not, both Aq-C${}^{\alpha }$ and PDEVA-5004 present thin and highly populated planes at all timesteps. In general, c/a is low ($\lesssim 0.3$) for all peaks at all timesteps when including up to $\sim 80 \%$ of satellites. This is already proving the oblate spatial distribution of the entire satellite population in both simulations. In particular, at all timesteps and in both simulations, there are planar structures compatible (in terms of c/a and ${f}_{\mathrm{sat}}$) with the M31 Ibata-Conn-14 plane and the MW classical plane. The GPoA value is also recovered in almost all timesteps. The strong consistency condition is met in many cases. For example, we can find very similar or higher-quality planar structures than that of the MW in Aq-C${}^{\alpha }$ ("bias" case) at ${T}_{\mathrm{uni}}$ = 8.6, 9.2 and 10.8 Gyr; and we find similar or higher-quality planar structures than that of M31 in PDEVA-5004 at ${T}_{\mathrm{uni}}$ = 4.9, 9.6, 10.8, 13.4 and 13.7 Gyr.

6.2.3. Evolution of Plane Quality with Cosmic Time

A more compact and clearer way of presenting the results on plane quality analysis showed in Figure 8 is looking at the "best" plane found at each timestep including a fixed fraction of the total number of satellites. This best plane is selected as the one with the lowest c/a at fixed ${f}_{\mathrm{sat}}$, which can be easily read from Figure 8 (note that this best-quality plane does not necessarily correspond to the same peak as ${f}_{\mathrm{sat}}$ changes).

In Figure 9, we show the properties of the "best" planes of satellites found at each timestep. In particular, we focus on c/a, Δrms height, and the inclination of the plane relative to the disk (latitude angle). Different shades of blue stand for planes with different ${f}_{\mathrm{sat}}$ = 30%, 50%, 70%, 90%. We show results with and without applying the observational obscuration bias in the left and right panels, respectively. In the right panels, green (magenta) lines mark the values of the universe ages where C₁ maxima (minima) appear in Figure 6. For comparison, the results for the best planes in the MW and M31 at z = 0 with the same ${f}_{\mathrm{sat}}$ are also shown as horizontal lines. These values, together with those of simulated results averaged over the last 1 Gyr analyzed, are given in Table 2. Note how at high ${f}_{\mathrm{sat}}$ (70% and 90% lines), M31 presents very large c/a and Δrms, due to the system configuration in two almost perpendicular planes (see Paper I).

Zoom In Zoom Out Reset image size

Table 2.  Plane Parameters of the Best (i.e., with Lowest c/a) Planes Found in Aq-C${}^{\alpha }$ and PDEVA-5004 Including a Fixed Fraction (X%) of Satellites

		Aq-C${}^{\alpha }$ (last 1 Gyr)			PDEVA-5004 (last 1 Gyr)			MW (z = 0)			M31 (z = 0)
		c/a	Δrms	${D}_{\mathrm{CG}}$	c/a	Δrms	${D}_{\mathrm{CG}}$	c/a	Δrms	${D}_{\mathrm{CG}}$	c/a	Δrms	${D}_{\mathrm{CG}}$
			(kpc)	(kpc)		(kpc)	(kpc)		(kpc)	(kpc)		(kpc)	(kpc)
Bias	30%	0.03	5.30	23.70	0.03	2.46	12.43	0.07 ± 0.001	9.54 ± 0.15	13.23 ± 0.15
	50%	0.07	12.70	23.77	0.06	6.27	11.96	0.10 ± 0.001	10.31 ± 0.12	14.61 ± 0.11
	70%	0.14	23.82	23.07	0.15	13.12	6.70	0.14 ± 0.002	12.67 ± 0.18	15.72 ± 0.15
	90%	0.26	41.05	12.60	0.25	21.11	4.13	0.21 ± 0.002	19.39 ± 0.19	10.46 ± 0.15
No bias	30%	0.04	6.72	27.92	0.03	3.45	7.76				6.24 ± 0.003	1.21 ± 0.35	6.17 ± 0.57
	50%	0.09	15.46	16.53	0.09	9.22	5.30				9.56 ± 0.002	2.05 ± 0.20	6.29 ± 0.58
	70%	0.17	26.02	17.73	0.15	13.14	5.83				25.83 ± 0.008	3.30 ± 0.83	2.62 ± 2.04
	90%	0.28	40.98	14.42	0.23	19.68	3.92				60.06 ± 0.009	8.02 ± 1.18	10.17 ± 6.27

Note. We show c/a, Δrms, and ${D}_{\mathrm{CG}}$. Values are averaged over the last Gyr of the corresponding analysis period (see Figure 9). Fractions for the MW and M31 have been calculated relative to a sample size of ${N}_{\mathrm{tot}}$ = 27 and 34 satellites, respectively. The results shown are the mean values and standard deviations resultant of 100 random realizations of radial distances within the observational uncertainties.

Download table as:  ASCIITypeset image

In terms of c/a and ${f}_{\mathrm{sat}}$, both Aq-C${}^{\alpha }$ and PDEVA-5004 simulations present high-quality planes. The best planes of satellites take c/a values that change with time, reaching at some timesteps, and particularly near the respective last periods analyzed, values compatible with those in the MW and M31 at z = 0 involving the same fraction of satellites.

We note that the fact that quality in our simulations increases toward low redshift is in contrast to Shao et al. (2019) findings, who report thinner planes of satellites at early times in EAGLE. On the other hand, Figure 9 shows that biased c/a results are systematically somewhat lower than non-biased ones in both simulations. This occurs because when applying the bias, we are removing a fraction of the volume where satellites can be considered, which prevents plane-thickening. This result may indicate that the quality of the MW's planar structure of satellites can appear artificially enhanced because of Galactic obscuration.

In terms of Δrms, PDEVA-5004 reflects similar results and the same fluctuation patterns seen with c/a. Especially at low redshifts, very low Δrms heights are found. In Aq-C${}^{\alpha }$, despite the low c/a values, we find larger Δrms values. This is because Δrms is a dimensional quantity that therefore depends on the overall size of the system at issue. This parameter very clearly decreases as the system evolves: Aq-C${}^{\alpha }$ is still settling its size until ${T}_{\mathrm{uni}}\sim 9$ Gyr. At the last moment of our analysis (${T}_{\mathrm{uni}}=11.5\,$ Gyr), the Δrms heights of planes are generally compatible with their observed counterparts at z = 0 (except for the "biased" (versus MW) results involving ${f}_{\mathrm{sat}}$ 90% and 70%, and the "non-biased" ones (versus M31) involving a 50%). However, we note that Aq-C${}^{\alpha }$ has still 2 Gyr to reach z = 0, and the system could still evolve toward a lower Δrms value.

In the third rows of each panel of Figure 9, we plot the angle formed at each timestep by the normal vector to the plane of satellites and the galaxy's disk plane (we use Galactic latitude angle, with $| b| \,\in \,$[0, 90°]). We can see that while low ${f}_{\mathrm{sat}}$ curves show a fluctuating behavior with ${T}_{\mathrm{uni}}$, the fluctuation level decreases as ${f}_{\mathrm{sat}}$ increases, and, finally, almost no fluctuations show up at ${f}_{\mathrm{sat}}$ = 90%. An important variation in this angle is an indication that the identities of the satellite members of planes with given ${f}_{\mathrm{sat}}$ have changed. Therefore, these results are indicating that the satellite members of the best-quality planes change quite a lot at low ${f}_{\mathrm{sat}}$, while at ${f}_{\mathrm{sat}}$ = 70% or even 50%, these identities are kept to an important extent.

Moreover, at times when c/a reach their minima (and the main peak strength C₁ reach their maxima) the latitude angle in both simulations is small and sometimes close to 0° (except for PDEVA-5004 at $z\approx 0$); that is, satellite planes are nearly perpendicular to the galaxy's disk. This result suggests that the best quality of satellite planes is, in many cases, reached at time intervals when no (or rather low) galaxy disk torques act upon the satellites belonging to the plane that best fits the satellite set (see, e.g., Danovich et al. 2015).²⁶

Finally, we focus on the distances ${D}_{\mathrm{CG}}$, or offsets, from the center of the galaxy to the previously presented best planes found in Aq-C${}^{\alpha }$ and PDEVA-5004. Contrary to other plane-identification methods used in simulation studies (see, e.g., Buck et al. 2015; Gillet et al. 2015; Ahmed et al. 2017), in the 4GND plot method, planes are not constrained to pass through the center of the main galaxy. Table 2 shows the offsets to the best planes averaged over the last 1 Gyr with ${f}_{\mathrm{sat}}$ = 30%, 50%, 70%, 90%. For comparison, the MW VPOS-3 (${f}_{\mathrm{sat}}$ = 24/27 = 88%) and M31 GPoA (${f}_{\mathrm{sat}}$ = 19/34 = 56%) planes present an offset of 10.4 and 1.3 kpc from the center of the MW and M31, respectively (see Table 3 in Pawlowski et al. 2013). Also, the collection of planes defined by the second peak in the 4GND plot of M31 presents distances ${D}_{\mathrm{CG}}$ between ∼15 and 35 kpc (see Paper I). Compared to these, the plane offsets measured in both simulations have reasonable values, passing close to the center of the main galaxy.

This section reveals that there are indeed preferential planar configurations of satellites at given moments in both Aq-C${}^{\alpha }$ and PDEVA-5004 simulations. These planes are thin and highly populated, compatible on average with all characteristics of the observed planar structures found in the MW and in M31 at z = 0, and even defining higher-quality planes at particular given times. A rough comparison to the planes of satellites reported in previous studies with zoom-in hydro-simulations that consider a moderately large sample of satellites (e.g., Gillet et al. 2015; Ahmed et al. 2017; Maji et al. 2017) focusing on the c/a, Δrms, and ${N}_{\mathrm{sat}}$ parameters when provided by the references, strongly suggests that the planes found in Aq-C${}^{\alpha }$ and PDEVA-5004 have, within some time intervals, a higher quality and reveal a higher degree of spatial ordering in the satellite distribution. However, this comparison is not completely unbiased due to the different types of simulations and methods for plane-identification used in our study and in others. First, only comparisons between zoom-in simulations that meet the conditions listed in Section 2 of this paper make sense: the dynamical effect of a massive, MW-like disk potential on satellite planes could be an important piece of the puzzle. Second, in the 4GND plot method, no priors are assumed: we do not choose the 11 most massive (or most luminous) satellites among the simulated satellite sample, planes are not required to pass through the center of the host galaxy, or to be thinner than a given Δrms thickness, etc.

One relevant feature of the main plane of satellites observed in the MW is that it presents a relatively high degree of coherent rotation within the plane (see Section 1). This means that the orbital angular momentum vectors (i.e., orbital poles, ${{\boldsymbol{J}}}_{\mathrm{orb}}$) of the constituent satellites are aligned with the normal to the plane. Orbital angular momentum is defined as ${{\boldsymbol{J}}}_{\mathrm{orb}}={\boldsymbol{r}}\times m{\boldsymbol{v}}$, where ${\boldsymbol{r}}$ and ${\boldsymbol{v}}$ are the position and velocity of the center-of-mass of the satellite relative to the center-of-mass of the host disk galaxy.

To study if the satellites included in the high-quality planes detected in the simulations with the extended 4GND plot method are co-orbiting within the plane, we first compute the ${{\boldsymbol{J}}}_{\mathrm{orb}}$ vectors of the satellites at each timestep and project them onto the sphere. Then, we quantify the clustering of ${{\boldsymbol{J}}}_{\mathrm{orb}}$ vectors around a given ${\boldsymbol{n}}$ direction (where ${\boldsymbol{n}}$ is the normal vector to a given plane), and we evaluate the fraction of satellites that are kinematically coherent (i.e., co-orbit) within the plane. To this end, we take this direction ${\boldsymbol{n}}$ as a reference axis and measure the angular distance DA²⁷ to each individual satellite orbital pole. In order to do this systematically at each timestep, we take, as the reference axis, the normal ${\boldsymbol{n}}$ to the best plane (i.e., with lowest c/a) including 50% of the total number of satellites at the respective timestep (see Figure 9).

This is exemplified in Figure 10, for the last timestep analyzed in each simulation. The x-axis shows $1-\cos (\mathrm{DA})$, therefore ranging from 0 to 1, and the y-axis shows the fraction of the total number of satellites with ${{\boldsymbol{J}}}_{\mathrm{orb}}$ enclosed by a certain angular distance DA from the reference axis ${\boldsymbol{n}}$. Since it is only possible to compare to MW data (no proper motion data is available for M31 satellites), the results shown include the observational obscuration bias. Indeed, the dashed line shows the MW case, where we use the latest available data for the confirmed MW satellites (Table 4 in Fritz et al. 2018 calculated from Gaia-DR2 data, or alternatively, for the satellites missing there, Table 4 in Pawlowski & Kroupa 2013).²⁸ For comparison, the dotted line illustrates the expectation from a uniform distribution of orbital poles, and a yellow vertical line is depicted at ${\rm{D}}{\rm{A}}=36\mathop{.}\limits^{{}^{\circ }}87$, which encloses 10%²⁹ of the sphere surface: this is the angle around the VPOS within which it is considered in Fritz et al. (2018) that MW satellites co-orbit (≳40% of MW satellites co-orbit; see Paper I). We see, at these example timesteps, that while for Aq-C${}^{\alpha }$ the fraction of co-orbiting satellites is ∼45%, for PDEVA-5004 it is much lower and below the MW fraction.

Zoom In Zoom Out Reset image size

We use the previous analysis to show in Figure 11 the fraction of co-orbiting satellites in the best-quality planes involving ${f}_{\mathrm{sat}}$ = 50% and 70% of satellites at each timestep. Figure 11 indicates that while in some cases there is consistency with the MW or even a higher degree of co-orbitation (particularly so in the Aq-C${}^{\alpha }$ case), in others, we can see that the fraction of co-orbiting satellites around the ${\boldsymbol{n}}$ direction is very low, despite these directions defining the highest-quality spatial planar arrangements found at that moment. Moreover, the abrupt changes in the fraction of co-orbitating satellites from one timestep to another is a consequence of the different identities of satellites constituting the best-quality planes at close times. Important differences in the fraction of co-orbiting satellites are also found between the 50% versus the 70% case at a same timestep: an indication that while one of the planes shows a kinematical coherence, the other (corresponding to a different density peak) does not.

Zoom In Zoom Out Reset image size

In this respect, it is interesting to compare the times when the C₁ maxima and minima occur (see Figure 6) with a measure of the co-orbitation of the involved satellites. Figure 11 indicates that, in both simulations, vertical green (magenta) bands do not always correspond to a high (low) fraction of co-orbiting satellites.

These results imply that, in general, the best planes in positions found with the 4GND plot method may just be fortuitous spatial alignments of satellites and, therefore, transient structures (see also Gillet et al. 2015; Buck et al. 2016; Ahmed et al. 2017; Fernando et al. 2017; Shao et al. 2019).

In order to efficiently detect kinematically coherent planes of satellites in these simulations, a deeper and more precise analysis of the persistence or not of good-quality planes of satellites across time is needed, which demands using the full six-dimensional phase-space information of satellites for plane searching. In a forthcoming paper (Santos-Santos et al. 2020, in preparation; Paper III), a new method is developed to address the plane of satellites persistence issue.

8.1. Quality Analysis in Terms of N_sat

8.2. Can the Peak Strength Be Used as a Measure of Plane Quality?

The green (magenta) vertical lines in Figure 9 mark the ${T}_{\mathrm{uni}}$ values where the main peak strength C₁ have maxima (minima) in Figure 6. A very relevant result is that maxima occur at ${T}_{\mathrm{uni}}$ values when the c/a value of the ${f}_{\mathrm{sat}}$ = 90% curve is at a minimum, that is, when the quality shows a maximum. And conversely, the magenta lines are close to maxima of c/a, that is, bad qualities. To find out whether this behavior keeps at other ${f}_{\mathrm{sat}}$ values and whether we can use peak strength C_p to measure quality at given ${f}_{\mathrm{sat}}$, we have calculated the main peak strength as a function of ${f}_{\mathrm{sat}}$, i.e., ${C}_{1}({f}_{\mathrm{sat}}),$ ³⁰ and compared it to the c/a of the respective best plane found with same ${f}_{\mathrm{sat}}$ at that timestep.

As an example, in Figure 12, we see for Aq-C${}^{\alpha }$ ("no bias" case, where all satellites are considered) that a correlation exists at all ${f}_{\mathrm{sat}}$. There is, however, an important dispersion particularly at low ${f}_{\mathrm{sat}}$. The same qualitative behavior is found for PDEVA-5004, and in the "bias" case for both simulations. Therefore, we conclude that although indicative, the peak strength C_p is not an accurate enough measure of the quality of its collection of planes.

Zoom In Zoom Out Reset image size

It is interesting to note that, when absolute numbers of satellites (${N}_{\mathrm{sat}}$, without ${N}_{\mathrm{tot}}$-matching) are used instead of fractions (${f}_{\mathrm{sat}}$) to analyze quality, the correlations shown in Figure 12 disappear or are very weak. This is a clear indication that the mutual relationships between quality and peak strength are best manifested when the quality analysis is made in terms of ${f}_{\mathrm{sat}}$, rather than ${N}_{\mathrm{sat}}$.

8.3. Does Radial Compactness of the Satellites Affect Quality?

In this subsection, we analyze the possible correspondence between plane quality and the radial compactness of a satellite system.

Figure 13 shows ${C}_{1,s}$, i.e., the contribution of each satellite s to the main peak, versus the radial distance of the satellite to the center of its host disk galaxy. The left panel corresponds to moments in the evolution of the respective simulations where C₁ is very high (${C}_{1}\sim 33$% for both simulations, see Figure 6), while the right panel corresponds to moments when C₁ is relatively lower (${C}_{1}\sim 10$% for both simulations).

Zoom In Zoom Out Reset image size

There is a clear correlation for satellites that are farther away to contribute more to a given peak (higher ${C}_{p,s}$). At fixed C₁, this effect is slightly intensified in systems where there is an agglomeration of satellites at shorter distances (more compact systems). This is expressed with systematically higher ${C}_{p,s}$ values but a similar slope in the correlation. See, for example, the left panel, where in PDEVA-5004—which is a more radially compact system than Aq-C${}^{\alpha }$—satellites show higher ${C}_{1,s}$ values.

On the other hand, if we compare PDEVA-5004 satellites in the left and right panels (i.e., timesteps where the C₁ is different but there is a similar radial compactness of satellites), we see that the ${C}_{1,s}$ is dramatically lower for the system with lower C₁.

Thus, while a spatial configuration with several central satellites and a few farther away ones (i.e., a compact system) increases the probability of four-galaxy-normals pointing in a similar direction, from these results, it is clear that this geometric effect is not the driving reason for having a high clustering of four-galaxy-normals or high C₁. We conclude that radial compactness does not have a determinant role in setting plane quality.

To address the so-called "Planes of satellites problem" (see discussion and references in Section 1), we have applied, for the first time, the four-galaxy-normal density (4GND) plot method (Pawlowski et al. 2013) to hydrodynamical simulations. An extension of the method, sketched in Santos-Santos et al. (2019, Paper I), is presented in detail and discussed in this paper.

Our choice has been zoom-in simulations because the higher resolution (as compared to large-volume simulations) one can reach by using this method allows MW-type galaxies with massive, extended disks to form, so that the dynamical effects of a live disk potential on satellite planes (torques, destruction of satellites on radial orbits) can be accounted for. Another advantage of zoom-in simulations is that they offer the possibility that these disk galaxies are surrounded by a number of resolved satellites high enough to be comparable to current samples of confirmed MW or M31 satellites (∼30, see the McConnachie 2012 database). Well-behaved disks and a large enough satellite number are conditions not recovered, at the moment, by larger-volume cosmological simulations.

The extension to the 4GND method is designed to identify, systematically catalog, and study in detail the quality of the predominant spatial planar configurations of satellites revealed by over-densities in the 4GND plots. It allows us to extract high-quality planes out of the number of combinations we can form with ${N}_{\mathrm{sat}}$ satellites from a sample of size ${N}_{\mathrm{tot}}$, with a low computational cost. Quality is evaluated through the outputs of a Tensor of Inertia analysis (ToI; Metz et al. 2007) using a normal-regression fitting technique. Quantitatively, planes (i.e., best-fitting solutions with high medium-to-long axes ratio, $b/a\sim 1$), have a good quality if they are populated relative to the sample size ${N}_{\mathrm{tot}}$ (high ${f}_{\mathrm{sat}}\equiv {N}_{\mathrm{sat}}/{N}_{\mathrm{tot}}$) and thin (low short-to-long axis ratio c/a). Being a two-parameter notion, the quality of two or more planes can be compared if one has lower c/a and higher ${f}_{\mathrm{sat}}$ than another, or if either ${f}_{\mathrm{sat}}$ or c/a are constant.

Density peaks are determined by local, isolated maxima in the 4GND plot. We have defined the peak strength, C_p, as the % of weighted four-galaxy-normals within 15° of the peak center. Peaks can be compared with each other across times through their respective strengths and with those of the MW and M31 satellite systems as well. Different satellites contribute differently to C_p. The satellite s contribution to peak p (i.e., ${C}_{p,s}$) is defined as the normalized weighted count of s contributions to four-galaxy-normals placed within 15° of the peak center. Satellites are ordered by decreasing C_p,s to a peak, and a plane is fitted to groups of increasing N_sat satellites following this order. This yields, for each density peak, a collection of planes.

In Paper I, we report on the application of the extended 4NDP method to the confirmed MW and M31 satellites. The method extension reveals a richer planar structure (i.e., a collection of planes rather than a unique plane per density peak), therefore allowing us to find planes of satellites around the MW and M31 with higher qualities than those previously reported with a given ${N}_{\mathrm{sat}}$. We find a second populated, high-quality plane around M31. Another important result, in view of the narrow range of satellite mass distributions that can be currently afforded in hydrodynamical simulations, is that satellite mass plays no role in determining a satellite's membership or lack thereof to the respective best-quality planes. This enables us to perform, through the extended 4GND plot method, comparisons between results from simulations and MW/M31 data (where the satellite mass range is wider).

In this paper, we present results of a detailed search for positional planar structures of satellites in two different (initial conditions, sub-grid modeling and numerical approaches) zoom-in cosmological hydrodynamical simulations of isolated disk galaxies that resemble the MW: Aq-C${}^{\alpha }$ and PDEVA-5004. They meet the conditions of having a central host galaxy with an extended disk, an overall quiet merger history after virialization, a numerous (∼30) satellite population, and more than 50 baryonic particles per satellite. No other particular selection criteria have been applied. In particular, we focus our analysis on the best-quality planes with a given satellite fraction ${f}_{\mathrm{sat}}$. We compare to the best-quality planes found in the MW and M31 systems at z = 0 with given ${f}_{\mathrm{sat}}$, i.e., those found from their so-called respective #1 Peaks as defined in Paper I.

The analysis goes from the halo virialization time to low redshifts, along the slow phase of host mass assembly. The halo mass-growth histories of Aq-C${}^{\alpha }$ and PDEVA-5004 present, per usual, two phases and no other relevant particularities. Satellite samples have been identified at z = 0.5. The number of satellites ${N}_{\mathrm{tot}}$ considered in the analysis at each timestep varies. In particular, satellites are not considered when they are orbiting beyond 450 kpc from the center of the host galaxy, after they have been accreted by the host galaxy, or if they are within the avoidance volume when we apply an observational obscuration bias to compare with the MW. The distributions of satellite radial distances to the center-of-mass of the main galaxy vary with time. Expressed in terms of ${f}_{\mathrm{sat}}$, Aq-C${}^{\alpha }$'s (PDEVA-5004's) radial distance distribution is very close to that of M31 (the MW), at particular times.

The analysis of density peaks reveals that C₁ varies with ${T}_{\mathrm{uni}}$, with time intervals where it is comparable or even higher than the MW value, ${C}_{1,\mathrm{MW}}=22.91 \% \pm 0.26 \%$. The number of strong peaks changes with time too, reaching values of 1–3 at most and only in those time intervals when C₁ is high. In this regard, density peaks found in simulations are similar to those observed in strength and number during given time intervals.

We find planar (i.e., high b/a and low c/a) configurations of satellites at all studied timesteps in both simulations. Indeed, no filamentary (i.e., $b/a\ll 1$) configurations have been detected in the period analyzed. The extended 4GND plot method allows us to identify, at many timesteps in both simulations, (i) planes of satellites with qualities that are compatible with the observed ones at z = 0 including a specific ${f}_{\mathrm{sat}}$, and also in some cases, (ii) planar structures that are compatible with the observed ones for all ${f}_{\mathrm{sat}}$.

We study the best-quality planes (i.e., with lowest c/a) including a fixed ${f}_{\mathrm{sat}}$ found at each timestep. In both simulations, their c/a values change with time, independently of the ${f}_{\mathrm{sat}}$ considered. Planes compatible with the observed ones in the MW and M31 at z = 0 are found at different timesteps or time intervals. Interestingly, these timesteps turn out to coincide with the time intervals where C₁ shows maxima. More specifically, a correlation has been found between C₁ and c/a at fixed ${f}_{\mathrm{sat}}$ but with important dispersion. Therefore, C₁ can be used as an estimation of plane quality but not to accurately measure it. Another interesting result is that the highest-quality planes are often close to perpendicular to the host disk plane. Ideally, perpendicular planes of satellites would not suffer torques from the disk.

No new information on quality is provided when using Δrms as a quality indicator once the satellite system size has settled. And no new information on quality (in terms of comparison to observational satellite planes) is obtained either when using absolute satellite numbers ${N}_{\mathrm{sat}}$ instead of fractions ${f}_{\mathrm{sat}}$ = $\tfrac{{N}_{\mathrm{sat}}}{{N}_{\mathrm{tot}}}$, with ${N}_{\mathrm{tot},\mathrm{obs}}$ matched to ${N}_{\mathrm{tot},\mathrm{sim}}$.

Interestingly, when the observational obscuration bias is applied, slightly higher peak strengths are measured, as well as somewhat lower c/a and Δrms values than in the non-biased case. The is due to solid-angle restrictions in the biased case, where satellites orbiting at low latitudes are neglected.

No clear, conclusive signal on the correspondence between plane quality and the radial compactness of a satellite system has been detected in this work. While satellites at larger distances from the host galaxy provide somewhat higher ${C}_{1,s}$ than nearby ones, the spatial satellite configurations that show the highest C₁ (and therefore highest overall plane qualities) are not those most compact.

We have further investigated if satellites composing the high-quality planes of satellites found with the extended 4GND plot positional method present a common orbitation within the plane they describe. We find that, in some cases, the fraction of co-orbiting satellites is very low, which we interpret as a sign of these positional planes consisting partly of interlopers: satellites that fall within the plane accidentally. Therefore, planes found with this method based on positional data in general do not constitute a kinematic unit and, in some cases, could be non-persistent in time. Because of this, the search for the physical reasons favoring or destroying positionally detected high-quality planes cannot be meaningfully addressed. In Santos-Santos et al. (2020, in preparation; Paper III), a new methodology is introduced where the full six-dimensional phase-space information of satellites is used, leading to the determination of persistent, kinematically coherent planes of satellites in both simulations.

Summing up, the application of the 4GND plot method (Pawlowski et al. 2013) with its extension presented here, (see also Santos-Santos et al. 2019, Paper I) to two zoom-in hydrodynamical simulations of disk galaxies, leads us to the following conclusions:

• 1.  
  Satellites are organized in planar configurations (not filamentary, i.e., $b/a\sim 1$) in both the Aq-C${}^{\alpha }$ and PDEVA-5004 simulations, at all timesteps analyzed. The plane short-to-long axis ratio c/a, and the plane population (${f}_{\mathrm{sat}}$) measure plane quality.
• 2.  
  The strengths of the strongest peaks in the 4GND plots (C₁) vary with time. Their values are, at all times, higher than that of the strongest peak in M31, and along given periods, also higher than that of the MW.
• 3.  
  The c/a ratios of the thinnest planes found including a fixed ${f}_{\mathrm{sat}}$ vary with cosmic time, and during some periods, reach a high quality. Along these good-quality time intervals, they are compatible with the observed planar structures found in the MW and in M31 at z = 0. These periods coincide with those when C₁ reaches maximum values. The timescale for these plane quality changes is $\sim 0.5\mbox{--}1\,\mathrm{Gyr}$.
• 4.  
  c/a and C₁ show correlations with increasing dispersion as ${f}_{\mathrm{sat}}$ decreases. C₁ can be used as a quality indicator, but not to accurately measure it.
• 5.  
  The application of the observational obscuration bias enhances plane quality, either measured by peak strength or plane thickness.
• 6.  
  The orientations of planes of satellites with respect to the disk of their host galaxy change with time. In most cases, planes are close to perpendicular to the disk during periods of good quality.
• 7.  
  The compactness of the distribution of satellite-host radial distances does not have a driving role at setting the quality of plane of satellites.
• 8.  
  In agreement with previous findings, the fraction of co-orbiting satellites found in high-quality positionally detected planes is rather low, suggesting that these planes do not represent a kinematic unit and are not persistent in time.
• 9.  
  The plane persistence issue in observations and simulations cannot be properly addressed unless the full six-dimensional phase-space information is considered. Such a methodology will be developed in a forthcoming paper.

The general conclusion of this paper is that even if two galaxies do not make a statistical sample, the fact that these two so different MW-type disk galaxies (whose selection method is quite general) do have, at given time intervals, high-quality positional satellite planes, would suggest that these planes can be expected not to be infrequent in ΛCDM L* disk galaxies in periods when they are free of major merger events in their assembly history.

We thank the anonymous referee for useful comments and suggestions that have helped improve this work. This work was supported through MINECO/FEDER (Spain) AYA2012-31101, AYA2015-63810-P, and MICINN/FEDER (Spain) PGC2018-094975-C21 grants. I.S.S. acknowledges support by the Arthur B. McDonald Canadian Astroparticle Physics Research Institute. This work used the Ragnar cluster funded by Fondecyt 1150334 and Universidad Andrés Bello, and Geryon cluster (Pontificia Universidad de Chile). We used a version of Aq-C-5 that is part of the CIELO Project run in Marenostrum (Barcelona Supercomputer Centre). This project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement No. 734374 (LACEGAL-RISE). I.S.S acknowledges funding from the same grant for a secondment at the Astrophysics group of Univ. Andrés Bello (Santiago, Chile), and from the Univ. Autónoma de Madrid for a stay at the Leibniz Institut fur Astrophysik Potsdam (Germany). I.S.S. thanks Dr. Patricia Tissera and Dr. Noam Libeskind for kindly hosting her.