Galactoseismology in cosmological simulations Vertical perturbations by dark matter, satellite galaxies, and gas (cid:63)

Context. Complex models recently became available for studying the dynamics of disk galaxies such as the Milky Way (MW). These models include the global dynamics from dwarf satellite galaxies, dark matter halo structure, gas infall, and stellar disks in a cosmological context. Aims. We use a MW model from a suite of high-resolution hydrodynamical cosmological simulations named GARROTXA to establish the relationship between the vertical disturbances seen in its galactic disk and multiple perturbations from the dark matter halo, satellites, and gas. Methods. We calculated the bending modes in the galactic disk in the last 6 Gyr of evolution. We computed the vertical acceleration exerted by dark matter and gas in order to quantify the impact of these components on the disk, and compared this with the bending behavior with Fourier analysis. Results. We ﬁnd complex bending patterns at di ﬀ erent radii and times, such as an inner retrograde mode with high frequency and an outer slower retrograde mode excited at di ﬀ erent times. The amplitudes of these bending modes are highest during the early stages of formation of the thin disk (20kms − 1 ) and reach up to 8.5kms − 1 in the late disk evolution. We ﬁnd that the infall of satellite galaxies leads to a tilt of the disk, and produces strong anisotropic gas accretion with a misalignment of 8 ◦ with subsequent star formation events and supernovae, creating signiﬁcant vertical accelerations on the disk plane. The misalignment between the disk and the inner stellar and dark matter triaxial structure, which formed during the ancient assembly of the galaxy, also leads to a strong vertical acceleration of the stars. We also ﬁnd dark matter subhalos that temporally coincide with the appearance of bending waves in certain periods. Conclusions. We conclude that several agents trigger the bending of the stellar disk and its phase spirals in this simulation, including satellite galaxies, dark subhalos, misaligned gaseous structures, and the inner dark matter proﬁle. These phenomena coexist and inﬂuence each other, sometimes making it challenging to establish direct causality.


Introduction
The stars of the Milky Way's (MW) disk appear to be vertically disturbed.These perturbations can be detected in the stellar density distribution of structures like the warp, corrugations, and north-south asymmetries (Widrow et al. 2012;Williams et al. 2013;Schönrich & Dehnen 2018).More recently, using the highquality data of stars in the Solar Neighbourhood from the Gaia satellite (Gaia Collaboration et al. 2018;Brown et al. 2021), Antoja et al. (2018) found a new structure in the phase-space: the phase spiral.This structure is the result of a perturbation to the disk (external or internal) which vertically displaces stars and triggers a phase mixing process.To explain its origin, the focus has been put on the Sagittarius dwarf galaxy, which is a known recent perturber of the Milky Way (Binney & Schönrich 2018;Laporte et al. 2019;Bland-Hawthorn et al. 2019;Li & Shen 2020;Bland-Hawthorn & Tepper-García 2021;Hunt et al. 2021;Gandhi et al. 2021;Widmark et al. 2021).Recent observations have shown that the formation mechanism of the phase spiral is more complex than initially expected (Hunt et al. 2022;Antoja et al. 2022;Alinder et al. 2023;Frankel et al. 2023).These structures in the vertical phase space are products of the perturbations of the stellar disk, which can be seen as deformations in both the kinematics and positions with respect to the galactic plane.These deformations can also be detected as bending waves.This bending behavior has also been detected in other galaxies (Gómez et al. 2021;Nandakumar et al. 2022;Urrejola-Mora et al. 2022).
Several theoretical studies using N-body and analytical models have attributed vertical perturbations to interactions with satellites, (e.g.Gómez et al. 2013;Debattista 2014;de la Vega et al. 2015;D'Onghia et al. 2016;Laporte et al. 2018;Semczuk et al. 2020;Bennett & Bovy 2021), the dark matter wake as a consequence of the interaction between the halo and an external perturber (Weinberg 1998;Gómez et al. 2016;Laporte et al. 2018;Grand et al. 2022), and the interaction with dark matter sub-halos (Feldmann & Spolyar 2015;Chequers et al. 2018;Tremaine et al. 2023).Other amplification mechanisms of the vertical perturbations of the disk include gas infalls along warps (Khachaturyants et al. 2022) and misalignment of the plane of stars at different radii (Sellwood & Debattista 2022).Thus, vertically disturbed stars can have a diverse range of origins and a perturbed disk may be a common state throughout galactic disks.
Grion Filho et al. (2021) used a purely gravitational N-body simulation to calculate the accelerations produced by the disk, the Sagittarius-like dwarf galaxy, and dark matter halo to discern the dominant source of perturbations in the different parts of the disk and different moments of evolution.They found that there are moments when the dark matter halo's response to the passing of the satellite enhances its effect.However, to have a complete picture of the dynamical mechanisms taking place in the Galaxy, one should include the hydrodynamical processes.Cosmological hydrodynamical simulations take into account realistic satellite accretion into the host galaxy and hence multiple satellites can perturb the disk, whereas isolated N-body models only account for one perturbing satellite interacting with the disk.However, multiple satellites can lead to different effects when interacting with the disk and with each other (Rocha et al. 2012;Widrow et al. 2014;D'Onghia et al. 2016;Trelles et al. 2022).In Gómez et al. (2017) they used cosmological simulations to study vertical patterns in galactic disks.They found that approximately 70% of the simulated galaxies show clear vertical patterns produced by close satellite encounters, distant flybys of massive companions, accretion of misaligned cold gas, and re-accretion of cold gas from the progenitors of gas-rich major mergers.However, the spatial resolution of the simulation does not allow to study scales of the size of the solar neighborhood.The study of van de Voort et al. (2015) used a cosmological model to follow the evolution of highly misaligned gas, where the inner parts of the galaxy re-align with it faster, resulting in a warp at later times of evolution.Semczuk et al. (2020) used galactic disks from the Illus-trisTNG cosmological simulation to study tidally induced warps, finding that both satellite interaction and gas accretion may induce warps.
In our previous work (García-Conde et al. 2022; GC22 hereafter), we found phase spirals in a Milky Way-like disk for the first time in a cosmological simulation, namely the GARROTXA (Roca-Fàbrega et al. 2016) model.However, its exact origin was not clear, since the satellites in this model are less massive (∼ 10 9 M ⊙ at the moment of the pericenter) than the minimum mass reported in the bibliography to excite phase spirals, at least in the context of the Sagittarius dwarf galaxy and the Milky Way.For example, the models by Laporte et al. (2019) and Bland-Hawthorn & Tepper-García (2021) used a mass of 6 × 10 10 and 2 × 10 10 M ⊙ respectively.We suggested that other effects such as the interaction of the collective effects of several satellites, or the misalignment of the angular momentum of the different components of the simulation, especially the stellar disk and dark matter halo, or the combination of some of them can generate strong perturbations in the disk.Later on, Grand et al. (2022) Table 1.Main characteristics of the MW-sized model G.323 of GAR-ROTXA simulation at z= 0. From left to right, virial radius, virial mass, the total stellar particles mass, mass of the disk, the typical mass of the disk particles, and number of disk particles.

