Galaxy shape statistics in the effective field theory

Intrinsic galaxy alignments yield an important contribution to the observed statistics of galaxy shapes. The general bias expansion for galaxy sizes and shapes in three dimensions has been recently described by Vlah, Chisari & Schmidt using the general perturbative effective field theory (EFT) framework, in analogy to the clustering of galaxies. In this work, we present a formalism that uses the properties of spherical tensors to project galaxy shapes onto the observed sky in the flat-sky approximation, and compute the two-point functions at next-to-leading order as well as the leading-order three-point functions of galaxy shapes and number counts. The resulting expressions are given in forms that are convenient for efficient numerical implementation. For a source redshift distribution typical of Stage IV surveys, we find that nonlinear intrinsic alignment contributions to galaxy shape correlations become relevant at angular wavenumbers l ≳ 100.


Introduction
The advent of wide and fast galaxy surveys has triggered an era of precision cosmology from the large-scale structure, with this cosmological tool being used to elucidate the nature of the dark components of the Universe and the physics of the early Universe, among other goals. But in order for the inferences drawn from ever more powerful data sets to be robust, a similar progression must occur on the modelling side: our model for the observables needs to be accurate enough to guarantee no residual biases in the analyses.
Weak gravitational lensing is one of the phenomena observed in the large-scale structure which can yield information about the nature of dark matter and the origin of the accelerated expansion of the Universe. Gravitational lensing by the large-scale structure of the Universe distorts the shapes of background sources (galaxies) in a correlated manner. Although the signal is challenging to extract (amounting to percent-level correlated changes in galaxy ellipticities), it is now regularly and successfully extracted from multiple data sets, and a number of collaborations have presented cosmological constraints from this probe in the last decade [1][2][3][4][5][6][7].
One of the challenges in the modelling of galaxy shape correlations is the coherent distortion induced by gravitational interactions between galaxies and the large-scale distribution of matter in the Universe. Such physical distortions are known as "intrinsic alignments" and have been detected to high significance by multiple experiments, with an increasing body of work focusing on which, how and when galaxies align, and how these alignments give rise to a contamination of the weak lensing signal. On the other hand, their detection is on such firm grounds [8][9][10][11] that there are proposals for their use as a cosmological probe on their own right [12][13][14][15].

JCAP05(2021)061
So far, the modelling of intrinsic galaxy alignments has lagged behind similar efforts for other large-scale structure probes. The most commonly adopted model is the "linear alignment model" (LA) [16], which postulates a linear relation between the projected shapes of galaxies and the projected tidal field of the large-scale structure. This ansatz in any case provides a good description at large separations only. In the nonlinear regime, an excess above LA predictions is observed for early-type galaxies, which can be well-approximated by an ad hoc replacement of the linear matter power spectrum by its nonlinear analogue, while preserving the assumption of linear alignment [17]. An alternative for nonlinear modelling is the so-called "halo model", by which the small-scale regime is modelled through the distribution of aligned central and satellite galaxies within a halo [18][19][20]. On the other hand, disk galaxies are assumed to be subject to alignments via tidal torques, and a separate set of models has been constructed to describe this flavor of alignments [21][22][23]. Cosmological simulations serve as a validation tool for all these models and help constrain their free parameters in as much as the simulated galaxy population is representative of the observed one [24][25][26][27][28][29][30][31][32][33][34][35][36][37].
The increasing precision of galaxy shape measurements requires, however, new modelling tools if we are to ensure the success of weak lensing cosmology [38]. Recently, a new theory approach was put forward that unifies the two alignment mechanisms under the same framework by means of standard perturbation theory (SPT) [39][40][41]. This more rigorous approach expands over the existing models and has already been applied to observational data sets [5,42]. In [43], we presented a complete formulation that is based on a similar expansion in the context of the "effective field theory" (EFT) of the large-scale structure. Considering the symmetries of a trace-free tensor, we were able to identify all potential gravitational observables at a given order that contribute to describing any given biased tensorial tracer of the large-scale structure, including the intrinsic shapes of galaxies.
The EFT framework we presented in [43] focused on three-dimensional intrinsic shapes and their correlations. In this work, we connect the three-dimensional observables computed in that work to the projected ones from galaxy surveys, namely: quantities corresponding to the convergence field κ and shear fields γ 1 and γ 2 . We present results for the one-loop angular power spectrum and tree-level angular bispectrum for the convergence and shear auto-and cross-correlations. We emphasise, however, that the formalism derived in this work does not rely explicitly on the EFT expansion. Our findings can be applied to projecting the statistics of any three-dimensional scalar and tensorial fields regardless of how one chooses to model them. In this paper, we do not include nonlinear projection contributions, such as those from reduced shear. We will argue below that these are much smaller than the nonlinearities in intrinsic shape correlations which we treat here.
This work is organized as follows. In section 2, we summarise the formalism presented in Paper I [43]. We define the projected observables in section 3. We construct predictions for intrinsic alignment two-point and three-point correlations assuming the flat-sky approximation in section 4. We produce numerical predictions for the typical observables and redshifts of galaxy surveys in section 5. Conclusions are presented in section 6. The contribution of gravitational lensing to the two-point shape statistics are included assuming the flat-sky approximation in appendix A. Details of the projection calculation for the angular power spectrum are presented in appendix B and similarly for the angular bispectrum in appendix C.
Throughout the paper, we assume an Euclidean ΛCDM cosmology with Ω m = 0.272, Ω b = 0.0455, n s = 0.967, σ 8 = 0.807 and h = 0.704. We work under the assumptions of adiabatic Gaussian perturbations and General Relativity. It is straightforward to include Dirac delta function Shear derivative k 1···n ≡ k 1 + · · · + k n Sum notation Y (m) ij (k) Rank-two spherical tensor basis in Fourier space ij (r) Rank-two spherical tensor basis in configuration space the impact of primordial non-Gaussianity following previous work [13]. In perturbative calculations, we always make the usual approximation of setting the n-th order growth factor D (n) (τ ) to [D(τ )] n .

