The Origin of Kinematically Persistent Planes of Satellites as Driven by the Early Evolution of the Cosmic Web in ΛCDM

The initial conditions of this simulation come from the Aquarius Project (Springel et al. 2008), a selection of DMO MW-sized halos formed in a ΛCDM simulation run in a 100 h⁻¹ Mpc side cosmological box. A new resimulation of the so-called Aquarius-C halo (hereafter Aq-C^α ), including the hydrodynamic and subgrid models described in Pedrosa & Tissera (2015; see also Scannapieco et al. 2005, 2006), has been analyzed in this work. The initial mass resolution of the baryonic and dark matter particles, m_bar and m_dm, respectively, and the parameters of the cosmological model are given in Table 1 (C and S blocks).

Table 1. Some Parameters of the Cosmological Model (C Block), Simulation (S Block), Host Galaxy (HG Block), and LVs (LV Block) for both the Aq-C^α and the PDEVA-5004 Systems

Block	Parameter	Aq-Cα	PDEVA-5004
C	Ωm	0.25	0.28
	Ωb	0.04	0.04
	ΩΛ	0.75	0.72
	H0 [km s–1 Mpc–1]	73.0	70.0
S	ns	1.0	1.0
	σ8	0.9	0.81
	mbar [M⊙]	4.1 × 105	3.94 × 105
	mdm [M⊙]	2.2 × 106	1.98 × 106
		z = 0	z = 0
HG	rvir,z=0 [kpc]	241.26	181.43
	Mvir [M⊙]	1.82 × 1012	3.44 × 1011
	Mstar [M⊙]	8.6 × 1010	3.05 × 1010
	Mgas [M⊙]	7.4 × 1010	8.6 × 109
	Temp${}_{\mathrm{lim}}$ [K]	2.0 × 105	2.0 × 104
	Tvir [Gyr]	7	6
	Tta,halo [Gyr]	4	3
LV	zhigh	8.45	10.00
K = 10	RLV [kpc]	255.26	166.07
	MLV [M⊙]	3.15 × 1012	1.04 × 1012
K = 15	RLV [kpc]	382.89	249.10
	MLV [M⊙]	6.59 × 1012	3.17 × 1012
K = 20	RLV [kpc]	510.52	332.14
	MLV [M⊙]	8.11 × 1012	6.31 × 1012

Note. Gas particles whose temperature is higher than Temp${}_{\mathrm{lim}}$ [K] are taken as pressurized hot gas particles and have not been considered here to calculate LV deformations. See text for parameter definitions.

Download table as:  ASCIITypeset image

The halo turnaround and virialization happen at a Universe age of T_ta,AqC ≃ 4 Gyr and T_vir,AqC ≃ 7.0 Gyr, respectively. In this case, 25% of the halo mass is accreted after collapse. Properties of this host galaxy measured at the final redshift of z = 0 are given in Table 1 (host galaxy or HG block).

This system presents a quiet history from z ≈ 1.5 to z = 0, where no major mergers occur. By an age of the Universe of T_uni ∼ 11.5 Gyr (z = 0.15), the main galaxy captures a massive dwarf (M_bar = 5.02 × 10⁹ M_⊙) carrying its own satellite system (six members with the satellite identification criteria used in this paper). The capture has been analyzed in Paper III, where it was shown that it has no perturbing effects on the dynamical behavior of the rest of Aq-C^α 's satellite population. As a low-redshift event, this capture is beyond the scope of this paper.

The PDEVA-5004 system comes from a zoom-in resimulation made with the PDEVA code (Martínez-Serrano et al. 2008) of a halo identified in a ΛCDM run. Parameters characterizing the cosmological model, the run, the host galaxy at z = 0, and some characteristic timescales are given in Table 1. For more details, see Doménech-Moral et al. (2012) and Paper II and references therein.

The host halo turnaround and virialization events occur at Universe age T_ta,5004 ≃ 3 and T_vir,5004 ≃ 6 Gyr, respectively, a bit earlier than in the Aq-C^α system. Only 20% of the virial mass is assembled after T_vir,5004, and no major mergers show up in the slow phase.

