© The Authors 2025

1 Introduction

Cosmic shear (CS) is arguably the most direct way to probe the fluctuations of the matter field on cosmological scales (Mandelbaum 2018). The slight deformations to the baseline shapes of galaxies are not affected by the clustering properties of those galaxies, making these measurements independent of galaxy bias. As a result, the information on the amplitude of matter fluctuations, σ₈, that can be extracted does not have the typical degeneracy with linear galaxy bias. However, an outstanding issue is that the shapes of galaxies before shearing are not entirely random but may be correlated with the large-scale gravitational field. Typically known as galaxy intrinsic alignments (IAs), this effect is highly dependent on galaxy-formation physics and contaminates CS measurements (Joachimi et al. 2015; Kiessling et al. 2015; Kirk et al. 2015; Troxel & Ishak 2015). Therefore, extracting reliable information from CS requires accurate modeling of IAs (Secco et al. 2022; Aricò et al. 2023; Preston et al. 2024).

Intrinsic alignment modeling is often based on a perturbative expansion (Desjacques et al. 2018). Early works modeled galaxy shapes as being linearly correlated to the large-scale tidal field, with the linear-alignment parameter (c_K) left free to describe any galaxy population (Catelan et al. 2001; Hirata & Seljak 2004). In recent years, this field has seen many advancements, with the formulation of models that include higher orders in perturbation theory (Blazek et al. 2015), complete formulations using the effective field theory of large-scale structures both in Eulerian Vlah et al. 2020, 2021; Bakx et al. 2023) and Lagrangian (Chen & Kokron 2024) space, and hybrid prescriptions that use a bias expansion combined with the results of N-body simulations to produce models that extend accurately into the nonlinear regime (Maion et al. 2024). An integral part of all of these models is the “shape bias” parameters – free coefficients that encode the average response of a galaxy’s shape to the large-scale properties of the matter field.

Although the bias parameters of a certain population of objects cannot be easily predicted a priori, measuring their values can shed light on the formation mechanism of said objects. As an example, the detection of nonzero halo assembly bias for density bias parameters (Gao et al. 2005; Wechsler et al. 2006; Gao & White 2007; Angulo et al. 2008) has shown that dependences on halos’ accretion histories have to be included in predictions from excursion-set theory (Bardeen et al. 1986; Bond et al. 1991; Lacey & Cole 1993) to be in agreement with halo-bias measurements from N-body simulations (Dalal et al. 2008). Furthermore, it explained why simple halo occupation distribution models, which assume the galaxy content of a halo depends solely on its mass (Zheng et al. 2005), are inadequate in describing realistic measurements from simulations and galaxy surveys (Croton et al. 2007; Zehavi et al. 2018; Contreras et al. 2024). Measurements of shape biases are still in their infancy but should contribute similarly to the understanding of shape formation as they advance into a more mature stage.

A significant amount of research has been carried out to analyze the IAs of halos and galaxies, in both gravity-only and hydrodynamical simulations (Tenneti et al. 2014, 2015b,a, 2016, 2020; Chisari et al. 2015; Velliscig et al. 2015; Hilbert et al. 2017; Delgado et al. 2023), but with little focus on large-scale shape-bias parameters. The first direct measurements of shape bias parameters of simulated halos were performed by Stücker et al. (2021) and Akitsu et al. (2021) employing the “anisotropic separate universe” technique (Schmidt et al. 2018; Wagner et al. 2014). In a later work, Akitsu et al. (2023) employed the quadratic bias technique (Schmittfull et al. 2015) to measure second-order shape biases. These works have shown that there is a power-law scaling of the linear shape bias with halo mass and demonstrated its dependence on concentration for high-mass halos. A universal relation between the linear shape bias and b₁, the linear density bias, was also established. Furthermore, they have quantified the dependence of second-order shape-bias parameters with halo mass and studied their dependence on b₁. Nevertheless, many interesting questions remain unanswered.

The precise origin of halo-shape alignments with the large-scale gravitational field is still poorly understood. In one possible scenario, the large-scale gravitational forces progressively deform halos after their formation, producing an alignment that grows over time. Another possibility is that the shape is already encoded during the formation of halos – for example through the anisotropy of infall velocities as determined by the large-scale gravitational potential. It should be possible to distinguish between these scenarios through measurements of the bias parameters. Beyond this, it is important to understand how alignments depend quantitatively on the (often arbitrary) boundary definition of the considered system and whether the alignment strength is affected by secondary properties beyond halo mass – such as the spin or formation time.

While studies of halos in gravity-only simulations may provide a qualitative understanding of necessary model ingredients, the ultimate goal of theoretical alignment studies is an understanding of the shape biases of galaxies. Determining these values in realistic scenarios, such as high-resolution hydrodynamical simulations, would be a great advancement toward obtaining informative priors that could be applied to IA models in weak-lensing analyses (Zennaro et al. 2022; Ivanov et al. 2025; Ivanov et al. 2024; Shiferaw et al. 2024; Zhang et al. 2024). Additionally, the comparison of measurements of galaxy shape bias parameters obtained from hydrodynamical simulations to those measured from observations could be used as a way to constrain “subgrid” models of galaxy formation. To answer all these questions, we need a fast and precise estimator of shape biases.

In this manuscript, we generalize the theoretical formalism presented in Stücker et al. (2025) to develop estimators of shape biases. The estimators are based on measurements of correlations in the joint probability distribution of Lagrangian large-scale tides and the shape of a given object. Similar to the separate universe (SU) approach, this yields simple and robust measurements of these biases but with the additional benefit of low computational costs, making the method viable for use in large hydrodynamical simulations. We validate the estimators through convergence tests and a comparison with literature results. Finally, we provide detailed measurements of the shape biases of halos and discuss the theoretical implications for the modeling of IAs.

The article is structured as follows. In Sect. 2, we derive the new estimators based on considerations from probability theory. In Sect. 3, we present the set of simulations we used to validate our results. In Sect. 4, we test the robustness of our shape bias measurements and provide accurate fitting formulae. Finally, in Sect. 5 we summarize our findings and discuss future applications.

2 Probabilistic bias

Considering an object composed of a finite set of particles, N, one can define its shape tensor, 𝓘. This quantity is a rank-2 tensor, generally defined as $${𝓘}_{ij}=\frac{1}{N}\sum _{n=1}^N(x_i^{(n)}-{\overline{x}}_i)(x_j^{(n)}-{\overline{x}}_j),$$(1)

with the corresponding reduced shape tensor $${\widetilde{𝓘}}_{ij}=\frac{1}{N}\sum _{n=1}^N\frac{(x_i^{(n)}-{\overline{x}}_i)(x_j^{(n)}-{\overline{x}}_j)}{r_n^2},$$(2)

where $x_i^{(n)}$ represents the Eulerian position of each particle, n is an index running over the particles in the set, and i runs over the three Cartesian coordinates, ${\overline{x}}_i$ is the mean of $x_i^{(n)}$, and $r_n^2\equiv {\sum }_i{(x_i^{(n)}-{\overline{x}}_i)}^2$. We used the non-reduced shape definition for all results in this article, except when comparing different definitions.