Galaxy shapes and bias expansion
The three-dimensional shapes of galaxies in their rest frame are spin-2 quantities that can be described in terms of a well-defined expansion in local gravitational observables. In Paper I [43], we assumed that, for each galaxy α, its light distribution is described by the symmetric where the Dirac delta function ensures the total momentum conservation, and is a consequence of the statistical translation invariance.

JCAP05(2021)061
second-moment tensor I ij (r α ) (i.e., an ellipsoid) and defined by the shape fluctuation field 2 where the trace-free tensor g ij describes galaxy shape perturbations, and the scalar field δ s , galaxy size perturbations. We presented a new EFT expansion for the former in Paper I [43], while the latter is described by the usual EFT of scalar biased tracers [44,45] in complete analogy to the case of galaxy number counts. This formalism allowed us to provide expressions for the two-point correlations between the intrinsic galaxy shape field and other scalar fields up to next-to-leading order, for which we drew analogies with the EFT application to scalar quantities, such as the number density of biased tracers [44,45]. Moreover, we presented three-point functions for these tracers at leading order. The EFT approach also shows that one has to allow for distinct bias parameters in the bias expansion for the trace (size) and the trace-free (shape) parts. Explicitly, we expand the number counts δ n , size fields δ s , and shape fields g ij in terms of local gravitational observables as follows where tr and TF denote the trace and trace-free components of a tensor, respectively, and O are the corresponding (renormalized) bias parameters. The expansion is not unique, but if the physical processes that determine the properties of the tracers are local, reference [44] showed that one possible complete basis is the one comprised of all scalar combinations of a set of operators Π [n] , defined recursively starting from Note that eq. (2.3) contains the leading gravitational observables at a given spacetime position r, τ : namely, the matter density perturbation δ m , and the scaled tidal field D ij δ m . In the case of shapes, the expansion must account for all possible trace-free tensor combinations, which in general results in more contributions at a given order than in the scalar (density) case. The computation of correlations between biased tensorial and scalar fields can be simplified by considerations of isotropy. To this end, in Paper I [43] we proposed to decompose any given tensorial field in spherical tensors [46] of multipole order = 2, whose transformation properties under rotation are known. We thus decompose the shape tensor field into these basis functions in terms of helicity m = 0, ±1, ±2 defined with respect to the wavevector k: 2 Different normalizations of the Sij field are possible, though this does not affect the form or range of validity of the EFT expansion.
and Y (0) = 1 is the single scalar helicity-0 mode (explicitly written in this form for the purpose of symmetry). Hereafter, we will drop the tensor indices on the basis vectors and use the boldface symbol Y in order to emphasize that we are dealing with a basis = 2 tensor, while for the = 0 component we will use the explicit form given in eq. (2.2). Note that the spherical tensor components S (0) 0 (k) and S (m) 2 (k) can be directly obtained from S ij by projections using the basis tensors.
To construct the basis explicitly, we start by defining an orthonormal basis of threedimensional Euclidean space comprised of tensors wheren can be chosen as the line of sight direction. With this basis defined, we are able to build our helicity basis, As a consequence, the orthonormal tensor basis functions can be defined as where the normalization is chosen as are trace-free, and complex , while a similar relation holds for the basis {e ± , e 0 }.

Two-point functions in 3D
In Paper I [43], we considered the auto-correlations of scalar and trace-free tensor fields, and their cross-correlation, defined as and we obtained explicit expressions for these by using the EFT expansion of the shape field (eq. (2.2)) and projecting out the trace and trace-free part of the power spectra. To simplify the calculation, we made use of the spherical tensor basis given in eq. (2.4) and obtained these power spectra by combination of the ones corresponding to the spherical tensor field components S (m) 2 , which are invariant quantities that also transform as spherical tensors. For projecting onto the flat sky (section 4), we can directly operate with the power spectra of the S (m) 2 components, which we proceed to summarize in what follows. Due to statistical isotropy and homogeneity, the component power spectra are given by where , = 0, 2 represent the trace and trace-free components, respectively, and m, m are helicity components. Different helicity components do not correlate. This is a consequence of statistical isotropy. Because the shape tensor field is real, and we assume parity invariance, this implies that the power spectra of the helicity components must satisfy: . As a consequence, the shape-shape tensor power spectrum can be decomposed into six independent contributions, where the last term is symmetrized over (q) and (−q) contributions, and the curly brackets in the subscripts denote symmetrization over the two pairs of indices.
We are now able to take the trace or trace-free components to obtain expressions for the power spectra defined in eq. (2.8), for which we obtain Explicit expressions for the spherical component power spectra up to one-loop order in perturbation theory were provided in section 5 of Paper I [43] relying on the EFT expansion. Nevertheless, the decomposition in spherical tensors that we presented in Paper I [43] and summarized here is valid nonlinearly for any order in PT.

Three-point functions in 3D
We now consider three tensor fields, each with a trace and trace-free part. The bispectrum of these fields is defined by Given the decomposition of S ij into spherical tensor components (eq. (2.4)), we can also define the bispectrum of a combination of the components as As a result, the decomposition of the bispectrum in terms of the Y (m) basis tensors is then