The identification of satellites has been done, in both simulations, at two different times, i.e., z = 0 and z = 0.5 (T_uni ∼ 8.6 Gyr), in order to include satellites that may end up accreted ¹⁴ by the central disk galaxy and do not survive until z = 0. We define satellite galaxies as those objects with stars (M_* > 0) that are bound to the host galaxy within any radial distance. To ensure objects are bound, we have computed their orbits back in time. The friends-of-friends algorithm and the SubFind halo finder (Springel et al. 2001) have been used to set structures and substructures in the Aq-C^α system (see also Dolag et al. 2009), while PDEVA-5004 satellites were selected using IRHYS. ¹⁵ By tracing the particle IDs across snapshots, we have followed individual satellites back in time to times even before they were assembled as bound structures. The total number of satellites is 34 (35) in Aq-C^α (PDEVA-5004). Of these, 32 (26) survive until z = 0. Satellites will be addressed throughout the paper with an identification code (see, for example, in Figures 1 and 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Satellites in Aq-C^α and PDEVA-5004 span baryonic mass ranges of M_bar = 8.5 × 10⁶−9.9 × 10⁸ and M_bar = 3.9 ×10⁷–1.8 × 10⁸ M_⊙, respectively (we recall that we impose that satellites are resolved with a minimum of 50 baryonic particles). The satellite mass distributions are addressed in Papers II, and III, where it is shown that, at least in these two simulations and within the satellite mass range available here, the baryonic mass is not a satellite property that determines whether or not a satellite belongs to the corresponding positional or persistent plane. This is an important result given our limited mass range due to the current computational possibilities.

Satellites present a diversity of orbital histories in both simulations. While a small fraction of satellites end up accreted by the central galaxy's disk (two in the Aq-C^α system, nine in PDEVA-5004), most of them have regular orbits with stable apocentric and pericentric distances. Some of them are backsplash satellites. Finally, a few satellite cases have just been captured by the main galaxy's halo and have not yet had time to complete their first pericentric passage. Such late satellite incorporations have only been found in the Aq-C^α system, where several satellites show first pericenters later than T_uni = 10 Gyr. Conversely, all of PDEVA-5004's satellites are fully incorporated to the system by that time.

Figure 1 shows the evolution of the orbital angular momentum of satellites, ${\vec{J}}_{\mathrm{orb}}$. Specifically, the top panels show the specific ${\vec{J}}_{\mathrm{orb}}$ vector moduli or magnitudes (i.e., normalized by their baryonic mass), hereafter sJ_orb, for the Aq-C^α system; see the figure caption. We can see that, for T_uni lower than ∼4 Gyr, most satellites' sJ_orb increases. As mentioned in Section 1, the current paradigm sets that galactic systems acquire their angular momentum $\vec{J}$ at high redshift, prior to the system turnaround. After this moment, the lever arm becomes too low for an effective angular momentum gain, and sJ_orb stops increasing or even decreases. Most satellites in Figure 1 show this behavior at high redshift. The same is true for the PDEVA-5004 satellite system. Thus, the results of our simulations fit into the TTT scheme as far as the qualitative high-z performance of some satellite's sJ_orb is concerned. In the slow phase of mass assembly, sJ_orb is roughly conserved in many cases, while in others it is not.

The behavior of orbital poles after T_vir was studied in Paper III, (see the Aitoff projections in their Figures 2 and B1). It was shown that the orbital angular momentum directions of most KPP satellites present only modest changes with time after T_vir, while for non-KPP members, changes can be more relevant.

For the purposes of this paper, it is very relevant to analyze angular momentum conservation before T_vir. To have a deeper understanding of when pole conservation sets in, in the middle and bottom panels of Figure 1, we plot the orbital pole (polar and azimuthal angles with respect to axes that are kept fixed in time) evolution of KPP and non-KPP satellites within the T_ta,AqC ≤ T_uni ≤ 10 Gyr interval, i.e., an interval beginning at halo turnaround time and centered at T_vir. Interestingly, some KPP member satellites maintain their orbital pole directions before T_vir roughly constant, while others change them only smoothly. Orbital pole changes are more frequent in non-KPP member satellites. We can also see that orbital pole clustering sets in already in the fast phase of mass assembly, an issue to be discussed in more detail in the next sections.

As explained in Section 1, a kinematically coherent persistent plane (KPP) is a fixed subset of satellites whose orbital poles are conserved and remain clustered for a long period of cosmic time, defining a high-quality planar structure in positional space in terms of a tensor of inertia (TOI) analysis (Cramér 1999).

In Paper III, axes of maximum satellite co-orbitation, ${\vec{J}}_{\mathrm{stack}}$, were determined through the so-called "Scanning of Stacked Orbital Poles Method." ¹⁶ One (two) such axes have been found in the Aq-C^α (PDEVA-5004) systems, giving 13 satellites for the Aq-C^α system KPP and 10 and 7 satellites, respectively, for the two PDEVA-5004 KPPs.

As an illustration of these results, in Figure 2, we show the time evolution of $\alpha ({\vec{J}}_{\mathrm{orb},{\rm{i}}},{\vec{J}}_{\mathrm{stack}})$, the angle formed by the axes of maximum co-orbitation, and the orbital poles of each satellite member of the KPP planes. This figure helps us answer the following important question: when is orbital pole clustering established? It is worth mentioning that clustering can appear at times much earlier than T_vir, before all the satellite members of a KPP are bound to the main galaxy.

According to Fritz's co-orbitation criterion mentioned in Section 1, translated into clustering properties, the ith satellite orbital pole will be considered to be aligned to the $\vec{V}$ vector when the $\alpha ({\vec{J}}_{\mathrm{orb},{\rm{i}}},\vec{V})$ angle is smaller than α_co-orbit = 3687. Figure 2 indicates that a fraction of satellites are aligned to their corresponding ${\vec{J}}_{\mathrm{stack}}$ already at T_uni = 2 Gyr. For those that are not aligned that early, important changes in their pole direction occur in the ∼2–4 Gyr time interval for both simulations. Specifically, an alignment timescale, T_cluster,i (relative to a given $\vec{V}$ vector), can be defined for the ith satellite as the time when the $\alpha ({\vec{J}}_{\mathrm{orb},{\rm{i}}},\vec{V})$ angle reaches a cosine higher than 0.8. Correspondingly, a measure of the time when a given KPP sets in as a clustered bundle of orbital poles around $\vec{V}$, T_cluster,V , is given by the median of the T_cluster,i values for each of its satellite members. An application of this protocol to data shown in Figure 2 gives results for the median timescales (and their corresponding differences with the 25th and 75th percentiles) of pole alignments with the ${\vec{J}}_{\mathrm{stack}}$ axes reported in Table 2, T_{cluster,Jstack} entry.

Table 2. Summary of Timescales Highlighted throughout the Paper (Values Correspond to Ages of the Universe, T_uni in Gyr)

Timescale	Aq-Cα		PDEVA-5004
T ${}_{\mathrm{dir},\,{{\rm{e}}}_{3}}$ (5)	∼2.0		∼0.5
T ${}_{\mathrm{shape},\,{{\rm{e}}}_{3}}$ (6)	∼4.5		∼3.5
	KPP	nKPP	KPP1	KPP2	nKPP
Tcluster, Jstack (2)	${3.7}_{-1.5}^{+0.9}$	⋯	${2.2}_{-0.2}^{+0.2}$	${3.5}_{-1.1}^{+1.0}$	⋯
T ${}_{\mathrm{align},\,{{\rm{e}}}_{i}}$ (7)	${4.5}_{-1.2}^{+2.5}$	⋯	${2}_{-0.2}^{+0.2}$	${3.5}_{-0.8}^{+0.5}$	...
Tdist, plane (10)	∼3.5	∼3.5	∼4.5	∼3.0	∼3.0
Tsat, infall	${7.2}_{-3.9}^{+1.8}$	${8.8}_{-4.5}^{+1.1}$	${6.4}_{-1.2}^{+1.4}$	${4.5}_{-0.3}^{+0.7}$	${4.3}_{-0.7}^{+0.5}$
Tsat, apo1	${1.7}_{-0.2}^{+1.3}$	${3.8}_{-2.1}^{+0.9}$	${2.5}_{-0.3}^{+0.5}$	${2.4}_{-0.2}^{+0.3}$	${2.2}_{-0.2}^{+0.2}$

Note. If a timescale is indicated in a specific figure, the figure number is indicated in parentheses. Timescales in the first block apply to the LVs of each simulation. Times in the second block concern the satellite populations and are given by medians and 25th–75th percentiles when possible. For definitions, see text and timescale summary in Section 8.3.

Download table as:  ASCIITypeset image

We see that clustering sets in very early, especially in the case of the PDEVA-5004 system. To understand why clustering sets in that early, we need to unveil its causes by analyzing LV deformations resulting from mass flows and their consequences.

The kinematic morphological parameter κ_rot is first defined as the fraction of the kinetic energy of a given galaxy coming from ordered motions (i.e., in-plane and close to circular) relative to its disk axis; see Sales et al. (2012). Quantitatively, this definition is

$$\begin{eqnarray}&&{\kappa }_{\mathrm{rot}}=\displaystyle \sum _{i}{w}_{m,i}{({v}_{i,\phi }/{v}_{i})}^{2},\end{eqnarray} \tag{ 1 }$$

where w_m,i = m_i /M_s is the mass weight of the ith constituent element (with m_i its mass and M_s the mass of the system), and v_i,ϕ and v_i are the tangential velocity relative to the system center of velocity in the plane normal to the fixed axis and the modulus of the velocity of the ith constituent element, respectively. κ_rot is a kinematic indicator of the morphology of a galaxy, such that those with κ_rot lower than 0.5 are considered to be a spheroid (dispersion-dominated galaxy), while galaxies with higher values are rotation-dominated.

This definition can be easily extended to the set of constituent elements of a given object or system, here satellite sets, relative to a fixed axis, here the ${\vec{J}}_{\mathrm{stack}}$ axes. Specifically, we aim at studying the evolution of the κ_rot parameter of KPPs as satellite sets versus that of sets of satellites outside these structures. In this kinematic analysis, the mass of each satellite is irrelevant, so we treat them as point sources, thus dropping the weighting term in Equation (1).

Results for KPP planes are given in Figure 3, where we show the median values (with their 25th–75th percentiles) for the ${({v}_{i,\phi }/{v}_{i})}^{2}(t)$ functions (note that the weighted mean is κ_rot). A first remarkable result is the behavior of κ_rot as a function of time, with largely constant median values (fluctuations in time are a result of low number statistics combined with parameter peaks at apocenter and pericenter). We see that satellites in KPP planes show high κ_rot values, and thus they behave as morphological disks from a kinematic perspective as well. Conversely, except for a few peaks, these medians are at any time lower than 0.5 for satellites outside KPPs. Finally, we see that high κ_rot parameter values set in early for satellites that will be KPP members, after they increase at high redshift. A rough estimation of a timescale for the establishment of the maximum ordered disky motion, ${T}_{\max ,\,\mathrm{rot}}$, allowed by Figure 3, points toward T ${}_{\max ,\,\mathrm{rot}}\sim$ 5 and 4 Gyr for the Aq-C^α KPP and PDEVA-5004 KPP2 systems, respectively. In the case of the PDEVA-5004 KPP1 system, satellites have already shown high κ_rot values since very early times, such that ${T}_{\max ,\mathrm{rot}}$ is virtually <2 Gyr. Indeed, as seen in Figure 2, KPP1 satellites have had their orbital poles aligned since very high redshifts, defining an in-plane motion.

Zoom In Zoom Out Reset image size

Together with the results above on the timescales for clustering establishment, this result reinforces the idea that satellite kinematic coherence (that is, aligned orbital poles) and the appearence of disklike orbits with high κ_rot values (i.e., in-plane and circular motion) set in at very early times, when the protosatellites are but part of the galaxy-to-be evolving environment, moving within it. Some answers to this possibility will be given in the next sections.

As said above, taking the LV as a whole, the ${I}_{{ij}}^{{\rm{r}}}$ eigenvectors, ${\hat{e}}_{1}(z)$, ${\hat{e}}_{2}(z)$, and ${\hat{e}}_{3}(z)$, mark the directions of the major, intermediate, and minor axes of its inertia ellipsoid at redshift z. It is very important to quantify the changes in such directions as cosmic evolution proceeds. It is particularly important to find out whether or not the three eigendirections become fixed at a given time, say, ${T}_{\mathrm{dir},{{\rm{e}}}_{i}}$. Should this happen, mass rearrangements at LV scales after ${T}_{\mathrm{dir},{{\rm{e}}}_{i}}$ would be organized in terms of a "skeleton" or fixed preferred directions, with the ${\hat{e}}_{3}$ direction corresponding to that of maximum compression in the LV deformation, or the direction along which the overall mass flow has been maximum.

In Figure 5, we show the time evolution of A_i (z), the angle formed by the eigenvectors ${\hat{e}}_{i}(z)$ and ${\hat{e}}_{i}(z=0)$, with i = 1, 2, 3, for the Aq-C^α simulation. That is, we measure the deviations from the eigendirections at a given z with respect to the final eigenvectors for an LV scale of K = 15. Notice that only two out of the three A_i angles are independent in such a way that if, for instance, A₁ = 0, then A₂ = A₃.

Zoom In Zoom Out Reset image size

We see that, when K = 15, A_i (z) change at high redshift. A₃(z) smoothly vanishes by T_uni ∼ 2 Gyr ≡${T}_{\mathrm{dir},{{\rm{e}}}_{3}}$. Results for PDEVA-5004 with K = 15 are shown in Table 2. Changes in ${\hat{e}}_{3}$ and the other principal directions become unimportant very early, except for a small change of 5% at most in ${\hat{e}}_{1}$ and ${\hat{e}}_{2}$ occurring around T_uni ∼ 6 Gyr. That is, the LV deformations get their three eigendirections fixed at high redshift, defining the freezing-out timescale T_freeze = 4 (2) Gyr for the Aq-C^α (PDEVA-5004) simulations well before the systems enter the slow phase of mass assembly.

It is important to figure out whether or not this behavior depends on the LV scale. In the Appendix (Figure A1), we show that the evolution of the principal directions for an LV scale of K = 20 for the Aq-C^α simulation makes it essentially invariant under this change in scale. The same results are obtained for the PDEVA-5004 simulation. That being said, we proceed with our analysis using a fiducial K = 15 value.

The shape evolution of the LVs is presented in Figure 6. We show the principal axes as a function of Universe age (results for Aq-C^α in the top panel and PDEVA-5004 in the middle panel) and their b/a and c/a ratios, color-coded by the Universe age (bottom panel; see sidebar).

Zoom In Zoom Out Reset image size

A remarkable result is the continuity of the a(t), b(t), and c(t) functions for all the LVs, with no mutual exchange of their respective eigendirections across evolution; i.e., the local skeleton is continuously built up, consistent with Hidding et al. (2014) and Robles et al. (2015). We see that at high redshift, the principal axes have very similar lengths, as expected for a sphere. Then, for both LVs in these plots, the a principal axis monotonously grows, while c decreases very rapidly and then the change slows down. Also, some periods when the change is very slow occur in PDEVA-5004. The b(t) axis, corresponding to the ${\hat{e}}_{2}$ principal direction, shows a decreasing behavior from T_uni ∼ 4 Gyr onward in the Aq-C^α system, remaining almost constant until that moment. b(t) shows some periods where it is almost constant in the PDEVA-5004 simulations as well.

The overall shape deformation of these LVs is well quantified through the evolution of the c/a and b/a ratios and the triaxiality parameter, T. The bottom panel of Figure 6 shows the c/a versus b/a diagram, where different shape specifications mark their corresponding parameter spaces, and the T = 1.0, 0.7, and 0.3 isotriaxiality curves have been drawn. We see that in both simulations, after a rapid evolution from a spherical shape (b/a ≃ c/a ≃ 1) toward more flattened structures (c/a always decreases, while b/a is constant along some periods), shape changes slow down as they increase their prolateness, with the Aq-C^α system reaching a final shape more prolate than the PDEVA-5004 system. Indeed, Figure 4 explicitly shows how a planar structure in Aq-C^α is clearly formed by z ∼ 1.0. In the Aq-C^α simulation, we note a fast change in b(t) by T_uni = 4 Gyr, marking a shape deformation from oblate to triaxial (see the "knee" feature in the bottom panel).

To better quantify how quickly the a(t), b(t), and c(t) functions change, their time derivatives are also plotted in the upper and middle panels of Figure 6. We see fast changes in the c(t) principal axis up to T_uni ∼ 4.5 Gyr (∼3.5 Gyr) for Aq-C^α (PDEVA-5004), marked by arrows in the top and middle panels as ${T}_{\mathrm{shape},{{\rm{e}}}_{3}}$ (see also Table 2), and then the decrement rate becomes almost constant and very low for PDEVA-5004 between T_uni ∼ 5 and 8 Gyr. A rapid decrement in c(t) reflects the strength of the early mass inflows in the ${\hat{e}}_{3}$ direction. It is worth mentioning that any anisotropic mass inflow implies a mass rearrangement and, consequently, a change in the principal axis values. Figure 6 informs us when these anisotropic mass flows become unimportant.

An important point is the possible scale dependence of results on LV shape evolution. Our analyses indicate that no remarkable, qualitative changes show up for either simulation at the larger scales tested here in the evolution of their respective principal axis lengths or the derivatives of the a(t), b(t), and c(t) functions.

An illustration of the LV global shape evolution for the Aq-C^α simulation just described is provided by Figure 4. This figure shows that by redshift z ∼ 1, a flattened structure, normal to the ${\hat{e}}_{3}$ principal direction (hereafter the ${\hat{e}}_{3}$-structure or, in this case, the ${\hat{e}}_{3}$-wall), clearly stands out for the first time in this plot. As explained in Section 6.2, from this time onward, vertical flows of LV material onto this flattened structure weaken to a great extent, while mass motions within it still occur, as can be seen in the projections on the X–Y plane of the LV evolution. These mass motions lead, by z ∼ 0.5, to the appearance of a prominent filament parallel to the ${\hat{e}}_{1}$ principal direction, where most mass now piles up at the expense of the wall-like structure population.

We note that the ${\hat{e}}_{3}$-structure is not a simple sheet, as the caustics appearing in the ZA, but the result of a complex history of mass (smaller-scale CW elements) incorporations. Similarly, the filament in the ${\hat{e}}_{3}$-structure is not a simple filament. In turn, this predominant filament also disappears later on in favor of halos (see z = 0 panels of Figure 4) that form a prolate configuration; see Figure 6, bottom panel.

The important point here is that the LV shape evolution expresses the CW unfolding at its scale, with all its complexity and diversity. The same is true in the case of the PDEVA-5004 system, with the difference that in this case the LV global shape is triaxial along almost all the evolution, leading to a prolate ${\hat{e}}_{3}$-structure (see bottom panel in Figure 6).

We first consider how the orbital poles of individual satellites are oriented relative to the LV principal directions, ${\hat{e}}_{1}$, ${\hat{e}}_{2}$, and ${\hat{e}}_{3}$.

The time development of the alignments of orbital poles with the principal axes is given in Figure 7 for both the Aq-C^α (upper block of panels) and PDEVA-5004 (bottom block of panels) systems. In this figure, panels in the first, second, and third columns stand for the angles formed by the satellite orbital poles and the principal directions ${\hat{e}}_{1}$, ${\hat{e}}_{2}$, and ${\hat{e}}_{3}$ of the LV-reduced TOI. Thin lines correspond to individual satellite pole orientations, colored according to the satellite identity as given in the sidebars. Cyan, orange, and purple thick lines represent the median values of the angle set at each time step for ${\hat{e}}_{1}$, ${\hat{e}}_{2}$, and ${\hat{e}}_{3}$ alignments, respectively, while the shaded areas correspond to the respective 25th and 75th percentiles.

Zoom In Zoom Out Reset image size

According to Figure 7, satellites in the Aq-C^α KPP orbit on planes close to normal to the direction of maximum global compression of the LV deformation, ${\hat{e}}_{3}$. That is, they move approximately within the flattened structure the initially spheroidal LV is deformed into along evolution. Two satellites are already aligned by T_uni ∼ 2. For those that are not, changes in their pole directions leading to improved alignments mostly occur between T_uni ∼ 2 and 4.5 Gyr.

A clear alignment signal (i.e., small angles) stands out for the PDEVA-5004 system, where the KPP1 satellite orbital poles tend to be close to parallel to the ${\hat{e}}_{1}$ axis. Indeed, the median of $\cos ({\vec{J}}_{\mathrm{orb}},{\hat{e}}_{1})$ is ≃0.9 after T_uni ≃ 2 Gyr. Therefore, after that, KPP1 satellites tend to orbit on planes close to normal to the ${\hat{e}}_{1}$ direction. On the other hand, satellite members of the KPP2 plane tend to have their poles aligned with the ${\hat{e}}_{3}$ directions. Thus, KPP2 satellites tend to orbit on planes close to normal to the direction of maximum global compression for the LV under consideration. It is worth mentioning that the alignment improves (the shaded area becomes narrower) at lower redshifts for two out of three KPPs we have identified.

Most satellites outside any kinematically coherent plane show no alignment with either the ${\hat{e}}_{1}$, ${\hat{e}}_{2}$, or ${\hat{e}}_{3}$ directions. They also show a larger spread in their angle values than those in the KPPs.

To extract information on the overall timescale for orbital pole alignment with the principal directions, we look for the Universe age when the median values of the $\alpha ({\vec{J}}_{\mathrm{orb}},{\hat{e}}_{i})$ angles (Figure 7), with i = 3 for the Aq-C^α and PDEVA-5004 KPP2 systems and i = 1 for the PDEVA-5004 KPP1 system, reach a value below α_co-orb = 3687. Results are given in Table 2, entry ${T}_{\mathrm{align},{{\rm{e}}}_{i}}$. We note the high dispersion of this value for the Aq-C^α KPP. T ${}_{\mathrm{align},{{\rm{e}}}_{i}}$ also gives a measure of the timescale for the clustering of orbital poles. Results for ${T}_{\mathrm{align},{{\rm{e}}}_{i}}$ are consistent with the corresponding T_{cluster,Jstack}, i.e., the timescale for the setting in of orbital pole clustering measured through alignments with ${\vec{J}}_{\mathrm{stack}}$ (see Table 2).

Another important point concerning satellite orbital pole alignments is how the probabilities of alignment with the axes of maximum co-orbitation, on the one hand, and the principal directions, on the other hand, are related to each other. For the PDEVA-5004 simulation, the respective numbers of orbital poles aligned with the ${\hat{e}}_{1}$ and ${\hat{e}}_{3}$ principal directions are eight (all of them aligned with the ${\vec{J}}_{\mathrm{stack}}$ axis defining the KPP1) and nine (including the seven satellite members of the KPP2 group aligned with the ${\vec{J}}_{\mathrm{stack}}$ axis defining it). Thus, a close relationship has been found in this case. In the case of the Aq-C^α system, we have 13 satellites in the KPP aligned with the ${\vec{J}}_{\mathrm{stack}}$ axis and 11 aligned with the ${\hat{e}}_{3}$ eigendirection. The number of satellites aligned with both of them at the same time is nine. This gives the conditioned fractions F(${\vec{J}}_{\mathrm{stack}}$ ∣ ${\hat{e}}_{3}$) = 0.82, F(${\hat{e}}_{3}$ ∣ ${\vec{J}}_{\mathrm{stack}}$) = 0.69, and F(${\hat{e}}_{i}$ ∣ no ${\vec{J}}_{\mathrm{stack}}$) = 0.10. Thus, we have found a relatively high (low) probability of a KPP satellite having its pole aligned with the ${\vec{J}}_{\mathrm{stack}}$ axis, in case it is (is not) aligned with the ${\hat{e}}_{3}$ or ${\hat{e}}_{1}$ principal directions.

Summing up, our results in this subsection indicate that the physical processes leading to the local CW development, in the case of these two zoom-in simulations, might have a significant impact on the dynamics of satellites, shaping their trajectories before the satellites are gravitationally bound to the central galaxy. Indeed, the same processes that cause the evolution of the LV principal directions across time induce the clustering of satellite orbital poles along some specific directions, hence contributing to the formation of kinematically persistent structures.

To characterize KPP orientations relative to the principal directions, either their normals, ${\vec{n}}_{\mathrm{KPP}}$, or the respective axes of maximum co-orbitation, ${\vec{J}}_{\mathrm{stack}}$, can be used. These are not equivalent analyses, as the normals ${\vec{n}}_{\mathrm{KPP}}$ are returned by the TOI analyses of the KPP satellite positions, while ${\vec{J}}_{\mathrm{stack}}$ determination uses the full six-dimensional phase space information. The angles between both vectors are given in Figure 8 (thick magenta lines), where we see that these vectors are highly aligned, except for the PDEVA-5004 KPP2 plane, where the alignment is more noisy. We also see that, as expected, ${\vec{n}}_{\mathrm{KPP}}$ and ${\vec{J}}_{\mathrm{stack}}$ are not perfectly parallel.

Zoom In Zoom Out Reset image size

Information on the alignments among the principal directions of LVs and co-orbitation axes, ${\vec{J}}_{\mathrm{stack}}$, is given in Figure 8 (thin solid lines) from T_uni = 7 Gyr onward for both the Aq-C^α and the PDEVA-5004 systems; see legends. As expected from the previous subsection, the ${\vec{J}}_{\mathrm{stack}}$ axis forms a small angle with the ${\hat{e}}_{3}$ direction in the case of Aq-C^α KPP and PDEVA-5004 KPP2 along most of the period analyzed. In other words, the direction of maximum co-orbitation in these two cases is close across time to that of maximum compression of the LV of reference here. The ${\vec{J}}_{\mathrm{stack}}$ axis for the PDEVA KPP1 plane, in turn, is close to the ${\hat{e}}_{1}$ direction.

Information on KPP orientations relative to the LV principal directions can be obtained from their respective normal vectors as positional planes, ${\vec{n}}_{\mathrm{KPP}}$. Figure 8 shows the angles between ${\vec{n}}_{\mathrm{KPP}}$ and the LV principal directions (thin dashed lines; see legends). Results for the PDEVA-5004 system show that the KPP1 (KPP2) structures, when considered as positional planes, are aligned with the ${\hat{e}}_{1}$ (${\hat{e}}_{3}$) principal directions. For the KPP identified in Aq-C^α , its normal vector aligns with ${\hat{e}}_{3}$ at high z and close to z = 0. In between, the alignment dims, and the angles both vectors form are similar, along some time intervals, to that ${\vec{n}}_{\mathrm{KPP}}$ forms with the ${\hat{e}}_{1}$ principal direction.

Put together, the results shown in Figure 8 indicate that the local environment flattening correlates with satellite kinematics rather than with their space positions. Indeed, kinematically coherent KPP satellites are characterized by their respective ${\vec{J}}_{\mathrm{stack}}$ axes, axes that highly align with the principal directions of the LVs, while this alignment gets worse when the positional planes (as characterized by their respective normals, ${\vec{n}}_{\mathrm{KPP}}$) are considered.

As suggested by Figure 7, the ${\hat{e}}_{3}$-structure plays an important role as a driver of orbital pole alignments in the two simulations analyzed here. Its formation and fate have been addressed in Section 6.3, in terms of the mass flows feeding it and, later on, piling up mass in a predominant filament that finally fades away in favor of halos. Let us now analyze the possible links between these processes and satellite plane formation. We do so by following the overall mass flows feeding the ${\hat{e}}_{3}$-structure from either side, using the satellites selected in this work as markers of the flows. Note that, in the following, we will use the term "${\hat{e}}_{3}$-plane" (in contrast to "${\hat{e}}_{3}$-structure") when referring to the mathematical plane that is perpendicular to the ${\hat{e}}_{3}$ direction and contains the center of mass of the LV.

Figure 9 shows some snapshots of the joint CW and KPP (proto)satellite ¹⁷ evolution in physical coordinates in a reference frame with the z-axis aligned with ${\hat{e}}_{3}$, the y-axis aligned with ${\hat{e}}_{2}$, and the x-axis aligned with ${\hat{e}}_{1}$, centered on the host galaxy formation site. We see that satellite-to-be mass elements (red points in the figure) come from small, dense subvolumes that at high redshift are small walls or filaments. Satellites are eventually formed by feeding from these mass elements. In some cases, these mass elements disappear in favor of collapsed satellites, which in these cases are essentially born in isolation. For example, the red structure located at Z ∼ 180 and X ∼ −200 kpc (in the first panel) evolves into a unique, isolated satellite. In other cases, satellites follow the filament they are formed from in a wide sense. These snapshots illustrate the idea that the global effects of the forces and torques produced by the whole CW on a particular (proto)satellite drive its incorporation to the ${\hat{e}}_{3}$-wall in the Aq-C^α simulation, with the orbital pole alignment results shown in Figure 7.

Zoom In Zoom Out Reset image size

The question arises of why some satellites are not part of KPPs. As already mentioned, in Figure 9, we see that different satellites reach the ${\hat{e}}_{3}$-wall under different circumstances; some flow through filaments along with mass elements that eventually feed into the formation of the ${\hat{e}}_{3}$-structure, others are part of a small-scale wall that merges with the ${\hat{e}}_{3}$-structure in a parallel manner, and finally, others reach the ${\hat{e}}_{3}$-structure freely, having cannibalized along the way the mass of the CW element they were initially embedded in. The natural evolution of the CW dynamics therefore leads to satellite trajectories with different characteristics as they reach the ${\hat{e}}_{3}$-structure (e.g., distance to the ${\hat{e}}_{3}$-plane, angle with respect to the plane), which in turn affects how the orbital pole alignment for each particular satellite comes about.

Qualitatively, we observe that satellites with trajectories that are initially closer to the ${\hat{e}}_{3}$-plane tend to run obliquely relative to it and are smoothly captured by the wall. Some satellites are almost in-plane from high redshift. On the other hand, satellites coming from further away fall onto the ${\hat{e}}_{3}$-plane at higher velocities, more perpendicularly to it, and they tend to maintain a velocity component normal to the ${\hat{e}}_{3}$-plane, such that their orbital poles are contained within this plane. In addition, satellites with a low impact parameter (i.e., low minimal distance with respect to the proto-host-galaxy center of mass) never suffer important specific angular momentum sJ_orb gains (see Figure 1), resulting in small apocenter orbits, with a high probability of suffering disturbing effects from the central regions of the system and so changing their orbital pole directions.

We have studied the evolution of distances from each satellite to the ${\hat{e}}_{3}$-plane. The left panel in Figure 10 shows the time evolution of the median distances for the different satellite populations in Aq-C^α , together with the corresponding 25th–75th percentiles (in comoving coordinates, i.e., removing the effects of the background Universe expansion). Two different regimes are clearly visible in either the KPP or non-KPP populations. At early times, distances decrease rapidly, an indication of strong mass flows in the ${\hat{e}}_{3}$ principal direction, feeding the ${\hat{e}}_{3}$-structure. In a second phase, median distances show fluctuations with essentially constant amplitude (in physical coordinates), except for nonlinear effects around ∼10 Gyr due to a pericenter accumulation event, affecting mainly the non-KPP population behavior. This constancy indicates that satellites do not feel the background expansion anymore. The two regimes define a separation interval of time, T_dist,plane, whose specific values depend on the satellite population; see below. For the sake of clarity, we define this value as the time when the distance curve reaches a minimum for the first time. This behavior, namely, a rapid distance decrement leading to a system with roughly constant size, decoupled from the background expansion, is reminiscent of collapse events suffered by halos as described by the spherical collapse model (Padmanabhan 1993) plus the ensuing violent relaxation (Lynden-Bell 1967) process leading to an equilibrium configuration. In the absence of any theory or model to describe the statistical fate of the particles involved in the collapse toward a wall (see footnote 12), we focus on the empirical characteristics just mentioned of the ${\hat{e}}_{3}$-structure behavior in a time interval of around ∼4 Gyr, and due to the reminiscences found with some characteristics of halo collapse in three dimensions, we will hereafter refer to this event as "${\hat{e}}_{3}$-structure collapse."

Zoom In Zoom Out Reset image size

A second interesting result Figure 10 reports on is that KPP and non-KPP satellites sample the two phases mentioned above differently. Differences before and along the ${\hat{e}}_{3}$-structure collapse are of particular interest here. Non-KPP satellites come from further away than KPP satellites before T_dist,plane (as Figure 11 illustrates), and, on average, non-KPP satellite members reach the ${\hat{e}}_{3}$-plane later than KPP satellite members as well. The specific values of T_dist,plane for the two satellite populations are given in Table 2.

Zoom In Zoom Out Reset image size

Additionally, by following (proto)satellite spatial trajectories, we have found that many non-KPP satellites trace mass flows that are mostly perpendicular to the ${\hat{e}}_{3}$-plane, converging onto the plane and feeding it, while KPP satellites' trajectories tend to run more obliquely relative to the ${\hat{e}}_{3}$-plane. Figure 11 shows some illustrative examples of these two different satellite trajectory behaviors.

As seen in Figure 7, it takes some time to have the orbital poles of KPP satellites aligned with the ${\hat{e}}_{3}$ directions (indeed, ${T}_{\mathrm{align},\,{{\rm{e}}}_{3}}\simeq$ 4.5 Gyr; see Section 7.1). To reach the aligned configuration that keeps roughly stable in time, a necessary condition is that the ${\hat{e}}_{3}$ principal direction is frozen. This happens very early (at ${T}_{\mathrm{dir},{{\rm{e}}}_{3}}\sim$ 2 Gyr; see Table 2). A second condition is that those not-yet-in-plane incoming satellites are placed in-plane.

Within this scheme, KPP satellite members would be those placed in-plane from very high redshift plus those having their orbital poles sucessfully bent while they are incorporated into the ${\hat{e}}_{3}$-structure. This "pole-bending" effect comes about as a result of the particular dynamical evolution of the CW, where satellites are just tracers of its flowing mass elements. Indeed, the trajectories of the KPP satellites shown in Figure 11 are evidently bent toward the ${\hat{e}}_{3}$-wall (approximately the X–Y plane) when collapsing onto it at T_dist,plane (see crosses).

Non-KPP satellites can have different origins: (i) those coming from further away than KPP satellites, infalling onto the ${\hat{e}}_{3}$-plane at late times and thus not having had enough time to complete a full orbit around the host yet; (ii) those with low impact parameters and hence small apocenters, like satellite #343 in Figure 11; and (iii) those with almost perpendicular infall onto the ${\hat{e}}_{3}$-plane but with a larger impact parameter, such that their orbital poles are parallel to the plane and aligned with the ${\hat{e}}_{2}$ axis (3 out of 21 non-KPP satellites, not classified as a second plane due to their low number); see satellite #347 in Figure 11.

Some considerations are in order concerning the effects that satellite motions within the ${\hat{e}}_{3}$-structure may have on satellite pole alignments. As already mentioned, between T_dist,plane and T_vir, the alignments between KPP satellite poles and the ${\hat{e}}_{3}$ direction improve somewhat, but to a much lesser extent than during the previous period of ${\hat{e}}_{3}$-structure formation. Thus, it would seem that, after T_dist,plane, the secondary collapse phase leading to the formation and evolution of this prolate structure does not have relevant effects on the alignments. Indeed, as already mentioned, these alignments keep overall constant until z ∼ 0 (see Figure 7), in contrast to early CW walls and filaments, that mostly vanish in favor of halos (see Figure 4).

In the case of PDEVA-5004, the right panel of Figure 10 indicates that the median distances to the ${\hat{e}}_{3}$-plane for the different satellite populations show similar trends as those mentioned above for the Aq-C^α simulation, i.e., two clearly distinguished phases, the first of them of rapid distance decrease. The main dissimilarity with the Aq-C^α case is that in PDEVA-5004, the ${\hat{e}}_{3}$-structure is a prolate structure and not a wall.

Another particularity of satellites in PDEVA-5004 is that, in many cases, protosatellite trajectories are initially approximately parallel to the ${\hat{e}}_{3}$ principal direction and have had their respective orbital poles aligned since very high redshift with the ${\hat{e}}_{1}$ principal direction. This is the origin of the KPP1 satellite system.

The KPP1 and KPP2 satellite groups show different T_dist,plane values (see Table 2): KPP1 satellite members reach, on average, the ${\hat{e}}_{3}$-plane ∼2 Gyr later than KPP2 satellites, during a phase of the evolution when dynamical events happen rapidly. Again, KPP2 satellite members are those whose trajectories have been succesfully bent by T_dist,plane ∼ 3 Gyr (and by the same reasons advocated previously, i.e., the CW dynamics) or those that are already aligned at high redshift. In addition, KPP1 satellite members come from further away than KPP2 members, with their poles already clustered in the ${\hat{e}}_{1}$ direction, and this clustering is maintained across ${\hat{e}}_{3}$-structure collapse.

As for non-KPP satellites, similarly to the Aq-C^α case, most of them reach the ${\hat{e}}_{3}$-plane with a low impact parameter relative to the host center or have been recently captured by the host.

Different timescales for the galaxy–satellite system formation and evolution have emerged in this work, summarized in Table 2. Apart from the halo turnaround, T_ta,halo, and virialization, T_vir, timescales (Table 1) coming from the spherical collapse model, we have pointed to and made specific definitions for timescales relative to the LV evolution, such as the Universe age when

• i)  
  the ${\hat{e}}_{3}$ principal vector, corresponding to the direction of the dominant compression flow of matter, gets its direction fixed, ${T}_{\mathrm{dir},{{\rm{e}}}_{3}}$ (see Figure 5); and
