# Shaping the Milky Way: The Interplay of Mergers and Cosmic Filaments

Arpit Arora¹,² ID, Nicolás Garavito-Camargo³,⁴ ID, Robyn E. Sanderson² ID, Martin D. Weinberg⁵ ID, Michael S. Petersen⁶ ID, Silvio Varela-Lavin⁷ ID, Facundo A. Gómez⁷ ID, Kathryn V. Johnston⁸ ID, Chervin F. P. Laporte⁹,¹⁰,¹¹ ID, Nora Shipp¹ ID, Jason A. S. Hunt¹² ID, Gurtina Besla⁴ ID, Elise Darragh-Ford¹³ ID, Nondh Panithanpaisal¹⁴,¹⁵ ID, and Kathryne J. Daniel⁴ ID

The EXP collaboration

¹Department of Astronomy, University of Washington, Seattle, WA 98195, USA; arora125@sas.upenn.edu

²Department of Physics & Astronomy, University of Pennsylvania, 209 S. 33rd Street, Philadelphia, PA 19104, USA

³Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

⁴Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

⁵Department of Astronomy, University of Massachusetts, Amherst, MA 01003-9305, USA

⁶Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

⁷Departamento de Física y Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena, Chile

⁸Department of Astronomy, Columbia University, 550 West 120th Street, New York, NY 10027, USA

⁹LIRA, Observatoire de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CY Cergy Paris Université, CNRS, 92190 Meudon, France

¹⁰Institut de Ciencies del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franqués 1, E-08028 Barcelona, Spain

¹¹Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

¹²School of Mathematics & Physics, University of Surrey, Stag Hill, Guildford, GU2 7XH, UK

¹³Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA

¹⁴Carnegie Observatories, 813 Santa Barbara St., Pasadena, CA 91101, USA

¹⁵TAPIR, California Institute of Technology, Pasadena, CA 91125, USA

Received 2025 April 28; revised 2025 June 3; accepted 2025 June 7; published 2025 July 24

## Abstract

The large-scale morphology of Milky Way (MW)—mass dark matter halos is shaped by two key processes: filamentary accretion from the cosmic web and interactions with massive satellites. Disentangling their contributions is essential for understanding galaxy evolution and constructing accurate mass models of the MW. We analyze the time-dependent structure of MW-mass halos from zoomed cosmological-hydrodynamical simulations by decomposing their mass distribution into spherical harmonic expansions. We find that the dipole and quadrupole moments dominate the gravitational power spectrum, encoding key information about the halo’s shape and its interaction with the cosmic environment. While the dipole reflects transient perturbations from infalling satellites and damps on dynamical timescales, the quadrupole—linked to the halo’s triaxiality—is a persistent feature. We show that the quadrupole’s orientation aligns with the largest filaments, imprinting a long-lived memory on the halo’s morphology even in its inner regions (~30 kpc). At the virial radius, the quadrupole distortion can reach 1–2 times the spherical density, highlighting the importance of environment in shaping MW-mass halos. Using multichannel singular spectrum analysis, we successfully disentangle the effects of satellite mergers and filamentary accretion on quadrupole. We find that, compared to isolated MW–LMC simulations that typically use a spherical halo, the LMC-mass satellite induces a quadrupolar response that is an order of magnitude larger in our cosmological halo. This highlights the need for models that incorporate the MW’s asymmetry and time evolution, with direct consequences for observable structures such as disk warps, the LMC-induced wake, and stellar tracers—particularly in the era of precision astrometry.

Unified Astronomy Thesaurus concepts: Galaxy dark matter halos (1880)

## 1. Introduction

The formation and evolution of galaxies are intimately tied to the morphology of their dark matter (DM) halos, which assemble hierarchically through anisotropic accretion along cosmic filaments and mergers with smaller halos (e.g., S. D. White & M. J. Rees 1978; G. R. Blumenthal et al. 1984; V. Springel et al. 2005). The surrounding large-scale environment dictates the angular momentum and accretion history of these halos, imprinting morphological and kinematic signatures that persist over cosmic time (D. J. Eisenstein & A. Loeb 1995; D. Aubert et al. 2004; O. Hahn et al. 2007; C. A. Vera-Ciro et al. 2011).

The shapes of Milky Way (MW)—mass halos, in particular, are complex and evolve in response to both internal and external processes. In their inner regions ~30 kpc, baryonic contraction and mergers tend to produce oblate or prolate shapes (J. Dubinski & R. Carlberg 1991; S. Kazantzidis et al. 2004; S. Bryan et al. 2013; V. P. Debattista et al. 2013; C. Vera-Ciro & A. Helmi 2013; J. Bovy et al. 2016). Beyond those radii, the same halos become markedly triaxial, retaining memory of their filamentary assembly and anisotropic infall (B. Allgood et al. 2006; C. A. Vera-Ciro et al. 2011; G. Despali et al. 2014; M. Valluri et al. 2021; J. Baptista et al. 2023). Past major mergers, tidal interactions, and ongoing cosmic accretion continuously distort the outer halo, breaking any perfect symmetry (e.g., J. Bailin & M. Steinmetz 2005; P. Bett et al. 2007; M. D. Schneider et al. 2012; M. Bonamigo et al. 2015; F. A. Gómez et al. 2017; A. Arora & L. L. Williams 2020; J. J. Han et al. 2022). A key question is whether—and how—we can recover a halo’s filamentary

<figure><img src="image_1.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
“memory” once an LMC-mass satellite has fallen in. Most MW-mass galaxies host or have hosted such massive companions (e.g., G. Besla et al. 2007; A. Robotham et al. 2012), so understanding the superposition of long-lived filamentary distortions and transient merger signatures is crucial—both for interpreting external analogs (e.g., SAGA; M. Geha et al. 2017) and for piecing together our Galaxy’s own history.

The MW provides a unique opportunity to study the shape and evolution of DM halos owing to the wealth of available stellar kinematic data and to the fact that the LMC’s orbit and mass (~10%–20% of the MW; G. Besla et al. 2007, 2010; M. Salem et al. 2015) are exceptionally well constrained. However, there exists a circular problem: accurate orbital modeling of stellar tracers requires a well-constrained Galactic potential, but the shape of the potential itself is inferred from the motions of these very tracers. Traditionally, the Galactic potential has been inferred by fitting stellar orbits using static, symmetric halo models (J. Binney & S. Tremaine 2008). However, these assumptions are increasingly recognized as oversimplified, leading to biases in mass estimates (F. A. Gómez et al. 2015; D. Erkal et al. 2019; E. Vasiliev et al. 2021; S. Reino et al. 2022) and orbital reconstructions of stellar tracers (e.g., R. D’Souza & E. F. Bell 2022; I. B. Santistevan et al. 2024). An accurate gravitational potential is critical for (1) constraining the MW’s DM distribution—including its total mass, radial profile, and triaxiality—which remains degenerate with halo shape assumptions (A. J. Deason et al. 2019; M. Cautun et al. 2020); (2) reproducing orbits of satellites and stellar streams, whose dynamics encode the nature of DM (e.g., K. V. Johnston et al. 2008; A. P. Cooper et al. 2010); and (3) testing $\Lambda$CDM predictions for halo triaxiality and merger-driven evolution.

High-precision astrometric data from Gaia (Gaia Collaboration et al. 2018) have revealed that the MW’s DM halo is in disequilibrium, induced by recent massive mergers such as the LMC, with a mass of ~10%–20% of the MW’s total mass (G. Besla et al. 2007), in the outer halo (N. Garavito-Camargo et al. 2019, 2021; M. S. Petersen & J. Peñarrubia 2020, 2021; C. Conroy et al. 2021; M. Cavieres et al. 2025) and the Sagittarius dwarf galaxy in the Galactic disk (K. V. Johnston et al. 1995; C. F. Laporte et al. 2018; J. A. Hunt et al. 2021). Recent work has incorporated time-dependent effects in MW potential models (e.g., D. Erkal et al. 2019; N. Shipp et al. 2021; E. Vasiliev et al. 2021; S. Lilleengen et al. 2023), particularly the influence of the LMC, but often assumes an initially symmetric MW halo. However, to fully understand the global structure of the MW halo and reconstruct orbits (B. Lowing et al. 2011; J. L. Sanders et al. 2020; A. Arora et al. 2024b; J. Nibauer & A. Bonaca 2025), it is essential to capture its intrinsic shape and evolution as well (see J. A. Hunt & E. Vasiliev 2025 for a recent review).

A key challenge in studying halo structure is disentangling the contributions of filamentary accretion and satellite perturbations, both of which distort the halo in complex, highly nonlinear ways. Filamentary accretion imprints large-scale anisotropies on the halo over extended timescales, while interactions with massive satellites induce more rapid, localized distortions that propagate throughout the halo. These competing effects, along with the halo’s own response—which depends on the perturber’s mass, velocity, and orbital phase and the halo’s internal structure—create feedback that can amplify or suppress

distortions (e.g., F. A. Gómez et al. 2016; C. F. Laporte et al. 2018; N. Garavito-Camargo et al. 2019). This raises a fundamental question for near-field cosmology: once we subtract the well-known LMC signals, can we still recover the halo’s primordial filamentary imprint? Answering it demands a framework capable of systematically capturing both transient and long-lived deformations across different spatial scales.

Basis function expansions (BFEs) provide a powerful approach to tackle this problem by decomposing the mass distribution into a sum of basis functions (M. Clutton-Brock 1972), often expressed in terms of spherical harmonics (L. Hernquist & J. P. Ostriker 1992; M. D. Weinberg 1999). This formulation is particularly well suited for nearly spherical systems, where perturbations naturally manifest as angular variations in the potential. In this framework, the dipole term represents the first-order deviation from symmetry, capturing large-scale displacements of the halo’s center of mass or coherent bulk flows induced by external perturbers. The quadrupole term, in contrast, reflects the halo’s intrinsic triaxial shape and/or time-dependent deformations (e.g., mergers, tidal interactions). Even in the absence of external perturbations, a triaxial halo inherently exhibits a nonzero quadrupole moment owing to its anisotropic mass distribution. Higher-order terms further resolve finer structural details and localized perturbations.

BFEs have been successfully applied to model evolving halo potentials with and without LMC-mass satellites (e.g., B. Lowing et al. 2011; E. Vasiliev 2013, 2023; J. L. Sanders et al. 2020; N. Garavito-Camargo et al. 2021; A. Arora et al. 2022; M. S. Petersen et al. 2022a, 2022b), allowing systematic tracking of time-dependent deformations across different spatial scales. Perturbations from massive satellites can induce long-lived dipole distortions in otherwise spherical halos (M. D. Weinberg 2023). Indeed, there is growing evidence that the MW halo exhibits a strong dipolar response to the infalling LMC, as seen in constrained MW-LMC simulations (N. Garavito-Camargo et al. 2021; E. Vasiliev et al. 2021; S. Lilleengen et al. 2023), cosmological zoom-in simulations of MW-LMC analogs (e.g., A. Arora et al. 2024a; N. Garavito-Camargo et al. 2024), and even observational data (C. Conroy et al. 2021; M. S. Petersen & J. Peñarrubia 2021; A. Byström et al. 2024; V. Chandra et al. 2024; R. Yaaqib et al. 2024). Moreover, the LMC-induced dipole must be accurately modeled and subtracted before one can robustly infer the halo’s shape from stellar kinematics.

In this work, we investigate the gravitational structure of MW-mass halos from the FIRE-2 suite of zoomed cosmological hydrodynamical simulations (P. F. Hopkins et al. 2018; A. Wetzel et al. 2023) using time-evolving BFE models, simultaneously analyzing the temporal-spatial evolution of each harmonic. We focus on the dominant quadrupole mode of the DM halo. The quadrupole captures the overall elongation and asymmetry of the halo, encoding the cumulative impact of both its large-scale structure and dynamical response to satellites. Unlike higher-order terms, it is robust against small-scale noise while still encoding global evolution. In a companion paper (E. Darragh-Ford et al. 2025, in preparation), we systematically compare MW-LMC analogs from the MW-est suite (D. Buch et al. 2024) to quantify the dipolar (bulk) response of MW-mass halos to massive satellites. Additionally, we apply a novel data-mining approach introduced in M. D. Weinberg & M. S. Petersen (2021) and M. S. Petersen et al. (2022b) to disentangle the effects of filamentary
accretion from those induced by satellite perturbations on the quadrupole moment. This distinction is crucial for determining whether halo distortions arise primarily from smooth cosmological accretion or discrete merger-driven perturbations, thereby providing a more complete picture of halo evolution.

The paper is organized as follows: In Section 2, we provide an overview of the simulations (Section 2.1) and basis function expansion techniques (Sections 2.2 and 2.3), and we introduce the spectrum analysis technique (Section 2.4). In Section 3, we describe the gravitational power spectrum of these halos, identifying the dominant dipole and quadrupole distortions and how they evolve in time (Sections 3.2 and 3.3). In Section 4, we quantify the connection between quadrupolar distortions and filamentary accretion; in Section 5, we separate these large-scale effects on the quadrupole gravitational structure from perturbations caused by satellite mergers. Finally, we discuss our findings and conclusions in Section 6.

