The Energy Shear of Protohaloes: Marcello Musso, Giulia Despali, and Ravi K. Sheth

Download as pdf or txt
Download as pdf or txt
You are on page 1of 9

MNRAS 000, 1–9 (2024) Preprint 31 May 2024 Compiled using MNRAS LATEX style file v3.

The energy shear of protohaloes


1⋆ 2,3,4 5,6
Marcello Musso , Giulia Despali , and Ravi K. Sheth †
1 Departamento de Fı́sica Fundamental and IUFFyM, Universidad de Salamanca, E-37008 Salamanca, Spain
2 Dipartimento di Fisica e Astronomia, Alma Mater Studiorum Università di Bologna, via Gobetti 93/2, I-40129 Bologna, Italy
3 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, I-40129 Bologna, Italy
4 INFN-Sezione di Bologna, Viale Berti Pichat 6/2, I-40127 Bologna, Italy
5 Center for Particle Cosmology, University of Pennsylvania, 209 S. 33rd St., Philadelphia, PA 19104, USA
6 The Abdus Salam International Center for Theoretical Physics, Strada Costiera, 11, Trieste 34151, Italy
arXiv:2405.20207v1 [astro-ph.CO] 30 May 2024

Accepted XXX. Received YYY; in original form ZZZ

ABSTRACT
As it collapses to form a halo, the shape of a protohalo patch is deformed by the initial shear field. This deformation
is often modeled using the ‘deformation’ tensor, constructed from second derivatives of the gravitational potential,
whose trace gives the initial overdensity. However, especially for lower mass protohalos, this matrix is not always
positive definite: One of its eigenvalues has a different sign from the others. We show that the evolution of a patch is
better described by the ‘energy shear’ tensor, which is positive definite and plays a direct role in the evolution. This
positive-definiteness simplifies models of halo abundances, assembly and of the cosmic web.
Key words: large-scale structure of Universe

1 INTRODUCTION of fastest infall (van de Weygaert & Bertschinger 1996): i.e.,


in these models, an initial sphere becomes increasingly non-
Understanding how dark matter halos form and cluster – how
spherical if the deformation tensor is anistropic. However, if
their abundance and spatial distribution evolves – potentially
one uses the actual positions of the particles that make up
unlocks a wealth of information about cosmology and pro-
the patch to define the initial inertia tensor, then one finds
vides the scaffolding on which models of galaxy formation
that the initial shapes are not spherical: the longest axis of
are built. In addition, the intrinsic distribution of halo shapes
the inertia tensor is very well-aligned with the direction of
is a leading systematic for constraining cosmological models
fastest infall (Porciani et al. 2002a; Despali et al. 2013; Lud-
with the next generation of galaxy cluster surveys, so halo
low et al. 2014). This alignment has a transparent physical
formation has been the subject of significant work.
interpretation: it allows particles from further away to reach
Numerical simulations have shown that the protohalo
the halo center at about the same time as particles that were
patches, which evolve into virialized halos at later times, form
initially less distant, but were collapsing towards the center
from special regions in the initial fluctuation field. Following
more slowly (also see Musso & Sheth 2023). The residual
Bardeen et al. (1986), most models of halo formation asso-
misalignment of the two tensors is instead what gives rise
ciate protohalo patches with peaks in the smoothed overden-
to torques/angular momentum (White 1984; Porciani et al.
sity field, which then ‘collapse’ due to gravitational instabil-
2002b; Cadiou et al. 2021).
ity. The simplest models assumed the patch and its collapse
were both spherically symmetric (Gunn & Gott 1972), with The word ‘collapse’ suggests that the object shrinks, in co-
later work exploring ellipsoidal rather than spherical collapse moving coordinates, as it evolves. Stated more carefully, the
(Bond & Myers 1996; Sheth et al. 2001). principal axes of the inertia tensor are expected to shrink as
These models typically distinguish between the ‘shape’ and they are squeezed by the velocity flows that define the defor-
the ‘deformation’ or ‘tidal’ tensor of the initial patch. The lat- mation tensor. Squeezing suggests that the principal axes of
ter uses the anisotropy of the second spatial derivatives of the the deformation tensor all have the same sign. Indeed, some
potential to infer how the initial shape is deformed. However, have used the requirement that all three axes have the same
to approximate the initial shape, it is usual to use the tensor sign as the guiding principle from which to build a model of
of second spatial derivatives of the density field. Since the halo abundances (Lee & Shandarin 1998). However, simula-
two tensors are correlated, this proxy for the protohalo shape tions have shown that this is simply not true for a significant
predicts preferential alignment of the shortest axis (corre- fraction of lower mass protohalo patches (Despali et al. 2013).
sponding to the steepest descent direction) with the direction Some problems also arise when working with peaks in the
matter overdensity field – the trace of the deformation tensor
– to identify protohalo centers. There is no guarantee that the
⋆ E-mail: mmusso@usal.es density peak coincides with the position onto which the local
† E-mail: shethrk@upenn.edu gravitational flow converges, which is a natural choice to iden-

© 2024 The Authors