Now, we can define the normalized and trace-free part of ℐ: $$I_{ij}(\mathbf{x})=\frac{{\mathcal{T}}_{ij}(\mathbf{x})}{{\mathcal{I}}_0}-\frac{1}{3}{\delta }_{ij}^K,$$(3)

in which ℐ₀ is the trace of ℐ_ij, and ${\delta }_{ij}^K$ is the Kronecker delta symbol.

The tensor I describes the shape of an object, factoring out its physical size. Being a symmetric tensor, it can be diagonalized to give three eigenvalues, I_1,2,3, such that ${\sum }_{i=1}^3{I}_i=0$. Each of these measures the relative extension of the object in each direction of an orthogonal basis formed by its eigenvectors. The goal of this work is to understand the response of this shape to changes in the large-scale gravitational field.

Attempting to determine which large-scale properties could generate an effect, we find that the simplest object we can build, allowed by the equivalence principle and symmetry of the problem, is the Lagrangian tidal field: $$T_{ij}(\mathbf{q})={\mathrm{\partial }}_i{\mathrm{\partial }}_j{\mathrm{\Phi }}_G(\mathbf{q}),$$(4)

where Φ_G is the displacement potential, q is a position in Lagrangian coordinates, and the derivatives are taken with respect to this coordinate system. For convenience, we also defined the traceless part of the tidal field: $$K_{ij}(\mathbf{q})=({\mathrm{\partial }}_i{\mathrm{\partial }}_j-\frac{1}{3}{\delta }_{ij}^K){\mathrm{\Phi }}_G.$$(5)

Assuming for now that the tidal field is the only property of the Lagrangian matter field affecting the shapes of gravitationally bound objects, we can now develop a theoretical framework to quantify its effect.

A key assumption of the forthcoming calculations is the validity of the peak-background split (PBS; Bardeen et al. 1986). That is, we assumed that the criterion for forming a halo is defined relative to small-scale quantities; large scales only affect this process by introducing modulations to the distribution of small-scale perturbations. SU simulations provide an exact implementation of this scenario (Gnedin et al. 2011; Wagner et al. 2014; Schmidt et al. 2018; Stücker et al. 2021; Akitsu et al. 2021), and it is interesting to consider them for a moment.

Consider a simulation where a constant perturbation, T₀, to the tidal field is introduced. Upon identifying a set of objects (e.g., halos or galaxies) in these simulations, one can compute their average shape as a function of the large-scale tide T₀ (see Figure 1 for an illustration of this procedure). This yields a functional relation ⟨I|T₀⟩_g – where our notation indicates an expectation value that has been conditioned on T₀ and that has been taken over galaxies (or halos). This function can be expanded in a perturbative series, $${⟨\mathbf{I}\mid {\mathbf{T}}_0⟩}_g={\mathbf{C}}_{K,1}{\mathbf{T}}_0+{\mathbf{T}}_0^T{\mathbf{C}}_{K,2}{\mathbf{T}}_0+\cdots ,$$(6)

where C_K,1 is a rank 2 tensor, C_K,2 is of rank 4 and each product corresponds to a contraction over two indices. Therefore, by running SU simulations at several values of T₀, one can compute derivatives of ⟨I|T₀⟩ and, hence, extract generalized tensorial bias parameters: $${\mathbf{C}}_{K,n}={\frac{{\mathrm{\partial }}^n⟨\mathbf{I}\mid {\mathbf{T}}_0⟩}{\mathrm{\partial }{\mathbf{T}}_0^n}|}_{{\mathbf{T}}_0=0}.$$(7)

SU simulations therefore provide a well-defined route to compute shape bias parameters. Unfortunately, they have the down-side of being costly due to the need to run several simulations at different values of T₀. We now show an alternative route to measuring these bias parameters that is fully consistent with the SU method, using quantities that can be obtained from only a single simulation.

We began by considering the linear field of a simulation in a finite volume, smoothed with a sharp filter in k-space such that only modes with k < k_d are retained. This linear field is characterized at first order by its tidal field, which we denote as T.

If we average the shape of galaxies that inhabit regions of the same tidal field T, we obtain the conditional expectation value ⟨I|T⟩_g. Note that this function is slightly different than ⟨I|T₀⟩_g, since the tidal field measured at the finite scale k_d carries a different amount of information than the idealized infinite scale field T₀. However, using probability theory, we can emulate a change in the field at infinite scales, which modifies the distribution of the field at finite scales.

First of all, consider that if we have an expectation value of some variable h that is conditioned on two probability variables ⟨h|x, y⟩, then we can obtain the expectation value that is only conditioned on the variable x through a weighted average as $$⟨h\mid x⟩=\frac{\int ⟨h\mid x,y⟩p(y\mid x)\mathrm{d}y}{\int p(y\mid x)\mathrm{d}y},$$(8)

where the denominator can be omitted if p(y|x) is a normalized distribution. Additionally, we need to consider that the distribution of interest is not the (Gaussian) conditional distribution at random locations – let us refer to it as p(T|T₀) – but instead it is the distribution at the locations of galaxies p_g(T|T₀). The two are related through a bias function: $$p_g\left(\mathbf{T}\mid {\mathbf{T}}_0\right)=f(\mathbf{T},{\mathbf{T}}_0)p\left(\mathbf{T}\mid {\mathbf{T}}_0\right),$$(9)

where f(T, T₀) describes the relative excess probability of finding a galaxy in a region that has the tidal field T at finite scales and T₀ at infinite scales (see Stücker et al. 2025, for a detailed explanation), $$f(\mathbf{T},{\mathbf{T}}_0)=\frac{p(g\mid \mathbf{T},{\mathbf{T}}_0)}{p(g)},$$(10)

or, through Bayes’ theorem, $$f(\mathbf{T},{\mathbf{T}}_0)=\frac{p(\mathbf{T},{\mathbf{T}}_0\mid g)}{p(\mathbf{T},{\mathbf{T}}_0)}.$$(11)

Under the PBS assumption the value of T₀ is irrelevant if the smaller scale field T is known, so that f(T, T₀) = f(T). Putting things together, we find $${⟨\mathbf{I}\mid {\mathbf{T}}_0⟩}_g=\frac{\int {⟨\mathbf{I}\mid {\mathbf{T}}_0,\mathbf{T}⟩}_{g}f(\mathbf{T})p\left(\mathbf{T}\mid {\mathbf{T}}_0\right)d\mathbf{T}}{\int f(\mathbf{T})p\left(\mathbf{T}\mid {\mathbf{T}}_0\right)d\mathbf{T}}.$$(12)

We can further simplify this expression by using that, under the PBS assumption, ⟨I|T₀, T⟩_g = ⟨I|T⟩_g. Further, we defined the “renormalized bias function”: $$\begin{array}{rl}F\left({\mathbf{T}}_0\right)& =\int f(\mathbf{T})p\left(\mathbf{T}\mid {\mathbf{T}}_0\right)d\mathbf{T}\\ & =⟨f(\mathbf{T}+{\mathbf{T}}_0)⟩,\end{array}$$(13)