• ii)  
  the rapid high-redshift decrement of the minor principal axis, c(t), stops, ${T}_{\mathrm{shape},{{\rm{e}}}_{3}}$ (see Figure 6).

We have also focused on timescales involving the (proto)satellites of the different populations, in which case, when possible, we show the population median values and the 25th–75th percentiles of the Universe age when

• iii)  
  the ith satellite becomes aligned with the axis ${\vec{J}}_{\mathrm{stack}}$ of maximum satellite co-orbitation determining the KPP plane the satellite belongs to, T_{cluster,Jstack} (see Figure 2);
• iv)  
  the ith satellite becomes aligned with the ${\hat{e}}_{3}$ principal vector (for satellites in KPP or KPP2) or the ${\hat{e}}_{1}$ eigenvector (satellites in KPP1); see Figure 7 and ${T}_{\mathrm{align},{{\rm{e}}}_{i}}$ entry in Table 2; and
• v)  
  there is a broad minimum for the first time in the median vertical distances to the ${\hat{e}}_{3}$-plane of the different satellite populations shown in Figure 10, T_dist,plane. From this age onward, the median distances stay roughly constant (in physical coordinates) and are much lower than the distances before this age is reached.

For the sake of the discussion in this subsection, we now report on two more timescales involving satellite populations as well, i.e., their medians and percentiles of the Universe age when

