Introduction
The shapes of dark matter halos and galaxies in the Universe are sensitive to the distribution of large-scale dark matter density fluctuations. The spatial correlations of shapes – their ‘intrinsic alignments’ – have been the subject of interest to cosmologists for decades [1 ] . They represent both a significant contaminant to weak lensing surveys [2 ] , as these correlations are inevitably present in any measurement of cosmic shear or ‘galaxy–galaxy lensing’, but they are also a cosmological signal of interest. As galaxy redshift and imaging surveys have become more powerful, large catalogs of galaxies with spectroscopic redshifts and high-quality shape measurements have been produced. 3D catalogs of galaxy shapes are promising probes of the BAO feature [3 ] , can help pin down the nature of dark energy [4 ] , and are potential probes of interesting new physics such as anisotropic non-gaussianities or parity violation [5 , 6 , 7 ] .
Beyond spectroscopic samples, Rubin Legacy Survey of Space and Time [8 ] , the Roman Space Telescope [9 , 10 , 11 ] and Euclid [12 , 13 ] space telescope will produce extremely dense samples of high-precision galaxy shapes. While these surveys will deliver shape catalogs of unprecedented quality, the degree to which intrinsic alignments can potentially contaminate weak lensing surveys is still in question. The depth to which modern lensing surveys will image galaxies (maximizing number density) far exceeds our spectroscopic follow-up capabilities, which are an important ingredient in data-driven alignment mitigation strategies. There have been many proposals to mitigate this effect, including selecting on subsets of galaxies which have smaller known alignment amplitudes (although only consistent with zero given measurement uncertainties) [14 ] , calibration and extrapolation of scaling relations from spectroscopic samples [15 ] , or alternative samples which can have their alignment amplitude directly calibrated from dedicated observing campaigns [16 ] .
So far, these efforts have predominantly focused on isolating and calibrating the linear alignment amplitude of shapes with the dark matter distribution – governed by its correlation with the large-scale tidal field. However, the symmetries of the Universe allow a much greater phenomenology of dependence of shapes, motivating ‘non-linear’ alignment models which include additional degrees of freedom such as contributions from tidal torquing, density-weighting, and so on. As imaging surveys become more powerful, our sensitivity to these additional alignment contributions becomes larger, and current direct calibration methods are not necessarily optimal to infer their signatures. This is because direct calibration relies on two-point functions, which are relatively insensitive to these contributions (they are ‘1-loop’ in the parlance of perturbation theory). From the perspective of analytic theory, many of the proposed models in the literature are not complete. EFT methods [17 , 18 , 19 , 20 , 21 ] are rigorous completions of the original crop of phenomenological models [1 , 22 ] . They have been tested rigorously for two point functions and, partially, for three point functions [23 , 24 , 25 , 26 ] . Notably, the EFT formulations of intrinsic alignments suffer from a large amount of ‘nuisance’ parameters – both due to the large number of degrees of freedom and the completely uncertain time evolution of the signal for galaxies. Employing maximally conservative priors on the amplitudes of these effects often significantly degrades constraining power on physics of interest. As we see it, the current status of intrinsic alignments presents some challenges: 1) direct calibration cannot extend to samples of galaxies representative of Roman, LSST, especially for their higher-order alignments and 2) it is unclear whether the full basis of operators in the EFT is needed to model galaxy alignments, and 3) to what degree do alignment models capture the correlation of shapes with the underlying dark matter density v.s. simply over fit to correlation functions.
The question of whether an EFT-based model of large-scale structure summary statistics adequately captures the full correlations or simply fits to said summaries is not unique to intrinsic alignments, it is a concern for all tracers of large-scale structure. To go beyond agreement at the level of summary statistics, there has been significant recent development in field-level perturbative models. The key principle of field-level modeling is that, by describing the entirety of the distribution of the tracer in question, all summary statistics are naturally well-described. Pioneering work in perturbative field-level models for the case of galaxy clustering was carried out in [27 ] , who developed a perturbative field-level forward model and applied it to study the distribution of halo number density fluctuations (see also [28 , 29 , 30 , 31 ] for perturbative field-level models that use a different construction). This construction has since been extended to redshift space [32 ] , galaxy densities in hydrodynamical simulations [33 , 34 ] , the HI distribution at low redshifts and during reionization [35 , 36 ] , and the Lyman-α \alpha α [37 ] . Beyond testing the validity of perturbative models in describing a tracer, by finding the set of EFT coefficients which minimize the field-level error, this approach is also a powerful estimator of said EFT coefficients. Measuring how EFT coefficients vary across different physical simulations of said tracers allows us to construct a mapping between model-dependent parameters in simulations and the general EFT space, understanding how different physics impacts EFT coefficients [38 , 39 , 40 , 41 ] . The measured values of coefficients can be used to calibrate priors for subsequent analyses in real data, enhancing the constraining power of EFT-based analyses [42 ] .
The focus of this paper is to develop the field-level formalism for the effective theory of intrinsic alignments, building on the recent development of a Lagrangian theory for IA [20 ] , whose power spectrum was calculated at one-loop order. We aim to test the ability of the EFT in describing the full distribution of correlations of dark matter halo and galaxy shapes, beyond linear theory, applying it to the large Quijote suite [43 ] of dark matter simulations as well as the recently publicly-released FLAMINGO suite of cosmological hydrodynamical simulations [44 , 45 ] . A visual summary of the field-level distribution of halo shapes, and best fits from the EFT model, is shown in Figure 1 . The use of this field-level model for alignments will also allow us to extract the values of the coefficients in the theory at high statistical significance even when limited by the volumes available, in complementary fashion to bispectrum-based measurements [23 ] , as well as “probabilistic estimators” [46 ] , which have been previously used in large-volume dark matter only simulations. The techniques, simulations and samples we utilize in this work are distinct from past investigations of intrinsic alignments in FLAMINGO, such as [24 , 47 , 48 , 49 ] . The past works have mostly focused on massive halos which host luminous red galaxies, as well as emphasizing the largest volume FLAMINGO box, at V = ( 2.8 G p c ) 3 V = (2.8\mathrm{Gpc})^3 V = ( 2.8 Gpc ) 3
Figure 1. A visual summary of the field-level shape distribution. Visualization of the halo density (top left) and halo shape helicity fields (top right) for halos in M h ∈ [ 10 12.8 , 10 13.8 ] M ⊙ h − 1 M_h \in[10^{12.8}, 10^{13.8}]M_\odot h^{-1} M h ∈ [ 1 0 12.8 , 1 0 13.8 ] M ⊙ h − 1 ≈ 100 M p c h − 1 \approx100\mathrm{Mpc}h^{-1} ≈ 100 Mpc h − 1 R = 7 M p c h − 1 R = 7\mathrm{Mpc}h^{-1} R = 7 Mpc h − 1 m = { 1 , 2 } m = \{1,2\} m = { 1 , 2 } k k k k = 1 / R k = 1/R k = 1/ R
This paper is structured as follows: in § II we review the theory of spin-2 biased tracers of the density field, constructing the set of 11 operators at cubic order, introducing the shifted basis, as well as the helicity decomposition we will use throughout this work to simplify measurements. We then discuss the field-level model of spin-2 tracers, and the associated transfer functions in § III. In § III C we discuss degeneracies in the spin-2 expansion that are different from those of scalar biases. Two key results from this section are the re-definition of the basis into the spin-2 analog of the Galileon basis often used in the scalar bias expansion, as well as uncovering a previously unknown redundancy that eliminates one operator, s a b 2 s^{2}_{ab} s ab 2
Technical details are relegated to a series of appendices: Appendix A discusses the construction of the orthogonalized transfer functions used in this work from first principles, Appendix B validates the spin-2 basis and helicity decomposition numerically against predictions from Gaussian field correlators, Appendix C discusses how we independently recreate corrupted data from the FLAMINGO simulation, Appendix D validates the analytic expressions for transfer function and parameter covariance matrices derived in this work, Appendix E shows an extended set of transfer functions for various halo masses that are not shown in the main work, and finally Appendix F derives the conversion between the shape bias coefficients used in this work and the alignment parameters customarily adopted in the weak lensing community.
Spin-2 theory of biased tracers – Recap and conventions
We begin by reviewing how we will define the distribution of shapes in this work, the correlations we will study, as well as the construction of the EFT model used in this work. Halo and galaxy shapes are spin-2 tracers of large-scale structure. In § IV we specify our exact scheme for how we construct a field of halo shapes from a catalog. For now, suppose a spin-2 field of shapes, I a b ( x ) I_{ab}(x) I ab ( x )
I a b ( x ) = δ a b K 3 δ I ( x ) + γ a b ( x ) . (1) I_{ab}(x)=\frac{\delta^{K}_{ab}}{3}\delta_I(x)+\gamma_{ab}(x). \tag*{(1)} I ab ( x ) = 3 δ ab K δ I ( x ) + γ ab ( x ) . (1) The trace of an object’s quadrupole tensor is related to its total size, and so δ I ( x ) \delta_I(x) δ I ( x ) size fluctuation field . Correspondingly, γ a b ( x ) \gamma_{ab}(x) γ ab ( x ) shape fluctuation field . It is γ a b \gamma_{ab} γ ab γ a b \gamma_{ab} γ ab γ a b \gamma_{ab} γ ab ∂ a ∂ b Φ \partial_a\partial_b\Phi ∂ a ∂ b Φ
∂ a ∂ b Φ = δ a b K 3 δ l i n + s a b , (2) \partial_a\partial_b\Phi=\frac{\delta^{K}_{ab}}{3}\delta_{\mathrm{lin}}+s_{ab}, \tag*{(2)} ∂ a ∂ b Φ = 3 δ ab K δ lin + s ab , (2) where δ l i n \delta_{\mathrm{lin}} δ lin s a b s_{ab} s ab
s a b = ( ∂ a ∂ b ∂ 2 − δ a b K 3 ) δ l i n . s_{ab}=\left(\frac{\partial_a\partial_b}{\partial^2}-\frac{\delta^{K}_{ab}}{3}\right)\delta_{\mathrm{lin}}. s ab = ( ∂ 2 ∂ a ∂ b − 3 δ ab K ) δ lin . Note that we have not specified the reference frame in which the expansion is constructed – either the late-time Eulerian frame with coordinates x x x q q q δ \delta δ [1 , 50 ] and ‘tidal alignment and tidal torquing’ (TATT) [2 ] , can be derived from counting powers of ∂ a ∂ b Φ \partial_a\partial_b\Phi ∂ a ∂ b Φ γ a b \gamma_{ab} γ ab ( a , b ) (a,b) ( a , b ) ∂ a ∂ b Φ \partial_a\partial_b\Phi ∂ a ∂ b Φ
γ a b ( x ) = c s s a b ( x ) . (3) \gamma_{ab}(x)=c_s s_{ab}(x). \tag*{(3)} γ ab ( x ) = c s s ab ( x ) . (3) where c s c_s c s linear alignment coefficient . It is often defined as ‘A 1 A_1 A 1 c s c_s c s A 1 A_1 A 1 s a b s_{ab} s ab linear tidal field gives the ‘linear alignment model’, but the traceless shear operator D a b \mathcal{D}_{ab} D ab nonlinear density field. This is equivalent to stating that the shapes of galaxies are immediately aligned and correlated with the present-day tidal field. The next tensor that one can construct from ∂ a ∂ b Φ \partial_a\partial_b\Phi ∂ a ∂ b Φ
∂ a ∂ c Φ ∂ c ∂ b Φ = δ a b δ 2 9 + 2 3 s a b δ + ( s 2 ) a b . \partial_a\partial_c\Phi\partial_c\partial_b\Phi= \frac{\delta_{ab}\delta^2}{9} + \frac{2}{3}s_{ab}\delta+ (s^2)_{ab}. ∂ a ∂ c Φ ∂ c ∂ b Φ = 9 δ ab δ 2 + 3 2 s ab δ + ( s 2 ) ab . which generates two new traceless tensor dependences – ( δ s ) a b (\delta s)_{ab} ( δ s ) ab ( s 2 ) a b (s^2)_{ab} ( s 2 ) ab [51 ] . The fields can also be generated from taking products of the trace and trace-free components of ∂ a ∂ b Φ \partial_a\partial_b\Phi ∂ a ∂ b Φ
γ a b ( x ) = c s s a b ( x ) + c δ s ( δ s ) a b ( x ) + c s 2 ( s 2 ) a b ( x ) . (4) \gamma_{ab}(x) = c_s s_{ab}(x) + c_{\delta s}(\delta s)_{ab}(x) + c_{s^2}(s^2)_{ab}(x). \tag*{(4)} γ ab ( x ) = c s s ab ( x ) + c δ s ( δ s ) ab ( x ) + c s 2 ( s 2 ) ab ( x ) . (4) There is one additional trace-free tensor that can be generated from the current fields at second order. Note that the following zero-mean, second-order, scalar field can be written
G 2 ( x ) = Tr [ ( s 2 ) a b ] ( x ) − 2 3 δ 2 ( x ) . (5) \mathcal{G}_2(x) = \operatorname{Tr}[(s^2)_{ab}](x) - \frac{2}{3}\delta^2(x). \tag*{(5)} G 2 ( x ) = Tr [( s 2 ) ab ] ( x ) − 3 2 δ 2 ( x ) . (5) This field in the context of scalar biasing theory is called the second-order Galileon field. Applying the tidal operator to this Galileon, D a b G 2 \mathcal{D}_{ab}\mathcal{G}_2 D ab G 2 D a b \mathcal{D}_{ab} D ab [17 ] .
Extending this construction to third order in ∂ a ∂ b Φ \partial_a\partial_b\Phi ∂ a ∂ b Φ
1 s t o r d e r : s i j 2 n d o r d e r : t i j , s i j 2 , δ s i j 3 r d o r d e r : L i j ( 3 ) , ( s t ) i j , δ t i j , s i j 3 , s 2 s i j , δ 2 s i j , δ s i j 2 . \begin{aligned}
1^{\mathrm{st}}\ \mathrm{order}:&\quad s_{ij}\\
2^{\mathrm{nd}}\ \mathrm{order}:&\quad t_{ij},\ s^2_{ij},\ \delta s_{ij}\\
3^{\mathrm{rd}}\ \mathrm{order}:&\quad L^{(3)}_{ij},\ (st)_{ij},\ \delta t_{ij},\ s^3_{ij},\ s^2s_{ij},\ \delta^2s_{ij},\ \delta s^2_{ij}.
\end{aligned} 1 st order : 2 nd order : 3 rd order : s ij t ij , s ij 2 , δ s ij L ij ( 3 ) , ( s t ) ij , δ t ij , s ij 3 , s 2 s ij , δ 2 s ij , δ s ij 2 . The Galileon term is normalized as t i j = 2 7 D i j G 2 = 2 7 D i j ( s 2 − 2 3 δ 2 ) t_{ij} = \frac{2}{7}\mathcal{D}_{ij}\mathcal{G}_2 = \frac{2}{7}\mathcal{D}_{ij}(s^2 - \frac{2}{3}\delta^2) t ij = 7 2 D ij G 2 = 7 2 D ij ( s 2 − 3 2 δ 2 ) velocity shear tensor in the literature [52 ] . It is equivalently defined by the trace-free part of the gradient of the second-order displacement computed in lagrangian perturbation theory (LPT) – t a b ∝ T F [ ∂ a Ψ b 2 ] t_{ab} \propto\mathrm{TF}[\partial_a\Psi^{2}_b] t ab ∝ TF [ ∂ a Ψ b 2 ] L i j ( 3 ) L^{(3)}_{ij} L ij ( 3 ) L a b ( 3 ) ∝ T F [ ∂ a Ψ b 3 ] L^{(3)}_{ab} \propto\mathrm{TF}[\partial_a\Psi^{3}_b] L ab ( 3 ) ∝ TF [ ∂ a Ψ b 3 ] ∇ 2 s a b \nabla^2s_{ab} ∇ 2 s ab q q q Ψ n \Psi_{\mathrm{n}} Ψ n
Before continuing onward let us make some remarks about this basis. When studying the IA power spectrum at one-loop order, the cubic basis can be further simplified. Operators which are only built from powers of s i j s_{ij} s ij δ \delta δ 1 L i j ( 3 ) L^{(3)}_{ij} L ij ( 3 ) ( s t ) i j (st)_{ij} ( s t ) ij δ t i j \delta t_{ij} δ t ij [17 , 20 ] due to a constraint that appears at the tadpole momentum configuration of the one-loop diagram. However, all three are non-trivial at the field level. Another remark we make, but which we will expand on in § III C, is that there are in fact only 10 independent operators in the cubic shape bias basis, instead of the usual 11 that are written down. A constraint that was written down in the one-loop context in Ref. [17 ] extends to the entire field level, which leads to finding the s 2 s i j s^2s_{ij} s 2 s ij
Having enumerated the standard bias shape basis, we write the complete perturbative expansion for the shape fluctuation field in terms of the allowed operators2
γ a b = ∑ O ∈ { s , δ s , s 2 , … } c O O a b + ϵ a b (6) \gamma_{ab} = \sum_{\mathcal{O}\in\{s,\delta s,s^2,\ldots\}} c_{\mathcal{O}}\mathcal{O}_{ab} + \epsilon_{ab} \tag*{(6)} γ ab = O ∈ { s , δ s , s 2 , … } ∑ c O O ab + ϵ ab (6) where all 10 operators receive a coefficient c O c_O c O ϵ a b \epsilon_{ab} ϵ ab shape stochasticity field , which by definition is uncorrelated with any of the deterministic contributions. The symmetries of the EFT constrain the correlations functions of ϵ a b \epsilon_{ab} ϵ ab k → 0 k \to0 k → 0
⟨ ϵ a b ( k ) ϵ c d ( k ′ ) ⟩ = ( 2 π ) 3 δ D ( k + k ′ ) ( δ a c K δ b d K + δ a d K δ b c K − 2 3 δ a b K δ c d K ) P ϵ , (7) \left\langle\epsilon_{ab}(\boldsymbol{k})\epsilon_{cd}(\boldsymbol{k}') \right\rangle= (2\pi)^3\delta_{\mathrm{D}}(\boldsymbol{k}+\boldsymbol{k}')\left(\delta_{ac}^{\mathrm{K}}\delta_{bd}^{\mathrm{K}}+\delta_{ad}^{\mathrm{K}}\delta_{bc}^{\mathrm{K}}-\frac{2}{3}\delta_{ab}^{\mathrm{K}}\delta_{cd}^{\mathrm{K}}\right)P_{\epsilon}, \tag*{(7)} ⟨ ϵ ab ( k ) ϵ c d ( k ′ ) ⟩ = ( 2 π ) 3 δ D ( k + k ′ ) ( δ a c K δ b d K + δ a d K δ b c K − 3 2 δ ab K δ c d K ) P ϵ , (7) where P ϵ P_{\epsilon} P ϵ P ϵ = σ γ 2 / n ˉ P_{\epsilon}=\sigma_{\gamma}^{2}/\bar{n} P ϵ = σ γ 2 / n ˉ σ γ 2 \sigma_{\gamma}^{2} σ γ 2 n ˉ \bar{n} n ˉ
Equation (6) can be defined either in late-time, Eulerian coordinates x \boldsymbol{x} x q \boldsymbol{q} q Ψ ( q ) \boldsymbol{\Psi}(\boldsymbol{q}) Ψ ( q )
γ a b ( x ) = ∫ d 3 q δ D ( x − q − Ψ ( q ) ) R a c ( q ) R b d ( q ) γ c d L ( q ) , (8) \gamma_{ab}(\boldsymbol{x})=\int\mathrm{d}^{3}q\,\delta_{\mathrm{D}}(\boldsymbol{x}-\boldsymbol{q}-\boldsymbol{\Psi}(\boldsymbol{q}))R_{ac}(\boldsymbol{q})R_{bd}(\boldsymbol{q})\gamma_{cd}^{L}(\boldsymbol{q}), \tag*{(8)} γ ab ( x ) = ∫ d 3 q δ D ( x − q − Ψ ( q )) R a c ( q ) R b d ( q ) γ c d L ( q ) , (8) where R a c ≡ δ a c + ∂ a Ψ c R_{ac}\equiv\delta_{ac}+\partial_{a}\Psi_{c} R a c ≡ δ a c + ∂ a Ψ c γ L \gamma^{L} γ L proto-halos [20 ] . A common assumption that is made in alignment theories is to assume “passive evolution” models of halo shapes. In the passive alignment picture, the distortion matrix R a c R_{ac} R a c [55 ] . When passive alignment is assumed, predictions can be made between the coefficients of the Lagrangian bias expansion and the Eulerian bias expansion, similar to the ‘co-evolution’ relations that have been derived for the bias of galaxy density fluctuations [56 ] .
The “active" v.s. “passive" distinction does not apply to the basis fields of the bias expansion themselves, which are defined in the Eulerian or Lagrangian frames and connected by a displacement vector. In this paper, we adopt the “shifted operator" basis used to study field-level perturbative EFT, first defined in [27 ] . Shifted operators are defined by treating the Zeldovich displacement non-perturbatively, but treating higher-order contributions in perturbation theory. In Fourier space a shifted shape operator is defined as
O ~ a b ( k ) = ∫ d 3 q e − i k ⋅ q e − i k ⋅ Ψ 1 ( q ) O a b ( q ) , (9) \widetilde{\mathcal{O}}_{ab}(\boldsymbol{k})=\int\mathrm{d}^{3}q e^{-i\boldsymbol{k}\cdot\boldsymbol{q}}e^{-i\boldsymbol{k}\cdot\boldsymbol{\Psi}_{1}(\boldsymbol{q})}\mathcal{O}_{ab}(\boldsymbol{q}), \tag*{(9)} O ab ( k ) = ∫ d 3 q e − i k ⋅ q e − i k ⋅ Ψ 1 ( q ) O ab ( q ) , (9) where a tilde applied to an operator denotes it is in the shifted basis.
We implement the tensor basis at the field level as follows: given the linear initial conditions used to seed a simulation, δ I C ( q ) \delta_{\mathrm{IC}}(\boldsymbol{q}) δ IC ( q ) s a b ( q ) s_{ab}(\boldsymbol{q}) s ab ( q ) D a b D_{ab} D ab δ \delta δ s a b s_{ab} s ab D a b D_{ab} D ab s a b s_{ab} s ab
s a b 3 → s a b 3 − 7 15 σ 2 s a b , δ 2 s a b → δ 2 s a b − σ 2 s a b , δ s a b 2 → δ s a b 2 s t a b → s t a b − 8 315 σ 2 s a b , δ t a b → δ t a b − 8 105 σ 2 s a b , (10) \begin{aligned}
s_{ab}^{3} &\to s_{ab}^{3}-\frac{7}{15}\sigma^{2}s_{ab}, \qquad\delta^{2}s_{ab}\to\delta^{2}s_{ab}-\sigma^{2}s_{ab}, \qquad\delta s_{ab}^{2}\to\delta s_{ab}^{2} \\
st_{ab} &\to st_{ab}-\frac{8}{315}\sigma^{2}s_{ab}, \qquad\delta t_{ab}\to\delta t_{ab}-\frac{8}{105}\sigma^{2}s_{ab}, \tag*{(10)}
\end{aligned} s ab 3 s t ab → s ab 3 − 15 7 σ 2 s ab , δ 2 s ab → δ 2 s ab − σ 2 s ab , δ s ab 2 → δ s ab 2 → s t ab − 315 8 σ 2 s ab , δ t ab → δ t ab − 105 8 σ 2 s ab , (10) and L a b ( 3 ) L_{ab}^{(3)} L ab ( 3 ) L = 1000 h − 1 M p c L=1000\,h^{-1}\mathrm{Mpc} L = 1000 h − 1 Mpc N = 1024 N=1024 N = 1024 N = 512 N=512 N = 512 N = 256 N=256 N = 256 k N y q ≈ 0.8 h M p c − 1 k_{\mathrm{Nyq}}\approx0.8\,h\mathrm{Mpc}^{-1} k Nyq ≈ 0.8 h Mpc − 1 [57 ] , except with perturbative displacements in lieu of the full N N N
The two-point correlations of shape operators are in principle rank 4 tensors. In Ref. [17 ] it was shown that tensor correlators simplify significantly if the fields are initially decomposed into the helicity basis of tensor spherical harmonics
O a b ( k ) = ∑ m = − 2 2 Y a b m ( k ^ ) O m ( k ) . (11) \mathcal{O}_{ab}(\boldsymbol{k})=\sum_{m=-2}^{2}Y_{ab}^{m}(\hat{\boldsymbol{k}})\mathcal{O}^{m}(\boldsymbol{k}). \tag*{(11)} O ab ( k ) = m = − 2 ∑ 2 Y ab m ( k ^ ) O m ( k ) . (11) Each m m m
⟨ O m ( k ) O m ′ ′ ( k ′ ) ⟩ = ( 2 π ) 3 δ m m ′ K δ D ( k + k ′ ) P O O ′ m ( k ) . (12) \left\langle\mathcal{O}_m(\mathbf{k})\mathcal{O}'_{m'}(\mathbf{k}') \right\rangle= (2\pi)^3\delta^K_{mm'}\delta^D(\mathbf{k}+\mathbf{k}')P^{m}_{\mathcal{O}\mathcal{O}'}(k). \tag*{(12)} ⟨ O m ( k ) O m ′ ′ ( k ′ ) ⟩ = ( 2 π ) 3 δ m m ′ K δ D ( k + k ′ ) P O O ′ m ( k ) . (12) In the helicity basis, shape power spectra are simplified to five scalar spectra. In the absence of parity violation, the positive and negative m m m [7 ] , and only three spectra are required to describe the two-point statistics of the shape field. The spherical tensors are defined in Fourier space in the following way: For every k \mathbf{k} k n ^ \hat{\mathbf{n}} n ^ { k ^ , e ^ 1 , e ^ 2 } \{\hat{\mathbf{k}},\hat{\mathbf{e}}_1,\hat{\mathbf{e}}_2\} { k ^ , e ^ 1 , e ^ 2 }
e ^ 1 = k × n ^ ∣ k × n ^ ∣ , e ^ 2 = k ^ × e ^ 1 . (13) \hat{\mathbf{e}}_1 = \frac{\mathbf{k}\times\hat{\mathbf{n}}}{\left|\mathbf{k}\times\hat{\mathbf{n}}\right|}, \qquad\hat{\mathbf{e}}_2 = \hat{\mathbf{k}}\times\hat{\mathbf{e}}_1. \tag*{(13)} e ^ 1 = ∣ k × n ^ ∣ k × n ^ , e ^ 2 = k ^ × e ^ 1 . (13) From this orthonormal basis we construct the helicity basis, composed of raising and lowering vectors { k ^ , e ^ + , e ^ − } \{\hat{\mathbf{k}},\hat{\mathbf{e}}_+,\hat{\mathbf{e}}_-\} { k ^ , e ^ + , e ^ − } e ^ ± = 1 2 ( e ^ 1 ± i e ^ 2 ) \hat{\mathbf{e}}_\pm= \frac{1}{\sqrt{2}}(\hat{\mathbf{e}}_1 \pm i\hat{\mathbf{e}}_2) e ^ ± = 2 1 ( e ^ 1 ± i e ^ 2 ) k ^ \hat{\mathbf{k}} k ^ e ^ ± \hat{\mathbf{e}}_\pm e ^ ±
Y ℓ = 2 , a b m = 0 ( k ) = 3 2 ( k ^ a k ^ b − 1 3 δ a b ) , Y ℓ = 2 , a b m = ± 1 ( k ) = 1 2 ( k ^ a e ^ b ± + k ^ b e ^ a ± ) , Y ℓ = 2 , a b m = ± 2 ( k ) = e ^ a ± e ^ b ± . (14) Y^{m=0}_{\ell=2,ab}(\mathbf{k}) = \sqrt{\frac{3}{2}}\left(\hat{k}_a\hat{k}_b-\frac{1}{3}\delta_{ab}\right), \qquad Y^{m=\pm1}_{\ell=2,ab}(\mathbf{k}) = \sqrt{\frac{1}{2}}\left(\hat{k}_a\hat{e}^{\pm}_b+\hat{k}_b\hat{e}^{\pm}_a\right), \qquad Y^{m=\pm2}_{\ell=2,ab}(\mathbf{k}) = \hat{e}^{\pm}_a\hat{e}^{\pm}_b. \tag*{(14)} Y ℓ = 2 , ab m = 0 ( k ) = 2 3 ( k ^ a k ^ b − 3 1 δ ab ) , Y ℓ = 2 , ab m = ± 1 ( k ) = 2 1 ( k ^ a e ^ b ± + k ^ b e ^ a ± ) , Y ℓ = 2 , ab m = ± 2 ( k ) = e ^ a ± e ^ b ± . (14) To construct these basis tensors on the lattice we pick a random vector n \mathbf{n} n k \mathbf{k} k 3 (13) and (14) using the resulting array of basis vectors. Projecting the shifted operators with the Y ℓ = 2 , a b m Y^{m}_{\ell=2,ab} Y ℓ = 2 , ab m
We show validation of our Lagrangian basis implementation, and the field-level helicity decomposition, against analytic theory predictions in Appendix B. Note: Henceforth, every index pair a b ab ab c d cd c d
The transfer functions of shapes
A catalog of tracer shapes and their positions can be used to generate the shape fluctuation field γ a b ( k ) \gamma_{ab}(\mathbf{k}) γ ab ( k ) (6) ) into the helicity basis by projecting the tensor expression with Y a b m ∗ ( k ) Y^{m*}_{ab}(\mathbf{k}) Y ab m ∗ ( k )
γ m ( k ) = ∑ O c O O ~ m ( k ) + ϵ m ( k ) . (15) \gamma_m(\mathbf{k}) = \sum_{\mathcal{O}} c_{\mathcal{O}}\widetilde{\mathcal{O}}_m(\mathbf{k})+\epsilon_m(\mathbf{k}). \tag*{(15)} γ m ( k ) = O ∑ c O O m ( k ) + ϵ m ( k ) . (15) For any set of coefficients c O c_{\mathcal{O}} c O ϵ m ( k ) \epsilon_m(\mathbf{k}) ϵ m ( k ) c O c_{\mathcal{O}} c O ∣ γ m ( k ) − ∑ O c O O ~ m ( k ) ∣ \left|\gamma_m(\mathbf{k})-\sum_{\mathcal{O}}c_{\mathcal{O}}\widetilde{\mathcal{O}}_m(\mathbf{k})\right| γ m ( k ) − ∑ O c O O m ( k ) power spectrum
⟨ ϵ m ( k ) ϵ m ′ ( k ′ ) ⟩ = ( 2 π ) 3 δ D ( k + k ′ ) δ m m ′ K P e r r m ( k ) . (16) \left\langle\epsilon_m(\mathbf{k})\epsilon_{m'}(\mathbf{k}') \right\rangle= (2\pi)^3\delta^D(\mathbf{k}+\mathbf{k}')\delta^K_{mm'}P^m_{\mathrm{err}}(k). \tag*{(16)} ⟨ ϵ m ( k ) ϵ m ′ ( k ′ ) ⟩ = ( 2 π ) 3 δ D ( k + k ′ ) δ m m ′ K P err m ( k ) . (16) There are, in total, three error power spectra P e r r m ( k ) P^m_{\mathrm{err}}(k) P err m ( k ) m m m ± m \pm m ± m P e r r ( k ) P_{\mathrm{err}}(k) P err ( k )
γ m ( k ) = c s s ~ m ( k ) + ϵ m ( k ) . (17) \gamma_m(\mathbf{k}) = c_s\widetilde{s}_m(\mathbf{k})+\epsilon_m(\mathbf{k}). \tag*{(17)} γ m ( k ) = c s s m ( k ) + ϵ m ( k ) . (17) In order to achieve the best possible field-level fit to the measured shape field, we can promote the coefficient c s c_s c s transfer function c s ; m ( k ) c_{s;m}(k) c s ; m ( k )
P e r r m ( k ) = ⟨ ∣ γ m ( k ) − c s ; m ( k ) s ~ m ( k ) ∣ 2 ⟩ . (18) P^m_{\mathrm{err}}(k)=\left\langle\left|\gamma_m(\mathbf{k})-c_{s;m}(k)\widetilde{s}_m(\mathbf{k})\right|^2\right\rangle. \tag*{(18)} P err m ( k ) = ⟨ ∣ γ m ( k ) − c s ; m ( k ) s m ( k ) ∣ 2 ⟩ . (18) Minimizing this error power spectrum for each k k k
c ^ s ; m ( k ) = ⟨ s ~ m ( k ) γ ~ m ∗ ( k ) ⟩ ⟨ s ~ m ( k ) s ~ m ∗ ( k ) ⟩ . (19) \hat{c}_{s;m}(k)=\frac{\langle\tilde{s}_m(k)\tilde{\gamma}^{*}_m(k)\rangle}{\langle\tilde{s}_m(k)\tilde{s}^{*}_m(k)\rangle}. \tag*{(19)} c ^ s ; m ( k ) = ⟨ s ~ m ( k ) s ~ m ∗ ( k )⟩ ⟨ s ~ m ( k ) γ ~ m ∗ ( k )⟩ . (19) Several comments are in order. First, the transfer functions of spin-2 tracers are m m m m m m ⟨ s ~ m ( k ) s ~ m ∗ ( k ) ⟩ \langle\tilde{s}_m(k)\tilde{s}^{*}_m(k)\rangle ⟨ s ~ m ( k ) s ~ m ∗ ( k )⟩ m = 1 m=1 m = 1 m = 2 m=2 m = 2 [20 ] .
The second important observation is that the error field defined as ϵ m = γ m − γ m m o d e l \epsilon_m=\gamma_m-\gamma^{\mathrm{model}}_m ϵ m = γ m − γ m model γ m m o d e l \gamma^{\mathrm{model}}_m γ m model (7) , which predicts a uniform m m m P e r r = σ γ 2 / n ˉ P_{\mathrm{err}}=\sigma_\gamma^2/\bar{n} P err = σ γ 2 / n ˉ P e r r P_{\mathrm{err}} P err k → 0 k\to0 k → 0 [58 ] . The former contribution can be removed by explicitly including the higher-order bias operators in the field-level shape model, which we discuss now.
This extension to higher-order operators is implemented by performing a helicity-dependent orthogonalization of the operators,
γ m ( k ) = ∑ O C O ; m ( k ) O ~ m ⊥ + ϵ m ( k ) , (20) \gamma_m(k)=\sum_O C_{O;m}(k)\widetilde{O}^{\perp}_m+\epsilon_m(k), \tag*{(20)} γ m ( k ) = O ∑ C O ; m ( k ) O m ⊥ + ϵ m ( k ) , (20) such that every new operator that is included is by construction uncorrelated with previous members of the basis. An advantage of performing this orthogonalization in the helicity basis is that the process becomes fully analogous to the standard approach of scalar fields from Ref. [27 ] , but with five copies for each operator. The field-level decomposition of galaxy shapes is thus described by three transfer functions per operator,
C ^ O ; m ( k ) = ⟨ O ~ m ⊥ γ m ⟩ ( k ) ⟨ ∣ O ~ m ⊥ ∣ 2 ⟩ ( k ) , (21) \hat{C}_{O;m}(k)=\frac{\langle\widetilde{O}^{\perp}_m\gamma_m\rangle(k)}{\langle|\widetilde{O}^{\perp}_m|^2\rangle(k)}, \tag*{(21)} C ^ O ; m ( k ) = ⟨ ∣ O m ⊥ ∣ 2 ⟩ ( k ) ⟨ O m ⊥ γ m ⟩ ( k ) , (21) after averaging over ± m \pm m ± m
Transfer functions of shifted shape fields
We can glean some more intuition for the shape transfer functions, and their m m m
We begin by recalling how perturbation theory can be used to understand corrections to the field-level galaxy bias expansion. Ref. [27 ] expresses the halo density field in shifted operator form as
δ h ( k ) = ∫ q e − i k ⋅ Ψ 1 [ 1 + δ h L + i k a ( 1 + δ h L ) ∂ b Ψ 1 a Ψ 1 b − ( 1 + δ h L ) ∇ ⋅ Ψ ˉ − Ψ ˉ ⋅ ∇ δ h L ] e − i k ⋅ q , (22) \delta_h(\mathbf{k})=\int_{\mathbf{q}}e^{-i\mathbf{k}\cdot\mathbf{\Psi}_1}\left[1+\delta_h^L+i k_a(1+\delta_h^L)\partial_b\Psi_1^a\Psi_1^b-(1+\delta_h^L)\nabla\cdot\bar{\mathbf{\Psi}}-\bar{\mathbf{\Psi}}\cdot\nabla\delta_h^L\right]e^{-i\mathbf{k}\cdot\mathbf{q}}, \tag*{(22)} δ h ( k ) = ∫ q e − i k ⋅ Ψ 1 [ 1 + δ h L + i k a ( 1 + δ h L ) ∂ b Ψ 1 a Ψ 1 b − ( 1 + δ h L ) ∇ ⋅ Ψ ˉ − Ψ ˉ ⋅ ∇ δ h L ] e − i k ⋅ q , (22) where Ψ ˉ = Ψ 2 + Ψ 3 + ⋯ \bar{\mathbf{\Psi}}=\mathbf{\Psi}_2+\mathbf{\Psi}_3+\cdots Ψ ˉ = Ψ 2 + Ψ 3 + ⋯ ( 1 + δ h ( x ) ) d 3 x = ( 1 + δ h L ( q ) ) d 3 q (1+\delta_h(\mathbf{x}))d^3x=(1+\delta_h^L(\mathbf{q}))d^3q ( 1 + δ h ( x )) d 3 x = ( 1 + δ h L ( q )) d 3 q γ ( x ) \gamma(\mathbf{x}) γ ( x )
γ a b ( x ) = ∫ d 3 q δ D ( x − q − Ψ ( q ) ) γ a b L ( q ) , (23) \gamma_{ab}(\mathbf{x})=\int d^3q\,\delta_D(\mathbf{x}-\mathbf{q}-\mathbf{\Psi}(\mathbf{q}))\gamma^L_{ab}(\mathbf{q}), \tag*{(23)} γ ab ( x ) = ∫ d 3 q δ D ( x − q − Ψ ( q )) γ ab L ( q ) , (23) where γ a b L ( q ) \gamma^L_{ab}(\mathbf{q}) γ ab L ( q ) c O c_O c O [55 ] , and which is equivalently formulated as a conservation ansatz
γ a b ( x ) d 3 x = γ a b L ( q ) d 3 q . (24) \gamma_{ab}(\mathbf{x})d^3x=\gamma^L_{ab}(\mathbf{q})d^3q. \tag*{(24)} γ ab ( x ) d 3 x = γ ab L ( q ) d 3 q . (24) From the above equation, ‘co-evolution relations’ for the coefficients can be derived [55 , 59 ] , but it is not clear to which degree they should hold in practice, as they neglect the contribution from active R a b ( q ) R_{ab}(\mathbf{q}) R ab ( q ) s ~ \tilde{s} s ~ (17) , neglecting the stochasticity term for now. After shifting, we have
γ ~ a b ( k ) = c s s ~ a b ( k ) = c s ∫ d 3 q e − i k ⋅ q e − i k ⋅ Ψ 1 ( q ) s a b ( q ) . (25) \tilde{\gamma}_{ab}(\boldsymbol{k})=c_s\tilde{s}_{ab}(\boldsymbol{k})=c_s\int\mathrm{d}^3q e^{-i\boldsymbol{k}\cdot\boldsymbol{q}}e^{-i\boldsymbol{k}\cdot\boldsymbol{\Psi}_1(\boldsymbol{q})}s_{ab}(\boldsymbol{q}). \tag*{(25)} γ ~ ab ( k ) = c s s ~ ab ( k ) = c s ∫ d 3 q e − i k ⋅ q e − i k ⋅ Ψ 1 ( q ) s ab ( q ) . (25) Expanding the Ψ 1 ( q ) \boldsymbol{\Psi}_1(\boldsymbol{q}) Ψ 1 ( q )
γ ~ a b ( k ) = c s [ s a b − Ψ 1 ⋅ ∇ s a b + δ s a b ] ( k ) . (26) \tilde{\gamma}_{ab}(\boldsymbol{k})=c_s\left[s_{ab}-\boldsymbol{\Psi}_1\cdot\nabla s_{ab}+\delta s_{ab}\right](\boldsymbol{k}). \tag*{(26)} γ ~ ab ( k ) = c s [ s ab − Ψ 1 ⋅ ∇ s ab + δ s ab ] ( k ) . (26) Thus, the linearly biased shifted tracer picks up a density-weighted contribution and therefore a m ≠ 0 m\ne0 m = 0 4
Instead of expanding in the Ψ 1 \boldsymbol{\Psi}_1 Ψ 1 (24) allows for an analogous expression to Equation (22) to be derived for galaxy shapes. It is given by (now i , j i,j i , j a , b a,b a , b [20 ] )
γ a b ( k ) = ∫ q e − i k ⋅ Ψ 1 [ γ ˉ a b L ( q ) + i k i γ ˉ a b L ∂ j Ψ 1 i Ψ ˉ j − γ ˉ a b L ∇ ⋅ Ψ ˉ − Ψ ˉ ⋅ ∇ γ ˉ a b L ] e − i k ⋅ q . (27) \gamma_{ab}(\boldsymbol{k})=\int_{\boldsymbol{q}}e^{-i\boldsymbol{k}\cdot\boldsymbol{\Psi}_1}\left[\bar{\gamma}^{L}_{ab}(\boldsymbol{q})+ik_i\bar{\gamma}^{L}_{ab}\partial_j\Psi_1^i\bar{\Psi}^j-\bar{\gamma}^{L}_{ab}\nabla\cdot\bar{\boldsymbol{\Psi}}-\bar{\boldsymbol{\Psi}}\cdot\nabla\bar{\gamma}^{L}_{ab}\right]e^{-i\boldsymbol{k}\cdot\boldsymbol{q}}. \tag*{(27)} γ ab ( k ) = ∫ q e − i k ⋅ Ψ 1 [ γ ˉ ab L ( q ) + i k i γ ˉ ab L ∂ j Ψ 1 i Ψ ˉ j − γ ˉ ab L ∇ ⋅ Ψ ˉ − Ψ ˉ ⋅ ∇ γ ˉ ab L ] e − i k ⋅ q . (27) The lack of a ‘1+’ implies that the second term in square brackets will always produce an operator larger than one-loop. The third and fourth terms will only get contributions from Ψ 2 \boldsymbol{\Psi}_2 Ψ 2 ∇ ⋅ Ψ 3 \nabla\cdot\boldsymbol{\Psi}_3 ∇ ⋅ Ψ 3
γ ˉ a b L ∇ ⋅ Ψ 2 → 3 14 c s L G 2 s a b , \bar{\gamma}^{L}_{ab}\nabla\cdot\boldsymbol{\Psi}_2\to\frac{3}{14}c_s^L\mathcal{G}_2s_{ab}, γ ˉ ab L ∇ ⋅ Ψ 2 → 14 3 c s L G 2 s ab , Ψ ˉ ⋅ ∇ γ ˉ a b L → c s L Ψ 2 ⋅ ∇ s a b . \bar{\boldsymbol{\Psi}}\cdot\nabla\bar{\gamma}^{L}_{ab}\to c_s^L\boldsymbol{\Psi}_2\cdot\nabla s_{ab}. Ψ ˉ ⋅ ∇ γ ˉ ab L → c s L Ψ 2 ⋅ ∇ s ab . Both of these terms contribute at one-loop order with the same-amplitude for its coefficient, but note that G 2 s a b \mathcal{G}_2s_{ab} G 2 s ab δ 2 s a b \delta^2s_{ab} δ 2 s ab s a b 3 s_{ab}^3 s ab 3 Ψ 2 ⋅ ∇ s a b \boldsymbol{\Psi}_2\cdot\nabla s_{ab} Ψ 2 ⋅ ∇ s ab c s L c_s^L c s L
From the perturbative decomposition of the shape field in the shifted basis, we compute transfer functions for shape operators. We first project the shape tensor into the helicity basis
γ m ( k ) = Y m a b , ∗ ( k ^ ) γ a b ( k ) , (28) \gamma_m(\boldsymbol{k})=Y_m^{ab,*}(\hat{\boldsymbol{k}})\gamma_{ab}(\boldsymbol{k}), \tag*{(28)} γ m ( k ) = Y m ab , ∗ ( k ^ ) γ ab ( k ) , (28) from which we get transfer functions for every m m m (19)
C s ; m ( k ) = ⟨ s ~ m ∑ O c O O ~ m ⟩ ⟨ s ~ m s ~ m ⟩ = c s L + α m k 2 + c s L ⟨ s ~ m 2 ∣ s ~ m ⟩ ⟨ s ~ m s ~ m ⟩ + c δ s L ⟨ [ δ s ~ ] m ∣ s ~ m ⟩ ⟨ s ~ m s ~ m ⟩ + c t L ⟨ t ~ m ∣ s ~ m ⟩ ⟨ s ~ m s ~ m ⟩ + c s L ⟨ [ S ~ 3 m ] ∣ s ~ m ⟩ ⟨ s ~ m s ~ m ⟩ . (29) \begin{aligned}
C_{s;m}(k)=\frac{\langle\tilde{s}_m\sum_O c_O\tilde{O}_m\rangle}{\langle\tilde{s}_m\tilde{s}_m\rangle} \\
&=\ c_s^L+\alpha_m k^2+c_s^L\frac{\langle\tilde{s}_m^2|\tilde{s}_m\rangle}{\langle\tilde{s}_m\tilde{s}_m\rangle}+c_{\delta s}^L\frac{\langle[\widetilde{\delta s}]_m|\tilde{s}_m\rangle}{\langle\tilde{s}_m\tilde{s}_m\rangle}+c_t^L\frac{\langle\tilde{t}_m|\tilde{s}_m\rangle}{\langle\tilde{s}_m\tilde{s}_m\rangle}+c_s^L\frac{\langle[\tilde{\mathcal{S}}_{3m}]|\tilde{s}_m\rangle}{\langle\tilde{s}_m\tilde{s}_m\rangle}. \tag*{(29)}
\end{aligned} C s ; m ( k ) = ⟨ s ~ m s ~ m ⟩ ⟨ s ~ m ∑ O c O O ~ m ⟩ = c s L + α m k 2 + c s L ⟨ s ~ m s ~ m ⟩ ⟨ s ~ m 2 ∣ s ~ m ⟩ + c δ s L ⟨ s ~ m s ~ m ⟩ ⟨[ δ s ] m ∣ s ~ m ⟩ + c t L ⟨ s ~ m s ~ m ⟩ ⟨ t ~ m ∣ s ~ m ⟩ + c s L ⟨ s ~ m s ~ m ⟩ ⟨[ S ~ 3 m ] ∣ s ~ m ⟩ . (29) where S ~ 3 , a b ≡ s ~ a b 2 L P T − s ~ a b Z A ≈ Ψ 2 ⋅ ∇ s a b + 3 14 G 2 s a b \tilde{\mathcal{S}}_{3,ab}\equiv\tilde{s}^{2\mathrm{LPT}}_{ab}-\tilde{s}^{\mathrm{ZA}}_{ab}\approx\boldsymbol{\Psi}_2\cdot\nabla s_{ab}+\frac{3}{14}\mathcal{G}_2s_{ab} S ~ 3 , ab ≡ s ~ ab 2 LPT − s ~ ab ZA ≈ Ψ 2 ⋅ ∇ s ab + 14 3 G 2 s ab α m k 2 \alpha_m k^2 α m k 2 k k k S ~ 3 \tilde{\mathcal{S}}_3 S ~ 3 ⟨ S ~ 3 s ⟩ = 0 \langle\tilde{\mathcal{S}}_3s\rangle=0 ⟨ S ~ 3 s ⟩ = 0 5
Example of the helicity-dependence of transfer functions
As discussed, the shifted Lagrangian tidal tensor generates a density-weighted contribution to the final shape anisotropy field. This produces a helicity-dependence of the transfer functions. In addition, the m = 2 m=2 m = 2 k → 0 k \to0 k → 0
C s ; m ( k ) ⊃ ⟨ [ δ s ~ ] m ∣ s ~ m ⟩ ( k ) ⟨ s ~ m s ~ m ⟩ ( k ) . (30) C_{s;m}(k) \supset\frac{\langle[\widetilde{\delta s}]_m|\widetilde{s}_m\rangle(k)}{\langle\widetilde{s}_m\widetilde{s}_m\rangle(k)}. \tag*{(30)} C s ; m ( k ) ⊃ ⟨ s m s m ⟩ ( k ) ⟨[ δ s ] m ∣ s m ⟩ ( k ) . (30) On general scaling grounds we expect this transfer function will not asymptote to a constant at low-k k k ( δ s ) m (\delta s)_m ( δ s ) m ϵ i ± \epsilon_i^\pm ϵ i ± p \boldsymbol{p} p α ≡ ϵ a ± p ^ a \alpha\equiv\epsilon_a^\pm\hat{p}_a α ≡ ϵ a ± p ^ a
( δ s ) m = 2 ( k ) ⊃ ∫ p α 2 δ p δ k − p . (\delta s)_{m=2}(\boldsymbol{k}) \supset\int_{\boldsymbol{p}} \alpha^2\delta_{\boldsymbol{p}}\delta_{\boldsymbol{k}-\boldsymbol{p}}. ( δ s ) m = 2 ( k ) ⊃ ∫ p α 2 δ p δ k − p . For the linearized shift, the higher-helicity contributions necessarily come from
s ~ a b ( k ) ⊃ i k ⋅ [ Ψ s a b ] ( k ) \widetilde{s}_{ab}(\boldsymbol{k}) \supset i\boldsymbol{k}\cdot[\boldsymbol{\Psi}s_{ab}](\boldsymbol{k}) s ab ( k ) ⊃ i k ⋅ [ Ψ s ab ] ( k ) ⊃ k ⋅ [ ∫ p p p 2 S a b ( k − p ) + k − p ( k − p ) 2 S a b ( p ^ ) ] δ k − p δ p . \supset\boldsymbol{k}\cdot\left[\int_{\boldsymbol{p}}\frac{\boldsymbol{p}}{p^2}S_{ab}(\boldsymbol{k}-\boldsymbol{p})+\frac{\boldsymbol{k}-\boldsymbol{p}}{(\boldsymbol{k}-\boldsymbol{p})^2}S_{ab}(\hat{\boldsymbol{p}})\right]\delta_{\boldsymbol{k}-\boldsymbol{p}}\delta_{\boldsymbol{p}}. ⊃ k ⋅ [ ∫ p p 2 p S ab ( k − p ) + ( k − p ) 2 k − p S ab ( p ^ ) ] δ k − p δ p . Applying Y a b m = 2 ( k ) Y^{m=2}_{ab}(\boldsymbol{k}) Y ab m = 2 ( k )
s ~ m = 2 ( k ) ∝ ∫ p α 2 k 2 ( k − p ) 2 δ k − p δ p , (31) \widetilde{s}_{m=2}(\boldsymbol{k}) \propto\int_{\boldsymbol{p}}\frac{\alpha^2 k^2}{(\boldsymbol{k}-\boldsymbol{p})^2}\delta_{\boldsymbol{k}-\boldsymbol{p}}\delta_{\boldsymbol{p}}, \tag*{(31)} s m = 2 ( k ) ∝ ∫ p ( k − p ) 2 α 2 k 2 δ k − p δ p , (31) In the low-k k k s ~ m = 2 ∝ k 2 \widetilde{s}_{m=2}\propto k^2 s m = 2 ∝ k 2 ⟨ ( δ s ) s ~ ⟩ m = 2 ∝ k 2 \langle(\delta s)\widetilde{s}\rangle_{m=2}\propto k^2 ⟨( δ s ) s ⟩ m = 2 ∝ k 2 k → 0 k\to0 k → 0 ⟨ ∣ s ~ m = 2 ∣ 2 ⟩ ∝ k 4 \langle|\widetilde{s}_{m=2}|^2\rangle\propto k^4 ⟨ ∣ s m = 2 ∣ 2 ⟩ ∝ k 4 diverges as k − 2 k^{-2} k − 2
Beyond asymptotic scalings, the k k k Ψ 1 \boldsymbol{\Psi}_1 Ψ 1 (30) is the Fourier transform of configuration-space correlation functions:
⟨ ( δ s ~ ) a b ∣ s ~ c d ⟩ = ∫ q e − i k ⋅ q i k i ( U i ( q ) C a b c d ( q ) + E a b ( q ) B i c d ( q ) ) . \langle(\widetilde{\delta s})_{ab}|\widetilde{s}_{cd}\rangle=\int_{\boldsymbol{q}}e^{-i\boldsymbol{k}\cdot\boldsymbol{q}}ik_i\left(U_i(q)C_{abcd}(q)+E_{ab}(q)B_{icd}(q)\right). ⟨( δ s ) ab ∣ s c d ⟩ = ∫ q e − i k ⋅ q i k i ( U i ( q ) C ab c d ( q ) + E ab ( q ) B i c d ( q ) ) . When projected onto the m = 2 m=2 m = 2 6
⟨ ( δ s ~ ) m = 2 ∣ s ~ m = 2 ⟩ = ∫ q e − i k ⋅ q i k [ − 4 35 ( J 1 J 4 − ( 35 J 6 + 14 J 8 + J 9 ) U ( q ) ) P 1 ( μ ) + 8 45 ( J 1 J 4 − ( 9 J 8 + J 9 ) U ( q ) ) P 3 ( μ ) − 4 63 ( J 1 J 4 − J 9 U ( q ) ) P 5 ( μ ) ] . \langle(\widetilde{\delta s})_{m=2}|\widetilde{s}_{m=2}\rangle=\int_{\boldsymbol{q}}e^{-i\boldsymbol{k}\cdot\boldsymbol{q}}ik\left[-\frac{4}{35}\left(\mathcal{J}_1\mathcal{J}_4-(35\mathcal{J}_6+14\mathcal{J}_8+\mathcal{J}_9)U(q)\right)P_1(\mu)
+\frac{8}{45}\left(\mathcal{J}_1\mathcal{J}_4-(9\mathcal{J}_8+\mathcal{J}_9)U(q)\right)P_3(\mu)
-\frac{4}{63}\left(\mathcal{J}_1\mathcal{J}_4-\mathcal{J}_9U(q)\right)P_5(\mu)\right]. ⟨( δ s ) m = 2 ∣ s m = 2 ⟩ = ∫ q e − i k ⋅ q ik [ − 35 4 ( J 1 J 4 − ( 35 J 6 + 14 J 8 + J 9 ) U ( q ) ) P 1 ( μ ) + 45 8 ( J 1 J 4 − ( 9 J 8 + J 9 ) U ( q ) ) P 3 ( μ ) − 63 4 ( J 1 J 4 − J 9 U ( q ) ) P 5 ( μ ) ] . where P ℓ ( μ ) P_\ell(\mu) P ℓ ( μ ) ℓ \ell ℓ (30) was computed in resummed form in [20 ] , however the implementation in spinosaurus is numerically unstable at low-k k k m = 2 m=2 m = 2
⟨ s ~ a b ∣ s ~ c d ⟩ ⊃ − 1 2 ∫ q e − i k ⋅ q [ 1 2 k i k j A i j ( q ) C a b c d ( q ) + k i k j B i a b ( q ) B j c d ( q ) ] . \langle\widetilde{s}_{ab}|\widetilde{s}_{cd}\rangle\supset-\frac{1}{2}\int_{\boldsymbol{q}}e^{-i\boldsymbol{k}\cdot\boldsymbol{q}}\left[\frac{1}{2}k_i k_j A_{ij}(q)C_{abcd}(q)+k_i k_j B_{iab}(q)B_{jcd}(q)\right]. ⟨ s ab ∣ s c d ⟩ ⊃ − 2 1 ∫ q e − i k ⋅ q [ 2 1 k i k j A ij ( q ) C ab c d ( q ) + k i k j B iab ( q ) B j c d ( q ) ] . This expression admits a similar decomposition into standard LPT configuration-space correlators when the m = 2 m=2 m = 2 Figure 2 , we show the contribution of Equation (30) computed analytically, as well as measured empirically from the correlations of shifted tensor operators on the lattice. The measured and calculated ratios are in good agreement, and validate that at low-k k k c δ s c_{\delta_s} c δ s k k k m = 2 m=2 m = 2
Figure 2. Calculation of the δ s \delta_s δ s m = 2 m=2 m = 2 (30) , while the orange shows the measurement of this ratio from field-level realizations of shifted operators.
In general, for a purely quadrupolar operator ( s a b , t a b ) (s_{ab},t_{ab}) ( s ab , t ab ) m m m L a b ( 3 ) L_{ab}^{(3)} L ab ( 3 ) m = { 0 , ± 1 } 7 m=\{0,\pm1\}^{7} m = { 0 , ± 1 } 7 7 m = 2 m=2 m = 2
Introducing Tidal Galileons
In the EFT for spin-2 tracers, as in the scalar case, there exist certain degeneracies between the operator basis that hold at the field level. The first one is the degeneracy between cubic operators and the linear tidal field at the level of the one-loop power spectrum, which are understood as a renormalization of the linear bias [61 ] . We have already fully taken this into account for our cubic operators by subtracting from the Lagrangian operators their proper overlap with s a b s_{ab} s ab (10) .
However, there exist other degeneracies at the field level that show up already for the quadratic operators. This was first discussed in [23 ] , where it was shown that for the m = 1 m=1 m = 1
( s 2 ) m = 1 ( k ) = 1 3 ( δ s ) m = 1 ( k ) = 1 3 ( Ψ ⋅ ∇ s ) m = 1 , (s^2)_{m=1}(\mathbf{k})=\frac{1}{3}(\delta s)_{m=1}(\mathbf{k})=\frac{1}{3}(\mathbf{\Psi}\cdot\nabla s)_{m=1}, ( s 2 ) m = 1 ( k ) = 3 1 ( δ s ) m = 1 ( k ) = 3 1 ( Ψ ⋅ ∇ s ) m = 1 , and so this channel cannot disentangle contributions from these operators. A related question that one can ask, then, is whether the angular content of tensor operators in the low-k k k k → 0 k\to0 k → 0
Since we are interested in the k → 0 k \to0 k → 0 Ψ 1 \mathbf{\Psi}_1 Ψ 1
O a b ( k ) = ∫ p 1 , … , p n ( 2 π ) 3 δ D ( k − p 1 … n ) F a b ( p 1 , … , p n ) δ p 1 ⋯ δ p n . (32) \mathcal{O}_{ab}(\mathbf{k}) = \int_{\mathbf{p}_1,\ldots,\mathbf{p}_n} (2\pi)^3\delta_D(\mathbf{k}-\mathbf{p}_{1\ldots n})F_{ab}(\mathbf{p}_1,\ldots,\mathbf{p}_n)\delta_{\mathbf{p}_1}\cdots\delta_{\mathbf{p}_n}. \tag*{(32)} O ab ( k ) = ∫ p 1 , … , p n ( 2 π ) 3 δ D ( k − p 1 … n ) F ab ( p 1 , … , p n ) δ p 1 ⋯ δ p n . (32) For quadratic operators, the k → 0 k \to0 k → 0 F a b F_{ab} F ab p , − p \mathbf{p},-\mathbf{p} p , − p k k k S a b ( k ) = k ^ a k ^ b − δ a b / 3 S_{ab}(\mathbf{k}) = \hat{k}_a\hat{k}_b-\delta_{ab}/3 S ab ( k ) = k ^ a k ^ b − δ ab /3 D a b D_{ab} D ab [ s 2 ] a b [s^2]_{ab} [ s 2 ] ab
s a c s c b ( k → 0 ) = ∫ p S a c ( p ) S c b ( p ) δ p δ p = 1 3 ∫ p S a b ( p ) δ p δ p = 1 3 ( δ s ) a b ( k → 0 ) (33) s_{ac}s_{cb}(k\to0)=\int_{\mathbf{p}}S_{ac}(\mathbf{p})S_{cb}(\mathbf{p})\delta_{\mathbf{p}}\delta_{\mathbf{p}}=\frac{1}{3}\int_{\mathbf{p}}S_{ab}(\mathbf{p})\delta_{\mathbf{p}}\delta_{\mathbf{p}}=\frac{1}{3}(\delta s)_{ab}(k\to0) \tag*{(33)} s a c s c b ( k → 0 ) = ∫ p S a c ( p ) S c b ( p ) δ p δ p = 3 1 ∫ p S ab ( p ) δ p δ p = 3 1 ( δ s ) ab ( k → 0 ) (33) So we see that in the k → 0 k \to0 k → 0 s a b 2 , δ s a b s^2_{ab},\delta s_{ab} s ab 2 , δ s ab all helicities as opposed to the prior m = 1 m=1 m = 1 δ 2 \delta^2 δ 2 s 2 s^2 s 2 ⟨ δ 2 G 2 ⟩ → 0 \langle\delta^2G_2\rangle\to0 ⟨ δ 2 G 2 ⟩ → 0 k → 0 k\to0 k → 0 s 2 , δ s s^2,\delta s s 2 , δ s tidal Galileon
G a b s = s a b 2 − ( δ s ) a b / 3 , \mathcal{G}^{s}_{ab}=s^2_{ab}-(\delta s)_{ab}/3, G ab s = s ab 2 − ( δ s ) ab /3 , which is the contribution from the tidal torquing operator explicitly deprojected from the density-weighting contribution. The parameter for the tidal Galileon is the same as the usual tidal torquing operator, c G s = c s 2 c_{\mathcal{G}s}=c_{s^2} c G s = c s 2 c δ s → c δ s − c G s / 3 c_{\delta s}\to c_{\delta s}-c_{\mathcal{G}s}/3 c δ s → c δ s − c G s /3 G s \mathcal{G}^s G s m = 0 m=0 m = 0
G m = 0 s ( k ) = 7 8 t m = 0 ( k ) (34) \mathcal{G}^{s}_{m=0}(\mathbf{k})=\frac{7}{8}t_{m=0}(\mathbf{k}) \tag*{(34)} G m = 0 s ( k ) = 8 7 t m = 0 ( k ) (34) so there is a degeneracy with the velocity shear t a b t_{ab} t ab m = 0 m=0 m = 0 t t t
G m = 2 s ( k ) = − s ~ m = 2 ( k ) , (35) \mathcal{G}^{s}_{m=2}(\mathbf{k})=-\tilde{s}_{m=2}(\mathbf{k}), \tag*{(35)} G m = 2 s ( k ) = − s ~ m = 2 ( k ) , (35) so rotating the quadratic terms into the tidal Galileon basis introduces to the m = 2 m=2 m = 2 s ~ \tilde{s} s ~ 8 c s c_s c s k → 0 k\to0 k → 0 ( δ s ) a b (\delta s)_{ab} ( δ s ) ab P e r r P_{\rm err} P err
Do similar degeneracies exist for the cubic fields? For the cubic scalar bias operators, the mode coupling kernels necessarily decompose cubic operators in this limit into a universal form:
F i ( 3 ) ( k − p 1 − p 2 , p 1 , p 2 ) p 1 , p 2 → ∞ = A i + B i f ( p 1 , p 2 ) , F_i^{(3)}(\mathbf{k}-\mathbf{p}_1-\mathbf{p}_2,\mathbf{p}_1,\mathbf{p}_2)^{p_1,p_2\to\infty}=A_i+B_i f(\mathbf{p}_1,\mathbf{p}_2), F i ( 3 ) ( k − p 1 − p 2 , p 1 , p 2 ) p 1 , p 2 → ∞ = A i + B i f ( p 1 , p 2 ) , where f ( p 1 , p 2 ) f(\mathbf{p}_1,\mathbf{p}_2) f ( p 1 , p 2 ) [62 ] . There is a richer degeneracy structure for the spin-2 case in this same kinematic limit. For the cubic shape operators, there are four distinct angular structures in their mode-coupling kernels:
S a b = ∑ i S a b ( p i ) ; S a b ( μ ) = ∑ i μ j k 2 S a b ( p i ) ; χ a b = ∑ i < j μ i j [ p ^ i p ^ j ] a b ; χ a b ( μ ) = ∑ c y c μ i k μ k j [ p ^ i p ^ j ] a b , (36) S_{ab}=\sum_i S_{ab}(p_i);\quad S^{(\mu)}_{ab}=\sum_i\mu_{jk}^2S_{ab}(p_i);\quad\chi_{ab}=\sum_{i<j}\mu_{ij}[\hat{p}_i\hat{p}_j]_{ab};\quad\chi^{(\mu)}_{ab}=\sum_{\rm cyc}\mu_{ik}\mu_{kj}[\hat{p}_i\hat{p}_j]_{ab}, \tag*{(36)} S ab = i ∑ S ab ( p i ) ; S ab ( μ ) = i ∑ μ j k 2 S ab ( p i ) ; χ ab = i < j ∑ μ ij [ p ^ i p ^ j ] ab ; χ ab ( μ ) = cyc ∑ μ ik μ k j [ p ^ i p ^ j ] ab , (36) where i , j , k i,j,k i , j , k i , j , k i,j,k i , j , k L ( 3 ) L^{(3)} L ( 3 ) k k k δ 2 s a b \delta^2s_{ab} δ 2 s ab
δ 2 s a b ( k → 0 ) = ∫ p 1 , p 2 , p 3 ( 2 π ) 3 δ D ( p 123 ) [ 1 3 ∑ i = 1 3 S a b ( p i ) ] δ p 1 δ p 2 δ p 3 (37) \delta^2s_{ab}(k\to0)=\int_{\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3}(2\pi)^3\delta_D(\mathbf{p}_{123})\left[\frac{1}{3}\sum_{i=1}^{3}S_{ab}(p_i)\right]\delta_{p_1}\delta_{p_2}\delta_{p_3} \tag*{(37)} δ 2 s ab ( k → 0 ) = ∫ p 1 , p 2 , p 3 ( 2 π ) 3 δ D ( p 123 ) [ 3 1 i = 1 ∑ 3 S ab ( p i ) ] δ p 1 δ p 2 δ p 3 (37) δ 2 s a b ( k → 0 ) ≡ 1 3 S a b . (38) \delta^2s_{ab}(k\to0)\equiv\frac{1}{3}S_{ab}. \tag*{(38)} δ 2 s ab ( k → 0 ) ≡ 3 1 S ab . (38) and for the additional three Zel’dovich fields we have
s 2 s a b → S a b ( μ ) / 3 − S a b / 9 s^2s_{ab} \to S^{(\mu)}_{ab}/3 - S_{ab}/9 s 2 s ab → S ab ( μ ) /3 − S ab /9 δ s 2 → χ a b / 3 − 2 S a b / 9 \delta s^2 \to\chi_{ab}/3 - 2S_{ab}/9 δ s 2 → χ ab /3 − 2 S ab /9 s a b 3 → χ a b ( μ ) / 3 − χ a b / 3 + S a b / 9. s^3_{ab} \to\chi^{(\mu)}_{ab}/3 - \chi_{ab}/3 + S_{ab}/9. s ab 3 → χ ab ( μ ) /3 − χ ab /3 + S ab /9. This decomposition hides an additional degeneracy that exists at the field level: the operators s 2 s a b s^2s_{ab} s 2 s ab s a b 3 s^3_{ab} s ab 3 3 × 3 3 \times3 3 × 3 A A A
A 3 − ( tr A ) A 2 + 1 2 ( ( tr A ) 2 − tr ( A 2 ) ) A − det ( A ) I 3 = 0 , (39) A^3 - (\operatorname{tr} A)A^2 + \frac{1}{2}\left((\operatorname{tr} A)^2 - \operatorname{tr}(A^2)\right)A - \det(A)I_3 = 0, \tag*{(39)} A 3 − ( tr A ) A 2 + 2 1 ( ( tr A ) 2 − tr ( A 2 ) ) A − det ( A ) I 3 = 0 , (39) where the 0 on the right-hand side is understood as a matrix of zeros. Applying this to A = s a b A = s_{ab} A = s ab tr ( s a b ) = 0 \operatorname{tr}(s_{ab}) = 0 tr ( s ab ) = 0
s a b 3 − 1 2 s 2 s a b − det ( s ) δ a b K = 0. (40) s^3_{ab} - \frac{1}{2}s^2s_{ab} - \det(s)\delta^K_{ab} = 0. \tag*{(40)} s ab 3 − 2 1 s 2 s ab − det ( s ) δ ab K = 0. (40) The trace of this equation implies s 3 / 3 = det ( s ) s^3/3 = \det(s) s 3 /3 = det ( s )
s a b 3 − s 3 3 δ a b K = 1 2 s 2 s a b , (41) s^3_{ab} - \frac{s^3}{3}\delta^K_{ab} = \frac{1}{2}s^2s_{ab}, \tag*{(41)} s ab 3 − 3 s 3 δ ab K = 2 1 s 2 s ab , (41) and so the s a b 3 s^3_{ab} s ab 3 s 2 s a b s^2s_{ab} s 2 s ab 9
The non-Zel’dovich operators do not generate new angular dependencies, and in fact suffer from an analogous ‘full degeneracy’ that { s 2 , δ s } \{s^2,\delta s\} { s 2 , δ s }
[ s t ] a b ( k → 0 ) = δ t a b ( k → 0 ) / 3 , [st]_{ab}(k \to0) = \delta t_{ab}(k \to0)/3, [ s t ] ab ( k → 0 ) = δ t ab ( k → 0 ) /3 , and its low-k k k
δ t a b ( k → 0 ) ∝ S a b ( μ ) − S a b , \delta t_{ab}(k \to0) \propto S^{(\mu)}_{ab} - S_{ab}, δ t ab ( k → 0 ) ∝ S ab ( μ ) − S ab , which is strictly contained within the s 2 s a b s^2s_{ab} s 2 s ab δ 2 s a b \delta^2s_{ab} δ 2 s ab 10 δ t \delta t δ t s t st s t s 2 s a b → G 2 s a b s^2s_{ab} \to\mathcal{G}_2s_{ab} s 2 s ab → G 2 s ab ( G 2 s a b , δ t a b , s t a b ) (\mathcal{G}_2s_{ab},\delta t_{ab},st_{ab}) ( G 2 s ab , δ t ab , s t ab ) k k k δ t , s t \delta t,st δ t , s t
G a b ( 3 a ) ≡ s t a b − 1 3 δ t a b , \mathcal{G}^{(3a)}_{ab} \equiv st_{ab} - \frac{1}{3}\delta t_{ab}, G ab ( 3 a ) ≡ s t ab − 3 1 δ t ab , G a b ( 3 b ) ≡ δ t a b − 2 7 G 2 s a b , \mathcal{G}^{(3b)}_{ab} \equiv\delta t_{ab} - \frac{2}{7}\mathcal{G}_2s_{ab}, G ab ( 3 b ) ≡ δ t ab − 7 2 G 2 s ab , such that they have no overlap with the four angular structures. These tidal Galileons are also, by construction, normal-ordered with the linear density field, i.e. their P 13 P_{13} P 13 c s c_s c s G a b s \mathcal{G}^s_{ab} G ab s
s 2 s a b → G 2 s a b s^2s_{ab} \to\mathcal{G}_2s_{ab} s 2 s ab → G 2 s ab δ s a b 2 → δ G a b s \delta s^2_{ab} \to\delta\mathcal{G}^s_{ab} δ s ab 2 → δ G ab s s a b 3 → G a c s s c b , s^3_{ab} \to\mathcal{G}^s_{ac}s_{cb}, s ab 3 → G a c s s c b , however eliminating s 2 s a b s^2s_{ab} s 2 s ab G 2 s a b \mathcal{G}_2s_{ab} G 2 s ab t a b t_{ab} t ab L a b ( 3 ) L^{(3)}_{ab} L ab ( 3 ) L a b ( 3 ) L^{(3)}_{ab} L ab ( 3 ) m = 0 m = 0 m = 0 m = 1 m = 1 m = 1 k → 0 k \to0 k → 0 δ s a b 2 \delta^2_{s ab} δ s ab 2 δ G a b s \delta G^s_{ab} δ G ab s ( G s ) a b (G^s)_{ab} ( G s ) ab P 33 P_{33} P 33 P r r m P^{m}_{rr} P r r m [20 ] :
1 s t o r d e r : s a b 2 n d o r d e r : δ s a b , G a b s , t a b 3 r d o r d e r : δ 2 s a b , δ G a b s , ( G s ) a b , G a b ( 3 a ) , G a b ( 3 b ) , L a b ( 3 ) (42) \begin{aligned}
1^{\mathrm{st}}\ \mathrm{order}:&\ s_{ab}\\
2^{\mathrm{nd}}\ \mathrm{order}:&\ \delta s_{ab},G^s_{ab},t_{ab}\\
3^{\mathrm{rd}}\ \mathrm{order}:&\ \delta^2s_{ab},\delta G^s_{ab},(G^s)_{ab},G^{(3a)}_{ab},G^{(3b)}_{ab},L^{(3)}_{ab} \tag*{(42)}
\end{aligned} 1 st order : 2 nd order : 3 rd order : s ab δ s ab , G ab s , t ab δ 2 s ab , δ G ab s , ( G s ) ab , G ab ( 3 a ) , G ab ( 3 b ) , L ab ( 3 ) (42) This basis isolates the contributions to low-k k k M \mathcal{M} M ( δ t , s t , L ( 3 ) ) (\delta t,st,L^{(3)}) ( δ t , s t , L ( 3 ) )
Extending the relation between the transfer function measured for one orthogonalized field and the shape bias coefficients which contribute to it from the linear case (Equation (29) ) to higher orders leads to the construction of a mixing matrix
M m ; i j ( k ) = ⟨ O ~ i ⊥ O ~ j ⟩ ( k ) ⟨ ∣ O ~ i ⊥ ∣ 2 ⟩ ( k ) (43) \mathcal{M}_{m;ij}(k)=\frac{\langle\widetilde{O}^{\perp}_{i}\widetilde{O}_{j}\rangle(k)}{\langle|\widetilde{O}^{\perp}_{i}|^2\rangle(k)} \tag*{(43)} M m ; ij ( k ) = ⟨ ∣ O i ⊥ ∣ 2 ⟩ ( k ) ⟨ O i ⊥ O j ⟩ ( k ) (43) which describes correlations between orthogonalized and non-orthogonalized fields. The mixing matrix for transfer functions is shown in Figure 3 . Several of the analytic degeneracies we have discussed in this section are visible in this figure: the low-k k k m = 2 m = 2 m = 2 s ~ \tilde{s} s ~ k k k t m = 0 t_{m=0} t m = 0 G m = 0 s G^s_{m=0} G m = 0 s k k k t a b t_{ab} t ab δ 2 s a b \delta^2s_{ab} δ 2 s ab δ G s \delta G^s δ G s k k k Figure 5 the curves corresponding to the C t ; m = 0 ( k ) C_{t;m=0}(k) C t ; m = 0 ( k ) C t ; m = 2 ( k ) C_{t;m=2}(k) C t ; m = 2 ( k ) k k k [27 , 63 ] , the influence of higher-loop contributions for higher-order/helicity transfer functions makes such an interpretation difficult. This complex structure motivates us to construct the full mixing matrix from lattice correlators of shifted fields.
Figure 3. Entries of the mixing matrix Equation (43) , for each helicity. The ∝ k − 2 \propto k^{-2} ∝ k − 2 M s , δ s ( m = 2 ) \mathcal{M}^{(m=2)}_{s,\delta s} M s , δ s ( m = 2 ) V = 50 ( G p c / h ) 3 V = 50\;(\mathrm{Gpc}/h)^3 V = 50 ( Gpc / h ) 3
Figure 5. Full set of transfer functions, C ^ i ; m ( k ) \hat{C}_{i;m}(k) C ^ i ; m ( k ) M ∈ [ 10 13.5 , 10 14 ] M ⊙ h − 1 M \in[10^{13.5},10^{14}]M_\odot h^{-1} M ∈ [ 1 0 13.5 , 1 0 14 ] M ⊙ h − 1
Finally, let us note that, the last three operators in the last line of Eq. (42) , one linear combination of coefficients is redundant, leaving two independent combinations at the level of the one-loop shape power spectrum. However, they are not redundant at the level of the shape trispectrum, which we effectively probe in our field-level analysis [64 , 65 , 66 ] .
Estimating coefficients from transfer functions
Having settled on a suitable basis which is minimally degenerate, we turn to measuring the coefficients of this basis. The transfer functions are related to the shape bias coefficients through a generalized least-squares (GLS) linear regression problem
C ⃗ m ( k ) = M m ( k ) ⋅ c ⃗ L , (44) \vec{C}_{m}(k)=\mathcal{M}_{m}(k)\cdot\vec{c}_{L}, \tag*{(44)} C m ( k ) = M m ( k ) ⋅ c L , (44) where note that the mixing matrix M m ( k ) \mathcal{M}_{m}(k) M m ( k ) M \mathcal{M} M 1 1 1 (27) ) modifies the matrix. The lowest-order inclusion comes from the S 3 S_3 S 3
M i j → M i j + ⟨ S ~ 3 ∣ O ~ i ⊥ ⟩ ⟨ ∣ O ~ i ⊥ ∣ 2 ⟩ δ 0 j K , (45) \mathcal{M}_{ij}\rightarrow\mathcal{M}_{ij}+\frac{\langle\widetilde{S}_{3}|\widetilde{O}^{\perp}_{i}\rangle}{\langle|\widetilde{O}^{\perp}_{i}|^2\rangle}\delta^{K}_{0j}, \tag*{(45)} M ij → M ij + ⟨ ∣ O i ⊥ ∣ 2 ⟩ ⟨ S 3 ∣ O i ⊥ ⟩ δ 0 j K , (45) and in the case of higher-loop shifts (with coefficients that are not c s c_s c s s a b s_{ab} s ab
C i ; m ( k ) = ∑ j c j M m ; j ( k ) + α m ; i k 2 + β m ; i k 4 , (46) C_{i;m}(k)=\sum_j c_j\mathcal{M}_{m;j}(k)+\alpha_{m;i}k^2+\beta_{m;i}k^4, \tag*{(46)} C i ; m ( k ) = j ∑ c j M m ; j ( k ) + α m ; i k 2 + β m ; i k 4 , (46) and the α m ; i , β m ; i \alpha_{m;i},\beta_{m;i} α m ; i , β m ; i M \mathcal{M} M
By construction, the transfer functions have no cross-covariance across different helicities and operators. We also neglect k k k k k k (44) with weights given by a diagonal precision matrix (Equation (48) ). While for large ensembles of dark matter-only simulations, such as the Quijote suite, we can empirically measure the scatter of transfer functions, this is not possible for the FLAMINGO suite. We thus implement an analytic approximation to the transfer function variance.
For a given field-level realization, the shifted operator fields are fixed. It is the statistics of ϵ m \epsilon_m ϵ m ϵ m \epsilon_m ϵ m 11
σ 2 ( C i ; m ) ( k ) = a m N m o d e s ( k ) P e r r , m ( k ) ⟨ ∣ O ~ i ; m ⊥ ∣ 2 ⟩ ( k ) , (47) \sigma^2(C_{i;m})(k)=\frac{a_m}{N_{\mathrm{modes}}(k)}\frac{P_{\mathrm{err},m}(k)}{\langle|\widetilde{O}^{\perp}_{i;m}|^2\rangle(k)}, \tag*{(47)} σ 2 ( C i ; m ) ( k ) = N modes ( k ) a m ⟨ ∣ O i ; m ⊥ ∣ 2 ⟩ ( k ) P err , m ( k ) , (47) where a m = { 1 , 1 / 2 , 1 / 2 } a_m=\{1,1/2,1/2\} a m = { 1 , 1/2 , 1/2 } m = { 0 , 1 , 2 } m=\{0,1,2\} m = { 0 , 1 , 2 } N m o d e s ( k ) = 4 π ( k / k F ) 2 N_{\mathrm{modes}}(k)=4\pi(k/k_{\mathrm{F}})^2 N modes ( k ) = 4 π ( k / k F ) 2 Δ k = k F = 2 π / L b o x \Delta k=k_{\mathrm{F}}=2\pi/L_{\mathrm{box}} Δ k = k F = 2 π / L box (47) against the scatter in the transfer functions of halo shapes in Appendix D. Given weights from Equation (47) , the linear regression formulation of the perturbative field-level analysis immediately yields an analytic expression for shape bias coefficients
c ^ L = ( M T ⋅ W ⋅ M ) − 1 ⋅ M T ⋅ W ⋅ C , (48) \hat{c}_L=(\mathcal{M}^{T}\cdot W\cdot\mathcal{M})^{-1}\cdot\mathcal{M}^{T}\cdot W\cdot C, \tag*{(48)} c ^ L = ( M T ⋅ W ⋅ M ) − 1 ⋅ M T ⋅ W ⋅ C , (48) with W ( k ) ≡ [ σ 2 ( C i ; m ) ( k ) ] − 1 W(k)\equiv[\sigma^2(C_{i;m})(k)]^{-1} W ( k ) ≡ [ σ 2 ( C i ; m ) ( k ) ] − 1 [67 ] . Both C C C M \mathcal{M} M np.linalg.lstsq routine12
F = M T ⋅ W ⋅ M , (49) F=\mathcal{M}^{T}\cdot W\cdot\mathcal{M}, \tag*{(49)} F = M T ⋅ W ⋅ M , (49) Uncertainties for the field-level shape bias coefficients follow directly from inverting this matrix, while correctly marginalizing over the impact of including the large number of broadband terms in this fit. This analytic parameter covariance is compared to the ensemble Quijote covariance in Appendix D and shown to agree qualitatively and quantitatively.
In Ref. [27 ] the components of the operator mixing matrix Equation (43) were calculated in perturbation theory assuming low-k k k 13 M i j ; m ( k ) \mathcal{M}_{ij;m}(k) M ij ; m ( k )
The fiducial analysis choice we make in this work is to fit the linear and quadratic transfer functions to k max = 0.4 h M p c − 1 k_{\max}=0.4\,h\mathrm{Mpc}^{-1} k m a x = 0.4 h Mpc − 1 k max = 0.07 h M p c − 1 k_{\max}=0.07\,h\mathrm{Mpc}^{-1} k m a x = 0.07 h Mpc − 1 m ≠ 0 m\ne0 m = 0 L ( 3 ) L^{(3)} L ( 3 ) c i c_i c i k max = 0.4 h M p c − 1 k_{\max}=0.4\,h\mathrm{Mpc}^{-1} k m a x = 0.4 h Mpc − 1 χ 2 / d o f ≈ 1 \chi^2/\mathrm{dof}\approx1 χ 2 / dof ≈ 1 k max k_{\max} k m a x χ 2 \chi^2 χ 2
Simulations and shape definitions
In this work we study the shapes of particles bound to dark matter halos across two sets of simulations: the dark matter-only Quijote suite of simulations and the FLAMINGO suite of hydrodynamical simulations. We restrict ourselves to z = 0.5 z=0.5 z = 0.5
For Quijote, we use halo catalogs derived from MPI-Rockstar [68 ] , which uses a phase space-based approach to halo finding and calculates the inertia tensors of output halos. We run MPI-Rockstar on the “High Resolution” suite, which has N = ( 1024 ) 3 N=(1024)^3 N = ( 1024 ) 3 L = 1000 h − 1 M p c L=1000\,h^{-1}\mathrm{Mpc} L = 1000 h − 1 Mpc reduced and iterative inertia tensor [69 ] . That is, for every bound particle within a sphere of radius R v i r R_{\rm vir} R vir
I a b r e d = 1 N ∑ n ( x n − x ˉ ) a ( x n − x ˉ ) b ∣ x n − x ˉ ∣ 2 , (50) I^{\mathrm{red}}_{ab}=\frac{1}{N}\sum_n\frac{(x_n-\bar{x})_a(x_n-\bar{x})_b}{\lvert x_n-\bar{x}\rvert^2}, \tag*{(50)} I ab red = N 1 n ∑ ∣ x n − x ˉ ∣ 2 ( x n − x ˉ ) a ( x n − x ˉ ) b , (50) where x ˉ \bar{x} x ˉ x n x_n x n (50) are used to define a subset of particles, and the reduced tensor is re-computed. MPI-Rockstar by default will compute this procedure until convergence, or up to 10 iterations. We also compute non-iterative shape tensors (using only the first spherical pass), for both the reduced inertia tensor and the ‘standard’ inertia tensor
I a b = 1 N ∑ n ( x n − x ˉ ) a ( x n − x ˉ ) b , (51) I_{ab}=\frac{1}{N}\sum_n(x_n-\bar{x})_a(x_n-\bar{x})_b, \tag*{(51)} I ab = N 1 n ∑ ( x n − x ˉ ) a ( x n − x ˉ ) b , (51) as well as iterative ‘standard’ tensors. That is, we will look at four definitions of halo shapes in Quijote, spanned by the 2x2 matrix of ‘Reduced / Standard’ and ‘Iterative / non-Iterative’. We use 50 Quijote-HR boxes of the 100 available in the suite. The resulting volume of V = 50 ( G p c h − 1 ) 3 V=50(\mathrm{Gpc}\,h^{-1})^3 V = 50 ( Gpc h − 1 ) 3
In Figure 4 we show the m = 0 m=0 m = 0 m = 2 m=2 m = 2 ∝ k − 2 \propto k^{-2} ∝ k − 2 m = 2 m=2 m = 2
Figure 4. Linear transfer functions, C ^ s ; m ( k ) \widehat{C}_{s;m}(k) C s ; m ( k ) m = 0 m=0 m = 0 m = 2 m=2 m = 2 c s c_s c s
Turning to the FLAMINGO suite, in this work we restrict ourselves to a subset of the entire suite of simulations. The use of FLAMINGO allows us to establish a baseline between the shapes of CDM halos (and compare to Quijote), but also allows us to study the shapes of galaxies in a hydrodynamical simulation (defined as the inertia tensors of stars bound to a halo). Since the FLAMINGO simulation adopts a different cosmology, a different halo finder, and slightly different algorithms for determining iterative halo shapes we do not necessarily expect 1:1 agreement with the Quijote results. Specifically, we will measure the shapes of halos and galaxies from the L1_M9_DMO, L1_M9 and L1_M8 simulations. L1 denotes these simulations correspond to a box length of L = 1000 M p c L = 1000\mathrm{Mpc} L = 1000 Mpc L = 681 h − 1 M p c L = 681 h^{-1}\mathrm{Mpc} L = 681 h − 1 Mpc m g a s ≈ 10 9 M ⊙ m_{\mathrm{gas}} \approx10^9M_\odot m gas ≈ 1 0 9 M ⊙ m C D M = 5.65 × 10 9 M ⊙ m_{\mathrm{CDM}} = 5.65 \times10^9M_\odot m CDM = 5.65 × 1 0 9 M ⊙ { m g a s , m C D M } = { 1.34 , 7.06 } × 10 8 M ⊙ \{m_{\mathrm{gas}},m_{\mathrm{CDM}}\} = \{1.34, 7.06\} \times10^8M_\odot { m gas , m CDM } = { 1.34 , 7.06 } × 1 0 8 M ⊙ 14 15 [70 ] + SOAP [71 ] (FLAMINGO’s default halo finder and code to compute halo properties), finding strong numerical agreement and validating our reconstruction.
We generate FLAMINGO initial conditions by running monofonIC [72 , 73 ] using the publicly available config files used to run the simulations. The cosmology is the same as the fiducial DES Y3 cosmology, and the phases are set by panphasia [74 ] using the identical descriptor. Two key differences are that we set z s t a r t = 0.5 z_{\mathrm{start}} = 0.5 z start = 0.5 δ c b \delta_{\mathrm{cb}} δ cb N c e l l = 900 N_{\mathrm{cell}} = 900 N cell = 900
Given a catalog of halo positions and inertia tensors, we produce the trace-free shape fluctuation field, γ a b ( x ) \gamma_{ab}(\boldsymbol{x}) γ ab ( x )
γ a b ( x ) = I a b ( x ) − ⟨ I a b ⟩ Tr ⟨ I a b ⟩ , (52) \gamma_{ab}(\boldsymbol{x}) = \frac{I_{ab}(\boldsymbol{x})-\langle I_{ab}\rangle}{\operatorname{Tr}\langle I_{ab}\rangle}, \tag*{(52)} γ ab ( x ) = Tr ⟨ I ab ⟩ I ab ( x ) − ⟨ I ab ⟩ , (52) which is consistent with normalizations of halo shape fields from prior works [23 ] . In Figure 1 we show both the halo number density for a given mass bin, as well as the inverse Fourier-transformed helicity fields for m = { 0 , 1 , 2 } m = \{0,1,2\} m = { 0 , 1 , 2 } m = 0 m=0 m = 0 m = { 1 , 2 } m=\{1,2\} m = { 1 , 2 }
Alignments of dark matter halos
We begin by studying the field-level distribution of 3D shapes of dark matter halos. Halos are selected by mass, in bins of width 0.5 dex, with lower bin edges log 10 ( M min / M ⊙ h − 1 ) = { 12.5 , 13.0 , 13.5 , 14.0 , 14.5 } \log_{10}(M_{\min}/M_\odot h^{-1}) = \{12.5, 13.0, 13.5, 14.0, 14.5\} log 10 ( M m i n / M ⊙ h − 1 ) = { 12.5 , 13.0 , 13.5 , 14.0 , 14.5 } M ∈ [ 10 13.5 , 10 14 ] M ⊙ h − 1 M \in[10^{13.5}, 10^{14}]M_\odot h^{-1} M ∈ [ 1 0 13.5 , 1 0 14 ] M ⊙ h − 1 Figure 5 as an illustrative example. Transfer functions for all other mass bins are shown in Appendix § 0.5.5 . We first present the error power spectra of the best-fit transfer functions for each simulation, and compare the performance of the field-level EFT against popular IA models (the ‘nonlinear linear alignment’ model and ‘tidal alignment and tidal torquing’ models). We will then discuss the best-fit coefficients extracted from these measurements. Finally, we will use the large ensemble of Quijote simulations to compare the parameters inferred from transfer functions to the parameters inferred from fitting the four power spectra of halo shapes (the ⟨ γ m = 0 δ h ⟩ \langle\gamma_{m=0}\delta_h\rangle ⟨ γ m = 0 δ h ⟩
Error power spectra
Having fit transfer functions to the shape fluctuation field which minimize the error power spectrum, we construct the helicity-dependent residual field ϵ m ( k ) \epsilon_m(\boldsymbol{k}) ϵ m ( k ) ϵ m \epsilon_m ϵ m σ γ 2 \sigma_\gamma^2 σ γ 2
σ γ 2 = 1 5 N h ∑ n = 0 N b γ i j , n γ n i j , (53) \sigma_\gamma^2 = \frac{1}{5N_h}\sum_{n=0}^{N_b}\gamma_{ij,n}\gamma^{ij}_n, \tag*{(53)} σ γ 2 = 5 N h 1 n = 0 ∑ N b γ ij , n γ n ij , (53) where γ n \gamma_n γ n [20 ] :
P e r r m = 0 ( k ) = A 0 σ γ 2 n ˉ + k 2 A 21 + ⋯ (54) P_{\mathrm{err}}^{m=0}(k) = A_0 \frac{\sigma_\gamma^2}{\bar{n}} + k^2 A_{21} + \cdots \tag*{(54)} P err m = 0 ( k ) = A 0 n ˉ σ γ 2 + k 2 A 21 + ⋯ (54) P e r r m = 1 ( k ) = A 0 σ γ 2 n ˉ + k 2 ( A 21 + A 22 ) + ⋯ (55) P_{\mathrm{err}}^{m=1}(k) = A_0 \frac{\sigma_\gamma^2}{\bar{n}} + k^2(A_{21} + A_{22}) + \cdots \tag*{(55)} P err m = 1 ( k ) = A 0 n ˉ σ γ 2 + k 2 ( A 21 + A 22 ) + ⋯ (55) P e r r m = 2 ( k ) = A 0 σ γ 2 n ˉ + k 2 ( A 21 + 4 A 22 ) + ⋯ (56) P_{\mathrm{err}}^{m=2}(k) = A_0 \frac{\sigma_\gamma^2}{\bar{n}} + k^2(A_{21} + 4A_{22}) + \cdots \tag*{(56)} P err m = 2 ( k ) = A 0 n ˉ σ γ 2 + k 2 ( A 21 + 4 A 22 ) + ⋯ (56) Any significant low-k k k
In Figure 6 we show the error power spectra, normalized by the shape noise amplitude, for the five halo mass bins in question. The reported uncertainties are the uncertainties on the mean from averaging over N = 50 N=50 N = 50 n ˉ P e r r ≈ 1 \bar{n}P_{\mathrm{err}} \approx1 n ˉ P err ≈ 1 [75 ] . The presence of sub-Poisson stochasticity hints at the imprint of the ‘halo exclusion’ mechanism manifesting itself in the alignments of halo shapes. The amplitude is non-monotonic with halo mass, which is the same behavior observed for halo clustering – there exists an interplay between the amplitude of the exclusion contribution and the amplitude of the shot noise. For the two lowest mass bins (densest samples), the error power spectra appear to increase at low-k k k P e r r P_{\mathrm{err}} P err
Figure 6. Error power spectra of Quijote halo shapes, re-scaled by the expected shape noise amplitude σ γ 2 / n ˉ \sigma_\gamma^2/\bar{n} σ γ 2 / n ˉ k k k k k k
Comparison with Popular IA Models
It is also interesting to ask the question of whether this full field-level EFT is required in order to describe the halo shape field. To this end, we compare our field-level pipeline to two standard models commonly adopted in the intrinsic alignment literature – the nonlinear alignment model (NLA) and TATT (tidal alignment and tidal torquing) [1 , 22 , 50 ] . The NLA model is defined by defining the halo shape field as being 100% correlated (up to stochasticity) with the nonlinear tidal field of matter
γ a b ( x ) = c S ( ∂ a ∂ b ∂ 2 − δ a b K 3 ) δ m N L ( x ) + ϵ a b ( x ) . (57) \gamma_{ab}(x) = c_{\mathrm{S}}\left(\frac{\partial_a\partial_b}{\partial^2} - \frac{\delta_{ab}^{K}}{3}\right)\delta_m^{\mathrm{NL}}(x) + \epsilon_{ab}(x). \tag*{(57)} γ ab ( x ) = c S ( ∂ 2 ∂ a ∂ b − 3 δ ab K ) δ m NL ( x ) + ϵ ab ( x ) . (57) It is immediately apparent that this model cannot work at the field level – it predicts zero power for m = 1 m = 1 m = 1 m = 2 m = 2 m = 2
c s D a b δ m ( k ) → C s ; m ( k ) D a b δ m ( k ) c_sD_{ab}\delta_m(\mathbf{k}) \to C_{s;m}(k)D_{ab}\delta_m(\mathbf{k}) c s D ab δ m ( k ) → C s ; m ( k ) D ab δ m ( k ) and find the Eulerian transfer function which minimizes the m = 0 m = 0 m = 0
( δ s ) a b ( x ) , ( s 2 ) a b ( x ) , (\delta_s)_{ab}(\mathbf{x}), (s^2)_{ab}(\mathbf{x}), ( δ s ) ab ( x ) , ( s 2 ) ab ( x ) , and the ‘velocity shear’ contribution t a b ( x ) t_{ab}(\mathbf{x}) t ab ( x ) R R R R = 0 h − 1 M p c R = 0\,h^{-1}\mathrm{Mpc} R = 0 h − 1 Mpc R = 20 h − 1 M p c R = 20\,h^{-1}\mathrm{Mpc} R = 20 h − 1 Mpc [23 , 56 ] ). This implementation is consistent with how galaxy surveys such as DES and HSC currently implement their intrinsic alignment models [76 ] , where the linear contribution is sourced from the nonlinear matter density field but the higher order terms are computed perturbatively. Each of these additional Eulerian operators also have their shape biases promoted to transfer functions.
Figure 7 shows the comparison between NLA, TATT and the field-level EFT for the M h ∈ [ 10 13.5 , 10 14.0 ] M ⊙ h − 1 M_h \in[10^{13.5}, 10^{14.0}]M_\odot h^{-1} M h ∈ [ 1 0 13.5 , 1 0 14.0 ] M ⊙ h − 1
Figure 7. Comparison of cross-correlation coefficients achieved for Eulerian IA models and the 10-parameter field-level EFT. Each model is slightly horizontally offset so the points do not overlap. The m = 1 m=1 m = 1 m = 2 m=2 m = 2 k k k k k k r m r_m r m
r m ( k ) = P γ , f i t m ( k ) P γ m ( k ) P f i t m ( k ) . (58) r_m(k) = \frac{P^{m}_{\gamma,\mathrm{fit}}(k)}{\sqrt{P^{m}_{\gamma}(k)P^{m}_{\mathrm{fit}}(k)}}. \tag*{(58)} r m ( k ) = P γ m ( k ) P fit m ( k ) P γ , fit m ( k ) . (58) For m = 0 m = 0 m = 0 m = 1 , 2 m = 1, 2 m = 1 , 2 k ≳ 0.1 h M p c − 1 k \gtrsim0.1\,h\mathrm{Mpc}^{-1} k ≳ 0.1 h Mpc − 1
The bottom panel shows the per-wavenumber improvement in cross-correlation coefficient Δ r m \Delta r_m Δ r m m = 1 , 2 m = 1, 2 m = 1 , 2 R = 20 h − 1 M p c R = 20\,h^{-1}\mathrm{Mpc} R = 20 h − 1 Mpc Δ r m = 1 , 2 \Delta r_{m=1,2} Δ r m = 1 , 2 k > 0.2 h M p c − 1 k > 0.2\,h\mathrm{Mpc}^{-1} k > 0.2 h Mpc − 1
Since the addition of t a b t_{ab} t ab Δ r m \Delta r_m Δ r m m = 0 m = 0 m = 0 Δ r m = 0 \Delta r_{m=0} Δ r m = 0 k k k k > 0.2 h M p c − 1 k > 0.2\,h\mathrm{Mpc}^{-1} k > 0.2 h Mpc − 1 Δ r m = 0 \Delta r_{m=0} Δ r m = 0
The upshot of our comparison is that full EFT outperforms all popular models at the field level, and this result is not driven by additional parameters. Indeed, having a correct implementation of the displacement is crucial for an optimal cross-correlation coefficient. The non-linear displacement of the EFT fields is captured by the Zel’dovich field, which does not involve any free parameters. At the same time, the filtered and unfiltered TATT models have a very different performance as well, even though they are characterized by the same set of parameters. Our results thus strongly motivate the use of the full EFT calculations to model the intrinsic alignments in the cosmic shear data.
The error power spectra obtained from the field-level EFT show the transfer functions are adequately capturing the structure of halo shape correlations. We now turn to extract the shape bias coefficients from the transfer functions of dark matter halo shapes. We focus on the coefficients { c i } \{c_i\} { c i } (44) . We restrict the design matrix to include the fiducial choice of scales discussed in § IIID: k max = 0.4 h M p c − 1 k_{\max} = 0.4\,h\mathrm{Mpc}^{-1} k m a x = 0.4 h Mpc − 1 k max = 0.07 h M p c − 1 k_{\max} = 0.07\,h\mathrm{Mpc}^{-1} k m a x = 0.07 h Mpc − 1 χ 2 / d o f \chi^2/\mathrm{dof} χ 2 / dof χ 2 \chi^2 χ 2 ± 0.05 \pm0.05 ± 0.05 d o f = 920 \mathrm{dof} = 920 dof = 920 Figure 8 . Note that the lowest mass bin has the most heavy-tailed distribution – this is the bin with the lowest shape noise, and so transfer functions are measured at higher significance than for the more massive bins. The fit to the mean transfer function has a very high χ 2 \chi^2 χ 2 χ 2 / d o f ∼ 1.5 \chi^2/\mathrm{dof} \sim1.5 χ 2 / dof ∼ 1.5 k max = 0.05 h M p c − 1 k_{\max} = 0.05\,h\mathrm{Mpc}^{-1} k m a x = 0.05 h Mpc − 1 k max = 0.2 h M p c − 1 k_{\max} = 0.2\,h\mathrm{Mpc}^{-1} k m a x = 0.2 h Mpc − 1
Figure 8. Distribution of χ 2 \chi^2 χ 2 χ 2 \chi^2 χ 2
In Figure 9 we plot the linear and quadratic shape bias coefficients which result from the transfer function fit, across Quijote and FLAMINGO -DMO, for the four shape definitions considered in this work. The first panel shows the scaling of the linear alignment amplitude, c s c_s c s c s c_s c s [48 ] . We also observe broad consistency between FLAMINGO and Quijote alignment amplitudes, despite vastly different halo finding algorithms and different definitions for how bound particles are used to compute halo shape. The shape bias of dark matter halos is small compared to their clustering bias, ∣ c s ∣ ≲ 0.4 \lvert c_s\rvert\lesssim0.4 ∣ c s ∣ ≲ 0.4 Lagrangian shape biases, and so comparison with past measurements of higher-order shape biases (such as [23 ] ) should be done with an appropriate conversion. This does not affect the linear coefficient but does affect higher-order coefficients.
Figure 9. Measurement of linear and quadratic IA coefficients from field-level transfer functions, as a function of halo mass and halo shape definition. The uncertainties are estimated from the ensemble of N = 50 N = 50 N = 50 y y y
The density-weighted alignment coefficient c δ s c_{\delta_s} c δ s c s c_s c s c s ≈ − 0.1 c_s \approx-0.1 c s ≈ − 0.1 c s c_s c s c s ( M h ) c_s(M_h) c s ( M h ) c s 2 c_s^2 c s 2 c t c_t c t c s 2 c_s^2 c s 2 c t c_t c t [23 ] , who infer the trends in Lagrangian shape biases from deviations of their measured Eulerian biases from the co-evolution prediction. A careful comparison between their (Eulerian) bispectrum method and our Lagrangian field-level technique is of interest, but is not straightforward in practice. One effect that makes such comparison difficult is the use of an explicit cutoff in [23 ] , which significantly affects shape bias parameter measurements. The cutoff effect can be mitigated by an explicit computation of the renormalization-group evolution of bias parameters, which we defer to future work.
Cubic alignment coefficients
We now turn to the cubic alignment coefficients inferred from our GLS fits. To our knowledge this is the first time cubic coefficients are measured in the literature – in Ref. [20 ] their posteriors were consistent with zero when fit to the helicity auto-spectra of halos with a volume of V = 8 ( h − 1 G p c ) 3 V = 8\,(h^{-1}\mathrm{Gpc})^3 V = 8 ( h − 1 Gpc ) 3 Figure 10 . Their behavior is qualitatively similar to what is observed in the linear and quadratic sector – increasing ∣ c s ∣ |c_s| ∣ c s ∣ M h ∈ [ 10 13.5 , 10 14.0 ] M ⊙ h − 1 M_h \in[10^{13.5},10^{14.0}]M_\odot h^{-1} M h ∈ [ 1 0 13.5 , 1 0 14.0 ] M ⊙ h − 1 M h ∈ [ 10 14.0 , 10 14.5 ] M ⊙ h − 1 M_h \in[10^{14.0},10^{14.5}]M_\odot h^{-1} M h ∈ [ 1 0 14.0 , 1 0 14.5 ] M ⊙ h − 1 P e r r P_{\rm err} P err k k k c δ 2 s c_{\delta^2s} c δ 2 s c δ G s c_{\delta Gs} c δ G s c G s s c_{Gss} c G ss
Figure 10. Constraints on cubic shape biases and impact on P e r r P_{\rm err} P err Left: The values of cubic shape bias parameters, as a function of c s c_s c s Right: The auto contributions of the three invariant sub-spaces of cubic parameters, { C i } \{C_i\} { C i } k 2 k^2 k 2 k 4 k^4 k 4
To investigate this, we study the spectrum of the cubic parameter covariance matrix. Specifically, given the 10x10 parameter covariance matrix, C C C
C = U Λ U T , C ⋅ u k = λ k u k , (59) C = U\Lambda U^T,\quad C \cdot u_k = \lambda_k u_k, \tag*{(59)} C = U Λ U T , C ⋅ u k = λ k u k , (59) where U ≡ u i k U \equiv u_{ik} U ≡ u ik C C C λ \lambda λ Λ \Lambda Λ k k k
S N R k = ∣ u i k c i ∣ λ k . (60) \mathrm{SNR}_k = \frac{|u_{ik}c_i|}{\sqrt{\lambda_k}}. \tag*{(60)} SNR k = λ k ∣ u ik c i ∣ . (60) In this decomposition, we find that for all mass bins, all eigenmodes are measured with SNR > 1 > 1 > 1 M h ∈ [ 10 14.0 , 10 14.5 ] M ⊙ h − 1 M_h \in[10^{14.0},10^{14.5}]M_\odot h^{-1} M h ∈ [ 1 0 14.0 , 1 0 14.5 ] M ⊙ h − 1 B 0 = span { u 0 } B_0 = \operatorname{span}\{u_0\} B 0 = span { u 0 } B 1 = span { u 1 , u 2 , u 3 } B_1 = \operatorname{span}\{u_1,u_2,u_3\} B 1 = span { u 1 , u 2 , u 3 } 16 B 2 = span { u 4 , u 5 } B_2 = \operatorname{span}\{u_4,u_5\} B 2 = span { u 4 , u 5 } u 0 ⋅ c i ≈ c δ 2 s − c δ G s / 3 u_0 \cdot c_i \approx c_{\delta^2s}-c_{\delta Gs}/3 u 0 ⋅ c i ≈ c δ 2 s − c δ G s /3 c i c_i c i O G s \mathcal{O}_{Gs} O G s − δ 2 s a b / 3 -\delta^2s_{ab}/3 − δ 2 s ab /3 B 2 B_2 B 2 c G ( 3 a ) c_{\mathcal{G}(3a)} c G ( 3 a ) c G ( 3 b ) c_{\mathcal{G}(3b)} c G ( 3 b ) c L 3 c_{L3} c L 3
In general, we can define projector matrices Π a b B = ∑ k ∈ B U a k U b k \Pi^B_{ab} = \sum_{k\in B} U_{ak}U_{bk} Π ab B = ∑ k ∈ B U ak U bk k k k Figure 10 we show the contributions of these sub-spaces to the auto-spectrum of halo shapes for the m = 0 m=0 m = 0 M h ∈ [ 10 13.5 , 10 14.0 ] M ⊙ h − 1 M_h \in[10^{13.5},10^{14.0}]M_\odot h^{-1} M h ∈ [ 1 0 13.5 , 1 0 14.0 ] M ⊙ h − 1 B B 17 ( k ) B^{17}_B(k) B B 17 ( k )
P B m ( k ) = Π a b B Π d c B a c a d P b c m ( k ) P_B^m(k) = \Pi^B_{ab}\Pi^B_{dc}a_ca_dP^m_{bc}(k) P B m ( k ) = Π ab B Π d c B a c a d P b c m ( k ) Note this is an illustration of how these combinations of coefficients contribute to the low-k k k ⟨ Π a b B c a O b γ h ⟩ \langle\Pi^B_{ab}c_aO_{b\gamma}^h\rangle ⟨ Π ab B c a O bγ h ⟩ B 0 B_0 B 0 ∼ 2.5 % \sim2.5\% ∼ 2.5% B 1 B_1 B 1 k k k B 0 B_0 B 0 B 2 B_2 B 2 k 2 k^2 k 2
A particular correlation between operators that is large (and can be seen explicitly in Figure 17 ) is that of the cubic Galileons c G ( 3 b ) c_{\mathcal{G}(3b)} c G ( 3 b ) c G ( 3 a ) c_{\mathcal{G}(3a)} c G ( 3 a ) k 0 k^0 k 0 s t a b st_{ab} s t ab δ t a b \delta t_{ab} δ t ab k 2 k^2 k 2
Figure 17. Comparison between empirical measurements of field-level IA parameters and the analytic expression given by generalized least-squares, for one of the Quijote + rockstar halo samples.
G ( 3 b ) , G ( 3 a ) → G ( 3 a ) + 1 6 G ( 3 b ) , \mathcal{G}^{(3b)},\mathcal{G}^{(3a)} \to\mathcal{G}^{(3a)}+\frac{1}{6}\mathcal{G}^{(3b)}, G ( 3 b ) , G ( 3 a ) → G ( 3 a ) + 6 1 G ( 3 b ) , in which case one component scales as O ( k 2 ) O(k^2) O ( k 2 ) O ( k 4 ) O(k^4) O ( k 4 ) B 2 B_2 B 2 c L 3 c_{L3} c L 3 k 4 k^4 k 4 m = 0 m=0 m = 0
In conclusion, we are able to measure all six coefficients across mass bins in the field level analysis for halos, however a spectral decomposition of the marginalized parameter covariance matrix implies that three parameter sub-spaces dominate the overall field-level fit, two which contribute as constant low-k k k k 2 k^2 k 2 k k k
Comparison with fits to power spectra
Are the best-fit values for { c i } \{c_i\} { c i }
d = ( P γ δ ( k i ) , P γ γ m ( k ) ) , (61) \mathbf{d}=\left(P_{\gamma\delta}(k_i),P^m_{\gamma\gamma}(k)\right), \tag*{(61)} d = ( P γ δ ( k i ) , P γ γ m ( k ) ) , (61) where P γ δ ≡ ⟨ γ m = 0 δ N L ′ ⟩ ( k ) P_{\gamma\delta}\equiv\langle\gamma^{m=0}\delta_{\mathrm{NL}}'\rangle(k) P γ δ ≡ ⟨ γ m = 0 δ NL ′ ⟩ ( k ) m = 0 m=0 m = 0 P γ δ P_{\gamma\delta} P γ δ P γ γ m = 0 P^{m=0}_{\gamma\gamma} P γ γ m = 0 k k k m = 0 m=0 m = 0 m = { 1 , 2 } m=\{1,2\} m = { 1 , 2 }
− ln L ∝ ∑ k ( r T C − 1 r ) m = 0 ( k ) + ∑ m = 1 2 ∑ k ( P ^ γ γ m − P ‾ γ γ m σ γ γ , m ) 2 ( k ) , (62) -\ln L \propto\sum_k \left(\mathbf{r}^{\mathsf{T}}\mathbf{C}^{-1}\mathbf{r}\right)_{m=0}(k)+\sum_{m=1}^{2}\sum_k\left(\frac{\widehat{P}^{m}_{\gamma\gamma}-\overline{P}^{m}_{\gamma\gamma}}{\sigma_{\gamma\gamma,m}}\right)^2(k), \tag*{(62)} − ln L ∝ k ∑ ( r T C − 1 r ) m = 0 ( k ) + m = 1 ∑ 2 k ∑ ( σ γ γ , m P γ γ m − P γ γ m ) 2 ( k ) , (62) where the m = 0 m=0 m = 0 k k k
r ( k ) = ( P ^ γ δ − P ‾ γ δ , P ^ γ γ m = 0 − P ‾ γ γ m = 0 ) ( k ) . \mathbf{r}(k)=\left(\widehat{P}_{\gamma\delta}-\overline{P}_{\gamma\delta},\widehat{P}^{m=0}_{\gamma\gamma}-\overline{P}^{m=0}_{\gamma\gamma}\right)(k). r ( k ) = ( P γ δ − P γ δ , P γ γ m = 0 − P γ γ m = 0 ) ( k ) . The predictions for the shape-matter and shape-shape spectra follow from the bias expansion
P ˉ γ δ ( k ) = ∑ i c i ⟨ O ~ i m = 0 δ m ⟩ ( k ) + α k 2 P s ~ δ ( k ) + a 2 k 2 , \bar{P}_{\gamma\delta}(k)=\sum_i c_i\langle\widetilde{O}_i^{m=0}\delta_m\rangle(k)+\alpha k^2P_{\tilde{s}\delta}(k)+a_2k^2, P ˉ γ δ ( k ) = i ∑ c i ⟨ O i m = 0 δ m ⟩ ( k ) + α k 2 P s ~ δ ( k ) + a 2 k 2 , P ˉ γ γ m ( k ) = ∑ i , j c i c j ⟨ O ~ i m O ~ j m ⟩ + β k 2 P s ~ s ~ m ( k ) + N 0 + N m ( k ) , \bar{P}_{\gamma\gamma}^{m}(k)=\sum_{i,j}c_ic_j\langle\widetilde{O}_i^m\widetilde{O}_j^m\rangle+\beta k^2P_{\tilde{s}\tilde{s}}^{m}(k)+N_0+N_m(k), P ˉ γ γ m ( k ) = i , j ∑ c i c j ⟨ O i m O j m ⟩ + β k 2 P s ~ s ~ m ( k ) + N 0 + N m ( k ) , where α \alpha α β \beta β ∇ 2 s a b \nabla^2s_{ab} ∇ 2 s ab a 2 a_2 a 2 N 0 N_0 N 0 N m N_m N m (54) . While the shifted operator cross-spectra ⟨ O ~ i m O ~ j m ⟩ \langle\widetilde{O}_i^m\widetilde{O}_j^m\rangle ⟨ O i m O j m ⟩ N = 50 N=50 N = 50 P γ δ P_{\gamma\delta} P γ δ [20 ] , but the use of Zel’dovich-only displacements for the basis implies the coefficient values should not be directly compared (as explicitly shown in Equation (27) , higher-order displacements modify coefficient values).
In Figure 11 we show all four shape noise-subtracted spectra for Quijote, for five mass bins. The uncertainties shown are uncertainties on the mean, averaged over the full 50 ( G p c / h ) 3 (\mathrm{Gpc}/h)^3 ( Gpc / h ) 3 [20 ] , and despite the large volume the residuals are closely scattered around zero with no large excursions. This full-freedom fit is done through markov chain monte carlo using emcee [77 ] , in the Lagrangian shape bias basis. The shaded bands indicate the 68% posterior predictive distributions from fitting to the full volume. We fit to slightly larger scales of k max = 0.25 h M p c − 1 k_{\max}=0.25\,h\mathrm{Mpc}^{-1} k m a x = 0.25 h Mpc − 1
Figure 11. Fits to halo helicity power spectra and cross-spectra with model. Solid lines are the result from MCMC fits with a strict one-loop model. Dashed lines keep the biases fit to values from transfer functions, and refit the nuisance coefficients. The dotted lines are fits where a subset of the cubic coefficients and the nuisance parameters are free to vary. The shape–matter cross spectrum is normalized so its amplitude matches the m = 0 m=0 m = 0 k = 0.1 h M p c − 1 k=0.1\,h\mathrm{Mpc}^{-1} k = 0.1 h Mpc − 1 N = 50 N=50 N = 50
We compare the full-freedom MCMC fit to maximum-likelihood fits where subsets of the shape bias parameters are constrained to parameters obtained from the transfer function fits. The dashed lines are fits where all 10 { c i } \{c_i\} { c i } k ∼ 0.15 h M p c − 1 k\sim0.15\,h\mathrm{Mpc}^{-1} k ∼ 0.15 h Mpc − 1 4 σ 4\sigma 4 σ
In summary, we have carried out a comprehensive examination of the spin-2 field-level bias expansion when applied to dark matter halo shapes across two suites of simulations – FLAMINGO-DMO and Quijote-HiRes – and four different definitions of halo shape. We find that the field-level model produces error power spectra that are flat at low-k k k c i c_i c i c s c_s c s c g s c_{g_s} c g s c t c_t c t M h ∈ [ 10 13.5 , 10 14 ] M ⊙ h − 1 M_h\in[10^{13.5},10^{14}]M_\odot h^{-1} M h ∈ [ 1 0 13.5 , 1 0 14 ] M ⊙ h − 1
Galaxy shape alignments in hydrodynamic simulations
Having studied and calibrated the field-level alignment machinery against dark matter halos, we turn to measurements of the transfer functions of galaxy shapes in hydrodynamical simulations. We split galaxies in FLAMINGO into two separate categories: The first split is motivated by measurements of direct intrinsic alignments in spectroscopic galaxy surveys, such as BOSS and DESI. Alignments have been detected in spectroscopic catalogs of LRGs for nearly two decades [78 ] , primarily with galaxies at z ≲ 0.5 z\lesssim0.5 z ≲ 0.5 [15 ] , and potentially new samples which can be used to calibrate photometric surveys [16 ] . Reproducing LRG-like and ELG-like samples within a hydrodynamic simulation provides an interesting point of comparison to current data. In FLAMINGO, we define LRGs as galaxies which have a low specific star formation rate (sSFR < 0.1 G y r − 1 < 0.1\,\mathrm{Gyr}^{-1} < 0.1 Gyr − 1 M ⋆ > 1.8 × 10 11 M ⊙ M_\star> 1.8 \times10^{11} M_\odot M ⋆ > 1.8 × 1 0 11 M ⊙ n ˉ L R G = 6.8 ( 7.8 ) × 10 − 4 ( h M p c − 1 ) 3 \bar{n}_{\mathrm{LRG}} = 6.8(7.8) \times10^{-4}\,(h\mathrm{Mpc}^{-1})^3 n ˉ LRG = 6.8 ( 7.8 ) × 1 0 − 4 ( h Mpc − 1 ) 3 z ≈ 0.5 z \approx0.5 z ≈ 0.5 ≥ 0.1 G y r − 1 \ge0.1\,\mathrm{Gyr}^{-1} ≥ 0.1 Gyr − 1 M ⋆ > 10 11 M ⊙ M_\star> 10^{11} M_\odot M ⋆ > 1 0 11 M ⊙ n ˉ E L G = 4.0 ( 1.6 ) × 10 − 4 ( h M p c − 1 ) 3 \bar{n}_{\mathrm{ELG}} = 4.0(1.6) \times10^{-4}\,(h\mathrm{Mpc}^{-1})^3 n ˉ ELG = 4.0 ( 1.6 ) × 1 0 − 4 ( h Mpc − 1 ) 3 [79 ] . The decrease in the number density of ELGs from going to a higher mass resolution is driven by a lower star formation rate at fixed stellar mass in the higher-resolution FLAMINGO box.
The second set of sample selections are based on a morphological split of the galaxies. We sub-select catalogs into ‘spiral’ and ‘elliptical’ galaxies (following [80 ] ), by splitting galaxies along their fraction of baryonic kinetic energy in rotation, κ c o r o t \kappa_{\mathrm{corot}} κ corot κ c o r o t > 0.5 \kappa_{\mathrm{corot}} > 0.5 κ corot > 0.5 > 0.1 G y r − 1 > 0.1\,\mathrm{Gyr}^{-1} > 0.1 Gyr − 1 κ c o r o t < 0.5 \kappa_{\mathrm{corot}} < 0.5 κ corot < 0.5 < 0.1 G y r < 0.1\,\mathrm{Gyr} < 0.1 Gyr M ⋆ > 2 × 10 11 M_\star> 2 \times10^{11} M ⋆ > 2 × 1 0 11 5.6 ( 6.4 ) × 10 − 4 ( h M p c − 1 ) 3 5.6(6.4) \times10^{-4}\,(h\mathrm{Mpc}^{-1})^3 5.6 ( 6.4 ) × 1 0 − 4 ( h Mpc − 1 ) 3 6.0 ( 20 ) × 10 − 4 ( h M p c − 1 ) 3 6.0(20) \times10^{-4}\,(h\mathrm{Mpc}^{-1})^3 6.0 ( 20 ) × 1 0 − 4 ( h Mpc − 1 ) 3
We define galaxy shapes as the inertia tensor of the star particles bound to a given subhalo, once again reconstructing the SOAP inertia tensors for missing subhalos. Having defined the galaxy sample we proceed to make a second quality cut on their shape tensors. We find that certain galaxy shapes (when using iterative shape definitions) have nearly-singular inertia tensors. These discrepant galaxies contribute more than 40% of the variance in the shape tensor for the sample, inflating the shape noise (through σ γ 2 \sigma_\gamma^2 σ γ 2 σ γ 2 \sigma_\gamma^2 σ γ 2
In Table I we report the inferred values of linear and quadratic shape biases for the FLAMINGO sample, from the field-level methodology developed in this work, using the same scale cuts and broad-band choices which were calibrated on the Quijote dark matter halo ensemble. Uncertainties on the measured transfer functions are given by Equation (47) , and the error bars on parameters come from Equation (49) . The error power spectra for all galaxy samples, shape definitions and helicities are consistent with the shape noise expectation, and this is explicitly shown for the standard iterative shapes in Figure 12 . In Figure 13 , we show the c i ( c s ) c_i(c_s) c i ( c s ) c δ 2 s c_{\delta2s} c δ 2 s Table 1 .
Figure 12. Error power spectra of galaxy shape samples from FLAMINGO considered in this work, normalized to their shape noise expectation.
Figure 13. Shape bias coefficients for the four L1_m8 galaxy samples extracted from FLAMINGO in this work, compared to the dark matter halo samples. Uncertainties shown are 2 σ 2\sigma 2 σ
Table 1. Linear and quadratic tidal-Galileon shape-bias coefficients for the FLAMINGO galaxy samples at high (L1_M8) and normal (L1_M9) resolution, for iterative, unreduced stellar inertia tensors (sshape_iter). The reported uncertainties are marginalized over the 31 cubic and broad-band parameters. FLAMINGO analyses have dof = 568
For galaxies living in massive halos (LRGs / ellipticals), the inferred shape biases are similar to each other in terms of central values and uncertainties. Linear alignments are found to be within c s ∈ [ − 9.49 , − 8.34 ] × 10 − 2 c_s \in[-9.49, -8.34] \times10^{-2} c s ∈ [ − 9.49 , − 8.34 ] × 1 0 − 2 σ c s ∼ 4 × 10 − 3 \sigma_{c_s} \sim4 \times10^{-3} σ c s ∼ 4 × 1 0 − 3 ∼ 20 σ \sim20\sigma ∼ 20 σ c δ s c_{\delta s} c δ s > 40 σ > 40\sigma > 40 σ
The alignments of LRGs in FLAMINGO were previously studied in [47 , 48 ] using the L = 2.8 G p c L = 2.8\mathrm{Gpc} L = 2.8 Gpc Figure 2 of [48 ] , we take their fiducial value of A I A = 5 ± 0.2 A_{\mathrm{IA}} = 5 \pm0.2 A IA = 5 ± 0.2 z = 0 z = 0 z = 0 10 2 c s ≈ − 11.05 ± 0.44 10^2 c_s \approx-11.05 \pm0.44 1 0 2 c s ≈ − 11.05 ± 0.44 − 9.4 × 10 − 2 -9.4 \times10^{-2} − 9.4 × 1 0 − 2 [47 ] ). What’s notable is the size of the uncertainty: the 2.8 G p c 2.8\ \mathrm{Gpc} 2.8 Gpc ∼ 20 × \sim20\times ∼ 20 × S N R ( 2.8 G p c ) ∼ 11 / 0.44 ∼ 25 \mathrm{SNR}^{(2.8\mathrm{Gpc})} \sim11/0.44 \sim25 SNR ( 2.8 Gpc ) ∼ 11/0.44 ∼ 25 S N R ( 1.0 G p c ) ∼ 9.25 / 0.46 ∼ 20 \mathrm{SNR}^{(1.0\mathrm{Gpc})} \sim9.25/0.46 \sim20 SNR ( 1.0 Gpc ) ∼ 9.25/0.46 ∼ 20
Turning to the blue samples, the statistical power of the field-level analysis allows us to make, for the first time, precise measurements of the alignments of blue galaxies in the L1_M8 FLAMINGO box. For ELGs, their low number density do not allow for a robust detection of linear alignment or tidal galileon contribution, but the central value is comparable to the spiral result. Direct measurements of intrinsic alignments for ELGs in DESI report non-detections (Table 3 of [15 ] ). with σ ( A 1 ) ∼ 1 \sigma(A_1) \sim1 σ ( A 1 ) ∼ 1 σ ( A 2 ) ∼ 4 \sigma(A_2) \sim4 σ ( A 2 ) ∼ 4 σ ( b T A ) ∼ 1 \sigma(b_{TA}) \sim1 σ ( b T A ) ∼ 1 A 1 A_1 A 1 b T A b_{TA} b T A σ ( A 1 δ ) ∼ 1 \sigma(A_{1\delta}) \sim1 σ ( A 1 δ ) ∼ 1
A 1 = − 0.52 ± 0.56 , A 1 δ = − 1.01 ± 0.55 , A 2 = 0.15 ± 0.27. \begin{aligned}
A_1 &= -0.52 \pm0.56, \\
A_{1\delta} &= -1.01 \pm0.55, \\
A_2 &= 0.15 \pm0.27.
\end{aligned} A 1 A 1 δ A 2 = − 0.52 ± 0.56 , = − 1.01 ± 0.55 , = 0.15 ± 0.27. We thus see that our inference is consistent with the current non-detections of alignment in ELG for these coefficients. Our measurements have smaller error bars than what the two-point function analysis in data is able to achieve, despite the larger effective volume of the ELG sample of V e f f E L G = 2 G p c 3 V_{\mathrm{eff}}^{\mathrm{ELG}} = 2\mathrm{Gpc}^3 V eff ELG = 2 Gpc 3 [81 ] . We find mild evidence (3.5 σ 3.5\sigma 3.5 σ 10 2 c δ s = 1.47 ± 0.42 10^2 c_{\delta s} = 1.47 \pm0.42 1 0 2 c δ s = 1.47 ± 0.42 18
More interestingly, we now turn to the denser spiral galaxy sample. Here, we find a > 4 σ > 4\sigma > 4 σ linear c s c_s c s [49 , 82 , 83 ] . However, in contrast with past work we note the sign of spiral alignment amplitude for L1_m8 spirals does not depend on the shape definition used. We also detect non-zero density-weighting contribution (> 13 σ > 13\sigma > 13 σ and non-zero tidal Galileon term (> 6 σ > 6\sigma > 6 σ
A 1 = − 0.76 ± 0.16 , A 1 δ = − 0.77 ± 0.16 , A 2 = 0.34 ± 0.08. \begin{aligned}
A_1 &= -0.76 \pm0.16, \\
A_{1\delta} &= -0.77 \pm0.16, \\
A_2 &= 0.34 \pm0.08.
\end{aligned} A 1 A 1 δ A 2 = − 0.76 ± 0.16 , = − 0.77 ± 0.16 , = 0.34 ± 0.08. It is appealing to ask the question of what amount of contamination relative to a cosmic shear data vector do these values of A i A_i A i n ( z ) n(z) n ( z )
The application of field-level perturbative IA modelling to hydrodynamical simulations, or simulations which use empirical prescriptions to imbue alignments [84 ] , paves the way for the calibration of simulation-based priors (SBPs) for these probes in much the same way as has recently been proposed for galaxy densities and the Lyman-α \alpha α [63 , 85 , 86 , 87 ] . The utilities of SBPs for alignments are two-fold: they will enhance the constraining power from growth-of-structure constraints using alignments, and will also be important in bounding the allowed contamination of IA in cosmic shear measurements in the absence of direct calibration [88 ] .
Conclusions
In this work, we have constructed the Lagrangian spin-2 bias expansion at the level of the full real-space 3D shape fluctuation field, including all operators to cubic order. This is in contrast with past works that focused on summary statistics calculated within this model (or the accompanying Eulerian formulation of the EFT of IA). We stress-tested the field-level EFT’s ability to describe the fluctuation field of dark matter halo and galaxy shapes in numerical simulations, finding a higher cross-correlation coefficient between the field-level EFT and halo shapes than any of the popular IA models used in the literature (NLA or TATT). Additionally, we find that for the volume of a single Quijote realization the distribution of reduced χ 2 \chi^2 χ 2
We also found the field-level formulation of shape statistics presents some additional subtleties not present in the scalar case – namely, the transfer functions measured (C m ; i ( k ) C_{m;i}(k) C m ; i ( k ) k k k
We then showed that the same procedure used to stress-test the field-level EFT can be re-cast as a formal generalized least-squares regression problem. This led to an estimator for the shape bias coefficients of this expansion, as well as its accompanying uncertainties, which cancels cosmic variance in the measurement. Notably, compared to past field-level EFT publications, the GLS-based estimator for coefficients allows us to characterize the uncertainty and covariance on measured parameters with a simple analytic expression. This analytic expression is compared to the ensemble covariance matrix measured from 50 boxes and found to be in good agreement for both numerical value and the off-diagonal cross-correlation coefficients.
Applying this estimator to dark matter halo shapes in the Quijote simulation leads to precise measurements of all ten coefficients through cubic order, in a way that is consistent with coefficients inferred from fitting to the one-loop shape auto-spectrum and shape-matter cross-spectrum. We find that the higher order halo alignment coefficients generally grow with halo mass, and trends with the linear coefficient c s c_s c s c g s , c t c_{g_s}, c_t c g s , c t ∣ c s ∣ ≳ 0.1 \lvert c_s \rvert\gtrsim0.1 ∣ c s ∣ ≳ 0.1 k k k c g ( 3 a ) , c g ( 3 b ) , c L 3 c_{g(3a)}, c_{g(3b)}, c_{L3} c g ( 3 a ) , c g ( 3 b ) , c L 3
As a final application we turned to the case of analyzing samples of galaxies in the FLAMINGO hydrodynamic simulation. Two sets of sample selections were made, one focusing on spectroscopic samples and another focusing on a morphological classification. While more massive / red galaxies were found to be similar across mass resolutions, the higher-resolution L1_m8 box is needed to resolve the properties of blue galaxies which are relevant for studies of ELG alignment and cosmic shear contamination. Our estimator recovers alignment amplitudes for LRGs consistent with those found from higher-volume Quijote boxes at similar signal-to-noise ratio, despite the 20 × 20 \times 20 ×
There are several directions to which we could extend our work. Perhaps the most pressing direction would be to characterize the dependence of alignment coefficients on hydrodynamical prescription, but we find that the necessity of using L1_m8 hinders the study of alignments for the various FLAMINGO realizations. Of direct relevance to characterizing the impact of intrinsic alignments in weak lensing surveys is to extend our studies to multiple redshifts, as the time-dependence of the alignment signal is essentially unknown. Another pressing avenue is using paired measurements from halos and galaxies to understand whether a ‘halo occupation distribution-like’ construction can be made for shape galaxies [84 , 89 , 90 , 91 ] , from which high-fidelity halo alignment coefficient measurements could be used to build priors for galaxy alignments (as has been successfully done for scalar galaxy bias [92 ] ). Finally, we also note that one of the unique applications of galaxy shape correlations is the ability to probe anisotropic forms of primordial non-Gaussianity. The new operators induced in the shape bias expansion from these forms of PNG carry their own EFT coefficients, but field-level techniques can be used to measure these coefficients [93 , 94 ] . It is a straightforward extension of the formalism developed out in this work to develop an estimator for these PNG shape bias coefficients, conditioned on simulations which have anisotropic PNG [5 ] . It is our intention to return to many of these points in upcoming work.
Acknowledgments
We thank Beatriz Tucci, Noah Sailer, Martin White, Francisco Villaescusa-Navarro, Elisabeth Krause and Matias Zaldarriaga for helpful conversations. We thank Kazuyuki Akitsu for insightful conversations and code comparison. NK acknowledges support from a potpourri of funding sources including Edward and Kiyomi Baird, the Bershadsky Fund, and the Fund for Memberships in Natural Sciences of the Institute for Advanced Study. NK gratefully acknowledges the Pacific Postdoctoral Program at the Dark Universe Science Center, University of Washington, where part of this work was carried out. The Pacific Postdoctoral Program is supported by a grant from the Simons Foundation (SFI-MPS-T-Institutes-00012000, ML). This work was initiated at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452 and a grant from the Simons Foundation (1161654, Troyer). MMI is supported by the Office of High Energy Physics of U.S. Department of Energy under grant Contract Number DE-SC0012567. JMS acknowledges support provided by The Brinson Foundation through a Brinson Prize Fellowship grant. Calculations and figures in this work have been made using nbodykit [95 ] and the SciPy Stack [96 , 97 , 98 ] . This research has made use of NASA’s Astrophysics Data System and the arXiv preprint server. Large-language models provided by OpenAI under an educational license – of the GPT-5.6 and GPT-6 family – as well as a suite of Anthropic’s Claude models were used to aid in some of the analytic and numerical aspects of this work. These uses were independently cross-checked by the authors, who take full responsibility. No LLMs were used for the purposes of writing text.
First-principles transfer function and orthogonalization:
The starting point for the orthogonalization procedure in the scalar case is the construction of the operator correlation matrix (as discussed in [27 ] )
C i j ( k ) = ⟨ O i ( k ) O j ∗ ( k ) ⟩ ⟨ O i ( k ) O i ∗ ( k ) ⟩ ⟨ O j ( k ) O j ∗ ( k ) ⟩ (63) C_{ij}(k) = \frac{\langle\mathcal{O}_i(k)\mathcal{O}_j^*(k)\rangle}{\sqrt{\langle\mathcal{O}_i(k)\mathcal{O}_i^*(k)\rangle\langle\mathcal{O}_j(k)\mathcal{O}_j^*(k)\rangle}} \tag*{(63)} C ij ( k ) = ⟨ O i ( k ) O i ∗ ( k )⟩ ⟨ O j ( k ) O j ∗ ( k )⟩ ⟨ O i ( k ) O j ∗ ( k )⟩ (63) Let’s consider, generically, the extension of this to the helicity modes m m m
C i j m m ′ ( k ) = ⟨ O i m ( k ) O j m ′ ∗ ( k ) ⟩ ⟨ O i m ( k ) O i m ∗ ( k ) ⟩ ⟨ O j m ′ ( k ) O j m ′ ∗ ( k ) ⟩ (64) C_{ij}^{mm'}(k) = \frac{\langle\mathcal{O}_i^m(k)\mathcal{O}_j^{m'*}(k)\rangle}{\sqrt{\langle\mathcal{O}_i^m(k)\mathcal{O}_i^{m*}(k)\rangle\langle\mathcal{O}_j^{m'}(k)\mathcal{O}_j^{m'*}(k)\rangle}} \tag*{(64)} C ij m m ′ ( k ) = ⟨ O i m ( k ) O i m ∗ ( k )⟩ ⟨ O j m ′ ( k ) O j m ′ ∗ ( k )⟩ ⟨ O i m ( k ) O j m ′ ∗ ( k )⟩ (64) Since operators of different helicities in this basis are orthogonal, this matrix is diagonal in δ m m ′ K \delta^K_{mm'} δ m m ′ K
C i j m ( k ) = ⟨ O i ( k ) O j ∗ ( k ) ⟩ ⟨ O i ( k ) O i ∗ ( k ) ⟩ ⟨ O j ( k ) O j ∗ ( k ) ⟩ , (65) C^m_{ij}(k) = \frac{\langle\mathcal{O}_i(\mathbf{k})\mathcal{O}_j^*(\mathbf{k})\rangle}{\sqrt{\langle\mathcal{O}_i(\mathbf{k})\mathcal{O}_i^*(\mathbf{k})\rangle\langle\mathcal{O}_j(\mathbf{k})\mathcal{O}_j^*(\mathbf{k})\rangle}}, \tag*{(65)} C ij m ( k ) = ⟨ O i ( k ) O i ∗ ( k )⟩ ⟨ O j ( k ) O j ∗ ( k )⟩ ⟨ O i ( k ) O j ∗ ( k )⟩ , (65) where we’ve suppressed the helicity dependence of the operators since there’s no additional complexity. The orthogonalization procedure then follows in exact analogy to the scalar case. Due to the field-level degeneracies discussed in Section III C, for some of the lowest k k k C C C C ′ = C + λ I C' = C + \lambda I C ′ = C + λ I λ = 10 − 6 \lambda= 10^{-6} λ = 1 0 − 6 k k k m m m i , j i,j i , j
Validation of numerical Lagrangian basis and helicity decomposition
In this appendix we validate the numerical implementation of the shape bias basis, as well as the implementation of the helicity decomposition. Our numerical implementation on the lattice is based on the usual Lagrangian (non-galileon) basis, which we then rotate to the galileon representation before measuring field-level transfer functions.
In the absence of any displacements, the shape-shape power spectra of the Lagrangian bias basis assumes simple forms. The only non-zero spectra are the linear ( s a b , s c d ) (s_{ab},s_{cd}) ( s ab , s c d ) P 13 P_{13} P 13 s a b s_{ab} s ab t a b t_{ab} t ab L 3 a b L_{3ab} L 3 ab s a b s_{ab} s ab t a b t_{ab} t ab m = 0 m = 0 m = 0
Given a realization of the linear density field δ I C \delta_{\rm IC} δ IC [20 ] . For each helicity m m m spinosaurus. Strictly for the purpose of validation, we apply a spherical sharp-k k k ∣ k ∣ ≤ k U V = 1 h M p c − 1 |\mathbf{k}| \leq k_{\rm UV} = 1\,h\mathrm{Mpc}^{-1} ∣ k ∣ ≤ k UV = 1 h Mpc − 1 P ( k ) → erfc ( k − k U V Δ k ) P ( k ) / 2 P(k) \to\operatorname{erfc}\left(\frac{k-k_{\rm UV}}{\Delta k}\right)P(k)/2 P ( k ) → erfc ( Δ k k − k UV ) P ( k ) /2 Δ k = 0.05 h M p c − 1 \Delta k = 0.05\,h\mathrm{Mpc}^{-1} Δ k = 0.05 h Mpc − 1
The comparison between the lattice correlators and analytic predictions is shown in Figure 14 . We observe very close agreement for all helicities between Lagrangian correlators in the analytic and lattice realizations. Other checks such as orthogonality of fields for different m m m δ \delta δ
Figure 14. Comparison between analytic and lattice realizations of Lagrangian tensor correlators. The color indicates the helicity of the measurement, solid lines are the analytic prediction and the points are the measurements from ensemble of Quijote initial conditions. The final panel shows two cubic correlators cross-correlated with the linear density field, so there is only an m = 0 m = 0 m = 0
For the de-facto field-level calculations in this paper we do not implement this spherical sharp-k k k [99 ] code. A sharp-k k k k U V = 1 h M p c − 1 k_{\rm UV} = 1\,h\mathrm{Mpc}^{-1} k UV = 1 h Mpc − 1
Reconstructing FLAMINGO inertia tensors
The publicly available inertia tensors for FLAMINGO are incomplete. The public documentation states “When computing the inertia tensors, a filter was applied based on subhalo bound mass (2 × 10 11 M ⊙ 2 \times10^{11}M_\odot 2 × 1 0 11 M ⊙ 19 StellarInertiaTensorNoniterative). Selecting on subhalos with more than 300 star particles still results in 67.1% of subhalos not having a recorded inertia tensor. This poses a challenge in studying the stellar shapes of lower-mass halos, such as those that host ELGs or spiral galaxies.
To this end, we have re-implemented the halo shape calculation functionality of SOAP, the code used to compute shapes of halos in FLAMINGO, in order to recalculate missing shapes. Specifically, we aim to recompute the halo shapes associated with the BoundSubhalo attribute in SOAP catalogs, as opposed to shapes measured from various spherical aperture definitions. This is done by selecting halos with more than 20 star particles and recording their HaloCatalogueIndex. Given FLAMINGO particle membership files, for dark matter, baryon and star particles, we sub-select required particle IDs by selecting for particles whose GroupNr_bound are found in the set of halo catalog indices. From these matched particle catalogs we record the required variables to compute both κ c o r o t \kappa_{\rm corot} κ corot
From these sub-selected particle catalogs, we proceed to compute inertia tensors through SOAP’s built-in inertia tensor function20 Figure 15 we show the cumulative distribution functions of the absolute difference between the SOAP-stored inertia tensors and our independent reconstruction (for halos that have a non-zero tensor), normalized by the dispersion of each component in the original catalogs. We only show iterative, non-reduced tensors but the results are consistent across all definitions we have adopted. Our reconstruction of the original tensors are highly accurate, with ∼ 90 % \sim90\% ∼ 90% 1 σ 1\sigma 1 σ
Figure 15. Cumulative distribution functions for components of the difference between public FLAMINGO and our reconstructed shape tensors, for iterative dark matter and stellar halo shapes.
Part of the (very small) disagreements between halo shapes can be attributed to the compression algorithm applied to public SOAP catalogs from FLAMINGO. Compressing our reconstructed tensors with the same FMantissa9 filter resolves these discrepancies, but we opt to use the uncompressed reconstructed shape tensors in this work.
Analytic transfer function and parameter covariances
In the limit of a Gaussian-distributed stochasticity field ϵ m \epsilon_m ϵ m σ γ 2 / n ˉ \sigma_\gamma^2/\bar n σ γ 2 / n ˉ C i , m ( k ) C_{i,m}(k) C i , m ( k ) (47) . To see this, we follow the derivation of uncertainties in parameter covariances from Ref. [75 ] . Start from the estimator for the transfer function (we have suppressed the helicity index here, since all expressions are diagonal in m m m
C ^ i ( k ) = ∑ k γ ( k ) O i ⊥ ∗ ( k ) ∑ k ∣ O i ⊥ ( k ) ∣ 2 (66) \hat{C}_i(k) = \frac{\sum_k \gamma(k)\mathcal{O}_i^{\perp*}(k)}{\sum_k \left|\mathcal{O}_i^\perp(k)\right|^2} \tag*{(66)} C ^ i ( k ) = ∑ k O i ⊥ ( k ) 2 ∑ k γ ( k ) O i ⊥ ∗ ( k ) (66) where the k \mathbf{k} k γ ( k ) \gamma(\mathbf{k}) γ ( k )
C ^ i ( k ) = ∑ k , j [ C j ( k ) O j ⊥ ( k ) + ϵ ( k ) ] O i ⊥ ∗ ( k ) ∑ k ∣ O i ⊥ ∣ 2 ( k ) . (67) \hat{C}_i(k) = \frac{\sum_{\mathbf{k},j} \left[C_j(k) O_j^{\perp}(k) + \epsilon(k)\right] O_i^{\perp*}(k)}{\sum_{\mathbf{k}} \left|O_i^{\perp}\right|^2(k)}. \tag*{(67)} C ^ i ( k ) = ∑ k O i ⊥ 2 ( k ) ∑ k , j [ C j ( k ) O j ⊥ ( k ) + ϵ ( k ) ] O i ⊥ ∗ ( k ) . (67) Under the assumption that the transfer function does not vary drastically within the bin width (and since the transfer functions are generically smooth, this is reasonable), and that numerical violations of orthogonality are small (which they are to machine precision) the first term is simply C i ( k ) C_i(k) C i ( k )
σ 2 ( C ^ i ) ( k ) = ∑ k , k ′ O i ⊥ ( k ) O i ⊥ ∗ ( k ′ ) ⟨ ϵ ( k ) ϵ ( k ′ ) ⟩ ( ∑ k ∣ O i ⊥ ∣ 2 ( k ) ) 2 . (68) \sigma^2\left(\hat{C}_i\right)(k) = \frac{\sum_{\mathbf{k},\mathbf{k}'} O_i^{\perp}(k) O_i^{\perp*}(k') \left\langle\epsilon(k)\epsilon(k') \right\rangle}{\left(\sum_{\mathbf{k}} \left|O_i^{\perp}\right|^2(k)\right)^2}. \tag*{(68)} σ 2 ( C ^ i ) ( k ) = ( ∑ k O i ⊥ 2 ( k ) ) 2 ∑ k , k ′ O i ⊥ ( k ) O i ⊥ ∗ ( k ′ ) ⟨ ϵ ( k ) ϵ ( k ′ ) ⟩ . (68) Assuming a white error power spectrum, and using the relation ∑ k ∣ O i ⊥ ∣ 2 ( k ) = N m o d e s ( k ) ⟨ ∣ O i ⊥ ∣ 2 ⟩ ( k ) \sum_{\mathbf{k}} \left|O_i^{\perp}\right|^2(k) = N_{\mathrm{modes}}(k)\left\langle\left|O_i^{\perp}\right|^2 \right\rangle(k) ∑ k O i ⊥ 2 ( k ) = N modes ( k ) ⟨ O i ⊥ 2 ⟩ ( k ) (47) . The factor a m a_m a m ± m \pm m ± m Figure 16 we compare Equation (47) to empirical estimates of the scatter in transfer functions. The empirical measurements use transfer functions from N s i m s = 50 N_{\mathrm{sims}} = 50 N sims = 50 M h ∈ [ 10 12.5 , 10 13.0 ] M ⊙ h − 1 M_h \in[10^{12.5}, 10^{13.0}]\,M_\odot h^{-1} M h ∈ [ 1 0 12.5 , 1 0 13.0 ] M ⊙ h − 1 ∼ 20 % \sim20\% ∼ 20%
Figure 16. Ratio between the analytic (Equation (47) ) and empirical estimates for standard deviations of transfer functions. Each panel denotes a separate transfer function, while the colors denote the three different helicities. The empirical estimates are calculated using N = 50 N = 50 N = 50
Formulating the bias parameter estimation problem in the language of generalized least-squares also allows us to estimate the covariances of the IA shape biases analytically. For a generalized least squares problem, the parameter Fisher matrix is given
F = M T C − 1 M , (69) \mathbf{F} = \mathcal{M}^{T}\mathbf{C}^{-1}\mathcal{M}, \tag*{(69)} F = M T C − 1 M , (69) where M \mathcal{M} M (43) and C C C Figure 17 we show a comparison between the analytic, GLS-based parameter covariance and the empirical covariance from studying N s i m s = 50 N_{\mathrm{sims}} = 50 N sims = 50 M h ∈ [ 10 13.5 , 10 14.0 ] M ⊙ h − 1 M_h \in[10^{13.5}, 10^{14.0}]\,M_\odot h^{-1} M h ∈ [ 1 0 13.5 , 1 0 14.0 ] M ⊙ h − 1 amplitude of parameter error bars are in good agreement for this mass bin – the analytic estimator is within 20-25% of the empirical result for all IA parameters. In the middle and right panels of Figure 17 we show the correlation matrices between IA parameters for both estimates. While there are numerical differences, the qualitative agreement between both parameter covariances is strong. We can also see that the cubic parameters are relatively well-measured, although certain sets of cubic parameters are block-degenerate ({ c δ 2 s , c δ g s } \{c_{\delta^2s}, c_{\delta g_s}\} { c δ 2 s , c δ g s }
The results of this appendix justify the use of analytic expressions for the errors of the transfer function, as well as the analytic expression for the covariance of our field-level IA parameters, which will be used extensively when applied to Flamingo due to the lack of independent realizations.
Halo transfer functions
In this appendix we show the remaining, 10-operator transfer functions for all halo mass bins. Figures 18, 19, 20 and Figure 21 show transfer functions for the [ 10 12.5 , 10 13.0 ] [10^{12.5}, 10^{13.0}] [ 1 0 12.5 , 1 0 13.0 ] [ 10 13.0 , 10 13.5 ] [10^{13.0}, 10^{13.5}] [ 1 0 13.0 , 1 0 13.5 ] [ 10 14.0 , 10 14.5 ] [10^{14.0}, 10^{14.5}] [ 1 0 14.0 , 1 0 14.5 ] [ 10 14.5 , 10 15.0 ] [10^{14.5}, 10^{15.0}] [ 1 0 14.5 , 1 0 15.0 ] Figure 5 .
Figure 21. The same set of averaged transfer functions as shown in Figure 5 , but for the M ∈ [ 10 14.5 , 10 15.0 ] M ⊙ h − 1 M \in[10^{14.5}, 10^{15.0}]M_\odot h^{-1} M ∈ [ 1 0 14.5 , 1 0 15.0 ] M ⊙ h − 1
Converting c s c_s c s A 1 A_1 A 1
The intrinsic alignment literature conventionally reports as the linear amplitude a coefficient ‘A 1 A_1 A 1 c s c_s c s C ˉ 1 = 5 × 10 − 14 M ⊙ − 1 h − 2 M p c 3 \bar{C}_1 = 5 \times10^{-14}M_\odot^{-1}h^{-2}\mathrm{Mpc}^3 C ˉ 1 = 5 × 1 0 − 14 M ⊙ − 1 h − 2 Mpc 3 [48 ] , as well as the weak lensing community. The linear alignment amplitude A 1 A_1 A 1 [16 ] )
Figure 18. The same set of averaged transfer functions as shown in Figure 5 , but for the M ∈ [ 10 12.5 , 10 13.0 ] M ⊙ h − 1 M \in[10^{12.5}, 10^{13.0}]M_\odot h^{-1} M ∈ [ 1 0 12.5 , 1 0 13.0 ] M ⊙ h − 1
c s = − 2 A 1 C ˉ 1 ρ c r i t Ω m D ( z ) − 1 . (70) c_s=-2A_1\bar{C}_1\rho_{\mathrm{crit}}\Omega_mD(z)^{-1}. \tag*{(70)} c s = − 2 A 1 C ˉ 1 ρ crit Ω m D ( z ) − 1 . (70) Since ρ c r i t = 2.78 × 10 11 h 2 M ⊙ M p c − 3 \rho_{\mathrm{crit}}=2.78\times10^{11}h^2M_\odot\mathrm{Mpc}^{-3} ρ crit = 2.78 × 1 0 11 h 2 M ⊙ Mpc − 3 C ˉ 1 ρ c r i t ≈ 0.014 \bar{C}_1\rho_{\mathrm{crit}}\approx0.014 C ˉ 1 ρ crit ≈ 0.014 Ω m = 0.306 \Omega_m=0.306 Ω m = 0.306 A 1 = 5 A_1=5 A 1 = 5
c s ≈ − 4.3 × 10 − 2 D ( z ) ( Ω m 0.306 ) ( A 1 5 ) . (71) c_s\approx-\frac{4.3\times10^{-2}}{D(z)}\left(\frac{\Omega_m}{0.306}\right)\left(\frac{A_1}{5}\right). \tag*{(71)} c s ≈ − D ( z ) 4.3 × 1 0 − 2 ( 0.306 Ω m ) ( 5 A 1 ) . (71) As mentioned in the main text, we take A 1 ≈ 5 ± 0.2 A_1\approx5\pm0.2 A 1 ≈ 5 ± 0.2 z = 0 z=0 z = 0 [48 ] , as approximate uncertainties from fitting to P g g P_{gg} P g g P g I P_{gI} P g I ( 2.8 G p c ) 3 (2.8\mathrm{Gpc})^3 ( 2.8 Gpc ) 3 [47 ] they also provide a fitting function for the redshift evolution of linear alignments as
A 1 ( M , z ) = α M ( z ) ( M 10 13.5 M ⊙ ) β M ( z ) , (72) A_1(M,z)=\alpha_M(z)\left(\frac{M}{10^{13.5}M_\odot}\right)^{\beta_M(z)}, \tag*{(72)} A 1 ( M , z ) = α M ( z ) ( 1 0 13.5 M ⊙ M ) β M ( z ) , (72) with α M ( z ) = − 2.18 z 2 + 12.69 z + 6.64 \alpha_M(z)=-2.18z^2+12.69z+6.64 α M ( z ) = − 2.18 z 2 + 12.69 z + 6.64 β M ( z ) = 0.26 z 2 − 0.5 z + 0.88 \beta_M(z)=0.26z^2-0.5z+0.88 β M ( z ) = 0.26 z 2 − 0.5 z + 0.88 A 1 ( M ∗ , z = 0 ) = 5 A_1(M^\ast,z=0)=5 A 1 ( M ∗ , z = 0 ) = 5 M ∗ = 2.3 × 10 13 M^\ast=2.3\times10^{13} M ∗ = 2.3 × 1 0 13 A 1 ( M ∗ , z = 0.5 ) = 9.9 A_1(M^\ast,z=0.5)=9.9 A 1 ( M ∗ , z = 0.5 ) = 9.9 D ( z = 0.5 ) = 0.77 D(z=0.5)=0.77 D ( z = 0.5 ) = 0.77 100 c s ≈ − 11.05 ± 0.44 100c_s\approx-11.05\pm0.44 100 c s ≈ − 11.05 ± 0.44
To convert the measured c G s c_{Gs} c G s c δ s c_{\delta s} c δ s
Figure 19. The same set of averaged transfer functions as shown in Figure 5, but for the M ∈ [ 10 13.0 , 10 13.5 ] M ⊙ h − 1 M \in[10^{13.0},10^{13.5}]M_\odot h^{-1} M ∈ [ 1 0 13.0 , 1 0 13.5 ] M ⊙ h − 1
s ~ a b = s a b + δ s a b − Ψ ( 1 ) ⋅ ∇ s a b . (73) \widetilde{s}_{ab}=s_{ab}+\delta s_{ab}-\boldsymbol{\Psi}^{(1)}\cdot\nabla s_{ab}. \tag*{(73)} s ab = s ab + δ s ab − Ψ ( 1 ) ⋅ ∇ s ab . (73) Conversion to Eulerian amplitudes follows through writing Ψ ( 1 ) ⋅ ∇ s a b \boldsymbol{\Psi}^{(1)}\cdot\nabla s_{ab} Ψ ( 1 ) ⋅ ∇ s ab
s a b ( 2 ) = D a b δ ( 2 ) = s a b 2 + 2 3 δ s a b − 5 2 t a b − Ψ ( 1 ) ⋅ ∇ s a b , s_{ab}^{(2)}=\mathcal{D}_{ab}\delta^{(2)}=s_{ab}^{2}+\frac{2}{3}\delta s_{ab}-\frac{5}{2}t_{ab}-\boldsymbol{\Psi}^{(1)}\cdot\nabla s_{ab}, s ab ( 2 ) = D ab δ ( 2 ) = s ab 2 + 3 2 δ s ab − 2 5 t ab − Ψ ( 1 ) ⋅ ∇ s ab , which was originally derived in [55 ] . Collecting terms and equating the shifted-operator expansion to the Eulerian expansion at second order gives
s ~ a b = s a b ( E ) − G s s a b + 5 2 t a b . \widetilde{s}_{ab}=s_{ab}^{(E)}-G_{s}^{s}{}_{ab}+\frac{5}{2}t_{ab}. s ab = s ab ( E ) − G s s ab + 2 5 t ab . We can read off that the Lagrangian coefficients of the shape bias expansion are re-mapped to
c s E = c s L , c δ s E = c δ s L , c_{s}^{E}=c_{s}^{L},\qquad c_{\delta s}^{E}=c_{\delta s}^{L}, c s E = c s L , c δ s E = c δ s L , c G s E = c G s L − c s L , c t E = c t L + 5 2 c s L , c_{G s}^{E}=c_{G s}^{L}-c_{s}^{L},\qquad c_{t}^{E}=c_{t}^{L}+\frac{5}{2}c_{s}^{L}, c G s E = c G s L − c s L , c t E = c t L + 2 5 c s L , but note the remapping is different for the non-Galileon basis. The usual Eulerian TATT amplitudes A 1 δ , A 2 A_{1\delta}, A_{2} A 1 δ , A 2 [22 ]
C δ s = − A 1 δ C ~ 1 ρ c r i t Ω m D ( z ) − 1 (74) C_{\delta s}=-A_{1\delta}\widetilde{C}_{1}\rho_{\mathrm{crit}}\Omega_{m}D(z)^{-1} \tag*{(74)} C δ s = − A 1 δ C 1 ρ crit Ω m D ( z ) − 1 (74) C 2 = A 2 ( 5 C ~ 1 ρ c r i t Ω m , f i d ) Ω m 2 D ( z ) − 2 (75) C_{2}=A_{2}\left(\frac{5\widetilde{C}_{1}\rho_{\mathrm{crit}}}{\Omega_{m,\mathrm{fid}}}\right)\Omega_{m}^{2}D(z)^{-2} \tag*{(75)} C 2 = A 2 ( Ω m , fid 5 C 1 ρ crit ) Ω m 2 D ( z ) − 2 (75) Figure 20. The same set of averaged transfer functions as shown in Figure 5 , but for the M ∈ [ 10 14.0 , 10 14.5 ] M ⊙ h − 1 M \in[10^{14.0},10^{14.5}]M_\odot h^{-1} M ∈ [ 1 0 14.0 , 1 0 14.5 ] M ⊙ h − 1
These Eulerian ‘C’ coefficients are defined with respect to the projected fields, so there is a similar C i = c i / 2 C_i=c_i/2 C i = c i /2
A 1 δ ≈ − c δ s ( E ) D ( z ) 8.5 × 10 − 3 ( 0.306 Ω m ) (76) A_{1\delta} \approx-\frac{c_{\delta s}^{(E)}D(z)}{8.5 \times10^{-3}}\left(\frac{0.306}{\Omega_m}\right) \tag*{(76)} A 1 δ ≈ − 8.5 × 1 0 − 3 c δ s ( E ) D ( z ) ( Ω m 0.306 ) (76) ≈ − D ( z ) ( c δ s ( G ) − c G s / 3 + c s / 3 8.5 × 10 − 3 ) ( 0.306 Ω m ) , (77) \approx-D(z)\left(\frac{c_{\delta s}^{(G)}-c_{Gs}/3+c_s/3}{8.5 \times10^{-3}}\right)\left(\frac{0.306}{\Omega_m}\right), \tag*{(77)} ≈ − D ( z ) ( 8.5 × 1 0 − 3 c δ s ( G ) − c G s /3 + c s /3 ) ( Ω m 0.306 ) , (77) where c δ s ( E ) c_{\delta s}^{(E)} c δ s ( E ) δ s \delta s δ s A 2 A_2 A 2
A 2 ≈ D ( z ) 2 ( c G s − c s 4.25 × 10 − 2 ) ( 0.306 Ω m ) . (78) A_2 \approx D(z)^2\left(\frac{c_{Gs}-c_s}{4.25 \times10^{-2}}\right)\left(\frac{0.306}{\Omega_m}\right). \tag*{(78)} A 2 ≈ D ( z ) 2 ( 4.25 × 1 0 − 2 c G s − c s ) ( Ω m 0.306 ) . (78) 
