text:
# Galaxy Spin Statistics and Spin-Density Correlation

## 1 INTRODUCTION

The origin and evolution of the galaxy angular momentum, i.e., the galaxy spin has been the subject of many studies in the last century. Hoyle (1949) suggested an original idea that the origin of the rotational galaxy motion could be ascribed to the gravitational coupling with the surrounding galaxies. Sciama (1955) applied Hoyle’s idea to his theory for the formation of galaxies in a steady-state universe model.

It was Peebles (1969) who first quantitatively examined Hoyle’s idea in the gravitational instability picture. He argued that the shear effect due to the primordial tidal torquing from the neighbor matter distribution should be mainly responsible for the acquisition of the angular momentum by a proto-galaxy. He pointed out that the alternative models for the origin of the galaxy angular momentum such as the initial vorticity model and the primeval turbulence model proposed by von Wiezsacker (1951) and Gamow (1952) to a wrong prediction of too early formation of galaxies. Assuming a spherical symmetry of a proto-galaxy, he analyzed quantitatively the growth rate of the magnitude of the galaxy angular momentum in the frame of the linear perturbation theory, and drew a conclusion that the galaxy angular momentum grows as proportional to the second order perturbation ($`t^{5/3}`$ for a $`\mathrm{\Omega }=1`$ universe).

In contrast, White (1984) showed that the proto-galaxy angular momentum grows at first order ($`t`$ for a $`\mathrm{\Omega }=1`$ universe) unless the restrictive condition of the spherical symmetry is imposed on proto-galactic sites, which was originally contended by Doroshkevich (1970). He expanded Doroshkevich’s contention in detail by means of the linear perturbation theory described by the Zel’dovich approximation, and confirmed that the proto-galactic angular momentum is generated by the misalignment between the proto-galactic inertia tensor and the local gravitational shear tensor, and grows to first order during the linear phase. He confirmed his results by N-body simulations.

Heavens & Peacock (1988) analyzed the correlation of the galaxy angular momentum with the local density maxima in the linear regime, and concluded that the total angular momentum in the linear regime is almost independent of the height of the density peaks (see also Hoffman 1986, 1988). Catelan & Theuns (1996) extended the Heavens-Peacock works and calculated the expectation value of the angular momentum assuming an ellipsoidal proto-galaxy centered on a peak of the Gaussian density field using White’s formula.

While all these studies concentrated on the magnitude of the angular momentum, the total angular momentum, or even the fraction of virial energy in rotation, is very difficult to observe. On the other hand, the direction of the angular momentum, i.e., the galaxy spin axis can be measured only from the position angle on the sky and the projected axis ratio, which can be implemented for very large surveys. Therefore, the galaxy spin axis could provide more useful statistics that can be easily tested against real observational data.

Recently Lee & Pen (2000, hereafter LP00) pointed out that the 1st-order linear perturbation theory predicts preferential alignments of the galaxy spin axis along with the 2nd principal axis of the local gravitational shear tensor, and suggested a unique statistical model that uses the galaxy spin axis as a tool to reconstruct the initial density field. Their theory is based on two basic assumptions: First, the spin axis of a galaxy aligns well with that of the underlying dark halo. Second, the galaxy spin aligns with the 2nd principal axis of the local gravitational shear tensor to a detectable degree.

The first assumption is generally accepted as a reasonable working hypothesis for the spiral galaxies. A spiral galaxy is a highly flattened disk with its plane perpendicular to the direction of the underlying halo spin axis in most galaxy-formation theories (e.g., Mo, Mao, & White 1998). There is also an observational clue to this assumption. The galaxy spin can be observed to much larger radii than the galaxy radius through radio emission of the gas, and the spin direction of a spiral galaxy has been seen to change only very modestly as one moves to larger radii. This suggests that the spin axis of a galaxy is well correlated with that of the whole halo.

Meanwhile the second assumption should work subject to numerical testing. In fact several N-body experiments have already shown that the dark halos have preferred direction in the spin orientation, and that this preferential spin alignment has likely a primordial origin. Dubinski (1992) for the first time found this preferential spin alignment in N-body simulations of dark halos. By comparing simulations in the absence and presence of a cosmological tidal field, he investigated the effect of the initial shear on the spin orientation. He concluded that the numerically detected preferential spin alignment measured in dark halos resulted from the shear effect due to the linear regime tidal torque. LP00 directly calculated the correlation in direction between the Lagrangian and the Eulerian angular momentum of dark halos measured in N-body simulations, and confirmed that the linear theory prediction for the orientation of the halo angular momentum is quite a good approximation.

The central concept of LP00 is that one can use this preferential spin alignment, if it really exists, as a linking bridge between the initial matter distribution of the universe and the observable unit galaxy-spin field. They have provided a mathematical algorithm to reconstruct the initial shear and density fields from the observable galaxy spin axes. Conventionally the peculiar velocity or the weak lensing shear fields have been used to reconstruct the total mass density field. This new method for the density reconstruction using the galaxy spin field is believed to be advantageous for a couple of reasons: First, the galaxy spin axis is relatively easier to measure observationally. Second, it is free of the standard galaxy biasing. Third, it allows the reconstruction of the full three dimensional density field.

Very recently, numerical study of the galaxy spin or galaxy ellipticity alignment due to the local gravitational shear has become quite topical. The flurry of recent activities (Croft & Metzler 2000; Heavens, Refregier, & Heymans 2000; Catelan, Kamionkowski, & Blandford 2000; Crittenden et al. 2000a; Crittenden et al. 2000b) is motivated partly by the statistical search in blank fields for weak lensing signal, for which the intrinsic galaxy alignment plays a role of systematic error. The gravitational shear effect on the galaxy spin axis due to the initial tidal torquing is local, distinguished from the weak lensing shear which can change only the apparent orientation of the galaxy spin axis due to the distant intervening matter far from both the source and the observer. From here on, the intrinsic cosmic shear is referred to as the local gravitational shear.

Heavens et al. (2000) calculated the intrinsic ellipticity correlation of dark halos found in high-resolution N-body simulations. They implied that at small redshift the intrinsic ellipticity correlation due to the local cosmic shear effect dominates the correlation signal due to the weak lensing effect. Crittenden et al. (2000a) reanalyzed the results of Heavens et al. and demonstrated that the results of Heavens et al. from high-resolution N-body simulations in fact indicate stronger intrinsic spin alignments than that of LP00 from low-resolution simulations. Croft & Metzler (2000) also detected the intrinsic correlation of projected ellipticities of dark halos in high-resolution N-body simulations, and also showed that the correlation signal is not strongly affected by the resolution of the simulations. Actually they found by comparing two simulations of different resolutions that the simulations of higher resolution found more intrinsic correlations.

Observationally, the question of galaxy spin alignment has received periodic attentions. The history of observational search for galaxy alignment traces back to the 19th century, and has been marked by checkered records (Strom & Strom, 1978; Gregory et al., 1981; Binggeli 1982, ; Helou & Salpeter, 1982; Helou, 1984; Dekel 1985, ; Lambas et al., 1988; Flin 1988, ; Hoffman et al., 1989; Kashikawa & Okamura 1992, ; Muriel & Lambas, 2000; Godlowski, 1994; Han et al., 1995; Cabanela & Dickey, 1999). For a review of the history of the field, see Djorgovski (1987) and Cabanela & Aldering (1998). However, past observational searches for galaxy alignment suffered from the small sample sizes Cabanela & Dickey (1999). It was only very recently that positive and reliable signals of galaxy alignments have been detected from large galaxy samples. Pen, Lee, & Seljak (2000, hereafter PLS00) have reported a tentative detection of the intrinsic spin correlation signal from the Tully galaxy catalog. The observed signal turns out to be significant at the $`97\%`$ confidence level with the amplitude of order of $`1\%`$ at $`1h^1`$ Mpc, which is consistent with the theoretical predictions made by PLS00. Brown et al. (2000) also detected the intrinsic alignment in galaxy ellipticities using the SuperCOSMOS Sky survey data. They showed that their results agree well with the linear theory predictions on the galaxy preferential alignment (Crittenden et al. 2000a; LP00). The observed and simulated amplitudes of correlations are expected to be stronger than the weak lensing effects for surveys such as Sloan Digital Sky survey (SDSS), and thus a quantitative analysis of intrinsic alignments must be completed before one can attempt to measure weak lensing shears within SDSS. These positive observational results hint at a possible detection of the spin-density cross correlation signal which will be addressed here.

The theory proposed by LP00 that the linear shear and density fields can be reconstructed using the detectable intrinsic galaxy spin alignments is quite speculative, based on many simplifying assumptions. The idea of LP00 must go through thorough observational and numerical testings in the future. However, the recent observational detections of intrinsic galaxy alignment and the agreement of the strength of the observed signals with the theoretical predictions encourage us to have a prospect for the plausibility of our theory and its application to the real universe. If, as predicted, the intrinsic alignment signal indeed dominates the weak lensing signal in shallow surveys like SDSS, the extraction of intrinsic shear becomes more plausible than that of the weak lensing.

In this paper, we present the galaxy spin correlation statistics with technical details. In $`\mathrm{\S }2`$, we review the mathematical algorithms for the density reconstruction given by LP00 in greater detail for the reader’s thorough understanding of our previous and future works. In $`\mathrm{\S }3`$ we review the spin-spin correlation statistics, and provide an analogous spin-density correlation statistics. In $`\mathrm{\S }4`$ we compare the theoretical estimates given in $`\mathrm{\S }3`$ with the observed signals. In $`\mathrm{\S }5`$ the results are summarized and final conclusions are drawn. We relegate the detailed calculations and derivations to Appendices A - J.

## 2 DENSITY RECONSTRUCTION

In the standard gravitational instability picture, a proto-galaxy acquires its angular momentum from the local gravitational shears due to the tidal interaction with the surrounding matter. The angular momentum of this proto-galaxy gradually evolves till the proto-galactic region reaches the moment of recollapse. On recollapse, separated out from the rest of the universe, its angular momentum would be approximately conserved afterwards. In other words, the galaxy angular momentum is expected to preserves its initial dependence on the local shear tensor fairly well that has been acquired during the linear regime. It is worth mentioning that the galaxy merging or secondary infall does not break the dependence of the galaxy angular momentum on the initial shears since the total rotational angular momentum after merging or infall process is the result of the constituent orbital angular momentum of the galaxies combined, which depends on the initial shear tensors. Similarly, the impact parameter of a collision is also determined by the shear field. Thus, what changes by those processes is only the smoothing scale of the intrinsic shear which correlates with the galaxy angular momentum.

It is true that one cannot expect the linear theory to fully describe the evolution of the galaxy angular momentum. Nonlinear effects such as galaxy-galaxy interaction and etc. may modify the galaxy angular momentum during the subsequent evolutionary stages. Nevertheless recent numerical simulations have found that in fact the linear theory predictions for the direction of galaxy spins are in fairly good agreement with numerical results (Dubinski 1992; LP00). Thus, we base our study of the direction of the galaxy angular momentum on the linear perturbation theory.

White (1984) and Catelan & Theuns (1996) have shown that in the 1st order linear perturbation theory described by the Zel’dovich approximation, the galaxy angular momentum in Lagrangian space is expressed as

$$L_i(t)=S^2(t)\frac{dD(t)}{dt}ϵ_{ijk}T_{jl}I_{lk},$$
(1)
where $`S(t)`$ is the expansion factor, $`D(t)`$ describes the growing mode of the density perturbations, $`𝐈=(I_{lk})=(q_lq_kd^3𝐪)`$ is the inertia tensor of a proto-galactic site in Lagrangian space, $`𝐓=(T_{jl})=(_j_l\varphi )`$ is the local shear tensor defined as the second derivative of the gravitational potential, $`\varphi `$ smoothed on a galactic scale of $`R`$.

Rotating the frame into the principal axis of the local shear tensor, $`𝐓`$, we can reexpress equation (1) in terms of the three eigenvalues, $`\lambda _1,\lambda _2,\lambda _3`$ of $`𝐓`$ such that

$$L_1(\lambda _2\lambda _3)I_{23},L_2(\lambda _1\lambda _3)I_{31},L_3(\lambda _1\lambda _2)I_{12},$$
(2)
where the three eigenvalues are ordered to be $`\lambda _1>\lambda _2>\lambda _3`$. We note three important implications of equation (2). First, if a proto-galactic region is spherically symmetric (corresponding to $`I_{12}=I_{23}=I_{31}=0`$), then the region gains no angular momentum at 1st order. Second, if the principal axis of the inertia tensor, $`𝐈`$, is aligned perfectly with that of the shear tensor, $`𝐓`$, then the off-diagonal elements of the inertia tensor is zero in the shear principal axis frame, no angular momentum generated at 1st order, either. Third, if the proto-galactic region is non spherical, and the principal axis of the inertial tensor is misaligned with that of the shear tensor, then one can expect the region to acquire a net angular momentum vector with $`L_2`$ being dominant since $`\lambda _1\lambda _3`$ is always bigger than the other two differences. In other words, the direction of the proto-galactic angular momentum is on average preferentially aligned with the 2nd principal axis of the shear tensor.

For the ideal situation where the principal axis of the inertia tensor is totally independent of that of the shear tensor Catelan & Theuns (1996), one can expect the maximal preferential alignment of the galaxy angular momentum vector along with the 2nd principal axis of the intrinsic shear tensor since $`I_{23}^2=I_{13}^2=I_{12}^2`$. What has been found in LP00 numerical simulations is, however, far from being idealistic. The principal axis of the inertia tensor has turned out to be quite strongly correlated with that of the shear tensor. However, a slight but detectable misalignment between the two tensors has been detected by LP00. It means that in spite of the strong correlation between the inertia and shear tensors a net nonzero angular momentum at first order is indeed generated to a detectable degree with its axis preferentially aligned with the intermediate principal axis of the local shear.

The essence of this idea is well represented by the following simple equation (see Appendix A):

$$\widehat{L}_i\widehat{L}_j|\widehat{𝐓}=\frac{1+a}{3}\delta _{ij}a\widehat{T}_{ik}\widehat{T}_{kj}.$$
(3)
Here $`\widehat{𝐓}`$ is a unit traceless local shear tensor ($`\widehat{𝐓}\stackrel{~}{𝐓}/|\stackrel{~}{𝐓}|`$ where $`\stackrel{~}{T}_{ij}T_{ij}\delta _{ij}\mathrm{Tr}/3`$), $`\widehat{𝐋}`$ is a unit galaxy spin vector, and $`a`$ is a spin-shear correlation parameter introduced by LP00 to measure the strength of the correlation between the local shear and the galaxy spin axis. If $`a=0`$, $`\widehat{L}_i\widehat{L}_j|\widehat{𝐓}=\delta _{ij}/3`$, spins are randomly oriented without any correlation with the local shears. While if the inertia and shear tensors are mutually uncorrelated, and there is no nonlinear effects, then the value of $`a`$ is calculated to be $`3/5`$ (it was mistakenly cited as unity in LP00, see Appendix A). The real value of $`a`$ should be determined empirically by numerical simulations, since in the linear theory one cannot estimate the strength of the correlation between the inertia and shear tensors from 1st principles, and one expects non-linear effects to be important as well. LP00 suggested the formula for the estimation of $`a`$ from N-body simulations:

