A New Estimator for Phase Statistics

We introduce a novel statistic to probe the statistics of phases of Fourier modes in two-dimensions (2D) for weak lensing convergence field $\kappa$. This statistic contains completely independent information compared to that contained in observed power spectrum. We compare our results against state-of-the-art numerical simulations as a function of source redshift and find good agreement with theoretical predictions. We show that our estimator can achieve better signal-to-noise compared to the commonly employed statistics known as the line correlation function (LCF). Being a two-point statistics, our estimator is also easy to implement in the presence of complicated noise and mask, and can also be generalised to higher-order. While applying this estimator for the study of lensed CMB maps, we show that it is important to include post-Born corrections in the study of statistics of phase.


Introduction
The weak lensing surveys which include Canada-France-Hawaii Telescope(CFHTLS) 1 , PAN-STARRS 2 , Dark Energy Surveys (DES) 3 [1], Prime Focus Spectrograph 4 , WiggleZ 5 [2], BOSS 6 [3], KiDS [4] and Subaru Hypersuprimecam survey 7 (HSC) [5] are already providing important cosmological insights. The future large scale structure (LSS) surveys Euclid 8 [6], Rubin Observatory 9 [7] and Roman Space Telescope 10 list weak lensing as their main science driver and are expected to take us beyond the standard model of cosmology [8] by answering some of the most profound questions regarding the nature of dark matter and dark energy (equivalently the modified theories of gravity) [9,10] and nature of neutrino mass hierarchy [11].
Weak lensing observations target the low-redshift universe and small scales where the perturbations are in the nonlinear regime and statistics are non-Gaussian [12]. Many different estimators exist which probe the higher-order statistics of weak lensing maps [13]. These include the well known real-space one-point statistics such as the cumulants [14] or their two-point correlators also known as the cumulant correlators as well as the associated PDF [15] and the peak-count statistics [16]. In the harmonic domain the estimators such as the Skew-Spectrum [17], Integrated Bispectrum [18] kurt-spectra [19], morphological estimator [20], integrated tripsectrum [21], Betti number [22], extreme value statistics [23], positiondependent PDF [24], density split statistics [25], response function formalism [26] estimators for shapes of the lensing bispectrum [27] are some of the statistical estimators and formalism recently considered by various authors in the context of understanding cosmological statistics in general and weak lensing in particular.
The gravitational clustering in the quasilinear and nonlinear regime generates coupling of Fourier modes that results in correlation of their phases [28]. The power spectrum does not contain any phase information [29]. Notable initial work on statistical evolution of phases of the Fourier modes in gravitational clustering in the quasilinear regime can be found in [30,31]. These studies relied on a perturbative framework and were tested against numerical simulations [32]. A universal behavior in evolution of phases of the Fourier modes in the nonlinear regime was reported in [33].
In recent years an estimator known also as the line correlation function (LCF) to measure for three-point (third-order) phase correlations was introduced in [34] (also see [35]). This was used in many different contexts. The possibility of improving the cosmological constraints using phase correlations were studied in [36,37]. In the context of redshift-space distortions this was used in [38]. The growth rate of perturbation were probed in [39]. One of the motivation in this paper is to introduce the third-order phase correlation functions for projected surveys in general and weak lensing surveys in particular. We introduce third-order estimators using both two-and three-point statistics.
We note here in passing that in addition to the summary statistics listed above, in recent years many novel techniques have gained popularity. These include Bayesian hierarchical modelling, likelihood-free or forward modelling approaches. However, it is important to realise that many of these methods often rely on simulations that are based on lognormal approximation or high-order Lagrangian theories and are only approximate compared to accurate ray tracing simulations which can be rather expensive (however, see also [40,41]).
At leading order the LCF takes contributions from the bispectrum. However, at smaller separations, it also take contribution from higher-order statistics. The perturbative treatment breaks down at smaller separation. The LCF encodes information that is highly complementary to that contained in the power spectrum. Next, we will introduce a two-point statistics also known as cumulant-correlator which can probe phase bispectrum with a higher signalto-noise.
The LCF has also been employed to distinguish various morphological types of collapsed objects [35]. In addition to the LCF other real space triangular configurations (TCF) have been considered for the study of third-order phase statistics, e.g., [42] employed TCF to probe the characteristic scale of ionized regions during the epoch of reionization from 21cm interferometric observations. This paper is organised as follows. In §2 we introduce the weak lensing bispectrum. In §3 we briefly review the modelling of the third-order phase statistics. In §4 we discuss our results. The conclusions are drawn in §5.

