Orientations of Dark Matter Halos in FIRE-2 Milky Way–mass Galaxies

The shape and orientation of dark matter (DM) halos are sensitive to the microphysics of the DM particles, yet in many mass models, the symmetry axes of the Milky Way’s DM halo are often assumed to be aligned with the symmetry axes of the stellar disk. This is well motivated for the inner DM halo, but not for the outer halo. We use zoomed-in cosmological baryonic simulations from the Latte suite of FIRE-2 Milky Way–mass galaxies to explore the evolution of the DM halo’s orientation with radius and time, with or without a major merger with a Large Magellanic Cloud analog, and when varying the DM model. In three of the four cold DM halos we examine, the orientation of the halo minor axis diverges from the stellar disk vector by more than 20° beyond about 30 galactocentric kpc, reaching a maximum of 30°–90°, depending on the individual halo’s formation history. In identical simulations using a model of self-interacting DM with σ = 1 cm2 g−1, the halo remains aligned with the stellar disk out to ∼200–400 kpc. Interactions with massive satellites (M ≳ 4 × 1010 M ⊙ at pericenter; M ≳ 3.3 × 1010 M ⊙ at infall) affect the orientation of the halo significantly, aligning the halo’s major axis with the satellite galaxy from the disk to the virial radius. The relative orientation of the halo and disk beyond 30 kpc is a potential diagnostic of self-interacting DM, if the effects of massive satellites can be accounted for.


Introduction
Dark matter (DM) halos have the potential to serve as macroscopic laboratories that can constrain the microphysics of a DM particle (e.g., Bullock & Boylan-Kolchin 2017, and references therein).Like all galaxies, the Milky Way (MW) is embedded within a DM halo, but despite our close-up view, its structure is not well constrained outside the inner 20 kpc.Especially difficult to constrain are the approximate principal axis ratios of the halo and their orientation relative to the Galactic disk (Vera-Ciro et al. 2011;Hattori et al. 2021).In external galaxies, the geometry of the DM halo is often determined by analyzing weak lensing, polar rings, and the distribution of satellite galaxies.Within our galaxy, methods for constraining the geometry of the DM halo rely on modeling the dynamics of distant tracers, such as standard candle stars, tidal streams, globular clusters, or dwarf galaxies, and are thus limited by the depths of current surveys (Vera-Ciro et al. 2011;Sanderson et al. 2019;Hattori et al. 2021;Gaia Collaboration et al. 2021;Reino et al. 2021).
With the advent of new dynamical information about these tracers, from missions such as Gaia, JWST, Nancy Grace Roman Space Telescope, and Vera C. Rubin Observatory, we will be able to detect stellar tracers out to the edge of the MW's DM halo (Sanderson et al. 2017(Sanderson et al. , 2019) ) and resolve the stars in the faint outskirts of up to 100 MW-like galaxies nearby (Pearson et al. 2022).Therefore, now is the time to examine the geometry of DM halo outskirts in realistic simulations of MW analogs in order to make predictions for the shape and orientation of the halo beyond 20 kpc (e.g., Prada et al. 2019), as well as to investigate its sensitivity to the physics of the DM particle (e.g., Vargya et al. 2022).
Theoretical studies of the expected DM halo shape have been done using simulations of MW-like galaxies both with and without the incorporation of baryonic physics (Vera-Ciro et al. 2011;Prada et al. 2019) and by varying the DM model (Vargya et al. 2022).However, somewhat less attention has been paid to the orientation of the symmetry axes of the DM halo as a function of radius and time, especially in cosmological baryonic simulations where central galaxy formation has been shown to affect DM halo shape (Bailin & Steinmetz 2004;Allgood et al. 2006;Garrison-Kimmel et al. 2018;Prada et al. 2019).The orientations of the DM halo axes may play an important role in the kinematic twisting of the stellar components in the outer stellar halo of observed Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.galaxies, since the potential due to DM mainly drives stellar kinematics at those radii (Kormendy et al. 2009;Foster et al. 2016;Ene et al. 2018;Huang et al. 2018;Pulsoni et al. 2018).
Furthermore, the ongoing interaction of the MW with the Large Magellanic Cloud (LMC) is expected to influence the orientation of the DM symmetry axes, especially at intermediate radii (Garavito-Camargo et al. 2021;Vasiliev et al. 2021), yet studies of the geometry of halos in cosmological MW analogs to date do not include LMC-like interactions.The trajectory of the LMC in a twisted DM halo may also influence the dynamical evolution of recently discovered dwarf galaxies like the Antlia 2 dwarf satellite (Chakrabarti et al. 2019;Ji et al. 2021) and their impact on the gas disk (Levine et al. 2006) of the MW.
Previous studies using DM-only simulations (e.g., Allgood et al. 2006;Vera-Ciro et al. 2011) indicate that the symmetry axes of halos can twist substantially from the interior to the virial radius.However, models of the MW potential nearly always incorporate the assumption that the DM symmetry axes are aligned to the central galaxy at all radii, which can end up producing inconsistent results (Law et al. 2009;Debattista et al. 2013;Prada et al. 2019;Vasiliev et al. 2021).Although the halo and stellar disk axes are likely to be well aligned where stars and gas dominate the mass distribution, which includes the inner 20 kpc or so where the MW shape is currently well constrained, there is no reason to assume the halo axes should align to the stellar disk at larger distances or out to the virial radius.Bailin et al. (2005) find evidence of decoupling between the disk axes and the halo axes in hydrodynamical simulations beyond ∼0.1R vir .An exploration by Emami et al. (2021a) into MW-mass halos in Illustris TNG50 finds that there is a distinct population of galaxies with DM radial profiles that rotate more than 50°relative to a fixed angular momentum vector.From a sample of 25 galaxies, they find that 32% of their sample experienced such a twist (Emami et al. 2021a).Analyses of the MW's potential that have been modeled using the assumption of an aligned disk and halo (Cunningham et al. 2020) may thus be flawed if the outer halo is tilted significantly relative to the disk.
Additionally, there is observational evidence to suggest the presence of a twisted halo and its connection to merging satellite galaxies.Current observations of the Sagittarius Stream (a powerful tracer of the MW potential) suggest the MW disk is oblique to the outer DM halo (Debattista et al. 2013;Hattori et al. 2021;Vasiliev et al. 2021).Han et al. (2022) use the phase space of the Gaia-Sausage-Enceladus merger remnants in the stellar halo (out to 30 kpc) to fit the shape and orientation of the DM halo with an ellipsoid model.Their findings indicate the halo twists to about 30°above the disk plane, with the tilt pointing to the apocentric positions of the merger debris.This motivates not only the observed presence of a tilt in the MW DM halo, but also the role that tidal effects from massive satellites play in setting the orientation of the DM halo.Our work explores this latter scenario by taking into account how a massive gas-rich satellite (similar to the LMC) influences the symmetry axes of the DM halo.Similarly, in the Auriga simulations, simulated halos can become oblique to the symmetry axes of the central galaxy as a function of host distance and time (Prada et al. 2019).
However, it is reasonable to expect that the radial extent of the alignment between galaxies and their host halos would be somewhat dependent on the nature and strength of the stellar and active galactic nucleus feedback in the galaxy, which can couple small-scale regions in the central galaxy to much larger scales.It is therefore useful to study the obliquity of DM halos in simulations with a variety of feedback prescriptions, in the interest of isolating these effects.
In this paper, we use the Latte suite of 13 simulated galaxies (M * ∼ 5 × 10 10 M e , M halo ∼ 10 12 M e ; Wetzel et al. 2016), which uses substantially different feedback physics than the simulations mentioned above13 (Hopkins et al. 2018).We also investigate how an infalling LMC-mass galaxy can dynamically perturb the orientation of the DM halo by analyzing halo orientations throughout several LMC-like mergers with varying masses.Finally, we examine the differences in the relative orientations of the halo and disk when the same halo is simulated with elastic-scattering self-interacting DM (SIDM) from the same initial conditions.The SIDM model implemented is modeled as elastic hard-sphere scattering with a fixed cross section per unit mass, as compared to alternative models that include dissipative SIDM (dSIDM) or excited-state SIDM.We perform this analysis to determine whether next-generation surveys could distinguish between cold DM (CDM) and SIDM by constraining the orientation of the DM halo.
This paper is organized as follows.In Section 2, we discuss the suite of MW-mass simulations we analyze in this paper and the evolutionary histories of each simulation.In Section 2.1, we present a method used to define the symmetry axes of the DM halo and the stellar disk.In Section 2.2, we introduce a procedure to quantify the orientation between the disk and halo based on the halo and disk symmetry axes.In Section 2.3, we introduce our SIDM and DM-only (DMO) halos and outline procedural differences for calculating the orientations in those simulations.Section 2.4 outlines the process of calculating the triaxiality (shape) of the halo over time as a supplementary measurement to quantify merger effects.In Section 3, we introduce the halo-disk orientations at the present day (Section 3.1) and as a function of redshift (Section 3.2).Additionally, we explore the effects of LMC-mass satellites on orientation (Section 3.4).In Section 3.3, we perform our analysis of SIDM and DMO halo-disk orientations at the present day and present the differences in results between the simulations.In Section 3.5, we discuss the orientation of the halo short axis relative to the enclosing DM filament across our different simulations.We conclude our results and discussion in Section 4.