$$a=26\widehat{\lambda }_i^2\widehat{L}_i^2.$$
(4)
Here $`\{\widehat{\lambda }_i\}_{i=1}^3`$ are the three eigenvalues of the trace-free unit shear tensor, satisfying $`_i\widehat{\lambda }_i^2=1`$ and $`_i\widehat{\lambda }_i=0`$, while $`\widehat{𝐋}`$ is the unit angular momentum vector measured in the shear principal axis frame. Note that if equation (3) holds as a theoretical estimation for $`\widehat{L}_i\widehat{L}_j`$ , then equation (4) becomes optimal (see Appendix J). LP00 found $`a0.24`$ in their N-body simulations. For the detailed description of the measurement of $`a`$ from N-body simulations used by LP00, see $`\mathrm{\S }2`$ in LP00.

The numerical result of LP00 indicates that the present galaxy spin axes are indeed (weakly but detectably) correlated with the intrinsic shears even though the correlation is not very strong. It is worth noting that $`a`$ is a universal value, independent of scale. The spin-shear correlation parameter of $`a`$, by its definition, must be measured from the tidal shears smoothed on the same scale that $`\widehat{𝐋}`$ is defined on.

Given the detectable preferential alignment of the galaxy spin along the 2nd principal axis of the intrinsic shear, it is possible to reconstruct the shear field from the observable unit galaxy spins. Let us say that we have $`m`$ galaxies with measured unit spins, $`\widehat{𝐋}(𝐱_\gamma )`$ for $`\gamma =1,2,\mathrm{},m`$. Now, we would like to find the maximum likelihood value of the traceless shear tensor, $`\stackrel{~}{𝐓}`$ at each galaxy position. Using Bayes’ theorem, $`P(\stackrel{~}{𝐓}|\widehat{𝐋})=P(\widehat{𝐋}|\stackrel{~}{𝐓})P(\stackrel{~}{𝐓})/P(\widehat{𝐋})`$. An immediate complication arises: $`P(\widehat{𝐋}|\stackrel{~}{𝐓})`$ is a purely local process, independent of the events at any other point, while $`P(\stackrel{~}{𝐓})=P[\stackrel{~}{𝐓}(𝐱_1),\stackrel{~}{𝐓}(𝐱_2),\mathrm{}]`$ is a joint random process linking different points with one another. As we noted in LP00, the linear shear expectation $`\stackrel{~}{𝐓}|\widehat{𝐋}=0`$ since $`P(\widehat{𝐋}|\stackrel{~}{𝐓})`$ is an even function of $`\stackrel{~}{𝐓}`$. Since we are using only directions of the spins (which are more readily observed and predicted), we cannot recover the magnitude of the shear field. In other words, the shear field can be reconstructed up to the ambiguity of a multiplicative normalization constant. Thus we can arbitrarily normalize the shear field. Here we use the normalization constraint of $`\stackrel{~}{T}_{ij}(𝐱)\stackrel{~}{T}_{ij}(𝐱)𝑑𝐱=1`$. The nontrivial quadratic maximum likelihood value of the shear field with this constraint is given as the solution to the following eigenvector equation (Appendix D):

$$\stackrel{~}{\xi }_{ijlm}(𝐱_\alpha ,𝐱_\beta )\stackrel{~}{T}_{ij}(𝐱_\alpha )d^3𝐱_\alpha =\mathrm{\Lambda }\stackrel{~}{T}_{lm}(𝐱_\beta ).$$
(5)
$`\mathrm{\Lambda }`$ is the largest eigenvalue of the posterior correlation operator,
$`\stackrel{~}{\xi }_{ijlm}(𝐱_\alpha ,𝐱_\beta )\stackrel{~}{T}_{ij}(𝐱_\alpha )\stackrel{~}{T}_{lm}(𝐱_\beta )|\widehat{𝐋}`$. In the asymptotic case of $`a1`$ as in LP00 simulation results, $`\stackrel{~}{\xi }_{ijlm}(𝐱_\alpha ,𝐱_\beta )`$ is given (Appendix C) as

$$\stackrel{~}{\xi }_{ijlm}(𝐱_\alpha ,𝐱_\beta )=a\stackrel{~}{C}_{ijnk}(𝐱_\alpha 𝐱_\gamma )\stackrel{~}{C}_{lmok}(𝐱_\beta 𝐱_\gamma )\widehat{L}_n(𝐱_\gamma )\widehat{L}_o(𝐱_\gamma )d^3𝐱_\gamma ,$$
(6)
where $`\stackrel{~}{𝐂}`$ is a two-point covariance matrix of the traceless shear tensor defined as $`\stackrel{~}{𝐂}=(\stackrel{~}{C}_{ijkl})=\stackrel{~}{T}_{ij}(𝐱)\stackrel{~}{T}_{kl}(𝐱+𝐫)`$ (see Appendix B).

In Appendix D, we explain in detail using the Lagrange multiplier method that the eigenvector associated with the largest eigenvalue is indeed the maximum likelihood expectation value of the shear. It is worth mentioning that in practice the posterior correlation function is defined only accurately at each galaxy position, so the integral in equation (6) must be replaced by a sum over discrete galaxy positions. We have regarded small $`a`$ ($`a1`$) as the limit of small signal to noise. We can also find a general expression for the shear reconstruction in Fourier space (Appendix E):

$$\frac{\stackrel{~}{\xi }_{ijlm}(𝐤_\alpha ,𝐤_\beta )}{P(k_\alpha )P(k_\beta )}\stackrel{~}{T}_{ij}(𝐤_\alpha )d^3𝐤_\alpha =\mathrm{\Lambda }\stackrel{~}{T}_{lm}(𝐤_\beta ).$$
(7)
Here $`P(k)`$ is the density power spectrum. Note that equation (7) is the optimal-filtered version of equation (5), holding without the constraint of small $`a`$.

Now, the expected shear field given the unit spin field can be found as the eigenvector of $`\stackrel{~}{\xi }_{ijlm}`$ associated with the largest eigenvalue. LP00 suggested an effective power iteration scheme to estimate the largest eigenvector of $`\stackrel{~}{\xi }_{ijlm}`$: One starts with an initial guess $`\stackrel{~}{T}_{ij}^0`$, and defines an iteration such that