Weak Lensing Bispectrum
The projected weak lensing convergence κ is a line-of-sight integration of the underlying threedimensional (3D) cosmological density contrast δ. The κ(θ) at a position θ can be expressed as follows: Here d A (r) is the comoving angular diameter distance at a comoving distance r. The kernel w(r) encodes geometrical dependence; a is the scale factor, H 0 is the hubble constant and Ω M is the cosmological density parameter. d A (r) and d A (r s ) are comoving angular diameter distances at a comoving distances r and r s . We have assumed all sources to be at a single source plane at a distance r s . The power spectrum P κ (l) and the bispectrum B κ (l 1 , l 2 , l 3 ) for the convergence maps are defined through the following expression : The two dimensional Fourier transform of κ(θ) is denoted as κ(l) with l being the twodimensional wave vector. we denote the two-dimensional Dirac delta function as δ 2D . We will use P δ and B δ to denote the power spectrum and bispectrum respectively of the underlying three dimensional cosmological density contrast δ. In the tree level standard perturbation theory the bispectrum B δ (k 1 , k 2 , k 3 ) can be expressed in terms of the kernel F 2 (k 1 , k 2 ) and power spectrum as follows [28]: The perturbative bispectrum as B PT δ which depends on the linear power spectrum P L . Here, k i represents 3D wave vectors and their moduli are represented as k i . In the flat-sky approximation, the convergence power spectrum P κ and bispectrum B κ can be expressed in terms of the underlying power spectrum P δ and B δ using the following line-of-sight integrations: , ; r . (2.4b) In the highly nonlinear regime many different halo-model based fitting functions have been proposed. We will use the most recent fitting function presented in Ref. [43] known to be more accurate compared to the previous fitting functions. We will use these expressions in our calculation for the Line Correlation Function (LCF) for the κ field. Notice that bispectrum only represents the leading contribution. Higher-order correction also get contributions from higher-order statistics such as the trispectrum which we have ignored in our study, For the validation of our theoretical results we use simulation that adopted cosmological parameters consistent with the WMAP 9 year result Ω m = 1 − Ω Λ = 0.279, Ω cdm = 0.233, Ω b = 0.046, h = 0.7, σ 8 = 0.82 and n s = 0.97

Third-order Phase Statistics in Projection
We will introduce the two-point correlation function ξ(θ) of the weak lensing convergence κ through the following expression: Here P l is the Legendre polynomial and l represents the angular harmonics. Also, J 0 represents the Bessel Function of order zero of the first kind. The angular power spectrum C l of convergence κ is identical to P κ (l) for the flat-sky approximation. We will specialise the discussion to weak lensing convergence κ in the following section, For the purpose of discussion here κ is a generic two-dimensional (2D) field defined over the celestial sphere. We will consider a flat patch of the sky for our discussion. The position vector θ is defined using the polar angle φ and Cartesian unit vectors i and j: Next we will consider two third-order statistics in projection.