Methods
Our analysis makes use of the Latte suite of FIRE-2 cosmological zoom-in simulations (Wetzel et al. 2016;Hopkins et al. 2018).These simulations capture the dynamical feedback between the live baryonic and DM components of MW-mass galaxies as they emerge from the cosmic web.The Latte suite consists of MW-mass galaxies (∼10 12 M e ) with different initial conditions and merger histories.We use the fiducial-resolution Latte simulations with a resolution of 7100 M e per initialized gas particle and 3.5 × 10 4 M e per DM particle.
The FIRE simulations offer several important advantages in the study of halo orientations.Each galaxy is initialized with different initial conditions selected from a large cosmological box, which provides a diversity of halo assembly histories.The simulations have intervals of 20-25 Myr between snapshots, which is essential for resolving the halo response, since typical dynamical times in the halo are ∼200-500 Myr.The spatial resolution is sufficiently high to account for dynamical torques and the production of wakes.The feedback model implemented in FIRE leads the global structure of the central galaxy at z = 0 to be consistent with observational distributions (Hopkins et al. 2018;Sanderson et al. 2020;Bellardini et al. 2021Bellardini et al. , 2022)).We thus expect our results to be relevant for present and future studies of the shape of the MW halo.
We analyze in detail the DM halos of four Latte galaxies (m12f, m12i, m12m, and m12w) described in Table 1.These four simulations are chosen since they possess LMC-like satellites.The criterion for "LMC-like" is defined in the following section.Additionally, two of the LMC analogs (m12f and m12w) used in this analysis were first presented in Samuel et al. (2021).The visualizations14 of these galaxies in the Appendix illustrate the visible differences between the central galaxies due to their diverse formation histories.
Given that LMC analogs can powerfully influence halo dynamics (Garavito-Camargo et al. 2021), we categorize and trace an "LMC analog" in each simulated halo by selecting the most massive and most recent/ongoing merger with pericenter distance similar to the LMC.We define LMC-like satellites as simulated dwarf galaxies that meet the following criteria: (a) it reaches at least a total mass of 4 × 10 10 M e or a stellar mass of 5 × 10 8 M e ; (b) the satellite crosses into the virial radius of the host; and (c) is not destroyed after 6 Gyr after the simulation begins.We refer to the interaction as a "major merger" if the LMC analog has an infall mass-to-central galaxy bound mass ratio of M LMC /M MW > 1/3 and we refer to it as a "minor merger" otherwise.
All four galaxies have either a major or minor merger with a simulated LMC analog across their evolutionary histories.These simulated mergers allow us to test the effect of an LMC-like merger on the DM symmetry axes as a function of redshift and mass ratio: 1. m12i is a relatively stable galaxy with no recent major mergers, which makes it an excellent baseline condition at z = 0 for comparison to simulations with recent major mergers.It has an LMC analog that has a minor merger with the host at z = 0.6.2. m12f has a recent major merger (  M M LMC 1 3 MW ) involving an LMC-like satellite with stellar mass of ∼10 9 M e .3. m12w has a major merger with an LMC companion, but has few stellar streams at present. 4. m12m is a galaxy with a minor LMC merger 3 MW ) at relatively early times.It also has the most massive stellar disk, which forms very early in the simulation.
The properties of the LMC analogs selected for this work are detailed in Table 2.We emphasize that the distinction in LMC mass ratios is critical, because the mass ratio at pericenter is a key indicator with which a potential can be accurately modeled in terms of axisymmetry, and it may impact the galaxy-halo alignment (Arora et al. 2022).In reality, Arora et al. (2022) find that the MW-LMC mass ratio is high enough such that loworder axisymmetric multipole modeling is insufficient in describing the galactic potential.Transitively, we pay attention to the mass ratios of the simulated MWs and LMCs to probe mass-dependent dynamic effects on the DM symmetry axes.