R vir
154 6.5 × 10 11 6 × 10 10 2.6 × 10 10 10 3 − 10 4 2 × 10 6 used a high-resolution simulation of a Milky Way-like galaxy from the Auriga suite of cosmological simulations and resolved phase spirals in the disk.The main formation mechanism of the phase spiral in that simulation was found to be the dark matter halo response to the passage of a satellite galaxy.Differently from our case, they have a dominant source of perturbation, so they can identify the main cause of the triggering of the phase spiral.
Here, we further analyze the vertical behavior of the disk of the GARROTXA model.In this work, we identify bending modes (as calculated in Widrow et al. 2014) and we study their behavior using Fourier decomposition.We calculate the vertical acceleration of the dark matter particles onto the disk and its relation with the passing satellites and the inner structure of the halo itself.We also analyze the vertical gravitational acceleration applied by the gas to understand the relationship between its misalignment with respect to the disk.Thus, we present a more realistic picture of the satellites, halos, gas, and disk interactions.
In this paper, we start in Sec. 2 by describing the simulation and its different structures.In Sec. 3 we describe the behavior of the disk and quantify its disturbances, focusing on the vertical component.In Sec. 4 we calculate the vertical accelerations (a Z ) applied by dark matter and gas onto the disk and identify all satellites in the system.In Sec. 5 we link the measured bending modes and disk tilting to the different typology of interactions occurring in our model, characterized in the previous sections, and we find many intertwined agents involved in the bending of the disk.In the last section, we discuss the results and we present our conclusions.

The GARROTXA Simulation
In this work, we use model G.323 of the GARROTXA simulation (Roca-Fàbrega et al. 2016), which is similar (although at the lower end of the predicted virial mass range) to the Milky Way in terms of disk size and rotation curve.The spatial resolution of this model is 109 pc, with a minimum dark matter particle mass of 10 5 M ⊙ , a typical mass of disk stellar particles of 3.5×10 3 M ⊙ , and a minimum time-step of 10 3 yr.At z ∼ 0, we have ∼ 2 × 10 6 stellar particles within the virial radius.The summary of these characteristics is shown in Table 1.
The analysis presented here focuses on data encompassing the last 6.3 Gyr (from 6.3 Gyr in lookback time to the present, equivalent to scale factors of a = 0.6 to a = 1.0, respectively) of the evolution of the simulated galactic system.This is the period where we find a well-defined stellar disk.In this period, we have 190 snapshots available.This allows us to undertake an exhaustive temporal evolution study of the disk's kinematics and density structures.To read and analyze the simulation we use the software yt (Turk et al. 2011).
The main dark matter structures of this model at the mentioned period are a halo with a virial radius R vir = 154 kpc1 with a total mass of 6.5 ×10 11 M ⊙ .The inner parts of the dark matter halo are dominated by a central triaxial substructure with an ellipsoidal shape extending to approximately 10 kpc in its major axis and 5 kpc in its minor axis.This structure was formed at z∼ 3 (11.6Gyr in lookback time) as a result of a 1:1 major merger (Roca-Fàbrega et al. 2016).The oldest stellar particles, both belonging to the primordial galaxies (pre-merger) and born during this major merger, show properties similar to an old elliptical galaxy embedded in a tri-axial dark matter halo, which is not rotationally supported.After this 1:1 major merger, a new stellar disk starts to form but is perturbed by another merger, which is the last major merger of the galactic system studied here.This last major merger occurs at z∼1.5 (∼1:10 merger at 9.55 Gyr of lookback time) and sets the more recent components of the structure of the galaxy that is a thick-old stellar disk and a much younger thin stellar disk.
At each snapshot, we apply a centering and alignment of the main galactic system.The center of mass is calculated iteratively, in spheres of shrinking radii, as in Power et al. (2003).After this centering, we compute the disk plane by calculating the angular momentum of stars inside a 20% of R vir , and we then apply a change of basis to align this with the Z axis.This places most of the galactic disk in a plane oriented with Z = 0.The X-axis is the projection of the original X-axis from the cosmological box onto the galactic plane whose normal vector is L.Then, we take all stellar particles born after z∼ 1.5 inside a fixed cylindrical section with Z < |3| and R < 7 kpc and apply the same process which further aligns the inner galaxy's L with the Z axis.We use galactocentric cylindrical coordinates (R, ϕ, Z) where ϕ and V ϕ are negative in the direction of rotation.The direction of ϕ = 0 is oriented towards the positive X axis.
In Fig. 1 we show the different structures of the stellar component of the model at the present time (z = 0) in its edge-on and face-on projections (top and central rows, respectively) and the face-on view colored by age (bottom panels).We separate stellar particles depending on whether they were born after (left column) or before the last major merger mentioned above (middle and right columns).In practice, we take a time of 9.05 Gyr (i.e.500 Myr after the merger) to ensure that we are selecting the particles born in the thin-young stellar disk.
Regarding the stellar particles born before the merger, they are distributed in the above-mentioned ellipsoid (right column) akin to the one in dark matter, and a structure that today composes an old, dynamically hot disk (central column).To separate this old disk from the ellipsoid, we calculate the angle α between the vector of Z = 0 plane and the angular momentum of the stars, considering stars with cos(α) < 0.7 as part of the ellipsoid.
For our study, we use only particles born after the last major merger, thus belonging to the 'young' disk (left panels) to track any signs of vertical perturbations.In the inner galactic regions (<5 kpc), there is an oval structure whose major axis is oriented at the same angle at all times and has an A 2 /A 0 average amplitude of 0.35.This elongated structure has a quadrupole signal in the velocity V R , but since we do not detect a pattern speed, we do not consider it to be a bar.This structure contains stellar particles of the same age as the disk.The semi-major axis of this young elongated structure is oriented as the old-stellar and dark matter inner triaxial structure's minor axis.There is also a one-armed structure, or lopsided mode, in the young disk stellar component.We come back to this aspect in Appendix B, although for this study we focus mostly on the vertical disturbances.Finally, for the following analyses, we further select stellar particles with |Z| < 2.5 kpc and R < 25 kpc, similar to GC22.In summary, in this study, we select particles from the disk based on their age (born after 9 Gyr in lookback time) and their location (aforementioned region).After this selection, we end up with ∼ 1.39 × 10 6 stellar particles in the disk.