$`\stackrel{~}{T}_{ij}^{n+1/2}`$ $`=`$ $`{\displaystyle \stackrel{~}{\xi }_{ijlm}\stackrel{~}{T}_{lm}^nd^3𝐤},`$
$`\stackrel{~}{T}_{ij}^{n+1}`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{T}_{ij}^{n+1/2}}{\sqrt{(\stackrel{~}{T}_{ij}^{n+1/2})^2d^3𝐤}}}+\stackrel{~}{T}_{ij}^{n1}.`$ (8)
Sufficiently large number of iterations converges the testing vector to the solution, i.e., the eigenvector associated with the largest eigenvalue with a small fractional estimation error proportional to $`(\mathrm{\Lambda }_1/\mathrm{\Lambda }_0)^m`$, where $`m`$ is the number of iterations, $`\mathrm{\Lambda }_0`$ and $`\mathrm{\Lambda }_1`$ are the largest and second largest eigenvalues respectively. Appendix F gives a general proof for the power iteration.

In order to find the expected shear field by solving equation (7) using the above iteration method, one has to know the power spectrum of the mass density field beforehand. Here we describe how one can actually determine the slope of the linear power spectrum in deriving the shear field by equation (7): From an observed set of $`N`$ galaxies and with an initial guess for the mass power spectrum, we construct a posterior shear correlation function (6), which is a $`5N\times 5N`$ matrix. From this posterior shear correlation function, one can construct a weighed posterior shear correlation function given in equation (7), whose largest eigenvector is the shear field to be reconstructed by the iteration method described above. Here the largest eigenvalue $`\mathrm{\Lambda }`$ is the likelihood. We can iterate this procedure itself to measure a self-consistent power spectrum by varying the power spectrum to maximize $`\mathrm{\Lambda }`$. In other words, at the same time when one reconstructs the initial shear field, one can also measure the slope of the initial power spectrum by finding such power spectrum as maximize $`\mathrm{\Lambda }`$ in equation (7). Note, however, that one can recover the slope but not the amplitude of the power spectrum since the likelihood $`\mathrm{\Lambda }`$ in (7) is independent of a multiplicative constant of the power spectrum.

The final step is the reconstruction of the density field, $`\delta (𝐱)`$ given the traceless shears, $`\stackrel{~}{𝐓}(𝐱_i)`$ reconstructed at each galaxy position $`𝐱_i`$. First we consider a orthonormal parametrization of the density and the five free components of the traceless shear tensor such that

$`y_0`$ $`=`$ $`{\displaystyle \frac{\delta }{\sqrt{3}}},`$
$`y_1`$ $`=`$ $`{\displaystyle \frac{(3\sqrt{3})\stackrel{~}{T}_{11}+2\sqrt{3}\stackrel{~}{T}_{22}+(3\sqrt{3})\stackrel{~}{T}_{33}}{6}},`$
$`y_2`$ $`=`$ $`{\displaystyle \frac{(3\sqrt{3})\stackrel{~}{T}_{11}+2\sqrt{3}\stackrel{~}{T}_{22}+(3\sqrt{3})\stackrel{~}{T}_{33}}{6}},`$
$`y_3`$ $`=`$ $`\sqrt{2}\stackrel{~}{T}_{12},y_4=\sqrt{2}\stackrel{~}{T}_{23},y_5=\sqrt{2}\stackrel{~}{T}_{31}.`$ (9)
$`𝐲`$ in fact is a orthonormal vector-representation of the full shear, $`𝐓`$ in terms of trace and traceless parts. Therefore the mutual correlation between the six components of $`𝐲`$ at the same position is always zero.

Reconstructing $`\delta (𝐱)`$ given $`\widehat{𝐓}(𝐱_i)`$ amounts to finding $`y_0(𝐱)|y_1(𝐱_i),y_2(𝐱_i),\mathrm{},y_5(𝐱_i)`$. Since a linear combination of the Gaussian variables is also Gaussian, $`𝐲`$ is a Gaussian variable, and the covariance matrix of $`𝐲`$, say, $`𝐕`$ ($`V_{ij}y_iy_j`$) can be obtained by the linear transformation of the shear two-point correlations (eq. \[B3\]). We obtain the following expression for $`y_0(𝐱)|y_1(𝐱_i),y_2(𝐱_i),\mathrm{},y_5(𝐱_i)`$ (see Appendix G):

$$y_0(𝐱)|y_1(𝐱_i),y_2(𝐱_i),\mathrm{},y_5(𝐱_i)=\frac{U_{0\nu }y_\nu }{U_{00}},$$
(10)
where $`𝐔𝐕^1`$, and the Greek index, $`\nu `$ goes from $`1`$ to $`5`$. Equation (10) allows us to reconstruct the density field at an arbitrary spatial position $`𝐱`$ once the traceless shear field is reconstructed at each galaxy position $`𝐱_i`$. Note that the only mathematical complication that arises in the reconstruction algorithm is a matrix inversion. Therefore it is computationally tractable, involving only linear algebra. It is worth mentioning that although the density field is supposed to be reconstructed in Lagrangian space, the galaxy spins are measured in Eulerian redshift space. We can regard this displacement between the Eulerian and Lagrangian spaces as noise, and convolve simply the two-point density correlation function, $`\xi (r)`$ with a Gaussian filter with a peculiar velocity dispersion $`\sigma _v=150`$ km/s for spiral galaxies (see Davis, Miller, & White 1997).

Fig.1 shows the accuracy of the reconstructed density field. It plots the cross correlation coefficient $`r=\delta _r(k)\delta _m(k)/\sqrt{P_r(k)P_m(k)}`$ between the reconstructed density field $`\delta _r`$ and the true density field $`\delta _m`$. We have implemented the algorithm with the realistic value of the correlation parameter of $`a=0.24`$ and simulated two million sample galaxy spins, and find a good ability of the algorithm to reconstruct the density field. As was shown in PLS00 and Crittenden et al. (2000a, 2000b), the intrinsic galaxy alignment is a much stronger effect than weak lensing for shallow surveys such as SDSS. Given that it was hoped that weak lensing power spectra could be measured, and thus the projected mass power spectrum, we would clearly expect that this stronger effect of intrinsic alignments would thus allow a better reconstruction of the density field in the source plane. The advantage now is that the full three dimensional shear field can be reconstructed, not just a two dimensional projected field as would be the case for the weak lensing. The above algorithm described in this section enables us to achieve this goal.

## 3 SPIN-DENSITY CROSS CORRELATION

In order for the algorithm given in $`\mathrm{\S }2`$ to be applied to the real universe, it is indispensable to have a nonzero shear-spin correlation parameter, $`a`$. Unfortunately, however, it is quite hard to measure the value of $`a`$ directly from real observational data since it requires us to know the initial shear field beforehand (see eq. ). An alternative simpler way to detect the intrinsic spin alignment is to investigate the spatial spin-spin correlation. The local shear effect is due to the surrounding matter distribution. But the matter in the universe is spatially correlated in the standard structure formation scenario. Consequently, if galaxy spins are indeed correlated with the linear shear tensor by equation (3) with a nonzero value of $`a`$, the spatial shear correlation must induce a spatial spin-spin correlation with themselves.

Since galaxy formation is an unsolved problem, it is difficult to make an accurate quantitative evaluation of the expected level of the spin-spin correlation signal. We can make at most approximate analytic estimates for the order of magnitude of the spin correlation and its qualitative behavior. PLS00 have attempted to estimate the expected strength of the galaxy spin-spin correlation (see eq.  in PLS00) using the 1st-order perturbation theory and the numerical normalization amplitude of $`a=0.24`$. Appendix H lays out the detailed derivation of the spin-spin correlation function presented in PLS00. PLS00 then measured the spin correlation signal directly from the observed spiral galaxies of the Tully catalog. The observed signal turned out to be significant at $`97\%`$ confidence level, and the amplitude of the signal is of order of $`1\%`$ at a separation of $`1h^1`$ Mpc, in agreement with the PLS00 theoretical estimates.

The consistent results of the observed spin-spin correlation signal with the theoretical predictions motivate us to consider a correlation of galaxy spin field with the density field. Galaxy spin alignment with the local gravitational shear field might result in the correlation of galaxy spins with the directional geometry of the nearby galaxy distribution. Perhaps the simplest statistic to observe is the correlation between the spin axis $`\widehat{𝐋}`$ and the spatial direction $`\widehat{𝐫}`$ to the nearest neighbor. Therefore, we first define a simple nontrivial spin-direction cross correlation function analogous to the spin-spin correlation function given by PLS00 such that :

$$\omega (r)|\widehat{𝐋}(𝐱)\widehat{𝐫}(𝐱)|^2\omega _0,$$
(11)
where $`\omega _0`$ is the value of $`\omega (r)`$ for the case of no correlation: $`\omega _0=1/3`$ for the three dimensional case while $`\omega _0=1/2`$ for the two dimensional case. Note that the two vectors, $`\widehat{𝐋}`$ and $`\widehat{𝐫}`$, are both defined at the same galaxy position, $`𝐱`$.

For a galaxy pair at $`𝐱`$ and $`𝐱+𝐫`$, let us consider the density and shear fields smoothed on two different top-hat scales, say $`R`$ and $`R^{}`$, where $`R`$ is the top-hat galactic radius, while $`R^{}`$ is the minimum top-hat radius that encloses the galaxy pair such that $`R^{}=R+r`$. In order to avoid confusion about the smoothing scale, in this section we use an explicit notation of $`\delta _R`$ and $`\widehat{𝐓}^R`$ respectively to represent the density and unit traceless shear fields smoothed on a scale of $`R`$, while $`\delta _R^{}`$ and $`\widehat{𝐓}^R^{}`$ for the density and unit traceless shear fields smoothed on a scale of $`R^{}`$. A simple directional vector that one could form is the gradient of the smoothed density field, $`\delta _R^{}(𝐱)`$. Assuming that galaxies form on peaks of the density field, one expects two neighboring galaxies to sit at the ends of a ridge connecting the two galaxies. A peak is by definition a location where the gradient is zero. If one considers the gradient halfway between the two peaks, one expects it to generically be near a saddle point, where again the gradient is zero. And in between, the gradient would be expected to point in a direction perpendicular to the peak separation $`𝐫`$.

Instead, one would expect that the direction of the galaxy separation vector correlates with the major principal axis of the local gravitational shear tensor smoothed on a scale of the galaxy separation. If we neglect the other two principal axes (as is often the case in principal component analysis), we can relate $`\widehat{r}_i\widehat{r}_j`$ to $`\widehat{𝐓}`$ up to a considerable ambiguity. But, it results in a trivial spin-direction correlation: $`\omega (r)\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{ij}=0`$ due to the even probability distribution of $`\widehat{𝐓}`$. We note, however, that the principal axis of a shear tensor is the same as that of its square, so we shall instead relate $`\widehat{r}_i\widehat{r}_j`$ to $`\widehat{T}_{ik}^R^{}\widehat{T}_{kj}^R^{}`$ analogous to equation (3) such that

$$\widehat{r}_i\widehat{r}_j|\widehat{𝐓}^R^{}=\frac{1b}{3}\delta _{ij}+b\widehat{T}_{ik}^R^{}\widehat{T}_{kj}^R^{},$$
(12)
where we introduce a new quantity, a direction-shear correlation parameter of $`b`$, to measure the strength of the correlation between the unit galaxy separation and the major axis of the local shear tensor. A careful reader may have noticed the difference of the sign ahead of the correlation parameters between equations (3) and (12). This sign difference arises because the direction of each alignment with the unit shear tensor is different. $`\widehat{𝐫}`$ is aligned with the major principal axis of $`(\widehat{T}_{ik}\widehat{T}_{kj})`$ while $`\widehat{𝐋}`$ is aligned with the minor principal axis. Let $`\lambda _1,\lambda _2,\lambda _3`$ be the three eigenvalues of the shear tensor, $`𝐓`$ with the order of $`\lambda _1>\lambda _2>\lambda _3`$. Then the three eigenvalues of the unit traceless shear tensor $`\widehat{𝐓}`$ is nothing but $`\widehat{\lambda }_1,\widehat{\lambda }_2,\widehat{\lambda }_3`$ with $`\widehat{\lambda }_i=\lambda _i\mathrm{Tr}/3`$, $`\mathrm{Tr}=_{i=1}^3\lambda _i`$. Obviously the order is the same: $`\widehat{\lambda }_1>\widehat{\lambda }_2>\widehat{\lambda }_3`$. Thus, the principal axes of $`𝐓`$ and $`\widehat{𝐓}`$ coincide. But if we consider the square of the shear, $`(\widehat{T}_{ik}\widehat{T}_{kj})`$, the eigenvalues are given as $`\widehat{\lambda }_1^2,\widehat{\lambda }_2^2,\widehat{\lambda }_3^2`$ but with the order of $`\widehat{\lambda }_1^2\widehat{\lambda }_3^2>\widehat{\lambda }_2^2`$. Thus, the intermediate axis of the shear tensor becomes the minor axis of $`(\widehat{T}_{ik}\widehat{T}_{kj})`$. It explains why $`\widehat{r}_i\widehat{r}_j|\widehat{𝐓}`$ is positively proportional to $`\widehat{T}_{ik}\widehat{T}_{kj}`$ (apart from the shear-independent constant) while $`\widehat{L}_i\widehat{L}_j|\widehat{𝐓}`$ is negatively proportional to $`\widehat{T}_{ik}\widehat{T}_{kj}`$. Also note that in equation (3) $`\widehat{𝐓}`$ is smoothed on the top-hat galactic scale of $`R`$ while in equation (12) $`\widehat{𝐓}`$ is smoothed on the minimum enclosing top-hat radius of $`R^{}`$ since the galaxy separation vector can be defined for a galaxy pair not for one galaxy.

This direction-shear correlation parameter of $`b`$ can be also determined in N-body simulations in principle. We suggest the following formula for estimation of $`b`$ in simulations:

$$b=\sqrt{2}\widehat{\lambda }_i\widehat{r}_i^2,$$
(13)
where $`\widehat{𝐫}`$ is the unit separation vector in the shear principal axis frame. In practice, each $`\widehat{r}_i`$ is obtained by measuring the separation vector of each closest galaxy pair and projecting the separation vector into the $`i`$-th principal axis of the local shears smoothed on the mean galaxy separation. Again equation (13) becomes optimal if equation (12) holds as a theoretical estimation formula for $`\widehat{r}_i\widehat{r}_j`$ (see Appendix J). We have found the average value of $`b=0.29\pm 0.01`$ from the same N-body simulation results that LP00 used for the measurement of $`a`$. It is worth mentioning, however, that the galaxy distribution is known to have a correlation function significantly different from that of the matter, measuring the value of $`b`$ requires a quantitative galaxy formation model beforehand. Thus, with having no quantitative galaxy formation model, equation (13) provides only a qualitative approximation for the magnitude of $`b`$.

With the similar method that we have used for the spin-spin correlation, one can find an analytic estimates of $`\omega (r)`$ (Appendix I) such that:

$$\omega (r)=A\frac{\delta _R^{}\delta _R^2}{(\sigma _R^{}\sigma _R)^2},$$
(14)
where the amplitude of $`A`$ depends on the correlation parameters, $`a`$ and $`b`$. It has the value of $`ab/6`$ and $`5ab/24`$ for the three and the two dimensional cases respectively. Here $`\delta _R\delta _R^{}`$ is the auto correlation of the density field smoothed on two different scales of $`R`$ and $`R^{}`$, and $`\sigma _R`$ and $`\sigma _R^{}`$ are the corresponding rms density fluctuations:

$`\delta _R(𝐱)\delta _R^{}(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^3}}{\displaystyle W_{th}(kR)W_{th}(kR^{})P(k)4\pi k^2𝑑k},`$
$`\sigma _R^2=\delta _R^2(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^3}}{\displaystyle W_{th}^2(kR)P(k)4\pi k^2𝑑k},`$
$`\sigma _R^{}^2=\delta _R^{}^2(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^3}}{\displaystyle W_{th}^2(kR^{})P(k)4\pi k^2𝑑k},`$ (15)
where the top-hat window function is given as $`W_{th}(kR)=3[\mathrm{sin}(kR)kR\mathrm{cos}(kR)]/(kR)^3`$.

To find a closed analytic form of $`\omega (r)`$, we can replace the top-hat filter with the Gaussian filter. Using a Gaussian filter of $`W_G(kR)=\mathrm{exp}(k^2R^2/2)`$, and a power-law power spectrum of $`P(k)=k^2`$, we find

$$\omega (r)=A\frac{2R^{}R}{R^2+R_{}^{}{}_{}{}^{2}}.$$
(16)
Equation (16) says that for neighboring galaxies, $`|\omega (r)|`$ decreases as $`r^1`$, less rapidly than the spin-spin correlation that decreases as $`r^2`$ (see PLS00). Note that in Lagrangian space the galaxy separation of $`r`$ cannot decrease below $`3R`$ since the top-hat radius enclosing a galaxy pair must be at least three galactic scale radius of $`R`$. Thus we assign $`\omega (r)`$ a constant value of $`\omega (3R)`$ for $`r3R`$.

## 4 SIGNAL FROM THE REAL UNIVERSE

The unique and advantageous feature of the galaxy spin statistics presented in $`\mathrm{\S }2`$ and $`\mathrm{\S }3`$ is that it is a readily testable theory against real observational data since it deals not with the magnitude but with the direction of a galaxy spin. The spin axis of a spiral galaxy can be easily determined by the information of the position angle (PA) and axial ratio ($``$): A spiral galaxy is a thin disk with a circular face-on shape, and its spin vector is perpendicular to the plane of the disk. Therefore, the apparent axial ratio gives the magnitude of the radial component of a spin vector, while the position angle determines the relative magnitude of the tangential components of the spin vector lying in the plane of the sky.

In order to apply observational tests to our theory, the most suitable dataset should be a large sample of spiral galaxies at low redshift with the information of position angle and axial ratio ($``$). The low redshift condition is required since at high redshift the weak lensing shear effect on the apparent orientation of the spin axis is dominant (Jain & Seljak, 1997; Wittman et al., 2000; Heavens et al., 2000). B. Tully (2000, private communication) generously has provided such a galaxy catalog: The Tully galaxy catalog is a compilation of 35674 nearby galaxy properties over the whole sky with median redshift of $`6740`$ km/s. Among the total 35674 Tully galaxy properties, 12122 galaxies are identified as spirals.

In measuring the spin-direction correlation, we consider all the 35674 galaxies in the Tully catalog to calculate the direction vectors, while we used only the spiral galaxies to measure the spin vectors. As mentioned in PLS00, we suspect that the shape-shape correlation of galaxies might cause a potential problem as a form of $``$-related systematic errors. The $``$-related systematic errors are involved in the measurement of the axial ratio, $``$, found in the Tully catalog, caused presumably by the deviation of the shape of spiral galaxies from a perfect ellipse, finite thickness of galaxies, and etc. For the detailed description of the Tully catalog and the data analysis, see $`\mathrm{\S }3`$ of PLS00.

An easy way to avoid any false signal from the $``$-related systematic errors is to measure the two dimensional spin-direction correlation. We project the three dimensional unit spin vector, $`\widehat{𝐋}`$ and unit separation vector, $`\widehat{𝐫}`$ onto the plane of the sky to obtain the two dimensional unit spin vector, $`\widehat{𝐒}=𝐒/|𝐒|`$, $`𝐒=\widehat{𝐋}(\widehat{𝐋}\widehat{𝐱})\widehat{𝐱}`$, and two dimensional separation vector, $`\widehat{𝐭}=𝐭/|𝐭|`$, $`𝐭=𝐫(\widehat{𝐫}\widehat{𝐱})\widehat{𝐱}`$. Now, the two dimensional spin-direction correlation is given by $`\omega _{2D}(r)=|\widehat{𝐒}\widehat{𝐭}|^21/2`$. Note that the projection of the spin vector onto the plane of the sky amounts to setting $`=0`$, so $`\omega _{2D}(r)`$ is free of the $``$-related systematic errors.

For the three dimensional spin-direction correlation, we use an effective redistribution-method to deal with the $``$-related systematic errors. We first bin the separation of every galaxy pair. At each bin we uniformly redistribute the given $``$ of each galaxy spin in the range of $`[0,1)`$ and the radial component of the given $`\widehat{𝐫}`$ in range of $`(1,1)`$ (the radial direction is along the line-of-sight at each galaxy position). We expect that the uniform redistribution of $``$ and $`\widehat{𝐫}`$ of galaxy pairs belonging to each bin eliminates effectively the systematic bias and false signal. Now we renormalize the spin and the separation vectors after the uniform redistribution of $``$ and $`\widehat{𝐫}`$ at each bin, and calculate the three dimensional spin-direction correlation, $`\omega _{3D}(r)=|\widehat{𝐋}\widehat{𝐫}|^21/3`$.

Fig. 2 plots the resulting observed signal (filled squares) vs. the galaxy separation, $`r=cz`$ (km/s) with error bars, and compares the observed signals with the theoretical estimates. Regarding the theoretical curves in Fig. 2, we convolve the Lagrangian correlation (dashed line) by a Gaussian filter with $`\sigma _v=150`$ km/s to obtain the Eulerian correlation (solid line) the observed signal is measured in Eulerian redshift space. For a detailed description of the convolution procedure, see also $`\mathrm{\S }3`$ in PLS00. The error bars are obtained from the experiment with the 500 sets of 12122 random two dimensional unit spins. We first generate the 12122 random spin vectors, and calculate the spin-direction correlation. We repeat this process 1000 times with different sets of random spin vectors, and compute the standard deviation of the spin-direction correlations. The solid line is the theoretical predictions given by equation (14) with the normalization amplitudes of $`a=0.24`$ and $`b=0.3`$ for the case of a power-law spectrum of $`P(k)=k^2`$.

Although the observed spin-direction correlation is fairly consistent with the theoretical estimates qualitatively, the signal is quite weak, and the error bars are still large. We expect that larger surveys like SDSS will make a precision measurement of the spin-density correlation signal in the near future.

## 5 SUMMARY AND CONCLUSIONS

We have presented the technical formalism in which we discuss the intrinsic galaxy spin correlation. We have shown how the intrinsic spin correlation is related to the initial potential and density fields, and how the problem can be inverted to derive the power spectrum and density field up to a multiplicative constant from the observable orientation of galaxy spins, as originally claimed by LP00. Since the intrinsic galaxy alignments are expected to dominate the weak lensing signal for shallow surveys such as SDSS, our algorithm for the density reconstruction by the intrinsic galaxy spin alignment should be more viable than the one by the weak lensing shear effect.

The formalism also allows us to address the issue of the spin-direction correlation, which we have estimated theoretically, and measured in the observational catalog. Although the observed signal is reasonably consistent with the theoretical estimates, the signal is quite weak, and the error bars are still large. We encourage future works on the spin-direction correlation with larger surveys, which will make a precision measurement of the spin-direction correlation signal.

We are very grateful to B. Tully for his catalog. We also thank U. Seljak for useful discussions. This work has been supported by Academia Sinica and partially by NSERC grant 72013704 and computational resources of the National Center for Supercomputing Applications.

## Appendix A SHEAR-SPIN CORRELATION

In order to find the expectation value of the unit galaxy angular momentum product given the unit shear tensor, $`\widehat{L}_i\widehat{L}_j|\widehat{𝐓}`$, let us first consider $`L_iL_j|𝐓`$. If the local shear tensor, $`𝐓`$ and the inertia tensor, $`𝐈`$ are mutually uncorrelated, the ensemble average of equation (1) over all orientations of the inertia tensor gives

$$L_iL_j|𝐓=ϵ_{iab}ϵ_{jcd}T_{ak}T_{cl}I_{kb}I_{ld}.$$
(A1)
where the time-dependent proportionality constant in equation (1) is set to be unity at present epoch.

From the statistical isotropy of the inertia tensor, we have the inertia tensor correlation $`I_{kb}I_{ld}=\gamma (\delta _{kb}\delta _{ld}+\delta _{kl}\delta _{bd}+\delta _{kd}\delta _{bl})/3`$ where $`\gamma (I_{11}^2+I_{22}^2+I_{33}^2)/3`$ is the proportionality constant. Here we stress that the density reconstruction algorithm depends not on the magnitude but on the direction of $`𝐋`$. Thus, the overall proportionality constants that arise in the middle of our derivations can be always set to unity and we set $`\gamma =1`$ hereafter. Of course, the renormalized angular momentum does not have the same magnitude as the original angular momentum. Furthermore, it does not have the dimension of the angular momentum anymore. But this renormalized vector does have the same direction as the original angular momentum, which is all that matters. Hereafter, we will use this kind of renormalization frequently by setting any proportionality constant to be unity whenever it does not affect the direction of the angular momentum regardless of the dimensionality.

Using $`I_{kb}I_{ld}=(\delta _{kb}\delta _{ld}+\delta _{kl}\delta _{bd}+\delta _{kd}\delta _{bl})/3`$, equation (A1) can be rewritten as

$$L_iL_j|𝐓=\frac{ϵ_{iab}ϵ_{jcd}\left(T_{ab}T_{cd}+T_{ad}T_{bd}\right)+\delta _{ij}|T|^2T_{ik}T_{kj}}{3}.$$
(A2)
On can verify that (A2) does not depend on the trace of the shear, so we rewrite it in terms of a traceless shear tensor $`\stackrel{~}{T}_{ij}=T_{ij}\delta _{ij}\mathrm{Tr}/3`$ such that,

$$L_iL_j|\stackrel{~}{𝐓}=\frac{2}{3}\delta _{ij}|\stackrel{~}{T}|^2\stackrel{~}{T}_{ik}\stackrel{~}{T}_{kj}.$$
(A3)
Note that $`L_iL_i|\stackrel{~}{𝐓}=|\stackrel{~}{𝐓}|^2`$ for this case of independent of $`𝐓`$ and $`𝐈`$. Since we are again only interested in the direction of the angular momentum vector, we may rescale equation (A3) to have the normalization constraint of $`L_iL_i|\widehat{𝐓}=1`$, dividing each side by $`|\stackrel{~}{𝐓}|^2`$ such that

$`L_iL_j|\widehat{𝐓}`$ $`=`$ $`{\displaystyle \frac{2}{3}}\delta _{ij}\widehat{T}_{ik}\widehat{T}_{kj}`$ (A4)
$`=`$ $`{\displaystyle \frac{1}{3}}\delta _{ij}+\left({\displaystyle \frac{1}{3}}\delta _{ij}\widehat{T}_{ik}\widehat{T}_{kj}\right).`$
The first term in the RHS of equation (A4) corresponds to the stochastic sources uncorrelated with the initial shear field. That is, the first term represents the ensemble average of $`L_iL_j`$ for the case that the direction of the angular momentum is completely random, having no correlation with the shear axis, due to the modification by the other stochastic sources such as nonlinear effect, mutual correlation between $`𝐈`$ and $`𝐓`$, and etc. While the second term in the RHS of equation (A4) corresponds to the deviation of the spin direction from the random average value due to it tendency to align preferentially with the intermediate axis of the shear tensor.

This linear theory prediction for $`L_iL_i|\widehat{𝐓}`$ holds provided that there is no nonlinear effects and that $`𝐓`$ and $`𝐈`$ are mutually independent. In practice, however, this condition is not guaranteed. Here we adopt the following simple assumption: the non-linear and stochastic effects are uncorrelated with the linear prediction, which adds noise to the unit spin vector. The nonlinear effect and the mutual dependence between $`𝐓`$ and $`𝐈`$ decreases the relative weight of the second shear-dependence term in equation (A4), which can be quantified by one parameter (say, $`c`$) such that

$$L_iL_j|\widehat{𝐓}=\frac{1}{3}\delta _{ij}+c\left(\frac{1}{3}\delta _{ij}\widehat{T}_{ik}\widehat{T}_{kj}\right).$$
(A5)
As the value of $`c`$ decreases, the first term of the RHS in equation (A5) dominates more, making the direction of the spin vector more random. the linear theory predictions with perfectly independent $`𝐈`$ and $`𝐓`$ corresponds to $`c`$ of unity, while the completely random direction of the angular momentum vector corresponds to $`c=0`$. This is the generalized quadratic relation we propose to express the correlation between the direction of the present galaxy angular momentum vector and the initial gravitational shear tensor. Retaining the framework of the linear perturbation theory, we treat the existence of nonlinear effect and the correlation between the shear and the inertia tensors as stochastic sources that tend to randomize the direction of the spin vector, decreasing the correlation of the angular momentum with the initial shear tensor.

Now, to find the expression for $`\widehat{L}_i\widehat{L}_j|\widehat{𝐓}`$, let us find the conditional probability density function, $`P(\widehat{𝐋}|\stackrel{~}{𝐓})`$. The conditional probability density function, $`P(𝐋|\stackrel{~}{𝐓})`$ is usually given as a Gaussian distribution (see Catelan & Theuns 1996) such that

$$P(𝐋|\stackrel{~}{𝐓})=\frac{|𝐐|^{1/2}}{\sqrt{(2\pi )^3}}\mathrm{exp}\left(\frac{𝐋^\mathrm{T}𝐐^1𝐋}{2}\right),$$
(A6)
where the covariance matrix $`𝐐`$ is defined in equation (A5). The conditional probability density distribution $`P(\widehat{𝐋}|\widehat{𝐓})`$ can be derived by integrating out $`P(𝐋|\stackrel{~}{𝐓})`$ over the magnitude of $`L=|𝐋|`$ such that

$$P(\widehat{𝐋}|\widehat{𝐓})=P(𝐋|\widehat{𝐓})L^2𝑑L=\frac{|\widehat{𝐐}|^{1/2}}{4\pi }\left(\widehat{𝐋}^\mathrm{T}\widehat{𝐐}^1\widehat{𝐋}\right)^{3/2},$$
(A7)
where $`P(\widehat{𝐋}|\stackrel{~}{𝐓})`$ is in fact equal to $`P(\widehat{𝐋}|\widehat{𝐓})`$. Here the unit covariance matrix, $`\widehat{Q}_{ij}`$ is given by equation (A5).

In the limit of $`c1`$, equation (A7) is simplified into

$$P(\widehat{𝐋}|\widehat{𝐓})=\frac{1}{4\pi }\left(1+\frac{3c}{2}[13\widehat{T}_{ik}\widehat{T}_{kj}\widehat{L}_i\widehat{L}_j]\right).$$
(A8)
With equation (A8) and the help of little algebra, it is straightforward to calculate the expectation value of $`\widehat{L}_i\widehat{L}_j|\widehat{𝐓}`$ in the limit of $`c1`$ such that

$`\widehat{L}_i\widehat{L}_j|\widehat{𝐓}`$ $`=`$ $`{\displaystyle \widehat{L}_i\widehat{L}_jP(\widehat{𝐋}|\widehat{𝐓})𝑑\widehat{𝐋}},`$ (A9)
$`=`$ $`\left({\displaystyle \frac{1}{3}}+{\displaystyle \frac{c}{5}}\right)\delta _{ij}{\displaystyle \frac{3}{5}}c\widehat{T}_{ik}\widehat{T}_{kj}.`$
Let us define a correlation parameter, $`a3c/5`$. Then, we finally get the desired expression :

$$\widehat{L}_i\widehat{L}_j|\widehat{𝐓}=\frac{1+a}{3}\delta _{ij}a\widehat{T}_{ik}\widehat{T}_{kj}.$$
(A10)
This equation says that for the ideal case of independent $`𝐈`$ and $`𝐓`$, the correlation parameter has the value of $`a=3/5`$ (corresponding to $`c=1`$) while for the random spins having no dependence on the shear tensor $`a=0`$ ($`c=0`$).

For practical purposes, it is also useful to have a similar expression to equation (A10) for the two dimensional unit spins. The two dimensional unit spins mean the galaxy spins projected onto the plane of sky, and normalized to have a unit magnitude. Let $`(\widehat{S}_1,\widehat{S}_2)=(\mathrm{cos}\varphi ,\mathrm{sin}\varphi )`$ be the two dimensional unit spins, and $`P(\widehat{S}_1,\widehat{S}_2|\widehat{𝐓})d\widehat{S}=P(\varphi |\widehat{𝐓})d\varphi `$ be the conditional probability density distribution of the two dimensional unit spins. Using the flat-sky approximation with $`\widehat{L}_3`$ in the line-of-sight direction, one can say $`\widehat{L}_1=\sqrt{1\widehat{L}_3^2}\mathrm{cos}\varphi `$, $`\widehat{L}_2=\sqrt{1\widehat{L}_3^2}\mathrm{sin}\varphi `$. Then, one can say $`P(\widehat{L}|\widehat{𝐓})=P(\varphi ,\widehat{L}_3|\widehat{𝐓})`$. Now, $`P(\varphi |\widehat{𝐓})`$ can be obtained by integrating out $`P(\varphi ,\widehat{L}_3|\widehat{𝐓})`$ over $`\widehat{L}_3`$ such that

$`P(\varphi |\widehat{𝐓})`$ $`=`$ $`{\displaystyle _1^1}P(\varphi ,\widehat{L}_3|\widehat{𝐋})𝑑\widehat{L}_3`$ (A11)
$`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _1^1}\left(1+{\displaystyle \frac{5a}{2}}[13\widehat{T}_{ik}\widehat{T}_{kj}\widehat{L}_i\widehat{L}_j]\right)𝑑\widehat{L}_3`$
$`=`$ $`{\displaystyle \frac{1}{2\pi }}\left\{1+5a\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{f_3}{2}}+f_1\mathrm{cos}^2\varphi +f_2\mathrm{sin}^2\varphi 2g_{12}\mathrm{cos}\varphi \mathrm{sin}\varphi \right)\right\},`$
where $`f_i=\widehat{T}_{ik}\widehat{T}_{ki}`$ for $`i=1,2,3`$, and $`g_{12}=\widehat{T}_{1k}\widehat{T}_{2k}`$.

Using equation (A11), it is straightforward to calculate $`\widehat{S}_i\widehat{S}_j|\widehat{𝐓}`$ such that

$`\widehat{S}_1^2|\widehat{𝐓}`$ $`=`$ $`{\displaystyle _0^{2\pi }}\mathrm{cos}^2\varphi P(\varphi |\widehat{𝐓})𝑑\varphi ,`$ (A12)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{5a}{8}}(f_1f_2),`$
$`\widehat{S}_2^2|\widehat{𝐓}`$ $`=`$ $`{\displaystyle _0^{2\pi }}\mathrm{sin}^2\varphi P(\varphi |\widehat{𝐓})𝑑\varphi ,`$ (A13)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{5a}{8}}(f_2f_1).`$
$`\widehat{S}_1\widehat{S}_2|\widehat{𝐓}`$ $`=`$ $`{\displaystyle _0^{2\pi }}\mathrm{cos}\varphi \mathrm{sin}\varphi P(\varphi |\widehat{𝐓})𝑑\varphi ,`$ (A14)
$`=`$ $`{\displaystyle \frac{5a}{4}}g_{12}.`$
## Appendix B SHEAR CORRELATIONS

Let us calculate the spatial shear correlation, $`C_{ijkl}=T_{ij}(𝐱)T_{kl}(𝐱+𝐫)`$. Using $`T_{ij}=_i_j\varphi `$ and $`\varphi =^2\delta `$, one can write

$`C_{ijkl}(𝐫)`$ $`=`$ $`T_{ij}(𝐱)T_{kl}(𝐱+𝐫)`$ (B1)
$`=`$ $`_{x_i}_{x_j}_𝐱^2\delta (𝐱)_{r_k}_{r_l}_𝐫^2\delta (𝐱+𝐫).`$
Replacing the ensemble average with the spatial average by the ergodic theorem, and applying the integration by parts, one can show that equation (B1) can be rewritten as

$`C_{ijkl}(𝐫)=_i_j_k_l_𝐫^4\xi (r)`$ (B2)
Using the identity relation, $`_𝐫^2=^r𝑑r^{}(1/r^2)^r^{}𝑑r^{\prime \prime }r^{\prime \prime 2}`$, and with the help of little algebra, equation (B2) can be arranged such that

$`C_{ijkl}(𝐫)`$ $`=`$ $`(\delta _{ij}\delta _{kl}+\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk})\left\{{\displaystyle \frac{J_3}{6}}{\displaystyle \frac{J_5}{10}}\right\}+(\widehat{r}_i\widehat{r}_j\widehat{r}_k\widehat{r}_l)\left\{\xi (r)+{\displaystyle \frac{5J_3}{2}}{\displaystyle \frac{7J_5}{2}}\right\}`$ (B3)
$`+(\delta _{ij}\widehat{r}_k\widehat{r}_l+\delta _{ik}\widehat{r}_j\widehat{r}_l+\delta _{il}\widehat{r}_k\widehat{r}_j+\delta _{jk}\widehat{r}_i\widehat{r}_l+\delta _{jl}\widehat{r}_i\widehat{r}_k+\delta _{kl}\widehat{r}_i\widehat{r}_j)\left\{{\displaystyle \frac{J_5}{2}}{\displaystyle \frac{J_3}{2}}\right\},`$
where $`\widehat{𝐫}=𝐫/r`$, $`J_nnr^n_0^r\xi (r^{})r^{n1}𝑑r^{n1}`$. The two-point covariance matrix of the traceless shears $`\stackrel{~}{𝐂}`$ can be also obtained by $`\stackrel{~}{C}_{ijkl}=C_{ijkl}\delta _{kl}C_{ijnn}/3\delta _{ij}C_{mmkl}/3+\delta _{ij}\delta _{kl}C_{mmnn}/9`$.

## Appendix C POSTERIOR CORRELATION FUNCTION

Let us consider the traceless posterior correlation function,
$`\stackrel{~}{\xi }_{ijlk}(𝐱_\alpha ,𝐱_\beta )\stackrel{~}{T}_{ij}(𝐱_\alpha )\stackrel{~}{T}_{lm}(𝐱_\beta )|\widehat{𝐋}(𝐱_\gamma )=\stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta |\widehat{𝐋}^\gamma `$ in the asymptotic limit of $`a1`$. Here $`𝐱_\alpha `$ and $`𝐱_\beta `$ are two fixed galaxy positions where we would like to reconstruct the shear field, while $`𝐱_\gamma `$ represents any arbitrary position of the given $`N`$ galaxies such that the index $`\gamma =1,2,\mathrm{},N`$ is dummy. Thus, it is in fact the expectation value of the quadratic shears given the whole galaxy spin field:

$`\stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta |\widehat{𝐋}^\gamma `$ $`=`$ $`{\displaystyle 𝑑\stackrel{~}{𝐓}^\alpha 𝑑\stackrel{~}{𝐓}^\beta \stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta P(\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta |\widehat{𝐋}^\gamma )},`$ (C1)
$`=`$ $`{\displaystyle 𝒟\stackrel{~}{𝐓}^\gamma 𝑑\stackrel{~}{𝐓}^\alpha 𝑑\stackrel{~}{𝐓}^\beta \stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta P(\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta ,\stackrel{~}{𝐓}^\gamma |\widehat{𝐋}^\gamma )},`$
$`=`$ $`{\displaystyle 𝒟\stackrel{~}{𝐓}^\gamma 𝑑\stackrel{~}{𝐓}^\alpha 𝑑\stackrel{~}{𝐓}^\beta \stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta P(\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta ,\stackrel{~}{𝐓}^\gamma )\frac{P(\widehat{𝐋}^\gamma |\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta ,\stackrel{~}{𝐓}^\gamma )}{P(\widehat{𝐋}^\gamma )}},`$
$`=`$ $`{\displaystyle 𝒟\stackrel{~}{𝐓}^\gamma 𝑑\stackrel{~}{𝐓}^\alpha 𝑑\stackrel{~}{𝐓}^\beta \stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta P(\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta ,\stackrel{~}{𝐓}^\gamma )P(\widehat{𝐋}^\gamma |\stackrel{~}{𝐓}^\gamma )},`$
where $`𝒟\stackrel{~}{𝐓}^\gamma _{\gamma =1}^Nd\stackrel{~}{𝐓}^\gamma `$, $`P(\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta ,\stackrel{~}{𝐓}^\gamma |\widehat{𝐋}^\gamma )=P(\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta ,\stackrel{~}{𝐓}^\gamma )P(\widehat{𝐋}^\gamma |\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta ,\stackrel{~}{𝐓}^\gamma )/P(\widehat{𝐋}^\gamma )`$ by Bayes’s theorem.

Here we use the approximation of $`P(\widehat{𝐋}^\gamma )=P(\widehat{𝐋}^1,\mathrm{},\widehat{𝐋}^N)_{\gamma =1}^NP(\widehat{𝐋}^\gamma )=1/(4\pi )^N`$. This approximation can be justified as follows: Let $`P(\widehat{𝐋}^1,\widehat{𝐋}^2)`$ be the joint probability distribution of the galaxy angular momentum. One can show that this joint probability distribution can be written as $`P(\widehat{𝐋}^1,\widehat{𝐋}^2)=_{\gamma =1}^2P(\widehat{𝐋}^\gamma )+𝒪(a^2)`$ (see Appendix H). Thus, in the limit of $`a1`$, this approximation holds at first order of $`a`$. In this asymptotic limit, we also have $`P(\widehat{𝐋}^\gamma |\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta ,\stackrel{~}{𝐓}^\gamma )=P(\widehat{𝐋}^\gamma |\stackrel{~}{𝐓}^\gamma )`$ using that the correlation of $`\widehat{𝐋}`$ with $`\stackrel{~}{𝐓}`$ at different points is also $`𝒪(a^2)`$ which we neglect. The constant, $`_{\gamma =1}^NP(\widehat{𝐋}^\gamma )=1/(4\pi )^N`$ is rescaled to unity in equation (C1) since the exact value of the overall constant is irrelevant to the shear reconstruction (any positive proportionality constant can be rescaled to unity).

Furthermore, equation (A8) says that in this asymptotic limit of $`a1`$, $`\widehat{𝐋}`$-dependent part of $`P(\widehat{𝐋}^\gamma |\stackrel{~}{𝐓}^\gamma )`$ is given as

$$P(\widehat{𝐋}^\gamma |\stackrel{~}{𝐓}^\gamma )=a\stackrel{~}{T}_{nk}^\gamma \stackrel{~}{T}_{ko}^\gamma \widehat{L}_n^\gamma \widehat{L}_o^\gamma ,$$
(C2)
apart from a proportionality constant ($`c=5a/3`$). Here the $`\widehat{𝐋}`$-independent part of $`P(\widehat{𝐋}^\gamma |\stackrel{~}{𝐓}^\gamma )`$ is ignored since it does not affect the shear-reconstruction, either. Inserting equation (C2) into equation (C1) gives

$`\stackrel{~}{\xi }_{ijlm}(𝐱_\alpha ,𝐱_\beta )`$ $`=`$ $`a{\displaystyle 𝒟\stackrel{~}{𝐓}^\gamma 𝑑\stackrel{~}{𝐓}^\alpha 𝑑\stackrel{~}{𝐓}^\beta \stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta \stackrel{~}{T}_{nk}^\gamma \stackrel{~}{T}_{ko}^\gamma P(\stackrel{~}{𝐓}^\alpha ,\stackrel{~}{𝐓}^\beta ,\stackrel{~}{𝐓}^\gamma )\widehat{L}_n^\gamma \widehat{L}_o^\gamma },`$ (C3)
$`=`$ $`a{\displaystyle 𝑑𝐱_\gamma \stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta \stackrel{~}{T}_{nk}^\gamma \stackrel{~}{T}_{ko}^\gamma \widehat{L}_n^\gamma \widehat{L}_o^\gamma },`$
$`=`$ $`a{\displaystyle 𝑑𝐱_\gamma \stackrel{~}{C}_{ijnk}(𝐱_\alpha 𝐱_\gamma )\stackrel{~}{C}_{lmok}(𝐱_\beta 𝐱_\gamma )\widehat{L}_n(𝐱_\gamma )\widehat{L}_o(𝐱_\gamma )}.`$
Here $`_\gamma \stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta \stackrel{~}{T}_{nk}^\gamma \stackrel{~}{T}_{ko}^\gamma =_\gamma \stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{nk}^\gamma \stackrel{~}{T}_{lm}^\beta \stackrel{~}{T}_{ko}^\gamma =_\gamma \stackrel{~}{C}_{ijnk}(𝐱_\alpha 𝐱_\gamma )\stackrel{~}{C}_{lmok}(𝐱_\beta 𝐱_\gamma )`$ by the Wick theorem. We ignore the other term $`\stackrel{~}{T}_{ij}^\alpha \stackrel{~}{T}_{lm}^\beta \stackrel{~}{T}_{nk}^\gamma \stackrel{~}{T}_{ko}^\gamma `$ since this term does not depend on the distance of $`𝐱_\alpha 𝐱_\gamma `$ (or $`𝐱_\beta 𝐱_\gamma `$), having no contribution to the shear reconstruction through the galaxy spins. In the continuum limit, the sum is replaced by the integration over $`𝐱_\gamma `$.

## Appendix D MAXIMUM LIKELIHOOD EXPECTATION VALUE

In this appendix, we provide a general argument that the maximum likelihood expectation value of a Gaussian random field can be given as the eigenvector associated with the maximum eigenvalue of the corresponding covariance matrix.

Let $`𝐯`$ be a Gaussian random field. Then, the probability distribution $`P(𝐯)`$ is Gaussian proportional to $`\mathrm{exp}(𝐯^T𝐀^1𝐯/2)`$ where $`𝐀`$ is the covariance matrix of $`𝐯`$. Provided that $`𝐀`$ is positive definite, the maximum likelihood value of $`𝐯`$ must be the one that maximizes $`\mathrm{exp}(𝐯^T𝐀^1𝐯/2)`$ or equivalently the one that minimizes $`(𝐯^T𝐀^1𝐯)/2`$. There is an obvious trivial solution, $`𝐯=0`$ for all points, which is of course not the solution to be sought for.

A nontrivial solution can be found by imposing a constraint. Let us choose a quadratic constraint of $`𝐯^T𝐯=1`$. Then, using the Lagrange multiplier method, we can say that the solution, i.e., the maximum likelihood value of $`𝐯`$ under this constraint should satisfy the following equation:

$$\frac{\delta }{\delta 𝐯}\left(\frac{𝐯^T𝐀^1𝐯}{2}\frac{\lambda }{2}[𝐯^T𝐯1]\right)=0,$$
(D1)
where $`\lambda `$ is a Lagrange multiplier.

Solving the above equation gives

$$𝐀^1𝐯=\lambda 𝐯,𝐯^T𝐯=1.$$
(D2)
Equation (D2) says that the solution to equation (D1) is the eigenvector of $`𝐀^1`$ with the associated eigenvalue of $`\lambda `$. Thus, the eigenvector of $`𝐀^1`$ associated with the smallest eigenvalue minimizes $`(𝐯^T𝐀^1𝐯)/2`$, since $`𝐯^T𝐀^1𝐯=𝐯^T\lambda _{min}𝐯=\lambda _{min}`$. But, the eigenvector of $`𝐀^1`$ associated with the eigenvalue of $`\lambda `$ is also the eigenvector of $`𝐀`$ itself associated with the eigenvalue, $`1/\lambda `$. Therefore, the eigenvector of $`𝐀^1`$ associated with the smallest eigenvalue, $`\lambda _{min}`$ is in fact the eigenvector of $`𝐀`$ associated with the largest eigenvalue, $`1/\lambda _{min}\mathrm{\Lambda }_{max}`$.

Hence, one can say that the maximum likelihood expectation value of $`𝐯`$ is in fact the eigenvector of the positive definite covariance matrix of $`𝐀`$ associated with the largest eigenvalue.

## Appendix E INVERSION THEOREM

This Appendix is devoted fully to prove equation (7), a nontrivial mathematical theorem (inversion theorem), which is at the core of our density reconstruction procedure.

The inversion theorem says the following: If a unit galaxy spin is related to a unit traceless intrinsic shear tensor by equation (3) with a nonzero value of $`a`$, then it is possible to invert the measurable unit galaxy spin field into the initial intrinsic shear field by equation (7). In other words, given the unit spin field, the expected intrinsic shear field is the solution to equation (7) as the eigenvector associated with the largest eigenvalue of the posterior correlation function defined in equation (6).

In order to prove this inversion theorem, we first prove the following three lemmas.

Lemma 1:

$$A_{lm}\underset{i,j}{}\widehat{T}_{ij}\widehat{T}_{jl}\widehat{T}_{im}\stackrel{~}{A}_{lm}=\frac{\widehat{T}_{lm}}{2}.$$
(E1)
In proving Lemma 1, we do not use the Einstein summation rule, so that the repeated indices do not mean the summation in the following proof (but the Einstein summation rule will be recovered after this Lemma 1). Let us first consider the off-diagonal elements, $`A_{lm}`$ with $`lm`$.

$`\stackrel{~}{A}_{lm}`$ $`=`$ $`A_{lm}={\displaystyle \underset{i,j}{}}\widehat{T}_{ij}\widehat{T}_{jl}\widehat{T}_{im}`$ (E2)
$`=`$ $`{\displaystyle \underset{j}{}}\widehat{T}_{lj}\widehat{T}_{jl}\widehat{T}_{lm}+{\displaystyle \underset{j}{}}\widehat{T}_{mj}\widehat{T}_{jl}\widehat{T}_{mm}+{\displaystyle \underset{il,im,j}{}}\widehat{T}_{ij}\widehat{T}_{jl}\widehat{T}_{im}`$
Note that the above equation is correct only in the three dimensional case where there is only one choice among $`1`$,$`2`$, $`3`$ for the dummy index $`i`$, if $`il`$ and $`im`$. So, in the final term of equation (E2), the index $`i`$ is not dummy. Since $`\mathrm{Tr}(\widehat{𝐓})=0`$, we have

$`{\displaystyle \underset{j}{}}\widehat{T}_{ij}\widehat{T}_{jl}`$ $`=`$ $`\widehat{T}_{il}\widehat{T}_{mm}+\widehat{T}_{im}\widehat{T}_{ml},`$
$`{\displaystyle \underset{j}{}}\widehat{T}_{mj}\widehat{T}_{jl}`$ $`=`$ $`\widehat{T}_{ii}\widehat{T}_{ml}+\widehat{T}_{im}\widehat{T}_{il}.`$ (E3)
Using the above equations, one can say

$`{\displaystyle \underset{j}{}}\widehat{T}_{ij}\widehat{T}_{jl}\widehat{T}_{im}+{\displaystyle \underset{j}{}}\widehat{T}_{mj}\widehat{T}_{jl}\widehat{T}_{mm}`$ (E4)
$`=`$ $`\left(\widehat{T}_{il}\widehat{T}_{mm}+\widehat{T}_{im}\widehat{T}_{ml}\right)\widehat{T}_{im}+\left(\widehat{T}_{ii}\widehat{T}_{ml}+\widehat{T}_{im}\widehat{T}_{il}\right)\widehat{T}_{mm}`$
$`=`$ $`\left(\widehat{T}_{im}^2\widehat{T}_{ii}\widehat{T}_{mm}\right)\widehat{T}_{ml}`$
Thus, we have

$`\stackrel{~}{A}_{lm}`$ $`=`$ $`{\displaystyle \underset{j}{}}\widehat{T}_{lj}\widehat{T}_{jl}\widehat{T}_{lm}+\left(\widehat{T}_{im}^2\widehat{T}_{ii}\widehat{T}_{mm}\right)\widehat{T}_{ml}`$ (E5)
$`=`$ $`\left(\widehat{T}_{il}^2+\widehat{T}_{ll}^2+\widehat{T}_{lm}^2+\widehat{T}_{im}^2\widehat{T}_{ii}\widehat{T}_{mm}\right)\widehat{T}_{lm}`$
$`=`$ $`\left(\widehat{T}_{il}^2+\widehat{T}_{lm}^2+\widehat{T}_{im}^2+\widehat{T}_{ii}^2+\widehat{T}_{mm}^2+\widehat{T}_{ii}\widehat{T}_{mm}\right)\widehat{T}_{lm}`$
$`=`$ $`{\displaystyle \frac{1}{2}}\widehat{T}_{lm},`$
since $`\widehat{T}_{ll}^2=\widehat{T}_{ii}^2+\widehat{T}_{mm}^2+2\widehat{T}_{ii}\widehat{T}_{mm}`$ and $`|\widehat{𝐓}|^2=1`$ With the exactly same manner, one can also prove for the diagonal element, $`\stackrel{~}{A}_{ll}=\widehat{T}_{ll}/2`$.

Lemma 2:

$$\stackrel{~}{C}_{ijlm}(r)\stackrel{~}{T}_{ij}(𝐱)\stackrel{~}{T}_{lm}(𝐱+𝐫)\stackrel{~}{C}_{ijkl}(k)=\left(\widehat{k}_i\widehat{k}_j\frac{\delta _{ij}}{3}\right)\left(\widehat{k}_l\widehat{k}_m\frac{\delta _{lm}}{3}\right)P(k),$$
(E6)
where $`P(k)=|\delta _k|^2`$ is the density power spectrum.

By the convolution theorem, we have

$`\stackrel{~}{C}_{ijlm}(k)=\stackrel{~}{T}_{ij}(𝐤)\stackrel{~}{T}_{lm}^{}(𝐤)`$ $`=`$ $`\left(k_ik_j{\displaystyle \frac{\delta _{ij}}{3}}k^2\right)\left(k_lk_m{\displaystyle \frac{\delta _{lm}}{3}}k^2\right)|\mathrm{\Phi }|^2`$ (E7)
$`=`$ $`\left(\widehat{k}_i\widehat{k}_j{\displaystyle \frac{\delta _{ij}}{3}}\right)\left(\widehat{k}_l\widehat{k}_m{\displaystyle \frac{\delta _{lm}}{3}}\right)P(k)`$
since $`T_{ij}(k)=k_ik_j\mathrm{\Phi }(k)`$, $`\mathrm{Tr}(𝐓)=\delta `$, and $`\delta (k)=k^2\mathrm{\Phi }(k)`$.

Lemma 3:

$$\frac{\stackrel{~}{C}_{ijml}(k)}{P(k)}\stackrel{~}{T}_{lm}(k)=\frac{2}{3}\stackrel{~}{T}_{ij}(k).$$
(E8)
$`{\displaystyle \frac{\stackrel{~}{C}_{ijml}(k)}{P(k)}}\stackrel{~}{T}_{lm}(k)`$ $`=`$ $`\left(\widehat{k}_i\widehat{k}_j{\displaystyle \frac{\delta _{ij}}{3}}\right)\left(\widehat{k}_l\widehat{k}_m{\displaystyle \frac{\delta _{lm}}{3}}\right)T_{lm}(k)`$ (E9)
$`=`$ $`\left(\widehat{k}_i\widehat{k}_j{\displaystyle \frac{\delta _{ij}}{3}}\right){\displaystyle \frac{2}{3}}k^2|\mathrm{\Phi }(k)|^2`$
$`=`$ $`{\displaystyle \frac{2}{3}}\stackrel{~}{T}_{ij}(k).`$
Now, we are ready to prove the inversion theorem with the help of the above three lemmas. From here on, we will regard all the proportionality constants as unity. We will use equation (3) as a theoretical estimation formula for $`\widehat{L}_i\widehat{L}_j`$, discarding the shear independent $`\delta _{ij}`$-term in equation (3) since it does not affect the shear inversion:

Inversion Theorem:

$$\widehat{L}_i(𝐱)\widehat{L}_j(𝐱)=\widehat{T}_{il}(𝐱)\widehat{T}_{lj}(𝐱)\frac{\stackrel{~}{\xi }_{abcd}(𝐤_\alpha ,𝐤_\beta )}{P(k_\alpha )P(k_\beta )}\stackrel{~}{T}_{ab}(𝐤_\alpha )𝑑𝐤_\alpha =\stackrel{~}{T}_{cd}(𝐤_\beta ),$$
(E10)
where

$$\stackrel{~}{\xi }_{abcd}(𝐱_\alpha ,𝐱_\beta )=\stackrel{~}{C}_{abil}(𝐱_\alpha 𝐱_\gamma )\stackrel{~}{C}_{cdjl}(𝐱_\beta 𝐱_\gamma )\widehat{L}_i(𝐱_\gamma )\widehat{L}_j(𝐱_\gamma )𝑑𝐱_\gamma .$$
(E11)
By the convolution theorem,

$$\stackrel{~}{\xi }_{abcd}(𝐤_\alpha ,𝐤_\beta )=\stackrel{~}{C}_{abil}(𝐤_\alpha )\stackrel{~}{C}_{cdjl}(𝐤_\beta )\widehat{L}_i(𝐤_\alpha +𝐤_\beta 𝐤^{})\widehat{L}_j(𝐤^{})𝑑𝐤^{}.$$
(E12)
Now, by Lemma 3, we have

$`{\displaystyle \frac{\stackrel{~}{\xi }_{abcd}(𝐤_\alpha ,𝐤_\beta )}{P(k_\alpha )P(k_\beta )}\stackrel{~}{T}_{ab}(𝐤_\alpha )𝑑𝐤_\alpha },`$ (E13)
$`=`$ $`{\displaystyle \frac{\stackrel{~}{C}_{cdjl}(𝐤_\beta )}{P(k_\beta )}\frac{\stackrel{~}{C}_{cdil}(𝐤_\alpha )}{P(k_\alpha )}\stackrel{~}{T}_{ab}(𝐤_\alpha )\widehat{L}_i(𝐤_\alpha +𝐤_\beta 𝐤^{})\widehat{L}_j(k^{})𝑑𝐤^{}𝑑𝐤_\alpha },`$
$`=`$ $`{\displaystyle \frac{\stackrel{~}{C}_{cdjl}(𝐤_\beta )}{P(k_\beta )}\stackrel{~}{T}_{il}(𝐤_\alpha )\widehat{L}_i(𝐤_\alpha +𝐤_\beta 𝐤^{})\widehat{L}_j(𝐤^{})𝑑𝐤^{}𝑑𝐤_\alpha }.`$
Let $`H_{ij}(𝐱)=\widehat{L}_i(𝐱)\widehat{L}_j(𝐱)=\widehat{T}_{in}(𝐱)\widehat{T}_{nj}(𝐱)`$. Then, by the convolution theorem, $`H_{ij}(𝐤_\alpha +𝐤_\beta )=\widehat{L}_i(𝐤_\alpha +𝐤_\beta 𝐤^{})\widehat{L}_j(𝐤^{})𝑑𝐤^{}`$. So equation (E13) is written as

$$\frac{\stackrel{~}{C}_{cdjl}(𝐤_\beta )}{P(k_\beta )}\stackrel{~}{T}_{il}(𝐤_\alpha )H_{ij}(𝐤_\alpha +𝐤_\beta )𝑑𝐤_\alpha =\frac{\stackrel{~}{C}_{cdjl}(𝐤_\beta )}{P(k_\beta )}\stackrel{~}{T}_{il}(𝐤_\alpha )H_{ij}(𝐤_\alpha +𝐤_\beta )𝑑𝐤_\alpha .$$
(E14)
Let us define $`A_{jl}(𝐱)H_{ij}(𝐱)\stackrel{~}{T}_{il}(𝐱)=|\stackrel{~}{𝐓}|\widehat{T}_{in}(𝐱)\widehat{T}_{nj}(𝐱)\widehat{T}_{il}(𝐱)`$. Then in Fourier space one can say $`A_{jl}(𝐤_\beta )=\stackrel{~}{T}_{il}(𝐤_\alpha )H_{ij}(𝐤_\alpha +𝐤_\beta )𝑑𝐤_\alpha `$ by the convolution theorem. Let us decompose $`𝐀`$ into the trace-free part and trace part such that $`A_{jl}(𝐱)=\stackrel{~}{A}_{jl}(𝐱)+\delta _{jl}\mathrm{Tr}(𝐀)/3`$. But we already know from Lemma 1, the trace-free part of $`\widehat{T}_{in}(𝐱)\widehat{T}_{nj}(𝐱)\widehat{T}_{il}(𝐱)`$ is given as $`\widehat{T}_{jl}(𝐱)`$ (apart from the proportionality constant). Therefore, $`\stackrel{~}{A}_{jl}(𝐱)=|\stackrel{~}{𝐓}|\widehat{T}_{jl}(𝐱)=\stackrel{~}{T}_{jl}(𝐱)`$. In Fourier space $`\stackrel{~}{A}_{jl}(𝐤_\beta )=\stackrel{~}{T}_{jl}(𝐤_\beta )`$. Thus, equation (E14) becomes

$`{\displaystyle \frac{\stackrel{~}{C}_{cdjl}(𝐤_\beta )}{P(k_\beta )}}{\displaystyle \stackrel{~}{T}_{il}(𝐤_\alpha )H_{ij}(𝐤_\alpha +𝐤_\beta )𝑑𝐤_\alpha }`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{C}_{cdjl}(𝐤_\beta )}{P(k_\beta )}}A_{jl}(𝐤_\beta )`$ (E15)
$`=`$ $`{\displaystyle \frac{\stackrel{~}{C}_{cdjl}(𝐤_\beta )}{P(k_\beta )}}\left(\stackrel{~}{A}_{jl}(𝐤_\beta )+{\displaystyle \frac{\delta _{jl}}{3}}\mathrm{Tr}(𝐀)\right)`$
$`=`$ $`{\displaystyle \frac{\stackrel{~}{C}_{cdjl}(𝐤_\beta )}{P(k_\beta )}}\stackrel{~}{T}_{jl}(𝐤_\beta ),`$
since $`\stackrel{~}{C}_{cdjl}\delta _{jl}=0`$ by equation (B3).