Determining DM Symmetry Axes
Figure 1 illustrates the process of determining the DM symmetry axes.We find the symmetry axes of the DM halo by diagonalizing the reduced moment-of-inertia tensor: Notes.All simulations have identical particle mass for DM (3.5 × 10 4 M e ) and initial particle mass for baryonic matter (7.1 × 10 3 M e ), as well as identical force softening parameters ( = 1.0 pc gas, min  , ò star = 4.0 pc, and ò dm = 40 pc).M 200 m and R 200 m : the virial mass and radius of the main halo at z ∼ 0. M *,90 and R *,90 : 90% of the stellar mass contained within 30 kpc of the galaxy and the spherical radius that encloses that mass.z 15 (z 2 ): redshift when the cumulative fraction of stars that formed in situ exceeded 0.5 when selecting stars at z = 0 within host-centric distances of 15 (2) kpc (Santistevan et al. 2020).z disk and t disk : redshift and cosmic time when the disk is roughly assembled.d 1 : the local scattering region radius, which is defined as the radius at which DM particles should have experienced at least a single self-interaction (Vargya et al. 2022).Assembly times for SIDM halos (i.e., where values replaced by †) are unable to be calculated, since star particle formation distances are not available for those simulations.where V is the set of position vectors u of each DM particle and where q = b/a (the ratio of the intermediate-axis to the majoraxis lengths) and s = c/a (the ratio of the minor axis to the major axis).We initialize q = s = 1 for the sake of simplicity, assuming that the axis orientation is independent of the axis ratios (Allgood et al. 2006;Vargya et al. 2022).
We calculate and diagonalize the reduced inertia tensor for each galaxy at different selection radii (i.e., selecting DM particles from a sphere of greater volume with each iteration).The eigenvectors from the resulting diagonalization are the symmetry axes (u, v, and w, for major, median, and minor axes) for the DM halo at that particular selection radius.The minor axis of the DM halo is defined to be the eigenvector with the smallest corresponding eigenvalue.The axis lengths of each symmetry axis (eigenvector) are defined to be the square root of each eigenvalue.
We emphasize that the particles are selected within ellipsoidal volumes rather than shells.An alternate method developed by Zemp et al. (2011) calculates local matter distributions by using ellipsoidal shells without distance-based weighting and omitting subhalo particles.We opt for volume selection, since the axes fitted to the particles within that volume can inform the force that stellar tracers are sensitive to.
Additionally, Emami et al. (2021a) utilize a volume-based method to determine their symmetry axes for DM halos.When comparing the eigenvalues of the volume-based method to a replicated analysis using a local shell-based method, they find that they are roughly similar, with the caveat that the shellbased method had more noise due to its sensitivity to local fluctuations in density (Emami et al. 2021b).Ultimately, they opt to use the volume-based approach, as it gives a better average to the shape and rotation of the halo.We adopt this approach in our work, as the simulations are at a higher resolution compared to TNG50 and the effect of noise due to local rotations induced by local density variations may be enhanced.
A mathematical motivation for the volume-based approach is that while Poissonʼs equation relates the potential and density of the halo, the intermediate derivative corresponds to the volume integral of the density.This quantity, in principle, is responsible for generating the gravitational force and is thus the quantity that dynamical galaxy models (and other processes that are sensitive to the mass density, such as lensing) are attempting to reproduce.
In this work, we opt to keep the subhalos within our selection sample.Simplified potential models of the MW do not attempt to model the influence of individual subhalos.To this end, this work is mainly geared toward detecting mismatches between this simplified model and a more realistic model due to the presence of massive subhalos.