## 2. Simulations and Methods

In this section, we outline the key methodologies used in this work. We begin by describing the galaxy simulations utilized in the analysis (Section 2.1). Next, we provide an overview of the BFE framework (Section 2.2), followed by a detailed discussion of how the expansions are fitted to the simulations and an assessment of the quality of the reconstructed density profiles (Section 2.3). Finally, we introduce the spectral analysis technique, used to analyze the BFE coefficient time series (Section 2.4).

### 2.1. Milky Way Analogs

We use the Latte suite of zoomed-in cosmological-baryonic simulations of isolated MW-mass galaxies, part of the Feedback In Realistic Environments (FIRE-2) project (A. R. Wetzel et al. 2016; A. Wetzel et al. 2023).<fnref n="16" /> These simulations employ the FIRE-2 physics model (P. F. Hopkins et al. 2018) implemented within the GIZMO code (P. F. Hopkins 2015), which combines a TREE+PM solver for gravity with a Lagrangian meshless-finite-mass hydrodynamics solver that adapts spatial resolution to local density. The FIRE-2 model incorporates detailed star formation and stellar feedback physics based on Starburst99 stellar evolution models (C. Leitherer et al. 1999) and assumes a ΛCDM cosmology consistent with Planck results (Planck Collaboration et al. 2016). For further details, see P. F. Hopkins et al. (2018) for the FIRE-2 model and A. Wetzel et al. (2023) for the specific implementation.

The Latte suite comprises seven simulated halos—m12b, m12c, m12f, m12i, m12m, m12r, and m12w—that resemble the MW in stellar mass, gas content, DM mass, and density profiles (S. Garrison-Kimmel et al. 2018; P. F. Hopkins et al. 2018; R. E. Sanderson et al. 2020). These halos have total masses of approximately (1–2) × 10¹² M☉ and are resolved with initial particle masses of $m_b = 7100 M☉$ for stars and gas and $m_{DM} = 35,000 M☉$ for DM in the zoomed-in region. Snapshots of the simulations are saved at intervals of approximately 25 Myr for the past 7.5 Gyr of evolution. While the Latte suite focuses on isolated halos, there exists the Elvis suite<fnref n="17" /> that

provides paired MW–Andromeda mass simulations, which better represent the real MW environment. However, analyzing the effects of massive neighbors on the host halo’s shape is beyond the scope of this work. Future studies should incorporate paired simulations to disentangle these influences.

The cosmological environments in the Latte suite ensure realistic accretion histories, where mergers and interactions with the surrounding large-scale structure play a critical role in establishing the density distributions of the halos, assembly histories, and angular momentum profiles (e.g., V. Springel et al. 2005; B. Allgood et al. 2006; G. Despali et al. 2014).

These zoomed-in cosmological simulations are run in a nonperiodic box frame on an expanding background. Since the zoom-in regions do not have zero total momentum, we recenter the simulation onto the host galaxy frame using the shrinking spheres method (C. Power et al. 2003) to find the galactic center at each time step applied to star particles. Subsequently, we rotate all snapshots of the simulation to align the galactic disk with the XY plane at the present day—the principal axes. We define the host-centered rotation in the principal axes as the galactocentric frame.

We identify DM subhalos in these simulations using the ROCKSTAR halo finder (P. S. Behroozi et al. 2012a). The merger trees linking these subhalos across the simulation snapshots are constructed using the consistent-trees algorithm (P. S. Behroozi et al. 2012b), as detailed in J. Samuel et al. (2021). The consistent-trees code reliably tracks subhalos with masses $\ge 10^6 M_{\odot}$ at the fiducial resolution of the FIRE-2 simulations (J. Samuel et al. 2021). N. Garavito-Camargo et al. (2024) present the key properties of the most massive mergers for these halos within the past 7.5 Gyr, including satellite and host masses at infall and pericenter, as well as pericentric distances. A. Arora et al. (2024a) and N. Garavito-Camargo et al. (2024) identify m12b as the best analog for an MW–LMC system based on its satellite-to-host mass ratio (1:8) and the orbital parameters—such as pericentric distance (~40 kpc) and orbital eccentricity (0.83)—of its most massive merging satellite.

### 2.2. Characterizing Halos with Basis Function Expansions

To analyze how the DM halo’s shapes and density evolve over time in our simulations, we expand each N-body snapshot onto a spherical harmonic basis (see below). BFEs provide a flexible and accurate means of capturing halo deformations (e.g., N. Garavito-Camargo et al. 2021; E. Vasiliev et al. 2021), including those induced by galactic evolution and satellite mergers. They have been shown to effectively reproduce orbits, even in the presence of massive satellites, in both constrained (N. Garavito-Camargo et al. 2021; S. Lilleengen et al. 2023; E. Vasiliev 2024) and cosmological simulations (B. Lowing et al. 2011; J. L. Sanders et al. 2020; A. Arora et al. 2022, 2024b). To model the roughly spherical DM halo shape, these expansions are constructed as bi-orthogonal functions of 3D spherical coordinates, where the angular dependence $(\theta, \phi)$ is modeled using spherical harmonics and the radial dependence $(r)$ is captured by orthogonal radial functions.

The radial functions are derived by numerically solving the Poisson equation, which is a specific form of the Sturm-Liouville equation (M. D. Weinberg 1999). This approach provides flexibility in capturing arbitrary spherical profiles without the need for a parameterized zero-order model. A key advantage of this empirical basis is its ability to accurately

<footnote n="16">¹⁶ These simulations are publicly available (A. Wetzel et al. 2023) at http://flathub flatironinstitute.org/fire.</footnote>

