# A Random Walk Model for Halo Triaxiality

Paul Menker,1,2 and Andrew Benson,2
1University of Southern California, 3551 Trousdale Pkwy, Los Angeles, CA
90089, USA
2Observatories, Carnegie Institution for Science, 813 Santa Barbara Street,
Pasadena, CA 91101, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We describe a semi-analytic model to predict the triaxial shapes of dark
matter halos utilizing the sequences of random merging events captured in
merger trees to follow the evolution of each halo’s energy tensor. When
coupled with a simple model for relaxation toward a spherical shape, we find
that this model predicts distributions of halo axis length ratios which
approximately agree with those measured from cosmological N-body simulations
once constrained to match the median axis ratio at a single halo mass. We
demonstrate the predictive and explanatory power of this model by considering
conditioned distributions of axis length ratios, and the mass-dependence of
halo shapes, finding these to be in good agreement with N-body results. This
model provides both insight into the physics driving the evolution of halo
triaxial shapes, and rapid quantitative predictions for the statistics of
triaxiality connected directly to the formation history of the halo.

###### keywords:

cosmology: theory, dark matter

††pubyear: 2015††pagerange: A Random Walk Model for Halo Triaxiality–A

##  1 Introduction

Dark matter halos, which form through the gravitational collapse of overdense
regions of the universe, have been found to be significantly non-spherical via
analysis of N-body simulations (Jing & Suto, 2002; Hopkins et al., 2005;
Allgood et al., 2006; Hayashi et al., 2007; Macciò et al., 2008; Vera-Ciro et
al., 2011; Despali et al., 2014; Bonamigo et al., 2015; Despali et al., 2017;
Vega-Ferrero et al., 2017). Such halos have been found to be well-
characterized as triaxial ellipsoids. Halos measured in N-body simulations are
found to be to mostly prolate, although with a broad distribution of shapes,
with lower mass halos being closer to spherical than higher mass ones (e.g.
Jing & Suto, 2002; Bonamigo et al., 2015). The triaxial shapes of halos have
also been found to exhibit weak trends with environment and halo formation
history, with higher mass halos in low density environments (e.g. voids and
walls of the cosmic web) being more prolate (e.g. Macciò et al., 2007; Lee et
al., 2017; Gouin et al., 2021; Lau et al., 2021; Hellwing et al., 2021). These
expectations from theoretical studies are broadly consistent with
observational studies (e.g. Morandi et al., 2011; Limousin et al., 2013; van
Uitert et al., 2017).

This triaxial shape has important consequences for the orbital evolution of
subhalos (and their associated satellite galaxies) leading to anisotropy in
the spatial distribution of subhalos around their host (Zentner et al., 2005;
Kuhlen et al., 2007), the strength of gravitational lensing exhibited by a
halo (Corless & King, 2007), cluster mass and concentration estimates
(Limousin et al., 2013), and two-point correlation functions (van Daalen et
al., 2012). Understanding the origins and statistical distribution of halo
shapes is therefore important to many astrophysical measurements.

In addition to these characterizations of halo triaxiality using measurements
from N-body simulations, there has been some analytic work aimed at
characterizing halo shapes, typically based upon the statistics of Gaussian
random density fields (Eisenstein & Loeb, 1995; Rossi et al., 2011).

Previous work, utilizing N-body methods, exploring the effects of mergers on
halo shapes (Moore et al., 2004; Drakos et al., 2019) has found the radial
mergers tend to produce prolate halos while mergers with higher angular
momentum tend to produce oblate halos, indicating that halo shape is directly
affected by the dynamics of halo mergers.

In this work, we construct a model for predicting dark matter halo shapes
using a different approach. We extend the semi-analytic framework of Benson et
al. (2020) and Johnson et al. (2021), in which key properties (spin and
concentration respectively) of halos are estimated by following the
accumulation of conserved quantities (angular momentum and energy
respectively) during each merger event in a halo’s formation history.
Specifically, in this work, we follow the accumulation of the energy tensor
during each merger to predict the triaxial shape of each halo. Unlike angular
momentum and the (scalar) energy, the energy tensor is not a conserved
quantity. However, we will show that it is approximately conserved during
mergers of halos to a sufficient level for this approach to succeed.

Benson et al. (2020) and Johnson et al. (2021) demonstrated that this merger
tree-based, semi-analytic framework can accurately predict the distribution of
halo spins and the distribution and mass-dependence of halo concentrations,
and provides a reasonable estimate of the auto-correlation over time of these
quantities as halos evolve. This approach has the advantage (compared to other
empirical and semi-analytic methods) that a halo’s properties are predicted
directly from its merger history. Since the statistical properties of such
merger histories depend fundamentally on the power spectrum of cosmological
density perturbations, and the cosmological model, this approach therefore
connects the statistics of these key physical properties of halos to the
underlying cosmological model, while also directly capturing correlations of
these properties with the specific formation history of a given halo.

In this work we will explore some of these same issues when this framework is
applied to the energy tensor to predict the triaxial shapes of halos. Our goal
is to both test the hypothesis that the distribution of halo shapes can be
directly related to halo merger histories, and to provide an approximate model
for halo shapes which connects, through the merger history, to fundamental
quantities such as the cosmological model and power spectrum. The remainder of
this paper is organized as follows. In §2 we describe how the existing
framework is applied to the case of the energy tensor and halo triaxiality. In
§3 we show the resulting distributions of halo triaxial shapes and compare to
measurements from N-body simulations. In §4 we discuss the successes and
limitations of this model, and in §5 we present our conclusions.

##  2 Methods

It is well-established (Navarro et al., 1997) that the density profiles of
spherically-averaged cold dark matter halos can be accurately described by the
Navarro-Frenk-White (NFW) profile given by:

$\rho(r)=\rho_{*}\left[\frac{r}{R_{\mathrm{s}}}\left(1+\frac{r}{R_{\mathrm{s}}}\right)^{2}\right]^{-1}.$
(1)

This distribution, which appears to hold at all scales (Wang et al., 2020),
describes a dense central cusp with $\rho\propto r^{-1}$, transitioning to a
more diffuse outer region with $\rho\propto r^{-3}$ at a characteristic scale
radius, $R_{\mathrm{s}}$. The parameter $\rho_{*}$ provides the overall
normalization of the density profile. The scale radius is often defined via
the relation $R_{\mathrm{s}}=R_{\mathrm{vir}}/c$ where $c$ is known as the
concentration parameter, and $R_{\mathrm{vir}}$ is the virial radius, defined
as the radius enclosing a certain mean density contrast111Throughout this work
we define $\Delta_{\mathrm{vir}}$ using the spherical collapse model (see, for
example, Eke et al., 1996, appendix A)., $\Delta_{\mathrm{vir}}$. Cold dark
matter halos are formed through a process of hierarchical merging (White &
Rees, 1978) of smaller halos, the remnants of which can persist for
cosmological timescales as ‘subhalos’ orbiting inside the host halo into which
they merged (Moore et al., 1999). A useful representation of this hierarchical
merging process is the “merger tree”—a graph that tracks the history of these
mergers, and which can be augmented by associating properties with each node
(halo) in the tree, such as halo mass, virial radius, and orbital properties.

The majority of prior (semi-)analytic work assumes a spherical halo (although
see Lee et al., 2005; Smith & Watts, 2005). It is well-established however
that cold dark matter halos formed in cosmological simulations are
significantly non-spherical, being better described as triaxial ellipsoids
(Frenk et al., 1988; Bett et al., 2007). In this paper, we utilize and extend
the merger tree-based approach of Benson et al. (2020) and Johnson et al.
(2021) to compute the triaxial shape of each halo in a merger tree by tracking
the evolution of the energy tensor of each halo. Specifically, we update the
energy tensor of each halo after each merging event in the tree, and find the
resulting triaxial axis ratios.

In the remainder of this section we present N-body simulations of merging
halos, and describe a simple model of how such mergers drive the evolution of
triaxiality. In §2.1 we present a suite of idealized merger simulations of
triaxial halos and examine the shapes of the resulting halos. Then in §2.2 we
define the energy tensor and how it is changed during halo mergers, and in
§2.3 describe how these results are applied to entire merger trees. In §2.4.1,
we describe the free parameters of our model— one which is new, and two which
are identical to, and constrained by, the parameters of the Johnson et al.
(2021) model. An important difference with respect to the prior works of
Benson et al. (2020) and Johnson et al. (2021) arises because the energy
tensor, unlike the angular momentum (Benson et al., 2020) and the scalar
energy (Johnson et al., 2021), is not a conserved quantity—we discuss this
issue in §2.4.2. In §2.5 we validate our model by comparing its predictions
with the results of our N-body simulations. Finally, in §2.6, we explain how
these equations are utilized to follow the energy tensor, and triaxial halo
shapes, through full merger trees.

###  2.1 N-body Simulations

To explore how mergers between halos drive the evolution of halo triaxiality
we perform a suite of N-body simulations of isolated, merging halo pairs using
the Gadget-2 code (Springel, 2005). Halos are modelled as triaxial
exponentially-truncated NFW profiles (Kazantzidis et al., 2006, with
$r_{\mathrm{decay}}=r_{\mathrm{vir}}$) with isotropic velocity dispersions. We
generalize the spherical NFW profile (eqn. 1) to a triaxial form, following,
for example, Jing & Suto (2002) we define a new coordinate, $m$, through the
relation

$m^{2}(\textbf{{x}})=\sum_{i}\frac{x_{i}^{2}}{a_{i}^{2}},$ (2)