Calculating the Obliquity Angles with Host Axes with Degeneracy Corrections
We define the obliquity angle of each simulated MW as the angle between the DM minor axis, determined as described in Section 2.1, and the short symmetry axis of the disk of the central galaxy, which is calculated in a similar manner to the DM, using the youngest 25% of in situ stars within a radius that contains 90% of the mass selected within 10 kpc of the halo center.In the galaxies we analyzed, this age selection corresponds to stars younger than 2-4 Gyr, depending on the formation history of each galaxy.Since the majority of the angular momentum of the galaxy is embedded within the gas   disk, tracing these stars will trace the gas disk closely, as they are expected to mostly originate from star formation within that gas, while we widen our selection so that we select enough stars to avoid shot noise and smooth out deviations from the youngest stars forming in spiral arms.We avoid adopting a circularity criterion to select the disk stars because our simulated galaxies possess complex structures such as bars or spiral arms that we wish to avoid by selecting young-but not solely the youngest-stars.
The obliquity angle is then determined from the dot product between the minor axis and the DM minor axis at a particular galactic radius and/or snapshot in time: The process of diagonalizing the moment-of-inertia tensor only determines the symmetry axes up to the sign of the eigenvectors, leading to a degeneracy that we need to account for.For degeneracy corrections, we choose the stellar disk angular momentum within r = 10 kpc at the present day as a reference vector.
To break these degeneracies, we invert the major (axis vector corresponding to length a) and minor (axis vector corresponding to length c) symmetry axes if the dot product between the DM minor axis and the disk angular momentum vector (the median angular momentum of the stellar disk within r = 10 kpc) is negative.
The angular momentum vector is defined as and the normalized angular momentum vector is For consistency, the minor axis is always defined as the cross product of the calculated major and median axis eigenvectors, such that our rotation tensor obeys the right-hand rule.
In Figure 2, the directions of the DM halo minor axes at different selection radii (blue arrow), the vector perpendicular to the stellar disk (red arrow), and the angular momentum vector (red) described in Section 2.2 are overlaid on the DM density map of m12i.The orientations of the DM halo at z = 0 are calculated based on the angle between the DM symmetry axes and the host symmetry axes, as shown in the figure.
For analyzing the orientation as a function of redshift, we select DM particles within initial spheres of radii of r = R *,90 , 5R *,90 , R 200 m , where R *,90 is the radius within which 90% of the star particles are contained and R 200m is the virial radius for each particular galaxy at a given redshift.We calculate the orientation relative to the galaxy axes at that particular redshift, as described in 2.1, with the short axis defining the orientation of the positive z-axis.
Additionally, we make note of the triaxiality as a function of redshift for each initial volume (see Section 2.4), as this can determine the axis of rotation of the symmetry axes.We caution that the shape measurements may be biased due to the presence of substructure, as mentioned in Section 2.1, which may dominate the shape measurement at the radius of the subhalos.

SIDM and DMO Simulations
Latte galaxies m12f, m12i, and m12m have additional simulations that are initialized with SIDM with baryons and CDM without baryons (Wetzel et al. 2016;Garrison-Kimmel et al. 2017, 2019;Sameie et al. 2021).We replicate the orientation calculation for the SIDM variants identical to the aforementioned process.However, this process needs to be tweaked for DMO simulations, given there is no stellar (baryonic) disk to serve as a reference frame.We define a pseudo-stellar disk by calculating the moment-of-inertia tensor at r = 10 kpc for each halo and using the short axis of the reduced tensor as the "stellar disk" reference vector.We caution that the orientation of the pseudo-stellar disk, i.e., the inner DM halo itself, would reflect stronger alignments to the outer halo, since the inner halo is no longer coupled to a distinct stellar disk.

Triaxiality and Determining Effect of LMC-like Mergers
We calculate the triaxiality (shape parameter) of the halo at the snapshot of interest.The triaxiality is defined as a function of the axis lengths (a, b, and c) that correspond to the major, median, and minor axes: We consider a halo to be prolate if 1 T > 2/3, triaxial if 2/3 > T > 1/3, and oblate if T < 1/3 (Allgood et al. 2006;Vera-Ciro et al. 2011;Prada et al. 2019;Vargya et al. 2022).
To constrain the triaxiality accurately, we choose (a, b, and c) using an iterative code that deforms the ellipsoid for a maximum of 100 iterations, as implemented in Vargya et al. (2022), at each snapshot in time (at different lengths).The iterative code optimizes the shape of the ellipsoid by recursively reshaping the selection geometry to the previously estimated ellipsoid shape.This process is repeated until either the maximum iteration count has passed or the ratios of the axis lengths converge where (a, b, and c) are the axis lengths.
This method is identical to the method described in Section 2.2, except for the triaxialities being determined by iterating over 100 selection ellipsoids, rather than the single iteration of the code used to determine the symmetry axes.We find that iterating the selection ellipsoid leads to a betterconstrained halo shape at the cost of poorly constrained obliquities in toy models.
We hypothesize that the minor axis of the main halo should be roughly orthogonal to the line joining the centers of the main galaxy and the LMC analog, while the major or intermediate axis of the main galaxy should point toward the LMC analog, since we expect the interaction to shift the DM particle density in this plane.Therefore, we also calculate the alignment between the host galaxy major axis and the direction to the LMC analog.Specifically, we calculate the alignment of the DM halo to the simulated LMC satellite, θ LMC , as the minimum angle between the major-axis vector a of the DM in the host galaxy and the position vector of the satellite x LMC , defined relative to the host galaxy center.We calculate θ LMC at three different scale lengths: at the stellar disk size R *,90 , at 5R *,90 , and at the virial radius of the halo R 200m :