2 M. Musso, G. Despali and R. K. Sheth
tify the geometric center of the protohalo patch. Moreover, is the peculiar kinetic energy tensor, and
some of the integrals required to self-consistently define peak Z
3
statistics diverge. Both these disadvantages are overcome if uij ≡ d3 rρ(r, t)(ri − rcm,i )(∇j ϕ − [∇j ϕ]cm,j ) , (7)
I V
one works instead with peaks in the energy overdensity field
(Musso & Sheth 2021), rather then in the matter overdensity. is the potential energy overdensity tensor, or the ‘energy
The main goal of the present work is to show that the energy shear’. The cosmological constant does not appear explicitly
analogue of the deformation tensor, which we refer to as the in equation (5); it only affects the perturbations through its
‘energy shear’ tensor, also has an appealing property: all its effect on the background evolution of ρ̄ and H.
eigenvalues have the same sign. Parameterizing the initial value of the comoving inertia
We motivate why this sign requirement is expected in Sec- tensor as Iij /a2 = (M/5)(AAT )ij , equation (5) is solved by
tion 2, and show tests in simulations in Section 3. A final  T

section discusses our findings and conclusions. Iij M T L L tr(AA )
≃ (AA )ij − D(z)(uij + uji ) (8)
a2 5 3
at first order in D(z), the growth function of matter density
2 THE ENERGY OVERDENSITY TENSOR perturbations, since kij is automatically of second order. In
this expression, uL ij is the linearized version of uij from equa-
In this Section, we study the conditions required for the tri- tion (7) divided by D, and is therefore time independent. For
axial collapse of protohalo patches of arbitrary shape, be- ease of notation, in the following we will drop the superscript
fore discussing how they relate to analytical models typically L and simply assume that all quantities are evaluated at the
based on spherically averaged quantities. lowest order in density perturbations, and rescaled by D(z).
This result can also be obtained directly by writing equa-
tion (1) in Lagrangian coordinates as
2.1 Evolution of the inertia tensor of protohaloes Z
The formation of a dark matter halo can be described as the Iij = a5 ρ̄ d3 q (qi + Ψi ) (qj + Ψj ) , (9)
V
collapse of all three axes of its inertia tensor
Z where r − rcm = a(q + Ψ(q, t)) and Ψ is the displacement
Iij ≡ d3 rρ(r, t)(ri − rcm,i )(rj − rcm,j ) , (1) relative to the center of mass from the initial comoving co-
V
ordinate q. At first order in Lagrangian perturbation theory,
where V is a freely evolving volume (in the physical coordi- since Ψ ≃ −D(∇ϕ(q) − [∇ϕ]cm ), one gets equation (8).
nates r) containing the conserved mass M , and rcm denotes While non-linear dynamics may modify the subsequent
its center of mass. Mass conservation guarantees that time evolution, linear theory is sufficient to predict if and when
derivatives “go through” the integration sign and the time each axis decouples from the Hubble flow and starts recol-
dependence of ρ(r, t). Thus, the evolution of the inertia ten- lapsing. Since non-linear corrections, even large ones, mani-
sor of a region comoving with the fluid (the so-called virial fest themselves at small scales, they will not be able to reverse
equation) reads an ongoing collapse. Hence, the condition for triaxial collapse
to happen is that (the symmetric part of) uij be positive def-
I¨ij 1 Λ inite: i.e., all its eigenvalues have the same sign.
= 2Kij + (Uij + Uji ) + Iij (2)
2 2 3 Taking the trace of equation (8) returns
(see e.g. Chandrasekhar 1969; Binney & Tremaine 1987, for
5I tr(AAT ) h ϵi
the Λ = 0 case), where Kij and Uij are the kinetic and po- RI2 ≡ ≃ a2 1 − 2D(z) , (10)
tential energy tensors of the body. They are defined as 3M 3 3

1
Z with ϵ ≡ tr(u), which corresponds (at first order in pertur-
Kij ≡ d3 rρ(r, t)(ṙi − ṙcm,i )(ṙj − ṙcm,j ) , (3) bations) to the result that the collapse of the inertial radius
2 V
Z RI is well described by spherical collapse with overdensity ϵ.
Uij ≡ − d3 rρ(r, t)(ri − rcm,i )(∇j Φ − [∇j Φ]cm,j ) , (4) This is the basis for identifying protohalos with local maxima
V of the energy overdensity (Musso & Sheth 2021), since these
where ∇Φ is the gravitational attraction due to matter, so are minima of the collapse time of RI .
that r̈ = −∇Φ + (Λ/3)r, and [∇Φ]cm is the one of the center Naively, for a halo to form, one might imagine that all the
of mass. three axes of the inertia tensor Iij should collapse at the same
To describe the evolution of Iij with respect to the back- time as RI . Equation (8) would then require that
ground, we split the matter gravitational potential as ∇Φ =
u + uT
4πGρ̄(r/3 + ∇ϕ), where ϕ is the potential perturbation obey- AAT ≃ tr(AAT ) , (11)
ing ∇2 ϕ = δm , and δm = (ρ/ρ̄) − 1 is the matter density 2ϵ
perturbation. We then rewrite equation (2) as in which case one would have
 ··
Iij
 ·
Iij 4πGρ̄ I M 2 u + uT  ϵ u + uT
+ 2H = 4kij − (uij + uji ) 2 , (5) Iij ≃ a tr(AAT ) 1 − 2D ≃ I, (12)
a2 a2 3 a 5 2ϵ 3 2ϵ
predicting that the eigenvectors of Iij and of uij + uji are
where a is the scale factor, H = ȧ/a is the Hubble parameter,
to be aligned, and their eigenvalues proportional. This pic-
I ≡ Ikk is the trace of Iij ,
ture is not entirely correct: we know in fact that protohalo
 ·  ·