<footnote n="17">¹⁷ The Elvis suite has a higher resolution, with $m_b = 3500–4000 M_{\odot}$ for stars and gas and $m_{DM} = 19,000 M_{\odot}$ for DM, compared to the Latte suite (A. Wetzel et al. 2023).</footnote>
represent both the inner and outer mass profiles of halos. Predefined bases, such as the Hernquist basis (L. Hernquist & J. P. Ostriker 1992), suffer from force biases due to mismatched potential gradients (e.g., arising from the Hernquist basis's $r^{-4}$ density falloff vs. the $r^{-3}$ decline of NFW-like halos). Such biases are particularly critical for studies extending to the halo outskirts ~$R_{vir}$ (M. S. Petersen et al. 2022b). By employing an empirical radial basis—free from assumptions about the density law—our method avoids these limitations, making it particularly effective in capturing the nonlinear evolution of the density field, including perturbations induced by satellite mergers and filamentary accretion.

The BFE represents the density field as

$$\rho_{\text{halo}}(\mathbf{x}, t) = \sum_{\ell m n} C_{\ell m n}(t) \varrho_{\ell n}(r) Y_{\ell}^{m}(\theta, \phi), \quad (1)$$

where $\varrho_{\ell n}(r)$ are the time-independent bi-orthogonal radial density functions, $Y_{\ell}^{m}(\theta, \phi)$ are the spherical harmonics, and $C_{\ell m n}(t)$ are the time-dependent weighting amplitude coefficients. Here $|m| \le \ell$, with $\ell$ denoting the angular degree and $n$ the radial index of the polynomials. A paired potential can be derived from the density field by solving the Poisson equation, $^{18}\Delta\Phi = 4\pi G\rho$, where the potential $\Phi(\mathbf{x}, t)$ is related to the density $\rho(\mathbf{x}, t)$ via the Green's function. The acceleration field can then be obtained by taking the gradient of the potential, $\mathbf{a} = -\nabla\Phi$.

We make use of the publicly available code EXP and the associated Python library pyEXP (M. S. Petersen & M. D. Weinberg 2025) to compute and analyze the BFE of the FIRE-2 Latte suite of simulations. A comprehensive review of the BFE mathematical background was presented in M. S. Petersen et al. (2022b).

## 2.3. Fitting BFE to Simulations

Here we provide details on how the BFE is fit to the zoomed cosmological halos in our simulations. To analyze the shape and evolution of the DM halo and its response to infalling satellites, we focus exclusively on the smooth DM distribution of the primary halo—excluding bound substructures such as subhalos and massive satellite debris. For each snapshot, we filter out all DM particles associated with subhalos with $M_{vir} \ge 10^6 M_{\odot}$ from the N-body particle data. Importantly, this filtering is applied at each snapshot, meaning that as subhalos tidally strip and their contributions become part of the main halo, they are no longer excluded once the subhalo catalog loses track of them. Additionally, we remove the debris from the most massive satellite in each simulation across all snapshots. Debris removal is accomplished by identifying the peak mass snapshot of each merging satellite using merger trees (see Section 2.1) and excluding all associated DM particles. As a result, our halos do not account for the growth driven by these most massive infalling satellites. For m12f, which has two massive satellites merging in succession, we remove the debris from both across all snapshots.

This procedure is not without limitations. First, subhalo tracking in ROCKSTAR can be unreliable, especially for

heavily stripped or disrupted subhalos, which may result in incomplete removal of substructure. This is particularly relevant for our most massive mergers, where tidal debris may be incorrectly classified as part of the main halo rather than the infalling satellite (B. Diemer et al. 2024; P. Mansfield et al. 2024). Additionally, by removing the anisotropically accreted debris from these satellites, we are reducing the imprint of their accretion history on the halo's global shape. However, this choice is justified, as our goal is to isolate the intrinsic response of the MW-mass halo to the infalling satellites, rather than the persistent anisotropic distortions left by their debris.

This approach minimizes the impact of substructure, ensuring that the BFE signal is not dominated by the hard-to-model density field of small subhalos, which could introduce noise in high-order radial polynomial coefficients. It ensures that the BFE captures the halo's intrinsic response to its environment and perturbations, while still accounting for halo growth—except for that driven by the most massive infalling satellites. These filtered halos are directly used for the BFE fits and to compute the radial basis functions.

We construct the radial basis functions using a spherically averaged density profile extending to 600 kpc (~$2R_{vir}$) averaged in time over the past 7.5 Gyr. This ensures that the radial functions remain fixed throughout the halo's evolution. A fixed basis is critical for studying the temporal evolution of spatial power: if the basis varied between snapshots, comparisons of structure across spatial scales would become ambiguous. In Appendix A, we demonstrate the stability of the radial basis over time, showing that it remains consistent in time except during massive mergers ($M_{sat}/M_{host} \ge 0.1$), further justifying this approach.

We truncate the expansion to $\ell_{max} \le 4$ for the angular components and $n_{max} \le 10$ for the radial components.<fnref n="19" /> Figure 1 illustrates the orthogonal radial basis constructed for m12i, with polynomials (color-coded) shown as basis potential (normalized with $GM_{tot}/r_{scale}$) as a function of radial distance from the galactic center. Each panel corresponds to a different value of $\ell$ shown up to the octupole ($\ell = 3$), while increasing $n$ is indicated by redder shades within each panel. The basis functions satisfy bi-orthogonality conditions given by

$$\int_{0}^{\infty} \varrho_{\ell n}(r) \psi_{\ell n'}(r) w(r) dr = \delta_{nn'}, \quad (2)$$

where $w(r) \equiv -4\pi Gr^2$ is the weighting function, $\psi_{\ell n}(r)$ are the time-independent radial potential functions, and $\delta_{nn'}$ is the Kronecker delta, ensuring bi-orthogonality for fixed $\ell$. The radial basis functions span scales from a few kiloparsecs (high-order polynomials) to hundreds of kiloparsecs (low-order terms), with node spacings tuned to resolve gravitational potential variations across the full spectrum of dynamically relevant scales—from small, localized perturbations to large-scale halo asymmetries. Finally, we fit these expansions to the recentered filtered N-body particle data for the DM halo at each snapshot for all the simulations, using particles within 600 kpc of the galactic center over the past 7 Gyr of evolution (roughly 350 snapshots).

<footnote n="18">¹⁸ The alternative notation $\nabla^2\Phi = 4\pi G\rho$ is occasionally used, though not strictly accurate! Mathematically, the Laplacian operator is written as $\Delta$, which by definition is the divergence of the gradient; $\nabla^2$ can be misinterpreted as taking the gradient twice, which yields a rank-2 tensor.</footnote>

<footnote n="19">¹⁹ While the expansion fits are computed up to $\ell_{max} \le 20$, we later show that most of the dynamical signatures and gravitational power are contained within 4.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Orthogonal radial basis functions derived by numerically solving the Poisson equation using the time-averaged spherically averaged density profile for m12i. Each panel corresponds to a different angular degree ($\ell = 0, 1, 2, 3$), and the radial basis functions are shown as normalized potentials as a function of radial distance. Polynomials are color-coded, with redder shades indicating higher radial orders ($n$). The basis function nodes span scales of a few kiloparsecs to hundreds of kiloparsecs, which gives resolving power for structure at all scales of interest in the problem.</figcaption></figure>

Figure 2 shows the spherically averaged density profile of m12w (solid black lines) and its reconstruction using the monopole ($\ell = 0$) term from the BFE (dashed red lines) over the past 7 Gyr of evolution. We highlight m12w as a critical test case owing to its extreme merger history: a 1:8 mass ratio satellite with eccentricity $e = 0.96$ undergoes a pericentric passage within 10 kpc at 8.7 Gyr. Such mergers induce rapid, nonadiabatic changes to the halo's potential, challenging the fidelity of basis function reconstructions.

The profiles are shown at 1 Gyr intervals, with the corresponding times labeled in the upper right corner of each panel. Subpanels below each main panel show the percentage difference between the reconstructed and true density profiles. Despite the merger's disruptive effects, reconstructions are accurate to within 5% out to the virial radius of the halo at all times. Larger discrepancies (8%–10%) are observed near the galactic center 2–3 kpc and beyond the virial radius at earlier times, reflecting the limitations of the monopole term in capturing steep central gradients and rapidly accreting outskirts. Notably, uncertainties remain stable even near the merger pericenter 8.7 Gyr, demonstrating the method's resilience to transient perturbations.

For comparison, simulations without significant mergers (e.g., m12i) maintain uncertainties within 5% across all radii and times (not shown). These results align with A. Arora et al. (2022, 2024b), who reported similar reconstruction errors in FIRE simulations using BFE with a grid-based radial basis, even during mergers. The consistency across halo histories underscores that a fixed radial basis robustly reconstructs time-evolving density profiles, enabling reliable temporal comparisons of expansions despite nonadiabatic events.

The time series of BFE coefficients encode rich dynamical information, including spatial and temporal correlations. Different dynamical events, such as the accretion of massive satellites or filamentary material, occur on distinct timescales

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Spherically averaged density profiles for m12w (solid black lines) and reconstructions using the monopole ($\ell = 0$) term from the BFE (dashed red lines) at 1 Gyr intervals over the past 7 Gyr. Times are labeled in the upper right corner of each panel. Subpanels show the percentage difference between the reconstructed and true profiles. Reconstructed profiles are accurate to within 5% errors out to the virial radius at all times.</figcaption></figure>

and can leave imprints at specific frequencies in the BFE coefficients. Spectral analysis techniques are well suited to disentangle these contributions and identify the dominant modes of evolution. Singular spectrum analysis (SSA) is a proven method for extracting correlated signals from time series by examining their covariance structure at successive time lags (N. Golyandina et al. 2001). By extending SSA to a multichannel framework, we can analyze the BFE coefficients —treated as multichannel time series—to uncover key modes of correlated evolution in the halo’s dynamics. This approach allows us to isolate and interpret the signatures of different physical processes shaping the halo over time.

## 2.4. Decomposing Halo Evolution and Structure with mSSA

Multichannel SSA (mSSA) is a technique for identifying correlated spatial–temporal patterns across multiple time series (N. Golyandina et al. 2001). Unlike traditional SSA, which focuses on a single time series, mSSA considers multiple channels simultaneously—in our case each coefficient series $C_{\ell nm}(t)$, enabling the discovery of patterns shared across different coefficients. This makes it particularly well suited for studying $C_{\ell nm}(t)$, which capture the dynamical evolution of a halo through their changing amplitudes.

The method involves constructing a covariance trajectory matrix that incorporates time-lagged copies of each coefficient series over a specified window length $L$, capturing both spatial and temporal correlations. For our analysis, we adopt half the length of the time series (maximum possible length), roughly 3.5 Gyr. The advantage of the longest window is that it is sensitive to slow trends (e.g., filamentary accretion and satellite orbits) while suppressing high-frequency noise. This effectively acts as a low-pass filter, prioritizing signals on cosmological timescales over transient perturbations. Singular value decomposition (SVD) is then applied to this matrix, decomposing the data into orthogonal principal components (PCs) ranked by their singular values. These modes represent the dominant correlated patterns, providing insight into the system’s evolving structure. For a detailed description of the mathematical implementation and applications of mSSA
+BFE, we refer the reader to M. D. Weinberg & M. S. Petersen (2021) and Appendix C of A. C. Johnson et al. (2023).

mSSA shares similarities with principal component analysis (PCA) but is distinguished by its focus on temporal patterns. While PCA identifies dominant spatial correlations from a covariance matrix built from a single field variable sampled across multiple spatial locations, SSA constructs its covariance matrix using a single time series sampled over temporal windows. The resulting temporal PCs capture the primary modes of variation across time, making SSA ideal for studying time-varying systems (M. D. Weinberg 2023). Notably, mSSA generalizes PCA by incorporating temporal lags: when the window length $L=1$, mSSA reduces to PCA. For BFE coefficient series—which encode spatial structure through spherical harmonics—mSSA can be interpreted as a hybrid spatiotemporal decomposition, resolving correlated patterns across both space (via harmonic modes) and time (via lagged covariance).

mSSA has proven highly effective in recent studies. M. D. Weinberg & M. S. Petersen (2021) demonstrated its utility in analyzing galactic bar formation, uncovering mode coupling and pattern-speed decay, as well as a slow retrograde oscillation of the disk with a gigayear-scale period in idealized MW simulations. Similarly, M. D. Weinberg (2023) used mSSA to identify dipole instabilities in NFW-like DM halos, while A. C. Johnson et al. (2023) revealed two distinct disk-halo dipole modes in an MW-like galaxy: one that grows and another that decays over time, driven primarily by the halo.

In this work, we apply mSSA to the BFE coefficients to reconstruct dynamically meaningful gravitational fields and isolate specific contributions to the halo's morphology. In Section 5, we show how mSSA disentangles the quadrupole response induced by the LMC analog from intrinsic halo triaxiality and filamentary influences in m12b.

### 3. Gravitational Power in the Harmonics

The gravitational dynamics of DM halos can be studied by decomposing the density and potential fields into harmonics using BFE. This approach provides the ability to calculate gravitational power in each harmonic ($\ell$) independently, enabling a detailed exploration of the distinct dynamical features, symmetries, and phenomena captured by individual modes. Each $\ell$ harmonic encodes distinct structural and dynamical signatures of the halo. For instance, the dipole term reflects asymmetries induced by merging satellites and the motion of the center of mass (M. D. Weinberg 1998), and the quadrupole mainly characterizes the halo's global shape, including triaxiality and deviations from spherical symmetry (C. A. Vera-Ciro et al. 2011; M. D. Schneider et al. 2012). It can also result from morphological distortions associated with the resonant response of the halo to a passing satellite (M. D. Weinberg 1998; J.-H. Choi et al. 2009; E. D'Onghia et al. 2010). Higher-order harmonics, though generally less significant in amplitude, reveal much more fine-grained features such as subhalo-induced perturbations or tidal streams (B. Lowing et al. 2011).

To quantify the gravitational contributions from these different harmonic orders, we compute the gravitational power for each mode. This metric captures the spatial distribution of the gravitational field across the harmonic spectrum and provides a tool to study its evolution over time. We define the gravitational power $P = -2*W$, where W is the energy in the

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Power in each harmonic averaged over the past 7.5 Gyr of the halo evolution relative to the monopole ($\ell=0$) for all simulations. The largest amplitudes (darker shades, and less negative values) of power are noted in the dipole ($\ell=1$) and quadrupole ($\ell=2$) modes, with the quadrupole power being systematically higher in all the simulations except m12w, which has a highly eccentric radial merger (with eccentricity $e=0.96$) with an LMC-mass satellite.</figcaption></figure>

gravitational field. For a DM halo with density, $\rho(x)$, and gravitational potential, $\Phi(x)$, $P$ is given by

$$P = - \int \rho(x) \Phi(x) dx. \quad (3)$$

Here the gravitational constant is absorbed into the weighting factor of the BFE (see Equation (2)). In the BFE framework, the power in a specific harmonic, $P_\ell$, is calculated as

$$P_{\ell} = \sum_{n,m} C_{\ell,n,m}^{2}, \quad (4)$$

where $C_{n,m}$ are the coefficients associated with the $\ell$th harmonic for all radial polynomial indices $n$ and azimuthal indices $m$.

Each harmonic has a distinct physical interpretation. The monopole term represents the spherically symmetric component of the density distribution, capturing the overall mass distribution and radial profile of the halo. The dipole term reflects asymmetries in the density field, often associated with the displacement of the halo center due to external perturbations or interactions with nearby structures. The quadrupole term encodes the halo's shape anisotropy and is sensitive to tidal forces, filamentary accretion, and interactions with the large-scale cosmic web. For a more detailed discussion of these modes and their implications, we refer the reader to our companion paper (E. Darragh-Ford et al. 2025, in preparation). We note that, as we are primarily interested in the response of the main halo, we have removed debris from the most massive satellite from the particle data (see Section 2.3). This removal systematically lowers the fractional gravitational power and asymmetry, ensuring that our analysis more realistically captures the shape and dynamics of the colder component of the main halo, leaving behind material that is strongly unmixed.

#### 3.1. Halo Response Captured in Dipole and Quadrupole

Figure 3 shows the time-averaged gravitational power in each harmonic mode ($\ell_{max} \le 4$) relative to the monopole term, across all simulated halos (with diverse assembly histories) in the Latte suite. The dipole ($\ell=1$) and quadrupole ($\ell=2$) consistently dominate the harmonic spectrum capturing the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Top row: fractional gravitational power as a function of time in the dipole (dotted) and quadrupole (solid) for m12b (blue, simulation with an LMC analog) and m12i (indigo). The quadrupole power exhibits long-term evolution, with a gradual decline over time, while the dipole power shows short-term variability on merger timescales, peaking after pericentric passages. Bottom row: total distance of the LMC-analog satellite from the halo center in m12b as a function of time. Dashed vertical lines mark the times of the first pericentric passages.</figcaption></figure>

most relevant dynamical features of the halo, with their amplitudes almost a magnitude higher than those of the higher-order modes ($\ell \ge 3$). These $\ell \ge 3$ poles remain relatively insignificant compared to the first two harmonics. The dipole naturally arises from the motion of the halo's center of mass, which can be perturbed by mergers, or external tidal forces. In contrast, the quadrupole reflects the halo's overall shape, including triaxiality and the dynamical effects of mergers and interactions that modify the halo's prolateness or oblateness.

The $\ell = 2$ systematically shows higher time-averaged power in all simulations except m12w, where the $\ell = 1$ and $\ell = 2$ exhibit comparable power amplitudes. The elevated $\ell = 1$ in m12r and m12w is a result of an LMC-mass satellite. In m12w, this satellite reaches a pericentric distance of 10 kpc, inducing a strong asymmetric response in the halo and driving a significant increase in the $\ell = 1$ power.

Overall, the $\ell = 1$ and $\ell = 2$ modes capture the most gravitationally relevant contributions to the halo's large-scale structure and dynamics (T. Quinn & J. Binney 1992; D. J. Eisenstein & A. Loeb 1995; D. Aubert et al. 2004; C. A. Vera-Ciro et al. 2011). These two modes are the primary drivers of the halo's response to mergers, asymmetries, and shape deviations, offering a comprehensive view of the halo's evolution. Given their dominance in the gravitational power spectrum, we focus our analysis on these two modes.

## 3.2. Transient Dipole and Long-term Quadrupole

Figure 4 illustrates the temporal evolution of gravitational power in the dipole (dotted) and quadrupole (solid) modes for two simulations: m12b (blue, with an LMC analog) and m12i (indigo). The top row shows the fractional power for each mode, while the bottom row tracks the total distance of the LMC-analog satellite from the center of the halo in m12b as a function of time. The dashed vertical line marks the time of first pericentric passage.

The dipole power exhibits transient variability on short timescales, reflecting merger-induced asymmetries. In m12b, it peaks immediately after the pericentric passage of the LMC

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Top row: fractional gravitational power in the quadrupole mode as a function of time for all simulations (color-coded) over the full 7 Gyr period (left panel) and zoomed into a ±1 Gyr window around the first pericentric passage of massive satellites (right panel). Dashed vertical lines mark the times of pericentric passages. Bottom row: distance of the most massive satellite from the halo center in each simulation (color-coded) as a function of time (left panel) and shifted to align the first pericentric passage time at 0 Gyr (right panel). The quadrupole power generally decreases as the halo evolves, but it shows transient increases near pericentric passages due to tidal wakes and direct contributions from the satellites.</figcaption></figure>

analog, as the halo responds to the strong perturbation. By contrast, m12i shows relatively lower dipole power after an initial peak from earlier mergers, consistent with the absence of ongoing strong interactions.

The quadrupole power, in contrast, evolves more gradually. It starts at a higher magnitude and steadily decreases over time as the halo dynamically settles. In m12b, the quadrupole also exhibits a secondary increase just before the LMC analog's pericentric passage, driven by the induced DM wake and torque from the satellite. This behavior highlights the quadrupole's sensitivity to long-term shape and structural changes. We note similar trends in $\ell = 1$ and $\ell = 2$ power in all simulations.

## 3.3. Declining Quadrupole across Halo Evolution

Figure 5 shows the temporal evolution of $\ell = 2$ power across all simulations in the top row, with the corresponding orbital histories of massive satellites shown in the bottom row. The left panels provide an overview of the full 7 Gyr period, while the right panels zoom in on the ±1 Gyr window around the first pericentric passage. The dashed lines indicate the times of pericentric passages for the satellites.

The quadrupole power here is primarily determined by the global shape of the halo and its triaxiality, gradually decreasing over time as the halo relaxes dynamically. This long-term decline reflects the halo's transition from an actively evolving asymmetric state to a stabler symmetric configuration (see Section 4 and Appendix B). However, short-term fluctuations in the quadrupole power are observed near the pericentric passages of satellites. These peaks are attributed to three mechanisms: (1) The dynamical friction wake induced by the satellite as it interacts with the host DM halo. (2) The gravitational influence of the satellite, which torques the host and induces a transient quadrupole moment. (3) The satellite's passage excites resonant oscillations in the halo's natural modes, leading to quadrupole distortions that decay over shorter timescales. While the torque mechanism aligns with the halo's principal axes, the resonant response arises from the halo's natural modes being perturbed by the satellite's gravitational field. Unlike the transient torque-driven signal,
resonant oscillations can contribute to the halo’s long-term shape evolution. This distinction is critical for disentangling short-term perturbations from sustained structural changes.

Our BFE models reconstruct the halo’s density field within a spherical boundary of 600 kpc (Section 2.3), with gravitational power computed using the full set of radial polynomials. To test sensitivity to radial truncation, we repeat the analysis for boundaries of 350 and 125 kpc. While the absolute quadrupole amplitude decreases slightly at smaller scales (e.g., 15% at 125 kpc), the long-term decline in quadrupole power and transient merger-driven enhancements persist across all radial limits. This consistency confirms that the observed quadrupole evolution reflects intrinsic dynamical processes, not artifacts of the chosen boundary.

To summarize, the dynamical state and shape of DM halos can be effectively studied by projecting the phase-space particle data into spherical harmonic expansions and computing their contributions to the total gravitational power. On average, the dipole and quadrupole modes dominate the harmonic spectrum. The dipole captures transient features on short timescales, while the quadrupole represents the large-scale shape and dynamical evolution of the halo. As we demonstrate in subsequent sections, the quadrupole’s orientation and amplitude are closely tied to the halo’s interaction with filamentary structures in the cosmic environment. Although both modes can exhibit peaks due to mergers, the quadrupole power consistently is higher and decreases over time, reflecting its long-term dynamical stability.

In the next sections, we will quantify how filamentary structures and satellites establish both the magnitude and orientation of the quadrupole. We will also investigate the reasons for its declining amplitude over time, providing insights into the processes that govern halo evolution and its interaction with the surrounding cosmic web.

## 4. Filaments’ Imprint on the Quadrupole

The shapes of DM halos are intrinsically set by their interactions with the surrounding large-scale structure, such as filaments and sheets of the cosmic web—a long-studied phenomenon in cosmology dating back to seminal work on structure formation (F. Hoyle 1953; P. J. Peebles 1969; A. Doroshkevich 1970), angular momentum acquisition (S. D. White 1984), and anisotropic accretion (D. J. Eisenstein & A. Loeb 1995). Modern simulations have reaffirmed and quantified this connection, demonstrating how filaments govern angular momentum and mass accretion (D. Aubert et al. 2004; C. A. Vera-Ciro et al. 2011). This alignment reflects the interplay between a halo’s mass assembly history and external tidal forces, with the large-scale environment leaving a lasting imprint on the halo’s morphology and quadrupole moments.

In this section, we explore the origin and evolution of the dominant quadrupole mode in these simulations. We investigate how the alignment and magnitude of the quadrupole moment correlate with the surrounding large-scale structure, specifically focusing on the filaments feeding into the halo, the weakening of these structures over time, and the evolution of the alignment.

### 4.1. Filaments Set the Quadrupole Amplitude and Orientation

We briefly describe the methodologies used to measure the quadrupole position angle and the orientations of major filaments, which enable the study of their alignment.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Quadrupole density field ($\rho_{\ell=2}$) in the XZ plane, referred to as the “Quadrupole-Environment View,” where the galactic disk is edge-on. Panels depict a time step with low fractional quadrupole power (left; m12b at 13.39 Gyr, $P_2/P_0 = 0.0006$) and high quadrupole power (right; m12r at 5.42 Gyr, $P_2/P_0 = 0.0127$). Red-blue open contours represent $\rho_{\ell=2}$ (positive and negative regions of the quadrupole), computed using the BFE formalism. The background shows the surface density of the environment with the main halo subtracted (diverging color map), highlighting filamentary structures within 1 Mpc. Green arrows indicate the quadrupole orientation, $\hat{n}_{\ell=2}$, while black arrows denote the dominant filament orientation, $\hat{n}_{\text{fila}}$. In the high-power epochs, the filaments are prominent and anisotropic, and the quadrupole is strongly aligned with the dominant filament (right panel), with a 3D dot product of 0.99. In contrast, the low-power epoch corresponds to broader, less defined filamentary structures, the quadrupole power is weaker, and the alignment between $\hat{n}_{\ell=2}$ and $\hat{n}_{\text{fila}}$ is lower (3D dot product = 0.34; left panel).</figcaption></figure>

1. *Quadrupole orientation* ($\hat{n}_{\ell=2}$). The quadrupole position angle is determined using the expansion coefficients from the spherical harmonic decomposition for $\ell=2$ and $m=0, 1, 2$ with the lowest-order radial polynomial ($n=0$). We construct the moment of inertia (I) tensor using ($C_{200}, C_{210}, C_{220}$), I’s eigenvalue decomposition provides the principal axes and the position angle, denoted as $\hat{n}_{\ell=2}$.

2. *Filament orientation* ($\hat{n}_{\text{fila}}$). Filament orientations are computed by identifying the location of the maximum-density/main filament in angular coordinates $(\theta, \phi)$ in the galactocentric frame for particles associated with filamentary structures. These particles are defined as those infalling toward the galactic center (negative radial velocity), outside the virial radius of the main halo, and with particles from the main halo and subhalos removed. The orientation of the main filament is represented by $\hat{n}_{\text{fila}}$.

The position angle of the quadrupole moment ($\hat{n}_{\ell=2}$) and the filament orientation ($\hat{n}_{\text{fila}}$) allow us to quantify their alignment and explore the correlation between the two.

To explore how the large-scale environment shapes the quadrupole mode, we examine the quadrupole density field, $\rho_{\ell=2}$, in the XZ plane (the galactic disk is edge-on) for two distinct epochs: one with low quadrupole power ($P_2/P_0 = 0.0006$) and another with high quadrupole power ($P_2/P_0 = 0.0127$). In Figure 6, the left panel corresponds to a low-power epoch (m12b at 13.39 Gyr), and the right panel corresponds to a high-power epoch (m12r at 5.42 Gyr).

The red-blue open contours in each panel represent the quadrupole density field, $\rho_{\ell=2}$, highlighting positive and negative contributions, computed using the BFE formalism. The background diverging color map<fnref n="20" /> shows the surface density of the environment with the main halo subtracted,

<footnote n="20">²⁰ Rather odd but conscious choice to highlight the structures.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time evolution of the quadrupole-environment view (as described in Figure 6) for the past 7 Gyr of the halo evolution. Each panel depicts a 1 Gyr time step for m12b (marked at the top of each panel). Red-blue open contours represent $\rho_{\ell=2}$ (positive and negative regions of the quadrupole) computed using the BFE formalism. The background surface density highlights filamentary structures within 1 Mpc, with the main halo subtracted. Green arrows indicate $\hat{n}_{\ell=2}$, while black arrows denote $\hat{n}_{\text{fila}}$. The trajectory of the infalling satellite is shown as a dashed black curve, with its current position marked by an open magenta star. The quadrupole initially aligns with the filaments but twists and weakens as the filaments diffuse over time, and a smaller quadrupole is induced during satellite interactions.</figcaption></figure>

revealing filamentary structures within 1 Mpc. The orientation of the quadrupole moment, $\hat{n}_{\ell=2}$, is indicated by the green arrow, while the orientation of the dominant filament, $\hat{n}_{\text{fila}}$, is shown by the black arrow. We introduce this visualization as the “Quadrupole-Environment View,” which will be used in subsequent analyses to investigate the interplay between quadrupole moments and filamentary structures across different simulations and epochs.

In the low-power epoch (left panel), the filamentary structures appear broad and diffuse, as we are in a late-stage evolution of the DM halo, with weaker anisotropies in the surrounding density field. The quadrupole power is multiple orders of magnitude lower, and the quadrupole orientation is weakly aligned with the dominant filament direction,<fnref n="21" /> with a 3D dot product of 0.34. In contrast, during the high-power epoch at higher redshifts (right panel), the filaments are well defined and anisotropic, and the quadrupole orientation is nearly perfectly aligned with $\hat{n}_{\text{fila}}$, with a 3D dot product of 0.99. The difference in quadrupole power between the low- and high-power epochs spans multiple orders of magnitude, underscoring the strong influence of the filamentary structures on both the amplitude and orientation of the quadrupole mode. This pattern is consistently observed across different simulations and times, reinforcing the idea that the galactic environment, particularly the filamentary structures, governs the quadrupole’s behavior.

Next, we investigate how the evolution of the filamentary environment influences and decreases the quadrupole power (Figure 5) and affects the overall shape of the halo.

Figure 7 shows the time evolution of the quadrupole-environment view (similar to Figure 6) for the m12b halo over the past 7 Gyr. Each panel corresponds to a snapshot taken at 1 Gyr intervals, highlighting the evolution of the quadrupole density ($\rho_{\ell=2}$), the surrounding filamentary structures, and the trajectory of an infalling satellite. Red-blue open contours represent $\rho_{\ell=2}$, while the background color map shows the main-halo-subtracted surface density of the environment. The green and black arrows indicate the orientations of the quadrupole ($\hat{n}_{\ell=2}$) and the dominant filament ($\hat{n}_{\text{fila}}$, respectively, and the dashed black curve traces the satellite’s trajectory, with its current position marked by a magenta star.

At early times (e.g., the 6.8 Gyr panel), the quadrupole moment has a higher absolute density value ($\rho_{\ell=2}$) and is closely aligned with the dominant filament, as indicated by the green and black arrows in the upper left corner of the panel. Over time, as the filaments broaden and diffuse, the quadrupole twists and changes direction, while its absolute density weakens significantly. This temporal evolution demonstrates the strong connection of the filamentary structures to the quadrupole’s amplitude and orientation.

Another notable feature is the interaction between the quadrupole and the infalling satellite. Following the satellite’s first pericentric passage at around 8.8 Gyr, a secondary, smaller quadrupole appears in the inner regions, becoming especially evident near the satellite’s second pericenter at 10.8 Gyr. This interaction not only perturbs the alignment of the main quadrupole but also temporarily enhances its amplitude.

These observations underline the importance of disentangling the quadrupole contributions from filamentary structures

<footnote n="21">²¹ Projecting from 3D to 2D affects the apparent alignment.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The absolute value of the dot product between the quadrupole position angle $\hat{n}_{\ell=2}$ and the maximum-density filament orientation $\hat{n}_{\text{fila}}$ for all simulations, shown as color-coded curves (colors are assigned based on the final letter of each halo name for consistency, e.g., m12r is red, m12i is indigo, except for m12w (green) and m12f (yellow)) as a function of time. Vertical dashed lines—two for m12f—also color-coded, mark the times of first pericentric passage of massive satellites. The dot product generally remains close to 1, with variations attributable to mergers or external influences, particularly in halos embedded within massive, broad filaments such as m12m (magenta) and m12c (cyan).</figcaption></figure>

and satellite interactions. A detailed separation of these effects is essential for understanding how the large-scale environment and local dynamics collectively shape the quadrupole moment and, by extension, the halo's triaxiality. In Section 5, we focus on isolating these contributions and assessing their individual impacts.

Finally, we show the time evolution of the quadrupole-environment view for the remaining halos in Appendix B. While our halos exhibit diverse configurations, such as differing numbers of filaments feeding them (e.g., m12i has two filaments, while m12r has three to four), or being embedded within massive, broad filaments (e.g., m12c and m12m), the trends observed here remain consistent. Specifically, the quadrupole weakens over time, diffuses alongside the filamentary structure, and adjusts to small perturbations caused by infalling satellites.

Figure 8 illustrates the temporal evolution of the alignment between the quadrupole orientation ($\hat{n}_{\ell=2}$) and the dominant filament orientation ($\hat{n}_{\text{fila}}$) for all simulations. The absolute value of their dot product is shown as color-coded curves, with each color corresponding to a specific halo (refer to Figure 5 for the legend). Vertical dashed lines, color-matched to the respective halo, indicate the first pericentric passage of massive satellites.

The dot product remains close to 1 across most simulations, indicating a small angle between the filament and quadrupole vectors and demonstrating their strong alignment. Notable deviations occur in halos embedded within broad, massive filaments, such as m12m (magenta) and m12c (cyan). In these cases, the isotropic feeding of material from the surrounding structure weakens the alignment but does not eliminate it. For m12i (indigo), the alignment improves over time as the halo recovers from an earlier perturbation caused by a merger, ultimately achieving excellent alignment at late times.

Simulations experiencing massive mergers show characteristic drops in the dot product after the first pericentric passage,

as seen in m12b (blue) and m12r (red). However, alignment often recovers in the absence of further significant perturbations, as evidenced by m12f (yellow). From Appendix B, it is also noticed that massive satellites preferentially fall along the major filaments feeding the halo (e.g., G. Tormen et al. 1997; N. I. Libeskind et al. 2014; E. Tempel et al. 2015).

To summarize, filaments play a crucial role in setting the quadrupole moment and, by extension, the shape of the halo (e.g., D. Aubert et al. 2004; J. Bailin & M. Steinmetz 2005; B. Allgood et al. 2006; C. A. Vera-Ciro et al. 2011). The strong alignment between the quadrupole and the filaments, even at later times, demonstrates that the global environment has a lasting influence on the halo evolution. While satellites contribute to the quadrupole, particularly during pericentric passages, their influence is subdominant at other times compared to the persistent and dominant role of filaments in determining the overall orientation and strength of the quadrupole moment.

This raises the question of how deep in the halo the influence of filaments extend. Can filamentary structures shape the inner regions of the halo, or are their effects confined to the outskirts? Addressing these questions requires examining the radial extent of filamentary influence, as well as their role in modulating the inner halo shape, which we will do in the next section.

## 4.2. Radial Extent of the Filaments Measured with Quadrupole

To measure the influence of the filaments throughout the halos, we compute the density of the quadrupole terms in the BFE relative to the spherical averaged density at a given radius. The fractional quadrupole density is computed as the maximum density in the quadrupole at a given radius and time, $\rho_{\ell=2}(r)_{\text{max}}$, normalized by the spherically averaged ($\rho_{\text{ell}} = 0(r, t)$) density at the same radius.

To calculate $\rho_{\ell=2}(r)_{\text{max}}$, we use the BFE to compute the quadrupole contribution to the density field on a uniform spherical grid with angles sampled evenly across the sphere by distributing points equidistantly along a spiral using the method described in C. Hüttig & K. Stemmer (2008). This approach avoids the biased sampling inherent in traditional spherical grids (e.g., uniform $\phi$ and $\cos\theta$ grids), which disproportionately concentrate points near the poles,<fnref n="22" /> leading to uneven angular resolution. The maximum quadrupole density at each radius represents the strongest aspherical contribution to the total density. The total density includes all multipoles so values of this ratio between the multiple and monopole contributions can exceed 1 for densities with strong $\ell > 0$ multipoles. This fractional quadrupole density serves as a proxy for the strength of the aspherical component set by the environment of the halo at different radii and epochs.

Figure 9 shows a heat map of the fractional quadrupole density as a function of spherical radius $r$ and time for all simulated halos. Values greater than 1 indicate regions where the quadrupole mode dominates over the spherically symmetric component, signaling highly triaxial deviations from spherical symmetry, while values near 0 correspond to spherical configurations. The top row shows simulations where the majority of the evolution occurs without massive mergers, while the bottom row depicts halos with significant satellite mergers. In each panel, satellite trajectories are

<footnote n="22">²² Traditional grids are area uniform but angle biased.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The maximum value of the absolute density in the quadrupole component $\rho_{\ell=2}$ at each spherical radius $r$, normalized by the spherically averaged density $\rho_{\ell=0}(r)$, shown as a function of radius and time across all simulations. The top row features simulations without massive mergers for most of their evolution, while the bottom row includes simulations with significant satellite mergers. Satellite trajectories are plotted as solid black lines, and the growing virial radius of the main halo is shown as a dashed black line. Contour levels are marked at 0.5 (yellow) and 1 (red), indicating regions of enhanced quadrupole density relative to the monopole. Fractional enhancements exceeding 1 are observed near the virial radius in most simulations, reflecting the strong influence of the large-scale filamentary environment. In all simulations, a persistent enhancement of at least 0.5 is observed as low as 30 kpc, highlighting the deep influence of the filaments on the inner halo. Localized increases in quadrupole density are also evident along satellite trajectories, corresponding to perturbations caused by mergers.</figcaption></figure>

marked by solid black curves, and the growing virial radius of the halo is shown as a dashed black line. Contour levels are drawn at fractional quadrupole densities of 0.5 (yellow) and 1 (red).

The $\rho_{\ell=2}(r)_{\text{max}}$ approaches values near 1 close to the virial radius (red contours) in most halos, signaling strong anisotropy dominated by filamentary accretion. Exceptions include m12m—embedded in a single massive, broad filament that accretes matter somewhat isotropically at the virial radius—and m12w, where the primary symmetry is not captured by a quadrupole (see Appendix B). In the other halos, early epochs (~6 Gyr, virial radius is ~150 kpc) exhibit high quadrupole amplitudes (~1) at the virial radius. These amplitudes persist as the halo grows, demonstrating the retention of a “filamentary memory”: structural imprints of anisotropic accretion that endure despite cosmic evolution. Simulations without major mergers (top row) show particularly stable alignment between the quadrupole and the virial radius, highlighting the virial radius as a dynamical transition boundary where isotropic accretion gives way to anisotropic filamentary influence.

The inner regions (~30 kpc) in all halos retain a persistent quadrupole (≥0.5; yellow contours), demonstrating the radial reach of filamentary influence into the inner halo morphology. However, in halos experiencing massive satellite mergers (bottom row), this filamentary memory is reconfigured by localized perturbations near the satellite trajectories. For example, in m12w, regions with higher fractional amplitudes correlate closely with the distances and times of satellite passages, particularly near pericenter, where the quadrupole enhancements are most pronounced. These perturbations form a secondary quadrupole signature superimposed on the dominant filamentary-driven mode. Importantly, the two contributions often interfere constructively, as massive satellites are typically accreted along the most massive filaments, reinforcing the overall quadrupole signal. Despite this, the satellite-induced contributions remain subdominant to the global filamentary mode. A similar increase in the amplitudes

of BFE coefficients during satellite mergers was noted by A. Arora et al. (2022).

The global influence of filaments evolves over cosmological timescales, while the localized effects of satellites are tied to their orbital dynamics. For example, Figure 5 shows quadrupole power peaking near satellite pericenters, while Figures 7 and 9 highlight sustained filamentary-driven quadrupole strength alongside localized increases near satellite trajectories. These distinctions in timescales allow us to disentangle the effects of mergers and filaments. In the next section, we use m12b as a case study to demonstrate how these signatures can be analyzed separately.

## 5. Decoupling the Effect of Filaments and Satellites on the Halo Shape

In this section, we apply multichannel SSA (mSSA) to the quadrupole coefficient time series in m12b—a halo that exhibits one of the strongest transient and merger-driven quadrupole signals in our sample—to identify temporally correlated gravitational signatures. We consider all the coefficients associated with the quadrupole term, including the full set of coefficients ($C_{\ell=2,m,n}$), where $m = 0, 1, 2$ and $n$ spans all the radial polynomials (Section 2.2). The inclusion of a broad range of polynomials enables the capture of distortions across multiple radial scales, from several hundred to a few kiloparsecs (see Figure 1). In Section 5.1 we discuss how to decompose the competing effects from the filaments and the satellites. In Section 5.2 we showed how once the PCs are identified one can reconstruct the density fields associated with the response of the halo to the filaments and to the satellite.

### 5.1. Spectral Analysis to Decouple the Quadrupole Response to the Filaments and to the Satellite

mSSA operates by constructing a time-lagged covariance matrix of the input series, capturing autocorrelations within the ensemble of BFE coefficients. Each time series is detrended by subtracting its mean and normalizing by its variance to avoid
amplitude bias (M. D. Weinberg & M. S. Petersen 2021). We use a window length *L* corresponding to half the total time series (roughly 3.5 Gyr), to prioritize low-frequency signals tied to filamentary accretion and satellite orbital timescales. The SVD of this matrix yields PCs, ranked by their singular values, which quantify the contribution of each PC to the total variance in the decomposition. In our analysis, the first six PCs account for ~95% of the total power, with the first 11 capturing ~99%, while PCs corresponding to numerical noise remain uncorrelated and exhibit low singular values.

A key advantage of mSSA is its ability to approximate
signals with time-dependent frequencies. Signals at different
characteristic frequencies are projected into multiple PC
groups, depending on the duration over which a given
frequency persists. For instance, a persistent signal with stable
frequency (e.g., filamentary accretion) will dominate fewer
PCs, as its covariance structure remains coherent across the
entire time series. In contrast, transient or modulated signals
(e.g., satellite interactions) require more PCs to capture their
time-dependent frequency shifts and localized perturbations.
Consequently, a single dynamical process may span multiple
PCs. By analyzing the frequency content and correlation
structure of the PCs, we identify two distinct groups: one
associated with filament-driven quadrupole variations (fila-
ment frequency group) and another corresponding to the
orbital response to the LMC analog (LMC frequency group).
The detailed methodology for this classification is described in
Appendix C.

We reconstruct the BFE coefficients corresponding to the
filaments and the satellite using the PC corresponding to each
group. We then compute the contribution to the quadrupole
from the filaments and the LMC. This procedure allows for a
direct comparison of the disentangled filamentary and satellite-
induced signatures. For details on the reconstruction proce-
dure, see M. D. Weinberg & M. S. Petersen (2021) and
A. C. Johnson et al. (2023).

Figure 10 shows the time evolution of the quadrupole power relative to the monopole for two reconstructed components: the filament-associated (cyan) and LMC-analog-associated (yellow) quadrupole response. The total quadrupole power from the original BFE expansion is shown in black, with vertical dashed lines marking the first and second pericentric passages of the LMC analog. The dashed-dotted gray curve shows the quadrupolar power for the MW-LMC constrained simulations from N. Garavito-Camargo et al. (2021), multiplied by a factor of 10, with time shifted to align the LMC’s first pericentric passage with that of the m12b LMC analog.

The mSSA decomposition reveals two dominant trends in
the quadrupole evolution. The first (cyan) corresponds to a
long-term, gradually declining quadrupole mode driven by the
filamentary structure of the halo (Section 4). This component
has a characteristic timescale of ~4 Gyr (refer to Appendix C),
reflecting the slow evolution of the large-scale halo anisotropy.
The second component (yellow) exhibits a shorter-timescale
oscillatory pattern, with peaks in quadrupole power occurring
before the pericenter passages. This mode, associated with the
LMC analog, has a timescale of ~1.5 Gyr (refer to
Appendix C) and is approximately a quarter of the total
quadrupole power. By the present day, both modes become
comparable in magnitude, with the satellite-associated comp-
onent exhibiting slightly higher power, reflecting its growing
dynamical influence even after the merger. This transient

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Reconstructed gravitational power in the quadrupole harmonics relative to the monopole as a function of time, derived from the BFE coefficients decomposed using mSSA. The cyan and yellow curves represent the contributions from the filament-associated and LMC-analog-associated modes, respectively. The dashed-dotted gray line represents the quadrupolar power in the MW-LMC constrained simulation from N. Garavito-Camargo et al. (2021), multiplied by a factor of 10, with time shifted to align the LMC’s first pericentric passage with that of the m12b LMC analog. Vertical dashed lines indicate the first and second pericentric passages of the LMC analog. This decomposition isolates the background halo response from the dynamical perturbations induced by the satellite.</figcaption></figure>

response captures the halo's reaction to the satellite's torque
and the formation of a dynamical wake.

Notably, the amplitude of the satellite-induced quadrupole
in m12b (and in other FIRE-2 analogs) is an order of
magnitude higher than in MW-LMC constrained simulations
(dashed-dotted gray curve), which assume an initially
spherical MW halo and LMC progenitor (N. Garavito-Cama-
rgo et al. 2021; S. Lilleengen et al. 2023). In those models, the
development of a quadrupole mode is driven by the satellite’s
orbital energy transfer, as the system first has to establish a
nonspherical response. In contrast, our simulations feature a
preexisting quadrupole structure, allowing the analog’s torque
to act on an already anisotropic halo, resulting in a much
stronger response.

The mSSA's ability to isolate these contributions provides a
framework for precise examination of how distinct gravita-
tional effects of large-scale structure and satellite perturbations
shape the evolving halo. This will allow us to compare the
response across halos that are in different environments or
simulated with different DM particle properties.

5.2. Decomposing the Density Field Response to the Filaments
and to the Satellite

Once the PCs for the filaments and the satellite are
identified, one can compute the contribution of each PC to
the quadrupole coefficients to obtain separate time-series
expansions for the filaments and the satellite. Figures 11 and
12 illustrate the spatial distribution of the quadrupole density
field at different times in the m12b halo, reconstructed using
the two identified frequency groups from the mSSA analysis.
These fields are computed by reconstructing the BFE
coefficients using the PC groups associated with the filament
frequency (left two columns) and the LMC frequency (right
two columns). Each row represents a different snapshot in
time, with Figure 11 covering the interval from 5.8 to 8.8 Gyr,
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Time evolution of the quadrupole-environment view (columns (1) and (3)) described in Section 4.1, reconstructed for the two identified frequency groups: filament frequency (first two columns) and LMC frequency (last two columns) PCs in the m12balo. Each row spans from 5.8 to 8.8 Gyr. The first and third columns show the global view (within 1 Mpc), while the second and fourth columns focus on the local region (within 100 kpc). Open contours represent the quadrupole density field (vertical color bar), with the LMC-analog trajectory (dashed black curve and a magenta star for the current position) overlaid for reference. In the global view, the background highlights filamentary structures with the main halo and LMC analog subtracted. In the local view, the background shows the main-halo- and LMC-analog-subtracted surface densities. The reported values in the local regions are scaled down by a factor of 100 for clarity. The dominant frequency quadrupole components align with large-scale filamentary structures, while the subdominant components capture the localized impact of the satellite's passage. A full time-evolving video of this figure is available at the Harvard Dataverse: doi:10.7910/DVN/QWFJZP.</figcaption></figure>

while Figure 12 covers the interval from 9.8 to 12.8 Gyr of the simulation. The LMC analog's trajectory (dashed black curve and a magenta star for the current position) is overlaid for reference in all panels. A full video illustrating the temporal evolution of this structure is available at the Harvard Dataverse: doi:10.7910/DVN/QWFJZP.

In each panel, the first and third columns display a global view of the quadrupole field within 1 Mpc, while the second and fourth columns provide a zoomed view of the inner

100 kpc, centered on the main halo. The quadrupole density field, reconstructed from the dominant frequency components in each PC group, is represented by open contours and is color-coded according to the vertical color bar, which denotes the 3D density values. The background in the global view panels (columns (1) and (3)) highlights the filamentary structures, with the contribution from the main halo and LMC analog subtracted to highlight the filaments. In contrast, the localized view panels (columns (2) and (4)) display the surface density
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Same as Figure 11, but for the final 4 Gyr of evolution in m12b. A full time-evolving video of this figure is available at the Harvard Dataverse: doi:10.7910/DVN/QWFJZP.</figcaption></figure>

of the main halo with the LMC analog subtracted, shown with a horizontal color bar. The reported density values in the local view panels are scaled down by a factor of 100 for clarity.

In the filament frequency panels, the global view (first column) reveals a strong alignment between the quadrupole structure and the surrounding filamentary environment at all times, indicating that the primary response of the quadrupole is dictated by the cosmic web. This is consistent with previous studies showing that large-scale tidal fields influence the morphology and orientation of DM halos over cosmological timescales (O. Hahn et al. 2007; N. I. Libeskind et al. 2013; S. Codis et al. 2015, 2018). As the filaments broaden and their density diminishes over time, the power in the quadrupole component correspondingly decreases (see Section 4.1),

reflecting the gradual weakening of external anisotropic forces (C. A. Vera-Ciro et al. 2011). The localized filamentary view (second column) shows that the quadrupole structure in the inner halo remains largely unaffected by the presence of the LMC analog, maintaining its alignment with the outer halo and the large-scale environment.

In contrast, the quadrupole structure reconstructed from the LMC frequency components exhibits a distinct response. The global view (third column) indicates only slight deviations in the outer halo, with the quadrupole retaining a residual memory of its initial orientation. This supports the idea that the outer shape is predominantly set by early accretion and large-scale tidal interactions, persisting long after the filamentary structures have weakened (J. Bailin & M. Steinmetz 2005;
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Temporal evolution of the absolute value of the dot product (alignment) between the PC group quadrupole position angles $\hat{n}_{\ell=2,\text{group kind}}$ for the filament (cyan) and LMC (yellow) groups with the maximum-density filament orientation $\hat{n}_{\text{fila}}$ and the LMC-analog position $\hat{n}_{\text{LMC}}$, respectively. The gray curve shows the alignment between $\hat{n}_{\text{fila}}$ and $\hat{n}_{\text{LMC}}$. Vertical dashed lines mark the first and second pericentric passages of the LMC analog. The filament-driven quadrupole remains nearly perfectly aligned with the large-scale filaments over time, maintaining a dot product close to unity. The LMC-driven quadrupole rapidly aligns as the LMC analog approaches first pericenter, reaching near-unity alignment and maintaining it thereafter.</figcaption></figure>

B. Allgood et al. 2006). Meanwhile, the localized LMC frequency view (fourth column) undergoes substantial torquing and modification as the analog approaches the galactic center. The density amplitude of the quadrupole field increases by an order of magnitude during the analog’s passage as also seen in idealized tailored N-body simulations (M. D. Weinberg 1998). Once the merger is complete and the wake is fully phase-mixed, the quadrupole stabilizes, reaching a new equilibrium state in which the local structure no longer undergoes substantial changes. However, the final alignment of the quadrupole remains distinct from that set by the filamentary structure, illustrating a nontrivial interaction between these two dominant modes that operate on different spatial and temporal scales.

Figure 13 quantifies the temporal evolution of the quadrupole orientation decomposed into contributions from the filamentary and LMC PC groups. The alignment here is shown as the absolute value of the dot product between position angles of each component (similar to Figure 8). The cyan curve tracks the alignment of the filament-dominated quadrupole mode with the direction of the maximum-density filament, while the yellow curve tracks the alignment of the LMC-driven quadrupole mode with the position of the LMC analog. The gray curve shows the alignment between the maximum-density filament and LMC-analog position angles. Vertical dashed lines mark the first and second pericentric passages of the LMC analog.

The filament-driven quadrupole remains consistently aligned with the large-scale filamentary structure, maintaining a dot product close to unity. Interestingly, we observe a slight decrease in the filament-driven quadrupole alignment immediately following the LMC analog’s second pericenter, but overall it remains well aligned with a dot product greater than 0.9. In contrast, the LMC-driven quadrupole starts not well aligned at early times, it is seen to be more aligned with the direction of the filaments (the gray curve traces the yellow curve well), and as the LMC analog approaches first pericenter it undergoes a rapid increase in alignment, reaching near unity

and maintaining it thereafter. This suggests that the satellite’s gravitational influence reorients the local quadrupole early in its interaction, after which the alignment remains stable.

These findings demonstrate the importance of both filamentary accretion and merging satellites in shaping the structure of DM halos. The interplay between these mechanisms provides critical insights into the evolution of nonaxisymmetric structures and their observable imprints in the stellar and gaseous components of galaxies. Specifically, the persistence of these structures and their alignment with the cosmic web offer a unique window into the disequilibrium dynamics of stellar halos, as observed in their kinematic and morphological signatures (K. V. Johnston et al. 2008; A. Pillepich et al. 2014; M. Bennett et al. 2022).

In summary, we demonstrate the utility of mSSA in conjunction with BFE to decouple and analyze the time-evolving quadrupole structure of an MW-mass halo. Our decomposition identifies two dominant components: a long-term, slowly evolving quadrupole mode aligned with the filamentary large-scale structure, and a transient, satellite-induced mode associated with the LMC analog. The filament-driven quadrupole gradually weakens over time as the filaments fade away, while the satellite-induced quadrupole undergoes significant short-term amplification and torquing during the LMC analog’s orbital evolution. The ability to isolate these components provides a clearer framework for understanding how global tidal fields and local satellite interactions shape halo morphology, with implications for the persistence of nonaxisymmetric features in galaxies.

## 6. Summary and Discussion

Understanding the DM halo shape and evolution is critical for both galaxy formation theory and galactic dynamics. From a galaxy formation perspective, the structure of halos is shaped by filamentary accretion and merging satellites, which regulate how gas collapses to form galaxies, influencing star formation, angular momentum buildup, and feedback processes. Our results provide a direct measurement of how these large-scale structures imprint long-lived distortions in the halo, which in turn affect the evolution of galaxies embedded within them.

On the other hand, galactic dynamics relies on accurate halo models to reconstruct the gravitational potential. The assumption of smooth, static, and spherical halos—implicitly adopting a shared principal-axis frame with the stellar disk—is commonly used to model the orbits of stars, satellite galaxies, and stellar streams. However, our findings challenge this assumption by demonstrating that asymmetric perturbations are an inherent and persistent feature of halos, which may systematically bias dynamical inferences of DM distributions. Indeed, both observations (J. J. Han et al. 2022; J. Nibauer & A. Bonaca 2025) and simulations show that DM halos are frequently tilted by >20° relative to their stellar disks—as shown in the TNG (J. J. Han et al. 2023b) and FIRE-2 Latte suite (used in this work; J. Baptista et al. 2023)—further underscoring the need to go beyond aligned, spherical models when interpreting dynamical data.

In this work, we have decomposed zoomed cosmological MW-mass halos into time-evolving spherical harmonic expansions to quantify their large-scale gravitational distortions. Two modes capture nearly all large-scale halo distortions (see Figure 3). The dipole encodes large-scale displacements of the halo’s center of mass, while the quadrupole measures the halo’s global triaxial shape and its
slow evolution under filamentary accretion. Our key conclu-
sions are as follows:

1. *Filament-driven quadrupole anchors halo morphology.* The quadrupole is set early by Mpc-scale filaments and aligns with them owing to anisotropic accretion except during major mergers (see Figure 8). Even as filaments dissipate, the halo retains a “memory” of its orientation (see Figures 7 and 15)—up to 50% density distortion at 30 kpc and 1–2× near the virial boundary (see Figure 9). The quadrupole amplitude declines gradually over several Gyr (see Figure 5).

2. *Dipole and quadrupole amplitudes reach maxima near pericenter.* The host's dipole peaks sharply at each satellite pericenter—reaching fractional density perturbations up to unity (Figure 4)—and damps on orbital timescales (e.g., 1–2 Gyr), and this bulk motion must be removed to recover the true intrinsic halo triaxiality (quadrupole). The quadrupole exhibits localized boosts from dynamical friction wakes and direct torquing prior to the satellite's pericenter (Figure 5), superimposed on a secularly declining baseline quadrupole tied to the halo's shape and cosmic environment.

3. *Halo response at pericenter*, and this bulk motion must be removed to recover the true intrinsic quadrupole. The quadrupole then shows localized boosts from wakes and torques (Figure 5) on top of a long-lived baseline.

4. *The MW's present-day halo shape and response.* Since the LMC crossed its first pericenter, both the dipole and quadrupole amplitudes in our Galaxy's halo likely remain near their peak values—consistent with idealized *N*-body simulations of the MW–LMC system. Accurately modeling and subtracting this ongoing dipolar motion is therefore essential before one can robustly infer the halo's intrinsic shape.

5. *The satellite-induced response is boosted in cosmological halos.* In our zoomed cosmological halos, an LMC-mass satellite induces a quadrupole power at pericenter that is an order of magnitude larger than in constrained simulations (see Figure 10) that assume a spherically symmetric MW (N. Garavito-Camargo et al. 2021; S. Lilleengen et al. 2023). This satellite-induced quadrupole response is fundamentally sensitive to the true shape of the halo (see Figures 11 and 12).

6. *mSSA disentangles filamentary and satellite quadrupole.* mSSA on the quadrupole coefficients cleanly isolates the slow, filament-driven quadrupole from the fast, satellite-driven response (see Figures 10 and 13). Models that assume spherical halos fail to capture this superposition, systematically misrepresenting the dynamical evolution of triaxial MW-mass halos.

Accounting for the MW’s full assembly history and triaxiality —captured by our time-evolving quadrupole—has direct consequences for modeling our Galaxy’s observable structures. Triaxial halo asymmetries and past mergers can seed vertical disturbances in the disk, such as warps, corrugations, and bending waves (M. D. Weinberg 1998; M. D. Weinberg & L. Blitz 2006; F. A. Gómez et al. 2013, 2016; C. F. Laporte et al. 2018, 2019; J. A. Hunt et al. 2021; R. J. Grand et al. 2023; J. J. Han et al. 2023a). Notably, the larger satellite response in an asymmetric halo can exert a stronger gravitational torque on

the disk, potentially resolving the discrepancy in which constrained MW–LMC simulations yield disk warps that are too small—despite accurately reproducing their spatial locations (C. F. Laporte et al. 2018; I. A. Stelea et al. 2024).

In the stellar halo, observed LMC-induced wake amplitudes (C. Conroy et al. 2021; J. A. Amarante et al. 2024; M. Cavieres et al. 2025) exceed predictions from constrained DM-only spherical halo simulations (N. Garavito-Camargo et al. 2019). The triaxial halos observed in cosmological simulations enhance these wake amplitudes owing to stronger gravitational coupling between the satellite and halo, naturally explaining the observed excess. Moreover, halo tracers such as satellites and stellar stream orbits are widely used to constrain the Galactic potential and the LMC’s total mass (K. V. Johnston et al. 1999; E. Patel et al. 2018; D. Erkal et al. 2019; K. Malhan & R. A. Ibata 2019; N. Shipp et al. 2021; S. E. Koposov et al. 2023). While recent stream-based models incorporate the LMC’s time-dependent influence within a symmetric MW halo (D. Erkal et al. 2019; N. Shipp et al. 2021; S. Lilleengen et al. 2023), our work shows that neglecting the halo triaxiality introduces systematic offsets in inferred masses. Including this asymmetry (E. Vasiliev et al. 2021; J. Nibauer & A. Bonaca 2025) is essential for accurate orbit modeling and Galactic potential recovery.

The upcoming decade of deep photometric and kinematic data from Rubin/LSST (Ž. Ivezić et al. 2019) and the Roman Space Telescope (D. Spergel et al. 2015) will push well beyond Gaia's reach. These facilities will map halo tracers with precise proper motions and radial velocities out to—and beyond—the MW's virial radius. With this unprecedented data set, we will have the opportunity to directly test whether these asymmetries manifest in the orbits of distant halo tracers. By comparing kinematics to predictions from static versus dynamically evolving halo models, we can take the next step toward fully characterizing the morphology of the MW's DM halo and its response to cosmic accretion.

Acknowledgments

A.A. and R.E.S. were supported in part by Simons Foundation grant 1018462, NSF grants AST-2007232 and AST-2307787, and NASA grant 19-ATP19-0068. N.G.-C. and K.J.D. acknowledge support provided by the Heising-Simons Foundation grant No. 2022-3927. R.E.S. and M.D.W. were supported in part by grant No. NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP). R.E.S. was also supported in part by a Sloan Fellowship. M.S.P. is supported by a UKRI Stephen Hawking Fellowship. F.A.G. acknowledges support from the FONDECYT Regular grant 1211370, by the ANID BASAL project FB210003, and from the HORIZON-MSCA-2021-SE-01 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement No. 101086388. K.V.J. was supported by Simons Foundation grant 1018465. C.F.P.L. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 852839). J.A.S.H. acknowledges the support of a UKRI Ernest Rutherford Fellowship ST/Z510245/1. G.B. is supported by NSF CAREER award AST-1941096 and NASA ATP award 80NSSC24K1225. The EXP collaboration gratefully acknowledges the support of the Simons Foundation.