But, by Lemma 3, equation (E15) is equal to $`\stackrel{~}{T}_{cd}(𝐤_\beta )`$, which finally proves the inversion theorem.

## Appendix F POWER ITERATION

We will provide a general proof for the power iteration scheme in this appendix.

Let us assume that we have a real symmetric positive definite $`n\times n`$ matrix, $`𝐀`$, and we seek for the eigenvector associated with the largest eigenvalue of $`𝐀`$. Let us say, $`𝐯_1,𝐯_2,\mathrm{},𝐯_n`$ are the the $`n`$ eigenvectors of $`𝐀`$ with the associated eigenvalues of $`a_1,a_2,\mathrm{},a_n`$ respectively (here we assume $`a_1a_2\mathrm{}a_n0`$). If $`n`$ is not too large, then we can always find the eigenvectors along with the associated eigenvalues by solving the eigenvector equation, $`\mathrm{𝐀𝐯}_i=a_i𝐯_i`$ numerically. However, in the case $`n`$ is very large, finding the maximum eigenvector by solving the eigenvector equation could be inefficient from a practical point of view since the computational time to solve the eigenvector equation could be too long. The power iteration scheme that we describe and prove here is a practical method to make a fast estimate of the eigenvector, $`𝐯_1`$ associated with the largest eigenvalue, $`a_1`$ fast without solving the eigenvector equation for the case of large $`n`$.