• vi)  
  the ith satellite distance to the host center of mass is, for the first time, smaller than the respective virial radii at that age (infall time), T_sat,infall; and
• vii)  
  the distance from the ith (proto)satellite to the (proto)halo center of mass reaches the first maximum (in physical coordinates), T_sat,apo1. This is essentially a turnaround timescale for the ith satellite, marking the beginning of its decoupling from the expansion of the background Universe.

As discussed in the previous subsections, in the physical process behind satellite orbital pole alignment with the ${\hat{e}}_{3}$ principal direction, a timescale stands out in the two zoom-in simulations studied: T_dist,plane, the collapse timescale for the ${\hat{e}}_{3}$-structure. Its values are ∼3.5 Gyr for the KPP system in the Aq-C^α simulation and ∼3 Gyr for the KPP2 system in PDEVA-5004. According to Table 2, T_dist,plane is roughly coeval to the clustering timescales, or the Universe age when KPPs are established: T_{cluster,Jstack} or ${T}_{\mathrm{align},{{\rm{e}}}_{3}}$. This is expected, because an aligned satellite orbits within the ${\hat{e}}_{3}$-structure; thus, its distance to the ${\hat{e}}_{3}$-plane must be low. This table also indicates that T_dist,plane is coeval to ${T}_{\mathrm{shape},{{\rm{e}}}_{3}}$, a timescale marking—within the accuracy of their determinations—the Universe age when the initially strong mass flows normal to the ${\hat{e}}_{3}$-structure, which feed it and eventually cause its collapse, slow down.