A.A. and N.S. acknowledge support from the DiRAC Institute in the Department of Astronomy and the eScience
Institute, both at the University of Washington. N.G.-C. and K. J.D. respectfully acknowledge that the University of Arizona is on the land and territories of Indigenous peoples. Today, Arizona is home to 22 federally recognized tribes, with Tucson being home to the O'odham and the Yaqui. The University strives to build sustainable relationships with sovereign Native Nations and Indigenous communities through education offerings, partnerships, and community service.

A.A. would like to thank the FIRE collaboration, the anonymous referee, and the Scientific Computing Team at the Flatiron Institute for access to data, computing resources, and valuable discussions that shaped this work. A.A. also expresses sincere gratitude to ChatGPT for its language assistance (and for writing its own acknowledgment).

*Software:* IPython (F. Perez & B. E. Granger 2007), Numpy (C. R. Harris et al. 2020), Scipy (P. Virtanen et al. 2020), Gizmo Analysis (A. Wetzel & S. Garrison-Kimmel 2020a), pyEXP (M. S. Petersen et al. 2022b; M. S. Petersen & M. D. Weinberg 2025), Rockstar (P. S. Behroozi et al. 2012a), Halo Analysis (A. Wetzel & S. Garrison-Kimmel 2020b), Matplotlib (J. D. Hunter 2007), CMasher (E. van der Velden 2020).