Settling of the Galaxy Orientation and Formation of the Galactic Thin Disk
Given the moment of inertia of the stellar component of the galaxy as a reference frame for calculating the halo position angle, it is of interest to know when the thin stellar disk begins to form and if it forms before or after the stellar component vector becomes influenced by the orbital vector of the LMC analog.The coherence of the thin disk appears when the starforming H I gas becomes rotationally supported within the galaxy.We quantify the rotational support strength of the cold H I gas as 〈v f 〉/σ v,f , which is the ratio between the average azimuthal velocity and the azimuthal velocity dispersion of the H I gas.This ratio describes the rotational support of gas particles within the inner galaxy (r < 0.05R vir ) as the mean rotational motion of the cold gas against the random azimuthal motions.We define the time at which the thin disk has or begins to appear coherent as t disk , which is determined by the transition time when roughly the rotational support ratio attains half the maximum mean rotational velocity of the cold gas.

Halo Orientation and Shape at z = 0
In the upper panel of Figure 3, the DM halo symmetry axes have greater obliquities at greater radii.We note that m12f is more highly aligned to the stellar disk, as its orientation at r = 100 kpc is θ < 10°compared to the θ ∼ 20°orientation of the other galaxies, possibly due to the merger happening edgeon to the stellar disk.Most notably, we find that m12w diverges in orientation significantly at r 100 kpc.This may be attributable to the aforementioned massive satellites that are infalling within 100-200 kpc of the halo center.Overall, these findings are consistent with the tentative evidence of a radially varying halo orientation in the MW, as found in Vasiliev et al. (2021).
To quantify at what particular radius the halo diverges, we calculate the fraction of the virial radius at which a particular halo tilts more than 20°relative to the disk.We present the divergence radii of these galaxies in Table 3. Halos that meet this divergence criterion diverge at about 0.2-0.4R 200m (Table 3).In comparison to a study of halo orientations (relative to the angular momentum of the stellar disk), Prada et al. (2019) find that a plurality of the simulated Aquarius halos aligns well to the stellar disk at all radii (most not exceeding 20°-30°in obliquity).Prada et al. (2019) also find that some halos experience extreme twisting (some being near perpendicular to the disk even at 0.5R 200m ).Although our analysis makes use of a much smaller halo sample, we observe similar trends: m12f is a prime example of a well-aligned halo at all radii and m12w is a prime example of extreme twisting.This trend is also observed in the TNG50 MW-like galaxies from Emami et al. (2021a), where the angle between the angular momentum vector of the stellar disk and minor axis of the DM halo grows and exceeds 20°at between 10 and 100 kpc (∼70 kpc, on average).Three of the FIRE galaxies in Table 3 diverge at similar radii.Compared to the measured disk-halo tilt in Han et al. (2022), none of the four FIRE galaxies exceed ∼30°tilt at 30 kpc.
The lower panel of Figure 3 shows the triaxiality of the DM halo as a function of radius at z = 0 for the four galaxies.m12f and m12m are either oblate or triaxial within 100 kpc, while m12i and m12w extend into prolate configurations farther from the galactic center.In particular, between 100 and 150 kpc, m12f experiences a dip in triaxiality.Further investigation reveals the presence of an infalling satellite within that range with a bound mass of 1.4 × 10 11 M e .Similarly, m12w has a triaxiality shift between 100 and 200 kpc, due to several satellites being within the vicinity (the most massive having a bound mass of 8.1 × 10 10 M e ).The shapes of the halos are possibly biased by the presence of these subhalos, and they may not reflect the true shape of the local mass distribution.Furthermore, the inner halos of all these galaxies are oblate, which is consistent with the distribution of the inner DM halo being influenced by the baryonic disk.However, recent phasespace analysis of the Helmi stream suggests that the inner halo may have a prolate shape, although action calculations may be invalid if the stream substructures are in resonance (Dodd et al. 2022;Craig et al. 2023).They acknowledge that the potential model used in their analysis ignores possible dynamical influences from massive satellites, e.g., introducing reflex motion of the MW toward Sagittarius (Vasiliev et al. 2021).

Halo Orientations across Redshift
In Figure 4, we show the orientations of the DM halo over time with respect to the evolving galaxy axes.We observe central radii aligning well to the galaxy symmetry axes as the galaxies approach z = 0. Furthermore, at the virial radius, we see that the minor axes are less aligned to the galaxy as the simulations approach z = 0.
In the outer halo (r = 5R *,90 , R 200 m ), we observe the halos with a major merger quickly aligning toward the disk as the LMC analog reaches pericenter with the halo.While the trend is present in the major-merging halos, the minor-merging halos have relatively smooth transitions.
The outer halo does not shift its density to be consistent with the inner halo; however, the inner halo orients the disk toward the LMC analog as it approaches pericenter.In Figure 5, we show the time-evolving orientations of the halo relative to the present-day axes.At around pericenter, we observe that the inner halo aligns to the present-day disk within 1-2 Gyr of pericenter of a major-merging satellite.This follows from the quick response of the inner halo to dynamical perturbations.In our CDM sample, the average dynamical timescale within the present-day divergence radius is roughly 2-4 Gyr, which is fast enough to allow for the inner halo/disk to respond to the LMC analog.Although this response is noticeable in the minormerging m12i, it is comparably weaker than the more massive satellite mergers.
These results suggest that massive LMCs (∼10 11 M e ) play a critical role in determining the orientation of the MW's disk at  ) and disk size (R *,90 ).The only SIDM/DMO halo that diverged was m12m (SIDM).
the present day.This result is consistent with the findings in Santistevan et al. (2021), where the orientation of the disk is heavily influenced by the most recent massive gas-rich merger.