On the other hand, the infall of satellites onto the halo happens after their capture by the ${\hat{e}}_{3}$-structure (particularly so for the Aq-C^α simulation), while the satellites reach first apocenter at times much earlier than this (see Table 2 T_sat,infall and T_sat,apo1 entries). Thus, none of these processes seem to be related to the origin of KPPs.

Alignments with the principal directions set in very early for KPP members in both simulations. For PDEVA-5004 KPP2 satellites, this happens earlier than for satellites in the Aq-C^α system. As for PDEVA-5004 KPP1 satellites, they come with the main mass flow normal to the ${\hat{e}}_{3}$-plane; therefore, their clustering (defined by the ${\hat{e}}_{1}$ direction) is already established at very high redshift.

To finish, let us mention the timescale for the establishment of the maximum ordered disky motion (circular and in-plane) of satellites in the two simulations analyzed here, ${T}_{\max ,\mathrm{rot}}$ (Section 4.2). We see that the morphological kinematic κ_rot parameter reaches its maximum shortly after T_dist,plane for both the KPP and the KPP2 satellite samples, in coherence with our previous interpretations in this section. To complete the scheme, for the KPP1 satellite members, the κ_rot parameter values have been high since very high redshift, in coherence with the very early alignments of their orbital poles with the ${\hat{e}}_{1}(t)$ principal direction.