## Appendix A
### Temporal Stability of Radial Basis Functions

In this appendix, we demonstrate that the choice of a time-averaged density profile for constructing the radial basis functions is well justified. Specifically, we show that the radial polynomials computed at each time step remain largely unchanged over the course of the simulation, reinforcing the validity of using a fixed, time-averaged basis.

A full mathematical description of the construction of adaptive basis functions and the solution to the Sturm-Liouville equation can be found in M. D. Weinberg (1999). In brief, these basis functions are derived by numerically solving the Poisson equation for an orthogonal potential-density pair,

using the empirically derived spherically averaged density profile of the halo at each time step. This procedure ensures that the basis functions capture the structure of the gravitational potential in a dynamically consistent manner.

To assess the stability of the basis, we compute the radial basis functions at individual snapshots and compare them over time. Figure 14 illustrates the resulting radial basis functions for the monopole (top row) and quadrupole (bottom row), computed from 5.8 to 13.8 Gyr at 2 Gyr intervals (columns). The functions are plotted as normalized potentials, color-coded by radial order (n), with redder shades representing higher-order polynomials.

We observe that the overall shapes and amplitudes of the polynomials remain remarkably stable over time. While minor variations are present in the inner regions owing to the deepening of the potential with ongoing accretion, these changes occur gradually on 2–4 Gyr timescales and do not significantly alter the global structure of the basis. This stability supports our approach of using a mean density profile to construct a fixed set of radial basis functions, as the basis functions computed at any given snapshot are highly similar to those derived from the time-averaged profile (see Figure 1).