Line Correlation Function in Projection
We will denote the LCF as L 2 (θ) which is defined through the following expression as a function of the angular scale θ is The LCF is a angle-averaged three-point collapsed correlation function where the three-points are in a collinear configuration. The two outer points are equidistant from the central point situated at θ 0 at a distance θ. A s represents the survey area. Assuming isotropy and homogeneity the estimator L(θ) depends only on the separation angular scale θ and does not depend on the position angle θ 0 . We consider the following convention for the Fourier transform ) map for the same source redshift. We use the publicly available maps discussed in [44].
The two-dimensional wave vector is denoted by l. Then The W (l) represents the smoothing window. We will not consider observational mask. The real-space statistics can be estimated simply by avoiding the masked region. The three-point correlation function of (θ) can be expressed as follows: The LCF in terms of the bispectrum can be expressed as follows: The three-point correlation function of the phase can be written in terms of the convergence bispectrum B κ (l 1 , l 2 , l 3 ) Using Eq.(3.7) in Eq.(3.3) we arrive at the following expression: We have used the following Bessel's first integral 11 (Eq.(71) of Mathsworld) to reduce Eq.(3.8a).
Here, J n denotes the Bessel functions of the first kind of order n.
In addition to the LCF of phases introduced in Eq.(3.8a)-Eq.(3.8b) we will also consider the following associated estimator: The derivation of Eq.(3.10b) follows the same steps as the derivation of Eq.(3.8b). The Effective Field Theory (EFT) provides a framework to extend the SPT results to smaller scales. Including the additional counter-terms from EFT in the expression of the kernel F 2 will extend the validity of results based on SPT to smaller angular scales not just for gravity induced non-Gaussianity but also for primordial non-Gaussianity [45].

Cumulant Correlators of Phases in Projection
The line correlation function considers a collapsed configuration of a triangle with equidistant points from the central location. The other collapsed configuration that is commonly used in 11 Bessel's first integral We use the fitting function proposed by [43] for the bispectrum. The simulation results show an ensemble average computed from 10 all-sky maps using publicly available software TreeCorr [46].
the literature was introduced in [47] in the context of density contrast in 3D and are known as the cumulant correlators. The two-to-one correlator defined below is of the lowest-order in the family of cumulant correlators and can be constructed cross-correlating a squared (θ) map with itself: Here S is the skew-spectrum of the (θ) map constructed from its bispectrum B Here µ denotes the cosine of the angle betwwn l 1 and l 2 i.e. µ = (l 1 · l 2 )/(l 1 l 2 ) with l i = |l i |.
Here the quantity in parentheses is the well-known Wigner-3j symbol. For more detailed derivation of flat-sky vs. all-sky correspondence see Appendix- §A. In Appendix- §B. we have presented generalisation to higher-order. The dashed-lines correspond to results computed using perturbation theory and the solid lines are obtained using a non-perturbative fitting function. We use the fitting function proposed by [43] for the bispectrum. The simulation results show an ensemble average computed using 10 all-sky maps. The three-point correlation function L 2 (θ) was computed from the phase maps using publicly available software TreeCorr [46].
The interest in cumulant correlators stems from the fact that they are two-point correlations but carry information about three-point statistics. They can be generalised easily to higher order. The information content is in general different for various triangular configurations. For an equilateral configuration see [42].

