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

Kinematically persistent planes (KPPs) of satellites are fixed sets of satellites co-orbiting around their host galaxy, whose orbital poles are conserved and clustered across long cosmic time intervals. They play the role of “skeletons,” ensuring the long-term durability of positional planes. We explore the physical processes behind their formation in terms of the dynamics of the local cosmic web (CW), characterized via the so-called Lagrangian volumes (LVs) built up around two zoom-in, cosmological hydro-simulations of Milky Way–mass disk galaxy + satellites systems, where three KPPs have been identified. By analyzing the LV deformations in terms of the reduced tensor of inertia (TOI), we find an outstanding alignment between the LV principal directions and the KPP satellites’ orbital poles. The most compressive local mass flows (along the eˆ3 eigenvector) are strong at early times, feeding the so-called eˆ3 -structure, while the smallest TOI axis rapidly decreases. The eˆ3 -structure collapse marks the end of this regime and is the timescale for the establishment of satellite orbital pole clustering when the Universe is ≲4 Gyr old. KPP protosatellites aligned with eˆ3 are those whose orbital poles are either aligned from early times or have been successfully bent at eˆ3 -structure collapse. KPP satellites associated with eˆ1 tend to have early trajectories already parallel to eˆ3 . We show that KPPs can arise as a result of the ΛCDM-predicted large-scale dynamics acting on particular sets of protosatellites, the same dynamics that shape the local CW environment.