The aim of this paper is to make a contribution to understanding the origin or physical processes behind the formation of persistent, kinematically coherent planes of satellites (KPPs), i.e., sets of satellites with fixed identities co-orbiting around their host galaxy, whose orbital poles are conserved and clustered across long cosmic time intervals and whose positions form good-quality positional planes.

Paper III identified such KPPs in their analyses of two cosmological, zoom-in hydrodynamical simulations where a system of some 30–35 satellites orbit around an MW-mass type galaxy with an extended thin gaseous and stellar disk. The two simulations differ in their initial conditions, the subgrid physics, and the methods to integrate both the gravitational and the hydrodynamical equations. In Paper III, simulations are analyzed from halo virialization time T_vir to z = 0. In both simulations, a relatively high fraction of the satellites have been found to be kinematically organized (a maximum of ∼60% and ∼80%, along some time intervals, in the Aq-C^α and PDEVA-5004 systems, respectively).

The specific aim of this paper is to elucidate which physical processes cause the satellite orbital angular momentum (${\vec{J}}_{\mathrm{orb}}$) direction clustering at early times. Using the same two zoom-in simulations as in Paper III but extending the analyzed period back in time until z_high = 8.45 and 10.00 for the Aq-C^α and PDEVA-5004 systems, respectively, we focus on the overall evolution of the CW around the satellite-to-be and galaxy-to-be objects from z_high until z = 0. We therefore follow the progenitors of the low-redshift satellites and host systems well within the fast phase of the system assembly and within the local (i.e., around the forming system) CW they are embedded in.