where the second line assumes the use of a sharp k-space filter to define T (Stücker et al. 2025). We then find $${⟨\mathbf{I}\mid {\mathbf{T}}_0⟩}_g=\frac{1}{F\left({\mathbf{T}}_0\right)}\int ⟨\mathbf{I}\mid \mathbf{T}{⟩}_{g}f(\mathbf{T})p\left(\mathbf{T}\mid {\mathbf{T}}_0\right)d\mathbf{T}$$(14) $$=\frac{1}{F\left({\mathbf{T}}_0\right)}\int \int \mathbf{I}\frac{p(\mathbf{I},\mathbf{T}\mid g)}{p(\mathbf{T}\mid g)}f(\mathbf{T})p\left(\mathbf{T}\mid {\mathbf{T}}_0\right)d\mathbf{I}d\mathbf{T}$$(15) $$=\frac{1}{F\left({\mathbf{T}}_0\right)}\int \int \mathbf{I}p(\mathbf{I},\mathbf{T}\mid g)\frac{p\left(\mathbf{T}\mid {\mathbf{T}}_0\right)}{p(\mathbf{T})}d\mathbf{I}d\mathbf{T}$$(16) $$=\frac{1}{F\left({\mathbf{T}}_0\right)}{⟨\mathbf{I}\frac{p\left(\mathbf{T}\mid {\mathbf{T}}_0\right)}{p(\mathbf{T})}⟩}_g,$$(17)

where in the second line we have expanded ⟨I|T⟩_g as an integral, in the third line we used Bayes’ theorem f(T) = p(g|T)/ p(g) = p(T|g) /p(T), and we rewrote the double integral over I and T weighted by p(I, T|g) as an average over all galaxies. Equation (17) simply tells us that one can recover the SU average ⟨I|T₀⟩_g by calculating the weighted average of shapes – where the weight corresponds to the relative increase in the probability of finding each galaxy’s environment if a large-scale tidal field T₀ is applied.

For the case of a sharp k-space filter, it is p(T|T₀) = p(T − T₀) so that after substituting into Eq. (7) we find $${\mathbf{C}}_{K,1}={⟨\frac{1}{p(\mathbf{T})}\mathbf{I}\otimes (-{\mathbf{p}}^{(1)}(\mathbf{T})-p{\mathbf{F}}^{(1)}(0))⟩}_g$$(18) $${\mathbf{C}}_{K,2}={⟨\frac{1}{p(\mathbf{T})}\mathbf{I}\otimes ({\mathbf{p}}^{(2)}+2{\mathbf{p}}^{(1)}\otimes {\mathbf{F}}^{(1)}-p{\mathbf{F}}^{(2)}+p{\left({\mathbf{F}}^{(1)}\right)}^2)⟩}_g,$$(19)

where we write the derivatives with respect to T in boldface to remark that this derivative is a 2k-rank tensor in which k is the order of the derivative.

The procedure outlined in this section defines an estimator of the shape biases, which we show to be unbiased and which has some notable properties:

• It can be defined on an object-by-object basis, by evaluating Eq. (18) without the averaging brackets.

• This estimator does not rely on the calculation of n-point functions, making it simpler to evaluate for any selection of objects, even in complex geometries.

• The re-normalization procedure is conceptually very simple in this picture.

Let us expand a bit more on this final point. When measuring bias parameters using n-point functions computed with perturbation theory, an issue appears that these measurements depend on the cutoff scale, Λ, used to divide large-from and small-scale modes (see Assassi et al. 2014; Desjacques et al. 2018; and Rubira & Schmidt 2024 for a detailed explanation). Similarly, our individual bias estimates also depend on the scale at which we damped our density field, k_d. However, after we average those estimates weighted by the probability density function (PDF) of the matter perturbations, this dependence should be approximately removed. This procedure should work well as long as it is a good approximation that the distribution of T_ij is a multivariate Gaussian, and that its cross-correlations with other observables are small. An issue appears of course when we make k_d similar to 1/R_*, where we denote by R_* the typical Lagrangian scale of halos. At these scales, higher-derivative operators become relevant, and one needs to consider their cross-correlation with T_ij to properly estimate the values of bias-parameters.

The properties mentioned above are also present in the estimators derived in the following sections of this manuscript.

Fig. 1 Illustration of halo shapes in a simulation slice and their relation to the underlying density and tidal fields. Left: density field linearly extrapolated to z = 0 (background), with a large set of small arrows superimposed. The sizes of these small arrows are scaled to the local projected tidal field’s largest eigenvalue, and their direction is that of the corresponding eigenvector. Middle: halos (colored ellipses) centered at their Lagrangian positions. Their shapes are those of the projected shape tensor estimated for the FoF group; their color corresponds to the value of their three-dimensional shape-alignment bias. Right: same as the middle panel except now the halos are at their final Eulerian positions.

2.1 Linear bias

The derivation leading to Eq. (18) has been left so far in a very general form; we now connect it to more common definitions of shape biases found in the literature.

From Eq. (7) taken at n = 1, we can see that the values assumed by the rank-4 tensor C_K,1 are constrained by the symmetries of the right-hand side, that is, C_K,1 has to be isotropic, symmetric over its first two indices (given the symmetry of ⟨I|T₀⟩) and the last two as well (given the symmetry of T). In Appendix B we summarize the findings of Stücker et al. (2025), who derived an orthogonal basis for all isotropic rank-4 tensors that respect the imposed symmetries, 𝕍₂₂ = span({J₂₂, J₂₌₂}), in which the tensors J₂₂ and J₂₌₂ are defined by $${\mathbf{J}}_{22}=S_{22}\left({\delta }_{ij}^K{\delta }_{kl}^K\right)$$(20) $${\mathbf{J}}_{2=2}=S_{22}\left({\delta }_{ik}^K{\delta }_{jl}^K\right)-\frac{1}{3}{S}_{22}\left({\delta }_{ij}^K{\delta }_{kl}^K\right),$$(21)

where S₂₂ is a double symmetrization operator, defined via its action on an arbitrary rank-4 tensor M, given by $$S_{22}\left(M_{ijkl}\right)=\frac{1}{4}(M_{ijkl}+M_{jikl}+M_{ijlk}+M_{jilk}).$$(22)

It is useful to keep in mind that the notation we employ gives information on the result of contracting these tensors with any arbitrary rank-2 tensors A and B. Contracting J₂₂ with A and B simply gives us the product of their traces – each subscript “2” without any connection symbol represents taking the trace of a rank-2 tensor. As for J₂₌₂, contracting it with A and B gives the trace of their matrix product; each “–” symbol connecting the subscripts represents summation over one repeated index; in this case, the first one is due to the matrix product, and a second one must appear to form the trace of this product.

Using this basis, one can uniquely decompose C_K,1 as a linear combination, $${\mathbf{C}}_{K,1}=c_{J_{22}}{\mathbf{J}}_{22}+c_{J_{2=2}}{\mathbf{J}}_{2=2},$$(23)