Results and Discussion
In this section we will discuss the results of our numerical investigations for various estimators and compare them against theoretical predictions.
Maps: To validate our analytical results we use the publicly available all-sky weak lensing maps generated by [44] 12 . The ray-tracing through N-body simulations were used to generate these maps. These simulation used 2048 3 particles to follow the evolution of gravitational clustering. To generate the convergence κ and the corresponding shear γ maps multiple lens planes were used; the source redshifts used were in the range z s = 0.05−5.30. We have chosen the maps with z s = 0.5, 1.0, 2.0 for our study. For CMB maps the lensing potentials were constructed using the deflection angles which were used to construct the lensing potentials and eventually the κ maps. These maps include post-Born corrections. Importance of post-Born terms in lensing statistics in the context of CMB studies were outlined [48]. However, recent studies have shown such corrections are not important in case of weak lensing statistics which 12 http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allsky raytracing/ We use the fitting function proposed by [43] for the bispectrum. We also show the results from perturbation theory (dashed-lines).
probe lower redshifts. The convergence maps were generated using an equal area pixelisation scheme in HEALPix 13 format [49].
The set of maps we use in this study are generated at N side = 4096 and were cross-checked against higher resolution maps constructed at a higher resolution N side = 8192, 16384 for consistency. These maps constructed at different resolution were found to be consistent with each other up to the angular harmonics ≤ 2000.
For our study, we have used high resolution maps N side = 4096. These maps were degraded to various lower resolution N side = 2048 In Figure 1 we show one such map for the source redshift z s = 0.5 (left panel) and the corresponding phase map (right panel). We analysed these maps using publicly available software TreeCorr 14 for the computation of two-and three-point correlation functions. The tso-point correlation functions are shown in Figure-2.
LCF for κ: The line correlation function Π 2 (θ) defined in Eq.(3.8b) is being plotted as a function of θ In Figure -3. From left to right we show results for source redshifts z s = 2.0, 1.0 and 0.5. We have presented the mean results estimated from ten all-sky maps. No noise was included in our study. Degraded maps with Nside = 2048 and 512 were respectively used to estimate the correlation function for the range θ s = 1 − 10 and 10 − 100 . The fitting function presented in [43] was used to compute the theoretical predictions. For the sake of completeness we have also shown the two-point correlation function of the κ maps in Figure  -2 for various redshifts, In this study we will see that at the low source redshift Post-Born correction plays a negligible role, but they play a significant role at higher redshift, e.g. in We use the fitting function proposed by [43] for the bispectrum (solid lines). The dots represent the cumulant correlators defined in Eq.(3.11) computed from one all-sky phase map by crosscorrelating the squared phase map 2 (θ) with itself (θ). We have used the publicly available software TreeCorr. Theoretical predictions correspond to the one computed using the fitting function proposed by [43] for the bispectrum. The deviations from the theoretical predictions at smaller angular scales is an effect of pixelisation.
case of lensing of CMB.
LCF for phase: In Figure -4 the LCF L ∈ (θ) is being plotted as a function of θ (in arcmin). The solid lines correspond to the ones defined in Eq.(3.10b) and the dashed-lines correspond to estimates from simulated all-sky weak lensing maps. Panels from left to right correspond to z s = 2.0, 1.0 and 0.5 respectively. The dashed-lines in each panel correspond to theoretical predictions computed using Born approximation. We have checked that the post-Born corrections do not make any appreciable difference. We have used ten realisations of all-sky maps to compute the numerical estimates. No noise was included. Degraded maps with Nside = 2048 and 512 were respectively used to estimate the correlation function for the range θ s = 1 − 10 and 10 − 100 . We use max = 2N side for our study. The theoretical results were computed using the fitting function presented in [43].
Skew-spectrum for phase: In addition to the LCF we have also introduced the skewspectrum estimator, not for κ(θ), but for the phase maps (θ) (denoted as S ) defined in Eq.(A.9b). For the bispectrum we have used the nonlinear fitting function [43] and for the perturbative calculations we have used the SPT result in Eq.(2.3b). We have presented three redshifts z s = 0.5, 1.0 and 2.0 in Fig.5 in different panels as indicated. The simulation results compare reasonably well against nonlinear predictions. We also see departure at high-which is a result of pixelisation. Our simulation results are derived using N side = 2048 which are generated by degrading maps originally created at N side = 4096. For PT results to be valid the maps need to be smoothed. We have also included z s = 1100.0 in our analysis which is relevant for CMB studies. The results are plotted in Fig.5 and for CMB lensing in Fig.7. We found that even with a single all-sky realisation we can estimate the S with a very high degree of accuracy. We have included realistic noise though it is expected from our previous study [17] inclusion of noise will not change our findings completely.
The importance of Post-Born (PB) corrections for CMB lensing lensing, was underlined in many recent studies, e.g., Ref. [17]. While such corrections do not make any significant contribution they do play important role for high redshift CMB studies. In Fig.7 we have shown the nonlinear results with and without PB corrections. Inclusion of Post-Born corrections is important to reproduce the simulation results. The maps used were of resolution N side = 2048. Higher-order generalisations of L 2 is presented in §B.
Cumulant correlators for phase: The cumulant correlators (CCs) carry equivalent information compared to the skew-spectrum S as they can be constructed from skew-spectrum. The theoretical cumulant correlator C 21 (θ) using Eq. (3.11). The numerical results were computed using the TreeCorr. The results are shown in Fig.6. This alows computation of CCs without broad binning as was the case for LCF. The deviation from theoretical prediction in the low-θ regime is related to the pixelisation effect. A comparison with the results presented for L 2 in Fig.4 confirms very high S/N in C 21 (θ). The generalisation to higher-order is straightforward and can by cross-correlating p-th power of (θ),i.e., p (θ) against its q-th power, i.e., q (θ 1 ) q (θ 2 ) .