JCAP05(2021)061
with 2 cyclic permutations as indicated. Statistical rotation and parity invariance give additional requirements that these bispectra need to satisfy, namely This significantly reduces the number of the independent bispectra of spherical tensor components in the expansion of eq. (2.14).
The EFT framework developed in Paper I [43] allowed us to find all the independent contributions to each bispectrum B (m 1 ,m 2 ,m 3 ) 1 2 3 . We refer the reader to that work for the explicit expansion of the three-dimensional bispectra of galaxy densities, sizes and shapes. Here, we will show how these bispectra, generically eq. (2.12), can be projected on the sky with the help of the spherical tensor decomposition of eq. (2.14).
The observables we are interested in are auto-and cross-correlators of galaxy densities, sizes and shapes. As mentioned in Paper I [43], one should note that helicity ±1 and ±2 can in principle contribute to any of the bispectra where at least one shape field is correlated. This is due to symmetry constraints being less stringent in the bispectrum, when compared to the power spectrum.

Statistics of shape fields projected on the sky
In imaging surveys we usually do not have access to the 3D shape of galaxies that we have thus far been describing. We can measure only 2D images, obtained by projections of the 3D shapes on the sky plane. Thus, in order to connect our theoretical models with the measured shape γ ij these sky-projections need to be taken into account. If we introduce the projector wherer is the direction pointing along the line of sight on the sky, we can write the total shape field on the sky as where γ G,ij denotes the weak lensing shear, and r s and r are respectively the redshift space and real space coordinates. 3 At linear order lensing shear can be simply added to the intrinsic shape, however this holds only at leading order. In general, the observed shape field is a nonlinear function of shear and intrinsic shape; we will discuss this issue below. We will not consider these nonlinear effects, but defer them to future work. Considering thus only the projections of the intrinsic shape, and neglecting redshiftspace distortions, we can define where TF here stands for the trace-free components of the projected, 2D part, of the 3D tensor and consequently we can define the total projecting operator 3 Note that the projectors Pij are invariant on the redshift space mapping, i.e.
since the directions of real and redshift space coordinates coincide.

JCAP05(2021)061
Projector operators have a simple idempotent property P ij P jk = P ik and are orthogonal tô r irj , but not to the isotropic term δ K ik . To obtain the Fourier transform of γ I,ij (r, z) field the following integral would have to be performed The mode coupling evidenced in the last equality is due to the fact that the projection operator breaks translation invariance. This issue is circumvented by performing the spherical harmonic decomposition on the sky, or the 2D Fourier transform in case when sufficiently small scales are considered (flat-sky approximation). However, let us first decompose the configuration space tensor field γ I,ij into the irreducible spherical tensor components, similar to what we did in Fourier space in section 2.1. We can thus introduce the helicity basis (r, m + , m − ) in analogy to eq. (2.5), where e x is a fixed unit vector that is not parallel tor, and construct out of it the rank two harmonic tensor basis M (defined in analogy to eq. (2.7)) these would be orthogonal to the total projection operator P ijkl given in eq. (3.4) and thus would not contribute to the decomposition of the projected shape γ I,ij . The traceless intrinsic shape field can then be decomposed in terms of two components on the sky An especially useful property of this basis is the realization that m ± basis vectors are unit eigenfunctions of the projectors P ij , i.e. we have which combined with eq. (3.3) gives the simple relation The components γ ±2 transform under rotation as helicity-two functions, and we have γ * If we integrate over the line of sight coordinate, introducing a window function W (χ) determined by the redshift distribution dN (z)/dz of galaxies, we have (3.11) where the variable χ is the comoving distance to redshift z. Figure 1. Coordinate setup in the flat sky approximation for the angular two-point functions. We use the coordinates labels where r = r n + κ, with r =n.r, and κ is the coordinate that lies in the observation plane, where we introduce the angular variable θ = κ/χ (χ is the comoving distance at redshift z). The observer is assumed to be at the point O, and the observation plane is assumed to be at the distance determined by redshift z.
In surveys that do not cover a wide area of the sky, it is often useful and simpler to work in the flat-sky approximation. Here, the area on the sky is assumed to be well approximated by a plane lying perpendicular to a fixed directionn. This allows us to setup the coordinate frame where r = r n + κ, with r =n.r, κ = r − r n, and θ = κ/χ. This coordinate setup is schematically shown in figure 1. In this setting the integral over the redshift is performed along the line of sightn and we have (3.12) and the spherical tensors M (±2) ij (n) are also defined relative to then direction. Since we have established a 2D planar coordinate θ, we can also introduce a corresponding 2D Fourier transform. This is again useful, in the flat sky approximation, since the statistical isotropy can serve as an homogeneity condition on the plane and we can thus use statistical translation invariance. We thus define a Fourier field on a plane (3.13) and the decomposition in spherical tensors defined in the plane orthogonal to the line of sightr is (3.14) The basis M (±2) spans a 2D Fourier space while orthogonal to then direction. We can think of it as being constructed from vector basis (n,ˆ + ,ˆ − ), analogously as we did in eq. (2.5), but withˆ ± referring to a fixed directionn. Each of the helicity components in this Fourier space is given by Note that it is not obvious that γ ±2 ( ), as defined above, will correspond to the small-angle limit of all sky treatment for high enough . However, this correspondence has been explicitly established in e.g. [47], by taking the small-angle limit of the all sky results, and we will return to it in our upcoming all-sky paper. Very frequently, different bases for decomposing the γ I,ij are used in the literature. Some of the most prominent are the polarization basis (E, B) used in the study of the CMB polarization [48][49][50], and what we label the Pauli basis (γ 1 , γ 2 ). The relation of our helicity basis to the polarization basis of galaxy shapes (γ E , γ B ) is a simple linear combination Since we deal with the simple linear combination of the fields, these relations hold in both all-sky and flat-sky, as well as in both angle and Fourier space of the flat-sky approximation. On the other hand Pauli basis can be obtained by similar linear combination and additional rotation around the azimuthal angle φ. We have 17) and the inverse is given by γ ±2 = e ±2iφ (γ 1 ± iγ 2 ), where φ = φ θ in real space and φ = φ in Fourier space, denoting the azimuthal angles of θ and , respectively. The Pauli basis can again be used either in all-sky or flat-sky, however it is more commonly used in the flat-sky case. As the name suggests, on the 2D sky in this basis can be organized as the linear combination where σ 1 and σ 3 are the Pauli matrices (see e.g. [51]). The relation of the Pauli basis and polarization basis is related by the simple plane rotation by a 2φ θ spin angle. Electric (E) and magnetic (B) components in angle space are then given by and equivalent expressions, as already noted, hold in the 2D Fourier space.