SIDM and DMO Comparisons
We extend our analysis to three available Latte galaxies (m12i, m12f, and m12m) initialized with SIDM with baryons and CDM-only simulations.Our results, shown in Figure 6, plot the present-day orientations for these different simulations.Compared to the original orientations at the present day (see Figure 3), the SIDM and CDM-only halos in m12i and m12m have smaller halo-galaxy obliquities as a function of radius.Furthermore, we find that the only halo that diverges (>20°) is the SIDM halo in m12m at 0.7 R 200 m .In m12m, we see that  the orientations between CDM, DMO, and SIDM are similar to each other.In Prada et al. (2019), they find that the DMO orientations are weakly correlated with each other at all radii.This finding appears to be consistent when comparing the DMO and CDM orientations.More insights are to be gained from future analyses of the correlation between the SIDM and DMO orientations.
In our analysis of SIDM galaxies, we make note of the orientations calculated within d 1 (the local self-interaction scattering region; see Table 1 for typical ranges of this value).This particular halo region is of importance, given that the physical differences between SIDM and CDM halos become significant due to SIDM particles having experienced at least one self-interaction.However, we observe that the differences in orientation between CDM and SIDM halos within d 1 are not apparent.This result indicates that the effects of self-interaction do not impact the orientation of the halo relative to the disk at small scales.
Given that the elastic SIDM orientations appear to be oriented closer to the disk at all radii, we expect this effect to be even stronger in alternative SIDM scenarios, such as the dSIDM model (Shen et al. 2021).Within the dSIDM model, DM can form dark disks or inner DM halos that behave as oblate rotators, which can more strongly force alignment with the baryonic disk (Shen et al. 2021(Shen et al. , 2022)).

Orientation of DM Halo Relative to LMC Infall
In Figure 7, we show the orientations of the DM halo relative to the LMC position as a function of time at different radii.We  analyze all four CDM+baryon galaxies with different LMC analog mergers.We indicate the time of the LMC pericenter, the total peak mass, and the different present-day scale lengths in Table 2.The major-merging ( > M M LMC 1 3 MW ) galaxies are m12f and m12w, while the minor-merging ( < M M LMC 1 3 MW ) galaxies are m12i and m12m.Additionally, in Figure 8, we present these orientations as a function of distance to the central galaxy.
We observe that the major merger in m12f corresponds to increased major-axis alignment to the LMC position over time, with a sudden period of misalignment when the LMC reaches its pericenter (at r = 36 kpc).In m12w, we see the major axis of the galaxy aligning well over time as the satellite approaches pericenter, with a short burst of misalignment following pericenter.The misalignment may be attributed to the shifting center of mass of the central galaxy as the LMC passes pericenter.Across both major-merging galaxies, their outer halos are strongly aligned to their respective LMC analogs (Figure 8).
In our minor-merging galaxies, m12i and m12m show a weak correlation between the halo orientation and direction of the merger as a function of time and distance (Figure 8).We hypothesize that due to the lower-mass LMC, the satellite applies less force on the host axes and does not influence the axes significantly.
Overall, the orientations of the outer halo in Figure 7 appear consistent with the major axis (a density dipole) aligning and tracking the LMC over time, but the mass must be  M 1 3 MW for this to happen.This is consistent with Vasiliev et al. (2021), where they find that the outer halo elongates toward the previous position of the LMC a few hundred million years ago in fitted static MW potentials.
Across both minor-and major-merging galaxy samples, we observe that the outer halo tends to be more well aligned to the infalling satellite position than the inner halo.These results are consistent with how simulated LMCs can effectively decouple the inner and outer DM halo, as per the findings of Garavito-Camargo et al. (2019).These results are similar to the findings in Shao et al. (2021), where they explored the twisting of the DM halo in MW-like galaxies in EAGLE simulations.They find that the outer halo (at R 200 m ) is more aligned to the vast orbital structure of satellites and that the presence of an LMCmass satellite does not affect the orientation of the outer halo itself, but can decouple the orientation of the central disk relative to the outer halo (Shao et al. 2021).However, they find that they cannot constrain a typical scale length ("twist radius") at which the inner and outer halo diverge, finding that this transition radius can vary between 30 and 150 kpc (Shao et al. 2021).
In Figure 9, we present the triaxiality as a function of time.The triaxialities for major-merging galaxies appear to be more erratic over time.On the other hand, minor-merging galaxies have a more gradual triaxiality transition as a function of z.The DM halo of m12i appears to have a stronger change in shape as the LMC reaches pericenter compared to the other minor merger in m12m.This is attributed to the extremely close pericenter (∼30 kpc) that the LMC analog makes with the host galaxy of m12i, while the merger in m12m has an orbital pericenter of ∼150 kpc.Additionally, the shape of m12i has been found to evolve distinctly from the other MW-mass galaxies in an analysis by Talei et al. (2022).They attribute this distinct shape evolution to its unique late-forming assembly history.