where the $x_{\mathrm{i}}$ are the position in a Cartesian coordinate system
aligned with the principle axes of the triaxial halo, and $a_{\mathrm{i}}$ are
the axis ratios of the triaxial halo, defined such that $a_{1}>a_{2}>a_{3}$,
and we set $a_{1}=1$ for convenience with no loss of generality. The density
at any point x is then simply $\rho(m[\textbf{{x}}])$ with $\rho(r)$ being the
spherical NFW profile as defined in equation (1).

Initial conditions for these triaxial halos are generated using the Agama code
(Vasiliev, 2019). We modify Agama to support the specific truncated, triaxial
NFW density profile described above, and approximate the corresponding
potential using Agama’s Multipole potential class with maximum angular orders
$l_{\mathrm{max}}=8$, $m_{\mathrm{max}}=6$. We then use Agama’s Schwarzchild
modeling capabilities to construct equilibrium particle distributions in a
variety of triaxial potentials. For each halo, we generate a library of
200,000 orbits, each integrated for 50 orbital periods, and then solve for the
weight of each orbit constrained to reproduce the target density profile and
isotropic velocity dispersion. From this library we then sample 2,000,000
particle coordinates in the 6-d phase space, which are then sub-sampled as
needed to construct our N-body halo realizations. For the subset of spherical
halos that we consider we also generate initial conditions using Eddington’s
formula (as implemented in Galacticus; Binney & Tremaine 2008, eqn. 4.46b). We
find that the results of our simulations when using spherical halos are
identical (to within statistical noise) for Agama or Galacticus-generated
initial conditions).

Table 1: Parameters of idealized merger simulations. In this table, “Label” is used to identify each simulation in the text, “$\xi$” is the mass ratio of the merging halos, “$\theta$” is the angle between their initial relative velocity vector and the line connecting their centers, “$c$” is the concentration of the halo, and “$a_{2}$” and “$a_{3}$” are the initial axis ratios of the halo. These last three parameters, $(c,a_{2},a_{3})$, are listed separately for the primary and secondary halos. |  |  | Primary | Secondary
---|---|---|---|---
Label | $\bm{\xi}$ | $\bm{\theta}$ | $\mathbf{c}$ | $\mathbf{a_{2}}$ | $\mathbf{a_{3}}$ | $\mathbf{c}$ | $\mathbf{a_{2}}$ | $\mathbf{a_{3}}$
a | 0.1 | 0 | 10.0 | 1.0 | 1.0 | 10.0 | 1.0 | 1.0
b | 0.1 | 0 | 9.0 | 1.0 | 1.0 | 12.0 | 1.0 | 1.0
c | 0.1 | $\pi/8$ | 10.0 | 1.0 | 1.0 | 10.0 | 1.0 | 1.0
d | 0.1 | $\pi/4$ | 10.0 | 1.0 | 1.0 | 10.0 | 1.0 | 1.0
e | 0.3 | 0 | 10.0 | 1.0 | 1.0 | 10.0 | 1.0 | 1.0
f | 1.0 | 0 | 10.0 | 1.0 | 1.0 | 10.0 | 1.0 | 1.0
g | 0.1 | 0 | 10.0 | 0.7 | 0.7 | 10.0 | 0.5 | 0.5
h | 0.1 | 0 | 10.0 | 1.0 | 0.7 | 10.0 | 1.0 | 0.4
i | 0.1 | 0 | 10.0 | 0.8 | 0.7 | 10.0 | 0.5 | 0.4
j | 0.3 | $\pi/4$ | 10.0 | 0.8 | 0.7 | 10.0 | 0.5 | 0.4
k | 0.5 | $\pi/8$ | 10.0 | 0.5 | 0.5 | 10.0 | 0.5 | 0.4
l | 0.2 | 0 | 10.0 | 1.0 | 1.0 | 10.0 | 0.5 | 0.5

We consider numerous combinations of primary and secondary halos, as detailed
in Table 1, spanning a range of mass ratios, angles of approach,
concentrations, and triaxialities. The primary halo always has a mass of
$10^{12}\mathrm{M}_{\odot}$. We use a particle mass of
$10^{7}\mathrm{M}_{\odot}$ and a softening length of 500 pc. The two halos are
initially separated by 300 kpc along the $x$-axis (approximately equal to the
virial radius of the primary halo). Both primary and secondary halos are
assigned randomly-selected orientations. Each pair of merging halos is then
allowed to evolve for 20 Gyr. We measure the energy tensor of the system (i.e.
including all particles, from both halos) at each snapshot of the simulation,
then find the eigenvectors and eigenvalues of the tensor allowing us to obtain
the components of the tensor in the coordinate system corresponding to the
principal axes of the system. We additionally determine axis ratios of the
merger product by measuring the inertia tensor, including all particles within
the approximate virial radius ($\approx R_{\mathrm{v,1}}(1+\xi)^{1/3}$ where
$R_{\mathrm{v,1}}$ is the pre-merger virial radius of the primary halo, and
$\xi$ is the mass ratio of the merging halos) of the merger product, and then
finding the eigenvalues of that inertia tensor.

|
---|---
Figure 1: Time evolution of the diagonal components of the energy tensor (and
scalar energy shown as $E/3$ for comparison with the expectation for the
diagonal components of the energy tensor for a spherical halo) in two N-body
simulations of mergers between idealized dark matter halos is shown by the
solid lines (numerical uncertainties on these lines—not shown for clarity—are
$\Delta E\approx 0.01\times
10^{16}\mathrm{M}_{\odot}\mathrm{km}^{2}\mathrm{s}^{-2}$ as determined from
bootstrap resampling). Left and right panels show results for simulations d
and g from Table 1 respectively. The merging halos are initially separated
along the $x$-axis and in the tangential merger case the initial velocity
vector lies in the $x$–$y$ plane. Subscript indices 0,1,2 correspond to the
$x,y,z$ axes respectively. Horizontal dashed lines indicate the initial value
of each component of the tensor. Dotted lines show the expectations from our
model, to be described in §2.2.

Figure 1 shows the time evolution of the diagonal components of the energy
tensor in simulations d and g from Table 1 in the left and right panels
respectively. The scalar energy is also shown as $E/3$ for comparison with the
expectation for the diagonal components of the energy tensor for a spherical
halo. (Subscript indices 0,1,2 correspond to the $x,y,z$ axes respectively.
Numerical uncertainties on these lines—not shown for clarity—are $\Delta
E\approx 0.01\times 10^{16}\mathrm{M}_{\odot}\mathrm{km}^{2}\mathrm{s}^{-2}$
as determined from bootstrap resampling.) Horizontal dashed lines indicate the
initial value of each component of the tensor. Dotted lines show the
expectations from our model, to be described in §2.2. We note that the total
energy (shown by the grey lines) is conserved throughout the simulation. There
are three distinct phases of evolution that can be identified in these
results. Between 0 and 2 Gyrs, corresponding to the violent merging phase, the
components of the energy tensor change rapidly. During this transitory phase
the potential is changing rapidly, and the remnant, merged halo has not yet
settled into a new equilibrium.

In the second phase, from around 2 Gyrs to 6 Gyrs, the remnant halo has
settled into a new equilibrium, and the components of the energy tensor are
relatively constant222In the left panel there are some fluctuations in the
values of the energy tensor components—particularly noticeable around 2 Gyr
where the $E_{00}$ and $E_{11}$ components switch values. This is a
consequence of the eigenvalue decomposition that we perform. Due to the
symmetry of the system and particle noise the labeling of axes (determined by
the rank order of eigenvalues) can switch. This does not affect the physical
attributes of the system, or our interpretation of the results.. Importantly,
they are constant at values close to their initial, pre-merger values. This is
consistent with our assumption that the energy tensor is approximately
conserved through the merger.

Finally, from about 6 Gyr onward, the components of the energy tensor slowly
converge towards $E/3$, indicating that the remnant halo is becoming more
spherical. This process of “ _sphericalization_ ”, in which the total energy
remains fixed, but the components of the energy tensor equalize over
timescales comparable to the age of the Universe is a consequence of the non-
conserved nature of the energy tensor. Interestingly, sphericalization appears
to happen on a timescale (around 5 Gyr) a few times longer longer than the
dynamical time of the halo (which is of order 1 Gyr). It is unclear if a large
merger always resets the process of sphericalization, or if they can speed up
the process in certain, coincidental circumstances.

The long-term (relative to the dynamical) stability of triaxial systems is
well-studied. de Zeeuw & Franx (1991), discussing equilibrium models for
elliptical galaxies, note that “the timescale for differential precession
between neighboring orbits is usually of the order 5–10 dynamical times and
becomes comparable to a Hubble time beyond a few optical radii.” Similarly,
Vasiliev & Athanassoula (2012) perform simulations of the evolution of cuspy
triaxial systems, finding that models with a weak ($\rho(r)\propto r^{-1}$)
density cusp are remarkably stable, while steeper cusps ($\rho(r)\propto
r^{-2}$) lead to evolution away from triaxiality due to chaotic diffusion of
orbits. Similarly, Merritt & Quinlan (1998), studying the effects of a central
black hole on the shape evolution of triaxial galaxies, found that shape
evolution occurs only on timescales of order hundreds of dynamical times in
the limit of small black hole masses. These results are consistent with
studies of cosmological N-body simulations which show that halos can maintain
triaxial shapes over cosmological timescales (Frenk et al., 1988; Dubinski &
Carlberg, 1991; Jing & Suto, 2002).