Angular power spectra
In this part, we explore the two-and three-point statistics in the flat-sky approximation. The coordinate set-up is as described in the earlier section and as suggested for the twopoint function in figure 1. These results agree with the corresponding all-sky result in the small-angle (high ) limit, as will be shown in an upcoming paper.

Two-point functions
We can evaluate the angular correlation function using eq. (3.15). We are interested in the following correlators: We use the s variable to label the helicity states ±2, and we will refer tos as the sign of the helicity state, i.e.s = s/|s|. Here and throughout, we assume galaxy number counts δ n as the scalar observable we are correlating with. The expressions however are applicable generally to any such scalar (e.g., galaxy sizes δ s ). We note again that in the flat-sky approximation statistical isotropy manifests itself as the translation invariance (and rotation invariance around the line of sight) of the correlators on the plane. This allows us to define simple spectral functions that depend only on the amplitude of the -mode on the plane The explicit calculation of these angular power spectra is given in appendix B. We use the standard approximation that the longitudinal modes k =n.k do not contribute noticeably to the angular spectra compared to the modes perpendicular to the plane, i.e. k k ⊥ . Note that this approximation is not strictly necessary, it however becomes accurate in the same regime where the flat-sky approximation does, and greatly simplifies the resulting expressions. Using the resulting expressions given in eq. (B.15) we obtain where W n and W g denote the window functions corresponding to number counts and shapes, respectively, and we have the normalization constants N 0,1,2 = 3 2 , 1 2 , 1 (related to the basis Y (m) ij in eq. (2.7)). Above, we allowed for the fact that the window functions for number density W n and for shapes W g can be different. Note that when using these flat-sky expressions, it is found to be a much better approximation to use → +1/2 in the arguments of P q ss ( /χ) power spectra above (see e.g. [52]). With this correction, the estimated error of the approximation is reduced from O( −1 ) to O( −2 ).

JCAP05(2021)061
We can also consider the angular spectra in different bases, introduced earlier. In the polarization basis, defined in eq. (3.16), angular power spectra are The cross correlators of B-modes with E-modes or a scalar field vanish as a consequence of statistical parity invariance (see e.g. Paper I [43]): C nB ( ) = C EB ( ) = 0. The angular power spectra in the helicity basis γ ± and the polarization basis are related by The Pauli basis is related to the polarization by simple rotation in the plane eq. (3.19). Angular power spectra in this basis are So far, all expressions included intrinsic contributions to shape correlations. The corresponding expressions including the leading lensing effect are given in appendix A.
Let us now return to the issue of nonlinear projection contributions. Since the shape measurement process is nonlinear, and the galaxy shape field is weighted by the number density of source galaxies, eq. (3.2) receives additional corrections. One such effect is that shapes are in fact measuring the reduced shear [53][54][55]. However there are also nonlinear effects which depend on the intrinsic shapes of galaxies [56][57][58]. Some examples of resulting contributions are where γ G denotes the lensing contribution to the shape distortion. On the other hand, in terms of purely intrinsic contributions, the effect of number-density weighting is already included in our shape bias expansion, as explained in Paper I [43]. The contributions in eq. (4.7) schematically add terms of the following form to angular shear correlations: That is, the convolution is performed after the projection. This is in contrast to the 1-loop shape contributions in eq. (4.4), which involve the projection of a 3D convolution integral.

JCAP05(2021)061
As shown in detail in [59], the latter contributions are much larger than the contributions from nonlinear projection effects. For this reason, we are justified in ignoring the latter here, even though they are formally of the same order in perturbations.
In terms of lensing contributions only, this hierarchy is well known: reduced-shear corrections to the shear power spectrum are much smaller than the corrections to the nonlinearities in the matter density field [55,57,58].