boundaries tend to follow equipotential surfaces (Musso &
Z
1 ri − rcm,i rj − rcm,j
kij ≡ d3 rρ(r, t) (6)
2 V a a Sheth 2023), which gives an independent prescription for Iij .

MNRAS 000, 1–9 (2024)


Protohalo shear 3
It must however be true to some degree, at least for the eigen- 2.3 Quantifying the amount of anisotropy
vectors, since equation (8) does show that the anisotropy of
The anisotropic strength of the protohalo environment can
uij favours infall from the direction of strongest compression,
be quantified by the magnitude of the traceless shear
with slower infall in the direction of weakest compression.
On the other hand, the mean deformation tensor (some- 3 3
qu2 ≡ tr(ūūT ) = (uij − ϵδij /3)(uji − ϵδji /3)
times also called shear, or tidal tensor when referring to its 2 2
traceless part) is defined as (λ1 − λ2 )2 (λ1 − λ3 )2 (λ2 − λ3 )2
= + + . (19)
1
Z 2 2 2
qij ≡ d3 r ρ(r, t) ∇i ∇j ϕ . (13)
M V (The suffix u is to distinguish it from the analogous variable
At leading order in perturbations, when ρ(r, t) ≃ ρ̄, qij that can be defined from the traceless part of the deforma-
is equal to the gradient of the center of mass acceleration tion tensor, which is often called q because it contributes a
∇i [∇j ϕ]cm , and its trace, the volume average of ∇2 ϕ, equals quadrupolar signature in perturbation theory.)
the mean matter overdensity. Note that although the trace of For a random Gaussian matrix, and therefore for random,
qij determines the evolution of the volume V , the tensor qij unconstrained positions in a Gaussian field, this random vari-
does not otherwise play a direct role in the evolution of Iij : ate would not be correlated with ϵ (the trace of uij ) and would
what matters is uij , rather than qij (cf. equation 8). follow a χ2 -distribution with five degrees of freedom (Sheth
& Tormen 2002). However, it is easy to see that simply re-
questing the energy shear to be positive definite induces a
2.2 Spherical initial volumes correlation between the two variables. In fact, since
When building models, one often approximates protohaloes ϵ2 − qu2 = 3(λ1 λ2 + λ1 λ3 + λ2 λ2 ) (20)
as spheres of the same mass. If V is a sphere of radius R
centered at x, containing mass M = (4π/3)ρ̄R3 , at leading and the right hand side is positive when all eigenvalues are
order in density perturbations its energy overdensity ϵR is positive, if uij is positive definite one has that ϵ > qu . Fur-
thermore, the limit can be saturated only if all three eigen-
d3 k
Z
ϵR (x) = δL (k) W2 (kR) eik·x (14) values are zero, which never happens in practice.
(2π)3 Some intuitive understanding of this correlation can be ob-
where δL (k) is the linear matter overdensity field normalized tained by noticing that if ϵ is small, then all eigenvalues are
to σ8 , and W2 (kR) ≡ 15j2 (kR)/(kR)2 (Musso & Sheth 2021). likely to be equally small (since they all have the same sign
Similarly, its energy shear is given by and can no longer cancel out from their sum). Hence, their
differences must be small as well. Below, we look for this
d3 k
Z
ki kj
uR,ij (x) = δL (k) 2 W2 (kR) eik·x , (15) correlation in the statistics of protohaloes.
(2π)3 k
whose trace is equation (14). Peaks in the energy overdensity
field are locations where ∇i ϵR = 0 and the Hessian 3 SIMULATION MEASUREMENTS
d3 k
Z
∇i ∇j ϵR = − δL (k) ki kj W2 (kR) eik·x (16) We now test our ansatz on the positivity of the eigenvalues
(2π)3
of the energy overdensity tensor in the protohaloes of two
is negative definite. simulations from the SBARBINE suite (Despali et al. 2016),
For a sphere, the mean mass overdensity δR , its deforma- namely Bice and Flora. Both simulations contain 10243 dark
tion tensor qR,ij (whose trace is δR ) and Hessian ∇i ∇j δR are matter particles in a periodic box of side Lbox = 125h−1 Mpc
obtained replacing W2 (kR) in equations (14), (16) and (15) and Lbox = 2h−1 Gpc, with initial conditions generated at z =
with the standard Top Hat filter W1 (kR) ≡ 3j1 (kR)/kR. 99. The corresponding particle masses are 1.55 × 108 h−1 M⊙
Hence, while the difference between quantities associated for Bice and 6.35 × 1011 h−1 M⊙ for Flora. The simulations
with energy rather than mass density is conceptually rather adopt a Planck13 background cosmology: Ωm = 0.307, ΩΛ =
important, in Fourier space the difference simply boils down 0.693, σ8 = 0.829 and h = 0.677.
to using a different, more UV-suppressed smoothing of the Haloes are identified at z = 0 using a Spherical Over-
overdensity field δ(k) than the usual Top-Hat one. density halo finder with a density threshold of 319 times
The distributions of these Gaussian variates are character- the background density, corresponding to the virial overden-
ized by the spectral moments sity. Our Flora sample contains 5378 haloes, distributed as
2
Z
dk k3 P (k) 2j 2 follows: all the 1378 halos more massive than 1015 h−1 M⊙ ,
σjn ≡ k Wn (kR) , (17) 2000 randomly selected haloes with masses between 1014 and
k 2π 2
1015 h−1 M⊙ , and 2000 randomly selected haloes with masses
where, following Musso & Sheth (2021),
between 4 × 1013 and 1014 h−1 M⊙ . Our Bice sample contains
Wn (x) ≡ (2n + 1)!!
jn (x)
. (18) 5000 randomly selected haloes with masses between 1011 and
xn 4×1013 h−1 M⊙ . For each halo identified at z = 0, we use ‘pro-
In particular, σ02 and σ22 are the standard deviations of ϵR tohalo’ to refer to the region occupied by that halo’s particles
and ∇2 ϵR , whereas σ01 and σ21 are those of δR and ∇2 δR . in the initial conditions.
For a homogeneous spherical density perturbation, and Since we are interested in initial quantities, we work with
only in this case, the two quantities are equal: ϵR = δR , and the first snapshot of the simulation. We measured each pro-
the traceless parts of uij and qij are both zero. This is be- tohalo’s center of mass position and velocity, qcm and vcm , by
cause only the k = 0 mode of δ(k) matters in this case, but averaging over all its particles. We then constructed estima-
both filters tend to 1 as k → 0, so they make no difference. tors for its inertia and potential energy overdensity tensor,