Conclusions and Future Prospects
The primary aim of this paper was to introduce statistics of Fourier phases in the context of weak lensing studies. We have introduced the three-point phase correlation LCF to probe the non-Gaussianity in weak lensing convergence maps used in the study of galaxy clustering L 2 . In addition, an associated three-point statistics Π 2 was also introduced for probing statistics of κ maps. We have used a set of state-of-the-art all-sky simulations of weak lensing convergence maps to test our theoretical results as a function of source redshift. We have generalised both Π 2 (θ) and L 2 (θ) to higher-order. Next, we have adopted the cumulant correlator typically used for 3D density field and 2D convergence maps for statistics phases. We showed that available theoretical models can reproduce the numerical results with reasonable accuracy and our results can be used to select the range to retain in order to maintain a given level of theoretical accuracy. While we have focussed on weak lensing convergence, the statistical estimators presented here will be useful in other areas cosmology, e.g., galaxy clustering (in real and redshift space) and in clustering of Lyman-α absorbers.
Several extensions are possible on the theoretical front. The Effective Field Theory (EFT) provides a framework to extend the validility domain of the standard perturbation theory (SPT). A formulation of the phase statistics in EFT will improve its domain of validity. It will also be interesting to formulate the phase-statistics directly for shear. Theoretical modelling of bispectrum generated by intrinsic alignment [50] will also be useful. The perturbative regime is not particularly affected by gas physics. Nevertheless, k-cut filtering can be included [51] to filter out particularly sensitive modes.
Finally, our study was performed in a rather idealized observational setting to establish a baseline. Realistically complex follow-up studies will be presented in future. In particular a Fisher-based analysis independently and jointly with power spectrum will be presented elsewhere.
implementation of the software TreeCorr. We would like to thank Alexander Eggemeier for many useful discussions.

A All-sky Expressions
Our primary aim in this appendix is to develop all-sky estimators presented in the text of the paper. Following [52][53][54] the spherical polar co-ordinate (θ, φ), radial co-ordinate r = 2 sin θ/2 ≈ θ Equivalently in the harmonic domain |l| = and φ l denotes the polar angle: The flat-sky κ(l) and its all-sky harmonic counterpart κ m are related by the following expression: The spherical harmonic basis Y m (θ, φ) and rhe flat sky basis are related by the following expressions: If we use the following expansion of the plane wave: Next, we consider the case of bispectrum and related estimator: Here the quantity in parentheses is the well-known Wigner-3j symbol which enforces the rotational invariance. It is only non-zero for the triplets ( 1 , 2 , 3 ) that satisfy the triangular condition and 1 + 2 + 3 is even. The reduced bispectrum b κ 1 2 3 is useful in directly linking the all-sky bispectrum and its flat-sky counterpart. For the convergence field κ, b κ 1 2 3 is defined through the following expression: The flat-sky bispectrum is similarly defined through: The flat-sky bispectrum B κ (l 1 , l 2 , l 3 ) is identical to the reduced bispectrum b κ 1 2 2 for high multipole [56]. This can be shown by using the following asymptotic relationship: Here the β functions represent the smoothing window functions and µ is the cosine of the angle between vector l 1 and l 2 i.e. µ = (l 1 · l 2 )/(l 1 l 2 ), where |l i | = l i . The skewness in terms of these expressions is given by:

B Higher-order LCFs
In this section we consider generalisation of the third-order LCFs to higher-order. The family of LCFs can be seen as a natural extension of two-point cumulant correlators often used in the literature. LCFs can be thought of as three-point generalisation of two-point cumulant correlators. The following statistics are relevant for the kappa field as well as the for the phases. In this section we generalise the line-correlation to higher order. The fourth-and fifth-order generalisations are as follows.