and, by contracting with the basis tensors and taking advantage of orthogonality, we can write the expressions for the bias parameters in Eq. (32), $$c_{J_{22}}=-{⟨\frac{1}{p(\mathbf{T})}(\frac{\mathrm{\partial }p}{\mathrm{\partial }\mathbf{T}}-p\frac{\mathrm{\partial }F}{\mathrm{\partial }\mathbf{T}}(0))\otimes {\mathbf{I}}_{ \cdot }^{(4)}\frac{{\mathbf{J}}_{22}}{{\parallel {\mathbf{J}}_{22}\parallel }^2}⟩}_{\mathrm{g}}$$(24) $$c_{J_{2=2}}=-{⟨\frac{1}{p(\mathbf{T})}(\frac{\mathrm{\partial }p}{\mathrm{\partial }\mathbf{T}}-p\frac{\mathrm{\partial }F}{\mathrm{\partial }\mathbf{T}}(0))\otimes {\mathbf{I}}_{ \cdot }^{(4)}\frac{{\mathbf{J}}_{2=2}}{{\parallel {\mathbf{J}}_{2=2}\parallel }^2}⟩}_{\mathbf{g}}.$$(25)

To evaluate these expressions, we needed to examine the form of the PDF for the Lagrangian tidal field p(T). Stücker et al. (2025) show that the PDF for the tidal field can be written as $$p(\mathbf{T})=N\exp [-\frac{1}{2}\mathbf{T}{\mathbf{C}}_T^{†}\mathbf{T}],$$(26)

where N is a normalization and ${\mathbf{C}}_T^{†}$ is the pseudo-inverse of the covariance of the tidal field, given by $${\mathbf{C}}_T^{†}=\frac{1}{{\sigma }^2}{\mathbf{J}}_{22}+\frac{15}{2{\sigma }^2}{\mathbf{J}}_{2=2},$$(27)

where σ is the variance of the density-field, damped at the scale of choice. From this expression, one can straightforwardly compute the first derivative of p with respect to the tidal field, given by $$\frac{1}{p}\frac{\mathrm{\partial }p}{\mathrm{\partial }\mathbf{T}}=-{\mathbf{C}}_T^{†}\mathbf{T}.$$(28)

Therefore, one can rewrite Eqs. (24) and (25) as $$c_{J_{22}}={⟨{\mathbf{C}}_T^{†}\mathbf{T}\otimes {\mathbf{I}}_{ \cdot }^{(4)}\frac{{\mathbf{J}}_{22}}{{\parallel {\mathbf{J}}_{22}\parallel }^2}⟩}_{\mathrm{g}}=\frac{1}{{\sigma }^2}⟨\mathrm{Tr}(\mathbf{T})\mathrm{Tr}(\mathbf{I}){⟩}_g=0$$(29) $$c_{J_{2=2}}={⟨{\mathbf{C}}_T^{†}\mathbf{T}\otimes {\mathbf{I}}_{ \cdot }^{(4)}\frac{{\mathbf{J}}_{2=2}}{{\parallel {\mathbf{J}}_{2=2}\parallel }^2}⟩}_{\mathrm{g}}=\frac{3}{2{\sigma }^2}⟨\mathrm{Tr}(\mathbf{K}\mathbf{I}){⟩}_g,$$(30)

in which the first term is equal to zero because we defined I to be trace-free; the first bias parameter, c_J2=2 is exactly equivalent to the more commonly defined c_K, that is, the response of galaxy shapes to perturbations in the traceless part of the tidal field.

2.2 Second-order biases

We can see from Eq. (6) that the rank-6 tensor C_K,2 is constrained by the symmetries to be isotropic and symmetric over its first, second, and third pairs of indices. As explained in Appendix B, the set of isotropic rank-6 tensors that respect the imposed symmetries can be generated by the following orthogonal basis: $${\mathbb{V}}_{222}=\mathrm{span}\left(\right\{{\mathbf{J}}_{222},{\mathbf{J}}_{2-2-2-},{\mathbf{J}}_{2=22},{\mathbf{J}}_{22=2},{\mathbf{J}}_{222}=\left\}\right).$$(31)

This implies that one can uniquely decompose C_K,2 as a linear combination of these tensors: $$\begin{array}{rl}{\mathbf{C}}_{K,2}=& c_{J_{222}}{\mathbf{J}}_{222}+c_{J_{22=2}}{\mathbf{J}}_{22=2}+c_{J_{2=22}}{\mathbf{J}}_{2=22}\\ & +c_{J_{222}}{\mathbf{J}}_{222}+c_{J_{2-2-2}}{\mathbf{J}}_{2-2-2-}.\end{array}$$(32)

By contracting with the basis tensor, one obtains the bias parameters $$c_{J_{2=22}}={⟨\frac{1}{p(\mathbf{T})}(p\frac{{\mathrm{\partial }}^{2}F}{\mathrm{\partial }{\mathbf{T}}^2}-2\frac{\mathrm{\partial }p}{\mathrm{\partial }\mathbf{T}}\frac{\mathrm{\partial }F}{\mathrm{\partial }\mathbf{T}}+\frac{{\mathrm{\partial }}^{2}p}{\mathrm{\partial }{\mathbf{T}}^2})\otimes {\mathbf{I}}_{ \cdot }^{(6)}\frac{{\mathbf{J}}_{2=22}}{{\parallel {\mathbf{J}}_{2=22}\parallel }^2}⟩}_{\mathrm{g}}.$$(33)

In the rest of this work, we assumed that the six-dimensional dot product is performed in such a way that I gets contracted with the last two indices of J_X and the first four are contracted with the derivatives of p and F. This is relevant because terms that have as their last index a disconnected “2” subscript, such as the one we have written, are proportional to the trace of I, which is zero. The same does not occur for terms for which the disconnected “2” subscript is contracted with the derivatives. We refrain from writing all of the terms for conciseness because they are completely analogous.

We now turn to the evaluation of the expression in parenthesis in Eq. (33). The second derivative of the PDF of the tidal tensor can be obtained straightforwardly $$\frac{1}{p}\frac{{\mathrm{\partial }}^{2}p}{\mathrm{\partial }{\mathbf{T}}^2}=\frac{1}{4}\left[{\mathbf{C}}_T^{†}(\mathbf{T}+{\mathbf{T}}^T)\right]\otimes \left[{\mathbf{C}}_T^{†}(\mathbf{T}+{\mathbf{T}}^T)\right]-{\mathbf{C}}_T^{†},$$(34)

and its first derivative has already been discussed in Eq. (28). As for F, the density bias function, we can write it explicitly as $$F(\mathbf{T})=1+{\mathbf{B}}_{T,1}\mathbf{T}+\frac{1}{2}\mathbf{T}{\mathbf{B}}_{T,2}\mathbf{T}+\frac{1}{6}{\mathbf{B}}_{T,3}\mathbf{T}\mathbf{T}\mathbf{T}+\cdots ,$$(35)

in which B_T,N are the generalized N-th order density bias parameters, which are defined as $${\mathbf{B}}_{T,1}=b_1{\mathbf{J}}_2$$(36) $${\mathbf{B}}_{T,2}=b_2{\mathbf{J}}_{22}+2b_{K^2}{\mathbf{J}}_{2=2}$$(37) $${\mathbf{B}}_{T,3}=b_3{\mathbf{J}}_{222}+b_{\delta K^2}{\mathbf{J}}_{22=2}+b_{K^3}{\mathbf{J}}_{2-2-2-}.$$(38)