Introduction
The discovery that satellites orbiting around the Milky Way (MW) and Andromeda (M31) have a highly anisotropic distribution has been long considered to be one of the most challenging issues for ΛCDM (see Bullock & Boylan-Kolchin 2017;Pawlowski 2018, for a review).
Most MW satellites define a "vast polar structure" (VPOS) of dwarf galaxies with respect to the Galactic disk (see Pawlowski et al. 2012Pawlowski et al. , 2013;;Kroupa 2015), with around 40% of their orbital poles aligned with the normal vector to the plane of satellites (Fritz et al. 2018; see Santos-Santos et al. 2020a, hereafter Paper I), suggesting that the VPOS is a robust positional structure and could be rotationally supported.
Satellites orbiting around M31 have also been discovered to be anisotropically distributed around it (Koch & Grebel 2006;McConnachie & Irwin 2006;Metz et al. 2007), with almost half the population (15 out of 27 satellites) forming a thin plane in positions, referred to as the "Great Plane of Andromeda," and almost edge-on from our perspective.Moreover, in Paper I, it has been shown that a second positional plane, roughly perpendicular to the former, shows up in M31.Finally, positionally flattened satellite structures have also been detected in other major galaxies beyond the Local Group (e.g., Chiboucas et al. 2013;Ibata et al. 2014Ibata et al. , 2015;;Tully et al. 2015;Müller et al. 2016Müller et al. , 2017Müller et al. , 2018Müller et al. , 2021;;Heesters et al. 2021;Martínez-Delgado et al. 2021;Paudel et al. 2021).
These recent observations have opened interesting debates on the issue of the positional planes of satellites.Indeed, different authors have studied ΛCDM simulations, finding that, while their presence is quite unusual (Libeskind et al. 2005(Libeskind et al. , 2009;;Lovell et al. 2011;Wang et al. 2013;Bahl & Baumgardt 2014;Cautun et al. 2015;Forero-Romero & Arias 2018), positionally detected planes of satellites can be found, in some cases even showing kinematic coherence similar to the planes found in the Local Group (Buck et al. 2015;Gillet et al. 2015;Ahmed et al. 2017;Maji et al. 2017;Garaldi et al. 2018;Shao et al. 2019;Samuel et al. 2021;Förster et al. 2022;Pham et al. 2023).
In many cases, however, these positional planes have been found to be unstable or transient structures; i.e., the membership of at least a fraction of satellites is lost within short timescales (Bahl & Baumgardt 2014;Buck et al. 2016;Lipnicky & Chakrabarti 2017;Maji et al. 2017;Zhao et al. 2023).In line with these results, Santos-Santos et al. (2020b, hereafter Paper II) have found that important fluctuations occur in the properties of positional planes as a function of time, presumably caused by the loss of a fraction of satellites that leave the planar structure on short timescales, while other transient satellites join it.These authors have also found that, at each simulation time step, only a fraction of satellites are in coherent co-orbitation11 within the planes, indirectly suggesting the possibility of a kinematic skeleton in positional planes.
This possibility has been explored by Santos-Santos et al. (2023, hereafter Paper III), who made a kinematic analysis of the satellite samples analyzed in Paper II, from the halo virialization time, T vir , onward.Specifically, two hydrodynamical, zoom-in MW-mass systems were studied.By focusing on the satellite orbital pole conservation and clustering, they identified the socalled kinematically persistent planes (KPPs) of satellites.These are groups of satellites whose identities are the same along extended time intervals and whose orbital poles are conserved and clustered in the same direction along them.In Paper III, it was numerically shown that KPP satellites, on the one hand, and satellite members of the thinnest positional planes (i.e., those with the lowest minor-to-major axes ratios c/a) among those with a fixed satellite fraction, on the other hand, share the same threedimensional space configuration.This result numerically proves that the positional planes include nonkinematically coherent satellites as well, i.e., satellites that are lost to the positional plane configuration and temporarily replaced by other transient-member satellites, as mentioned.In this way, KPPs play the role of a kind of skeleton, shaping long-lasting planar structures whose satellite membership fluctuates along time, except for the kinematically coherent ones.
We see that the key point of KPP structures is the clustering of a fraction of the host satellites' orbital poles and its persistence along time.But questions remain about what is causing this clustering, when this clustering is established, and whether or not the cosmic web (CW) development has some role in answering to the two previous questions.
Indeed, the processes behind the clustering of orbital polesand hence behind the origin of KPPs-remain unclear.Different approaches have been proposed in order to explain such phenomenon, such as group capture of satellites onto the central galaxy (Lynden-Bell & Lynden-Bell 1995;D'Onghia & Lake 2008;Li & Helmi 2008), satellites being formed from tidal dwarf galaxies in ancient gas-rich mergers (Kroupa et al. 2010;Hammer et al. 2013;Kroupa 2015), capture of satellites during host mergers (Angus et al. 2016;Smith et al. 2016), or the effect of aspherical halo tides at increasing the orbital pole collimation of satellites orbiting inside them or at populating KPPs (see, e.g., Shao et al. 2019;Wang et al. 2020).
Additionally, partly based on early observations of the alignment between the VPOS and the Large Magellanic Cloud (LMC; see Kunkel & Demers 1976;Lynden-Bell 1976), some authors (e.g., Garavito-Camargo et al. 2021;Samuel et al. 2021) have suggested that the LMC infall onto the MW might help to explain the existence of the VPOS (see, however, Pawlowski 2021).In Paper III, it was shown that the late infall of an LMC-like group of satellites is not needed in order to have kinematically coherent satellite planes.However, this presence might help to enhance the fraction of satellites in coherent coorbitation and, more particularly, the ratio of those rotating in one sense over those rotating in the contrary sense within KPPs.
Although the previously mentioned processes might have been operative along cosmic evolution, and they could have concurrently contributed to satellite plane formation, other authors' approaches have stressed the role of the CW evolution.It is known that the mass elements that currently form galaxies were organized at high redshift as a CW, whose emergence and evolution are analytically described via the Zel'dovich approximation (ZA; Zel'dovich 1970) and its extension to the adhesion model (see, e.g., Gurbatov et al. 1989;Kofman et al. 1992;Gurbatov et al. 2012, and references therein), as well as via numerical simulations (Cautun et al. 2014, and references therein).These works show that the morphology of the CW comes from a hierarchical, multiscale, anisotropic collapse, where large-scale flattened structures, frequently containing coplanar filaments (see, e.g., Aragón-Calvo et al. 2010), are but one of its elements.It is worth noting that when the collapse of a CW filament or a wall is mentioned,12 it is generally not a simple caustic formation that is meant but rather a multiscale, complex, nonrelaxed structure, made in turn of different smaller-scale CW elements and so on.We empirically know that, at a given scale, these are morphologically transient structures vanishing in favor of halos, where mass piles up.
There is increasing evidence that the CW morphological development shapes some halo and galaxy properties.Using numerical simulations, Libeskind et al. (2014) and Kang & Wang (2015) have shown that the major axes of dark matter (DM) halos are well aligned with the slowest-collapsing directions of the LSS density field Hessian (i.e., their major axes or filaments; see also Porciani et al. 2002a for protohalo alignments).Other studies reach similar conclusions analyzing simulations using different techniques (Vera-Ciro et al. 2011;Shao et al. 2016;Cataldi et al. 2023).See also Wang et al. (2020) for similar results obtained in the Sloan Digital Sky Survey DR12 data set analysis.
Another basic example is the role played by the CW in the acquisition of angular momentum by gas as it travels toward the galaxy formation region (Pichon et al. 2011;Codis et al. 2015;Kraljic et al. 2020).According to the so-called tidal torque theory (TTT; see, e.g., Peebles 1969;Doroshkevich 1970;White 1984;Schäfer 2009, for a review), angular momentum acquisition at very high redshift is the result of the misalignment between the inertia tensor and the shear tensor due to tidal forces.The extended TTT, revised to include the CW anisotropic configuration (Codis et al. 2015;Kraljic et al. 2020), predicts that massive galaxies have their spin perpendicular to filaments, while for low-mass galaxies, it is aligned with the filament.These results have been confirmed through simulations (e.g., Dubois et al. 2014;Welker et al. 2014) and observational data (e.g., Welker et al. 2020).Similar results involving halos are presented by Kang & Wang (2015) and references therein (Codis et al. 2012;Aragon-Calvo & Yang 2014;Welker et al. 2014).Alignment of halos with the filaments of the CW have been analyzed by Ganeshaiah Veena et al. (2018, 2019, 2021) in different largevolume simulations.
Mass flows within CW structures also have an impact on the distribution of matter around galaxies and halos, especially on satellite distributions.For example, some studies focus on the anisotropic character of subhalo/satellite capture along filaments (Benjouali et al. 2011;Libeskind et al. 2014;Wang et al. 2014;Kang & Wang 2015;Tempel et al. 2015;Dupuy et al. 2022), as well as its possible consequences for satellite systems' shapes (Tempel et al. 2015;Wang et al. 2020) or orbital coplanarity (Benjouali et al. 2011;Goerdt et al. 2013;Buck et al. 2015).Other authors analyze alignments between the principal directions of the inertia ellipsoids of satellite systems, on the one hand, and those of a few-R 200 -scale shells around their respective host galaxies (Shao et al. 2016) or the axes of slowest collapse in the matter distribution at larger scales (Libeskind et al. 2015(Libeskind et al. , 2019)), on the other hand.An analysis of dark matter only (DMO) simulations led Libeskind et al. (2012) to find satellite orbital pole alignments with the intermediate axis of the shear field tensor. 13 Welker et al. (2018) study alignments of satellite galaxies with filaments in their neighborhood in the Horizon-AGN simulation.
For our purposes here, the interesting alignments involve KPP orbital poles with directions characterizing large-scale mass flows converging toward their host, specifically, scales large enough to include high-redshift protosatellites before they are bound to the host galaxy and the mass surrounding them responsible for their torquing as well.These kind of alignments have not yet been fully addressed through hydrodynamic simulations and are the key to unveiling the physical origin of KPPs.Indeed, satellites were initially formed in connection with galaxies, hence following the same mass flows responsible for the formation of the latter (i.e., flows where satellites emerge essentially as the nodes of mass distributions at smaller scales).As mentioned above, flattened and elongated structures on ever larger scales appear at particular locations as the CW develops.Protosatellites emerging within these regions are expected to follow the mass flows causing them, traveling long distances before they reach the halo, and suffering (to different extents) the effects of forces and torques coming from the CW (proto)elements as they develop.These effects would give rise to different types of satellite orbital pole alignments with the directions of the main compression of flows and possibly to orbital pole clustering.This pole clustering would explain the formation of KPPs.
In this paper, this idea is explored by making use of cosmological simulations.Specifically, in this paper, we try to advance the quest for an answer to the following questions involving the kinematic structuring of satellites and their timescales: why and when was the clustering of KPP satellite orbital poles established, why are not all of the satellites involved in this clustering, and which role did the local CW development play therein.
Regarding the analysis of kinematic organization, our methodology here will be based on previous works focusing on the study of local CW developments around forming galaxies (Hidding et al. 2014;Robles et al. 2015).We aim at characterizing the local skeleton emergence by studying the shape deformation around galaxy-to-be objects, quantifying the timescales of deformation and the possible changes of the orientation of their principal directions and their freezing-out timescales.This information will be used along with the orbital angular momentum information of those satellites orbiting at the present time around these galaxies.By using the information about the kinematically coherent persistent planes formed around these systems (see Paper III), we try to clarify whether the satellites' anisotropic distribution forming satellite planes and the principal directions along which mass flows at larger scales are connected.
The paper is organized as follows.The simulations analyzed and their corresponding satellite samples are introduced in Sections 2 and 3, respectively.Section 4 is devoted to KPP satellite properties.The analysis of the mass density evolution around galaxy-to-be objects is addressed in Section 5 by introducing the Lagrangian volumes (LVs) in order to characterize the local CW.Section 6 reports on the specific properties of the LV evolution around the two galaxy systems analyzed in this paper.In Section 7, we study alignments of satellites (either individual poles or planes of kinematically coherent satellites) relative to the LVʼs principal directions across time.Results are discussed in Section 8. Finally, in the last section, we summarize our work and expose the conclusions reached.

Simulations
The simulated satellite samples analyzed in this paper are the same as studied in Papers II, and III, where the conditions to be met by galactic systems are discussed; the codes, runs, and satellite identification methods are described; and the relevant references are given.To guide the reader, a brief summary follows.
We study the planes of satellites orbiting around isolated, simulated host galaxies selected so that the following requirements are met: (a) the host galaxy at redshift z ∼ 0 is endowed with an extended (R ∼ 15 kpc) thin stellar and gaseous disk, (b) the assembly history of this central galaxy is free of major-merger events after halo virialization, (c) the system hosts a large (∼30) satellite population around the host, and (d) the simulation is run with a resolution high enough to ensure that the analysis of angular momentum conservation is made with sufficient accuracy.Thus, we require satellite objects to include more than 50 baryonic particles.
We have preanalyzed a set of different zoom-in cosmological hydrodynamical simulations, finding among them two that reach the previous prerequisites, the so-called "Aquarius-C α " resimulated halo and PDEVA-5004.The two simulations 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.
It is worth remembering that a standard two-phase process characterizes the halo mass growth: first a fast phase, where mass growth is largely provided by frequent merger activity, and then a slow phase, where growth rates and dynamical/ merging activity are low.The halo virialization time, T vir , roughly marks the separation between both phases.In the slow phase, the system formed by the halo and its bound satellites evolves independently from cosmic expansion.

Aq-C α
The initial conditions of this simulation come from the Aquarius Project (Springel et al. 2008), a selection of DMO MWsized halos formed in a ΛCDM simulation run in a 100 h −1 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. 2005Scannapieco et al. , 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).
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 9 M e ) 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 lowredshift event, this capture is beyond the scope of this paper.

PDEVA-5004
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.

Satellite Identification
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 accreted14 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.15By 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).
Satellites in Aq-C α and PDEVA-5004 span baryonic mass ranges of M bar = 8.5 × 10 6 −9.9 × 10 8 and M bar = 3.9 × 10 7 -1.8 × 10 8 M e , 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.