An alternative approach to mitigate fluctuations in the empirical density profile is to fit the profile with a smooth double power-law function and derive the basis functions from this fitted model. This method further stabilizes the basis, particularly in cases where transient features introduce short-lived variations in the empirical density profile. However, we find that, even without such smoothing, the empirical density profile itself produces basis functions that are highly stable over time.

Across all our simulated halos, we find that the basis remains stable, even in the presence of minor fluctuations. However, during pericentric passages of massive satellites—such as in m12w—the basis exhibits transient deviations. Even

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Temporal evolution of orthogonal radial basis functions for m12i, computed by solving the Poisson equation using the empirical density profile at each time step. Columns represent different epochs from 5.8 to 13.8 Gyr, while rows correspond to the monopole (top) and quadrupole (bottom) radial functions. The basis functions are shown as normalized potentials as a function of radial distance, color-coded by radial order (n), with redder shades indicating higher-order polynomials. The overall shapes and amplitudes of the basis functions remain stable over time, with minor variations in the inner regions due to the deepening of the potential. These results justify the use of a time-averaged spherical density profile to construct a fixed basis for the expansion, ensuring consistency in time-series analysis of the expansion coefficients.</figcaption></figure>
in these cases, the overall structure of the basis is retained outside of these brief perturbations.

Thus, using a time-averaged density profile to derive a fixed basis is a well-motivated and robust choice. It ensures that the expansion coefficients are evaluated consistently over time and can be analyzed through mSSA to reconstruct physically meaningful gravitational structures.