Vertical bending modes of the disk
In this section, we focus on the Z and V Z distributions of the disk through its bending modes, similarly to Widrow et al. (2014), as a way to quantify the disk's vertical deformation.These modes can lead to phase spirals (Bland-Hawthorn & Tepper-García 2021; Darling & Widrow 2019), and therefore our additional motivation is to be able to explain those found in GC22.By identifying the timing and spatial location of the perturbations, we hope to establish a link with the underlying phenomenon, whether internal or external, that cause this deformation.
To calculate the vertical modes, we divide the X-Y plane into 70 bins from -25 to 25 kpc in both coordinates X and Y (bin size of ≈ 0.7 kpc).This binning allows us to have enough stellar particles (an average of 350 particles per bin) to perform the following calculations.In each bin, we calculate the linear regression of V Z as a function of Z.The slope and the intercept correspond to the amplitudes of the breathing (A) and bending (B) amplitude respectively: An example of this method can be seen in Fig. 2, where we show the face-on projection of the galactic disk at two different lookback times of evolution corresponding to 2.33 Gyr (left columns) and 0.99 Gyr (right columns).We show an overdensity plot of the disk particles with a Gaussian filter (top panels), the bending (middle panels) and breathing (bottom panels) amplitudes of the young disk computed with Eq. 1.In the top panels, we can see a one-armed density structure, which is leading and retrograde.In the central panels, we show how strong bending waves appear at specific times while they can be weaker in others (right vs left panels).In the breathing modes (bottom panels), we see a prominent mode m = 2, which follows the oval structure at the inner parts of the stellar disk, and it is present all through the recent galaxy evolution.This is consistent with previous studies, where this breathing behavior is linked to the bar or other non-axisymmetric internal structures (Widrow et al. 2014;Faure et al. 2014;Monari et al. 2015).
Here we mostly focus on the bending modes, while the breathing modes and their relation to the density patterns (i.e. the central oval structure and the retrograde one-armed pattern) are briefly examined in Appendix B for the sake of completeness.
To visualize the evolution of the bending parameters over time, we take the values of B in the X-Y grid and apply a Fourier analysis that we evaluate in a polar grid, dividing the disk into radial sections with 1 kpc of width.For each ring, we have N bins from the initial X-Y grid, each of which is located in an az-imuthal position ϕ i , and has a bending amplitude (weight in the Fourier decomposition) of W i = B i .Then, we perform a Fourier decomposition following the next expression2 : We estimated the amplitudes for each mode, A m , at each radius and time.We focus on the results for m = 1 and m = 2 in the rest of the analysis.We also analyzed higher modes, but, even though at some times they are comparable with the amplitudes of modes m = 1 and m = 2, they are never dominant.We also study the phase at which the maximum amplitude for each mode is located at each radius, following the expression: In Fig. 3 we show the amplitude and phases of the bending modes m = 1 (left panels) and m = 2 (right panels) as a function of radius and time.The bending pattern seems to be very complex.Over the timespan studied, at R > 7, there is an overall dominance of red-black colors (corresponding to a phase of 0-100 • ), indicating that the bending is prone to be excited at this angle.At the beginning of the evolution (between 6.3 and 5 Gyr), there is a highly perturbed disk, as indicated by the higher amplitudes of the bending modes (top row of panels), which can be up to 20 km s −1 at external radii.At these times, the phase of this bending mode m = 1 (bottom left panel) at radii R > 7 shows a slow retrograde pattern ( φ > 0 in our reference system).From 2 Gyr to the present, we see that the amplitude of the bending mode m = 1 increases, with the phase oriented as in previous times, reaching its maximum at approximately 1.1 Gyr, with a value of 7.5 km s −1 between 10 and 15 kpc, and 8.5 km s −1 between 15 and 20 kpc with a leading shape.
At inner radii (< 7 kpc), we generally see a fast retrograde pattern ( φ > 0, in the bottom left panel) that is excited at certain times (top left panel).We observe that its amplitude is maximum at 3.5 Gyr (with 3.5 km s −1 of amplitude).
The m = 2 mode of the bending shows high amplitudes at earlier times (top left panel), and generally at higher radii, has prograde motion ( φ < 0).From 2.5 Gyr to 1.8 Gyr we see an increase at intermediate regions of the disk, although the bending mode m = 1 generally dominates at these times and radii (we further discuss the behavior of the bending modes amplitudes in Sec. 6).

Perturbation agents: Dark matter, satellite galaxies and gas
To quantify the impact that the dark matter and gas have on the disk's vertical direction, we have calculated the vertical gravitational accelerations a Z they produce on the galactic disk which is on the Z = 0 plane at every snapshot (see Sec. 2).This is computationally faster than calculating accelerations on all disk particles.This mesh has 100 × 100 bins in the X and Y directions and extends to R = 40 kpc.We have used the following expression: where the vertical acceleration at each bin (a z,bin ) is the sum of all individual vertical accelerations from all the mass elements N part (dark matter particles or gas cell) of each component (dark matter and gas).m i is the mass of each mass element, and r i,bin is the total distance between the bin and the mass element.The relation r z;i,bin /r i,bin takes the vertical component of the acceleration.Some mass elements can be very close to the Z = 0 plane, producing large accelerations, and increasing the noise in the Fourier analysis.To mitigate this, we establish a limit of maximum acceleration in a bin, corresponding to the acceleration an element can apply if it is situated at a vertical distance of ±54 pc to the plane, which is half of the minimum size of the gas cell (equivalent to the spatial resolution of the simulation: 109 pc).If an acceleration in a bin exceeds this limit, it is set to this maximum.
In Fig. 4 we show a Z in the Z = 0 plane at two different time periods (three leftmost panels vs. the three rightmost) coinciding with two instances when there are at least one of the satellites near its pericentre.In the top row, we show the acceleration by the dark matter particles, where we can see the effects of the satellite, crossing the plane between the first panel (negative Z, negative vertical accelerations) and the next one (positive Z, positive vertical accelerations).The three rightmost panels show another interaction, but in this case, the satellite is not crossing the plane at the moment of minimum distance.However, we can see the vertical accelerations it induces in the outer parts of the grid (upper right regions in red), which are outside of the disk.
In the bottom row of Fig. 4, we show the acceleration exerted by the gas.We can see that it is stronger at the disk region, where it is more concentrated.At certain times (e.g. at 4.09 Gyr), the map appears to be divided into one half of the gaseous component inducing positive acceleration and the other, a negative one, suggesting the possibility of a tilted or warped gaseous disk.This will be further analyzed in Sec.4.3.
Once the calculation of a Z is carried out, we apply the same Fourier decomposition as in previous sections.This allows us to better visualize the evolution of the vertical accelerations with time, understand the effects that these have on the disk dynamics, and compare them with the results from the previous sections.In the next subsections, we present the accelerations from the different components.

Accelerations by Dark Matter
Here we first take the approach of considering the impact of the dark matter as a whole, without excluding any particles belonging to satellites.This will include the effects of the global structure of the halo, its satellites, and associated unbound material, and the possible wakes excited by these satellites.The results are shown in Fig. 5, where we display the amplitude of the m = 1 and m = 2 modes and their respective phases of maximum amplitude, ϕ m .At inner radii (< 10 kpc, and mostly in the m = 1 mode), the particles that are likely to dominate the acceleration are the ones from the gravity of the inner structure of the dark matter halo by proximity, over the particles situated at higher Fig. 4. Vertical acceleration of dark matter (top panels) and gas (bottom panels).Black dot-dashed circles are situated at 20 kpc, illustrating the size of the stellar disk of the galaxy.The dark matter acceleration at the top three leftmost panels is dominated by the crossing of dark matter sub-halos DM1119 and DM2595.The effect of gas at these times (bottom panels) is bimodal.The dark matter's a Z of the right block of panels is dominated by the satellite galaxy "Escarabajo", which applies positive a Z in the positive x -positive y region of the mesh.radius, like the ones from the satellites.Indeed, in Sec. 5 we show that the inner dark matter structure is triaxial and slightly tilted in relation to that of the disk, most of the time, producing a constant acceleration.This tilt is especially relevant at certain moments, such as 1.5 Gyr approximately, which coincides with the moments of higher accelerations at these radii.
In the outer parts of the disk, the acceleration produced by the dark matter component (> 20 kpc) comes mainly from the dark matter of satellites orbiting the system, especially when they cross the galactic plane.In the following Sec.4.2 we study in detail the isolated effects of the satellites.