Orbital Properties
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,  J orb .Specifically, the top panels show the specific  J orb vector moduli or magnitudes (i.e., normalized by their baryonic mass), hereafter sJ orb , for the Aq-C α system; Note.Gas particles whose temperature is higher than Templim [K] are taken as pressurized hot gas particles and have not been considered here to calculate LV deformations.See text for parameter definitions.
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  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.

Orbital Pole Clustering Timescales
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,  J stack , were determined through the so-called "Scanning of Stacked Orbital Poles Method."16One (two) such axes have been found Figure 1.The components of the orbital angular momentum vector,  J orb , of each satellite as a function of time for satellites in the Aq-C α simulation, relative to axes that are kept fixed along cosmic evolution.The top panels show the evolution of the specific  J orb vector moduli or magnitudes (hereafter sJ orb ) from high redshift to z = 0.The second and third rows show the  J orb directions, represented by the polar (θ) and azimuthal (f) angles in spherical coordinates.Their evolution is given for T ta,AqC = 4 Gyr T uni 10 Gyr, an interval beginning at halo turnaround and centered at the halo virialization time T uni ∼ 7 Gyr (marked with vertical lines).The satellite samples are divided according to belonging (left panels) or not (right panels) to kinematically coherent, persistent planes (see Table 1, second column in Paper III).Different colors and line types stand for the satellite IDs, as coded in the sidebar on the legends.
As an illustration of these results, in Figure 2, we show the time evolution of ) , 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 ) angle is smaller than α co-orbit = 36°.87. Figure 2 indicates that a fraction of satellites are aligned to their corresponding  J 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  V vector), can be defined for the ith satellite as the time when the ) 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  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  J stack axes reported in Table 2, T cluster,Jstack entry.
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.