Figure 2: Axis ratios of the merged halo as a function of time in our
idealized N-body simulations. Each panel shows results from a different
simulation, with parameters as specified in the lower left corner of the
panel. (Specifically, $M_{2}/M_{1}$ gives the mass ratio of the merging halos,
$\theta$ the angle between the initial velocity vector of the halos and the
line joining their centres, $c$ is the concentration of the primary or
secondary halo as indicated by the subscript, and $\bm{a}$ is a vector
consisting of $(a_{2},a_{3})$ of the primary or secondary halo as indicated by
the subscript.) Green and blue solid lines show $a_{2}$ and $a_{3}$ measured
from the simulation, while the corresponding dotted lines show results from
our model (to be described in detail in §2.5).

Figure 2 shows axis ratios, $a_{2}$ and $a_{3}$, measured from six of our
simulations. Each panel shows the parameters of the simulation (mass ratio,
angle of approach, concentrations and axis ratios of primary and secondary
halos). Solid lines show the axis ratios measured directly from each
simulation. Dotted lines indicate expectations from our model, to be discussed
in §2.5.

The axis ratios of the simulated merger remnants show the same transitory
phase (0–2 Gyr) as was apparent in the energy tensors (Figure 1), during which
time the axis ratios show large oscillations. After this initial phase, the
axis ratios evolve more gradually, generally trending toward $a_{2}=a_{3}=1$
as the halo sphericalizes over time. An interesting exception is shown in the
lower left panel, which the axis ratios actually decrease between 10 and 16
Gyr, before increasing again. The cause of this long-timescale behaviour in
this simulation is unclear—examining the simulation snapshots suggests that
fallback of ejected particles may contribute to this effect, although it is
not clear why this does not affect other simulations.

###  2.2 Energy Tensor Evolution

We define the energy tensor of a self-gravitating system (in this work, a dark
matter halo) as

$\textbf{{E}}=\textbf{{K}}+\textbf{{W}},$ (3)

where K is the kinetic energy tensor defined as (Binney & Tremaine, 2008, eqn.
4.240b):

$K_{\mathrm{ij}}=\frac{1}{2}\int\mathrm{d}^{3}x\rho(\textbf{{x}})\overline{v_{\mathrm{i}}v_{\mathrm{j}}},$
(4)

and W is the potential energy tensor defined as (Binney & Tremaine, 2008,
2.19):

$W_{\mathrm{ij}}=-\int\mathrm{d}^{3}x\rho(\textbf{{x}})x_{\mathrm{i}}\frac{\partial\Phi}{\partial
x_{\mathrm{k}}}.$ (5)

Systems in equilibrium obey a tensor virial theorem,
$2\textbf{{K}}+\textbf{{W}}=0$ (Binney & Tremaine, 2008, eqn. 4.248). The
trace of the energy tensor is equal to the sum of the scalar potential and
kinetic energies, and is therefore a conserved quantity (in a time-independent
system). However, the energy tensor itself is not a conserved quantity.
Nevertheless, we will show in §2.5 that it is approximately333The conservation
is approximate during the merger due to interactions with the environment,
including loss of mass from the system during the merger. A similar break down
of precise conservation was also seen by Johnson et al. (2021) in the scalar
energy (see §2.4.1 for details). We will show in §3.1 that the same
corrections for this non-conservation that were found by Johnson et al. (2021)
also apply here. conserved during mergers, although over cosmological
timescales the energy tensor evolves to become more isotropic (see §2.4.2).
Therefore we proceed, for now, under the assumption that E is constant for a
given pair of merging halos and the resulting merged system.

Considering two dark matter halos, labelled by superscript (1) and (2), that
are about to merge, we can write the total energy tensor of the system,
$\textbf{{E}}^{\mathrm{(tot)}}$, as:

$\textbf{{E}}^{\mathrm{(tot)}}=\frac{\textbf{{W}}^{\mathrm{(1)}}}{2}+\frac{\textbf{{W}}^{(2)}}{2}+\textbf{{V}}^{\mathrm{(12)}}+\textbf{{K}}^{\mathrm{(orb)}},$
(6)

where we have used the tensor virial theorem to write
$\textbf{{E}}^{\mathrm{(i)}}=\textbf{{W}}^{\mathrm{(i)}}/2$ for the energy
tensor of each individual halo, and where $\textbf{{V}}^{\mathrm{(12)}}$
represent the potential energy tensor due to the gravitational interaction of
the two halos immediately prior to merging, and
$\textbf{{K}}^{\mathrm{(orb)}}$ is the kinetic energy tensor associated with
the orbital motion of the two halos immediately prior to merging.

Our choice to set $\textbf{{E}}^{\mathrm{(i)}}=\textbf{{W}}^{\mathrm{(i)}}/2$
differs from the approach taken in Johnson et al. (2021). In that work the
energy of each halo was found by summing potential and kinetic terms, with the
kinetic energy found by solving the Jeans equation in spherical symmetry. As
shown by Cole & Lacey (1996), this approach leads to virial ratios, $2K/|W|$,
in good agreement with N-body simulations, even though cosmological N-body
halos are found to deviate somewhat (10–20%) from $2K/|W|=1$ (i.e. they are
not in virial equilibrium). Johnson et al. (2021) further show that this
approach results in good agreement with the total energies of cosmological
halos.

Unfortunately, it is not possible in general to solve the Jeans equations for
triaxial systems of the type considered in this work444van de Ven et al.
(2003) find general solutions for the Jeans equation for triaxial Stäckel
models, but we are not aware of other solutions.. As such, we resort to
assuming that the (tensor) virial theorem holds, and thereby estimating the
total energy tensor from the potential energy tensor, without the need to
compute the kinetic energy tensor. Given the fact that cosmological halos are
not precisely in virial equilibrium this will introduce some error, at the
10–20% level, into our model. Given the other simplifying assumptions made we
do not consider this to be a critical problem.

The potential energy tensor, W, can be found analytically for a triaxial NFW
potential. Defining (Binney & Tremaine, 2008, eq. 2.117),

$\psi(m)=\int_{0}^{m}2R_{\mathrm{s}}\frac{\rho^{*}}{(1+\frac{m^{\prime}}{R_{\mathrm{s}}})^{2}}\mathrm{d}m^{\prime}$
(7)

where the triaxial coordinate $m^{2}$ is defined in equation (2), we find

$W^{\prime}_{ii}=-\pi^{2}a_{2}^{2}a_{3}^{2}Ga_{i}^{2}A_{i}\int_{0}^{cR_{\mathrm{s}}}[\psi(m)-\psi(cR_{\mathrm{s}})]^{2}\mathrm{d}m.$
(8)

This gives

$\textbf{{W}}^{\prime\mathrm{(i)}}=-\gamma\frac{G\tilde{M}^{2}}{R_{\mathrm{s}}}\begin{pmatrix}A_{1}&0&0\\\
0&a_{2}^{2}A_{2}&0\\\ 0&0&a_{3}^{2}A_{3}\\\ \end{pmatrix},$ (9)

(Chandrasekhar, 1987, eqn. 93) where

$\gamma=\frac{9}{4}\left[1-\frac{1}{(1+c)^{2}}-\frac{2\log(1+c)}{1+c}\right],$
(10)

and

$\tilde{M}=\frac{4\pi}{3}a_{2}a_{3}\rho^{*}R_{\mathrm{s}}^{3},$ (11)

and the virial radius is $cR_{\mathrm{s}}$. $\textbf{{W}}^{\prime}$ indicates
components of the tensor defined in the coordinate system corresponding to the
principal axes of the triaxial halo, and $\tilde{M}$ can be thought of as a
scale mass.

The $A_{\mathrm{i}}$ terms are defined as (Chandrasekhar & Lebovitz, 1962,
§II):

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\frac{2}{\sin^{3}\phi\sin^{2}\theta}[F(\theta,\phi)-E(\theta,\phi)],$
$\displaystyle A_{2}$ $\displaystyle=$
$\displaystyle\frac{2}{\sin^{3}\phi\sin^{2}\theta\cos^{2}\theta}$
$\displaystyle\,\,\,\left[E(\theta,\phi)-F(\theta,\phi)\cos^{2}\theta-\frac{a_{3}}{a_{2}}\sin^{2}\theta\sin\phi\right],$
$\displaystyle A_{3}$ $\displaystyle=$
$\displaystyle\frac{2}{\sin^{3}\phi\cos^{2}\theta}\left[\frac{a_{2}}{a_{3}}\sin\phi-E(\theta,\phi)\right],$
(12)

where $E()$ and $F()$ are the incomplete elliptic integrals

$E(\theta,\phi)=\int_{0}^{\phi}(1-\sin^{2}\theta\sin^{2}\phi^{\prime})^{1/2}\mathrm{d}\phi^{\prime},$
(13)

and

$F(\theta,\phi)=\int_{0}^{\phi}(1-\sin^{2}\theta\sin^{2}\phi^{\prime})^{-1/2}\mathrm{d}\phi^{\prime},$
(14)

with

$\sin\theta=\sqrt[]{\frac{1-a_{2}^{2}}{1-a_{3}^{2}}},\quad\mathrm{and}\quad\cos\phi=a_{3}.$
(15)