**Appendix B**

**Quadrupole-environment Evolution for All Simulations**

Figure 15 presents the time evolution of the quadrupole-environment view for all simulations, complementing the examples shown in Section 4.1 and Figure 6 (view description) and the detailed evolution highlighted for m12b in Figure 7. Each row represents a different simulation, while columns

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Time evolution of the quadrupole-environment view (as described in Figure 6) for the past 7 Gyr of the halo evolution. Each row corresponds to a different simulation, and each column shows a snapshot at 2 Gyr intervals, from 6.8 to 12.8 Gyr. The quadrupole density strength, $\rho_{\ell=2}$, is represented by red-blue open contours, indicating positive and negative regions, respectively, computed using the BFE formalism. The background surface density (diverging color map) highlights filamentary structures within 1 Mpc, with the main halo subtracted. Green arrows indicate the orientation of the quadrupole moment, $\hat{n}_{\ell=2}$, while black arrows denote the dominant filament orientation, $\hat{n}_{\text{fila}}$. Dashed black lines trace the trajectories of massive infalling satellites, with their current positions marked by open magenta stars. The fractional quadrupole gravitational power is marked in lower left corner in each panel. The time evolution illustrates the strong correlation between the halo quadrupole and its filamentary environment, as well as the influence of massive satellites.</figcaption></figure>
show snapshots at 2 Gyr intervals over the past 7 Gyr of evolution. The corresponding fractional gravitational power in the quadrupole is marked in the lower left corner.

We note that the simulations exhibit a wide variety of filamentary environments. For example, m12i is influenced by two dominant filaments, while m12r features three to four distinct anisotropic filaments feeding the halo. In contrast, m12m is embedded in a massive, broad filament that affects the halo nearly isotropically, whereas m12c is surrounded by a single dominant filament. Despite these differences, the amplitude of the quadrupole moment decreases consistently as the filamentary structures broaden and their surface densities decrease over time.

The evolution of the quadrupole moment varies depending on the filamentary configuration. In m12m, the gradual fading of the surrounding broad filament corresponds to a significant reduction in the quadrupole density amplitude. In m12r, the persistent presence of multiple distinct filaments maintains a stronger and more dynamic quadrupole structure over time. These trends emphasize that the halo's quadrupole responds directly to changes in its filamentary environment, with alignment and amplitude varying based on the filaments' density distributions and spatial configurations.

Interactions with infalling satellites also contribute to variations in the quadrupole, particularly near their pericentric passages. For example, m12b and m12r show temporary deviations in quadrupole alignment and amplitude during such events. However, these effects are secondary compared to the dominant influence of the filamentary environment; the halo tends to retain the memory of its original quadrupole orientation (e.g., J. Dubinski & R. Carlberg 1991; J. Bailin & M. Steinmetz 2005; B. Allgood et al. 2006; C. A. Vera-Ciro et al. 2011) even after satellite perturbations. Notably, massive satellites predominantly infall along major filaments (e.g., G. Tormen et al. 1997; D. Aubert et al. 2004; N. I. Libeskind et al. 2014; E. Tempel et al. 2015), further reinforcing the connection between these structures and the halo's evolution.

These results confirm and extend the findings in Section 4.1, demonstrating that the primary driver of halo shape and evolution is the global filamentary structure feeding the halo. While satellite interactions introduce localized perturbations, the filamentary environment predominantly determines the long-term evolution of the quadrupole and, by extension, the triaxiality of the halo.

# Appendix C
## Finding mSSA Principal Component Groupings

In this appendix, we detail the methodology used to group PC from the SVD of the lagged trajectory matrix in mSSA. Since the input coefficient time series are often interdependent, a single dynamical signature can be spread across multiple PCs. Additionally, varying frequency modes may be distributed among different PCs, necessitating a systematic approach to identifying physically meaningful groupings.

The number of PCs required to describe a dynamical signal is tied to its temporal complexity. For instance, linear or exponential trends (e.g., gradual filamentary growth) are captured by a single PC, as their evolution depends only on the immediate past. Periodic signals (e.g., satellite orbital motion), however, require groups of PCs to encode amplitude and phase, analogous to sine and cosine components. Nonstationary signals—such as those involving abrupt

changes (e.g., pericenter passages) or frequency modulation —demand additional PCs to account for their evolving structure.

One of the most fundamental yet tedious methods for grouping PCs involves reconstructing gravitational fields from individual PCs and visually inspecting their spatial and temporal behavior. However, a more quantitative grouping can be achieved by applying the following criteria:

1. *Singular Values*. PCs associated with the same dynamical process tend to have comparable singular values.

2. *Time Evolution*. PCs that exhibit similar time-dependent behavior can be identified by plotting their amplitudes over time (top panels of Figure 16).

3. *Frequency Analysis*. Applying a discrete Fourier transform (DFT) to the PCs allows us to examine their dominant frequency components and identify those that share similar power spectra (bottom left panel of Figure 16).

4. *w-Correlation Matrix*. The weighted correlation matrix, as introduced in M. D. Weinberg & M. S. Petersen (2021), quantifies the correlation between different PCs, helping to establish natural groupings (bottom right panel of Figure 16).

For the case of m12b, we analyze the quadrupole ($l = 2$) coefficients, including all azimuthal modes ($m = 0, 1, 2$) and the full range of radial orders ($n$). We focus on the first 12 PCs, as they collectively account for 98% of the total variance, with the singular values dropping exponentially beyond this set. In our analysis, the filament-driven quadrupole aligns with a quasi-linear trend, requiring few PCs, while the satellite-induced signal spans multiple PCs owing to its quasi-periodic orbital motion and transient interactions with the halo.

Figure 16 illustrates the process used to define PC groupings. The top row presents the PC amplitudes (color-coded) over time. The left panel shows the four most dominant PCs (PC0–PC3), while the right panel displays PC4–PC11. These groupings are made based on their temporal signatures and amplitudes. The bottom row provides further insight: the left panel shows the DFT power spectrum for each PC, revealing characteristic timescales, while the right panel displays the w-correlation matrix, quantifying the correlation between different PCs. From this analysis, we identify two primary groupings:

1. *Filamentary structure group*. PC0–PC3 exhibit low-frequency variations with characteristic timescales of ~4 Gyr, consistent with the slow evolution of the large-scale halo anisotropy. These components also share power at a secondary ~3 Gyr timescale. The w-correlation matrix confirms their mutual correlation.

2. *LMC-analog response group*. PC4–PC11 contain power at higher frequencies, peaking at ~1.5 Gyr, indicative of the satellite-induced perturbations. These PCs show moderate correlations in the w-correlation matrix.

3. *Uncorrelated component*. As an example, PC12 exhibits three-orders-of-magnitude-lower power at very high frequencies and appears uncorrelated with other PCs, suggesting that it primarily represents noise or a higher-order transient effect.

In Section 5, we utilize these identified filamentary and LMC-associated frequency groups to reconstruct and analyze the gravitational structures driving halo evolution.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** PC grouping for the quadrupole time series in m12b, based on mSSA decomposition. Top row: PC amplitudes as a function of time (half-window length shown). The left panel presents PC0–PC3, which exhibit similar long-timescale evolution, while the right panel shows PC4–PC11, dominated by higher-frequency variations. Bottom row: the left panel displays the DFT power spectrum for each PC, highlighting characteristic timescales; the right panel shows the w-correlation matrix, quantifying the correlation between different PCs. These analyses define two primary groupings—one associated with filamentary structure (PC0–PC3) and another with the LMC-analog response (PC4–PC11). PC12, which exhibits low correlation and significantly lower power, remains ungrouped.</figcaption></figure>