Let us start with an initial arbitrary vector, $`𝐮^0`$. We can construct a new vector, $`𝐮^1`$ out of $`𝐀`$ and $`𝐮^0`$ such that $`𝐮^1=\mathrm{𝐀𝐮}^0`$. Now, using the eigenvectors of $`𝐀`$ as a basis, we can expand $`𝐮^0`$ such that $`𝐮^0=_{i=1}^nb_i𝐯_i`$. So, we can write $`𝐮^1`$ such that $`𝐮^1=𝐀_{i=1}^nb_i𝐯_i=_{i=1}^na_ib_i𝐯_i`$. Iterating this process $`m`$ times leads to a $`m`$-th vector, $`𝐮^m`$ such that $`𝐮^m=_{i=1}^na_ib_i𝐯_i`$. Since $`a_1`$ is the largest eigenvalue, the first component proportional to $`a_1^m`$ dominates. Thus, if we iterate this process sufficiently large times, $`𝐮^m`$ converges effectively to the eigenvector associated with the largest eigenvalue. After $`m`$ iterations, the fractional error caused by approximating $`𝐯_1`$ by $`𝐮^m`$ is proportional to $`(a_1/a_2)^m`$, which goes to zero as $`m`$ becomes large.

Here, the key assumption made for this power iteration to function is that $`𝐀`$ is positive definite, which guarantees $`a_i>0`$ for all $`i=1,n`$. However, even for the case one has a matrix which is not positive definite so that not all eigenvalues are positive, one can still use the power iteration to find the maximum eigenvector as far as the largest eigenvalue is positive. It can be made by inserting secondary steps between each iteration such that