By following satellites back in time, we find that, in most cases, their progenitors gain a specific orbital angular momentum magnitude as predicted by the TTT from T_uni ∼ 2 to ∼4 Gyr. The directions of the orbital angular momentum (i.e., the so-called orbital poles) of KPP satellite progenitors are also conserved from very early, while this is not the case for many satellites outside KPPs (Figure 1). Our analysis here indicates that clustering of orbital poles occurs already at high redshift; see Table 2. An analysis by means of the morphological kinematic κ_rot parameter (see, for example, Sales et al. 2012) indicates that, from very early times, the collections of KPP satellites represent, kinematically, "disky" systems, with a high fraction of their kinetic energy coming from in-plane and almost circular motion within KPPs, while systems outside these structures are more spheroidal (Figure 3).

To elucidate how the aforementioned clustering came about, we analyze satellite pole evolution as part of the CW they are embedded in through the LV deformation method (see Robles et al. 2015). For each simulation, we mark the particles that at z_high are within a sphere of radius R_LV = K · r_vir,z=0/(1 + z_high) (K = 15 and 20) centered at the protogalaxy center and follow their trajectories forward in time up to z = 0. The volumes these particles span at each time are the so-called LVs. We analyze the evolution of their principal directions ${\hat{e}}_{i}(t)$, with i = 1, 2, 3, and their principal axes a(t) > b(t) > c(t) through the reduced TOI method (Cramér 1999).

The a(t), b(t), and c(t) functions inform us about the shape deformations of the LV, including how quickly they happen. The general result is that while a(t) grows, c(t) decreases, in most cases monotonously, but with some stagnation periods (Figure 6). The b(t) axis keeps roughly constant in PDEVA-5004 and decreases in Aq-C^α after T_uni ∼ 4 Gyr. As for axis ratios, c(t)/a(t) decreases rapidly up to T_uni ∼ 4 Gyr in the Aq-C^α LV, causing a quick flattening of its initially spherical shape into a wall-like CW structure (Figure 4). Then, the decrement of the b(t)/a(t) ratio takes over, transforming the LV shape from oblate to triaxial and finally prolate. Ratio changes in the PDEVA-5004 LV are such that its shape is always triaxial. The principal directions also freeze out, meaning that the direction of overall maximum compression ${\hat{e}}_{3}$ does not change (within a threshold) after its freezing out at ${T}_{\mathrm{dir},{{\rm{e}}}_{3}}$, and, consequently, from this moment onward, the overall maximum compression takes place along a fixed direction (Figure 5). All three principal directions freeze out very early (T_freeze ≃ 4 and 2 Gyr for Aq-C^α and PDEVA-5004, respectively). In this way, we witness the emergence of a kind of global LV "skeleton," such that for T_uni > T_freeze, overall anisotropic mass rearrangements of LV particles preferentially occur along fixed directions.

To elucidate the role that the local CW development around the forming system has in driving pole clustering, we have analyzed the alignments between the LV principal directions and the satellite orbital poles across cosmic time (Figure 7). We find that, for KPP satellites in the Aq-C^α system, a clear alignment signal stands out with the ${\hat{e}}_{3}$ axis after ${T}_{\mathrm{align},{{\rm{e}}}_{3}}\sim$ 4.5 Gyr (see Table 2); i.e., the orbital poles of KPP satellites tend to be parallel to the LV's direction of maximum overall compression. Some satellites in the KPP show alignments as early as at T_uni ∼ 2 Gyr. For those that do not, the orbital poles are bent efficiently between T_uni ∼ 2 and 4.5 Gyr. This is roughly coeval to the timescale for ${\hat{e}}_{3}$-structure collapse. Satellites outside KPP structures do not show any particular alignments with any LV principal direction.

In the PDEVA-5004 system, the pole alignment signal of KPP2 satellites with the ${\hat{e}}_{3}$ axis after ${T}_{\mathrm{align},{{\rm{e}}}_{3}}\sim$ 3.5 Gyr is even clearer than in the Aq-C^α case. KPP1 satellite poles have been well aligned with the ${\hat{e}}_{1}$ axis since at least T_uni ≃ 2.0 Gyr. No alignment signals show up for satellites not belonging to KPPs.

We study the evolution of protosatellite mass elements in relation to the evolving local CW (Figures 9 and 10). Our findings show that KPP satellites aligned with ${\hat{e}}_{3}$ present closer, more oblique trajectories relative to the ${\hat{e}}_{3}$-plane, allowing for "pole bending" by the same forces and torques that drive the evolution of the local CW dynamics (Figure 11). Satellites with poles aligned with the ${\hat{e}}_{1}$ axis show trajectories parallel to ${\hat{e}}_{3}$ and are not bent. Finally, non-KPP satellites can present different origins, with many showing low impact parameters relative to the host center and hence a high probability of suffering disturbing phenomena.

To our knowledge, this paper represents the first time that the effects of the CW as a driver of the satellite orbital pole organization into KPPs is analyzed in some detail through numerical simulations and where comparison with previous results, including observational ones, can be easily made.

A few works have attempted to address the characteristics and evolution of co-orbiting satellite planes or their origin as connected to the large-scale structure they are embedded in.