By construction, $\textbf{{W}}^{\prime\mathrm{(i)}}$ is diagonal in a
coordinate system aligned with the principal axes of the triaxial halo.
Generically, though, the two coordinate systems aligned with the principal
axes of our two halos will not coincide with each other. We therefore chose to
rotate the potential energy tensor $\textbf{{W}}^{\prime\mathrm{(2)}}$ into
the coordinate system corresponding to the principal axes of halo (1). We
assume no correlation between the principal axes of halos (1) and (2). Thus,
$\textbf{{W}}^{\prime\mathrm{(1)}}$ and $\textbf{{W}}^{\prime\mathrm{(2)}}$
are both generated as in equation 9, and we then set
$\textbf{{W}}^{\mathrm{(1)}}=\textbf{{W}}^{\prime\mathrm{(1)}}$ and
$\textbf{{W}}^{\mathrm{(2)}}=\textbf{{R}}\textbf{{W}}^{\prime\mathrm{(2)}}\textbf{{R}}^{\mathrm{T}}$,
where R is a pitch-yaw-roll rotation matrix chosen such that the principal
axes of halo (2) are isotropically distributed relative to those of halo (1).

We now define a coordinate system with the first axis aligned with the vector
connecting the centres of our two halos, and the remaining axes perpendicular
to this with a rotation chosen uniformly at random. In this coordinate system
the interaction tensor, $\textbf{{V}}^{\mathrm{(12)}}$, can be written as

$V^{\prime}_{\mathrm{ij}}=-\frac{GM_{1}M_{2}}{c^{\mathrm{(1)}}R_{\mathrm{s}}^{\mathrm{(1)}}}\begin{pmatrix}1&0&0\\\
0&0&0\\\ 0&0&0\\\ \end{pmatrix},$ (16)

where
$R_{\mathrm{v}}^{\mathrm{(1)}}=c^{\mathrm{(1)}}R_{\mathrm{s}}^{\mathrm{(1)}}$
is both the virial radius of halo 1, and the distance between halos at time of
merging. In this same coordinate system, the components of the kinetic energy
tensor, $\textbf{{K}}^{\mathrm{(orb)}}$, are trivially related to the radial
and tangential velocities of halo (2). We must then rotate the components of
these tensors into the coordinate system aligned with the principal axes of
halo (1). The components of the required rotation matrix, Q, can be found from
the dot products of the basis vectors of our primed coordinate system,
$\hat{\textbf{{x}}}^{\prime}$ and the coordinate system of halo (1),
$\hat{\textbf{{x}}}^{\mathrm{(1)}}$,

$Q_{\mathrm{ij}}\equiv\langle\hat{x}_{\mathrm{i}}|\hat{x}^{\prime}_{\mathrm{j}}\rangle,$
(17)

which can be computed explicitly.

###  2.3 Applications to Merger Trees

In §2.2, we described the effect of single merger between two halos on the
energy tensor of the merged halo. However, in the formation process of any
current halo, a very large number of mergers will have taken place. To account
for this, we create merger trees, which represent a ledger of every merging of
two progenitor halos. The merger tree records the mass of each progenitor
halo, and the time and ordering of every merging event between progenitor
halos. We can then apply the model of §2.2 to each merging event in order,
recording the resulting values of $a_{2}$ and $a_{3}$ along with the full
energy tensor for each halo. This energy tensor is then used in the next
merger event involving the resulting merged halo. In this way, our model is
applied recursively through the tree, until we reach the final, $z=0$ halo and
compute its triaxial shape.

Fully-specifying this application of our model to a merger tree requires
additionally specifying how to assign orbits to merging halos, how to assign
energy tensors to halos in the tree that have no progenitors, and how to deal
with subresolution accretion. These aspects are described in the remainder of
this section.

In §2.4 we will discuss additional physics beyond that described in §2.2,
including possible non-conservation of the energy tensor through
sphericalization and interactions with the environment.

####  2.3.1 Merging Halo Orbits

As in Benson et al. (2020) and Johnson et al. (2021), we assign orbital
velocities to merging halo pairs by drawing values at random from the joint
distribution of these velocities measured from cosmological N-body
simulations. Specifically, we use the distribution functions measured by Jiang
et al. (2015). As in Benson et al. (2020) we introduce a weak correlation
between the orbital angular momentum of merging secondary halos and the spin
vector of the primary halo, $\bm{\lambda}_{\mathrm{P}}$, such that the angle
between these two vectors is distributed as

$P(\cos\theta)=\frac{1}{2}(1+\eta|\lambda_{\mathrm{P}}|\cos\theta),$ (18)

where the parameter $\eta$ controls the strength of the correlation. Benson et
al. (2020) found a best fit value of $\eta=(6.6^{+11.3}_{-4.5})\times 10^{-1}$
when constraining their model to match the distribution of halo spin
parameters measured in N-body simulations. We use a value of $\eta=0.16$
(consistent with the value from Benson et al. 2020 within the uncertainties)
in this work based on an updated fit of the Benson et al. (2020) model to data
from the MultiDark Planck simulation suite (Riebe et al., 2013; Klypin et al.,
2016) which constrains $\eta<0.65$ at 65% confidence with a best-fit value of
$\eta=0.16$.

####  2.3.2 Initial Conditions

Due to the finite resolution of our merger trees there are always some halos
which have no progenitors. Therefore, for these halos we can not determine
their energy tensor based on an application of equation (6). Instead, we
assume that these halos are spherically-symmetric, and set the scale radius of
the NFW profile using the concentration model of Diemer & Joyce (2019). The
merger trees used in this work have a mass resolution equal to $10^{-4}$ times
the mass of the final, $z=0$ halo. We have checked that this resolution is
sufficient that the statistics of the $z=0$ halo triaxial shapes are
insensitive to our choice of initial conditions in these unresolved progenitor
halos.

####  2.3.3 Subresolution accretion

The finite resolution of merger trees also means that not all growth in halo
mass is due to merging events captured in the tree—there is some sub-
resolution accretion which also leads to mass growth.

To account for this subresolution accretion we simply assume that the energy
tensor of the halo grows in proportion to the characteristic virial energy of
the primary halo,
$E_{\mathrm{vir}}\propto-\mathrm{G}M^{2}_{\mathrm{vir}}/R_{\mathrm{vir}}$,
where $M_{\mathrm{vir}}$ and $R_{\mathrm{vir}}$ are the virial radius and mass
of the halo respectively. Given the resolution of our merger trees ($10^{-4}$
times the $z=0$ halo mass) we find that this subresolution accretion has only
minimal effects on our results.

###  2.4 Beyond Energy Tensor Conservation

So far, we have considered the idealized case in which the energy tensor is
perfectly conserved during, and between merger events. In this subsection we
explore processes which break this assumption, and extend our model to account
for them. These additions do not change how our model is applied to a merger
tree, as outlined in §2.3.

####  2.4.1 Corrections due to interactions with the environment

The model described in §2.2 correctly describes the merger of two smooth, NFW,
isolated halos if all components of the energy tensor are conserved, and their
gravitational interaction can be well-approximated by the interaction of two
point masses. However, as was discussed in §1, this is not the case. Even in
the prior work of Benson et al. (2020) and Johnson et al. (2021), which worked
with quantities that are precisely conserved in isolated systems (angular
momentum and energy), we found the need to introduce terms breaking this
conservation. First, cosmological halos, are not perfectly smooth, nor
perfectly described by a triaxial NFW profile, and nor is their interaction
precisely that of two point masses. These are approximations made in our
model, and we expect some order unity differences between our model
predictions and direct N-body simulations as a result. Furthermore, in
practice, pairs of merging halos are not isolated systems—they interact with
their environment, exchanging angular momentum, energy, and mass. We briefly
review the corrections introduced by Johnson et al. (2021) to account for this
below.

These corrections will be included both when comparing our model to idealized
halo merger simulations in §2.5, and when calibrating our model to
cosmological simulations in §3.1. While the “environment” of a merging pair of
halos differs in these two cases—being a vacuum for the case of the idealized
mergers, and a cosmological distribution of matter and halos in the
cosmological case—interactions with the environment happen nevertheless.
Specifically, even for the idealized simulations particles can be scattered
onto unbound (or close to unbound) orbits, and lost from the system to the
environment, carrying away energy. For this reason, we apply the model
described below for non-conservation effects in both cases, although we allow
for the possibility that the quantitative details will differ by retuning the
parameters of the model in each case.

In addition to these environmental exchanges, in the current work we must
account for the fact that the energy tensor is not conserved even in an
isolated, post-merger system, for both idealized and cosmological
calculations. We describe how this is accounted for in §2.4.2.

To incorporate the effects of interactions with the environment we employ the
same model as Johnson et al. (2021). Specifically, we modify equation (6) to
become

$\textbf{{E}}^{\mathrm{(tot)}}=\frac{\textbf{{W}}^{\mathrm{(1)}}}{2}+\left[\frac{\textbf{{W}}^{(2)}}{2}+\textbf{{V}}^{\mathrm{(12)}}+\textbf{{K}}^{\mathrm{(orb)}}\right]\left(1+\frac{M_{2}}{M_{1}}\right)^{-\gamma}(1+b).$
(19)

The new term multiplies the kinetic, potential, and interaction energy of halo
(2) but does not affect halo (1). This is because it is only the _added_
energy of the merger that we are modifying. This term contains two free
parameters, $\gamma$ and $b$, both of which we set to their best-fit values
from Johnson et al. (2021): $\gamma=1.518$ and $b=0.673$.

The form of this correction is somewhat arbitrary. In Johnson et al. (2021),
we hypothesized that the failure of (scalar) energy conservation would depend
on the mass ratio of the merging halos as the mass fraction which becomes
unbound and lost from the system (thereby directly leading to a failure of
energy conservation in the final halo) is strongly mass-dependent. This
motivates the form of the $(1+M_{2}/M_{1})^{\gamma}$ term—in the limit of
$M_{2}\ll M_{1}$ it tends to $1$ (i.e. perfect conservation) while allowing
for either an increase or decrease in the energy for cases where $M_{2}\sim
M_{1}$. The precise form was chosen as a simple parameterization which
captures this desired behaviour. Any mass-independent failure of conservation
is then accounted for by the $(1+b)$ term. Here we implicitly assume that this
mass-independent failure of conservation is also independent of all other halo
properties, and so this term is a constant. The justification for this choice
is that it is the simplest possible assumption, and allowed for a successful
match to N-body measurements of halo concentrations in Johnson et al. (2021).