Its derivatives are given by $$\frac{\mathrm{\partial }F}{\mathrm{\partial }\mathbf{T}}(0)={\mathbf{B}}_{T,1}$$(39) $$\frac{{\mathrm{\partial }}^{2}F}{\mathrm{\partial }{\mathbf{T}}^2}(0)={\mathbf{B}}_{T,2}.$$(40)

Not all of these tensors provide relevant contributions when substituted into Eq. (33). Notice that, according to the convention adopted in this work, the last two indices of these rank-6 tensors are contracted with I, and therefore J₂₂₂ and J₂₌₂₂ are guaranteed to be zero since they would extract the trace of I, which is null. Furthermore, J₂₌₂₂ and J₂₂₂₌ result in identical results due to the symmetry in permutation of the first and second pairs of indices in Eq. (33). Therefore, the relevant tensors, which define our second-order bias parameters, are $${\mathbf{J}}_{2=22},{\mathbf{J}}_{2-2-2-}.$$(41)

Substituting the results of Eqs. (28), (34), (39), and (40) into Eq. (33), we can work out the expressions for each of these bias parameters. These calculations are cumbersome; therefore, we developed a semi-numerical method for computing them. The code used for these calculations is made public in a GitHub repository¹. We simply provide the results obtained from those calculations, $$c_{\delta K}=\frac{3}{2}{⟨\delta \frac{\mathrm{Tr}(\mathbf{K}\mathbf{I})}{{\sigma }^4}⟩}_g-b_1{c}_K$$(42) $$c_{K\otimes K}=\frac{135}{7}{⟨\frac{\mathrm{Tr}(\mathbf{K}\mathbf{K}\mathbf{I})}{{\sigma }^4}⟩}_g.$$(43)

It is interesting to spend a moment on Eq. (42). Notice that this estimator, unlike the two others, depends on quantities that are defined in an ensemble average, notably c_K and b₁, implying that, with this definition, it is not possible to estimate this bias parameter on an object-to-object basis.

The usual bias expansion for IAs at second order has three operators, δK_ij, (K ⊗ K)_ij, and t_ij, defined by $$t_{ij}(\mathit{q})=-\frac{4}{21}[\frac{{\mathrm{\nabla }}_i{\mathrm{\nabla }}_j}{{\mathrm{\nabla }}^2}-\frac{1}{3}{\delta }_{ij}^K]({\delta }^2-\frac{3}{2}{K}^2).$$(44)

The first two operators are included in our analysis, and we derived estimators for their corresponding bias parameters, but the same is not true for t_ij. We can see from its expression in Eq. (44) that it cannot be written linearly in terms of T_ij T_lm, due to the presence of the inverse Laplacian term $\frac{1}{{\mathrm{\nabla }}^2}$, and therefore its PDF is not Gaussian, and the derivative of the PDF with respect to t_ij cannot be built as a linear combination of the tensor $$C_{2,ijkl}=\frac{{\mathrm{\partial }}^{2}p}{\mathrm{\partial }{T}_{ij}\mathrm{\partial }{T}_{kl}}.$$(45)

Therefore, to derive an estimator for the corresponding bias parameter, c_t, we would have to compute the joint PDF of K_ij and t_ij. Even if this could in principle be done, we choose to defer it to future work as it would not add significantly to the results presented here.

2.3 Linear density bias

We are also interested in measuring the linear density bias of halos, b₁, to understand its coevolution with the shape bias parameters. The method we used for measuring this parameter is very similar in spirit to the ones described in this article, and is discussed extensively in Stücker et al. (2025). We refer the reader to that article for more information and suffice with quoting the equation that defines the estimator for b₁, $$b_1={⟨\frac{\delta }{{\sigma }^2}⟩}_g.$$(46)

3 Simulations

The simulations employed in this work are gravity-only N-body simulations run using a modified version of L-GADGET3 (Springel 2005; Angulo et al. 2021), with parameters summarized in Table 2. The larger simulation, which we refer to as the “Planck” simulation, has a cosmology compatible with that found by Planck Collaboration VI (2020). As for the two smaller ones, they have cosmological parameters specifically designed to reduce the errors of emulators in cosmological parameter space (Angulo et al. 2021). The initial conditions are created at z = 49 using fixed amplitudes (Angulo & Pontzen 2016), and then evolved to z = 0. At each output redshift, the code stores properties of friends-of-friends (FoF) and SUBFIND groups, which are interpreted to be halos and subhalos, respectively. Properties of these halos and subhalos are then computed from these particle distributions, using many different definitions and criteria, which we make explicit in the following subsection.

3.1 Shape definitions

Halos can be defined in many different ways, and the resulting measurements of quantities such as mass or shape are highly dependent on the precise definition used. We list in Table 1 a few definitions of shapes that are used throughout this manuscript. We generally used the FoF definition for the halo shapes and M_200,b as the mass of our halos. This choice was made for convenience since these properties were natively computed in our simulations, while the other ones had to be recomputed in postprocessing.

3.2 Shape convergence

Another important matter is to determine the minimum number of particles a halo should have so that its shape-measurement is well converged. We estimated this limit by measuring shape bias parameters for the two versions of the Narya simulation with different mass resolutions, described in Table 2. The results of these measurements can be seen in Fig. C.1 and show that value of shape-biases are well converged for halos with M_h > 10¹² M_⊙h⁻¹. These halos are resolved with ~270 particles in the low-resolution Narya simulation. Consequently, in the remainder of this manuscript, we only analyze halos resolved with at least 300 particles.

3.3 Error calculation

Before entering into the description of our results, we describe how we estimated uncertainties in our measurements of the bias parameters. We adopted the jackknife technique, dividing our simulation volume into N_J = 4³ = 64 subregions, and then evaluated N_J by means of the bias parameters, each time leaving out one of the sub-volumes. We computed the covariance and the mean of our estimate as $$\mathbf{C}=\frac{1}{N_J-1}\sum _n({\mathbf{b}}_n-{\mathbf{b}}_0)\otimes ({\mathbf{b}}_n-{\mathbf{b}}_0)$$(47) $${\mathbf{b}}_0=\frac{1}{N_J}\sum _n{\mathbf{b}}_n,$$(48)

where b_n represents a vector of the three bias parameters averaged in subregion n. Finally, the uncertainty on c_K was computed as ${\sigma }_{c_K}=\sqrt{{\mathbf{C}}_{00}}$. We represent these uncertainties as error bars that extend up and down from the points with a length of σ_c, except in Fig. 2 where we display these errors with a region filled one σ_c upward and one downward of the average curve.

Fig. 2 Halo bias parameters as a function of mass in the Planck simulation. Different colors indicate the damping scale (kd), measured in [hMpc−1], used to smooth the Lagrangian fields employed in the measurement of the bias parameters, and the colored regions indicate the 1σ interval around the mean value. Values presented here were measured for halos at z = 0, and the measurements at higher redshifts are even less dependent on kd, so any conclusion on the convergence of the estimators obtained here are also valid at other times.

4 Results

In this section we present the results of applying our shape bias estimators to halos in the simulations described in the previous section. Before presenting the results, we emphasize that we can estimate c_K and c_K⊗K for each individual object. This gives us great flexibility since once this set of biases is measured, we can average them over any quantity of interest, provided it is a small-scale property of the halos. To avoid confusion, we state clearly in the figure caption the quantity used for binning since it might not be the one being displayed.