$$𝐮^{m+1/2}=\mathrm{𝐀𝐮}^m,𝐮^{m+1}=\frac{𝐮^{m+1/2}}{|𝐮^{m+1/2}|}+𝐮^{m1},$$
(F1)
with the assumption that not all eigenvalues are negative. After the first iteration, we have $`𝐮^1=_{i=1}^n[1+a_i/(_{j=1}^na_j^2b_j^2)]b_i𝐯_i`$. If $`a_i<0`$, then $`1+a_i/(_{j=1}^na_j^2b_j^2)<1`$. Thus, this refined power iteration effectively suppresses the eigenvectors associated with the negative eigenvalues and converges $`𝐮^m`$ to $`𝐯_1`$.

## Appendix G DENSITY RECONSTRUCTION ALGORITHM

Let $`𝐲(y_0,y_1,\mathrm{},y_n)`$ is a Gaussian random vector of length $`n+1`$. The conditional expectation value of the first component given the rest components, $`y_0|y_1,\mathrm{},y_n`$ can be calculated such that :

$`y_0|y_1,\mathrm{},y_n`$ $`=`$ $`{\displaystyle y_0P(y_0|y_1,\mathrm{},y_n)𝑑y_0},`$
$`=`$ $`{\displaystyle y_0\frac{P(y_0,y_1,\mathrm{},y_n)}{P(y_1,\mathrm{},y_n)}𝑑y_0}.`$
Here, since $`𝐲`$ is Gaussian, $`𝐲^{}(y_1,\mathrm{},y_n)`$ is also Gaussian:

$$P(y_0,\mathrm{},y_n)=\frac{1}{\sqrt{(2\pi )^{n+1}|𝐕|}}\mathrm{exp}\left(\frac{𝐲^T𝐕^1𝐲}{2}\right),$$
(G2)
$$P(y_1,\mathrm{},y_n)=\frac{1}{\sqrt{(2\pi )^n|𝐕^{}|}}\mathrm{exp}\left(\frac{𝐲^T𝐕^1𝐲^{}}{2}\right),$$
(G3)
where $`𝐕`$ is the $`(n+1)\times (n+1)`$ covariance matrix for $`𝐲`$, and $`𝐕^{}`$ is the $`n\times n`$ covariance matrix for $`𝐲^{}`$ such that $`V_{\mu \nu }=V_{\mu \nu }^{}`$ for $`\mu ,\nu =1,\mathrm{},n`$, and $`|𝐕|`$, $`|𝐕^{}|`$ represent the determinants of $`𝐕`$, $`𝐕^{}`$ respectively. In this Appendix, the Greek indices $`\mu ,\nu ,\tau `$ run from $`1`$ to $`n`$.