Three-point functions
We turn next to the angular bispectrum calculation in the flat-sky approximation. The derivation details are presented in appendix C, and here we give a summary of the results. For some notable earlier work that used the perturbation theory to obtain the projected bispectrum statistics and taking intrinsic alignments into account, see e.g. [53,[60][61][62][63].
In analogy to how we define the flat-sky angular power spectrum in eq. (4.9) we can define the angular bispectrum. We thus have To obtain the relation of these angular bispectra to the 3D spherical tensor bispectra in section 2.3 in the flat-sky approximation we again use the approximation that the longitudinal modes k =n.k contribute negligibly to the perpendicular modes k ⊥ that lay in the plane of the sky (see figure 1). We thus assume that k k ⊥ . We again note that relaxing this approximation would give rise to the subleading corrections in variables. However, with this approximation we ensure translation invariance in the plane manifested as the Dirac delta's on the left hand side of definitions in eq. (4.9). We obtain the results given in eq. (C.6) that we review here. The density-density-density angular bispectrum is given by where we use the shorthand notation˜ = χ. This result is the standard clustering bispectrum [53]. The density-density-shape angular bispectrum is given by (4.11) In the polarization basis introduced in eq. (3.16), density-density-shape angular bispectrum can be split into E and B components yielding This decomposition has already been shown in Paper I [43] (see also [41] for similar results).

JCAP05(2021)061
The density-shape-shape angular bispectrum in helicity basis can compactly be written as (1,0) 022 It is important to notice that if we consider the cross-bispectrum (correlating at least two different shape tracers) we have four independent angular bispectra, since B n−+ can be different from B n+− . In the case when only auto-bispectra (same tracers) are considered this simplifies and we have B n−+ = B n+− . This can be seen if we consider the difference of the two angular bispectra. Using the explicit form above we have (4.14) , which vanishes if the auto-correlation bispectrum is considered. This situation is intrinsically different from the case of power spectrum where C +− = C −+ (see eq. (4.3)) was guaranteed even for the cross correlations of different shape tracers. The reason for this lies in the fact that the imposed symmetries (statistical isotropy and parity invariance) are more constraining for the two-point functions, leaving less functional freedom than in case of three-point functions.
We can transform these results into the E and B basis to get We again note that B nEB and B nBE are equal in the case of tracer auto-correlations, while in general they can be different when cross-correlating different tracers.

JCAP05(2021)061
Finally, we turn to the correlations of the three shape fields. The general expression in terms of the helicity basis is given by As in case of the density shape-shape bispectra, we see that in the case where we are considering the single tracer auto-correlations, we have B ++− = B +−+ = B −++ , and similarly for the permutations of B +−− . Transforming these results into the E and B basis, we get Where, to obtain the values, for other cross-corrlation components, like B EBE etc., we need to take permutations of the q = 1 helicity index in B (q 1 ,q 2 ,q 3 ) 222 bispectra in the expressions above. In case of auto correlations these are all equivalent again.
In Paper I [43], we presented the results for the 3D bispectra B at tree level (LO -leading order) PT. We also showed explicitly the results for the density-densityshape angular bispectra, demonstrating that B nnE both B nnB contributions arise already at leading PT order (see also [41]). For density-shape-shape and shape-shape-shape angular -15 -JCAP05(2021)061 bispectra, the situation is different. At tree level, only one nonzero helicity q can contribute. This is because at leading order in PT, the Π ij (r) field given in eq. (2.3) has only scalar contributions. The contributions that survive at leading order for all bispectra are then listed in eq. (C.7) of appendix C. We see that the B nnE and B nnB expressions given in eq. (4.12) remain unchanged, while for the rest the number of contributions is substantially reduced. The expressions given in eq. (4.15) at tree level thus reduce to where B nBE is again equal to B nEB in case of auto-correlations, while in case of crosscorrelations in the expression above we need to replace the B  , while the rest of the angular bispectra, involving more then one γ B field are zero at tree-level.

Results for one-loop power spectrum and tree-level bispectrum
In this section we show the angular power spectrum and bispectrum of intrinsic galaxy shapes, based on the results obtained in Paper I [43]. As we indicated in the previous section we show the results using the flat-sky approximation. The primary motivation for this lies in the realization that on the scales where one-loop terms are relevant, flat-sky is a good approximation to the all-sky result. We will address the return to the all-sky results in the subsequent paper.
Before we continue, we would like to point out the issue of the redshift dependence of the free bias parameters. Generally, the bias parameters are expected to evolve on a Hubble time scale, so that one can write If the redshift distribution considered is not too wide, then one can approximately work with constant bias parameters (the correction is of order (∆z) 2 where ∆z is the width of the redshift distribution). We will however not explore this issue further in this paper and we thus assume that bias parameters are independent of redshift. Furthermore, in order to show some concrete results we assume a simple weight function W g (χ) that corresponds to the redshift distribution of galaxies with shapes expected for a next generation wide and deep survey. In other words, e.g. [38], where α = 1.24, β = 1.01 and z 0 = 0.51. This specific example corresponds to the so-called "gold sample" of galaxies with shapes expected to be observed with the Rubin Observatory's Legacy Survey of Space and Time [64,65]. Note that the results shown here for shape correlations are for intrinsic alignments alone, i.e. without the contribution of the lensing signal. For simplicity, we assume that the same sample is used for both shapes and clustering, i.e. W (χ) ≡ W n (χ) = W g (χ). This is generally not the case in galaxy surveys, which tend to focus on either spectroscopic or color-selected clustering samples to achieve better redshift accuracy. Nevertheless, we adopt this simple case here for illustration purposes. In figure 2 we show how these weight functions contribute to the angular power spectrum and bispectrum. We see that the weight function W 2 g for the power spectrum peaks around the redshift z ≈ 0.2, while in case of bispectrum the weight function W 3 g (the square of what is shown in the figure) diverges at low redshifts as ∝ z 3α−4 .
In figure 3 we show the various cross and auto angular power spectra for galaxy number density and E and B-modes of galaxy shape. Below each sub-plot, we show the ratio between the terms obtained from the EFT one-loop prediction and the linear power spectrum. The exception is the bottom right case, since C BB is expected to be null at linear order. In this case, we normalize it to the EE power spectrum. The first angular power spectrum (upper-left panel) is the galaxy number density-number density contribution C nn . We plot the EFT contributions to the angular power spectra as given in section 5.1 of Paper I [43]. We show three curves corresponding to: the linear bias with linear matter power spectrum supposed to be contrasted with the data or simulations, these parameters would need to be fitted for. The values of these bias parameters will depend on the specifics of the objects we look at and will thus encode the information on, e.g., galaxy type and mass, and other details of the relevant small-scale physics. However, we would not expect the parameters to differ by orders of magnitude from the unity values. Similarly to this we show the contributions to the number density-shape angular power spectrum C nE (upper-right panel). Notice that the b 1 term in case of the shape spectra corresponds to the linear alignment model from [16]. From the ratio sub-plot shown below the C nn prediction, and in general for all other correlations shown in figure 3, we see that the EFT contributions are scale-dependent and have the highest impact at the smallest separations between galaxies (large values), as expected. In practice, the EFT expansion should only be trusted as long as the new terms are a fraction of the linear one.
The bottom panels of figure 3 show the intrinsic shape auto-correlations: C EE (left) and C BB (right), using the same bias value prescription. Our predictions indicate that C BB , which is null in linear alignment model predictions, is present at the mildly-nonlinear level. This is also in line with SPT predictions [40,66]. However, because E and B-modes are connected by the same choice of biases, the bottom right panel of figure 3 can be interpreted as the typical order of magnitude of B-modes at mildly non-linear scales in comparison to the E-modes. We find these to be suppressed on large scales.
So far we did not add in the discussion the contributions of shot noise. In [43] we have shown that the leading large-scale contribution to cross shape-noise and auto shape-shape 3D power spectrum gives where P is a constant power spectrum. In appendix B we show that this implies that the shot noise contributions to the angular power spectrum are while C nE ( ) receives no constant large-scale contribution. We have also used the abbreviation W = dχ[W (χ)] 2 /χ 2 . This implies that the shot noise of the C EE and the C BB are the same. Given that the contribution from the long mode coupling shown in figure 3 is highly suppressed, as discussed in the paragraph above, we can expect the C BB contribution to be noise dominated on the large scales. This also implies that the difference of the C EE and C BB (which is equal to C +− , as shown in eq. (4.5)) are free of the constant shot noise contributions and receives only the 2 contributions, similar to C nE ( ).
In figure 4, we show a plausible range for fractional contribution of one-loop terms to each linear angular power spectrum, as a function of multipole and based on the variation of bias parameters. The variation of the bias parameters corresponds to the choice given in figure 3 of Paper I [43]. The three panels represent C nn , C nE and C EE . Note that in the figures we do not show the stochastic bias contributions. The contributions are contained until scales of l ∼ 10 2 and become comparable to the linear prediction for larger multipoles Since this paper is focused on intrinsic alignment contributions to galaxy shape correlations, we have not included the contribution from gravitational lensing here. Assuming that lensing and alignment effects are additive, which holds to leading order, it is straightforward to include gravitational lensing as described in appendix A. The result for C nE shown here Figure 4. Relative one-loop bias contributions compared to the dark matter one-loop power spectra for C nn , C nE , C EE . We again use b 1 = 1, b R * = −3/2 and for the rest of the biases we use the values |b n | = b 1 /4. The sign of bias values b n are adjusted so that the grey area represents the band of potential bias one-loop contributions where |b n | < b 1 /4. In all cases we neglect the stochastic bias contributions. Effectively this gives an estimate of the validity regime of linear alignment model, provided our bias coefficient choice is realistic. then functions as a contaminant to galaxy-galaxy lensing, while that for C EE is a contaminant to cosmic shear. Note however that there are also important cross-contributions of intrinsic alignments and lensing, which are included in the relations given in appendix A.
Next, in figure 5, we show the number density auto-and number density, number density and shape cross-bispectra. As mentioned, for all bispectra cases we adopt the simplification for the window function W g = W n , thus all the weights are proportional to W 3 g . We adopt the same bias values as in the case of power spectra above. The tree-level 3D bispectra needed for the computation of the angular bispectra here are explicitly given in section 5.2 of Paper I [43]. All panels are normalized by the maximal value attained over the configurations shown. In the upper panels we show the contributions for the number density auto-bispectrum B nnn for configurations with maximal wavenumber 1 = 30 (left panel) and 1 Figure 5. The angular number density and shape auto-and cross-bispectrum. Left panels show bispectra with maximal wavenumber 1 = 30 while on the right panels it is 1 = 150. Upper panels show the number density auto-bispectrum B nnn , while the middle and lower panels show the cross-bispectrum of two number density and one shape fields in the B nnE and B nnB cross-bispectra configurations, respectively. All bispecra are normalized to its maximal amplitude configuration value. The maximal amplitudes are, typically, an order of magnitude smaller in B nnB compared to B nnE , for both 1 = 30 and 1 = 150 cases. Values of the bias parameters are the same as in the case of power spectra.
both cases we see that the equilateral contributions are suppressed compared to the 2 + 3 ≈ 1 diagonal (flattened triangles), where the amplitude the bispectrum amplitudes reach their maximal values. Middle panel shows the number density, number density and E-mode crossbispectrum B nnE for the same maximal wavenumber 1 choices. We notice the similar relative configuration distribution as in the case of B nnn , with the equilateral contributions suppressed relative to the 2 + 3 ≈ 1 diagonal. The property that the bispectrum amplitudes peak for the flattened triangles, where 2 + 3 ≈ 1 , is inherited from the similar behavior of the 3D bispectrum of matter and biased tracers (see the right panels of figure 10 in [45]).
Similar to the E-mode, in the lower panels we show the number density, number density and B-mode cross-bispectrum B nnB . For the same choice of bias parameters as before, the maximal amplitude of this bispectrum is suppressed by an order of magnitude relative to the B nnn and B nnE . Additionally, relative to this maximal amplitude configuration, the rest of -21 -JCAP05(2021)061 the configurations are further suppressed more strongly than in the B nnn and B nnE cases. In other words, in the B nnB cross-bispectrum, most of the relevant contributions are roughly centered around the maximal amplitude configuration and they are again clustered close to the flattened triangle shapes, i.e. 2 + 3 ≈ 1 is diagonal. The latter feature can again be explained by the shape of the 3D bispectra.
Let us comment on the fact that the maximal B nnB contribution is only one order of magnitude suppressed relative to B nnE . This is in stark difference compared to the power spectrum case, where the only channel of comparison was C BB compared to C EE . There, in the absence of shot noise, the large suppression of C BB on large scales was caused by the fact that the leading order contribution was a one-loop result, versus the linear one for the C EE . In the case of B nnE and B nnB , this is no longer so, and both of these bispectra have non-vanishing tree-level contributions, so that the B-mode contribution is less suppressed. Moreover, we expect the similar behaviour to persist in the other cross-and auto-bispectra contributions, as shown in Eqs, (4.18) and (4.19). The difference comes in the B nBB , B EBB , and B BBB terms that do not have tree-level contributions and are thus expected to be suppressed on large scales (similar to the C BB case in the power spectrum case). We note that we did not take the shot noise contributions into account here. These could again significantly affect the discussion above, and we refer the reader to Paper I [43] for a more detailed discussion of the noise contributions to the 3D shape auto-and cross-bispectra.
Lastly, we remind the reader that similar statistics have been derived in the context of other intrinsic alignment models. Examples are the LA model [67], the nonlinear alignment model [17] or SPT [39][40][41]. The formalism presented in this work is applicable to any of these models, the EFT, and any spin-2 observable of the large-scale structure, as it provides a general framework for projected statistics of such fields. Nevertheless, we can make the following remarks about how these models can be compared. First, it is clear that the linear alignment model is simply the lowest order prediction from either the EFT and SPT. This has been seen to underpredict the level of alignment between galaxies at small scales [11,68], and hence the need for higher-order models. The EFT and SPT are alternative methods to reach quasi-linear scales. The EFT presents an advantage with respect to SPT in that all higherorder terms are consistently included and degeneracies between them can be straightforwardly resolved. The two approaches should be equivalent as long as all the terms contributing to a galaxy shape within SPT can be identified. The nonlinear alignment model is on the other hand an empirical approach to modelling intrinsic alignment power spectra which makes an adhoc replacement of the linear matter power spectrum by the nonlinear one in otherwise linear equations. Although this enhances the intrinsic alignment power at small scales in a manner that is more consistent with observations, the empirical nature of the model limits its applicability to deriving predictions for alignment observables. For example, the nonlinear alignment model predicts no B-modes whatsoever. More details regarding the comparison between different models can be found in Paper I [43].

Conclusions
In Paper I [43], we had presented an EFT expansion of galaxy intrinsic shapes, sizes and number density tracers that allows one to perform a complete prediction of three-dimensional Fourier-space statistics, i.e. power spectra up to the one-loop order and tree-level bispectra. In this work, we developed the projection formalism that relates these 3D correlators, which describe the physical dynamics of the observables, with the projected quantities that live on -22 -

JCAP05(2021)061
the 2D plane of the sky. We adopted the flat-sky approximation for computing the angular power and bispectra. This provides an accurate description of angular power spectrum on scales where our one-loop corrections apply. Furthermore, our results, which establish the connection between components of the 3D statistics and projected 2D statistics, are applicable to more general 3D correlators, and not just the ones predicted by the EFT of Paper I. Indeed, these projections are performed independently from the methods by which the 3D statistics is obtained and instead rely entirely on the underlying statistical symmetries (homogeneity, isotropy and parity invariance).
We presented results for the angular two-and three-point statistics of galaxy counts and shapes. The results are presented using three different bases; the helicity basis, which is the natural basis for applying the above-mentioned symmetry constraints; the electricmagnetic component basis, which is the commonly used basis in the field and in which we present our final results; and lastly, the Pauli basis that is also frequently used in the lensing community. The development of this formalism allows one to identify the forms of degeneracies taking place when projecting into the angular power spectra and bispectra, as well as which contributions survive at the tree-level bispectrum and one-loop power spectrum. In case of the power spectrum, the situation is relatively simple, and all number density and E-mode auto-and cross-correlations have contributions already at linear order. In contrast, the cross-correlations of these with the B-mode vanish to all orders as a consequence of the statistical parity invariance (in addition to the homogeneity and isotropy). The remaining Bmode auto correlation gets its leading order contribution from the one-loop mode-coupling terms. It is for this reason few orders of magnitude suppressed compared to the number density and E-mode correlators (on the scales of our interest) which, as mentioned, benefit in that respect from the leading linear order contribution. In the case of bispectra, the statistical symmetry considerations are less constraining allowing for all possible number density, E-mode and B-mode auto-and cross-correlations. However, many of these do not contribute in the tree-level PT, appearing only at higher-loop orders. This is typically the case in bispectra that correlate more than one B-mode. An interesting observation in case of B-mode contributions to the bispectrum is that it no longer is as strongly suppressed (as was the case in the power spectrum) compared to the other correlations. The reason for this is that all of the considered cross-correlations (B nnn , B nnE and B nnB ) have contributions starting at the tree-level in PT.
Assuming plausible bias parameter values, we made estimates of these correlators using the EFT results of Paper I [43], and we present several numerical results for the statistics discussed above. The current study is only qualitative and for more robust quantitative performance assessment, comparisons to simulations would be needed. On the other hand, our formalism can be readily applied to ongoing weak gravitational lensing surveys of the large-scale structure to model position-shape and shape-shape correlations, as well as threepoint statistics, in a consistent and complete manner up to quasi-linear scales. Moreover, we emphasize that its potential applications are even more general: it can be applied to any study of tensorial fields on the sky. In the future, we plan to present analogous allsky derivations that can improve the accuracy of our predictions at low values, as well as predictions for configuration space correlation functions in this regime.

A Gravitational lensing
We now provide expressions that include the gravitational lensing contribution to shapes, under the approximations stated after eq. (3.2): first, we assume that lensing simply adds to the intrinsic shape; second, we neglect all post-Born, reduced-shear and lensing-bias corrections, so that γ G,ij can be written as the second derivative of a lensing potential, which in turn is proportional to the gravitational potential. This means that γ G only receives helicity-0 contributions.
In the polarization basis, defined in eq. (3.16), the angular power spectra including lensing can then be written as 22 ( /χ)

( /χ).
Here, is the lensing kernel in the absence of curvature. P nm (k, z) is the cross-power spectrum between matter density and number counts n (which can be calculated up to 1-loop order using the EFT of biased tracers), while P mm (k, z) is the nonlinear matter power spectrum. P (0) 2m (k, z) is the cross-correlation between matter and the helicity-2 component of the intrinsic shape. This can be obtained in the EFT of shapes from P (0) 02 (k, z) by setting all the numbercount bias parameters to zero, except for b 1 which is set to 1. For explicit expressions for the EFT power spectra, we refer the reader to [43]. We defer the inclusion of the lensing contribution to bispectra statistics [e.g. 53, 60, 62] to future work.

B Angular power spectrum
Let us assume we have a real tensor field X ij (r) such that we can identify the trace part as the simple scalar overdensity δ(r), while the trace-free part can be decomposed into the helicity two components γ ±2 (r). In other words, we can write ij (r) = δ K ij . It is convenient to use this field in order to compute the relevant statistics. Moreover, we assume that the power spectrum of this field is of the form given in eq. (2.10), i.e. we have We are interested in angular statistics on the sky. Here, we focus on the flat-sky approximation results and thus we can introduce the field (see figure 1 for coordinate setup) in analogy to eq. (3.15), and wherê The flat-sky projected two-point function can be expressed as It is convenient to setup the coordinate frame so that k =n.k and k ⊥ = k − k n, and using the variables χ 2,1 = χ ± 1 2 ∆χ. This gives us To further simplify the calculation note that the leading wavenumber k modes that can contribute to the integral are those ones with small k =n.k, so that k 1/(χθ). Since we work in the small-angle approximation, we have 1/(χθ) 1/χ and thus modes with longitudinal wavenumber k much larger than 1/χ do not give rise to angular correlations because of cancellations along the line of sight. Only modes with k similar to 1/χ lead to angular correlations. The relevant transverse wavenumbers 1/(χθ) are thus much larger than the relevant longitudinal wavenumbers, and we can set the argument of the spherical components of the tensor power spectrum to from which it follows that

JCAP05(2021)061
We can now integrate over k which gives us the delta function that sets ∆χ = 0 and This allows us to compute the flat-sky projected angular power spectrum. Using the definition in eq. (B.3) and the expression above we have We also used the fact that Given the decomposition of P ijkl given in eq. (2.10) we have Using the definition of flat-sky angular power spectrum we arrive at the final expression: ij , introduced in eq. (2.7). Since all the components of the C s 1 s 2 angular power spectra are real, we have X s 1 ( 1 )X * s 2 ( 2 ) = X * s 1 ( 1 )X s 2 ( 2 ) . (B.16) As an example we can have a look at the leading noise contribution to the size-shape and shape-shape power spectrum given by as shown in [43], and where P is a constant in k. As shown there, this noise spectrum can be decomposed in the form given in eq. (2.10) with components P (0) 00 (k) = P Using the definition of the 3D bispectrum given in eq. (2.12), and the same set of approximations given in appendix B (e.g. k k ⊥ ), we have dk ,123 (2π) 3 e −i l χ l k ,l X ij k ,1 , 1 /χ 1 , z 1 X nm k ,2 , 2 /χ 2 , z 2 X rs k ,3 , 3 /χ 3 , z 3 (C.2) Using the definition of flat-sky angular bispectrum Using these results and eq. (2.10), one can also recover the angular power spectrum results given in eq. (4.3) and derived explicitly in appendix B. However, in this section we focus on the angular bispectrum results. Using the 3D bispectrum decompositions as given in section eq. (2.14) we have Where ({a, b, c}) notation is used and we understand that the summation over all the nontrivial index permutations should be performed. Above we also used the bispectrum parity relation given in eq. (2.15) to reduce the number of independent terms. Each of the cross angular bispectra has contributions from different helicities, which give rise to the terms in the integrand. Depending on the number of helicity two fields in the angular bispectrum number of terms in the integrand varies from three for the density-density-shear (002) bispectrum, six for the density-shear-shear (022) bispectrum and ten for the shear-shear-shear (222). These multiplicities correspond to the number of combinations with repetition for the helicity of the shear field.

JCAP05(2021)061
In Paper I [43] we gave the explicit bispectrum results for biased tracers up to the treelevel (leading order) PT. In that case only, few of the terms above survive given that leading order fields carry only helicity zero contribution. In other words, only single non-zero helicity terms survive and we have These results for B LO nns 3 agree with the projections used in Paper I [43] to decompose the bispectrum at tree-level in PT.