4.1 Validation of estimators

As described in Sect. 2, our bias estimators require smoothing the Lagrangian density field on an arbitrary scale k_d in Fourier space. An important validation is to make sure that our estimates are, in fact, independent of k_d if we make it sufficiently small (since bias parameters are formally scale-independent quantities). Figure 2 shows the measurements of c_K, c_δK and c_K⊗K averaged in bins of halo mass for increasing values of k_d. These measurements are displayed in the figure as shaded bands marking the 1σ region around the mean value of the biases. One can see that the measurements using k_d = 0.05 [hMpc⁻¹] and k_d = 0.1 [hMpc⁻¹] are compatible within 1σ over all masses, with the exception of the largest mass bin for c_K⊗K, for which these estimators are compatible within 1.5σ. Based on these checks, we adopted k_d = 0.1 [hMpc⁻¹] as the smallest scale at which our estimators are scale-independent and used it to estimate all bias parameters that appear in the following.

To validate our method, we compared the results to previous measurements in the literature (Stücker et al. 2021; Akitsu et al. 2021, 2023). In Fig. 3 we compare c_δK and c_K⊗K versus c_K for halos in different mass bins. Different colored lines represent different halo-shape definitions. Comparing the results for our standard definition, I_FoF, to the measurements of Akitsu et al. (2023) shows similar trends but a significant quantitative disagreement. However, the shapes of halos are expected to be very dependent on the halo definition used. Shape definitions based on spherical overdensity, such as ℐ_200,c, ℐ_500,c, and ℐ_200,b, produce more isotropic halos due to the fact that the particles belonging to the halo are preselected to be in a spherical region. On the other hand, the FoF shape definition allows the halo out-skirts to have varying shapes, which can heavily influence the measurement of ℐ_FoF.

If we compare the measurements of Akitsu et al. (2023) to our results for spherical boundary definitions – in particular the reduced inertia tensor ${\tilde{\mathcal{I}}}_{200,\mathrm{b}}$ – we find a good agreement. Since these authors use the Amiga’s halo finder (Knollmann & Knebe 2009), a spherical-overdensity based halo-finder, to define their halos and compute their bias parameters with the reduced shape tensor, this is the comparison that makes the most sense. The agreement is particularly good for ${\tilde{\mathcal{I}}}_{200,\mathrm{b}}$, which is unsurprising since those authors employ the same definition for the halo boundary. Therefore, we consider our measurements to be consistent with those reported by Akitsu et al. (2023), but we note that care must be taken when comparing different halo definitions – much more so than, for example, when comparing halo density bias parameters. As the tidal alignment of halos is not directly observationally relevant, the preferential definition is not clear, and we stuck to the FoF boundary for simplicity. However, we note that, for example for measuring the alignment of galaxies in simulations, it is advisable to use definitions that match observational criteria (Tenneti et al. 2015b).

4.2 Universal relations

As mentioned at the beginning of this section, our method has the advantage of being very flexible and allowing the measurement of shape bias parameters averaged over any quantity of choice. Figure 4 shows the shape bias parameters measured from the halos in the Planck and Narya-HD simulations, averaged directly in peak height ν = δ_c/σ(M, z), where σ(M, z) is defined as $$\sigma (M,z{)}^2=\int \frac{d^3\mathbf{k}}{(2\pi {)}^3}P(k,z)W^2(kR),$$(49)

where $R={\left(\frac{3M}{4\pi \overline{\rho }}\right)}^{1/3}, \overline{\rho }$ is the mean density of the universe, and W is a top-hat profile in Fourier space: $$W(kR)=\frac{1}{(kR{)}^2}(\frac{3 \sin (kR)}{kR}- \cos (kR)).$$(50)

This quantity is commonly used in excursion-set theory to measure how likely it is that the density in a region of radius R has crossed δ_c, the threshold for collapse. The relations in terms of ν are universal, that is, they do not depend on redshift or cosmology. Furthermore, we notice that the functional forms relating the shape biases to ν can be described by simple linear or quadratic functions, which get damped at very low ν, much simpler than those relating the shape biases to b₁. Fits for the relation of the IA biases with ν are provided in Table 3. For comparison, we also show these relations with b₁ in Fig. 5 with fits as listed in Table 4.

The fact that c_K depends strongly on ν is very revealing about the nature of alignments. Let us consider the hypothesis that halos get their alignments through the dynamical effect of the tidal field after formation. The tidal field may slightly distort the virialized halos – shifting them to an anisotropic equilibrium configuration. We refer to this as the “post-formation-response” scenario. Such a mechanism has been considered in several studies to explain the IA of galaxies (e.g., Tugendhat & Schäfer 2018; Zjupa et al. 2020; Ghosh et al. 2024), but it could also be a plausible effect for halos.

If post-formation-response is the primary alignment mechanism of halos, then the degree of alignment should depend primarily on the anisotropy of the gravitational field in the vicinity of the halo. For a spherical halo with mass profile M(r) plus a traceless large-scale tide K, the gravitational acceleration is given by $$\mathit{a}(\mathit{r})=4\pi G{\rho }_m\mathbf{K}\mathit{r}-\frac{GM(r)}{r^3}\mathit{r},$$

where the factor 4πGρ_m is needed to correct for the dimensionless form of K. If we consider an Navarro–Frenk–White halo at the virial radius r_200c, we can define the effective tidal tensor at the scale of the halo: $$\begin{array}{rl}{T}_{200\mathrm{c}}& =4\pi G{\rho }_{\mathrm{m}}\mathbf{K}-\frac{G{M}_{200\mathrm{c}}}{r_{200\mathrm{c}}^3}{\mathbf{J}}_2\\ & =4\pi G{\rho }_{\mathrm{m}}(\mathbf{K}-\frac{200}{3}{\mathbf{J}}_2),\end{array}$$(51)

where J₂ is the unit matrix. The isotropic component is sourced purely by the halo, and the anisotropic component purely by the large-scale tide. Therefore, the relative anisotropy in the gravitational field of an already formed halo is on the order 3K₁/200 (where K₁ is the largest eigenvalue of K). The anisotropy in the triggered density response should be proportional to this and should probably also not be larger in amplitude. Therefore, we expect a post-formation-response that is mass-independent and of order |c_K| ≲ 10⁻². This simple expectation clearly disagrees with the measured c_K relations, which exhibit a strong mass (or ν) dependence and a much larger degree of alignment up to |c_K| ≲ 0.5. We conclude, therefore, that post-formation-response is not the dominant mechanism for the alignment of halos (see Camelio & Lombardi 2015, for a similar consideration for the case of galaxies). The anisotropy must already be encoded during the formation and accretion processes.

We speculated about three possible scenarios that could explain how the ν dependence of c_K is encoded during formation. (1) The tidal field only significantly introduces anisotropy up to the formation of a halo and becomes irrelevant at later times. In this scenario, halos that form earlier have a lesser response because they have formed in a weaker tidal field. Combined with the fact that a smaller ν leads to an earlier formation, this could qualitatively explain the observed trend. However, it is unclear whether this could match the measurements quantitatively and whether it is consistent with the observed universality. (2) Halos at low ν (and therefore large variance) are born in significantly larger tidal fields than halos at high ν. The relative changes in the tidal tensor at the halo scale due to an imposed large-scale tidal field are, therefore, much smaller for small halos. (3) Halo shapes are not responding as a continuous deformation to the large-scale tidal field but are rather governed by discrete events – for example through the formation of larger-scale pancakes and filaments. The response of small halos may then be small because they are almost all embedded in larger-scale filaments, so it is unlikely that an increase in the tidal field triggers an additional filament formation. Such a process could be described through an excursion set formalism, which naturally leads to universality.