Accelerations by satellite dwarf galaxies and dark sub-halos
The GARROTXA model is a zoom-in cosmological simulation, thus, the halo of the main galactic system contains many substructures both in stars and dark matter.Here we used the Rockstar halo finder (Behroozi et al. 2012) to identify all sub-halos that may impact the disk's kinematics.With this halo finder, we identify all dark matter structures with M 200 > 10 6 M ⊙ (to avoid poorly resolved sub-halos and contamination by field halo dark matter particles), and, using the dark matter particle IDs we also find all gravitationally bound particles.If the sub-halos have stellar mass within 20% of the R 200 of the satellite (as calculated by Rockstar) we identify them as luminous.We use the value of M 200 right before the infall into the R Vir of the host galaxy as a reference of the size and mass of the sub-halo.Once the satellites have entered the virial radius of the host, it is highly difficult to determine their mass.In these cases, we use or indicate the bound mass given by Rockstar.
The complete data of the sub-halos (including both dark and luminous satellites) found by Rockstar can be found in Fig. C.1, in the Appendix C. By applying this software, we have identified new dark satellites and some dwarf galaxies compared to GC22.In the timespan studied, we found nearly 500 independent subhalos (luminous or dark) that enter the galaxy.As expected from the ΛCDM Universe, we find a large number of dark satellites, i.e. sub-halos with no stellar particles.For instance, inside the virial radius of the studied galaxy at redshift z = 0, we find 139 of the aforementioned sub-halos (those of lower mass get destroyed over time) without any stellar particle, compared to 10 dwarf galaxies.These dark satellites, although they are generally less massive than dwarf galaxies, can interact with the disk as luminous satellites would.
In Fig. 6 we show again the top left panel of Fig. 5, i.e., the amplitude of m = 1 mode of the dark matter vertical acceleration, but now showing also the total, not just vertical, acceleration from each sub-halo on the center of the galactic disk.The total acceleration by the most important dark satellites is shown with dashed lines, whereas the ones containing stellar particles are shown with solid lines.The radius and time at which each satellite crosses the plane are marked by a cross of the same color.We observe that the increase in the vertical acceleration in the regions from 25 to 40 kpc (background colors) is consis-tent with the close encounter with sub-halos both dark and luminous (lines and crosses).The constant disk crossing of the many dark and luminous satellites is permanently perturbing the disk.However, only those of higher mass, and those that get closer to the galactic plane show significant a Z (see background colors vs. lines).For instance, at 5.1 Gyr or 2.9 Gyr the dark satellites DM2512 and DM1747 (violet and yellow lines) cross the disk but we do not observe a strong increase in a Z (background colors).
In contrast, the strongest interactions that we observe are the following.First, a strong interaction occurs at 4.1 Gyr, and it is due to a dark satellite (DM1119, blue dashed curve) containing a mass of 1.1 × 10 9 M ⊙ .This interaction is fast, and the pericentre occurs when the object is crossing the plane, at 20 kpc.This satellite is accompanied by another dark satellite, DM2595, with a mass of 10 8 M ⊙ (orange), and the dwarf galaxy Arania's second pericenter, whose mass at this moment is of 1.3 × 10 9 M ⊙ (blue line with stars).
Secondly, another noticeable moment is the interaction occurring at 2.2 Gyr with a dwarf galaxy which we will refer to as "Escarabajo" (light green).This satellite starts its interaction with 1.6 × 10 9 M ⊙ , being a slower interaction than the previous example at 4.1 Gyr.In this case, the pericenter does not occur near the galactic plane yet we see a band of strong accelerations in the background colors.The black dashed line is the sum of the accelerations of those satellites whose mass does not exceed 10 8 M ⊙ or maximum acceleration does not exceed 2 km s −1 Myr s −1 .This shows the accelerations of different dark matter sub-haloes that, even though they are less massive, may have a relation to some of the increases seen in the a Z applied by the dark matter particles.However, the increase in acceleration by dark matter is not always related to passing satellites.We find that at 1.1 Gyr, the increase is not entirely explained by the passing of minor satellites, as in the case of Fig. 4. By examining the acceleration by dark matter in the X-Y projection, we find that this is related to the overall dark matter distribution.As will be discussed in Sec.5.2, at these times there are instabilities of the galaxy that reflect as a net acceleration onto the plane.Finally, from around 1 Gyr until the present we observe several sub-halos, which overall increase the accelerations in the intermediate and outer regions, like the pink and light green ones with pericenters occurring at 0.4 and 0.26 Gyr respectively.