MNRAS 000, 1–9 (2024)


4 M. Musso, G. Despali and R. K. Sheth

Figure 1. Left panel. Eigenvalues of the potential energy overdensity tensor uij . Only 5% of protohaloes have one negative eigenvalue; this
is our main result. Right panel. Eigenvalues of the mean deformation tensor. The lowest eigenvalue is negative in about 40% of the all
protohaloes, and in more than 50% of the lower mass ones, in stark contrast to the left panel.

7 1, Flora 1, Flora
2, Flora 2, Flora
6 3, Flora 2.0 3, Flora
1, Bice
5 1, Bice
2, Bice
3, Bice 1.5 2, Bice
4
P( i/ 02)
P( i/ 02)

3, Bice

3 1.0
2
0.5
1
0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0
i/ 02 i/ 02

Figure 2. Left panel. Distribution of the (normalized) eigenvalues λi of the potential energy overdensity tensor. Filled and open histograms
show results for all sampled protohaloes having masses larger than 4 × 1014 h−1 M⊙ (from Flora) and larger than 1011 h−1 M⊙ (from Bice),
respectively. Right panel. Same, but now for the eigenvalues λ̄i of the traceless part of the potential energy overdensity tensor. Filled and
open histograms are very similar, indicating that differences in the left panet are mainly driven by the mass-dependence of the trace.

defined in equations (1) and (7): initial displacement field Ψ, and estimate the mean deforma-
NH
tion tensor (13) as:
X
Iˆij ≡ (q (n) − qcm )i (q (n) − qcm )j (21) NG
1 X ∂Ψj
n=1 q̂ij ≡ , (23)
P NH NG ∂qi k
(n) (n) k=1
n=1 (q − qcm )i (v − vcm )j /f DH
ûij ≡ −3 PNH (n) , (22) where NG is the number of grid cells contained within the La-
n=1 |q − qcm |2
grangian protohalo volume. This means that we also include
where D is the ΛCDM density perturbation growth factor (at empty grid cells (i.e. those not occupied by halo particles)
the redshift z of the snapshot), f = d ln D/d ln a, and n runs that are located inside the Lagrangian volume, because they
over the NH protohalo particles. Note that Iij is not divided too contribute to the total potential field that acts on the
by NH . To obtain accelerations from velocities we used the protohalo. This slightly refines the measurement one obtains
Zel’dovich approximation, in which ∇ϕ ≃ v/Ḋ = v/f DH. if one uses the particle grid points only (see Despali et al.
For each protohalo we also measure the mean deforma- 2013, for more details).
tion tensor, defined in equation (13) as the average over the The left panel of Fig. 1 shows the eigenvalues λ1 , λ2 and λ3
protohalo particles of the second spatial gradient of the grav- of ûij + ûji as a function of halo mass. In the range of masses
itational potential at each position. Since we need to take one that we sampled, only one out of the 5378 Flora protohaloes
more derivative, instead of getting the accelerations from the had one slightly negative eigenvalue. For Bice protohaloes,
particle velocities (evaluated at irregular positions), we com- which extend to much lower masses, this fraction rises to 0.5%
pute the gradient of the initial displacements that were im- including only haloes mor massive than 1012 h−1 M⊙ , and to
posed to create the initial conditions with N-GenIC (Springel nearly 5% for all haloes in the sample. Nevertheless, this
et al. 2005), which are computed on the grid. In practice, we is substantially smaller than the fractions shown in Despali
select the grid points occupied by each protohalo, take their et al. (2013) for density-based statistics – i.e. the eigenvalues

MNRAS 000, 1–9 (2024)


Protohalo shear 5

15 15
10 10
5 5
0 0
15 Flora 15 Flora
P[cos( ij)]

P[cos( ij)]
Bice 10 Bice
10
5 5
0 0
15 15
10 10
5 5
0 0
0 1 0 1 0 1 0 1 0 1 0 1
cos( ij) cos( ij)
Figure 3. Cosine of the angle θij between the i-th eigenvector of the inertia tensor and the j-th one of the energy shear (left) or of the
deformation tensor (right). The alignment with the energy shear is stronger.

of q̂ij – which we also show in the right panel of Fig. 1. This