Orientation of the Halo Short Axis within Filaments
At the present day, all the galaxies in our analysis have their minor axes in the outer halo aligned roughly perpendicular to the enveloping filament (see Figure 2 for the case of m12i).(2011), where they identify a trend of halos with short axes (at the virial radius) aligning perpendicular to the filament direction.
This pattern of halo-filament orientation has a critical implication: infalling substructure mass will preferentially accrete along the major axis of the halo (Vera-Ciro et al. 2011).In the context of m12w, we observe significant galaxyhalo decoupling in the outer halo (likely due to the major merger), which indicates that the galaxy may become independently oriented to the filament.In Figure 13, we observe that m12w is located at the nexus of three filaments that are accreting massive subhalos near the central galaxy.All the other Latte galaxies appear to be roughly perpendicular to the main accreting filaments they are embedded within, but m12w experiences subhalo accretion from multiple directions, causing the DM halo axes to twist significantly as numerous infalling substructures enter the selection aperture between the twist radius (R ∼ 100 kpc) and the virial radius.Future analysis of the environment of the filamentary substructure and its relation to galaxy/outer-halo orientation in Latte galaxies may be critical to understanding the mechanisms behind this observed decoupling.

Formation of the Thin Disk and the Stability of the Galaxy Direction
Figure 10 presents the rotational support ratio 〈v f 〉/σ v,f of cold H I gas in each galaxy and the approximate time at which the thin disk appears (t disk ).We find that three of four galaxies in our CDM+baryons sample form their thin disks after the first pericentric passage of the LMC analog.The cold-gas component that eventually forms in the thin disk is deposited on the same orbital vector as the massive merging satellite (Santistevan et al. 2021), and we expect that the thin disk imprints the direction of the merger and, by loose correlation, the filament direction at the time of the merger.
From our results in Figures 4 and 5, it appears that the direction of the stellar component is stabilized due to the merger, and when the thin disk forms it is oriented along the direction of the mergerʼs angular momentum, while simultaneously the inner DM halo strongly aligns to this stellar axis configuration.
We believe that deviations or twists from the inner halo (which should be set by the stellar disk) to the outer halo are attributed to the evolution of the filament direction.For example, m12w is presently accreting massive substructure along three different filamentary directions-most likely causing the DM symmetry axes to distort as those substructures accrete into the moment-of-inertia aperture.

Conclusions
In this work, we have investigated the orientations of the DM halo symmetry axes as a function of radius and time for DM, SIDM, and DMO halos.Furthermore, we have explored how midscale to major mergers (1/10-1/3 M MW ) can influence these symmetry axes in terms of orientation relative to the stellar disk and halo alignment to the LMC.The main conclusions we can draw from this analysis are: 1.At z = 0, the orientations of the DM halo are well aligned to the stellar disk at central to intermediate radii (0.2-0.4 R 200m ) and diverge at greater radii (>0.4 R 200m ; Figure 3).These results are consistent with the fact that there is no a priori reason that the DM halo should be aligned with the disk axes at all radii.Galactic potential models that assume DM-stellar disk alignment may be insufficient in describing the true global potential in the outer halo.together with instruments such as Rubin can give us the necessary observations required to detect streams far out enough in the MW to detect divergences in orientation.3.In the presence of an infalling LMC analog, we find differential major axis-LMC alignment behavior depending on the mass of the analog satellite.In Figure 11, galaxies undergoing a major LMC merger (m12f, m12w) have their DM major axes align to the satellite as infall progresses.However, galaxies undergoing a minor LMC merger do not have an appreciable effect on the direction of the major axes.4. When viewing the alignment as a function of radius, we see that major-merging galaxies show a distinct radius alignment relationship (Figure 8).This apparent trend is not replicated in minor-merging galaxies.We note that this relationship may be modulated by the angle of LMC infall.We speculate that the lower-mass LMC satellites do not generate a significant amount of torque on the halo to strongly shift the DM density, leaving the axes roughly stationary (see Figure 13).5.The outer halo more closely aligns to the LMC analog position over time than the inner halo, showing the presence of the differential response of the halo relative to the LMC across galactic radii.
6.In simulations initialized with SIDM (DMO), the orientations to the stellar disk (pseudo-stellar disk) are more aligned compared to the orientations of CDM halos with baryons.We predict this effect to be stronger in dSIDM models.7. We find evidence that the orientation of the thin disk is set by the orbital direction of the merger, and that the inner DM halo simultaneously couples to the stellar axes with the same configuration.Deviations to these axes in the outer halo are believed to be caused by changes in the direction of substructure accretion, due to the evolution of filamentary structure.
An observational test of the halo orientation can be useful in explaining phenomena such as the figure rotation of the DM halo (Valluri et al. 2021).A conceivable observational method for determining the halo orientation may come from fitting the MW DM distribution function using data from RR Lyrae stars (Hattori et al. 2021).This method has been used to fit the halo shape within the disk using observational data from Gaia, but it assumes a perfectly aligned halo (Hattori et al. 2021).Extending this method to tracer stars beyond the baryonic potential of the disk could be useful in constraining the true MW orientation.
Next-generation observations of the MW will reach the orientation transition radius by enabling the detection of these faint tracers in the outer halo.In Figure 3, we indicate the different distance ranges for different instruments as well.The Gaia satellite, when observing giant branch stars, is limited to observations at 60 kpc (Sanderson et al. 2019).However, using ground-based asteroseismology measurements of M-giants, as detailed in Auge et al. (2020), we may be able to measure stellar (and inherently stream) features up to 200 kpc.We see that in m12w (a case of extreme obliquity), we may be able to resolve this orientation signal through observation.Additionally, the next-generation LSST and its instruments can provide us with up to 400 kpc (beyond the virial radius) of distance and will be more than capable of resolving the stream features necessary to detect divergences in the halo orientation (Sanderson et al. 2019).
Critical conclusions from future observations and simulations can help resolve theories on Modified Newtonian Dynamics (MOND), where an oblique halo relative to the disk challenges a requirement of MOND that the net potential of the halo and the disk should be well aligned (Debattista et al. 2013).Thus, continuing our investigations of simulated DM halos can provide us with powerful predictions of outer-halo dynamics in advance of future observational studies.