## ORCID iDs

Arpit Arora https://orcid.org/0000-0002-8354-7356

Nicolás Garavito-Camargo https://orcid.org/0000-0001-7107-1744

Robyn E. Sanderson https://orcid.org/0000-0003-3939-3297

Martin D. Weinberg https://orcid.org/0000-0003-2660-2889

Michael S. Petersen https://orcid.org/0000-0003-1517-3935

Silvio Varela-Lavin https://orcid.org/0000-0003-3213-8736

Facundo A. Gómez https://orcid.org/0000-0002-1947-333X

Kathryn V. Johnston https://orcid.org/0000-0001-6244-6727

Chervin F. P. Laporte https://orcid.org/0000-0003-3922-7336

Nora Shipp https://orcid.org/0000-0003-2497-091X

Jason A. S. Hunt https://orcid.org/0000-0001-8917-1532

Gurtina Besla https://orcid.org/0000-0003-0715-2173

Elise Darragh-Ford https://orcid.org/0000-0002-8800-5652

Nondh Panithanpaisal https://orcid.org/0000-0001-5214-8822

Kathryne J. Daniel https://orcid.org/0000-0003-2594-8052

Arora, A., Garavito-Camargo, N., Sanderson, R. E., et al. 2024a, ApJ, 974, 286

Arora, A., Sanderson, R., Regan, C., et al. 2024b, ApJ, 977, 23

Arora, A., Sanderson, R. E., Panithanpaisal, N., et al. 2022, ApJ, 939, 2

Arora, A., & Williams, L. L. 2020, ApJ, 893, 53

Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376

Bailin, J., & Steinmetz, M. 2005, ApJ, 627, 647

Baptista, J., Sanderson, R., Huber, D., et al. 2023, ApJ, 958, 44

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2012a, ApJ, 762, 109

Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2012b, ApJ, 763, 18

Bennett, M., Bovy, J., & Hunt, J. A. 2022, ApJ, 927, 131

Besla, G., Kallivayalil, N., Hernquist, L., et al. 2007, ApJ, 668, 949

Besla, G., Kallivayalil, N., Hernquist, L., et al. 2010, ApJL, 721, L97

Bett, P., Eke, V., Frenk, C. S., et al. 2007, MNRAS, 376, 215

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton, NJ: Princeton Univ. Press)

Blumenthal, G. R., Faber, S., Primack, J. R., & Rees, M. J. 1984, Natur, 311, 517

Bonamigo, M., Despali, G., Limousin, M., et al. 2015, MNRAS, 449, 3171

Bovy, J., Bahmanyar, A., Fritz, T. K., & Kallivayalil, N. 2016, ApJ, 833, 31

Bryan, S., Kay, S., Duffy, A., et al. 2013, MNRAS, 429, 3316

Buch, D., Nadler, E. O., Wechsler, R. H., & Mao, Y.-Y. 2024, ApJ, 971, 79

Byström, A., Koposov, S. E., Lilleengen, S., et al. 2024, arXiv:2410.09149

Cautun, M., Benitez-Llambay, A., Deason, A. J., et al. 2020, MNRAS, 494, 4291

Cavieres, M., Chanamé, J., Navarrete, C., et al. 2025, ApJ, 983, 83

Chandra, V., Naidu, R. P., Conroy, C., et al. 2024, arXiv:2406.01676

## References

Allgood, B., Flores, R. A., Primack, J. R., et al. 2006, MNRAS, 367, 1781

Amarante, J. A., Koposov, S. E., & Laporte, C. F. 2024, A&A, 690, A166
Choi, J.-H., Weinberg, M. D., & Katz, N. 2009, MNRAS, 400, 1247

Clutton-Brock, M. 1972, Ap&SS, 16, 101

Codis, S., Jindal, A., Chisari, N., et al. 2018, MNRAS, 481, 4753

Codis, S., Pichon, C., & Pogosyan, D. 2015, MNRAS, 452, 3369

Conroy, C., Naidu, R. P., Garavito-Camargo, N., et al. 2021, Natur, 592, 534

Cooper, A. P., Cole, S., Frenk, C., et al. 2010, MNRAS, 406, 744

Deason, A. J., Belokurov, V., & Sanders, J. L. 2019, MNRAS, 490, 3426

Debattista, V. P., Roškar, R., Valluri, M., et al. 2013, MNRAS, 434, 2971

Despali, G., Giocoli, C., & Tormen, G. 2014, MNRAS, 443, 3208

Diemer, B., Behroozi, P., & Mansfield, P. 2024, MNRAS, 533, 3811

D'Onghia, E., Springel, V., Hernquist, L., & Keres, D. 2010, ApJ, 709, 1138

Doroshkevich, A. 1970, Afz, 6, 581

Dubinski, J., & Carlberg, R. 1991, ApJ, 378, 496

D'Souza, R., & Bell, E. F. 2022, MNRAS, 512, 739

Eisenstein, D. J., & Loeb, A. 1995, ApJ, 439, 520

Erkal, D., Belokurov, V., Laporte, C., et al. 2019, MNRAS, 487, 2685

Gaia Collaboration, Brown, A., Vallenari, A., et al. 2018, A&A, 616, A1

Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2019, ApJ, 884, 51

Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2021, ApJ, 919, 109

Garavito-Camargo, N., Price-Whelan, A. M., Samuel, J., et al. 2024, ApJ, 975, 100

Garrison-Kimmel, S., Hopkins, P. F., Wetzel, A., et al. 2018, MNRAS, 481, 4133

Geha, M., Wechsler, R. H., Mao, Y.-Y., et al. 2017, ApJ, 847, 4

Golyandina, N., Nekrutkin, V., & Zhigljavsky, A. A. 2001, Analysis of Time Series Structure: SSA and Related Techniques (Boca Raton, FL: CRC press)

Gómez, F. A., Besla, G., Carpintero, D. D., et al. 2015, ApJ, 802, 128

Gómez, F. A., Grand, R. J., Monachesi, A., et al. 2017, MNRAS, 472, 3722

Gómez, F. A., Minchev, I., O'Shea, B. W., et al. 2013, MNRAS, 429, 159

Gómez, F. A., White, S. D., Grand, R. J., et al. 2016, MNRAS, 465, 3446

Grand, R. J., Pakmor, R., Fragkoudi, F., et al. 2023, MNRAS, 524, 801

Hahn, O., Porciani, C., Carollo, C. M., & Dekel, A. 2007, MNRAS, 375, 489

Han, J. J., Conroy, C., & Hernquist, L. 2023a, NatAs, 7, 1481

Han, J. J., Naidu, R. P., Conroy, C., et al. 2022, ApJ, 934, 14

Han, J. J., Semenov, V., Conroy, C., & Hernquist, L. 2023b, ApJL, 957, L24

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375

Hopkins, P. F. 2015, MNRAS, 450, 53

Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800

Hoyle, F. 1953, ApJ, 118, 513

Hunt, J. A., Stele, I. A., Johnston, K. V., et al. 2021, MNRAS, 508, 1459

Hunt, J. A., & Vasiliev, E. 2025, NewAR, 100, 101721

Hunter, J. D. 2007, CSE, 9, 90

Hüttig, C., & Stemmer, K. 2008, GGG, 9, 137

Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111

Johnson, A. C., Petersen, M. S., Johnston, K. V., & Weinberg, M. D. 2023, MNRAS, 521, 1757

Johnston, K. V., Bullock, J. S., Sharma, S., et al. 2008, ApJ, 689, 936

Johnston, K. V., Spergel, D. N., & Hernquist, L. 1995, ApJ, 451, 598

Johnston, K. V., Zhao, H., Spergel, D. N., & Hernquist, L. 1999, ApJL, 512, L109

Kazantzidis, S., Kravtsov, A. V., Zentner, A. R., et al. 2004, ApJL, 611, L73

Koposov, S. E., Erkal, D., Li, T. S., et al. 2023, MNRAS, 521, 4936

Laporte, C. F., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286

Laporte, C. F., Johnston, K. V., & Tzanidakis, A. 2019, MNRAS, 483, 1427

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Libeskind, N. I., Hoffman, Y., Forero-Romero, J., et al. 2013, MNRAS, 428, 2489

Libeskind, N. I., Knebe, A., Hoffman, Y., & Gottlöber, S. 2014, MNRAS, 443, 1274

Lilleengen, S., Petersen, M. S., Erkal, D., et al. 2023, MNRAS, 518, 774

Lowing, B., Jenkins, A., Eke, V., & Frenk, C. 2011, MNRAS, 416, 2697

Malhan, K., & Ibata, R. A. 2019, MNRAS, 486, 2995

Mansfield, P., Darragh-Ford, E., Wang, Y., et al. 2024, ApJ, 970, 178

Nibauer, J., & Bonaca, A. 2025, ApJL, 985, L22

Patel, E., Besla, G., Mandel, K., & Sohn, S. T. 2018, ApJ, 857, 78

Peebles, P. J. 1969, ApJ, 155, 393

Perez, F., & Granger, B. E. 2007, CSE, 9, 21

Petersen, M. S., & Peñarrubia, J. 2020, MNRAS: Letters, 494, L11

Petersen, M. S., & Peñarrubia, J. 2021, NatAs, 5, 251

Petersen, M. S., Peñarrubia, J., & Jones, E. 2022a, MNRAS, 514, 1266

Petersen, M. S., & Weinberg, M. D. 2025, JOSS, 10, 7302

Petersen, M. S., Weinberg, M. D., & Katz, N. 2022b, MNRAS, 510, 6201

Pillepich, A., Vogelsberger, M., Deason, A., et al. 2014, MNRAS, 444, 237

Planck Collaboration, Ade, P., Aghanim, N., et al. 2016, A&A, 594, A13

Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14

Quinn, T., & Binney, J. 1992, MNRAS, 255, 729

Reino, S., Sanderson, R. E., Panithanpaisal, N., Rossi, E. M., & Kuijken, K. 2022, MNRAS, 509, 5365

Robotham, A., Baldry, I. K., Bland-Hawthorn, J., et al. 2012, MNRAS, 424, 1448

Salem, M., Besla, G., Bryan, G., et al. 2015, ApJ, 815, 77

Samuel, J., Wetzel, A., Chapman, S., et al. 2021, MNRAS, 504, 1379

Sanders, J. L., Lilley, E. J., Vasiliev, E., Evans, N. W., & Erkal, D. 2020, MNRAS, 499, 4793

Sanderson, R. E., Wetzel, A., Loebman, S., et al. 2020, ApJS, 246, 6

Santistevan, I. B., Wetzel, A., Tollerud, E., et al. 2024, MNRAS, 527, 8841

Schneider, M. D., Frenk, C. S., & Cole, S. 2012, JCAP, 2012, 030

Shipp, N., Erkal, D., Drlica-Wagner, A., et al. 2021, ApJ, 923, 149

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv:1503.03757

Springel, V., White, S. D., Jenkins, A., et al. 2005, Natur, 435, 629

Stelea, I. A., Hunt, J. A., & Johnston, K. V. 2024, ApJ, 977, 252

Tempel, E., Guo, Q., Kipper, R., & Libeskind, N. 2015, MNRAS, 450, 2727

Tormen, G., Bouchet, F. R., & White, S. D. 1997, MNRAS, 286, 865

Valluri, M., Price-Whelan, A. M., & Snyder, S. J. 2021, ApJ, 910, 150

van der Velden, E. 2020, JOSS, 5, 204

Vasiliev, E. 2013, MNRAS, 434, 3174

Vasiliev, E. 2023, Galax, 11, 59

Vasiliev, E. 2024, MNRAS, 527, 437

Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279

Vera-Ciro, C., & Helmi, A. 2013, ApJL, 773, L4

Vera-Ciro, C. A., Sales, L. V., Helmi, A., et al. 2011, MNRAS, 416, 1377

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Weinberg, M. D. 1998, MNRAS, 297, 101

Weinberg, M. D. 1999, AJ, 117, 629

Weinberg, M. D. 2023, MNRAS, 525, 4962

Weinberg, M. D., & Blitz, L. 2006, ApJL, 641, L33

Weinberg, M. D., & Petersen, M. S. 2021, MNRAS, 501, 5408

Wetzel, A., & Garrison-Kimmel, S., 2020a GizmoAnalysis: Read and Analyze Gizmo Simulations, Astrophysics Source Code Library, ascl:2002.015

Wetzel, A., & Garrison-Kimmel, S., 2020b HaloAnalysis: Read and Analyze Halo Catalogs and Merger Trees, Astrophysics Source Code Library, ascl:2002.014

Wetzel, A., Hayward, C. C., Sanderson, R. E., et al. 2023, ApJS, 265, 44

Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJL, 827, L23

White, S. D. 1984, ApJ, 286, 38

White, S. D., & Rees, M. J. 1978, MNRAS, 183, 341

Yaaqib, R., Petersen, M. S., & Peñarrubia, J. 2024, MNRAS, 531, 3524