Evolution of the Kinematic Morphological Parameter
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 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,f 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  J stack axes.Specifically, we aim at studying the evolution of the κ rot parameter of KPPs as 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.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 ) 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, rot , allowed by Figure 3, points toward T max, rot ~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,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.
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., inplane 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.

Setting Up LVs: Method
To study the local CW evolution, we measure through TOI analysis the deformations of the LVs around the two host galaxies-to-be.We define these regions by marking the particles that at high redshift are within a spherical volume around the center of the galaxy progenitor and follow them forward in time (see Robles et al. 2015).
To build up LVs, the first step is halo selection at z = 0. Their virial radii r vir,z=0 are given in Table 1.Next, for each halo at z = 0, we have traced back to high redshift z high (see values in Table 1) all the particles inside the sphere defined by its respective r vir, z=0 .Using the position of these particles at z high , we have calculated a new center of mass  r c .Then, we have selected at z high all the particles enclosed by a sphere of radius R LV = K × r vir, z=0 /(1 + z high ), with K = 10, 15, and 20 (the motivation for this choice is discussed in the next subsection), around their respective centers  r c (see first row of Figure 4).
These particles sample mass flows shaping the CW elements as the Universe evolves.Particles follow geodesic trajectories until they possibly get stuck and begin the formation of, or are accreted onto, a CW structure element (i.e., a caustic).Our interest here focuses on the global deformations of LVs as an average of the fate of their constituent particles.Thus, we have followed these particles until z = 0; i.e., we have followed the evolution of the LVs from z high until z = 0. Note that, by construction, the mass of an LV is constant across evolution, as is the number of particles it is made of.
The choice of initially spherically distributed sets of particles aims to unveil the anisotropic nature of the local cosmological evolution, illustrated in Figure 4, where the LV corresponding to the Aq-C α LV at z high and its corresponding deformations until its final shapes and orientation at z = 0 are displayed.In this figure, we note that the LV has evolved into a highly irregular, anisotropic, and multiscale mass organization, including very dense subregions as well as other much less dense and even rarefied ones, with an overall flat structure from z ; 1 onward, corresponding to the formation of a large-scale sheet of the CW.It is very remarkable that filaments become coplanar, largely embedded into the flattening structure.
To quantify the local LV transformations illustrated in Figure 4, we have calculated, at different redshifts, the reduced inertia tensor, I ij r , of each LV relative to its center of mass: where r n is the distance of the nth LV particle to the LV center of mass and N is the total number of such particles.We note that the summation does not include hot gas particles, as their shapes are mostly driven by hydrodynamical/thermal pressure forces; see Robles et al. (2015).We have used the reduced tensor instead of the non-reduced tensor (Porciani et al. 2002b) to minimize the effect of substructure in the outer part of the LV (Gerhard 1983;Bailin & Steinmetz 2005).In addition, the reduced inertia tensor is invariant under LV mass rearrangements in radial directions relative to the LV center of mass; that is, characterizations of the LV shape would not be affected by these mass flows, hence making the I ij r tensor particularly suited to describe anisotropic mass deformations as those predicted by the ZA and the AM and observed in Figure 4.
In order to measure the LV shape evolution, first, we have calculated the principal axes of the inertia ellipsoid, a, b, and c, derived from the eigenvalues (λ i , with λ 1 λ 2 λ 3 ) of the I ij r tensor, so that a b c (see González-García & van Albada 2005;Robles et al. 2015), where M is the total mass of a given LV.Note that λ 1 + λ 2 + λ 3 = 2M, and this implies a 2 + b 2 + c 2 = 1.We denote the directions of the principal axes of inertia by e i ˆ, i = 1, 2, 3, where e 1 ˆcorresponds to the major axis, e 2 ˆto the intermediate axis, and e 3 ˆto the minor axis.The deformation of these LVs is conveniently described by the triaxiality parameter, T (Franx et al. 1991), defined as where T = 0 corresponds to an oblate spheroid and T = 1 to a prolate spheroid.An object with axis ratio c/a > 0.9 has a nearly spheroidal shape, while one with c/a < 0.9 and T < 0.3 has an oblate triaxial shape.On the other hand, an object with c/a < 0.9 and T > 0.7 has a prolate triaxial shape (González-García et al. 2009).