For example, Shao et al. (2019) used the EAGLE-100 volume to identify "MW-like-orbit" satellite planes, i.e., narrow planes formed by the 11 most massive satellites around MW-mass halos, where eight of them show a high degree of coplanarity (a small dispersion of their orbital poles) at z = 0. In that paper, their aim was to look for alignments with the principal directions of the host halos resulting from a TOI analysis. They noted that the degree of co-orbitation of their subsets of eight coherent satellites selected at z = 0 is better at present than at earlier times and suggest it is driven by halo torques after infall. Our findings are consistent with theirs in that the collimation and clustering of satellite orbital poles in two KPPs improves with time (see alignments with the ${\hat{e}}_{3}$ direction in both simulations; Figure 7). Finally, it is worth noting that they found a wide variety of times at which these eight MW-like-orbit satellites started to show co-orbitation, with some setting in early, while others do later on (see the wide range implied by percentiles in ${T}_{\mathrm{align},{{\rm{e}}}_{3}}$ in Table 2). The higher fraction of satellites that are established much later on relative to our results could come from other possible channels for orbital pole clustering enhancement, such as LMC-like group infall or tides from aspherical halos.

In a more general perspective, our results are also consistent with—and provide an explanation for—those from Libeskind et al. (2014) and Dupuy et al. (2022), who found a preferred direction of subhalo infall onto halos. These works show that subhalos are mainly incorporated onto halos along a direction that is contained within the plane orthogonal to the direction of fastest collapse and that aligns with the spines of filaments. Following these predictions, Libeskind et al. (2015, 2019) tested the possible alignments between the observed satellite planes in the local Universe and the principal directions of the cosmic density field as reconstructed from the CosmicFlows-2 peculiar velocity survey (Tully et al. 2013). Similarly, Xu et al. (2023) also suggest a possible connection between the presence of a rotating plane of satellites in TNG50 and the large-scale sheet structure it is embedded in.

Our work shows that, in the two simulations analyzed in this paper, kinematically coherent satellites in fact gain their common dynamics much before they reach the halo by tracing the mass flows as mass collapses into the CW elements (Zel'Dovich 1970; Shandarin & Zeldovich 1989). We emphasize, therefore, that here it is but the initial step of anisotropic mass collapse, together with the initial location of protosatellites relative to the skeleton of CW collapsed structures, which leaves an imprint on the dynamics of satellites.

Indeed, we find that—different from suggestions from previous works (Goerdt et al. 2013; Buck et al. 2015; Ahmed et al. 2017)—kinematic planes are not only driven by filamentary infall, as KPP satellites do not always reach their stationary positions via one or a few strong filaments.

We note that the principal directions of collapse identified in the previous works result from velocity shear tensor analyses at z = 0. Indeed, most methods to analyze the CW evolution and classify its elements (see, e.g., Hoffman et al. 2012; Cautun et al. 2013; Libeskind et al. 2018, and references therein) are local ones, where the tensorial tools used are defined on points, needing a smoothing procedure to enable calculation.

Conversely, to study the evolution of the local environment of sites where galaxy systems are to form, our perspective is rather global. We use here a simple method that tracks the average large-scale deformations of an initially spherical LV. The method provides the accumulated deformations the LV suffers along time intervals between the different snapshots the simulation provides. This method is simpler than the usual local, tensorial methods in that the LV evolution is described through three principal directions (that happen to freeze out at very early times here) and three time functions, the principal axes a(t), b(t), and c(t), from which only two are independent. The method accurately catches the evolution of the local environment of galaxy formation sites. In our analysis here, the method singles out an overall direction of maximum compression of the mass flows whose relevance here has already been mentioned. In practice, this direction of maximum compression is the same as those returned by the usual methods, as the velocity shear tensor analysis. ¹⁸ While the velocity shear tensor is a more appropriate scheme when trying to detect and classify CW structures (see Cautun et al. 2013; Libeskind et al. 2018), both formalisms allow one to identify the main directions of mass flows.

These are the conclusions of this work concerning the origin of KPPs, according to the two simulations analyzed in this paper.

• 1.  
  The formation of KPPs is closely related to the early anisotropic collapse of mass in ΛCDM that drives the evolution of the CW: the same physical processes behind the emergence of the CW at high redshift are behind the emergence of clustering of (proto)satellite orbital poles, that is, behind the formation of KPPs at high redshift. The initial location of protosatellites relative to the sites of early CW collapse also plays a role.
• 2.  
  A timescale stands out for the establishment of KPP satellite orbital pole clustering, T_dist,plane, the Universe age when satellites' distances to the plane defined by the direction of overall compression of the local mass distribution (i.e., ${\hat{e}}_{3}$ direction) become overall roughly constant. This is similar to a collapse event. This event occurs well before the median timescale for satellite infall onto their host halo.
• 3.  
  KPP member satellites characterized by pole alignments with the ${\hat{e}}_{3}$ principal direction are those already aligned at high redshift plus those whose orbital poles have been succesfully bent as they are incorporated into the ${\hat{e}}_{3}$-structure when it collapses. The latter's trajectories at high redshift tend to be oblique relative to the ${\hat{e}}_{3}$-plane. This is the case of the so-called Aq-C^α and PDEVA-5004 KPP2 groups. Additionally, in the case of Aq-C^α KPP satellites, they collapse earlier onto the ${\hat{e}}_{3}$-structure and are closer to the ${\hat{e}}_{3}$-plane at given times than non-KPP satellites.
• 4.  
  KPP satellite orbital pole alignments not only occur along the direction of maximum compression ${\hat{e}}_{3}$ but can also appear along the other principal directions. Such is the case of a sizable subset of satellites in the PDEVA-5004 simulation (the KPP1 system), whose orbital poles are aligned with the ${\hat{e}}_{1}$ direction, and of some satellites in the Aq-C^α simulation, with orbital poles aligned with the ${\hat{e}}_{2}$ eigenvector.
• 5.  
  In PDEVA-5004, the CW dynamics leads to a triaxial mass distribution (rather than to a wall-like structure as in Aq-C^α ). Those PDEVA-5004 satellites characterized by pole alignments with the ${\hat{e}}_{1}$ direction tend to have trajectories that, at high redshift, follow the direction of overall maximum compression. These satellites do not suffer from orbital pole bending when the ${\hat{e}}_{3}$-structure collapses.

We would like to emphasize that just two simulations have been analyzed. However, these two simulations differ in multiple aspects. In Section 2, we presented the differences between the two simulations regarding their hydrodynamics, and we pointed out that the initial conditions and the subresolution physics models differ as well. The possibility of KPP formation in both cases implies that their origin must be driven by the more fundamental common physical processes of structure evolution in a ΛCDM cosmological context and less dependent on the details of galaxy formation modeling. We have shown not only that KPPs can form in ΛCDM simulated disk galaxy systems but, importantly, that their existence is a natural consequence of ΛCDM's prediction for large-scale mass flows at high redshift shaping the local CW structure. In the two zoom-in simulations analyzed here, KPPs are the result of the same dynamics acting on protosatellites' mass elements placed at particular locations and/or endowed with particular kinematic characteristics.

We note that other channels for KPP enhancement at low redshift are possible as well, for example, through the late capture of a satellite with its own system of subsatellites (see Paper III). On the other hand, satellite interactions within the inner regions of halos could destroy kinematic coherence.

Finally, we want to stress that our scientific conclusions are an interpretation of results from only two simulations, which exhibit a correlation between LVs and KPPs. Further work involving extending our analysis to a broader sample of simulations is therefore needed in order to assess the frequency of finding KPPs and to robustly conclude that this is a generic feature of KPPs in ΛCDM. In particular, analyses of large-volume simulations of galaxy formation might shed light on how frequently the different channels for satellite plane formation appear throughout cosmic evolution.