Let $`𝐔𝐕^1`$. Then, we can rewrite

$`{\displaystyle \frac{𝐲^T𝐕^1𝐲}{2}}`$ $`=`$ $`{\displaystyle \frac{𝐲^T𝐔𝐲}{2}},`$ (G4)
$`=`$ $`{\displaystyle \frac{U_{00}y_0^2}{2}}U_{0\nu }y_\nu y_0{\displaystyle \frac{U_{\mu \nu }y_\mu y_\nu }{2}},`$
$`=`$ $`{\displaystyle \frac{U_{00}}{2}}\left(y_0+{\displaystyle \frac{U_{0\nu }y_\nu }{U_{00}}}\right)^2{\displaystyle \frac{1}{2U_{00}}}\left(U_{00}U_{\mu \nu }U_{0\nu }U_{0\nu }\right)y_\mu y_\nu ,`$
But, we have:

$$\frac{1}{U_{00}}\left(U_{00}U_{\mu \nu }U_{0\mu }U_{0\nu }\right)=V_{\mu \nu }^1.$$
(G5)
which can be proved by

$`V_{\mu \nu }^{}{\displaystyle \frac{(U_{00}U_{\nu \tau }U_{0\nu }U_{0\tau })}{U_{00}}}`$ $`=`$ $`V_{\mu \nu }^{}U_{\nu \tau }{\displaystyle \frac{V_{\mu \nu }U_{\nu 0}U_{0\tau }}{U_{00}}},`$ (G6)
$`=`$ $`\delta _{\mu \tau }{\displaystyle \frac{\delta _{\mu 0}U_{0\tau }}{U_{00}}}=\delta _{\mu \tau },`$
since $`\mu 0`$, i.e., $`\delta _{\mu 0}=0`$.

Therefore, we can express the conditional probability density distribution,
$`P(y_0|y_1,\mathrm{},y_n)`$ such that

$$P(y_0|y_1,\mathrm{},y_n)=\frac{1}{\sqrt{2\pi |𝐕||𝐕^{}|^1}}\mathrm{exp}\left\{\frac{U_{00}}{2}\left(y_0+\frac{U_{0\nu }y_\nu }{U_{00}}\right)^2\right\}.$$
(G7)
Thus, equation (LABEL:eqn:cond) can be simplified into

$`y_0|y_1,\mathrm{},y_n`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi |𝐕||𝐕^{}|^1}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}y_0\mathrm{exp}\left\{{\displaystyle \frac{U_{00}}{2}}\left(y_0+{\displaystyle \frac{U_{0\nu }y_\nu }{U_{00}}}\right)^2\right\}𝑑y_0,`$ (G8)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{|𝐕||𝐕^{}|^1}}}{\displaystyle \frac{U_{0\nu }y_\nu }{U_{00}^{3/2}}}.`$
To calculate $`|𝐕||𝐕^{}|^1`$, let us construct a $`(n+1)\times (n+1)`$ matrix, say $`𝐕^{\prime \prime }`$ from the $`n\times n`$ maxtrix $`𝐕^{}`$ such that

$$𝐕^{\prime \prime }\left(\begin{array}{cc}1,& \mathrm{𝟎}\\ \mathrm{𝟎},& 𝐕\end{array}\right).$$

Then, obviously $`|𝐕^{\prime \prime }|=|𝐕^{}|`$, i.e., $`|𝐕^{\prime \prime }|^1=|𝐕^{}|^1`$.
Thus one can say $`|𝐕||𝐕^{}|^1=|𝐕||𝐕^{\prime \prime }|^1=|𝐕||𝐕^{\prime \prime 1}|=|𝐕𝐕^{\prime \prime 1}|`$. But, $`𝐕𝐕^{\prime \prime 1}`$ is a $`(n+1)\times (n+1)`$ matrix such that

$$𝐕𝐕^{\prime \prime 1}=\left(\begin{array}{cc}V_{00},& 𝐕_{0\tau }\\ 𝐕_{\nu 0},& 𝐈\end{array}\right),$$

where $`𝐈`$ is a $`n\times n`$ identity matrix. Its determinant is straightforwardly calculated to be

$$|𝐕𝐕^{\prime \prime 1}|=V_{00}V_{0\nu }V_{0\tau }V_{\nu \tau }^1.$$
(G9)
Hence, we finally find equation (G8) equal to

$$y_0|y_1,\mathrm{},y_n=\frac{U_{0\nu }y_\nu }{U_{00}},$$
(G10)
since through equations (G5) and (G9),

$`U_{00}|𝐕||𝐕^{}|^1`$ $`=`$ $`U_{00}(V_{00}V_{0\nu }V_{0\tau }V_{\nu \tau }^1),`$ (G11)
$`=`$ $`U_{00}V_{00}V_{0\nu }V_{0\tau }\left(U_{00}U_{\nu \tau }U_{0\nu }U_{0\tau }\right),`$
$`=`$ $`1.`$
We emphasize that this expression for $`y_0|y_1,\mathrm{},y_n`$ given as equation (G10) here is equivalent to equation (10) in Betschinger (1987). Bertschinger’s formula is written as
$`y_0|𝐲^{}=V_{0\nu }V{}_{\nu \nu ^{}}{}^{1}y_{\nu ^{}}^{}`$. Let us reexpress the matrix $`𝐕`$ such that

$$𝐕=\left(\begin{array}{cc}V_{00},& 𝐕_{0i}\\ 𝐕_{i0},& 𝐕_{ij}\end{array}\right).$$

Then, one can say $`𝐕_{0i}^1/𝐕_{00}^1=|𝐀^{0i}|`$, where $`𝐀`$ represents a cofactor of $`𝐕`$. Now the
determinant of the cofactor can be expressed as a sum such that $`|𝐁|=_i𝐁_{0i}|𝐂^{0i}|`$ where $`𝐁_{ij}𝐀^{0i}`$ and $`𝐂`$ is the cofactor of $`𝐁`$. This is explicitly the same as Betschinger’s formula, when we expand each of Betschinger’s matrix elements as the determinants of cofactors, and sum it over $`𝐕_{0\nu }`$.

However, there is one advantage of our formula over Bertschinger’s: Our conditional expectation value is often computationally much cheaper, since for a translation invariant random field, $`𝐕^1`$ can be computed with the fast fourier transformation method (FFT), while that is not generally possible for Bertschinger’s formula.

## Appendix H SPIN-SPIN CORRELATION

Let us first consider the three dimensional spatial spin-spin correlation,
$`|\widehat{𝐋}(𝐱)\widehat{𝐋}(𝐱+𝐫)|^2\widehat{L}_i\widehat{L}_i^{}\widehat{L}_j\widehat{L}_j^{}`$. Before estimating it theoretically, it is instructive to understand the interpretation of $`\widehat{𝐋}`$ in this expression. In practice, $`\widehat{𝐋}`$ is the measured unit spin vector of an observed galaxy. In theory, however, there is no way to calculate the unit spin vector of an observed galaxy analytically since galaxy formation is still an unsolved problem. What one can do at most (and at best) theoretically is to use some analytic estimation formula for the expectation value of $`\widehat{𝐋}`$. Here we use equation (3) as an theoretical formula for $`\widehat{𝐋}`$ based on the linear perturbation theory.

In equation (3) we have taken the ensemble average over the inertia tensors to obtain a result that depends on shear tensors. Taking the ensemble average of equation (3) over the shear tensors as well would result in a trivial $`\widehat{L}_i\widehat{L}_j=\delta _{ij}/3`$.

$`\widehat{L}_i\widehat{L}_i^{}\widehat{L}_j\widehat{L}_j^{}`$ $`=`$ $`\widehat{L}_i\widehat{L}_j\widehat{L}_i^{}\widehat{L}_j^{},`$ (H1)
$`=`$ $`\left({\displaystyle \frac{1+a}{3}}\delta _{ij}a\widehat{T}_{ik}\widehat{T}_{kj}\right)\left({\displaystyle \frac{1+a}{3}}\delta _{ij}a\widehat{T}_{il}^{}\widehat{T}_{lj}^{}\right),`$ (H2)
$`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{a^2}{3}}+a^2\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{il}^{}\widehat{T}_{lj}^{}.`$ (H3)
We took the expectation value over the inertia tensors going from (H1) to (H2) assuming the inertia tensors at two different positions to be independent.

In equation (H3) it is formidable to calculate $`\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{il}^{}\widehat{T}_{lj}^{}`$ analytically since $`\widehat{𝐓}`$ is in general not a Gaussian random field while $`\stackrel{~}{𝐓}`$ is. We approximate $`\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{il}^{}\widehat{T}_{lj}^{}=(\stackrel{~}{T}_{ik}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{il}^{}\stackrel{~}{T}_{lj}^{})/(|\stackrel{~}{𝐓}|^2|\stackrel{~}{𝐓}^{}|^2)`$ by $`\stackrel{~}{T}_{ik}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{il}^{}\stackrel{~}{T}_{lj}^{}/|\stackrel{~}{𝐓}|^2^2`$, and apply the Wick theorem such that