LV Properties
As a simple characterization of the local CW dynamics around forming galaxy systems, we study the global evolution of LVs from their initial spherical shape to the structures they span at low redshift.Different episodes stand out: (i) orientation of principal directions and freezing-out timescales (local CW skeleton emergence), (ii) LV deformations, and (iii) characterization of the times when deformation slows down.They will be analyzed in turn.As mentioned in Section 1, the local CW dynamics (or LV deformations) are mainly driven by the surrounding mass density evolution.Thus, the scale the LV tracks must be big enough so that it not only traces back the whole host galaxy plus satellites system but also represents the local CW dynamics around them.We have tested the robustness of our results against changes in this scale by comparing our results using different LV radii R LV = K × r vir, z=0 /(1 + z high ), with K = 10, 15 (fiducial value), and 20 (see Appendix).The results of the respective analyses have been compared.Sizes and masses of the different LVs, taking into account their scale, are given in Table 1.

Evolution of the Principal Directions
As said above, taking the LV as a whole, the I ij r eigenvectors, e z 1 ˆ( ), e z 2 ˆ( ), and e z 3 ˆ( ), 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 dir,e i .Should this happen, mass rearrangements at LV scales after T dir,e i would be organized in terms of a "skeleton" or fixed preferred directions, with the 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 e z i ˆ( ) and e z ), 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 1 = 0, then A 2 = A 3 .
We see that, when K = 15, A i (z) change at high redshift.A 3 (z) smoothly vanishes by T uni ∼ 2 Gyr ≡T dir,e 3 .Results for PDEVA-5004 with K = 15 are shown in Table 2. Changes in e 3 ând the other principal directions become unimportant very early, except for a small change of 5% at most in e 1 ˆand e 2 ôccurring 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.

Eigenvalue and Principal Axes: Shape Evolution
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).
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 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 shape,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 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.