While so far a clear understanding of the formation of shapes is elusive, we note that there may be some hope to gain further insights: (a) The simple universal form of the response motivates the notion that a simple description exists – similar to the way that excursion sets offer a quite reasonable model for the density biases of halos. (b) All the abovementioned scenarios could be turned into simple models with quantitative predictions that can be easily tested. (c) The response of halos could be investigated in warm dark matter simulations. The smallest warm dark matter halos have a clear and well-defined formation time and size (Diemand et al. 2005; Ishiyama et al. 2010; Angulo et al. 2017; Ogiya & Hahn 2018; Ondaro-Mallea et al. 2024). It would be insightful to explore how such objects respond to tides before and after formation. (d) Some of the points above can be studied in anisotropic SU simulations. For example, for (2) to be correct, the differential response |∂/∂T₀⟨I|T₀⟩_g| at large values of the tidal field should be notably smaller than the response at T₀ = 0.

Fig. 3 Comparison of shape-bias measurements performed in the Planck simulation with a range of different halo definitions to the measurements reported in Akitsu et al. (2023). There is a large variability in the relations found, of which ${\tilde{\mathcal{I}}}_{200,\mathrm{b}}$ gives the best agreement with previous works.

Fig. 4 Measurements of shape bias parameters in the Planck simulation as a function of peak-significance ν = δc/σ(M, z). All measurements, regardless of redshift, seem to follow a common trend with ν. Solid black lines show the fits performed to approximate these common relations. These fitting functions are exceptionally simple – linear for cK and cK⊗K and quadratic for cδK – and provide good descriptions of the measured relations.

Fig. 5 Shape biases in the Planck simulation averaged over bins in ν = δc/σ(M, z), the peak-significance, as a function of linear density bias (b1), also averaged over bins of ν. Different colors represent the redshifts indicated in the legend. Each panel shows the orthogonal distance of the measurements to the fit, in units of the measurement error. The postulated functions provide a good description of the relations between the shape biases and b1.

4.3 Assembly bias

In addition to the dependence of halo IA bias on ν, it is interesting to explore whether splitting the same ν populations by some secondary property would result in significantly distinct signals. We explored this possibility by splitting our samples according to their spin, which is defined as (Bullock et al. 2001) $$\lambda =\frac{|\mathbf{L}|}{\sqrt{2}{v}_{200,c}{r}_{200,c}},$$(52)

where $$\mathbf{L}=⟨\mathbf{r}\times \mathbf{v}{⟩}_{\text{particles}},$$(53)

is the specific angular momentum vector of the object. We then split the population at each bin of ν into four quartiles of the spin distribution such that Q₁ contains the halos with the 25% smallest spins and Q₄ the 25% largest. Figure 6 shows the ratios between the biases of these subpopulations and the bias of the full sample. We can see that at low masses, the amplitude of $c_K^Q$ is anticorrelated with spin: lower-spin halos have a higher bias and higher-spin halos have a lower bias. This trend is inverted around ν ≈ 1.75, and beyond that the amplitude of c_K is larger for high-spin halos and smaller for low-spin ones. As far as we are aware, this is the first time that a secondary dependence of c_K on halo spin is observed.

It is well known that the linear density bias of halos b₁ has a secondary dependence on spin (Gao & White 2007), and hence it is interesting to probe whether the dependence of c_K with spin is simply a consequence of the dependence of b₁ with spin. With that purpose, in Fig. 7, we show c_K and b₁ averaged over ν bins and split into subpopulations by spin, as described earlier. One can clearly see the assembly bias in b₁ – measurements that have the same ν are connected by dashed lines, and they have different values of b₁ depending on the spin-quartile they belong to. Nevertheless, one can also identify an independent c_K assembly bias signal since c_K does not vary along the b₁ − c_K relation at fixed ν with varying spin.

5 Conclusions

In this article we have presented a new methodology for computing shape biases in cosmological simulations from correlations in the distribution of the initial tidal tensor and the final shapes of objects.

Our method is conceptually similar to the anisotropic SU approach (Stücker et al. 2021; Akitsu et al. 2021) and provides a robust measurement of the shape bias parameters that is independent of the number of considered parameters – in notable contrast to methods that need to fit a model to power spectra or even field-level measurements (Zennaro et al. 2022; Schmittfull et al. 2019). We have shown that for sufficiently large damping scales, k_d ≤ 0.1h/Mpc, our method is independent of the damping scale and matches the results from SU simulations. However, the new method allows us to perform these measurements at a much reduced computational cost for an arbitrary set of parameters from a single simulation – making it viable for application in modern hydrodynamical simulations.

The application of these estimators to simulated halos has allowed us to quantify the dependence of the bias parameters c_K, c_δK, and c_K⊗K on the peak height (ν) and linear density bias (b₁), and to provide fitting functions that accurately describe these relations. The strong mass dependence of the alignment response shows that the anisotropy must already be encoded during the formation and accretion history of objects. The post-formation-response – as it is considered in some models of galaxy alignment (e.g., Tugendhat & Schäfer 2018; Zjupa et al. 2020; Ghosh et al. 2024) – is too weak and too mass-independent to explain the linear large-scale alignment of halos.

Further, we have presented the first detection of the dependence of the tidal alignment response (c_K) on the spin of halos at a fixed peak height. This secondary dependence is rather complicated, leading to a greater response for higher spins at high masses but a lesser response for higher spins at lower masses. Understanding this relation may pose a significant challenge for analytical approaches to modeling alignments and tidal torquing (e.g., White 1984).

The response of halo shapes to the large-scale tidal field provides important qualitative insights into collisionless systems and the ingredients that are necessary to model IAs. However, quantitatively, it is far more relevant to understand the alignment of galaxies – in particular, to correctly interpret weak-lensing observations. While we have validated our method in the case of halos in this article, it can easily be applied to simulated galaxies as well – an avenue that we will explore in future studies.

Fig. 6 Measurements of linear shape bias in the Planck simulation, averaged over bins of ν, for halos selected in the lowest and highest quartiles in the distribution of spins. A clear assembly bias signal is detected: for low values of ν, objects with low (high) spin have an above (below) average cK value; this behavior is inverted around ν* = 1.75, as marked by the vertical black line.

Fig. 7 Linear shape bias as a function of linear density bias in the Planck simulation. Different colors represent the quartiles in the distribution of halo spins. One can clearly see a difference in the values of cK even at fixed b1, which indicates a nontrivial assembly-bias signal.

Acknowledgements

The authors thank Kazuyuki Akitsu for providing his measurements at our request. The authors acknowledge support from the Spanish Ministry of Science under grant number PID2021-128338NB-I00 and from the European Research Executive Agency HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Sklodowska-Curie grant agreement number 101086388 (LACEGAL).