$`\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{il}^{}\widehat{T}_{lj}^{}`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{T}_{ik}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{il}^{}\stackrel{~}{T}_{lj}^{}}{|\stackrel{~}{𝐓}|^2|\stackrel{~}{𝐓}^{}|^2}},`$ (H4)
$``$ $`{\displaystyle \frac{\stackrel{~}{T}_{ik}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{il}^{}\stackrel{~}{T}_{lj}^{}}{|\stackrel{~}{𝐓}|^2|\stackrel{~}{𝐓}^{}|^2}},`$
$`=`$ $`{\displaystyle \frac{9}{4\xi ^2(0)}}\left(\stackrel{~}{T}_{ik}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{il}^{}\stackrel{~}{T}_{lj}^{}+\stackrel{~}{T}_{ik}\stackrel{~}{T}_{il}^{}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{lj}^{}+\stackrel{~}{T}_{ik}\stackrel{~}{T}_{lj}^{}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{il}^{}\right),`$
since $`|\stackrel{~}{𝐓}|^2^2=4\xi ^2(0)/9`$ by $`\stackrel{~}{T}_{ij}\stackrel{~}{T}_{kl}=(3\delta _{ik}\delta _{jl}+3\delta _{il}\delta _{jk}2\delta _{ij}\delta _{kl})\xi (0)/45`$ Bardeen et al. (1986).

The $`r`$-independent term in equation (H4) is straightforwardly calculated such that

$$\frac{9\stackrel{~}{T}_{ik}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{il}^{}\stackrel{~}{T}_{lj}^{}}{4\xi ^2(0)}=\frac{1}{3}.$$
(H5)
The rest $`r`$-dependent two terms in equation (H4) can be calculated using the given $`\stackrel{~}{C}_{ijkl}`$ (Appendix B). We find

$`{\displaystyle \frac{9}{4\xi ^2(0)}}\left(\stackrel{~}{T}_{ik}\stackrel{~}{T}_{il}^{}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{lj}^{}+\stackrel{~}{T}_{ik}\stackrel{~}{T}_{lj}^{}\stackrel{~}{T}_{kj}\stackrel{~}{T}_{il}^{}\right)`$ $`=`$ $`{\displaystyle \frac{9}{4\xi ^2(0)}}\left(\stackrel{~}{C}_{ikil}\stackrel{~}{C}_{kjlj}+\stackrel{~}{C}_{iklj}\stackrel{~}{C}_{kjil}\right)`$ (H6)
$`=`$ $`{\displaystyle \frac{9}{4\xi ^2(0)}}\stackrel{~}{C}_{iklj}\stackrel{~}{C}_{kjil},`$
since $`\stackrel{~}{C}_{ikil}=\stackrel{~}{C}_{kjlj}=0`$. Now using equation (B3), one can show that

$`\stackrel{~}{C}_{iklj}\stackrel{~}{C}_{kjil}`$ $`=`$ $`{\displaystyle \frac{4}{9}}\xi ^2+{\displaystyle \frac{14}{5}}J_5^24J_3J_5+{\displaystyle \frac{14}{9}}J_3^2{\displaystyle \frac{8}{5}}J_5\xi +{\displaystyle \frac{8}{9}}J_3\xi .`$ (H7)
For the case of a power law spectrum of $`\xi (r)r^n`$, every term in equation (H7) is proportional to $`\xi ^2(r)`$ since $`J_3=\frac{3}{n+3}r^n\xi (r)`$ and $`J_5=\frac{5}{n+5}r^n\xi (r)`$. Thus, one can say that for the case of a power law spectrum, equation (H6) is proportional to $`\xi ^2(r)`$. Let us say equation (H6) can be written as $`\frac{1}{3}+\beta \frac{\xi ^2(r)}{\xi ^2(0)}`$ where $`\beta `$ is a proportionality constant to be determined. Note that it is true only for the case of a power-law spectrum.

In order to optimize the approximation used in equation (H4), we determine the optimal value of the proportionality constant by considering the real value of $`\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{il}^{}\widehat{T}_{lj}^{}`$ in the limit of $`r=0`$. In the asymptotic limit of $`r=0`$, we have

$$\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{il}^{}\widehat{T}_{lj}^{}=\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{il}\widehat{T}_{lj}=\frac{1}{2}\widehat{T}_{lj}\widehat{T}_{lj}=\frac{1}{2},$$
(H8)
by equation (E1). But we also have

$$\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{il}^{}\widehat{T}_{lj}^{}\frac{1}{3}+\beta \frac{\xi ^2(r)}{\xi ^2(0)}.$$
(H9)
By equating equation (H9) with $`r=0`$ to (H8), we find the best-approximation proportionality constant, $`\beta =1/6`$.

So far we have treated a galaxy as a point-like object. However, a real galaxy has a finite size with a typical one dimensional Lagrangian scale of $`0.55h^1`$ Mpc. Dealing with a real galaxy with a finite size amounts to replacing $`\xi `$ by $`\xi _R`$ (a top-hat convolved density correlation on a galaxy scale of $`R`$). Finally, we find the approximate estimation of the three dimensional spatial spin-spin correlation through equations (H3) to (H6) such that

$`\widehat{L}_i\widehat{L}_i^{}\widehat{L}_j\widehat{L}_j^{}`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{a^2}{3}}+a^2\widehat{T}_{ik}\widehat{T}_{kj}\widehat{T}_{il}^{}\widehat{T}_{lj}^{}`$ (H10)
$``$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{a^2}{3}}+a^2\left({\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{6}}\xi _R^2(r)\right)`$
$``$ $`{\displaystyle \frac{1}{3}}+{\displaystyle \frac{a^2}{6}}\xi _R^2(r),`$
with the normalization of $`\xi _R(0)=1`$.

Note that the only approximation made in the derivation of equation (H10) is to replace $`\widehat{T}\widehat{T}\widehat{T}^{}\widehat{T}^{}`$ with $`\stackrel{~}{T}\stackrel{~}{T}\stackrel{~}{T}^{}\stackrel{~}{T}^{}/|\stackrel{~}{T}|^2^2`$. We have tested the validity of this approximation by Monte-Carlo simulation, and found that it holds good for all value of $`r`$. We generated a three dimensional Gaussian shear field with the density correlation function $`\xi (r)r^1`$, and smoothed it with a top-hat window on a scale radius of two grids. We then generate the unit spins by multiplying the shear tensors with the uncorrelated inertia tensors with an approximate normalization. We then plot the LHS of equation (H10) using the randomly generated random fields shown as the solid line in Fig. 3, while the RHS derived from the density correlation function is shown as the dotted line. In this model we have used $`a=3/5`$. We note the excellent agreement between the model and the simulation results, showing that the application of Wick’s theorem to the unit shears in equation (H4) is well justified.

With the similar manner one can also find a two dimensional spin-spin correlation, $`|\widehat{𝐒}(𝐱)\widehat{𝐒}(𝐱+𝐫)|^2`$. For the case of a power law spectrum, all quadratic statistics in the spin should depend on the square of the density correlation function. It is reasonable to assume that the two dimensional spin-spin correlation should be also expressed as $`1/2+A\xi _R^2(r)`$ for powerlaw correlations. Here $`1/3`$ is replaced by $`1/2`$ since $`|\widehat{𝐒}(𝐱)\widehat{𝐒}(𝐱+𝐫)|^2=1/2`$ for the two dimensional random spins.

Through equations (A12) and (A13), we have $`\widehat{S}_1^2\widehat{S}_2^2|\widehat{𝐓}=5a(\widehat{T}_{1k}\widehat{T}_{k1}\widehat{T}_{2k}\widehat{T}_{k2})/4`$, $`\widehat{S}_1\widehat{S}_2|\widehat{𝐓}=5a\widehat{T}_{1k}\widehat{T}_{k2}/4`$. The comparison of this two dimensional expression with equation (3) indicates that the shear-spin correlation parameter $`a`$ should be replaced by $`5a/4`$ for the two dimensional case. In other words, the amplitude of $`A=a^2/6`$ for the three dimensional case should be replaced by $`A=(5a/4)^2/6`$ for the two dimensional case. Thus, we have

$$\widehat{S}_i\widehat{S}_j\widehat{S}_i^{}\widehat{S}_j^{}\frac{1}{2}+\frac{25}{96}a^2\xi _R^2(r).$$
(H11)
## Appendix I SPIN-DIRECTION CORRELATION

The same technique used in Appendix H can be also applied to calculate the spin-direction correlation, $`\widehat{L}_i\widehat{L}_j\widehat{r}_i\widehat{r}_j`$ such that

$`\widehat{L}_i\widehat{L}_j\widehat{r}_i\widehat{r}_j`$ $`=`$ $`\left({\displaystyle \frac{1+a}{3}}\delta _{ij}a\widehat{T}_{ik}^R\widehat{T}_{kj}^R\right)\left({\displaystyle \frac{1b}{3}}\delta _{ij}+b\widehat{T}_{ik}^R^{}\widehat{T}_{kj}^R^{}\right),`$ (I1)
$`=`$ $`{\displaystyle \frac{1}{3}}+{\displaystyle \frac{ab}{3}}ab\widehat{T}_{ik}^R\widehat{T}_{kj}^R\widehat{T}_{il}^R^{}\widehat{T}_{lj}^R^{}.`$ (I2)
Using the same approximation that is used in equation (H4) and verified by Monte-Carlo simulations (see Fig. 3), we have

$`{\displaystyle \frac{9\widehat{T}_{ik}^R\widehat{T}_{kj}^R\widehat{T}_{il}^R^{}\widehat{T}_{lj}^R^{}}{4(\sigma _R\sigma _R^{})^2}}`$ $`=`$ $`{\displaystyle \frac{9\stackrel{~}{T}_{ik}^R\stackrel{~}{T}_{kj}^R\stackrel{~}{T}_{il}^R^{}\stackrel{~}{T}_{lj}^R^{}}{4(\sigma _R\sigma _R^{})^2}},`$ (I3)
$`=`$ $`{\displaystyle \frac{9(\stackrel{~}{T}_{ik}^R\stackrel{~}{T}_{kj}^R\stackrel{~}{T}_{il}^R^{}\stackrel{~}{T}_{lj}^R^{}+\stackrel{~}{T}_{ik}^R\stackrel{~}{T}_{il}^R^{}\stackrel{~}{T}_{kj}^R\stackrel{~}{T}_{lj}^R^{}+\stackrel{~}{T}_{ik}^R\stackrel{~}{T}_{lj}^R^{}\stackrel{~}{T}_{kj}^R\stackrel{~}{T}_{il}^R^{})}{4(\sigma _R\sigma _R^{})^2}}.`$
Note that we calculate the shear correlation at the same position but smoothed on two different scales. Following the same logic explained in Appendix H, we have
$`9\stackrel{~}{T}_{ik}^R\stackrel{~}{T}_{kj}^R\stackrel{~}{T}_{il}^R^{}\stackrel{~}{T}_{lj}^R^{}/(4\sigma _R^2\sigma _R^{}^2)=1/3`$,
$`9(\stackrel{~}{T}_{ik}^R\stackrel{~}{T}_{il}^R^{}\stackrel{~}{T}_{kj}^R\stackrel{~}{T}_{lj}^R^{}+\stackrel{~}{T}_{ik}^R\stackrel{~}{T}_{lj}^R^{}\stackrel{~}{T}_{kj}^R\stackrel{~}{T}_{il}^R^{})/(4\sigma _R^2\sigma _R^{}^2)\delta ^R(𝐱)\delta ^R^{}(𝐱)^2/(\sigma _R^2\sigma _R^{}^2)`$. The proportionality constant can be easily obtained to be $`1/6`$ again by considering the limit of $`r=0`$ (i.e., $`R^{}=R`$).

So, finally we find

$$\widehat{L}_i\widehat{L}_j\widehat{r}_i\widehat{r}_j=\frac{1}{3}+\frac{ab}{6}\frac{\delta _R^{}\delta _R^2}{(\sigma _R^{}\sigma _R)^2}.$$
(I4)
For the two dimensional case, replacing $`1/3`$ with $`1/2`$, and $`a`$ with $`5a/4`$, we also find

$$\widehat{S}_i\widehat{S}_j\widehat{r}_i\widehat{r}_j=\frac{1}{2}+\frac{5ab}{24}\frac{\delta _R^{}\delta _R^2}{(\sigma _R^{}\sigma _R)^2}.$$
(I5)
## Appendix J CORRELATION PARAMETERS

In this final appendix, we provide the derivations of the optimal estimation formula (eq.  and eq. ) for the two correlation parameters ($`a`$ and $`b`$ respectively), and the involved error-bar formula as well.

Multiplying $`\widehat{\lambda }_i`$ to each side of equation (3), and using $`_i\widehat{\lambda }_i=0`$, $`_i\widehat{\lambda }_i^2=1`$, $`_i\widehat{\lambda }_i^4=1/2`$, $`_{i<j}\widehat{\lambda }_i\widehat{\lambda }_j=1/2`$, and $`_{i<j}\widehat{\lambda }_i^2\widehat{\lambda }_j^2=1/4`$, equation (4) is straightforwardly derived from equation (3). Here note that equation (3) is used as a theoretical estimation formula for $`\widehat{L}_i\widehat{L}_j`$. The error $`ϵ_a`$ involved in the measurement of the average value of $`a`$ is given as the standard deviation of $`a`$ for the case of no correlation between $`\widehat{𝐓}`$ and $`\widehat{𝐋}`$.

For the case of no correlation, $`a=0`$. So, $`ϵ_a=\sqrt{a^2}`$. Now, by equation (4) we have

$$ϵ_a=\sqrt{a^2}=\sqrt{(26\underset{i}{}\widehat{\lambda }_i^2\widehat{L}_i^2)^2}=\sqrt{\frac{4}{5}},$$
(J1)
since $`\widehat{L}_i^2=1/3`$, $`\widehat{L}_i^4=1/5`$ for each $`i=1,2,3`$ and $`\widehat{L}_i^2\widehat{L}_j^2=1/15`$ for each $`ij`$ if $`\widehat{𝐋}`$ is random. Thus, for the $`N_t`$ ensemble, we finally have $`ϵ_a=\sqrt{4/(5N_t)}`$.

One can derive an optimal formula for $`b`$ with the similar argument. But, there is one notable difference in deriving the optimal formula for $`b`$. For $`a`$, we have directly used equation (3) to derive equation (4) where each $`|\widehat{L}_i|^2`$ is weighed by the square of the eigenvalue, $`\widehat{\lambda }_i^2`$. It is adequate since the unit galaxy spin vector, $`\widehat{𝐋}`$, is expected to be aligned with the intermediate principal axis of the shear tensor which does not involve any degeneracy upon a sign change of the shear tensor ($`\widehat{𝐓}\widehat{𝐓}`$). Whereas, the unit galaxy separation vector, $`\widehat{𝐫}`$, is expected to align with the major principal axis of the shear tensor, and thus orthogonal to the minor principal axis. But the major and minor principal axes of the shear tensor are interchanged by a sign change of the shear tensor.

In order to measure the real correlation between $`\widehat{𝐫}`$ and the major principal axis of the shear tensor, we first define a secondary correlation parameter, $`d`$, by weighing each $`\widehat{r}_i`$ by $`\widehat{\lambda }_i`$ instead of $`\widehat{\lambda }_i^2`$ such that

$$d=\widehat{\lambda }_i\widehat{r}_i^2.$$
(J2)
In the following we show that in fact $`b=d/\sqrt{2}`$, and prove that the optimal estimation formula for $`b`$ is indeed equation (13).

Let $`𝐞^M`$ be the major eigenvector of $`\widehat{𝐓}^R^{}`$. Since $`\widehat{𝐫}`$ is supposed to be aligned with $`𝐞^M`$, one can model the unit galaxy separation vector as a mixture of the major eigenvector with a random component such that

$$\widehat{r}_i\widehat{r}_j|\widehat{𝐓}^R^{}=\gamma ^2e_i^Me_j^M+(1\gamma ^2)\frac{\delta _{ij}}{3}$$
(J3)
where $`\gamma ^2`$ measures the strength of the correlation between $`\widehat{𝐫}`$ and $`𝐞^M`$. Using $`𝐞^M=(1,0,0)`$ in the principal axis frame of $`\widehat{𝐓}^R^{}`$, one gets $`\widehat{r}_1^2|\widehat{𝐓}^R^{}=\frac{1}{3}+\frac{2\gamma ^2}{3}`$, $`\widehat{r}_2^2|\widehat{𝐓}^R^{}=\frac{1}{3}\frac{\gamma ^2}{3}`$, and $`\widehat{r}_2^2|\widehat{𝐓}^R^{}=\frac{1}{3}\frac{\gamma ^2}{3}`$. Now, using equation (J3) for the theoretical estimation of $`|\widehat{r}_i|^2`$ in equation (J2), we have

$$d=\gamma ^2\left(\frac{2}{3}\widehat{\lambda }_1\frac{1}{3}\widehat{\lambda }_2\frac{1}{3}\widehat{\lambda }_3\right).$$
(J4)
Approximating the three eigenvalues of the traceless shear tensor as $`\widehat{\lambda }_1=1/\sqrt{2}=\widehat{\lambda }_3`$ and $`\widehat{\lambda }_3=0`$ since $`\widehat{\lambda }_1^2=\widehat{\lambda }_3^2\widehat{\lambda }_2^2`$ and $`\widehat{\lambda }_i\widehat{\lambda }_i=1`$, and inserting these values of the three eigenvalues into equation (J4), we find $`\gamma ^2=d\sqrt{2}`$.

Now, let us model the unit spin vector as a mixture of the intermediate eigenvector of $`\widehat{𝐓}^R`$, $`𝐞^I`$, with a random component such that

$$\widehat{L}_i\widehat{L}_j|\widehat{𝐓}^R=\alpha ^2e_i^Ie_j^I+(1+\alpha ^2)\frac{\delta _{ij}}{3},$$
(J5)
where $`\alpha ^2`$ measures the correlation strength between $`\widehat{𝐋}`$ and $`𝐞^I`$. Using $`𝐞^I=(0,1,0)`$ in the principal axis frame of $`\widehat{𝐓}^R`$, one gets $`\widehat{L}_1^2|\widehat{𝐓}^R=\frac{1}{3}+\frac{\alpha ^2}{3}`$, $`\widehat{L}_2^2|\widehat{𝐓}^R=\frac{1}{3}\frac{2\alpha ^2}{3}`$, and $`\widehat{L}_2^2|\widehat{𝐓}^R=\frac{1}{3}+\frac{\alpha ^2}{3}`$. Inserting the above approximate values of the eigenvalues into equation (J5), we have $`\alpha ^2=a/2`$.

Extrapolating equation (J3) to the limit of $`R^{}=R`$, let us calculate $`\widehat{L}_i\widehat{L}_j|\widehat{𝐓}^R\widehat{r}_i\widehat{r}_j|\widehat{𝐓}^R`$ at this asymptotic limit of $`R^{}=R`$. Using equations (I4), (J3) and (J5) for the theoretical estimation formula, and given $`𝐞^M𝐞^I=0`$ due to the orthogonality of the eigenvectors, we have

$$\widehat{L}_i\widehat{L}_j|\widehat{𝐓}^R\widehat{r}_i\widehat{r}_j|\widehat{𝐓}^R^{}=\frac{1}{3}\frac{\alpha ^2\gamma ^2}{3}=\frac{1}{3}+\frac{ab}{6}.$$
(J6)
It shows that $`b=d\sqrt{2}`$ since $`\alpha ^2=\frac{a}{2}`$ and $`\gamma ^2=\sqrt{2}d`$.

The error $`ϵ_b`$ involved in the measurement of the average of $`b`$ is also given as the standard deviation of $`b`$ for the case of no correlation between $`\widehat{𝐓}`$ and $`\widehat{𝐫}`$. Now we have

$$ϵ_b=\sqrt{b^2}=\sqrt{2\underset{i,j}{}\widehat{\lambda }_i\widehat{\lambda }_j\widehat{r}_i^2\widehat{r}_j^2}=\sqrt{\frac{4}{15}}.$$
(J7)
Thus, for the $`N_t`$ ensemble, $`ϵ_b=\sqrt{4/(15N_t)}`$.