fraction amounts to about 40 % of all protohaloes, and nearly 30
half of those between 1011 h−1 M⊙ and 1012 h−1 M⊙ . Even at Flora, deformation tensor
Flora, energy tensor
masses close to 1014 h−1 M⊙ one of the eigenvalues may still 25 Bice, deformation tensor
be negative. This obvious difference, that the energy shear Bice, energy tensor
tensor is positive-definite even when the deformation tensor 20
P(matrix cosine)

is not, is the main result of this paper.


The left panel of Fig. 2 shows the distributions of the three 15
eigenvalues of ûij + ûji , scaled by σ02 defined in equation
(17), in the two simulations: filled histograms are for Flora
10
protohaloes more massive than 4 × 1014 h−1 M⊙ , empty ones
5
are for Bice protohaloes more massive than 1011 h−1 M⊙ . The
less massive objects clearly have smaller λi /σ02P
. This is not so 0
surprising, since we know that the trace, ϵ = i λi , depends 0.70 0.75 0.80 0.85 0.90 0.95 1.00
matrix cosine
on mass (Musso & Sheth 2021). More precisely, whereas the
unnormalized ϵ increases at smaller masses, ϵ/σ02 tends to de-
crease, which explains why all the λi /σ02 tend to be smaller. Figure 4. Distribution of the matrix cosine (equation 25) of the
To remove this dependence, the right panel of Fig. 2 shows – inertia tensor Iij with the energy shear tensor uij (orange) or
for the same objects – the distributions of the eigenvalues of the deformation tensor qij (blue). Filled histograms are for Flora
the traceless part of the potential energy overdensity tensor, (massive halos) and empty ones for Bice (smaller masses).

λ̄i ≡ λi − ϵ/3, (24)


inertia tensor with the energy shear, for each protohalo we
again normalized by σ02 (same color coding). Now the his- computed the matrix cosine
tograms for the two simulations (filled and empty) are very q
similar, showing that these distributions are nearly mass in- ˆ ≡ ûij Iˆji / ûkl ûlk Iˆmn Iˆmn ,
cos(û, I) (25)
dependent.
A detailed analysis of the mutual alignment of the three which quantifies with a single number the alignment between
eigenvectors of the inertia tensor with those of the energy the eigenvectors of the two matrices and the similarity of their
shear is shown in the panels on the left of Fig. 3. This align- eigenvalues. We repeated the same measurement for the de-
ment is quite strong, and stronger than the one with the formation tensor. A value of the cosine closer to 1 indicates
deformation tensor, shown for comparison in the panels on stronger alignment, and is what equation (12) would predict
the right. Interestingly, in both cases the alignment is better exactly. The orange histograms in Fig. 4 show that most mea-
in Bice (i.e. for smaller halo masses) than in Flora (larger sured cosines are greater than 0.9, although the alignment is
masses). slightly better for massive halos (filled histogram). For com-
To get a more concise description of the overlap of the parison, the blue histograms show the result of replacing û

MNRAS 000, 1–9 (2024)


6 M. Musso, G. Despali and R. K. Sheth

3.5 0.8
R1, Flora Flora
3.0 R2, Flora 0.7 Bice
R3, Flora 2
5
2.5 R1, Bice 0.6
R2, Bice
R3, Bice 0.5
2.0