LV Shape Evolution from the Perspective of CW Structure Formation and Dissolution
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 e 3 ˆprincipal direction (hereafter the e 3 ˆ-structure or, in this case, the 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 e 1 ˆprincipal direction, where most mass now piles up at the expense of the wall-like structure population.
We note that the 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 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 e 3 ˆ-structure (see bottom panel in Figure 6).

LV Alignments with Satellite Planes
In this section, we analyze alignments of satellites (either individual orbital poles or planes of kinematically persistent satellites) relative to the LV's principal directions across time from z high to z = 0.The robustness of our results against K changes is assured by the previous discussion on scale effects in Section 6.

Alignments with the Orbital Poles of Individual Satellites
We first consider how the orbital poles of individual satellites are oriented relative to the LV principal directions, e 1 ˆ, e 2 ˆ, and 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 e 1 ˆ, e 2 ˆ, and 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 e 1 ˆ, e 2 ˆ, and e 3 ˆalignments, respectively, while the shaded areas correspond to the respective 25th and 75th percentiles.
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, 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 e 1 ˆaxis.Indeed, the median of  J e cos , orb 1 ( ˆ) is ;0.9 after T uni ; 2 Gyr.Therefore, after that, KPP1 satellites tend to orbit on planes close to normal to the e 1 direction.On the other hand, satellite members of the KPP2 plane tend to have their poles aligned with the 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 e 1 ˆ, e 2 ˆ, or 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  J e , i orb a( ˆ) angles (Figure 7), with i = 3 for the Aq-C α and PDEVA-5004 KPP2 systems and i = 1 for the PDEVA-5004 KPP1 reach a value below α co-orb = 36°.87.Results are given in Table 2, entry T align,e i .We note the high dispersion of this value for the Aq-C α KPP.T align,e i also gives a measure of the timescale for the clustering of orbital poles.Results for T align,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  J 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 e 1 ˆand e 3 ˆprincipal directions are eight (all of them aligned with the  J stack axis defining the KPP1) and nine (including the seven satellite members of the KPP2 group aligned with the  J 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  J stack axis and 11 aligned with the e 3 ˆeigendirection.The number of satellites aligned with both of them at the same time is nine.This gives the conditioned fractions F(  J stack | e 3 ˆ) = 0.82, F(e 3 ˆ|  J stack ) = 0.69, and F(e i ˆ| no  J stack ) = 0.10.Thus, we have found a relatively high (low) probability of a KPP satellite having its pole aligned with the  J stack axis, in case it is (is not) aligned with the e 3 ˆor 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.

Global Alignments with KPP Planes
To characterize KPP orientations relative to the principal directions, either their normals,  n KPP , or the respective axes of maximum co-orbitation,  J stack , can be used.These are not equivalent analyses, as the normals  n KPP are returned by the TOI analyses of the KPP satellite positions, while  J 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,  n KPP and  J stack are not perfectly parallel.Information on the alignments among the principal directions of LVs and co-orbitation axes,  J 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  J stack axis forms a small angle with the 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  J stack axis for the PDEVA KPP1 plane, in turn, is close to the e 1 ˆdirection.Information on KPP orientations relative to the LV principal directions can be obtained from their respective normal vectors as positional planes,  n KPP .Figure 8 shows the angles between  n 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, ˆ(e 3 ˆ) principal directions.For the KPP identified in Aq-C α , its normal vector aligns with 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  n KPP forms with the 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  J 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,  n KPP ) are considered.

A CW Wall-like Structure in Aq-C α
As suggested by Figure 7, the 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 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 "e 3 ˆ-plane" (in contrast to "e 3 ˆ-structure") when referring to the mathematical plane that is perpendicular to the e 3 ˆdirection and contains the center of mass of the LV. Figure 9 shows some snapshots of the joint CW and KPP (proto)satellite17 evolution in physical coordinates in a reference frame with the z-axis aligned with e 3 ˆ, the y-axis aligned with e 2 ˆ, and the x-axis aligned with 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 e 3 ˆ-wall in the Aq-C α simulation, with the orbital pole alignment results shown in Figure 7.
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 e 3 ˆ-wall under different circumstances; some flow through filaments along with mass elements that eventually feed into the formation of the e 3 ˆ-structure, others are part of a small-scale wall that merges with the e 3 ˆ-structure in a parallel manner, and finally, others reach the 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 e 3 ˆ-structure (e.g., distance to the 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 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 e 3 ˆ-plane at higher velocities, more perpendicularly to it, and they tend to maintain a velocity component normal to the 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 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 e 3 ˆprincipal direction, feeding the 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 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 "e 3 ˆ-structure collapse."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 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 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.
Additionally, by following (proto)satellite spatial trajectories, we have found that many non-KPP satellites trace mass flows that are mostly perpendicular to the e 3 ˆ-plane, converging onto the plane and feeding it, while KPP satellites' trajectories tend to run more obliquely relative to the e 3 ˆ-plane.As seen in Figure 7, it takes some time to have the orbital poles of KPP satellites aligned with the e 3 ˆdirections (indeed,  T align, e 3 4.5 Gyr; see Section 7.1).To reach the aligned configuration that keeps roughly stable in time, a necessary condition is that the e 3 ˆprincipal direction is frozen.This happens very early (at T dir,e 3 ~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 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 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 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 e 3 ˆ-plane but with a larger impact parameter, such that their orbital poles are parallel to the plane and aligned with the 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 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 e 3 direction improve somewhat, but to a much lesser extent than during the previous period of 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).

PDEVA-5004: Two KPPs within a Prolate CW Structure
In the case of PDEVA-5004, the right panel of Figure 10 indicates that the median distances to the 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 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 e 3 ˆprincipal direction and have had their respective orbital poles aligned since very high redshift with the 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 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 e 1 Figure direction, and this clustering is maintained across e 3 ˆ-structure collapse.
As for non-KPP satellites, similarly to the Aq-C α case, most of them reach the e 3 ˆ-plane with a low impact parameter relative to the host center or have been recently captured by the host.

Timescales
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 e 3 ˆprincipal vector, corresponding to the direction of the dominant compression flow of matter, gets its direction fixed, T dir,e 3 (see Figure 5); and ii) the rapid high-redshift decrement of the minor principal axis, c(t), stops, T shape,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  J stack of maximum satellite co-orbitation determining the KPP plane the satellite belongs to, T cluster,Jstack (see Figure 2); ˆprincipal direction (two upper panels) or along the e 2 ˆprincipal direction at z = 0 (two bottom panels).A zoom of the respective central regions is also shown.Trajectories go from z high to z = 0 in comoving (box) coordinates with the host galaxy center as the origin of the reference system.The dots in the nonzoomed panels indicate the position of the protosatellite (see ID in legend) at z high , given in Table 1.Small crosses on the trajectories mark satellite positions at T dist,plane (i.e., the collapse timescale for KPP satellites onto the e 3 ˆ-structure).
iv) the ith satellite becomes aligned with the e 3 ˆprincipal vector (for satellites in KPP or KPP2) or the e 1 êigenvector (satellites in KPP1); see Figure 7 and T align,e i entry in Table 2; and v) there is a broad minimum for the first time in the median vertical distances to the 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 e 3 principal direction, a timescale stands out in the two zoom-in simulations studied: T dist,plane , the collapse timescale for the 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 align,e 3 .This is expected, because an aligned satellite orbits within the e 3 ˆ-structure; thus, its distance to the e 3 ˆ-plane must be low.This table also indicates that T dist,plane is coeval to T shape,e 3 , a timescale marking -within the accuracy of their determinations-the Universe age when the initially strong mass flows normal to the 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 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 e 3 ˆ-plane; therefore, their clustering (defined by the 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,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 e t 1 ˆ( ) principal direction.

Summary
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 (  J 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 e t i ˆ( ), 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 e 3 ˆdoes not change (within a threshold) after its freezing out at T dir,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 e 3 ˆaxis after T align,e 3 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 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 e 3 ˆaxis after T align,e 3 ~3.5 Gyr is even clearer than in the Aq-C α case.KPP1 satellite poles have been well aligned with the 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 e 3 ˆpresent closer, more oblique trajectories relative to the 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 e 1 ˆaxis show trajectories parallel to e 3 ând 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.

Some Comparisons with Other Works
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 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 align,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 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. (2015Libeskind et al. ( , 2019) ) tested the possible alignments between the observed satellite planes in the local Universe and 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. 18While 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.

Conclusions
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., 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.

KPP member satellites characterized by pole alignments
with the 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 e 3 ˆ-structure when it collapses.The latter's trajectories at high redshift tend to be oblique relative to the 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 e 3 ˆ-structure and are closer to the 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 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 e 1 ˆdirection, and of some satellites in the Aq-C α simulation, with orbital poles aligned with the e 2 êigenvector.
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 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 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 largevolume simulations of galaxy formation might shed light on how frequently the different channels for satellite plane formation appear throughout cosmic evolution.
the 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 Center, Spain), the NLHPC (funded by ECM-02) and Ladgerda cluster (Fondecyt 12000703).This project has received funding from the European Union Horizon 2020 Research and Innovation Programme under Marie Skłodowska-Curie Actions (MSCA) grant agreement No. 101086388-LACEGAL.

Appendix Evolution of the Principal Directions
In this Appendix, we test how the principal directions of compression behave when we change the LV scale, R LV = K × r vir,z=0 /(1 + z high ).In the top panel of Figure A1, we see that the timescale for principal direction fixing when using a scale of K = 20 (involving a volume increase of a ∼2.4 factor compared to the fiducial K = 15) behaves identically as for K = 15, as none of the principal directions change after T uni ∼ 4 Gyr ≡ T freeze for either K = 15 or K = 20.The same is true in the case of PDEVA-5004 after T uni ∼ 2 Gyr ≡ T freeze , except for the very small changes around T uni ∼ 6 Gyr.
To deepen these results, we have also analyzed how the principal directions of LVs at different scales are oriented with respect to each other.For the Aq-C α simulation, these principal directions after T uni ∼ 4 Gyr are essentially the same for either scale; see Figure A1, bottom panel.For PDEVA-5004, this is also true from the very beginning of our analysis, more accurately for the e 3 ˆdirections.Finally, it is worth noting that, however, when using a K = 10 value, only one axis freezes out that early, while for the other two, it takes a longer time, as they become frozen within 10% by T uni ∼ 6 Gyr.This is due to the fact that a shorter scale tracing of the density field around the galaxy-to-be formation site gives an LV dominated by denser mass elements, where more abrupt/complex dynamic processes take place.Therefore, K = 10 LVs at high z are dominated, in both simulations, by early activity in their central regions and are thus not suited for our purposes here.

Figure 2 .
Figure2.Alignment between the axes of maximum satellite co-orbitation,  J stack , and the orbital poles of satellites belonging to the KPP in Aq-C α (upper panel) and the two persistent planes KPP1 and KPP2 in PDEVA-5004 (middle and lower panels, respectively).

Figure 3 .
Figure 3. Behavior of the v v i i , 2 f() ratios as a function of cosmic time for KPP and non-KPP satellites in the Aq-C α and PDEVA-5004 systems; see legends.Lines are the median values at each simulation output time, and the shaded bands give the corresponding 25th and 75th percentiles.Velocity components are taken relative to the corresponding  J stack axes.

Figure 4 .
Figure 4. Shape evolution of the reference Aq-C α LV from z high = 8.45 to z = 0 (six snapshots).Two LV projections along fixed axes are shown for each snapshot, one along the z = 0 e 3 ˆprincipal direction onto the X-Y plane (first and third columns) and the other along the z = 0 e 2 ˆprincipal direction onto the X-Z plane (second and fourth columns).The snapshot redshift z and Universe age are given in each panel.This figure shows how the initially spherical mass distribution flattens with time in such a way that by z ∼ 1, a wall-like structure has emerged, and by z ∼ 0.5, mass piles up in a predominant filament within the X-Y plane.

Figure 5 .
Figure5.Evolution of the cosine of the angle A i formed by the eigenvectors e i (z) and e i (z = 0) with i = 1,2,3 for the Aq-C α simulation for the LV with a K = 15 scale.The upper horizontal axes give the redshift scale, while the lower ones stand for the Universe age T uni .

Figure 6 .
Figure 6.Properties of the LVs of the Aq-C α (upper panel) and PDEVA-5004 (middle panel) systems.Arrows mark the respective T shape,e3 timescales, when the rapid high-redshift changes of the c(t) principal axes slow down.The top and middle panels show the evolution of the principal axis lengths across time and their respective derivatives.The bottom panel gives the evolution of the axis ratios b/a and c/a for both simulations, with the Universe age color-coded according to the color bar.The c/a vs. b/a plane is split into three regions according to the values the T shape parameter takes on them; see the blue, green, and orange isotriaxiality curves.

Figure 7 .
Figure7.Alignment between the  J orb of satellites with respect to the principal directions of the LVs in Aq-C α (upper block panels) and PDEVA-5004 (lower block panels).Alignments are shown for satellites belonging to KPPs and outside these structures as well; see legends.Thin lines correspond to individual satellites, with colors and line types as encoded on the right of the panels.Thick lines are the medians at each T uni , and shaded bands mark the 25th-75th percentile range.

Figure 8 .
Figure 8. Global alignments of the principal directions with  J stack and the normal vectors to the KPP planes,  n KPP , along cosmic evolution.The upper panel shows the results for the Aq-C α simulation, while the middle and lower panels show the results for PDEVA-5004 KPP1 and KPP2, respectively.Different line types stand for the cosine of the different angles.Thick magenta lines:  J stack and  n KPP .Thin solid lines:  J stack and the principal directions.Thin dashed lines: n KPP and the principal directions.Cyan, orange, and purple colors stand for each of the corresponding principal directions, as usual (see legends).

Figure 9 .
Figure 9. Projections along the e 2 ˆprincipal direction of the density field around the host galaxy formation site for different values of the Universe age given in the top left corner of each panel (Aq-C α simulation in physical coordinates).Red points sample the (proto)satellite mass elements.Evolution goes from left to right panels, and the rightmost panel refers to the host halo virialization time.
Figure11shows some illustrative examples of these two different satellite trajectory behaviors.

10.
Median distances of the different satellite populations in the Aq-C α (left panel) and PDEVA-5004 (right panel) simulations to the corresponding plane normal to e 3 ˆthat contains the center of mass of the LV.Shaded areas stand for the corresponding 25th and 75th percentiles.

Figure 11 .
Figure 11.Trajectories of four satellites belonging to KPP (leftmost columns) and non-KPP (rightmost columns) populations in the Aq-C α system.Two projections are shown for each satellite, either along the e 3ˆprincipal direction (two upper panels) or along the e 2 ˆprincipal direction at z = 0 (two bottom panels).A zoom of the respective central regions is also shown.Trajectories go from z high to z = 0 in comoving (box) coordinates with the host galaxy center as the origin of the reference system.The dots in the nonzoomed panels indicate the position of the protosatellite (see ID in legend) at z high , given in Table1.Small crosses on the trajectories mark satellite positions at T dist,plane (i.e., the collapse timescale for KPP satellites onto the e 3 ˆ-structure).

Figure A1 .
Figure A1.Top panel: evolution of the cosine of the angle A i formed by the eigenvectors e i (z) and e i (z = 0) with i = 1, 2, 3 for the Aq-C α simulation for an LV K = 20 scale, that is, R LV =20 • r vir,z=0 /(1 + z high ).Bottom panel: relative orientations of the principal directions of LVs tracing the K = 15 and K = 20 scales.Upper horizontal axes give the redshift scale, while the lower ones stand for the Universe age T uni .

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

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