####  2.4.2 Sphericalization

Assuming that sphericalization of the shape of our halos is driven by
differential precession, we can estimate the timescale for this process. Erkal
et al. (2016) derive that the precession frequency in a flattened potential is
given by

$\omega=\Omega(1-1/q),$ (20)

where

$\Omega=\sqrt{\frac{\mathrm{G}M_{\mathrm{vir}}}{r_{\mathrm{vir}}^{3}}}$ (21)

is the characteristic orbital frequency at the virial radius, and we estimate
the flattening parameter $q$ as

$q=(a_{2}a_{3})^{-1/4}.$ (22)

(Note that Erkal et al. (2016) define the flattening in terms of the
potential, rather than in terms of the density profile as in this work. The
above expression accounts for the fact that the potential is less flattened
than the density) Between each merger event in our merger tree we allow the
eigenvalues of the energy tensor to evolve according to:

$\dot{\lambda}_{\mathrm{i}}=-\alpha\omega(\lambda_{\mathrm{i}}-\bar{\lambda}),$
(23)

where $\bar{\lambda}$ is the mean of the eigenvalues, and $\alpha$ is a
dimensionless parameter. This approach ensures that the total energy
($\sum\lambda_{\mathrm{i}}\equiv 3\bar{\lambda}$) is unchanging while driving
the system toward a spherical shape.