2)
P(qu2/ 02
P(Ri)

0.4
1.5
0.3
1.0
0.2
0.5 0.1
0.0 0.0
0.5 0.0 0.5 1.0 1.5 2.0 0 2 4 6 8 10
ratio of eigenvalues Ri qu2/ 2
02

Figure 5. Ratios of the i-th eigenvalue λi of the energy shear uij 2 /σ 2 for Flora (filled) and Bice (open)
Figure 6. Distribution of qu 02
to the i-th eigenvalue a2i of the inertia tensor Iij , both normalized protohaloes; smooth curve shows the expected distribution for un-
to their traces: Ri ≡ (λi /ϵ)/(a2i /I). Equation (12) would predict constrained positions (a χ2 distribution with 5 degrees of freedom).
this ratio to be 1, whereas on average R1 > 1 and R3 < 1, mean-
ing that protohaloes are slightly less elongated than this simplest
prediction.

with the deformation tensor q̂: clearly, the inertia tensor is


better aligned with the energy shear than with the deforma-
tion tensor.
The direct comparison of the eigenvalues a2i of the inertia
tensor and λi of the energy shear, each normalized to their
sum I and ϵ, is shown in Fig. 5, which displays their ratios
Ri ≡ (λi /ϵ)(I/a2i ). Equation (12) would predict that the two
tensors are proportional to each other, and therefore each
ratio should be equal to 1. This prediction is clearly not ver-
ified. Fig. (5) suggests instead that protohaloes are slightly
less elongated than they would have been if their inertia ten-
sor and energy shear were proportional.
The histograms in Fig. 6 show the distribution measured
in protohalos (in Flora and Bice respectively) of the shear
strength qu2 /σ022
defined in equation (19). This distribution is Figure 7. Correlation between the energy overdensity ϵ and the
broader than a χ25 , shown for comparison by the smooth solid traceless energy shear u of protohaloes. The color coding displays
curve, which would be expected if the eigenvalue of the energy the (logarithm of the) halo mass. Dashed line shows the limiting
overdensity tensor were free to take up any sign. Nevertheless, value (ϵ = qu ).
the it is still rather mass independent.
A scatter plot of ϵ vs qu (colored by protohalo mass) is dis-
played in Fig. 7, and shows a strong degree of (approximately homogeneous, in which case the two tensors are identical.
linear) correlation between the two variables. This behaviour Under this assumption, there is an exact analytical expres-
is also a consequence of the eigenvalues of uij being all posi- sion for the potential of the ellipsoid. Working with uij is on
tive, as argued in Section 2.3. the other hand completely general – albeit based on a per-
turbative expansion – and is therefore the correct choice for
describing the collapse of the inertia tensor Iij . The downside
is that uij has no closed analytical expression in terms of Iij
4 COMPARISON WITH PREVIOUS WORK
(to be expected, since it also depends on the density profile,
Our analysis bears some qualitative similarity with the results except for particular, highly symmetrical cases).
of Bond & Myers (1996), but with two important differences. Porciani et al. (2002a), and also later Ludlow et al. (2014),
First, they considered an initially spherical region, in which dropped the spherical protohalo assumption, but continued
case AAT in equation (8) would be diagonal and reduce to to assume that the deformation tensor provides the natural
R2 δij , where R is the Lagrangian radius of the sphere. While framework for understanding evolution. In particular, they
this approximation can be useful and significantly simplifies assumed that protohaloes are ellipsoidal with the longest axis
the calculations, protohalos are not spherical (or even ellip- being exactly aligned with the direction of maximum com-
soidal), but have rather irregular shapes. pression of the deformation tensor. While this assumption
Second, they focused on the deformation tensor qij , rather works well in practice, our analysis here suggests that it has
than on uij . This is correct only under the assumption (which no fundamental physical justification. Rather, its accuracy is
they make) that V is spherical (or at least ellipsoidal) and more of statistical origin (and is not as good as the one with

MNRAS 000, 1–9 (2024)


Protohalo shear 7
the energy shear). Furthermore, an exact alignment would We have thus added another ingredient to the recipe for
not generate angular momentum via gravitational coupling, characterizing protohaloes in terms of energy-related quanti-
since the torque induced by the external gravitational field ties:
on the protohalo is zero.
• the center of mass of a protohalo is much better identified
On the other hand, as shown by Musso & Sheth (2023),
by a local maximum of the energy overdensity field than of
protohaloes are well approximated by equipotential surfaces
the matter overdensity (Musso & Sheth 2021);
that deviate from spherical symmetry in response to the lo-
• protohalo shapes and boundaries can be described by
cal anisotropy of the infall potential. This gives an indepen-
equipotential surfaces, which further maximize the enclosed
dent prescription for determining the initial protohalo shape,
energy with respect to spherical configurations (Musso &
which then evolves under the action of uij . The symmetric
Sheth 2023);
part of uij determines the collapse and deformation, and the
• the energy overdensity tensor (whose eigenvalues must
antisymmetric one the torque (Musso & Sheth 2023), which
be all positive) does a much better job at predicting both the
can then generate angular momentum as described by Cadiou
collapse of the three axes and their initial orientation than
et al. (2021).
the deformation tensor (this work).
Despali et al. (2013) studied the eigenvalues of the de-
formation tensor, and noted that at lower masses (below More specifically, we have shown that the principal axes of
1013 h−1 M⊙ ) it is common for one of them to be negative. the inertia tensor of protohaloes are very well aligned with
In hindsight, this is not surprising since in a generic setup the eigenvectors of the energy tensor; the latter identify the
the average of ∇i ∇j ϕ (that is, qij ) describes the gradient of main directions of compression of the protohalo region. The
the center of mass acceleration (see discussion at the end of longest/shortest protohalo axis tends to align with the direc-
Section 2.1), which is not directly implicated in the collapse tion of maximum/minimum compression, corresponding to
dynamics. Rather, as we argue here, triaxial contraction fol- the largest/smallest eigenvalue of the energy tensor, for which
lows from the positivity of the eigenvalues of the energy shear infall is favoured/disfavoured. Moreover, no collapse should
uij , which our measurements in Fig. 1 confirm. be possible if one eigenvalue is negative. Our measurements
Working with uij has also practical advantages. For a strongly support this prediction: less than 5% of protohaloes
ΛCDM power spectrum the variance of ∇i ∇j ϵR remains fi- down to 1011 h−1 M⊙ have one negative eigenvalue (and even
nite, whereas that of ∇i ∇j δR diverges. For ϵR it is thus then, only slightly negative); none has two negative values.
possible to carry out the full calculation of the peak num- This is important for several reasons. First of all, it gives
ber density from first principles, without having to resort to an effective and unambiguous description of the protohalo
hand-waving changes of the smoothing filter like in a δR - environment: protohaloes can be identified with regions with
based approach (Musso & Sheth 2021). However, one might positive definite energy overdensity tensor, as this is what
also be interested in other statistics, like the number density triggers triaxial collapse. Secondly, while the trace of the ten-
of points with null determinant. These are the so-called criti- sor (the energy overdensity ϵ) sets the time of collapse and
cal events (Hanami 2001), where an eigenvalue changes sign, therefore determines the current halo mass, its traceless part
because – for instance – a peak and a saddle point merge and (the traceless energy shear) provides a set of variables that
disappear as the smoothing scale R changes. These events can be used to describe assembly bias. It can incorporate ef-
are Lagrangian proxies for halo mergers (Cadiou et al. 2020). fects of the anisotropy of the environment, (as done e.g. by
Computing their number density would require taking one Sheth & Tormen 2002; Musso et al. 2018, using the tidal de-
extra derivative, but the variance of the third derivative di- formation tensor, the analogous quantity for the density), and
verges, this time also for ϵR . On the other hand, if the eigen- would appear naturally in symmetry-based bias expansions
value changing sign is that of uij , the calculation remains (Sheth et al. 2013; Lazeyras et al. 2016; Desjacques et al.
well behaved. Defining critical events in terms of the energy 2018).
tensor, rather than the Hessian, would thus not only come In fact, in our formulation, the energy shear is expected to
with a stronger connection with dynamics, but also better source the first correction to the collapse time of RI , and
mathematical properties. hence to appear in any analytical expression of a critical
threshold for ϵ. Fig. 6 shows that, simply by working with
qu2 /σ02
2
, the mass-dependence can be nearly scaled-out; this
should vastly simplify the development of models of the ef-
5 CONCLUSIONS
fects of anisotropies on halo collapse. Furthermore, our anal-
Our measurements strongly support the conclusion that pro- ysis also shows (see Fig. 7) that the positivity of the eigenval-
tohalo patches can be identified with regions where the en- ues introduces a strong correlation between ϵ and the trace-
ergy overdensity tensor is positive definite. This fact is not less shear (which would be independent if uij were a generic
a simple statistical correlation, but a dynamical feature that Gaussian matrix, or equivalently, at random locations in the
follows directly from the evolution equations of the inertia Gaussian field). This can help to explain the scatter and the
tensor, which must collapse along three axes. On the other trend with mass of the measured values of ϵ in protohaloes
hand, the failure of the deformation tensor defined by the (Musso & Sheth 2021), thus enabling the construction of a
mass overdensity to predict triaxial collapse also has an ex- physically motivated model of the collapse threshold in the
planation. While this tensor is directly related to the change energy-based approach. This would provide a natural way for
of shape of a microscopic volume element, for a macroscopic the approach to include assembly bias effects.
body its physical interpretation is rather that of the gradient There are also a number of practical implications for ana-
of the center of mass acceleration: therefore, not a quantity lytical models of halo aboundance and of the cosmic web. For
directly related to the collapse dynamics. instance, halos are typically identified with peaks in the (mat-

MNRAS 000, 1–9 (2024)


8 M. Musso, G. Despali and R. K. Sheth
ter or) energy overdensity field. To distinguish between peaks as Figure 1 shows, the signatures of the potential energy and
and other stationary points, one requires that the Hessian of deformation tensors can be different, so the actual cosmic web
the field be positive definite (Bardeen et al. 1986). However, classifications will differ. Another obvious connection would
one can replace this constraint on the Hessian with the same be with the various works that have tried to predict the emer-
requirement on the energy overdensity tensor. Since the two gence of the cosmic web from the formation of caustics in
tensors are correlated, the Hessian would then acquire a mean the velocity field (Arnold et al. 1982; Hidding et al. 2014;
value proportional to −γuij , where γ denotes the normalized Neyrinck 2012; Feldbrugge et al. 2018). Unlike these works,
correlation coefficient between −∇2 ϵ and ϵ. Thus, it is likely which investigate the microscopic properties of the velocity
that if uij is positive definite, then so is the negative Hessian, field, the approach based on the energy shear attempts to
and the point is a peak. predict the evolution of macroscopic volumes. It is therefore
For a more quantitative grasp of the effect, note that re- more suitable to make analytical estimates of the statistics of
quiring positive definiteness shifts the mean of ϵ from zero objects.
to ⟨ϵ|+⟩ ≈ 1.65σ02 (see, e.g., Lee & Shandarin 1998). The Finally, as a practical matter, the measurement of uij in
trace of the Hessian x is commonly used to denote the peak simulations is very easy, since it only depends on positions
‘curvature’. If γ denotes the normalized correlation coefficient and velocities of protohalo particles. In contrast, qij requires
between x and ϵ, requiring positive definiteness of ϵ shifts the computing the gradient of the displacement, which needs to
mean of x to ⟨x|+⟩/σ22 ≈ γ ⟨ϵ|+⟩/σ02 ≈ 1.65γ; i.e., the trace be done at grid points and thus requires extra work.
of the Hessian is likely to be positive, so the Hessian itself is
likely to be positive definite, without explicitly requiring it
to be so.
ACKNOWLEDGEMENTS
Besides the Hessian, which requires two derivatives of the
field, some models of the cosmic web make use of third deriva- Thanks to Corentin Cadiou, Christophe Pichon and Dmitri
tives to identify the disappearance of peaks (or other station- Pogosyan for helpful discussions, and to the organisers and
ary points) via mergers (Hanami 2001; Cadiou et al. 2020). participants of the 2023 KITP Cosmic Web workshop. Our
In a ΛCDM scenario, the energy-based approach is attractive visit to KITP was supported by the National Science Foun-
because the statistics of its second derivatives depend on in- dation under PHY-1748958. GD acknowledges the funding by
tegrals that do not diverge. However, the statistics of its third the European Union - NextGenerationEU, in the framework
derivatives still diverge, so, to estimate abundances of merg- of the HPC project – “National Centre for HPC, Big Data and
ers will require using a different (essentially ad hoc) filter. Our Quantum Computing” (PNRR - M4C2 - I1.4 - CN00000013
results suggest that more principled progress can be made by – CUP J33C22001170001).
replacing the Hessian of the energy overdensity field with the
energy tensor. This would enable estimation of the relevant
statistics (which now have the same degree of convergence
DATA AVAILABILITY
in k-space as a gradient) without mathematical ambiguities.
I.e., a constraint based on actual collapse dynamics may be The simulation data and post-processed quantities used in
better suited for this task than a merely topological one, and, this work can be shared on reasonable request to the authors.
like for energy peaks, considering a better-motivated physical
quantity also solves the mathematical issues in the calcula-
tion of the statistics.
The idea of replacing the role of the Hessian with the en- REFERENCES
ergy overdensity tensor can be naturally extended to other Arnold V. I., Shandarin S. F., Zeldovich I. B., 1982, Geophysical
constituents of the cosmic web (Bond et al. 1996). The signa- and Astrophysical Fluid Dynamics, 20, 111
ture of the eigenvalues of the energy tensor could be used as Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ,
the fundamental marker that distinguishes between different 304, 15
cosmic web environments: the four possible different signa- Binney J., Tremaine S., 1987, Galactic dynamics. Princeton Uni-
tures of the potential energy tensor should be used to char- versity Press, NJ, USA
Bond J. R., Myers S. T., 1996, ApJS, 103, 1
acterize nodes, filaments, walls and voids. In all these cases
Bond J. R., Kofman L., Pogosyan D., 1996, Nature, 380, 603
the underlying physical interpretation of the energy tensor, Cadiou C., Pichon C., Codis S., Musso M., Pogosyan D., Dubois
based on predicting the collapse vs expansion of the various Y., Cardoso J. F., Prunet S., 2020, Mon. Not. Roy. Astron.
axes, would be the same. It would be interesting to inves- Soc., 496, 4787
tigate how well this approach could reproduce properties of Cadiou C., Pontzen A., Peiris H. V., Lucie-Smith L., 2021, Mon.
the cosmic web that are usually studies with density-based Not. Roy. Astron. Soc., 508, 1189
methods, such as its skeleton (Pogosyan et al. 2009) or con- Chandrasekhar S., 1969, Ellipsoidal figures of equilibrium. Yale
nectivity (Codis et al. 2018). University Press, New Haven, CT, USA
A classification of the constituents of the cosmic web based Codis S., Pogosyan D., Pichon C., 2018, Mon. Not. Roy. Astron.
Soc., 479, 973
on the dynamical role of the energy shear would connect with
Desjacques V., Jeong D., Schmidt F., 2018, Phys. Rep., 733, 1
the well-established literature on using the signature of the
Despali G., Tormen G., Sheth R. K., 2013, Mon. Not. Roy. Astron.
deformation tensor (i.e., the Hessian of the potential rather Soc., 431, 1143
than of the density) to characterize the cosmic web (e.g. Hahn Despali G., Giocoli C., Angulo R. E., Tormen G., Sheth R. K.,
et al. 2007; Feldbrugge et al. 2023, but also, e.g., Shen et al. Baso G., Moscardini L., 2016, MNRAS, 456, 2486
2006). Our work suggests that using the potential energy ten- Feldbrugge J., van de Weygaert R., Hidding J., Feldbrugge J., 2018,
sor instead should be more robust. Perhaps more importantly, JCAP, 05, 027

MNRAS 000, 1–9 (2024)


Protohalo shear 9
Feldbrugge J., Yan Y., van de Weygaert R., 2023, Mon. Not. Roy.
Astron. Soc., 526, 5031
Gunn J. E., Gott J. Richard I., 1972, ApJ, 176, 1
Hahn O., Porciani C., Carollo C. M., Dekel A., 2007, MNRAS,
375, 489
Hanami H., 2001, Mon. Not. Roy. Astron. Soc., 327, 721
Hidding J., Shandarin S. F., van de Weygaert R., 2014, Mon. Not.
Roy. Astron. Soc., 437, 3442
Lazeyras T., Musso M., Desjacques V., 2016, Phys. Rev., D93,
063007
Lee J., Shandarin S. F., 1998, ApJ, 500, 14
Ludlow A. D., Borzyszkowski M., Porciani C., 2014, MNRAS, 445,
4110
Musso M., Sheth R. K., 2021, Mon. Not. Roy. Astron. Soc., 508,
3634
Musso M., Sheth R. K., 2023, Mon. Not. Roy. Astron. Soc., 523,
L4
Musso M., Cadiou C., Pichon C., Codis S., Kraljic K., Dubois Y.,
2018, MNRAS, 476, 4877
Neyrinck M. C., 2012, Mon. Not. Roy. Astron. Soc., 427, 494
Pogosyan D., Pichon C., Gay C., Prunet S., Cardoso J. F., Sousbie
T., Colombi S., 2009, Mon. Not. Roy. Astron. Soc., 396, 635
Porciani C., Dekel A., Hoffman Y., 2002a, Mon. Not. Roy. Astron.
Soc., 332, 339
Porciani C., Dekel A., Hoffman Y., 2002b, MNRAS, 332, 325
Shen J., Abel T., Mo H., Sheth R. K., 2006, ApJ, 645, 783
Sheth R. K., Tormen G., 2002, MNRAS, 329, 61
Sheth R. K., Mo H., Tormen G., 2001, MNRAS, 323, 1
Sheth R. K., Chan K. C., Scoccimarro R., 2013, Phys. Rev. D, 87,
083002
Springel V., et al., 2005, Nature, 435, 629
White S. D. M., 1984, ApJ, 286, 38
van de Weygaert R., Bertschinger E., 1996, MNRAS, 281, 84

This paper has been typeset from a TEX/LATEX file prepared by


the author.

MNRAS 000, 1–9 (2024)

You might also like