Figure 1 .
Figure 1.The selection radius defines the enclosed spherical region where DM particles are chosen.The moment-of-inertia tensor is determined from the enclosed DM particles.Symmetry axes of the DM halo (labeled a, b,and c) are defined by the eigenvectors of the reduced moment-of-inertia tensor.The length of each symmetry axis is the square root of the corresponding eigenvalues.These eigenvectors and eigenvalues define the converged halo shape.

Figure 2 .
Figure 2. DM density plot of m12i with host minor axes overplotted at z = 0.The red arrow indicates the direction of the stellar disk/galaxy vector.The blue arrows indicate the minor axes of the halo, with the lengths scaled to the radii at which the principal axes were calculated.Each column shows a different projection of the halo.Brighter regions indicate higher densities of DM. Figure 13 replicates this density map for all galaxies.

Figure 3 .
Figure 3. Upper panel: DM host minor-axis orientation plots of m12f, m12i, m12m, and m12w as a function of radius at z = 0.The distance ranges for Gaia (for MSTO), ground-based M-giant asteroseismology, and LSST are shown with the different line styles.Lower panel: DM triaxiality plots of m12f, m12i, m12m, and m12w as a function of radius at a redshift of z = 0. Triaxiality is determined by the convergent ellipsoid at each radius.The horizontal lines demarcate the transitions between oblate (T < 1/3), triaxial (1/3 < T < 2/3), and prolate (T > 2/3) shape classifications.This plot is available as a function of time in Figure 9.

Figure 4 .
Figure 4.The misalignment angle of the DM halo minor axis relative to the stellar component over time at different scale lengths.The yellow solid line is the orientation of the DM halo at the edge of the stellar component of the galaxy (R *,90 ), the green dashed line is the orientation outside the galaxy but well within the DM halo (5 × R *,90 ), and the blue dashed-dotted line is the global orientation of the DM halo at R 200m .Orientations at each snapshot are calculated with respect to the galaxy vector at that particular snapshot.The vertical black line indicates when the LMC analog reaches pericenter for the particular galaxy and the red line indicates when the thin disk of the galaxy begins to form.

Figure 5 .
Figure 5.The same as Figure 4, but the orientations at each snapshot are calculated with respect to the stellar disk at the present day.

Figure 6 .
Figure6.Halo orientations relative to the stellar disk for different simulations at z = 0.Each galaxy begins with identical initial conditions, with the exception of how the DM interacts within the system.The gray line represents the orientation of the halo with CDM and baryons (see Figure3), and the light blue (dark blue) lines correspond to DMO (SIDM) halo orientations.Additionally, we plot distance ranges for Gaia (MSTO), M-giant asteroseismology, and LSST.The DMO pseudostellar disk axis is defined to be the angular momentum of the DM particles within 10 kpc.

Figure 7 .
Figure 7. Angle between the major axis and the LMC analog position.The black line indicates the time when the satellite reaches pericenter and the red line indicates when t bursty occurs.

Figure 8 .
Figure 8. Major-axis alignment to LMC as a function of the distance to the host.Different selection scale lengths for the DM halo are coded by color.The virial radius tends to be more closely aligned with the LMC than the inner halo, showing differential alignment to the LMC over time.

Figure 9 .
Figure 9. Halo triaxiality as a function of time at different scale lengths.Each galaxy's LMC pericenter is indicated with a vertical black line, and the time at which the disk is formed is shown with the vertical red line.

Figure 10 .
Figure10.〈v f 〉/σ v,f ratios for CDM+baryon galaxies.The vertical red line indicates the time when the thin disk becomes roughly assembled.The vertical black line indicates the time of the LMC analog pericenter.The reference frame is defined as the cylindrical principal axes, where the z-direction is defined by the minor axis of the moment of inertia of young stellar particles.

Figure 11 .
Figure 11.Position of the LMC satellite (circle) relative to the galaxy (origin; in the principal/disk axis frame) and the direction of the galaxy's major axis (crosshair, exaggerated scale).The color indicates the cosmic time at which the major axes and LMC positions were measured.

Figure 12 .
Figure 12.The same as Figure 11, except with the axes fixed in the latest tracked reference frame of the disk.

Figure 13 .
Figure 13.The DM densities within each Latte CDM galaxy.The red arrow indicates the direction of the stellar minor axis.The blue arrows indicate the directions of the DM halo minor axis, with the lengths scaled to the radii at which the axes are determined.The dashed green inner circle indicates a radius of 100 kpc.

Table 1
Properties of Simulations Used in This Work

Table 2
Properties of Identified LMC Satellites Used in This Work

Table 3
Approximate Radii at Which Each Latte Galaxy's Halo Orientation Exceeds 20°(Divergence Radius r div )