{'fraction_non_alphanumeric': 0.05131194588223236, 'fraction_numerical': 0.018564546712448112, 'mean_word_length': 4.082568209936837, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Galaxy imaging surveys as spin-sensitive detector for cosmological colliders

Galaxy imaging surveys provide us with information on both the galaxy distribution and their shapes. In this paper, we systematically investigate the sensitivity of galaxy shapes to new physics in the initial conditions. For this purpose, we decompose the galaxy shape function into spin components, and compute the contributions to each spin component from both intrinsic alignment and weak lensing. We then consider the angular-dependent primordial non-Gaussianity, which is generated by a non-zero integer spin particle when active during inflation, and show that a galaxy imaging survey essentially functions as a spin-sensitive detector of such particles in the early universe. We also perform a forecast of the PNG generated from a higher spin particle, considering a Rubin Observatory LSST-like galaxy survey.


Introduction
Over the past few decades, our understanding of the Universe has improved dramatically by virtue of various cosmological observations. In particular, the cosmic microwave background (CMB) observations have played an important role in establishing the concordance cosmology, i.e. ΛCDM model. The observed existence of dark matter and dark energy urges us to seek for a new idea beyond standard model of particle physics. In the next decade, their properties will be investigated more profoundly through many galaxy surveys such as the Rubin Observatory Legacy Survey of Space and Time (LSST) [1], Euclid [2], the Nancy Grace Roman Space Telescope, which is formerly known as WFIRST [3], and Spectro-Photometer for the History of the Universe, Epoch of Reionization, and Ices Explorer (SPHEREx) [4]. While the CMB only observes the two-dimensional surface, galaxy surveys chart the three-dimensional distribution of galaxies. Since enlarging the survey volume increases the constraining power, the upcoming galaxy surveys with a huge volume may enable a more precise measurement than the CMB. More importantly, by probing smaller-scale perturbations than the CMB, galaxy surveys allow us to explore the structure of the Universe on a broader range of scales. Through the galaxy imaging surveys, we can observe not only the galaxy distribution but also galaxy shapes. Similarly to the CMB, for which measuring the polarization provides complementary information to measuring the temperature fluctuation [5][6][7][8], measuring galaxy shapes supplies us with unique information about the structure of the Universe.
In particular, the cosmic shear caused by the weak lensing effect (see [9] for a review) is a powerful tracer of the dark matter distribution, which was reported e.g. in Dark Energy Survey [10], Kilo-Degree Survey [11], and Subaru HSC observation [12]. Broadly speaking, the primary target of galaxy imaging surveys has been measuring the weak lensing effect or equivalently the dark matter distribution. Nevertheless, the observed galaxy shear can be generated also by the external tidal fields, which is known as 'Intrinsic Alignment (IA)' [13][14][15][16] (for reviews, see Refs. [17][18][19][20]). The measurements of the cosmic shear have been used to constrain the cosmological parameters such as the equation of state for dark energy [21][22][23][24], the fluctuation amplitude σ 8 and the matter density Ω m0 [24][25][26][27][28], and parameters for a deviation from general relativity [29]. Although the IA has been usually treated as a contamination to the weak lensing contribution, the IA itself includes precious information on the model of the Universe [30].
One may wonder if measuring the IA can also be a powerful tool to seek a fossil of the primordial Universe. Inflation is the most successful scenario of the early universe, which consistently explains the seed of the fluctuations measured by the CMB and the large scale structure observations. A detection of the primordial gravitational waves (PGWs) generated during inflation can pin down the energy scale of the inflation, giving the benchmark value for building models of inflation. The measurement of the CMB B-mode [6,7,31] is a promising way to detect the PGWs. Meanwhile, a consequence of the PGWs is also encoded in the cosmic shear through the gravitational lensing [32] and also through the IA [33,34] (see also Ref. [35]).
Without a detection of the PGWs, the energy scale of the inflation is just bounded by above. According to the latest Planck data [36], the energy scale of inflation can be as high as H inf < 2.7 × 10 −5 M pl ∼ 6.6 × 10 13 GeV in 95%CL, which by far exceeds the accessible energy scale by any ground-based accelerator experiments. In such a high energy environment, an extremely massive particle, which cannot be explored by the conventional accelerator experiments, can be excited. If the particles had interacted with the inflaton, that drove the inflationary expansion, their imprints can be encoded in the primordial perturbations. The cosmological collider program [37][38][39] aims at probing massive particles which had existed in inflationary Universe through a precise measurement of the primordial perturbations. The information of the mass and spin for such heavy particles are encoded in the squeezed primordial non-Gaussianity (PNG) as the distinct oscillatory feature [37] and angular dependence [39], respectively. In particular, it was shown that a particle with a different spin leads to the PNG with a different angular dependence [39][40][41] (see Ref. [42] for a different approach).
Detecting an imprint of these massive particles can open a unique window to explore new physics beyond standard model such as string theory. One may explore an imprint of grand unification encoded in the angular dependent PNG generated by the predicted spin-2 Kaluza-Klein gravitons [43]. Going one step further, a typical prediction of string theory includes an infinite tower of higher spin particles, whose spins s are larger than 2 (see e.g., the review article [44] and references therein). The seminal works by Vasiliev [45,46] show that the no-go theorem for higher spin gauge theories known in the flat spacetime [47][48][49] can be evaded by introducing a negative or positive cosmological constant, including the de Sitter spacetime, which serves a good approximation of an inflationary spacetime. Probing a higher spin particle which might have existed during inflation is the most challenging goal of the cosmological collider program, which nevertheless can lead to a direct evidence of string theory.
The angular dependent PNG has been searched through CMB [50][51][52][53] and LSS [51,[54][55][56][57][58][59][60][61]. Since the squeezed PNG generated by a scalar field, the inflaton 1 or other spectator scalar fields, is typically bigger than the one by a non-zero spin particle, we need to establish an observational method which can selectively pick up the latter without being hidden by the former. In this regard, observing galaxy shapes from imaging surveys has a compelling advantage [54]. It is widely known that the angular independent squeezed PNG, whose amplitude is often expressed as f NL , causes the scale dependent bias of the galaxy number density [64][65][66][67][68][69]. Similarly, the angular dependent PNG which can be generated by a spin-2 particle induces a scale-dependent bias for the IA correlation of the 2nd shape moment, which is the usual cosmic shear [54][55][56]. Recently, this was confirmed in full N -body simulation by Ref. [70]. Remarkably, no matter how large f NL is, it does not contaminate the second shape moment at large scales where the non-linear evolution is negligible. Furthermore, even if there exists an infinite tower of higher spin particles that generate different types of angular dependent (squeezed) PNG, they do not contribute to the 2nd shape moment either. This leads us to an ambitious speculation that an ideal galaxy imaging survey that can separately observe different galaxy shape moments may enable us to measure the PNG generated by massive particles with different spins separately, playing the role of a spinresolved spectroscopy. In this paper, we confirm this speculation, generalizing the analysis in Ref. [54] to an arbitrary nth galaxy shape moment.
The goal of this paper has two folds; First, we compute various effects that contribute to the nth galaxy shape moment, assuming the ΛCDM cosmology with the adiabatic Gaussian initial condition, based on the formalism of Refs. [33,71] (see [72][73][74][75] for related approaches). This generalizes known results for the 2nd moment to all higher moments, similarly to the recent Ref. [76]. In order to find new physics by observing galaxy shapes, the imprint of new physics needs to dominate these contributions at least over a certain range of scales. Second, as an example of new physics, we consider the angular-dependent PNG generated by higher spin particles, showing that galaxy imaging surveys indeed work as spin-sensitive detectors of the cosmological collider, which can separately measure the PNG generated by particles with different spins.
The outline of this paper is as follows. In Sec. 2, we introduce the definition of the galaxy shape moment and its spin decomposition. In Sec. 3, we estimate various contributions to the galaxy shape moments in the ΛCDM model with the adiabatic Gaussian initial condition, providing a benchmark value for a signal of new physics to exceed for its discovery. In Sec. 4, we compute the scale-dependent bias generated by the angular dependent PNG from non-zero integer spin particles. In Sec. 5, we forecast the detectability of the PNG from spin-4 particles, considering the future galaxy imaging surveys. In this paper, we assume the Planck fiducial [77] flat ΛCDM cosmology with Ω b0 h 2 = 0.022, Ω CDM0 h 2 = 0.12, h = 0.67, n s = 0.9645, A s = 2.2 × 10 −9 and k p = 0.05 Mpc −1 . With Ω b0 and Ω CDM0 , Ω m0 is defined as Ω m0 ≡ Ω b0 + Ω CDM0 . We have used the public Boltzmann code CAMB [78] to calculate the transfer function and the matter power spectrum at present.

Spin decomposition of galaxy shapes
In this section, we introduce the galaxy shape function which we will then use as a window to explore new physics.

Definition and spin decomposition
First, let us define the shape function for a galaxy located at (τ, x) with τ being the conformal time as where x and y denote the 3D spatial coordinates for the centroid of the galaxy and those measured from the 3D centroid of light emission from the galaxy, respectively. The indices i 1 , · · · , i n run from 1 to 3 and are raised and lowered by the 3-dimensional Kronecker delta. We have introduced the weight by using the brightness of the emission, B(x, τ ). Here, R * is the size of the galaxy andB(x, τ ) is the total brightness of the galaxy, defined as B(x, τ ) = d 3 yB(x + y, τ ). As is clear from the definition, g i 1 i 2 ···in is the rank-n symmetric tensor in 3D. The centroid is defined so that g i (x, τ ) vanishes. For a localized galaxy, the integral in g i 1 i 2 ···in converges. There are two different origins for the deformation of the galaxy shape, The first term is the intrinsic galaxy shape deformation and the second term is the deformation due to the weak lensing, which is caused by the gravitational field between galaxies and us. We will further discuss each component in the next section.
In galaxy imaging surveys, we can only observe the galaxy shapes projected onto the 2D sky which is orthogonal to the line of sightn. Operating with the projection tensor, defined as withn being the unit vector along the line of sight, we obtain the galaxy shape on the projected plane as Over the area of an individual galaxy, which subtends a very small angle on the sky, we can approximate the line of sightn as constant, so that the sky can be approximated as a two-dimensional plane with coordinates θ, where we choose the origin to correspond to the centroid of the galaxy image. The surface brightness can be obtained by integrating the emission along the line of sight, where we ignore any absorption or other radiative transfer effects. However, our results below are based solely on symmetry considerations, so they continue to hold even in the presence of a more complicated local mapping from emissivity to surface brightness. In Eq. (2.5), using the source comoving distance χ, defined as we have introduced y ,θ i , and θ i as y ≡ y ·n, respectively. Here, the dependence on x·n was identified with τ dependence. We then obtain the projected galaxy moment as (2.8) which is normalized by the total intensityĪ, Since g i vanishes, so do I i , satisfying I i (θ, τ ) = 0. Next, let us introduce the spin±n components of I i 1 ···in . Here, we generalize the approach of [71]. An alternative, closely related approach to higher-order moments has recently been presented in Ref. [76]. It will turn out that the contribution of the PNG generated by a spin-n particle selectively appears in the spin±n components at large scales. Using the orthonormal basis, (n, e ι , e ψ ), let us introduce the spin±1 unit basis vectors as which satisfy The projected nth moment tensor I i 1 i 2 ···in can be expanded by a product of m ±i as m ±i 1 m ±i 2 · · · m ±in , (m ±i 1 m ∓i 2 · · · m ∓in + (perms)), · · · . (2.12) Apart from the first terms, m ±i 1 m ±i 2 · · · m ±in , all other terms, which contain both m + and m − , have trace part contributions, as one can confirm by using Eq. (2.11). Therefore, the traceless components of I i 1 i 2 ···in ,Ĩ i 1 i 2 ···in , which satisfy, for any p, q ∈ [1, n], can be expanded by the first terms of Eq. (2.12) as where the coefficients ±n γ are given by In the second equality, we inserted the trace part of I i 1 i 2 ···in , which does not contribute to Since m ± transform as spin±1 vectors under a rotation by an angle β, i.e. m ± → e ±iβ m ± , the n products of m ± transform as spin±n basis tensors, m ± · · · m ± → e ±inβ m ± · · · m ± , being invariant under the rotation by an angle β = 2πN/n where N is an integer.

Visual image of spin-n distortion
In the previous subsection, we argued that the traceless components,Ĩ i 1 ···in , can be expressed by the spin±n components, ±n γ. Before we calculate ±n γ, let us further discuss the spin decomposition to develop an intuitive understanding. Relevant discussions can be found e.g., in Refs. [76,79]. For this purpose, here we introduce the 2D coordinates on the projected plane along the directions of e ι and e ψ . In this paper, we express the 2D coordinates on the projected plane as a 1 , a 2 , · · · ∈ {ι, ψ}, distinguishing them from the 3D spatial coordinates, expressed as i 1 , i 2 , · · · = 1, 2, 3. Introducing where e a i is the ith component of e a , we can rewrite the spin±n components, ±n γ as where p and q run over p, q = 1, · · · , n, ±n γ can be given by a more compact expression as Notice that an index for the former runs a 1 , a 2 , · · · = ι, ψ and the one for the latter runs i 1 , i 2 , · · · = 1, 2, 3. As discussed above, the spin ±n components, ±n γ, can be expressed by the two traceless components ofĨ a 1 a 2 ···an , which are the only independent components ofĨ a 1 a 2 ···an . In fact, with a use of the symmetric and traceless condition (2.19), all other components ofĨ a 1 a 2 ···an can be expressed by these two components in the right hand side of Eq. (2.20) as for k = 0, 1, · · · , [ n− 1 2 ]. It will turn out that the imprint of the PNG generated by the spin-n particles is encoded in the spin-n component of the galaxy shape function without being mixed with the contributions of the PNG from particles with other spins.
Next, let us discuss a visual image of the spin-n distortion, described byĨ a 1 a 2 ···an . For this purpose, using the 2D polar coordinates (r, φ) with which θ ι ≡ e ι i θ i and θ ψ ≡ e ψ i θ i are given by we express the surface brightness I(θ a ) in the Fourier series expansion as Apart from θ a = (θ ι , θ ψ ), the surface brightness also depends on τ and the centroid, while we don't write them explicitly. Using Eqs. (2.16) and (2.23), we can express the projected 2D galaxy moment I a 1 ···an as an nth moment described by the integral of θ a 1 · · · θ an , which has been used, e.g., in Ref. [9]. The condition I a = 0 leads to c 1 = s 1 = 0. Further, we have s 0 (r) = 0. Notice that the presence of a non-zero c n or s n with n ≥ 2 indicates a deviation from the circular symmetric distribution I(θ a ) = I(|θ a |). We now show that the two independent spin-n components ofĨ a 1 a 2 ···an with n ≥ 2, I ιι···ιι andĨ ιι···ιψ , correspond to c n (r) and s n (r), respectively. As a simple exercise, let us start with n = 2, which corresponds to the usual cosmic shear. Inserting Eq. (2.24) intõ I a 1 a 2 = I a 1 a 2 − δ a 1 a 2 I bb /2, we obtaiñ where we have introduced C n ≡ π drr n+1 c n (r) , S n ≡ π drr n+1 s n (r) . (2.27) Extending this discussion to n ≥ 3, we obtaiñ where we can avoid a messy computation by virtue of the formulae (the derivation can be found in App. A.2) r n cos(nφ) = 2 n−1 [θ a 1 θ a 2 · · · θ an ] TL 2 a 1 =···=an=ι , (2.29) r n sin(nφ) = 2 n−1 [θ a 1 θ a 2 · · · θ an ] TL 2 a 1 =···=a n−1 =ι,an=ψ . (2.30) Thus, we find that the spin-n components of the galaxy shape function are directly related to a deviation from the circular symmetric distribution which is proportional to cos(nφ) or sin(nφ). As it should be, these terms remain invariant under the rotation by an angle ∆φ = 2πN/n where N is an integer. Figure. 1 shows the deviation of I(θ) from the circular symmetric distribution, described by c n (r) and s n (r) for n ≥ 2.

Angular power spectrum
After having described how an individual galaxy image is decomposed into spin-n components, we now generalize to spin-n fields on the skyĨ i 1 ···in (n), wheren is a unit vector denoting a location on the sky. This field can be estimated by dividing the survey area of a galaxy imaging survey into pixels, and measuring the average shape component of the many galaxies in each pixel.
Assuming only the adiabatic perturbation to the FLRW metric, and working to linear order in perturbations, we can relateĨ i 1 ···in to the matter density perturbation δ m (k) without loss of generality viã where z is the redshift, dN/dz is the redshift distribution function of galaxies, δ m (k) is the total matter density fluctuation at present, D(z) is the growth factor normalized as (1 + z) −1 during matter domination, µ =k ·n, and x = kχ(z). All physical effects are encoded in the kernel F (n) (k, µ, z) which will be derived in the next two sections. As discussed above, there are two different contributions to F (n) (k, µ, z), which stem from the intrinsic alignment and the weak lensing, respectively (at linear order in perturbations, there is no coupling between these two contributions). As a consequence of the global rotation symmetry, for scalar perturbations the tensor indices ofĨ i 1 ···in can be expressed only in terms of the (projected) wavenumbers k j 1 · · · k jn . In this paper, we further assume that the matter power spectrum preserves the global rotation symmetry as on +n γ, we obtain a scalar quantity γ (n) as where we have introducedk ± ≡ m i ∓k i . Using this expression, we compute the angular power spectrum defined as Following Ref. [88], let us define the E-mode and B-mode of the spin-n components as The matter density perturbation δ m (k) will be replaced with the matter power spectrum at present P m (k), which preserves the global rotation symmetry, after taking the correlation in the angular power spectrum (2.35). Then, since the contribution from each k to a (n) lm does not depend on the angular direction of k, we can simply compute the angular integral in a (n) lm , choosing k i = (0, 0, k) as Here, we used that a spin-s function s f with s > 0 which only depends on the polar angle θ satisfies [71,[89][90][91]ð where µ =k ·n = cos θ. We can further rewrite the angular integral, using where the argument ∂ x in F (n) means that the derivative operates on all x-dependent terms to the right in F (n) . Finally, using the partial wave expansion and the following formula where we use ∼ instead of the equality, since we have already employed the global rotation symmetry of the matter power spectrum, which becomes manifest only after taking the correlation. Using this expression, we find that whenĨ i 1 ···in is given by Eq. (2.31) and the matter power spectrum preserves the global rotation symmetry as given in Eq. (2.32), the B-mode distortion identically vanishes. The B-mode distortion can be generated from the primordial gravitational waves (PGWs) [32,34] and also from a violation of the global rotation symmetry in the PNG [56]. In Ref. [56], it was shown that the off-diagonal EB correlation (parity-even) can be generated from the violation of the global rotation symmetry and in Ref. [61], it was shown that the diagonal EB correlation (parity-odd) can be generated from the helical PGWs. (Recall the argument below Eq. (2.38).) Using Eq. (2.42), we obtain the angular power spectrum for the E-mode of the spin-n galaxy shape component as Here and hereafter, for our notational brevity, we drop the index E, since we only consider the E-mode. Unless stated, we use the linear spectrum P L (k) for the matter power spectrum P m (k). For a later use, we introduce which captures the rough structure of the integrand of C (n,n ) l for a given k and z.

Galaxy shapes as observable
In order to probe a consequence of new physics through shape observation, we first derive the power spectrum of the galaxy shape functionĨ i 1 ···in predicted in the ΛCDM cosmology with adiabatic Gaussian initial condition. The galaxy distribution is known to represent a biased tracer of the matter density distribution. Similarly, the galaxy shape represents a biased tracer of the tidal field. All the contributions toĨ i 1 ···in can be summarized as The observed galaxy shape can be decomposed into the intrinsic shape (as would be measured by an observer close to the galaxy) and the weak lensing shear due to the propagation of light through the inhomogeneous universe (projection effect). We further divide the former into the linear alignment (A), which behaves as a linearly biased tracer of the tidal field, and the non-linear alignment (B), which corresponds to the higher order correction to the linear alignment. In addition, since galaxy shapes are intrinsically non-spherical, there is a random (noise) contribution to the shape. For the 2nd moment, this is known as shape noise. For higher order shape contributions, the magnitude of the noise is uncertain. In the rest of this paper, we will ignore the noise in the galaxy shapes, considering a noise-free ideal observation. Weak lensing also adds contributions to the shape both at linear and non-linear orders.
In this section, we compute the scale dependence of each contribution one by one. In Sec. 3.1, we will discuss the bias expansion of an arbitrary moment galaxy shape function. In Sec. 3.2, using this expansion, we will compute the linear and non-linear alignment contribution. In Sec. 3.3, we briefly summarize the weak lensing contributions, deferring the detailed computation to App. B.

Bias expansion
The goal of the bias expansion of galaxy density and shapes is to capture the effect of longwavelength modes on these observables. By completely parametrizing this dependence with free parameters, it is possible to describe the large-scale statistics of galaxies and their shapes rigorously [92,93], even though the formation of galaxies is extremely complex and nonlinear. This approach works via an expansion in perturbations and derivatives, and thus applies only on scales where the matter density δ is much less than one, which can be phrased as k < k NL where k NL is the nonlinear wavenumber, and which are larger than the length scale R * on which galaxies form (k < k * = 1/R * ).
The position and shape of a galaxy at a given time τ depend on the trajectory of the non-linear matter density around the galaxy under consideration. The matter density is affected by the small scale physics such as the gas cooling, radiation and star formation, and the past trajectory of matter around the focused galaxy. However, when we consider a much larger scale than R * , which is the typical size of galaxies, the position and shape of galaxies are predominantly determined by the coarse-grained gravitational potential Ψ. In this paper, choosing the conformal Newtonian gauge, we express the line element for the perturbed FLRW spacetime as where a(τ ) is the scale factor and Ψ, Φ correspond to the gravitational potential and the curvature perturbation in this gauge. In the absence of the anisotropic pressure, the Einstein equation relates Ψ and Φ as Ψ + Φ = 0. Here, we set the background spatial curvature to 0 and ignore the tensor perturbation and the non-linear metric perturbations. Focusing on much larger scales than R * , we express the number density and shape of galaxies evaluated at the Eulerian spacetime coordinates (x, τ ) as where δ, K ij , and s i 1 ···in are defined as whereρ m denotes the background matter energy density and D i 1 i 2 is defined as On the second equality of Eq. (3.4), we used the Poisson equation. The symmetric and traceless tensor K i 1 i 2 is called the (second-rank) tidal tensor. The equivalence principle states that gravity can be locally eliminated, making the system equivalent to a local inertial system. This ensures that we can eliminate the gravitational potential Ψ and its derivative ∂ i Ψ by performing a local coordinate transformation to choose a free-falling observer frame. Therefore, in Eqs. (3.2) and (3.3), we did not include the gravitational potential and its first derivative. As a result, the energy density δ and the tidal tensor K i 1 i 2 give the leading order local contributions. Note that the galaxy density and shape in general depend on the matter density and tidal field on their entire past history. This can also be incorporated within the perturbative expansion [94]; in case of the galaxy density, contributions beyond the local dependence on δ and K ij appear at third order. The same holds for the n = 4 moment of galaxy shapes, while in case of the 2nd moment, there is a contribution at second order [93]. Since we are interested in the leading contribution to n = 4 moments, which is given by the second-order contribution, it is sufficient to work with the local, Eulerian density and tidal field in the bias expansion.
To determine the functional form of N and G i 1 ···in , we need to trace the history of the galaxy formation. Instead, here, along the line with Refs. [54,93,95,96], we expand δ n and the traceless part of g i 1 i 2 , using the dominant local contributions as X with X = δ, δ 2 , · · · and n = 0, 2 denote the bias parameters which are defined as the response of the mean nth moment galaxy shape to the change of X. For example, the linear bias parameter for the number density, b (0) δ , is given by the response of the mean galaxy density with respect to the background matter density as ∂ lnn g /∂ lnρ m and the one for the 2nd moment galaxy shape, b (2) K , is given by the response of the mean 2nd galaxy shape to the external tidal field K i 1 i 2 and so on. Here and hereafter we put a tilde to denote the traceless part of g i 1 ···in , i.e.g i 1 ···in ≡ [g i 1 ···in ] TL 3 . In Eqs. (3.8) and (3.9), the higher order terms in perturbation and the terms suppressed by |R * ∂| ∼ R * k with k being the Fourier mode are abbreviated.
Similarly, we can expandg i 1 ···in with n ≥ 3, using δ, K i 1 i 2 , and s i 1 ···in . For example, the 4th momentg ijkl up to the quadratic order in perturbation is given bỹ where the coefficient b K 2 describes the response to the quadratic external tidal field 2 . The terms with

Linear and non-linear intrinsic alignment
Using the bias expansion, discussed in the previous subsection, we can compute the linear and non-linear intrinsic alignment. The linear alignment term which contributes to F (n) , introduced in Eq. (2.31), is given by with k * ≡ 1/R * . We set k * as k * = 1 Mpc −1 for our numerical results. Notice that for n ≥ 3, the leading linear alignment term is a higher-derivative term and suppressed for k k * . Therefore, for n ≥ 3, the most significant contribution of the intrinsic alignment on large scales comes from the non-linear alignment.
Non-linearity in the galaxy alignment appears in two different ways. First, the relations (3.2) and (3.3) are already non-linear, as shown in Eqs. (3.8) and (3.9). Second, each argument of Eqs. (3.2) and (3.3), i.e. δ, K i 1 i 2 , · · · evolves non-linearly, which is the usual non-linear clustering effect. While we considerĨ i 1 ···in , which is given by projectingg i 1 ···in to the 2D sky, the loop integrals should be performed in the 3D space. The 3D loops are qualitatively different from the non-linear weak lensing effects, whose loop integrals are performed in the projected 2D space.
The non-linear alignment (NLA) for the 2nd moment was discussed, using simulation in Refs. [97][98][99], galaxy catalogs in Ref. [100] and bias expansion in Refs. [93,101]. Here, let us provide a crude estimation of (3D) loop integrals from non-linear alignment for galaxy shape moments. The scaling of loop contributions in k → 0 is known to be different between the unbiased and biased tracers. While the 1-loop of the former scales as k 2 P L (k) and k 4 , the one of the latter scales as P L (k) and k 0 . As discussed in Refs. [102,103], the scaling of the stochastic term for unbiased tracers, k 4 , can be derived from the mass and momentum conservation. Since the conservation is not imposed in our bias expansion, the stochastic terms with ∝ k 0 appear without being canceled for biased tracers. In Ref. [104], the 1-loop contribution to the power spectrum of the number density of halos was computed based on EFTofLSS. A representative term which scales as k 0 in k → 0 is given by where P hh is the halo auto-power spectrum andĉ 2 (k 1 , k 2 ) is the subkernel in the second order kernel for biased tracer and given bŷ (3.14) Correlation functions of shape momentsg i 1 ···in additionally have dimensionless tensor legs, which simply modifies O(1) factor. Loop computations for the 2nd moment can be found in Refs. [93,101]. Having considered this, we use to provide a crude estimation of the NLA contribution for the nth shape moments. This contribution is shown in Fig. 2 by brown dot-dashed lines. The bias parameter is set to |b (n) NLA | = 0.1 both for the 2nd and 4th shape moments. Since the 1-loop contribution which appears from the linear alignment term of g ijkl is suppressed by (k/k * ) 2 for the crosscorrelation and (k/k * ) 4 for the auto-correlation, the above contributions only appear from K 2 type term like the one in Eq. (3.10).
As is known, when the linear power spectrum is given by a simple power law spectrum with being k-independent, the scaling of the l-loop contribution can be estimated as 3 (see e.g., [92,95,106]) Here, k NL denotes the non-linear scale, at which the dimensionless matter power spectrum becomes unity. The 1-loop contribution shown in Fig. 2 follows this scaling around k NL , which explicitly corresponds to n L = 1.5 at k NL,z=0 = 0.25h Mpc −1 . For the actual Universe, since the power spectrum does not scale with a single power, the estimation becomes more complicated (see e.g., Ref. [107]). Since k * > k NL , these loop contributions dominate the linear alignment terms, suppressed by (k/k * ) 2 for the 4th moment as given in Eq. (3.12).
For the nth moment function with n = 2m (m = 1, 2, · · · ), the leading contributions, which are not suppressed by powers of k/k * , only appear from the non-linear bias expansion terms schematically in the formg i 1 ···i 2m ∼ (K) m , where the tensor indices of K ij are abbreviated. For a larger n, the non-linear alignment contributions without k/k * suppression start with higher loops. For instance, the auto-correlation ofg i 1 ···i 6 starts from 2loop contributions.

Weak lensing
In this subsection, we briefly summarize the weak lensing contribution to the galaxy shape function, whose detailed computation can be found in App. B. At leading order in the lensing deflection (Born approximation), the weak lensing contribution can be expressed by using the deformation matrix A ij [108] as Here,θ andθ s denote the 2D coordinates on the image plane and the source plane, respectively. We put a tilde to distinguish them from the offset from the centroid, introduced in the previous section. Because of the spatial inhomogeneity of Φ and Ψ, the deformation matrix A ij depends onθ.
The weak lensing contribution to the nth moment function,Ĩ WL i 1 ···in , should have n tensor indices which are projected to the 2D plane. For the 2nd moment, the lensing contribution, I WL ij is nothing but the deformation matrix A ij . Meanwhile, for n ≥ 3, the additional tensor indices are supplied either by acting with angular derivatives or by multiplying several A ij s. The former yields a suppression by powers of k/k * , while the latter contributes as twodimensional loop integrals.
The linear weak lensing (LWL) contribution, whose indices can thus only be supplied by operating with the spatial derivative operator, are suppressed as where we have dropped the redshift dependence. This result is in agreement with Eq. (84) of [76]. Therefore, similarly to the intrinsic alignment, for n > 2, the leading contribution for k k * stems from the loop contributions. Notice that since A ij , which is given by integrating along the line of sight, is a 2D object, loop integrals of the weak lensing, expressed by nonlinear terms of A ij , should be performed in the projected 2D space. Therefore, weak lensing loops are qualitatively different from the loop contributions in the intrinsic alignment, which are given by first computing loop integrals in 3D space and subsequently projecting it to the 2D space. As discussed in App. B, the lensing loops obtained after projection are much smaller than the latter for the auto-correlation of the 4th moment.

Summary
In order to probe an imprint of new physics encoded in galaxy shapes, we need to understand the contributions predicted in the concordance cosmology, which is ΛCDM cosmology, so in general relativity, with the Gaussian adiabatic initial condition. The sweet spot of hunting new physics is located at scales or angular multipoles where a signal of new physics can be exposed without being hidden by the contributions predicted in the standard cosmology.
In this section, estimating each contribution to the nth moment shape function, we have found that the situation is very different between n = 2, which has been widely investigated, and n ≥ 3. This is essentially because for n ≥ 3, both the intrinsic alignments and the weak lensing are suppressed by a factor of (k/k * ) 2 at linear order in perturbations. Then, the leading contributions appear from the loop contributions, especially those in intrinsic alignment. Meanwhile, as is widely known, the 2nd moment or the spin±2 component  includes the linear alignment and the linear weak lensing effect which are not suppressed at large scales. Figure 2 shows each contribution to the auto-correlation of the 2nd moment, P (2,2) , and 4th moment galaxy shape, P (4,4) , respectively, evaluated at z = 1.5. As the redshift becomes closer to z = 0, LA and l loop NLA in Fig. 2 change as D 2 (z) and D 2(l+1) (z), respectively. Since Fig. 2 shows results for the 3D power spectrum, it is mostly relevant for considering shape correlations in spectroscopic surveys, where the redshift space distortion is yet to be included. We will later perform the projection to angular correlations.

Intrinsic alignment from angular dependent PNG
As a candidate of new physics, in the rest of the paper, we consider the PNG generated by the higher spin particle excited in inflationary Universe. In Refs. [39,40,109], it was shown that the primordial bispectrum generated from the massive particles with non-zero spins exhibits the characteristic dependence on the angle between the wave vectors of the bispectrum. In this section, we show that such an angular dependent PNG generates an intrinsic distortion of the galaxy shape, extending the analysis in Ref. [54] to an arbitrary spin.

PNG from higher spin particles
Chen and Wang showed that the imprints of additional massive scalar particles can be encoded as an oscillatory contribution in the squeezed bispectrum [37]. This was extended to a general even spin particle by Arkani-Hamed and Maldacena in Ref. [39] and to a general integer spin particle, including odd spin particles, by Lee et al. in Ref. [40] (see also Refs. [41,110]). Such a higher spin particle is naturally predicted in string theory as a part of the infinite tower of the higher spin states [44] (see also Ref. [111]). Since the oscillatory feature in the PNG is characterized by the non-analytic scaling determined by the particle's spin s and mass M s , a precise measurement of the squeezed limit PNG may allow us to measure the mass and spin of the particles excited almost around the GUT scale. This program was dubbed the cosmological collider [39].
When the mass of the exchanged particle, M s , is much heavier than H inf , integrating out the massive particle yields a local contribution to the Lagrangian of the inflaton. Unless we know the original Lagrangian of the inflaton, we cannot selectively capture the information on the exchanged massive particle by analyzing the local contribution. On the other hand, when M s is of O(H inf ), the contribution of the particle exchange becomes non-local.
The squeezed bispectrum with the soft exchange of the even spin particles in the slowroll limit is given by [39] with k L /k S 1,k α ≡ k α /k α (α = L, S) and P s denotes the Legendre polynomials. For example, f s for a particle in the principal mass series is given by Reflecting the fact that the contribution of the massive particle is non-local, the bispectrum acquires the oscillatory contribution as a consequence of the non-analytic scaling. The contribution from the odd spin particles is cancelled in the leading order of k L /k S [39], but appears in the sub-leading order with (k L /k S ) 5/2 [40] (see also Ref. [109]).
In what follows, we show that the contribution of the spin-s particles to the PNG, described by Eq. (4.1), results in the generation of the intrinsic spin-s shape moment,Ĩ int i 1 i 2 ···is , keeping the spin-n shape moment with n = s intact as far as the non-linear evolution is negligible. To keep the generality, in the following, we use a generalized ansatz of the squeezed bispectrum (4.1) with s = 0, 1, 2, · · · , including the odd s and leaving f s as a general function of k L /k S . As was shown in Ref. [40], the sub-leading contribution from a spin-s particle includes P l with l = s. Therefore, at the sub-leading orders, the sth galaxy shape moment does not selectively single out the contribution of the spin-s particle.
When the function f s can be rewritten into the separable form as with a pivot scale k p , the squeezed bispectrum (4.1) modifies the local power spectrum at x as [95,96] where α L i 1 ···is is given by Here, we used Eq. (C.5) in App. C.1. These results are straightforwardly generalized to the case where the squeezed bispectrum is given by a sum of separable contributions. When f s is given by a linear combination of a power of k L /k S , e.g., Eq. (4.2), which is given by with ∆ s = 3/2 + iν s , each term can be rewritten as in Eq. (4.4).
In a simple setup, the squeezed PNG from massive particles is suppressed by the smallness of the coupling, the Boltzmann factor, e −Ms/T H with T H being the Hawking temperature T H = H inf /(2π), and the dilution factor (k L /k S ) 3/2 (see Eq. (4.2)). A larger coupling between the inflaton and the massive particles enhances the radiative corrections, which can lead to the break down of the EFT description. In Ref. [112], Kumar and Sundrum showed that the PNG mediated by heavy particles in several scenarios are too small to be detected, when we require the validity of the EFT description. They also showed that this difficulty can be circumvented in a curvaton scenario, where two different cutoff scales can be introduced. The Boltzmann suppression is not significant for a model where M s is around H inf as discussed in Refs. [43,109] by considering explicit models. In Ref. [113], Wang and Xianyu pointed out that a coupling between the inflaton and the massive sector is strictly restricted to evade the Boltzmann suppression without fine-tuning, since a large coupling leads to a large radiative correction to M s . This restriction has guided their attention to parity odd dimension 5 operators, which can relax the Boltzmann suppression [113,114] (see also Refs. [115,116]). Finally, the suppression by the dilution factor (k L /k S ) 3/2 for M s ≥ O(H inf ) can be evaded by lowering M s , while the straightforward attempt conflicts with the Higuchi bound for s ≥ 2 [117]. In Refs. [118,119], it was argued that the Higuchi bound and subsequently the suppression by the dilution can be avoided, when the de Sitter symmetry is explicitly broken. The model proposed by Kehagias and Riotto in Ref. [118] violates the global rotation symmetry.
In this paper, deferring an attempt to build a model that predicts a large angular dependent PNG to elsewhere, we simply take A s as a parameter. For our forecast, we consider both cases with and without the dilution, i.e. Re[∆ s ] = 3/2 and ∆ s = 0, respectively. The tightest constraint on the angular dependent PNG parameters has been obtained from the CMB observation by Planck satellite. The Planck result put the limit on A 2 for ∆ 2 = 0 as σ(A 2 ) ∼ 77 [50]. In Ref. [52], this analysis was extended to a more general angular dependence with even numbers of s, including the PNG generated by higher spin particles. However, the CMB constraint almost reaches the cosmic variance limit [120]. In Ref. [59], it was argued that from the bispectrum of the galaxy number density, A s can be constrained as σ(A 2 ) ∼ 15 for ∆ 2 = 0 and σ(A 4 ) ∼ 21 for ∆ 4 = 0. Ref. [60] discussed the constraint on A s for Re[∆ s ] = 3/2 by combining the galaxy power spectrum and bispectrum. These analyses require the bispectrum information, which is more complex than the power spectrum analysis, especially in galaxy surveys (e.g. the decomposition of anisotropic signals in the galaxy bispectrum [96,121,122]), and the angular dependent PNG signals are mixed up in the bispectrum. Here, one can find an advantage to use galaxy shape [54,56,70], whose observation enables us to pick up the imprints of particles with different spins separately as shown in this paper.

PNG contribution to galaxy shape function
Now, we are ready to calculate the imprint of the angular dependent PNGs encoded in the galaxy shape function.

Imprint of spin-4 particles
In Ref. [54], it was shown that the contribution of the PNG (4.1) in the second shape moment g i 1 i 2 , is only that generated by the spin-2 particles, as far as we consider the large scales where the non-linear secondary evolution is negligible. Here, we will show that similarly only the PNG generated by the spin-4 particle contributes to the 4th shape moment,g i 1 i 2 i 3 i 4 . Along the line with Ref. [54], we compute the contribution of the PNG tog i 1 i 2 i 3 i 4 by evaluating the correlation between δ andg i 1 i 2 i 3 i 4 .
In Sec. 3, various contributions to the 4th shape moment for the Gaussian initial condition were computed. The angular dependent PNG leads to an additional contribution to the correlation δ(x)g i 1 i 2 i 3 i 4 (y) . For example, when the PNG is given by Eq. (4.1) with ∆ 4 = 0, we obtain 3 ,sym and ξ δΦ being the cross-correlation between the linear δ and Φ. A more detailed computation of Eq. (4.8) can be found in App. C. Performing the Fourier transformation, we obtain where P m (k, z) is the linear matter power spectrum at z and M(k, z), which relates δ(k, z) to the primordial curvature perturbation Φ k as is given by with T (k) being the transfer function. Since Eq. (4.9) depends on δ 2 , which diverges when we send the UV cutoff to the infinity, we need to perform the renormalization to compute the observable effect. The renormalization proceeds analogously to the one for the PNG generated by the spin-0 and spin-2 particles discussed in Refs. [54,95,96]. As one can see in Eq. (4.5), the counter term for the contribution of the PNG generated by the spin-4 particle should be in the form, In the presence of the PNG generated by the spin-4 particle, the local matter density for the short mode k is modified as with The renormalized bias is then defined through the response of 4th moment to α b (4) which leads to The detailed computation is summarized in App. C.2. This result can be straightforwardly extended to the PNG with ∆ 4 = 3/2 ± iν 4 as Here and hereafter, we set the pivot scale k p to k * . For a different choice of k * , the corresponding b (4) NG and Θ 4 differ.

Imprint of particle with a general integer spin
We can also evaluate the imprint of the PNG generated by the spin-n particle with a general integer n, encoded in δ(x)g i 1 ···in (y) . The PNG from the spin-n particle selectively appears in the traceless part of the nth shape moment,g i 1 ···in . Therefore, even if A 0 , which corresponds to f NL , is much larger than A s with s ≥ 2, the contribution of A 0 does not contaminateg i 1 ···in with n ≥ 2 in the linear regime. Repeating a similar computation to App. C, we find that the spin-n contribution in the squeezed PNG (4.1) yields the additional intrinsic alignment contribution given by where F (n) PNG denotes the PNG contribution in Eq. (2.31) and Θ n is the phase determined for a given halo model, i.e. it is a function of k * . The leading PNG contribution appears from the quadratic terms in the bias expansion of the nth shape moment, which are accompanied with additional spatial gradient for n > 4. However, the spatial gradient does not yield an additional suppression by k/k * in Eq. (4.17), since it is replaced with (k S /k * ) n−4 , where k S is the short mode. Therefore, as written in Eq. (4.17), the leading PNG contribution scales as ∝ k Re[∆n] M −1 (k, z) for a general integer n. The coefficient C n is determined by conducting the renormalization. The counter term needed for the renormalization differs for a different n. For n = 0, 2, 4, C n is given by C 0 = 1/2, C 2 = 3 ([54]) and C 4 = 35/4 (Sec. 4.2.1).
Let us emphasize that the separability of the PNG contributions from different spins no longer holds, once the loop contributions become important, because the kernel functions also induce the angular dependence. Then, the PNG generated by the spin-s particles can contribute to the nth moment of the galaxy shape function only for n = s. Having considered this, in the next section, we explore whether there is a scale where the contribution from the PNG becomes dominant, keeping the late time non-linear contributions subdominant. ] NG A2 = 0.5 LA, b

Forecast on PNG from higher spin particles
In the previous section, we have computed the contribution of the PNG generated by a spin-n particle, to the nth moment galaxy shape function,Ĩ i 1 ···in . In this section, comparing it to other contributions predicted in ΛCDM cosmology with the adiabatic Gaussian initial condition, let us discuss whether we can observe the imprint of the higher spin particles from future observations.

Dominant contribution at different scales
Combining the results obtained in Sec. 3 and Sec. 4, we can evaluate which effect is dominant at each scale or at each multipole moment. Figure 3 compares the contributions of the PNG from the spin-2 particle to LA and NLA, computed in Sec. 3 for the 2nd galaxy moment. Here, considering the ΛCDM Universe, we set n L ∼ −1.7 and k NL,z=0 = 0.25h/Mpc [92]. For ∆ 2 = 0, the contribution of the PNG dominates at the large scales as pointed out in Ref. [54] (see also Ref. [70]). The colored region shows the range of k at which the PNG dominates the other contributions. Meanwhile, for Re[∆ 2 ] = 3/2, the contribution of the PNG is dominated by NLA at all scales with k < k eq , where the galaxy imaging surveys work as a spin-sensitive detector. When the perturbative expansion holds, satisfying b (2) NG A 2 10 4 (assuming b (2) NG = 0.1), there is no range of k where the PNG generated from the massive spin-2 particle with Re[∆ 2 ] = 3/2 dominates the 2nd moment of the galaxy shape function. The PNG from the spin-2 particle in the principle mass series, which is suppressed by the dilution, Re[∆ 2 ] = 3/2, exhibits the oscillatory resonance feature as shown in Eqs. (4.1) and (4.2). Since our purpose here is to compare the amplitudes of different contributions, for the illustrative purpose, the oscillatory contribution in Eq. (4.17) is ignored. The situation is different for the PNG generated from the spin-n particle with n > 2, because the linear alignment and the linear weak lensing forĨ i 1 i 2 ···in are both suppressed by (k/k * ) n−2 . Meanwhile, as shown in Eq. (4.17), the PNG can contribute toĨ i 1 i 2 ···in without being suppressed at large scales, especially for ∆ s = 0. Figure 4 compares the different contributions to the auto-correlation ofĨ i 1 i 2 i 3 i 4 . Since the dominant contamination of the 4th shape moment comes from NLA, i.e. the 3D loops under the Gaussian initial condition (even if we include the lensing contribution), here let us compare it to the PNG contribution. Around k NL , where the 1-loop NLA acquires the additional factor (k/k NL ) n L +3 , the PNG contribution dominates the NLA contribution in the range of k which satisfies Here, focusing on the amplitude, the oscillatory contribution (if any) is again ignored. On the other hand, for k < k eq , the leading contribution of NLA, which corresponds to the stochastic term, is roughly given by which is almost scale independent like the white noise. In practice, this contribution would be absorbed by the amplitude of the noise in the 4th moment, for which we do not have a reliable estimate. In our idealized setting then, the PNG contribution should exceed this to be detectable. The PNG contribution is z-independent while the linear alignment contribution varies as ∼ (1 + z) −2 and the NLA contribution depends on ∼ (1 + z) −4 . Therefore, a smaller amplitude of PNG can be detected for a deeper survey. This point is quantitatively analyzed in Fig. 7.

Angular power spectrum
Using Eq. (2.43), we can compute the angular power spectrum. Here, we focus on C (2,2) l and C (4,4) l , by which we can explore the imprint of the spin-2 and spin-4 particles, respectively. Inserting Eq. (4.17) into Eq. (2.43), we can compute the contribution of the PNG from spin-s particle to the angular power spectrum C (n,n ) l , where n and (or) n are (is) equal to s. The PNG contribution from the spin-2 particle was addressed in Ref. [54] (∆ 2 = 0) and in Ref. [56] (Re[∆ 2 ] = 3/2). In Fig. 5, we show the contributions of LA, NLA, WL, and PNG to C (2,2) l , including NLA, which was not taken into account in Refs. [54,56]. Assuming an LSST like lensing survey [1,123], we have used the function dN/dz given by with and Equation (5.6) can be derived by operating the source distribution function on the expression for a source at fixed redshift given in App. B 4 . Figure 6 compares different contributions to C  4 , which corresponds to the PNG from a spin-4 particle in the principle mass series. For a larger ν 4 , the net contribution of the PNG to C (4,4)PNG l becomes smaller, being smoothed out by the integration over k. This aspect is common with massive particles in the principal mass series with other spins, as discussed in Ref. [56].
Since the PNG contribution with Re[∆ s ] = 3/2 is almost constant, being independent of l, one may think it can be detected at high ls even if it is subdominant at lower ls. However, at high ls, because of the mixing due to the non-linear evolution, we cannot separately pick up the PNG contributions from particles with different spins. (Roughly speaking, using the redshift of the peak of dN/dz, we can estimate l NL ∼ k NL,z=0.5 χ(z = 0.5) ∼ 470 as the non-linear angular scale.) In the range where the non-linear evolution is negligible, while the PNG contribution does not become dominant for n = 2, it can be dominant for n = 4 in l 200 5 .

Forecast for future imaging survey
Finally, we investigate the detectability of the PNG generated by higher spin particles, considering the future galaxy survey, LSST. As discussed in Sec. 4.2, on large scales the PNG from each particle with spin-s only contributes to the sth galaxy moment of the galaxy shape function, leaving an imprint in the auto-correlation C (n, n) l with n = s and the crosscorrelations C (n, n ) l with n = s or n = s. Therefore, no matter how large the angular independent PNG f NL is, this does not disturb probing the signal of the higher spin particle encoded inĨ i 1 ···in with the corresponding value of n (see [55] for a demonstration including f NL and A 2 ). Furthermore, even if there exists an infinite tower of higher spin particles as predicted in string theory, we can single out the contribution of the spin-s particle by looking at the sth galaxy shape moment.
On the other hand, when the contribution of the loops becomes important, the nontrivial momentum dependence in the kernel functions allows the contribution of the PNG from the spin-n particle to contaminateĨ i 1 ···in with n = n. Therefore, in what follows, we focus on k < k NL , where the loop contribution remains subdominant.
As discussed in Sec. 3, since the dominant contamination to the s = 4 signal at k < k NL comes from NLA, in the following, we ignore the contributions of the linear alignment and weak lensing. The signal to noise ratio (S/N) for each redshift, which is not cumulative, is expressed by where f sky is the sky fraction to the galaxy imaging survey and C (4,4)NLA (l, z) denotes the NLA contribution which is crudely estimated by using Eq. (3.15). Here, we have discarded the contribution from the non-linear scales, setting l max = k NL (z)χ(z). For ∆ 4 = 0 (without dilution), S/N becomes independent of l max when we choose a sufficiently large l max , since the shape signal of PNG is significant at the large scales. On the other hand, for Re[∆ 4 ] = 3/2 (with dilution), since the signal dominates the shape noise only at the small scales, S/N depends on the choice of l max . However, since we cannot selectively single out the contribution of the spin-4 particle, we expect that increasing l max will not improve the actual S/N very much.
NG A 4 which is required for a detection with S/N = 5 by a survey covering a fraction f sky = 0.5 of the sky when including scales up to l max = k NL (z)χ(z). The blue solid line corresponds to ∆ 4 = 0 and the orange dotted line corresponds to ∆ 4 = 3/2. Using Eq. (5.8) we can compute the corresponding amplitude of the PNG which can be detected for a given S/N. Figure 7 shows the amplitude of b (4) NG A 4 which is detectable with S/N = 5 for a survey covering a fraction f sky = 0.5 of the sky up to redshift z. The increasing volume and number of modes available at higher redshift leads to improved constraints. Further, as discussed in Sec. 5.1, since the NLA contribution is smaller at higher redshifts while the PNG contribution remains constant, we find that the high redshift survey works in favor of PNG detection.

Conclusion
In this paper, we have explored the possibility that a new physics can be searched through future imaging surveys. For this purpose, we have established a way to compute the angular power spectrum for the nth moment galaxy shape function, performing the spin decomposition. This generalizes the conventional analysis for the cosmic shear, the spin-2 component, to an arbitrary spin component of the galaxy shape function.
In order to look for a signal of new physics, we need to understand the prediction in the concordance cosmology. We estimated the contribution to the nth shape moment through the linear and non-linear alignment and also the linear and non-linear weak lensing, assuming the ΛCDM cosmology with the adiabatic Gaussian initial condition. This provides a benchmark value which the signal of the new physics under consideration should exceed. For the 2nd moment, the usual cosmic shear, the dominant contamination to the signal of the new physics comes from the linear weak lensing and the linear alignment. Meanwhile, for the higher shape moment with n > 2, since these linear contributions are suppressed by a positive power of k/k * 1, the dominant contamination comes from loop contributions through the non-linear alignment. We have estimated these loop contributions, based on the effective field theory of large scale structure.
As an example of new physics, which can be searched through galaxy imaging surveys, we have studied the angular-dependent primordial non-Gaussianity, which can be generated through a non-zero integer spin particle during inflation. We have shown that an ideal galaxy imaging survey, which can measure the different spin components of the galaxy shape function independently, enables us to detect the imprints of particles with different spins separately, playing the role of a spin-resolved detector of the cosmological collider. This separation only works at smaller wavenumbers than k NL , since the angular dependence of the kernel functions mixes up contributions of different spins.
We have considered two different scaling behaviors of PNG, with and without the dilution, corresponding to ∆ 4 = 0 and Re[∆ 4 ] = 3/2, respectively. We have found that for ∆ 4 = 0, an LSST-like survey can detect the PNG of O(1) generated by a spin-4 particle. This is comparable to the expected constraints on local primordial non-Gaussianity, which probes the existence of an additional scalar (spin-0) particle. Meanwhile, detecting the PNG for Re[∆ 4 ] = 3/2 is rather challenging, mainly because the signal grows towards smaller scales (larger k) and is difficult to distinguish from nonlinear contributions arising from the dominant Gaussian perturbations.
In an actual observation, of course we also need to take into account the shape noise, which is ignored in this paper. In case of the 2nd moment, the amplitude of shape noise is well known, while we are not aware of reliable published shape noise estimates for the 4th moment studied here. Instead, we have used the constant, k 0 , limit of the nonlinear alignment loop contribution as a rough estimate, which is likely to underpredict the actual noise in the 4th moment. Furthermore, in order to probe the actual amplitude of the PNG, one needs to constrain the shape bias b (4) NG for the 4th moment, as the observed correlation functions only constrain b (4) NG A 4 . We can approach this issue from both an analytic (e.g. the peak theory [124]) and numerical approaches. In particular, using N -body simulations we can estimate the response of the 4th moment of halo shapes to PNG in analogy to Ref. [70]. These issues will be addressed in future work.

A Formulae
Here, we derive several formulae used in the main text.

A.1 Derivation of Eq. (2.41)
Using mathematical induction, we prove that for all integers s ≥ 0, We can explicitly confirm P (s) for s = 0, 1, 2. Given that P (s) is true for s = k, k − 1 (k ≥ 2), taking the derivative of P (k) with respect to x, we obtain where we used The left hand side of P (k + 1) is given by Using P (k), Eq. (A.1) and the Bessel differential equation: Thus, P (k + 1) is true, whenever P (k) and P (k − 1) are true. Hence, by the Principle of Mathematical Induction, P (s) is true for all integers s ≥ 0.
Then, we get

A.3 Symmetric traceless tensor in arbitrary dimension
Next, we consider an nth symmetric traceless tensor in d-dimension, expressing it as where I i 1 i 2 ···in is an arbitrary rank-n symmetric tensor andĨ i 1 i 2 ···in is its traceless part. In the square brackets of the right hand side, the indices I i 1 i 2 ···in are replaced with the contraction byĩ m (m = 1, · · · , k) and the replaced indices appear as the indices of the Kronecker delta. Imposing the traceless condition, we obtaiñ For n ≥ 2, d ≥ 2 and 0 ≤ k ≤ [n/2], 2k − (d + 2n − 4) ≤ 0. Then, we can obtaiñ For example, the coefficient for d = 2 is given bỹ .14) and the one for d = 3 is given byC , (A.23) where we have changed k − l → k. Rewriting the summation over k by using the following formulae: we arrive at Eqs. (A.18) and (A.19).

B Weak lensing
As discussed in Sec. 3, the galaxy shape moment also can be generated through the gravitational lensing. In this Appendix, we estimate the spin-n component of the distortion generated through the weak lensing,Ĩ WL i 1 i 2 ···in (θ), providing the detailed computation. As discussed around Eq. (2.2),Ĩ WL i 1 ···in is given by the difference between the observed galaxy shape function, I i 1 i 2 ···in , and the intrinsic one, I int i 1 i 2 ···in . The former is defined by using the coordinates on the image plane,θ i =θ i + θ i , and the latter is defined by using those on the source plane,θ s i =θ s i + θ s i . In Sec. B.1, solving the geodesic equation, we derive the relation between these two coordinates,θ i =θ i (θ s ). We put a bar on the centroid coordinates and a tilde on the coordinates measured from the common origin on each plane. Using this relation, we can express the deviation from the centroid on the image plane, θ i , using the coordinates on the source plane (→ Eq. (B.14)). Using this expression, in Sec. B.2, we computeĨ WL i 1 i 2 ···in (θ).

B.1 Solving the null geodesics
The relation between the intrinsic galaxy shape on the source plane and the apparent galaxy shape on the image plane can be computed by tracing the photon propagation in the perturbed FLRW spacetime. Using the affine parameter λ, the null geodesic equation, which describes the photon propagation, is given by where Γ µ αβ is the Christoffel symbol of the perturbed spacetime. The Greek indices µ, α, β run from 0 to 3. Solving the geodesic equation along the orthogonal direction ton with a use of the null condition we obtain the lens equation as Here, we have imposed the boundary condition asθ i (χ = 0) =θ i and θ i (χ) =θ i s . The projected coordinatesθ i andθ s i correspond to those on the image plane and the source plane, respectively. A detailed computation of Eq. (B.3) can be found e.g. in Ref. [71]. We also useθ s i to express the mapping of a coordinate on the image planeθ i to the corresponding one on the source plane asθ s i (θ).
Taking the partial derivative ofθ s i with respect toθ j , we obtain the deformation matrix A ij [108] as Because of the spatial inhomogeneity of Φ and Ψ, the deformation matrix A ij depends onθ.
These coordinates on the source plane and the image plane are visually explained in Fig. 8. The first square brackets describe the difference between the two coordinates on the source plane which are mapped fromθ i andθ i on the image plane, respectively. The second square brackets describe the difference between the coordinates on the source plane mapped from the apparent centroid,θ i , and the actual centroid,θ s i . Since we cannot directly measure the source distribution, the actual centroidθ s is left undermined. Using Eq. (B.3), we can compute the first square brackets in Eq. (B.5) minus θ i as Here we are evaluating Φ and Ψ along the unperturbed photon path. This is only correct at linear order in the deflection. If one includes the lower-order deflection in the argument of these metric perturbations and expands this effect as well, one obtains additional contributions which involve the integral over products of ∂ k ∂ l Φ and ∂ k ∂ l Ψ. These contributions, which are known as "post-Born corrections" (e.g., [125]), to the nth moment can be roughly approximated as being of order (δA (2) ) m , with m = n if n is even and m = (2n + 1)/2 if n is odd. This means that they are of the same order as the nonlinear lensing contribution we will discuss below, and hence we will not derive them in detail here.
Expanding δθ i with respect to θ i =θ i −θ i , we obtain The spatial variation of A ij yields the higher order terms of θ i , characterized by δA ij is nothing but the deviation of the deformation matrix from the unit matrix. The traceless part of δA (n) ii 1 ···i n−1 describes the spin-n deformation due to the gravitational lensing, whose visual image is discussed in Sec. 2.2.

B.2 Estimation of weak lensing effect
Using the formula derived in the previous subsection, let us compute the weak lensing contribution in Eq. (2.2), which can be rewritten as Similarly to Eq. (2.8), using the projected coordinates on the source plane, θ s i , the n-th moment intrinsic galaxy shape function is defined as The brightness theorem, which states that the lensing does not change the surface brightness, relates the apparent surface brightness I(θ) to the intrinsic surface brightness I int (θ s ) as I int (θ s ) = I(θ) or equivalently in our notation, First, let us compute the first square brackets of Eq. (B.10) by using θ i expressed in terms of θ s i . In the previous subsection, we derived the expression of δθ i (θ), solving the null geodesics. Using this expression, θ s i in the left hand side of Eq. (B.5) can be expressed in terms of θ i and the global shift between the apparent and actual centroids. Solving Eq. (B.5) recursively, we can express θ i in terms of θ s i and the contribution due to the global shift of the centroid as where, e.g., the leading contributions of δB (2) i 1 i 2 and F i , which is a function ofθ s (θ) −θ s , are given by δB (2) Since δA The leading lensing effect on the nth moment of a galaxy image can be computed by making the lensing operators derived above act on a perfectly symmetric (spherical) intrinsic image, whose moments are given by where N is a natural number and C 2N is of order the angular size of the image to the nth power; Ref. [76] also presents the effects of non-symmetric intrinsic images. In the following, we calculate the terms in the first brackets of Eq. (B.10) for n = 2N . For the moment, let us ignore the contributions due to the shift of the centroid, which turns out to be sub-leading. Inserting Eq. (B.14) into Eq.
with M = 1, 2, · · · . The factor (R * /χ) 2(M −1) of the first contribution appears, since I i 1 ···i 2N and I int i 1 ···i 2(N +M −1) carry different powers of R * /χ. In computing the terms in the second line, which is given by 2N product of δA (2) , we used δA (2) ipj δA (2) iqj P jj = δA (2) ipj δA (2)j iq . The ellipses stand for the trace part and other non-linear terms in perturbations. As discussed above, these include post-Born corrections which are of similar order of magnitude. While we did not write it explicitly, for M < N , the lacking indices should be supplied by the projection tensor P ij . Therefore, such terms do not contribute to the traceless component. Since the order of δA In Eq. (B.18), we have replaced ∂/∂θ i with kχ , where k is the Fourier mode of the 3D coordinates of the centroid on the image plane. Equation (B.19) indicates that at the linear order of the perturbation, the traceless component of the first brackets of Eq. (B.10) with n = 2N ≥ 4 is suppressed at least by (kR * ) n−2 , which is much smaller than 1 since we only consider larger scales than the typical size of galaxies. To be precise, the integrand of δA (2) in Eq. (B.19) is further suppressed by (χ /χ) 2(M −1) . For example, for the 4th moment, the tree-level diagram which appears by contracting the first contribution in Eq. (B.17) roughly scales as Meanwhile, for the 2nd moment with N = 1, the linear contribution simply gives the deformation matrix δA (2) i 1 i 2 . Similarly, all the non-linear contributions which include δA (n) i 1 ···in with n ≥ 3 are suppressed by positive powers of kR * . Therefore, the only contributions which are not suppressed by kR * 1 are products of δA (2) , i.e. the second contribution in Eq. (B.17). In the end of this subsection we will discuss the loop contributions from the products of δA (2) .
So far, we have not considered the contribution from the shift of the centroid. According to Ref. [81], the centroid shift can be estimated by (see also Ref. [126]) Repeating a similar argument, we find that the leading contribution of the centroid shift has more (kR * ) than the leading contribution and is suppressed more on the large scale. For this reason, we do not consider the contribution of the centroid shift. Next, let us discuss the first brackets of Eq. (B.10) for an odd number n. Repeating the same argument, we find that the odd nth moment is always suppressed by kR * , since the intrinsic nth moment of the galaxy shape function I int i 1 i 2 ···in (θ s ) vanishes for an odd n under the assumption of Eq. (B.16). In particular, the first brackets of Eq. (B.10) is suppressed at least by (kR * ) n−2 at the linear perturbation.
In summary, we showed that the only contributions that are not suppressed by kR * ( 1) are the products of δA (2) . A contraction of δA (2) yields weak lensing loops, computed in the projected 2D plane. As emphasized in Sec. 3.2, these 2D loops are qualitatively different from 3D loop contributions to the galaxy shape moments, which are projected into 2D after computing loop contributions in 3D (see, e.g. Eq. (3.15)). The 3D loop contributions are included as the non-linear alignment effect. For a comparison of 2D and 3D loops, here we estimate their contributions to the angular power spectrum of the 4th shape moment, C (4,4) l . Taking the flat-sky limit, the 2D weak lensing loops in C (4,4) l is roughly estimated as, C (4,4)2D-1L l ∼ (δA (2) δA (2) )(l)(δA (2) δA (2) where C γ (l) is the lensing shear for 2nd moment and l is a 2D vector on the Cartesian coordinate. We find that the contribution of 2D lensing 1-loop, given by Eq. (B.22), is much smaller than the one of the 3D 1-loop(NLA). Therefore, evaluating the forecast in this paper, we neglect the weak lensing contribution, whose linear contribution is suppressed by kR * and whose loop contributions are much smaller than the 3D loops, included in non-linear alignment effects.
In this Appendix, we have estimated the nth moment deformation of the galaxy shape due to the weak lensing, assuming the absence of the intrinsic deformation, i.e. assuming a circular intrinsic image. As was argued in Sec. 4.2, the angular dependent PNG can generate the intrinsic alignment of the galaxy shape. To estimate the weak lensing contribution more accurately by taking into account the intrinsic deformation, the assumption, (B.16), should be abandoned.

C PNG from spin-4 particle
In this Appendix, we compute the contribution of the PNG generated by a spin-4 particle, given in Eq. (4.1), to the 4th moment galaxy shape function.

C.1 Contributions to galaxy shape function
From Eq. (3.10), the cross-correlation of the matter density field withg ijkl is Let us compute δ(x) [K ij (y)K kl (y)] TL 3 , sym term in the presence of the non-Gaussian initial condition such as Eq. (4.1).
In the first line, we used the symmetric tensor property and in the second line, we used Eq.(4.1) and expanded B φ (k, p, |k + p|) and M(|k + p|) in the powers of q 1. We also applyp ∼p in the powers of q 1. Choosing the z direction of the polar coordinates along the direction ofk and integrating over the azimuthal angle ϕ, we obtain This can be understood by noticing that since the left hand side of Eq. (C.4) should be independent of ϕ, all the tensor indicies should be along the direction ofk, satisfying the symmetric traceless condition. The amplitude of the left hand side of Eq. (C.4) can be determined as follows. As shown in Ref. [127] (in this paper, the spatial dimension is d, so here we set d = 3), the contraction between [p i 1 · · ·p in ] TL 3 and [k i 1 · · ·k in ] TL 3 is given by Operating [k i 1 · · ·k in ] TL 3 on the left hand side of Eq. (C.4) and using Eq. (C.5), we find that the amplitude of the left hand side of Eq. (C.4) should be P n (µ), i.e.
Here, we used Using Eq. (C.6) in Eq. (C.3), we obtain where we haved used the orthonormality of the Legendre polynomials.

C.2 Renormalization
Next, we define the renormalized bias parameter b (4) NG following the discussions in Refs. [54,128].

C.2.1 Gaussian initial conditions
First, we consider the case of Gaussian initial conditions. Let us introduce the coarse-grained density field δ L and tidal field K L,ij with a coarse-graining scale R L . Then, using a functional F L,ijkl , we can formally expressg ijkl (x) as g ijkl (x) = F L,ijkl (δ L (x); K L,pq (x); δ s (x)), (C.9) with δ s (x) ≡ δ(x) − δ L (x) being the small-scale fluctuations on which in principle the 4th moment depends other than δ L and K L,ij . A formal Taylor expansion of F L,ijkl in δ L and K L,ij leads tog K 2 (R L ; δ s (x)) [K L,ij K L,kl ] TL 3 , sym (x) + O(δ 3 L , ∇ 2 δ L ), (C. 10) In general, the coefficient c K 2 depend on on the short modes δ s and the coarse-graining scale R L 7 . In the Gaussian case, however, the short modes and long modes are uncorrelated, so the c K 2 can be regarded as an effective constant. The spherical symmetry requires the expectation value ofg ijkl (x) to vanish. Let us consider the following modification of the tidal field: K L,ij (x, τ ) → K L,ij (x, τ ) + D(τ )β ij ; δ L (x, τ ) → δ L (x, τ ), (C.11) or equivalently the modification of the Newtonian potential, where β ij is a constant symmetric traceless tensor. One can interpret this as the leading observable effect of a potential perturbation in the k → 0 limit, as constant and pure-gradient potential perturbations can be removed by coordinate transformations. Alternatively, this effect can be realized in simulations by implementing an anisotropic expansion, roughly resembling a Bianchi I spacetime [129,130]. Under this modification, the expectation value ofg ijkl changes to Notice that g ijkl (x) does not vanish owing to the presence the preferred direction β ij . Then, the renormalized bias parameters are introduced via b (4) (C.14) Note that here the summation is not taken. This definition of the bias parameter b (4) K 2 is independent of R L by construction. The point is that the renormalized bias is defined as the response to the locally uniform transformation of the tidal field (second derivatives of the potential).

C.2.2 Non-Gaussian initial conditions with no scaling (local-type)
In the presence of the primordial non-Gaussianity (4.1), Eq. (C.8) yields δ(x) [K ij K kl ] TL 3 , sym (y) = 2 9 A 4 D ijkl ξ δΦ (|x − y|) δ 2 L . (C. 16) This expression strongly depends on the coarse-graining scale through δ 2 L , so the bias expansion Eq. (C.15) is not sufficient for non-Gaussian initial conditions. In this case, we have to take into account the dependence of g ijkl on the small-scale fluctuations δ s explicitly, since the primordial non-Gaussianity couples the long-modes with short-modes. To do this, let us introduce the parameter y ijkl where W L is an isotropic window function of the scale R L , K ij s (x ± r/2) ≡ D ij r δ s (x ± r/2) with D ij r being the derivative operator acting on r and W s (k) ≡ 1− W L (k). We now introduce explicitly the dependence of g ijkl (x) on y ijkl s (x): g ijkl (x) = F L,ijkl (δ L (x); K L,pq (x); y pqrs s (x)). (C. 19) Expanding F L,ijkl to linear order in y pqrs s leads to adding a term c (4) NG y ijkl s (x) (C. 20) to the expansion on the r.h.s. of Eq. (C.10). Then, we have additional contribution to matter-shape correlation, δ(x)g ijkl (y) =b NG should consist of R L -independent combination of these contributions. In other words, physically, b (4) NG should correspond to the response of the 4th moment to a specific (R L -independent) transformation of the density field as implied in the previous subsection.
In fact, under the initial condition given by Eq. (4.1) the local power spectrum is modulated like Eq. (4.5), which means that the long-wavelength potential perturbation leads an anisotropic modulation of the local initial matter power spectrum described by Note that we treat α L,pqrs as a locally constant number. Here, we assumed f 4 (k L /k S ) 1. This implies that in terms of the density field the presence of the anisotropic non-Gaussianity alters the local density field such that δ α (k S ) = 1 + 1 2 α L pqrs k p This is an anisotropic (hexadecapole), scale-independent rescaling of the density field. After this transformation, the expectation value of the 4th moment galaxy shape changes to To proceed the computation of Eq. (C.31), we use the following identity: where [δ i 1 i 2 · · · δ i n−1 in ] sym means to symmetrize the expression δ i 1 i 2 · · · δ i n−1 in in the indices i 1 · · · i n . For instance, the explicit expression for n = 4 case is given by This tells us that the linear response of the mean 4th moment of galaxies, through which we define the renormalized bias b (4) NG , is given by NG .

C.2.3 Non-Gaussian initial conditions with scaling
In this case, the modulation of the local initial matter power spectrum becomes P ini m,α (k S ; k L k p k L,pkL,qkL,rkL,s TL 3 Φ(k L )e ik L ·x . (C. 40) In terms of the density field, the presence of the anisotropic non-Gaussianity alters the local density field such that Obviously, in this case the local density field is modulated with the scale-dependence, unlike (C.28) where the local density field is rescaled uniformly. In this case, defining the renormalized bias with respect to the locally uniform modulation leads to the final expression in the Fourier space,