Appendix A Notation

In this section we define some conventions in notation to be used throughout this work. This is particularly relevant because we define some convenient ways to define operations between tensors that appear frequently in this paper. Table A.1 resumes the definitions for most of the symbols employed throughout this article. Table A.2 resumes several of the most important notation for tensor operations used throughout this manuscript.

Appendix B Isotropic tensors

In this appendix we condense the results of Stücker et al. (2025) by determining an orthogonal basis for all isotropic tensors of rank n that respect the following symmetry: let M_{i0i1...i2n,i2n+1...} be an isotropic tensor. Then we required M to be symmetric using the interchange i_2n ↔ i_2n+1. For tensors of rank 1, there is only one such tensor, the Kronecker delta symbol, $${\mathbf{J}}_2={\delta }_{ij}^K.$$(B.1)

For rank-4 tensors, there are two such tensors, $${\mathbf{J}}_{22}={\delta }_{ij}^K{\delta }_{kl}^K$$(B.2) $${\mathbf{J}}_{2=2}=S_{22}\left({\delta }_{ij}^K{\delta }_{kl}^K\right)-\frac{1}{3}{\mathbf{J}}_{22};$$(B.3) (B.4)

the subscripts indicate properties of these tensors, namely that by contracting J₂₂ with two rank-2 tensors A and B, one obtains the product of their traces – so there is no mixing of their components. By performing the same calculation with J₂₌₂, on the other hand, one obtains the trace of their matrix product. In fact, each “-” symbol in the subscripts represents one index contraction between tensors A and B. Finally, for rank 6 there are five different tensors $${\mathbf{J}}_{222}={\delta }_{ij}^K{\delta }_{kl}^K{\delta }_{mn}^K$$(B.5) $${\mathbf{J}}_{22=2}=S_{222}\left({\delta }_{ij}^K{\delta }_{kn}^K{\delta }_{ml}^K\right)-\frac{1}{3}{\mathbf{J}}_{222}$$(B.6) $${\mathbf{J}}_{2=22}=S_{222}\left({\delta }_{il}^K{\delta }_{kj}^K{\delta }_{mn}^K\right)-\frac{1}{3}{\mathbf{J}}_{222}$$(B.7) $${\mathbf{J}}_{222}=S_{222}\left({\delta }_{in}^K{\delta }_{kl}^K{\delta }_{mj}^K\right)-\frac{1}{3}{\mathbf{J}}_{222}$$(B.8) $${\mathbf{J}}_{2-2-2-}=S_{222}\left({\delta }_{ik}^K{\delta }_{lm}^K{\delta }_{nj}^K\right)-\frac{1}{9}{\mathbf{J}}_{222}-\frac{1}{3}({\mathbf{J}}_{22=2}+{\mathbf{J}}_{2=22}+{\mathbf{J}}_{222=}).$$(B.9)

Appendix C Shape convergence

In this appendix we present measurements of the shape-bias parameters obtained from two simulations, both run with Narya cosmology but at different resolutions. This test allows us to define a convergence criterion, which is used throughout the entire article. As can be seen from Fig. C.1, these bias parameters are in agreement with each other for M_h > 10¹²M_⊙h⁻¹. This mass corresponds to halos with roughly 270 particles in the low-resolution Narya simulation, so we defined our convergence criterion to be that halos must have more than 300 particles.

Fig. C.1 Comparison of measurements of shape-biases cK, cδK, and cK⊗K in two versions of the Narya simulation, with increasing resolution. For these measurements, shapes were computed using the FoF halo definition. One can see that the results generally match within the error bars, except for measurements of cK below M200,b = 1012M⊙h−1. The damping scale employed in measuring these bias parameters is kd = 0.1hMpc−1

Appendix D Comparison to the literature

Similar measurements to the ones reported in this article were performed by Akitsu et al. (2023), but for the Eulerian bias parameters. Nevertheless, in that work, they present expressions for translating the bias parameters from Eulerian to Lagrangian space. Therefore, when comparing our results to theirs, we always converted our measurements using Eq. (2.15) of Akitsu et al. (2023).

References

All Tables

All Figures

Fig. 1 Illustration of halo shapes in a simulation slice and their relation to the underlying density and tidal fields. Left: density field linearly extrapolated to z = 0 (background), with a large set of small arrows superimposed. The sizes of these small arrows are scaled to the local projected tidal field’s largest eigenvalue, and their direction is that of the corresponding eigenvector. Middle: halos (colored ellipses) centered at their Lagrangian positions. Their shapes are those of the projected shape tensor estimated for the FoF group; their color corresponds to the value of their three-dimensional shape-alignment bias. Right: same as the middle panel except now the halos are at their final Eulerian positions.

In the text

Fig. 2 Halo bias parameters as a function of mass in the Planck simulation. Different colors indicate the damping scale (kd), measured in [hMpc−1], used to smooth the Lagrangian fields employed in the measurement of the bias parameters, and the colored regions indicate the 1σ interval around the mean value. Values presented here were measured for halos at z = 0, and the measurements at higher redshifts are even less dependent on kd, so any conclusion on the convergence of the estimators obtained here are also valid at other times.

In the text

	Fig. 3 Comparison of shape-bias measurements performed in the Planck simulation with a range of different halo definitions to the measurements reported in Akitsu et al. (2023). There is a large variability in the relations found, of which ${\tilde{\mathcal{I}}}_{200,\mathrm{b}}$ gives the best agreement with previous works.
In the text

Fig. 4 Measurements of shape bias parameters in the Planck simulation as a function of peak-significance ν = δc/σ(M, z). All measurements, regardless of redshift, seem to follow a common trend with ν. Solid black lines show the fits performed to approximate these common relations. These fitting functions are exceptionally simple – linear for cK and cK⊗K and quadratic for cδK – and provide good descriptions of the measured relations.

In the text

Fig. 5 Shape biases in the Planck simulation averaged over bins in ν = δc/σ(M, z), the peak-significance, as a function of linear density bias (b1), also averaged over bins of ν. Different colors represent the redshifts indicated in the legend. Each panel shows the orthogonal distance of the measurements to the fit, in units of the measurement error. The postulated functions provide a good description of the relations between the shape biases and b1.

In the text

	Fig. 6 Measurements of linear shape bias in the Planck simulation, averaged over bins of ν, for halos selected in the lowest and highest quartiles in the distribution of spins. A clear assembly bias signal is detected: for low values of ν, objects with low (high) spin have an above (below) average cK value; this behavior is inverted around ν* = 1.75, as marked by the vertical black line.
In the text

	Fig. 7 Linear shape bias as a function of linear density bias in the Planck simulation. Different colors represent the quartiles in the distribution of halo spins. One can clearly see a difference in the values of cK even at fixed b1, which indicates a nontrivial assembly-bias signal.
In the text

Fig. C.1 Comparison of measurements of shape-biases cK, cδK, and cK⊗K in two versions of the Narya simulation, with increasing resolution. For these measurements, shapes were computed using the FoF halo definition. One can see that the results generally match within the error bars, except for measurements of cK below M200,b = 1012M⊙h−1. The damping scale employed in measuring these bias parameters is kd = 0.1hMpc−1

In the text