Acceleration by gas
We apply the same method as in Sec.4.1 to the gas cells of the simulation, whose results are shown in Fig. 7.In the bottom left panel, which shows the phase of the mode m = 1, we see a sharp transition between the outer parts and the inner parts of the disk.This transition is located between 5 and 10 kpc during the first snapshots and between 10 to 15 kpc from 5 Gyr on.This transition can also be seen as a sort of break in the amplitude.
At the inner regions, between 4.5 and 3.5 Gyr, we observe a retrograde rotating mode m = 1 that extends to 10 kpc approximately.This pattern has a period between 300 and 500 Myr and is similar to the one of the inner bending wave seen in the left panels of Fig. 3. Indeed, we see that at these times and radius, the gaseous disk can be misaligned with the stars by about 4 degrees (top panel of Fig. A.1), thus causing this specific mode m = 1 in acceleration.
At the outer parts (at radius of R > 15 kpc), we observe low amplitudes of the acceleration's mode m = 1 (top left panel).At these radii, there is a rather constant phase (bottom left panel), which corresponds to an approximate angle of -90 • (blue colors).In this same panel, from 2 Gyr on, this phase changes to red/black (corresponding to angles between 50 • and 100 • ), propagating towards the outer radius.This is accompanied by a higher amplitude (top left panel).From 1 to 0 Gyr the amplitude of m = 1 decreases, and the phase shows similar angles to previous instances.To study this in more detail we examined the gas distribution and its temperature in our model.
In Fig. 8 we show the gas in a zoom-out, face-on, projection (left column), a zoom-out, edge-on projection (central column), and a zoom-in, edge-on projection of the gaseous disk (right column) colored by temperature.In the first row (5.61 Gyr) we see a highly clumpy distribution of cold (blue) gas, stripped from Arania at its first encounter (at >6 Gyr).The right panel of this row (zoom-in) shows a deformed gaseous disk, highly tilted from the galactic plane.
Over time, as seen in the next zoom-out face-on panels rows, the gas falls into the galaxy preferably through a cold (blue) gas filament in the negative Y and positive X direction of -50 • approximately at radii of > 50 kpc, and -90 • approximately at smaller radii throughout almost all snapshots studied.This is especially well seen at the panels of 2.88, 2.56 and 1.88 Gyr (first column of the third, fourth and fifth row respectively, where the gas falls from the left parts of the plane, dominated by blue colors.Indeed, the latter direction of -90 • is compatible with the observed phase of the acceleration in the Fourier mode m = 1.This gas is a mixture of the one that has been stripped from satellites (Arania's at early times, mostly tidally in its first pericenter, and ram-pressure stripped from Escarabajo at later times) and the one that outflowed from the disk at an early time through supernovae winds and initially kept in the circumgalactic medium.This metal-rich gas has shorter cooling times and thus falls into the galaxy (Dekel & Birnboim 2006).
In particular, at 2.88 Gyr Escarabajo is already inside the virial radius and, as mentioned above, has suffered ram-pressure stripping due to the presence of a warm-hot gas corona around the galaxy.Gas stripped from Escarabajo mixes with the pre-existing inflowing gas (as seen at 2.88 Gyr and 2.56 Gyr in the third and fourth rows, respectively), cools down and, over the next half gigayear, enters the disk and increases the star formation (see figure 5 in GC22).Following this star formation event, from 2 Gyr on, we observe a change in the phase of the m = 1 mode (in Fig. 7), along with a higher amplitude of m = 1 in a Z which propagates from the inner parts towards the outer radii.At these times, the gaseous disk is tilted with respect to the stellar disk due to the recent interaction with the satellite Escarabajo.After examining the snapshots, we observe that there is also an anisotropic ejection of gas from the supernovae winds following the aforementioned star formation event in this direction which might be further perturbing the disk in this direction.As a consequence of these different phenomena, at 1.88 Gyr (right panel of 5th row in Fig. 8) the gaseous disk is tilted by about 4 • to 8 • (as seen in Appendix A, and in Sec. 5).Along the last gigayear of evolution, star formation is regulated by the galactic fountain leading to new supernovae events (although with weaker mass ejection than in the former case).

Triggering agents of bending deformation
In this section, we make a connection between the times of higher vertical acceleration (from both dark matter and gas) and the behavior of the disk in terms of bending modes calculated in previous sections.In Fig. 9 we take the data from Figs. 3, 5,  and 7, and get the mean amplitude of the Fourier modes within a radial bin (see top labels) to compare between them at specific radial regions.In panel A we show the mean amplitude of the m = 1 Fourier mode of the vertical acceleration produced by dark matter, between 20 and 40 kpc.In this region, the effects of satellites and dark matter sub-halos are dominating, as seen in Sec.4.1, and we see clear peaks at their moments of maximum approach.Panels B and F show the mode m = 1 of the acceleration by gas (orange line) and dark matter (blue line) at radii from 10 to 20 kpc and 0 to 10 kpc, respectively.We observe that the amplitude of the m = 1 mode of a Z applied by the gas is as significant as the one by the dark matter and even more dominant in the region from 10 to 20 kpc, due to its strong misalignment with the disk.At radii from 0 to 10 kpc, the m = 1 acceleration by the gas is comparable to that exerted by the dark matter halo.However, in terms of total vertical acceleration, the inner dark matter triaxial structure dominates due to the contribution of higher modes (Fig. 5).
In addition to this, we calculate the angle between the angular momentum vector of the stellar disk and that of the gas measured in spherical shells with radii between 10 to 20 kpc and between 0 and 10 kpc (dashed orange line in panels C and G, respectively).We observe that there is a correlation between the acceleration applied by gas and its angular misalignment with the disk.As mentioned in Sec.4.1, here we also obtained the tilt of the central dark triaxial structure with respect to the stellar disk by calculating its maximum height and corresponding radius, with the former being the Fourier amplitude of the vertical distribution of dark matter particles in a cylinder with |Z| < 10 kpc.Although this tilt is only of 0.5-1.25 • at times (blue lines in panels C and G), given the total mass of the triaxial inner dark matter structure this is enough to create large accelerations like the ones observed in Fig. 5.We also measure the changes in the angular momentum direction of the stellar disk in relation to the previous snapshot taking the cosmological box as a reference (panels D and H).Panels E and I show the mean amplitudes of the bending modes m = 1 (green lines) and m = 2 (purple lines), Fig. 7. Fourier decomposition of the acceleration of the gas cells (similar to Fig. 5) onto the Z = 0 plane (galactic plane).We show m = 1 (left panels) and m = 2, from which we have the amplitude of the mode (top panels) and their corresponding phase of the maximum amplitude (bottom panels).
for regions between 10 and 20 kpc and between 0 and 10 kpc, respectively.
The perturbed initial state of the galaxy and the successive interactions with satellites make it very difficult to pinpoint a moment in which the galaxy is at equilibrium along its history.However, in Sec. 4, we have described the key moments where interactions with satellites and dark matter sub-halos occur.Based on these, we divide our analysis into two intervals.The first is defined by the interaction with the dark satellites (at 4.1 Gyr), and the other is that of the dwarf satellite galaxy Escarabajo (at 2.2 Gyr).We study both interactions in more detail and their possible consequences on the disk.In Appendix D we briefly study how these variables correlate to each other more quantitatively.A more descriptive explanation based on that and Fig. 9 is presented below.

Multiple dark and luminous satellites and bending waves
At 4.1 Gyr, in Fig. 3, we can see that a relatively higher amplitude of bending mode m = 1 is excited in the inner (< 5 kpc) parts of the disk.In the external parts (> 10 kpc) the increase starts later (from 3.9 Gyr on approximately).Both increases are in relation to previous instances where the amplitude is lower. .In the case of outer radii, we observe a very slow bending mode m = 1 and prograde bending mode m = 2, which coincide with the pericenters of two of the dark satellites (DM1119 and DM2595) at 4.1 Gyr.These have, respectively, virial masses M 200 of 10 9 and 10 8 M ⊙ , plane-crossing positions at cylindrical galactocentric radii of 20 kpc and 15 kpc, and velocities V Z of 364 km s −1 and 386 km s −1 , respectively.Their gravitational acceleration can be seen in the first three top panels in Fig. 4, and in Fig. 6.The pericenters of these satellites are compatible with the increase in bending modes at outer radii.During this time, the second pericentric passage of the satellite Arania also takes place at 3.9 Gyr (blue line in Fig. 6), which at the time of this in-teraction has a bound mass M bound of 1.1 × 10 9 M ⊙ , as calculated by Rockstar (5×10 8 M ⊙ as calculated in GC22, by measuring the tidal radius at apocenters of the orbit), which may also reinforce the bending patterns.
The inner retrograde bending pattern does not have a clear origin.It seems to be present throughout all the timespan studied.The oscillation of the inner bending mode (panel I of Fig. 9) could be related to the inner density structures of the disk (see Appendix B).But, specifically from 4 Gyr, it has a higher amplitude with a maximum at 3.5 Gyr (top left panel of Fig. 3) and a more clear phase (bottom left panel of Fig. 3).In Sec.4.3 we found that the accelerations from the gas show a wave-like pattern with similar behavior in rotation and periodicity.Its amplitude (panel F) increases approximately 400 Myr before the interaction with the dark satellites DM1119 and DM2595 (blue and orange line respectively, at 4.1 Gyr), and other less massive sub-halos near the disk at that moment.Thus, it is unlikely that these sub-halos are the cause of the increase in inner bending.At these times, we also see a slight misalignment in the gas (panel G).Our hypothesis is that gas is being accreted and destabilizing both the stellar and gaseous disk, inducing these specific bending waves, or at least enhancing the already existing one.

Distant interaction triggering bending modes via disk misalignment
The dwarf satellite Escarabajo has a total mass of M 200 of 1.6 × 10 9 M ⊙ , and a stellar mass of 6 × 10 6 M ⊙ at its first infall, before its closest approach (2.2 Gyr).This interaction is slower (V Z = −50 km s −1 at the moment of the pericenter) than the ones mentioned in the previous subsection (see Sec. 5.1).Escarabajo's pericentric passage occurs at R = 40 kpc and at 16 kpc above the plane, thus far from it.We saw in Fig. 6 that the vertical acceleration by this satellite is mostly at the outer disk and barely reaches the inner parts.Even though the amplitude of the bending mode  m = 1 (panel E of Fig. 9) reaches its maximum a Gyr later, its increase starts approximately 300 Myr after the pericenter of Escarabajo, which is half of the rotation period at these outer parts (600 Myr at 20 kpc).Furthermore, we need to take into account that, this being a slower interaction, the acceleration is not only applied at the moment of the pericenter but also several Myr before and after its maximum approach.Thus, this initial increase is compatible with the passing of Escarabajo.
The angular momentum of the stars also starts to change at 300 Myr after Escarabajo's passing (as seen in panels D and H).This indicates that Escarabajo's passage pulls and tilts the disk, especially at outer radii, where self-gravity is smaller.It also slightly shifts from its alignment with the central triaxial structure of the dark matter halo, which also generates accelerations that prevail longer in time.This creates a misalignment between the disk and the central dark matter and old stars (panel G) triaxial structure, which also exerts an acceleration into this disk.This is reflected in the amplitude of mode m = 1 of the acceleration by dark matter (panels B and F, blue lines), which reaches its maximum at approximately 1.2 Gyr.
Furthermore, in Sec.4.3, we have seen how Escarabajo loses its gas and this mixes with previous existing inflows.We measure the cold gas inside Escarabajo's virial radius before its maximum approach and we find that this satellite contributes to 10 7 M ⊙ of cold gas (T < 10 4 K).At these times, the misalignment of the gas disk in relation to the plane is 8 degrees (panel C of Fig. 9) and has an associated increase in the vertical acceleration by the gas at radii of 10-20 kpc (panel B).This is a period of destabilization of both the gaseous and stellar disk as a consequence of the direct interaction with Escarabajo.Even though the infalling gas may not induce large perturbations by itself due to its mass, this new gas is settling in a different angle in respect to the disk.This new gas forms stars and the subsequent SNe bubbles further perturbs the disk.
The bending mode m = 1 finally reaches its maximum at approximately 1.1 Gyr.This perturbation is present across stellar populations as a whole, but we find that it is stronger in stellar particles from 4 to 6 Gyr of age, approximately.The phase of this bending mode (bottom left panel of Fig. 3) indicates that it is aligned with the major axis of the central triaxial dark matter structure, which at this moment is tilted in the opposite direction.Thus, the central triaxial structure of dark matter halo appears to be crucial to the behavior of the bending pattern once it has been triggered.In Sec. 6, we also discuss briefly the possible effects of the disk itself on the shape of the bending wave.
Being such a complex system implies the existence of interactions that are amplified (or dampened) by the different components and structures.Thus, the gas and especially the central triaxial dark matter structure seem to contribute to a general destabilization between different parts of the disk, leading to the amplification and prolongation in time of the bending mode.It is expected to have a delay in the large-scale effects caused by the "tug" of satellites if they induce other imbalances.In addition to this, there are also all the new passages of other, less massive, but non-negligible, satellites from 1 Gyr on (see Fig. 6).

Summary
The existence of vertical waves and vertical phase mixing in the MW and in other galaxies has motivated our exploration of the bending modes in the stellar disk of the GARROTXA cosmological simulation.Using Fourier decomposition, we analyzed these modes and the vertical gravitational accelerations produced by the dark matter halo, gas, and satellite galaxies.Our main findings on the study of the bending of the stellar disk are: 1.The bending behavior of the stellar disk is much more complex than previously seen in N-Body simulations with different regimes at different radii and times.Throughout the simulation, we find different m = 1 modes coexisting but being exited at different times: One is slower and extends from ∼7 kpc to 20 kpc.This external mode is generally aligned with the major axis of the inner structure of old stars and dark matter.We also find an inner, with much lower amplitude, retrograde mode with high frequency dominating at <5 kpc. 2. The amplitudes of the bending mode are high in the early stages of the thin disk formation (20 km s −1 ) and they can be up to 8.5 km s −1 in the late disk evolution.3. The slow pericentric passage of a satellite galaxy of mass M 200 of 1.6 × 10 9 M ⊙ triggers a disk tilting.When the stellar disk is tilted, it becomes misaligned with the inner dark matter triaxial structure (formed in an ancient major merger) which is less responsive to these perturbations as it is more massive.The effect on the disk bending is magnified by this new misalignment with the dark matter inner structure, leading to the highest vertical deformation in the last several Gyrs of evolution.4. When the satellite enters the virial radius of the host, the gas of the satellite is also stripped by ram-pressure and it slowly cools down and falls into the disk along with preexisting inflows.This cold gas enters the stellar disk nonisotropically, settles in a slightly tilted new gaseous disk, and, lately, induces new star formation events and the creation of supernovae bubbles.Even though these processes may have a weaker impact on the disk, they are compatible with the accelerations seen exerted by the gaseous disk.Thus, all these processes might magnify the initial bend.5.After one of the several interactions with dark sub-halos, the inner m = 1 retrograde bending mode increases its amplitude and shows a phase pattern that correlates with the one of the gas vertical acceleration.Although this amplitude is relatively lower than the external bending discussed above, this may suggest that the gas-stellar disk misalignment plays an important role in the observed inner disk bending.

Discussion
We thus identified the following agents that trigger the bending of the stellar disk: 1. Dark sub-halos: Dark satellites with masses ∼ 10 9 M ⊙ interact with the stellar disk on many occasions coinciding with a higher amplitude of the bending mode m = 1.The presence of such dark sub-halos is expected in the ΛCDM galaxy formation theories and their interaction with the disk has been mostly studied as a source of disk heating (see e.g.Benson et al. 2004;D'Onghia et al. 2016).Here we highlight the role of such structures as triggers of bending perturbations of the disk.2. Satellite galaxies: Satellite galaxies in the GARROTXA simulation at z> 2 have at most M 200 ∼10 10 M ⊙ before their first infall, and most of them do not reach inner regions of the galaxy.However, although their direct tidal influence in the inner stellar disk is not large, they lead to a cascade of events that trigger some perturbations as discussed below.3. Inner density profile: The inner tri-axial structure of the dark matter is an imprint of its ancient accretion history.At certain times, the largest vertical gravitational acceleration comes from a misalignment between the stellar disk and this dark matter ellipsoid.The misalignment is triggered by the tilting of the stellar disk, and the resulting bending can last for more than 1 Gyr.Similar phenomena have also been observed in Sellwood & Debattista (2022), where slow bending waves arise from misalignment between the inner and outer regions of the disk.We cannot discard an intrinsic coupling mode between the disk and halo contributing to this deformation.Furthermore, it is possible that the overall inner density profile of the disk also exerts a torque after the initial disk tilting, resulting in the leading shape of the warp (Shen & Sellwood 2006;Gómez et al. 2016), which is worth exploring in future analysis.4. Gas accretion and misaligned gaseous structures: The vertical accelerations from the gaseous structures, often misaligned with the stellar disk, are comparable to those from satellites and dark matter.The origin of this misalignment is the inhomogeneity of the cold flows, the slightly different angular momentum of the incoming gas with respect to the actual stellar disk, and the expansion of supernovae bubbles.
We speculate that all these phenomena may have a further impact on the bending behavior of the stellar disk, although we cannot measure in which exact proportion.The effect of extra-planar non-homogeneous distributed gas has also been detected as a source of vertical patterns and disk tilting in other models (e.g.Gómez et al. 2017;Khachaturyants et al. 2022).
Below we discuss other perturbation mechanisms that we have not looked into in detail but may need special consideration: 5. Wake: We have not found a dominant global asymmetry like the one expected when a strong wake is present in the Fourier analysis of the GARROTXA model.This could be because there is only one satellite galaxy with comparable mass to the Sagittarius Dwarf Galaxy, and its pericenter occurs at large radii from the main host's center.This satellite has a mass above 10 10 M ⊙ before entering the virial radius of the host for the first time, and it loses mass down to 4.7×10 9 M ⊙ at the moment of the first pericenter.Since the orbit of the satellite is mostly polar and slightly retrograde, we expect most of the torques to be in the vertical direction.However, in our case, both the collective and transient response to the wake would be weaker than expected for a Sagittarius-like satellite.Even if we have not quantified the exact torque applied by the wake, its effects are already captured by our calculation of accelerations by all dark matter particles.Furthermore, due to the dipolar nature of the halo's response, the wake effects are accounted for in our Fourier analysis.6. Streams: We also have studied the tidal streams from Arania after its first pericenter (both stellar and dark matter debris) and we have found that the maximum vertical acceleration applied by this material is negligible in comparison with the progenitor and the rest of the dark matter halo.However, in the case of more massive satellites, their debris could be significant enough to affect the general structure of the dark matter halo (Garavito-Camargo et al. 2021).7. Density structures of the disk: Other structures present in the GARROTXA model are a central stellar oval and a onearmed retrograde density pattern (see Appendix B).We have found some weak correlations between these stellar disk's non-axisymmetric structures and the breathing mode.Thus, we suspect that the physics of the vertical behavior in our model is also related to the density structure of the disc.However, we will study how these self-gravitating structures can perturb the disk in the future (see e.g.Grion Filho et al. 2021 for an analysis in a pure N-Body simulation).
In GC22 we detected an increase in the amplitude of the phase spirals in intermediate-age stars coinciding with the most important satellites' percienters.We observe that the periodicity of the thick bands present in Figure 5 of GC22 his of approximately 1.2 Gyr, which is roughly coincident with the period of the one-armed density mode (Appendix B).This suggests that the phase spiral may propagate as this density mode in our model.The observed phase spirals are related to the phase mixing of the bending modes of the disc studied here, even if these appear to be not very large compared to other models (amplitudes of 8.5 km s −1 in this model compared to amplitudes up to 30 km s −1 in Gómez et al. 2017).In Tremaine et al. (2023) they propose a mixed scenario with several large-scale perturbations and numerous small-scale perturbations.This is qualitatively compatible with our model in the sense of having multiple dark-(and not dark) satellites affecting the disk.However, due to the resolution limitations, and the several interactions occurring simultaneously, we could not accurately track the spiral pattern or measure the winding of the phase spiral and, except for some temporal coincidences, we cannot attribute them to a specific event.The conclusions of this work point also to a mix of many processes (not present in the model of Tremaine et al. 2023) such as the structure of dark matter halo, accretion history, gas content, etc.
We conclude that although further study is needed, the phase spirals we observed can be a consequence of the many perturbations that the stellar disk suffers directly or indirectly related to the satellite galaxies' infall but also gas accelerations and inner dark matter structure.

Caveats
There are also some caveats in our study that need to be considered.For our analysis, we define a new reference system at each timestep.However, it is not trivial to define a midplane with such a complex vertical structure, with different bending at different radii or for different populations.Variations in the determination of the galactic plane can lead to inaccuracies in the measurements of dynamical properties (Beane et al. 2019).
The application of the tools used in this work is necessary to detect and analyze the variety of processes taking place in cosmological simulations.They are also necessary in order to analyze a larger set of models, which is crucial for a more global understanding of disk dynamics in real galaxies.However, another caveat in our methodology is that our Fourier analysis does not always provide conclusive evidence to establish true causality.This likely comes from two aspects.Firstly, when there are many processes in place, and some of them are complexly interconnected, it is not straightforward to distinguish between real agents of perturbation or simple coexisting effects, or even to separate causes and their effects.Secondly, this limitation might also be related to the Fourier analysis.Localized accelerations like those from sub-halos and satellite galaxies will be somehow smoothed by the sine/cosine decomposition, especially when we only keep the series to the first orders (m = 1 and m = 2 modes).In the case of the internal non-axisymmetric dark matter structure, its total acceleration is not completely captured by the m = 1 mode, since it is also dominant in higher modes.Finally, since there are different disk regimes with radial limits that change with time, our analysis can be biased by exploring the disk in rings of fixed radius that might not be the natural ones.
To overcome these limitations, in the future, it might be worth exploring Basis Function Expansions or Multichannel Singular Spectrum Analysis for both the acceleration field and the disk kinematic perturbations.More sophisticated tools to link the different agents and effects might improve our correlation analysis (Appendix D).These techniques are starting to be ap-plied to disc dynamics in more controlled simulations (Weinberg & Petersen 2021;Garavito-Camargo et al. 2021;Johnson et al. 2023) and this might establish the foundations for its application in more complex cases such as the model used here.These expansions might be also useful for a more global study of the dynamics, not centered exclusively on the vertical or planar dimensions.
Also, even if the number of satellites in this model is compatible with the predictions for the Milky Way's mass (Bullock & Boylan-Kolchin 2017), a lack of interaction with massive satellites (M * >10 10 M ⊙ ) in this system that could be equivalent to the LMC and SMC makes it challenging to compare directly with the Milky Way at the present time.However, this exercise helps us to understand the evolution of disk galaxies in a global way, where past interactions can have a long-lasting impact on stellar disk morphology and kinematics.

Final conclusions and remarks
A general conclusion of our work is that even galaxies that are far from dense groups can be highly and complexly perturbed systems in a realistic cosmological context, never reaching dynamical equilibrium.We adopt the idea by Grion Filho et al. (2021) of the importance of a holistic approach to disc dynamics, to which we add here additional factors: that of the effects of the inner profiles of dark matter and old stars from previous major mergers, and the gas brought in by satellites and accreted through filaments.These factors are often overlooked by more simplistic models that therefore do not capture the full complexity of galactic systems.
It is undoubtedly a huge effort to bring cosmological simulations into the picture of disk dynamics.Although detecting and understanding the complex interplay between all existing phenomena and the different dynamical regimes in the disc in cosmological simulations is challenging, the analysis of coexisting perturbative processes will lead to more realistic galactoseismology.These figures serve as a visual aid to identify in which moments different variables increase or decrease in relation to each other, and, combined with Fig. 9, can help to better understand the causes and effects at play at the different time intervals discussed in Sec. 5.However, due to the small number of snapshots per window (average of 23 snapshots), we cannot extract accu-rate and definite conclusions based on these correlations on their own.
In the three figures we observe that the first gigayear of the studied period presents high correlations between all variables due to the highly perturbed state of the galaxy and the subsequent relaxation over the next gigayear.We have a special interest in the correlations with the bending m = 1 at external radius starting at 2 Gyr, since it is the higher bending mode we find over time and radius.In the bottom panel of Fig. D.1 we find a high correlation between this bending mode and the mode m = 1 of a Z applied by dark matter.We also find breathing modes being

Fig. 1 .
Fig. 1.Different structures of the stellar component of the snapshot at present times z = 0.The left panel shows stars born 500 Myr after the last major merger (at 9.55 Gyr in lookback time).The panels in the central column show the older, hot disk.Lastly, in the right column, we show the old stellar ellipsoid.The top two rows show the surface density of the stellar particles in the X-Z and X-Y projections, respectively.The last row shows the face-on projection colored by age.The disk rotates clockwise in this reference system.

Fig. 2 .
Fig. 2. Overdensity plot of the disk particles with a Gaussian filter σ = 6 (top panel), bending (middle), and breathing (bottom) maps at two different snapshots of this simulation in each column, 2.33 Gyr and 0.99 Gyr in lookback time.A strong bending mode can be seen at 0.99 Gyr (middle right panel).

Fig. 3 .
Fig. 3. Temporal evolution of the bending modes' Fourier amplitude A m (top panels) and phase ϕ m (bottom panels) of the m = 1 (left panels) and m = 2 (right panels) modes.The phase of mode m = 2 ranges between −90 • and 90 • due to the duplicity of the structures.The white horizontal grey-dashed lines mark the division of the disk at a radius of 10 and 20 kpc.The bending m = 1 dominates over m = 2 and shows higher amplitudes at outer radii (R> 7 kpc) excited mainly at the beginning of the timespan studied and at later times from 2 Gyr on.The inner parts show faster retrograde bending waves with the highest amplitude happening around 4 Gyr of lookback time.

Fig. 5 .
Fig.5.Fourier decomposition of the vertical acceleration by dark matter particles onto the Z = 0 plane.We show m = 1 (left panels) and m = 2, from which we have the amplitude of the mode (top panels) and the phase where the maximum amplitude is located (bottom panels).These are similar to 3, but extended to 40 kpc of radius.The inner acceleration is dominated by the inner dark matter structure.

Fig. 6 .
Fig.6.Total gravitational acceleration magnitude produced by each individual sub-halo on the center of the galactic system (colored lines, right y-axis) and the sum of acceleration by the rest of (less massive) sub-halos (black dashed line, right y-axis).Sub-halos containing stellar particles are indicated as solid lines, and those with pure dark matter are shown as dashed lines.We also marked the moments and positions (cylindrical radius, left axis) when the satellites ("x") cross the galactic plane.The grey "x" symbols correspond to the crossings of less massive satellites.Background color and left axis are equivalent to the top-left panel in Fig.5.Notice that, while the background color shows the vertical acceleration by all dark matter on the galactic plane Z = 0, lines show the total acceleration produced by sub-halos at the galactic center, not at the Z = 0 plane.The latter gives proportionally lower values of acceleration due to the difference of the distance between the satellite and the closest point of the plane vs to the center of the galaxy.For the lines, we also use the data of M bound for each satellite, re-evaluated at each apocenter.The masses indicated on the label are the first mass registered for each satellite, which corresponds to M 200 .In this figure, we see that the vertical acceleration from the dark matter halo between 25-40 kpc is correlated with the pericentric passages of several satellites, dark or luminous.The ones exerting more vertical acceleration are DM1119 and Escarabajo, as well as DM479 and DM2113 at later times.

Fig. 8 .
Fig. 8. Distribution and temperature of the gas in GARROTXA.Left and central columns: Zoom-out face-on (Y-X plane) and edge-on (Y-Z plane) projection of the gas temperature at six different instances of the simulation colored by temperature.The crosses and dashed lines indicate the position and direction of the satellites Arania and Escarabajo.Right column: Edge-on projection of the gaseous disk, colored by temperature.The dotted horizontal line represents the plane of the stellar disk.The main inflow is reaching the galaxy from the negative Y direction.

Fig. 9 .
Fig. 9. Top left panel (A) mean value of the amplitude of Fourier mode m = 1 of a Z by the dark matter (blue lines) at external radii (20-40 kpc).The vertical lines show the moments of maximum approach to the disk by those satellites of higher mass, and the labels indicate their M 200 at infall.The left panels show the calculated properties at radii from 10 to 20 kpc, whereas the right panels show between 0 to 10 kpc.Panels B and F show the mean value of the amplitude of mode m = 1 of the vertical acceleration of dark matter (blue lines) and gas (orange lines).Panels C and G show the angle difference between the angular momentum of stars and the angular momentum of gas (dotted orange lines, leftmost axis) and the tilt of the central dark matter triaxial structure with respect to the disk (dotted blue lines rightmost axis), denoted by θ gas and θ ellipsoid , respectively.Panels D and H show the angular momentum vector variation in relation to the previous snapshot divided by the timespan.Panels E and I show the mean amplitude of Fourier amplitudes for bending modes m = 1 (green lines) and m = 2 (purple lines).

Fig
Fig. B.1.Similar to Fig. 3 but for the stellar mass.

Fig. C. 1 .
Fig. C.1.Distance of satellites to the centre of the galaxy.The top panel shows dark satellites of 10 6 M ⊙ < M <10 7 M ⊙ .The second panel are satellites 10 7 M ⊙ < M <10 8 M ⊙ .Third one M > 10 8 M ⊙ .The last panel show all satellites with stellar mass.

Fig
Fig. D.1.Pearson's correlation coefficient in moving windows of 800 Myr of width for bending and breathing modes vs accelerations of gas and dark matter.We mark as circles the moments where Pearson's coefficient fulfils |r| > 0.4 and p − value < 0.05

Fig. D. 3 .
Fig. D.3.Same as Fig D.1 but with variables of density vs bending and breathing modes.
excited at these times, which also correlate with the acceleration by dark matter at these times.From 2.2 Gyr until 1.6 Gyr the density mode m = 1 correlates with the acceleration by dark matter and gas at more internal radii between 0 and 10 kpc (top panel of Fig. D.2, light blue and orange dots, respectively).At external radii (bottom panel of said figure) there is not a clear tendency.At these times and radii, in the top panel of Fig. D.3, we also see a correlation between the density mode m = 2 and the breathing mode m = 2 (dark green dots) and between the modes m = 1 of density and