We note that the form of equation (23 implies that halos will be driven toward
a spherical shape over cosmological timescales, unless $\alpha\ll 1$. This may
be in conflict with the results of prior work, such as that of Vasiliev &
Athanassoula (2012) discussed in §2.1.

However, because of the halo shape dependence in the precession timescale
(equations 20 and 22) the timescale for sphericalization becomes longer as the
halo becomes more spherical—as such, complete sphericalization is less rapid
than might be expected from the dynamical timescale and value of $\alpha$
alone. Second, in the cosmological case of most interest here, mergers will
continuously drive the halo away from spherical, such that what matters most
is the application of equation (23) during the short intervals between
frequent merger events—the long term behaviour implied by equation (23) is
less important to determining the results of our model.

Nevertheless, it should be kept in mind that the above model for
sphericalization is clearly an oversimplification. Based on the idealized
merger simulation results of §2.5 applying this model over times of 10 Gyr or
more (with no ongoing merger activity to drive the halo shape evolution) is
unlikely to be accurate.

###  2.5 Comparison with N-body Simulations

Having developed our model for the evolution of halo axis ratios after a
halo–halo merger, we can examine the accuracy of this model by comparing to
the idealized N-body simulations presented in §2.1.

In Figure 1, which shows the evolution of the components of the energy tensor
in two of these idealized simulations, the dotted lines show the predictions
of our model. Specifically, we used the eigenvalues of the energy tensor
measured in the initial conditions of each simulation as the starting point
for our model calculation. We then apply the model described in equation (23)
to account for the effects of sphericalization as the simulation progresses.
This gives the eigenvalues of the energy tensor at each timestep of the
simulation. We then solve equation (9) at each timestep to find the
corresponding axis ratios predicted by our model.

Our model has three free parameters, $(b,\gamma,\alpha)$. We choose the
optimal values of these parameters to minimize the median absolute deviation
between our model prediction and the simulation-measured axis ratios across
all of our idealized simulations, excluding the first 5 Gyr of each simulation
to remove the transitory phase where the eigenvalues of the energy tensor and
axis ratios show large fluctuations. We find optimal values of
$(b,\gamma,\alpha)=(0.295,1.175,2.145)$. We note that these values are optimal
for these idealized simulations, but will not necessarily be the optimal
choice for cosmological halos. Therefore, we will re-visit tuning of these
parameters when we apply our model to cosmological cases. For the values of
$(b,\gamma,\alpha)$ given above, the median absolute deviation between
simulation and model axis ratios ranges between 0.010 and 0.125, with a mean
of 0.042. Therefore, our model expectations are within $\Delta a\approx 0.04$
of the simulation axis ratios 50% of the time.

We begin by considering the left panel of Figure 1, which shows results from
simulation “d”, a 10:1 mass ratio merger with angle of approach $\theta=\pi/4$
between spherical halos. In this case, our model accurately captures the
evolution of the components of the energy tensor with time. In particular, the
non-conservation of these components due to sphericalization is clearly seen
in this figure. Recall that the horizontal dashed lines indicate the initial
values of each component, and so show the expectation if the energy tensor
were perfectly conserved, with no sphericalization. The N-body results (solid
lines) clearly deviate from this expectation, and this deviation is accurately
captured by our model. We find that the majority of our simulations show
similarly good agreement with our model.

There are, however, some exceptions. In the right panel of Figure 1 we show
results from simulation “g”, a radial 10:1 merger between two oblate halos. In
these case some sphericalization can once again be seen. However, here the
sphericalization occurs more slowly than predicted by our model, meaning that
at late times our model will predict a more spherical halo than was found in
the N-body simulation. This demonstrates that, while sphericalization seems to
occur in all of our simulations, our simple model (§2.4.2) does not fully
capture the details of this process, and so is an avenue for further
exploration.

To judge the validity of our model in matching the simulation results for the
energy tensor evolution we compute median absolute deviations (MADs) between
the diagonal components of the tensor as measured from the simulation and as
predicted by our model, normalized to $|E|/3$. For the optimal values of
$(b,\gamma,\alpha)$ given above we find that the MAD ranges from $0.03$ in the
best case (simulation “b” in Table 1), to $0.22$ in the worst case (simulation
“h”). Two thirds of the simulations (8 out of 12) have a MAD of $0.10$ or
smaller. In the majority of cases, our model assumptions (an energy tensor
that is conserved through the merger and which then undergoes slow
sphericalization) result in a match to the simulation results to within 10%.

We next examine our model predictions for the axis ratios directly. Figure 2
shows axis ratios measured from six of our idealized simulations as solid
lines, with dotted lines indicating the expectation from our model. The panels
of Figure 2 are arrange in order of increasing maximum absolute deviation
between the simulation and model axis ratios (from left to right, top to
bottom).

The upper left panel shows results for a 10:1 radial merger between two
spherical halos of differing concentrations. The high degree of symmetry in
this case implies that $a_{2}=a_{3}$ should be expected. This is correctly
predicted by our model, and is approximately true in the simulation also.
After the initial transitory phase our model accurately matches the results of
the simulation, including the gradual sphericalization. The middle left panel
shows a similar result, in this case for a 5:1 mass ratio and an oblate
secondary halo. Again our model accurate matches the N-body axis ratios and
the slow sphericalization up to around 15 Gyr. After this time the simulated
halo begins to become less spherical, causing a deviation from our model
expectation. As discussed in §2.1 the cause of this late-time behaviour
remains unclear, although examination of the simulation snapshots suggests
that re-accretion of ejected particles may play a role.

The top and middle right panels show cases with a lower degree of symmetry,
due to an angle of approach of $\theta=\pi/4$. The top right panel shows a
10:1 mass ratio merger between spherical halos, while the middle right panel
shows a 3:1 mass ratio merger between two triaxial halos. In these cases our
model is slightly offset, overestimating $a_{2}$ and $a_{3}$, but nevertheless
captures the rate of sphericalization, and the approximate difference between
the axis lengths, $a_{2}-a_{3}$.

Finally, in the bottom row of Figure 2 we show two cases where our model fails
more seriously. In the left panel we examine a 10:1 radial merger between two
prolate halos. Our model performs very well in matching the N-body results up
until a time of around 10 Gyr. After this time the simulation-measured axis-
ratios decreases for around 6 Gyr before beginning to climb again. Such
behaviour is not expected in our model. Once again, the cause of this
behaviour is not clear to us, but may be related to re-accretion of ejected
particles. In the bottom right panel we show a 10:1 radial merger between two
oblate halos. In this case sphericalization seems to cease after around 7 Gyr,
resulting in our model becoming a progressively worse match to the simulation
results.

As a further test of our sphericalization model we set $\alpha=0$ (thereby
resulting in no sphericalization), recompute the predictions from our model,
and compare these new, non-sphericalized predictions to the results of the
N-body simulations. We find that the median absolute deviations between
simulated and model axis ratios increase from the range 0.010–0.125 for
$\alpha=2.145$ to the range 0.048–0.463 for $\alpha=0$. Thus, ignoring
sphericalization results in discrepancies between our model and the
simulations increasing by a factor of around 4—indicating that
sphericalization is an important component of our model.

In summary, our model captures the axis ratios and their evolution during
sphericalization measured in N-body simulations in a majority of cases with
reasonable accuracy (typical median absolute deviations of $\Delta a\approx
0.04$). The most common failure occurs at late times (10 Gyr or more after the
merger) where some simulated halos show a a trend of becoming less spherical—a
behaviour not captured in our model. In a minority of simulations the
sphericalization process appears to be significantly slower than assumed in
our model, resulting in the degree of triaxiality of the halo being
underestimated by our model. Overall, while imperfect, this level of agreement
is sufficient to allow us to make predictions for the statistical properties
of axis ratios in ensembles of mergers with at least moderate accuracy.

###  2.6 Calibration and Validation

To calibrate and validate our model we compare to the results reported by
Bonamigo et al. (2015) from the Millennium XXL simulation, and by Schneider et
al. (2012) from the Millennium I & II simulations. We generate a distributions
of triaxial shapes by constructing 10,000 merger trees for halos with $z=0$
masses of $10^{10}$, $10^{11}$, $10^{13}$, $10^{13}$, $10^{14}$,
$10^{15}\mathrm{M}_{\odot}$ utilizing the approach of Parkinson et al. (2008)
with best-fit parameters given by Benson et al. (2019). In constructing these
trees we match cosmological and power spectrum parameters to the Millennium
XXL simulation (Angulo et al., 2012).

Each tree provides virial masses and radii for all halos, along with the
cosmic time at which each halo exists, orbital velocities for each merger, and
progenitor-descendent relationships between halos. For each tree, we
initialize scale radii and shapes of all halos with no progenitors (i.e. those
just above the resolution limit of $10^{-4}M_{0}$, where $M_{0}$ is the mass
of the $z=0$ halo) as described in §2.3.2. The energy tensor of each of these
initial halos is then computed.

Working forward in time through the tree, at each merger event in the tree we
compute the energy tensor of the merged halo using equation (19), accounting
for any accreted mass as described in §2.3.3. The axis ratios,
$a_{\mathrm{i}}$, are then determined from the eigenvalues of the energy
tensor555This must be done numerically as the potential energy tensor depends
on the $a_{\mathrm{i}}$ in a non-trivial way. In practice we tabulate the
relation between the eigenvalues of the potential energy tensor and the axis
ratios to allow rapid solution.. Between each merger event in the tree we
apply our sphericalization procedure as defined in equation (23). At the end
of this process we have determined the energy tensor, and axis ratios for
every halo in the tree. For this work we examine only the properties of the
final, $z=0$ halo in the tree, which we have verified are converged with
respect to the tree mass resolution.

Bonamigo et al. (2015) and Schneider et al. (2012) both excluded unrelaxed
halos with a large offset (5% and 7% of the virial radius respectively)
between the position of the most bound particle and the centre of mass. We
cannot apply precisely this selection to our halos, but instead, when
comparing results with Bonamigo et al. (2015) and Schneider et al. (2012), we
exclude any tree which experienced a major (mass ratio of 4:1 or greater)
merger within the last 1 Gyr of evolution as these late merging halos may not
have fully virialized by $z=0$. Other cutoff times between 0.5 and 2 Gyrs were
tested, but the results were similar in all cases.

##  3 Results

We begin, in §3.1, by constraining the value of our sphericalization
parameter, $\alpha$, by matching to the distribution of $a_{3}$ found in
N-body simulations. With the calibrated model we then examine the resulting
behaviour for a single tree and compare the distribution of halo shapes with
those found in N-body simulations.

###  3.1 Calibration

We find that, for $\alpha\gtrapprox 10$, all $a_{i}\approx 1$, that is halos
are driven to be almost perfectly spherical, confirming our expectation that
$\alpha\sim 1$ will be required. To constrain the value of $\alpha$ we tune it
to find the best match to the distribution of $a_{3}$ measured from
cosmological N-body simulations by Bonamigo et al. (2015). Specifically, we
applied our model to a sample of 10,000 merger trees, each with $z=0$ halo
mass of $10^{14}\mathrm{M}_{\odot}$, using values of $\alpha$ ranging from 1.5
to 2.5 chosen based on a preliminary examination of the resulting $a_{3}$
distribution. For each value of $\alpha$ we computed the median $a_{3}$ across
our 10,000 trees and compared computed a goodness of fit measure
$\chi^{2}=[\bar{a}_{3}(\alpha)-\bar{a}_{3}^{*}]^{2}$ where
$\bar{a}_{3}(\alpha)$ is the median value of $a_{3}$ computed from our model
for a given $\alpha$, and $\bar{a}_{3}^{*}=0.53$ is the median $a_{3}$
measured by Bonamigo et al. (2015) for halos of this mass. We then seek the
value of $\alpha$ that minimizes $\chi^{2}$. We find an optimal value of
$\alpha=1.75$.

In the remainder of this paper we therefore use
$(b,\gamma,\alpha)=(0.673,1.518,1.75)$.

###  3.2 Individual Tree Behavior

Figure 3: The evolution with time of the axis lengths, $a_{i}$, in the main
branch of a single, representative merger tree of mass
$10^{14}\mathrm{M}_{\odot}$ at $z=0$ as described by our model. Note that
$a_{1}=1$ by definition. Major merger events are clearly seen as sharp
declines in $a_{2}$ and $a_{3}$, which then slowly drift toward 1 as a result
of sphericalization and multiple, smaller mergers.

Figure 3 shows the evolution of the axial ratios along the main branch of a
single, representative tree of mass $10^{14}\mathrm{M}_{\odot}$ at $z=0$. By
definition, $a_{1}$ is always 1—because the coordinate system is rotated after
every merger to align with the principle axes of the halo, $a_{1}$ is always
along the largest axis and should not be thought of as a particular direction
in space. The other axis lengths, $a_{2}$ and $a_{3}$, show large
variations—these are strongly correlated both because an individual merger
event will affect both axes, and because they are both normalized to the
principle axis length. After each major merger, the halo becomes very prolate,
but sphericalization and the collective effects of many smaller mergers drive
the halo back toward a more spherical configuration. Major mergers become more
rare as the halo grows in mass and the timescale for growth (and, therefore,
for changes in triaxial shape) increases as the halo evolves due to the
decreasing density (and dynamical timescale) of the halos.

###  3.3 Unconditioned Distributions of Axis Ratios

|
---|---

Figure 4: Probability distribution functions for the axis ratios $a_{2}$ (left
panel) and $a_{3}$ (right panel) derived by applying our model to 10,000
realizations of merger trees with mass $10^{14}\mathrm{M}_{\odot}$ at $z=0$.
In the left panel, for $a_{2}$, we overlay the corresponding result measured
from cosmological N-body simulations by Schneider et al. (2012). In the right
panel, for $a_{3}$, we overlay N-body results from Bonamigo et al. (2015) and
Schneider et al. (2012)

.

Figure 5: Probability distribution function for the triaxiality parameter
($T\equiv[1-a_{2}^{2}]/[1-a_{3}^{2}]$) derived by applying our model to 10,000
realizations of merger trees with mass $10^{14}\mathrm{M}_{\odot}$ at $z=0$.
N-body simulation results are taken from Schneider et al. (2012).

Figure 4 shows the distributions of $a_{2}$ and $a_{3}$ values (left and right
panels respectively) derived by applying our model to 10,000 realizations of
merger trees with mass $10^{14}\mathrm{M}_{\odot}$ at $z=0$. For the $a_{2}$
distribution, we overlay the corresponding result from Schneider et al.
(2012)666Schneider et al. (2012) uses a halo mass definition based on a
density threshold of 200 times critical density, instead of the spherical
collapse-motivated threshold used in this work and by Bonamigo et al. (2015).
In making the comparison in Figure 4 we convert from this 200 times critical
mass definition to a spherical collapse-motivated definition assuming NFW
density profiles and the concentration model of Diemer & Joyce (2019).. For
the $a_{3}$ distribution, we overlay N-body results from both Schneider et al.
(2012) and Bonamigo et al. (2015). The agreement in the median value of
$a_{3}$ with this N-body measurement is a result of our use of this statistic
to tune the value of the sphericalization parameter, $\alpha$. Our predicted
distributions for $a_{2}$ and $a_{3}$ show clear peaks with extended tails to
low values. Comparing the distribution of $a_{3}$ with that from Bonamigo et
al. (2015) we find that our model predicts a broader distribution of
$a_{3}$—in particular is has many more halos with very small values of
$a_{3}\lessapprox 0.3$. Our distribution of $a_{3}$ has 25 and 75 percentiles
of $0.41$ and $0.64$, while that measured by Bonamigo et al. (2015) has $0.47$
and $0.61$. Schneider et al. (2012) found an $[25,50,75]$ percentiles of
$[0.46,0.52,0.58]$, which is narrower than Bonamigo et al. (2015)’s result,
but matches their median of 0.54 quite well. This discrepancy will be
discussed further in sections 4 and 5. For the $a_{2}$ distribution, we found
a median of 0.69 as opposed to the median of $0.67$ from the N-body
simulations of Schneider et al. (2012). Our 25 and 75 percentiles are $0.57$
and $0.79$ respectively, while Schneider et al. (2012) found values of $0.58$
and $0.78$. Therefore, this work overestimates the median value slightly and
over predicts low $a_{2}$ values, but otherwise agrees quite well with
simulation. It should also be noted that the distribution of $a_{2}$ was
measured from a different simulation to which our model has not been tuned.
Given that there are significant differences between axis ratios measured from
various simulations (as will be discussed in §3.4), this agreement is evidence
for the robustness of our model.

Finally, figure 5 shows the distribution of the triaxiality,

$T=\frac{1-a_{2}^{2}}{1-a_{3}^{2}},$ (24)

of this model and compares it with the same parameter measured from N-body
simulations by Schneider et al. (2012). The overall shape of the predicted
distribution is similar to that found by Schneider et al. (2012), although the
peak is offset to slightly smaller values of $T$.

###  3.4 Conditioned Distributions and Mass Dependence

Figure 6: The median value of $a_{2}$, $\bar{a}_{2}$, conditioned on $a_{3}$
for halos of mass $10^{14}\mathrm{M}_{\odot}$ at $z=0$. Blue points show
results from this work found by subdividing our sample of 10,000 merger trees
into bins of $a_{3}$. Black points show the medians reported by (Bonamigo et
al., 2015, their Figure 12) with error bars indicating the inter-quartile
range in both data sets. Bonamigo et al. (2015) show results for a variety of
mass intervals. Here we show their result for halo masses between $10^{13.4}$
and $10^{14.2}\mathrm{M}_{\odot}$.

In Figure 6 we examine the correlation between axis ratios predicted by our
model by considering the median value of $a_{2}$ conditioned on the value of
$a_{3}$, $\bar{a}_{2}(a_{3})$. This is shown by the blue points, which were
computed by subdividing our sample of 10,000 merger trees into bins of
$a_{3}$, and then finding the median $a_{2}$ in each bin. For this calculation
we again used trees with $z=0$ halo mass of $10^{14}\mathrm{M}_{\odot}$. For
comparison we show N-body results from Bonamigo et al. (2015) as the black
points. Bonamigo et al. (2015) show that the relation is largely independent
of halo mass—nevertheless, we plot their result for halos in the mass range
$10^{13.4}$ to $10^{14.2}\mathrm{M}_{\odot}$ to most closely match our halo
masses. Our model matches the N-body results of Bonamigo et al. (2015) very
well for $a_{3}>0.3$. At lower value it somewhat overestimates the value of
$\bar{a}_{2}(a_{3})$.

Figure 7: The median value of $a_{3}$, $\bar{a}_{3}$, as a function of $z=0$
halo mass. Black triangles show results from this work found by applying our
model to 10,000 merger tree realizations in each of six different halo mass
bins. The various N-body curves are taken from Bonamigo et al. (2015; Figure
13).

Finally, in Figure 7, we plot the median $a_{3}$, $\bar{a}_{3}$, as a function
of $z=0$ halo mass. Black triangles show results from this work found by
applying our model to 10,000 merger tree realizations in each of six different
halo mass bins. Cyan squares show results from Bonamigo et al. (2015) while
lines show a compilation of other N-body results (Bonamigo et al., 2015;
Allgood et al., 2006; Bett et al., 2007; Muñoz-Cuartas et al., 2011; Schneider
et al., 2012; Despali et al., 2014; Jing & Suto, 2002) gathered by Bonamigo et
al. (2015). Statistical uncertainties on the results from this work are at the
sub-percent level—small compared the the difference between N-body results
seen here. While the N-body results generally have small statistical
uncertainty777Bonamigo et al. (2015) do not show error bars on their points,
but the bin-to-bin variations suggest statistical variations of no more than
$\pm 0.02$ in $\bar{a}_{3}$., there is a significant variation between the
various N-body-derived results888This variation is likely due, at least in
part, to the poor convergence properties of halo shapes in N-body simulations
as shown, for example, by Mansfield & Avestruz (2021; Figure 3)—with
significant differences between different simulations persisting even in halos
resolved with $10^{5}$ particles.. Nevertheless, all show the same general
trend of decreasing $\bar{a}_{3}$ with increasing mass. Our model predicts a
similar trend and lies within the spread of the N-body results.

##  4 Discussion

Our model for the triaxial shape of dark matter halos functions by tracking
the energy tensor of a halo during each merger event in its formation history,
under the assumption that the energy tensor is approximately conserved. Minor
corrections, which break perfect conservation of the energy tensor, are
applied during each merger, and a slow sphericalization is allowed to occur
between each merger event. We verify our model by comparing its expectations
to idealized N-body simulations of merging, triaxial halo pairs, finding that
in the majority of cases it matches the measured axis ratios with a median
absolute deviation of $\Delta\approx 0.04$.

Applying this model to a merger tree, we are able to track the energy tensor
and axis length ratios along the main branch of the tree999We in fact track
these along all branches, although the results may be affected by resolution
in lower mass branches., and record how the triaxial shape of halos changes
over time.

We find that statistics of halo triaxial shape are very sensitive to the
sphericalization parameter, $\alpha$, which controls the timescale for
sphericalization. A value of $\alpha\sim 1$ is found to result in a good match
to the properties of triaxial halos found in cosmological N-body simulations,
suggesting that the orbital precession frequency (equation 20) is the
appropriate timescale for sphericalization. However, some of our idealized
N-body experiments show more complex sphericalization behaviour. This is an
area where improved understanding of the relaxation of triaxial halos could
significantly improve our model.

Figure 4 demonstrates that our model produces reasonable agreement with the
distribution of axis length ratios measured in cosmological N-body simulation.
The medians of the distributions from our model and the N-body simulation
match quite well (once we have tuned the sphericalization parameter,
$\alpha$), and the general shape is similar. However, as one can see from the
$a_{2}$ and $a_{3}$ histograms (Figure 4), our model produces a significantly
more broad distribution than that from the N-body simulation. This is
particularly apparent for very low $a_{3}$ values. As can be seen in Figure 3,
low $a_{2}$ and $a_{3}$ values occur after a major merger. Therefore, it is
possible that our model does not fully capture the dynamics of the largest
mergers. For example, the $(1+M_{2}/M_{1})^{-\gamma}$ term in equation (19),
introduced to account for the possibility of mass/energy ejection from the
system during the merger, may be insufficient to capture the behaviour of the
most massive mergers. Similarly, it is possible that the underlying assumption
that the energy tensor is approximately conserved may simply fail for
cosmological major mergers (although it holds well for mergers of idealized
halos as discussed in §2.4.2). A second possible cause for the excess of low
$a_{2}$ and $a_{3}$ values is that our simple criterion for excluding
unvirialized halos (based on the time since the last major merger) is
insufficient to match in detail the relaxation criteria applied by Bonamigo et
al. (2015) to their sample. Other possible causes for this discrepancy include
the limitations of the simple “$(b,\gamma)$” model101010We have used this
model unchanged from the form used by Johnson et al. (2021) who applied this
to a similar model for halo concentrations based on the scalar energy of
halos. Since the trace of the energy tensor is the scalar energy, we expect
that applying this same model here should ensure that we similarly match the
behaviour of the trace of the energy tensor, but there is no guarantee that it
will give the correct behaviour in a more general sense. for non-conservation
of the energy tensor in mergers, and our simple model for
sphericalization111111As already discussed, it is clear that the precession
frequency given by equation (20) is not the correct timescale for
sphericalization. It is also possible that equation (23) is not the correct
functional form for the evolution of the eigenvalues of the energy tensor
during sphericalization., given by equation (23). As can be seen in Figure 1,
the dynamics of mixing are quite complex, and sphericalization does not
necessarily turn on immediately in the same way as Figure 3.

These limitations do not invalidate the underlying intent of the model, namely
to describe the qualitative evolution of triaxial halo shapes measured in
N-body simulations, and to provide an approximate quantitative description of
them. They do, however, suggest a rich variety of improvements one could
explore in future works.

To explore the predictive power of our model, once optimized to match
unconditioned distributions of axis length ratios, we have considered the
triaxiality parameter (Figure 5), conditioned distributions (Figure 6), and
mass-dependent distributions (Figure 7). The distribution of triaxiality
parameters is of the same shape as that found in N-body simulations, although
offset to somewhat lower values, while we find that our model’s predictions
for conditioned distributions and mass dependence to both agree extremely well
with measurements from N-body simulations. The median $a_{2}$ conditioned on
$a_{3}$ shown in Figure 6 suggests that our model correctly captures the
correlated evolution of $a_{2}$, $a_{3}$ during mergers—except perhaps at
$a_{3}<0.3$ where we slightly overestimate $\bar{a}_{2}(a_{3})$. Given the
heavy tail of low $a_{3}$ values in 4, this is perhaps due to the
overpopulation of low $a_{3}$ halos in our model. Mass-dependence, which can
only arise from the differences encoded into the structure of the merger
trees, suggests that our model is correctly capturing the physics that drives
triaxial evolution during halo assembly. Improvements such as those discussed
above could likely improve the agreement at $M=10^{15}\mathrm{M}_{\odot}$, but
given the scatter between different N-body determinations at lower mass
scales, we are hesitant to overfit to a single data set.

##  5 Conclusion

The triaxial shapes of dark matter halos are important due to their
consequences for the orbital evolution of subhalos, the strength of
gravitational lensing produced by a halo, and many other factors. This work
builds on two previous papers (Benson et al., 2020; Johnson et al., 2021)
which track halo spin and concentration using a similar approach. In this
paper we expand their random-walk methodology to gain insight into the
physical mechanisms governing the evolution of triaxial halo shapes. We find
that the evolution of these triaxial halo shapes is driven by a combination of
mergers and sphericalization. After calibration of our model to match the
median $a_{3}$ in cluster halos, we find that these features are sufficient to
explain conditioned distributions of axis ratios and mass dependence measured
in cosmological N-body simulations, although the model underpredicts the width
of the distributions of $a_{2}$ and $a_{3}$. While our model does not produce
high-precision matches the N-body measurements, this was not the goal.
Instead, we aimed to explore whether thinking about triaxial halo shape
evolution in this way is a productive avenue of exploration. Future
improvements, if desired, should focus on generalizing the correction factors
(for environmental interactions, and sphericalization) employed here.

There are a number of other avenues that could be explored in future work.
Firstly, this work connects triaxial shape evolution directly to the growth
and merger history of a halo. Connecting this back to the previous work with
this framework would allow for predictions for the joint distributions of halo
spin, concentration, and triaxial shape and how these grow together. Second,
our model only considers the global triaxial shape of halos, while radial
variations are known to exist. For example Vera-Ciro et al. (2011) found that
each radial shell of a halo retains a memory of its shape at the time at which
it was assembled. Our model applies to halos, but could be adapted for
subhalos. For these structures there are additional drivers of triaxial shape
evolution, notably tidal forces from their host halo, which can cause subhalos
to become more prolate (Kuhlen et al., 2007), but likely also lead to some
additional sphericalization due to tidal heating effects. Finally, it may be
interesting to extend our model, currently developed for cold dark matter, to
cases where dark matter has self-interactions which is known to affect halo
triaxiality (Peter et al., 2013; Tulin & Yu, 2018). Such self-interactions
could lead to more rapid thermalization and sphericalization which could be
incorporated into our sphericalization model.

## Acknowledgements

We thank Xiaolong Du, Fangzhou Jiang, Ethan Nadler, and Shengqi Yang for
helpful discussions and comments on this work. The CosmoSim database used in
this paper is a service by the Leibniz-Institute for Astrophysics Potsdam
(AIP). The MultiDark database was developed in cooperation with the Spanish
MultiDark Consolider Project CSD2009-00064. The authors gratefully acknowledge
the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the
Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu) for
funding the MultiDark simulation project by providing computing time on the
GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
  * Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
  * Benson et al. (2019) Benson A. J., Ludlow A., Cole S., 2019, MNRAS, 485, 5010
  * Benson et al. (2020) Benson A., Behrens C., Lu Y., 2020, Monthly Notices of the Royal Astronomical Society, 496, 3371
  * Bett et al. (2007) Bett P., Eke V., Frenk C. S., Jenkins A., Helly J., Navarro J., 2007, MNRAS, 376, 215
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press, http://adsabs.harvard.edu/abs/2008gady.book.....B
  * Bonamigo et al. (2015) Bonamigo M., Despali G., Limousin M., Angulo R., Giocoli C., Soucail G., 2015, MNRAS, 449, 3171
  * Chandrasekhar (1987) Chandrasekhar S., 1987, Ellipsoidal Figures of Equilibrium. Dover
  * Chandrasekhar & Lebovitz (1962) Chandrasekhar S., Lebovitz N. R., 1962, ApJ, 136, 1037
  * Cole & Lacey (1996) Cole S., Lacey C., 1996, MNRAS, 281, 716
  * Corless & King (2007) Corless V. L., King L. J., 2007, MNRAS, 380, 149
  * Despali et al. (2014) Despali G., Giocoli C., Tormen G., 2014, MNRAS, 443, 3208
  * Despali et al. (2017) Despali G., Giocoli C., Bonamigo M., Limousin M., Tormen G., 2017, MNRAS, 466, 181
  * de Zeeuw & Franx (1991) de Zeeuw T., Franx M., 1991, ARA&A, 29, 239
  * Diemer & Joyce (2019) Diemer B., Joyce M., 2019, ApJ, 871, 168
  * Drakos et al. (2019) Drakos N. E., Taylor J. E., Berrouet A., Robotham A. S. G., Power C., 2019, MNRAS, 487, 993
  * Dubinski & Carlberg (1991) Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496
  * Eisenstein & Loeb (1995) Eisenstein D. J., Loeb A., 1995, ApJ, 439, 520
  * Eke et al. (1996) Eke V. R., Cole S., Frenk C. S., 1996, MNRAS, 282, 263
  * Erkal et al. (2016) Erkal D., Sanders J. L., Belokurov V., 2016, MNRAS, 461, 1590
  * Frenk et al. (1988) Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 327, 507
  * Gouin et al. (2021) Gouin C., Bonnaire T., Aghanim N., 2021, A&A, 651, A56
  * Hayashi et al. (2007) Hayashi E., Navarro J. F., Springel V., 2007, MNRAS, 377, 50
  * Hellwing et al. (2021) Hellwing W. A., Cautun M., van de Weygaert R., Jones B. T., 2021, Phys. Rev. D, 103, 063517
  * Hopkins et al. (2005) Hopkins P. F., Bahcall N. A., Bode P., 2005, ApJ, 618, 1
  * Jiang et al. (2015) Jiang L., Cole S., Sawala T., Frenk C. S., 2015, MNRAS, 448, 1674
  * Jing & Suto (2002) Jing Y. P., Suto Y., 2002, ApJ, 574, 538
  * Johnson et al. (2021) Johnson T., Benson A. J., Grin D., 2021, The Astrophysical Journal, 908, 33
  * Kazantzidis et al. (2006) Kazantzidis S., Zentner A. R., Kravtsov A. V., 2006, The Astrophysical Journal, 641, 647
  * Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
  * Kuhlen et al. (2007) Kuhlen M., Diemand J., Madau P., 2007, ApJ, 671, 1135
  * Lau et al. (2021) Lau E. T., Hearin A. P., Nagai D., Cappelluti N., 2021, MNRAS, 500, 1029
  * Lee et al. (2005) Lee J., Jing Y. P., Suto Y., 2005, ApJ, 632, 706
  * Lee et al. (2017) Lee C. T., Primack J. R., Behroozi P., Rodríguez-Puebla A., Hellinger D., Dekel A., 2017, MNRAS, 466, 3834
  * Limousin et al. (2013) Limousin M., Morandi A., Sereno M., Meneghetti M., Ettori S., Bartelmann M., Verdugo T., 2013, Space Sci. Rev., 177, 155
  * Macciò et al. (2007) Macciò A. V., Dutton A. A., van den Bosch F. C., Moore B., Potter D., Stadel J., 2007, MNRAS, 378, 55
  * Macciò et al. (2008) Macciò A. V., Dutton A. A., van den Bosch F. C., 2008, MNRAS, 391, 1940
  * Mansfield & Avestruz (2021) Mansfield P., Avestruz C., 2021, MNRAS, 500, 3309
  * Merritt & Quinlan (1998) Merritt D., Quinlan G. D., 1998, ApJ, 498, 625
  * Moore et al. (1999) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19
  * Moore et al. (2004) Moore B., Kazantzidis S., Diemand J., Stadel J., 2004, MNRAS, 354, 522
  * Morandi et al. (2011) Morandi A., Pedersen K., Limousin M., 2011, ApJ, 729, 37
  * Muñoz-Cuartas et al. (2011) Muñoz-Cuartas J. C., Macciò A. V., Gottlöber S., Dutton A. A., 2011, MNRAS, 411, 584
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Parkinson et al. (2008) Parkinson H., Cole S., Helly J., 2008, MNRAS, 383, 557
  * Peter et al. (2013) Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS, 430, 105
  * Riebe et al. (2013) Riebe K., et al., 2013, Astronomische Nachrichten, 334, 691
  * Rossi et al. (2011) Rossi G., Sheth R. K., Tormen G., 2011, MNRAS, 416, 248
  * Schneider et al. (2012) Schneider M. D., Frenk C. S., Cole S., 2012, J. Cosmology Astropart. Phys., 2012, 030
  * Smith & Watts (2005) Smith R. E., Watts P. I. R., 2005, MNRAS, 360, 203
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Tulin & Yu (2018) Tulin S., Yu H.-B., 2018, Phys. Rep., 730, 1
  * van Daalen et al. (2012) van Daalen M. P., Angulo R. E., White S. D. M., 2012, MNRAS, 424, 2954
  * van de Ven et al. (2003) van de Ven G., Hunter C., Verolme E. K., de Zeeuw P. T., 2003, MNRAS, 342, 1056
  * van Uitert et al. (2017) van Uitert E., et al., 2017, MNRAS, 467, 4131
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 482, 1525
  * Vasiliev & Athanassoula (2012) Vasiliev E., Athanassoula E., 2012, MNRAS, 419, 3268
  * Vega-Ferrero et al. (2017) Vega-Ferrero J., Yepes G., Gottlöber S., 2017, MNRAS, 467, 3226
  * Vera-Ciro et al. (2011) Vera-Ciro C. A., Sales L. V., Helmi A., Frenk C. S., Navarro J. F., Springel V., Vogelsberger M., White S. D. M., 2011, MNRAS, 416, 1377
  * Wang et al. (2020) Wang J., Bose S., Frenk C. S., Gao L., Jenkins A., Springel V., White S. D. M., 2020, Nature, 585, 39
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
  * Zentner et al. (2005) Zentner A. R., Kravtsov A. V., Gnedin O. Y., Klypin A. A., 2005, ApJ, 629, 219

##  Appendix A Test of Non-conservation

Figure 8: Probability distribution functions for the axis ratio $a_{3}$
assuming a conserved energy tensor. Results were derived by applying our model
to 10,000 realizations of merger trees with mass $10^{14}\mathrm{M}_{\odot}$
at $z=0$, and are compared to N-body simulation results from Bonamigo et al.
(2015).

A key assumption of this work is that the halo’s energy tensor is not
precisely conserved during its evolution. There is strong justification for
this from N-body simulations, as discussed in §2.1. However, it is interesting
to explore what would happen to our model if these phenomena were ignored. As
discussed in §2.4, non-conservation involves 3 tunable parameters,
$(b,\gamma,\alpha).$ Figure 8 explores how essential these features are to the
model by setting $b=\gamma=\alpha=0.$ As expected, the resulting probability
distribution function of $a_{3}$ now strongly disagrees with that measured
from N-body simulations. Without these non-conservation effects, especially
sphericalization, the axis ratios are biased too small, resulting in halos
that are far too elongated. Clearly these corrections are a necessary
component of our model—Figure 8 reinforces the result of §2.1 that the
underlying phenomena of halo mergers are being successfully captured in this
model.