{'fraction_non_alphanumeric': 0.03706391999631319, 'fraction_numerical': 0.02018526199364026, 'mean_word_length': 4.411284289276808, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 33, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Information content of weak lensing power spectrum and bispectrum: including the non-Gaussian error covariance matrix

We address the amount of information in the non-Gaussian regime of weak lensing surveys by modelling all relevant covariances of the power spectra and bispectra, using 1000 ray-tracing simulation realizations for a Lambda-CDM model and an analytical halo model. We develop a formalism to describe the covariance matrices of power spectra and bispectra of all triangle configurations. In addition to the known contributions which extend up to six-point correlation functions, we propose a new contribution `the halo sample variance (HSV)' arising from the coupling of the lensing Fourier modes with large-scale mass fluctuations on scales comparable with the survey region via halo bias theory. We show that the model predictions are in good agreement with the simulation once we take the HSV into account. The HSV gives a dominant contribution to the covariance matrices at multipoles l>10^3, which arises from massive haloes with a mass of>10^14 solar mass and at relatively low redshifts z<0.4. Since such haloes are easily identified from a multi-colour imaging survey, the effect can be estimated from the data. By adding the bispectrum to the power spectrum, the total information content or the cumulative signal-to-noise ratio up to a certain maximum multipole of a few 10^3 is improved by 20--50 per cent, which is equivalent to a factor of 1.4--2.3 larger survey area for the power spectrum measurement alone. However, it is still smaller than the case of a Gaussian field by a factor of 3 mostly due to the HSV. Thus bispectrum measurements are useful for cosmology, but using information from upcoming surveys requires that non-Gaussian covariances are carefully estimated.

Sky Survey (LSST 6 ), the European Space Agency (ESA) Euclid satellite mission 7 , and the NASA Wide-Field Infrared Survey Telescope (WFIRST) satellite mission. 8 Among different cosmological probes, weak gravitational lensing or cosmic shear is recognized as one of the most promising methods for constraining the nature of dark energy, provided systematic errors are well under control (see Bartelmann & Schneider 2001;Schneider 2006;Hoekstra & Jain 2008, for reviews). The bending of light rays emitted from a distant galaxy due to the foreground mass distribution causes the image to be distorted. The distortion signal is too weak for us to measure in single galaxies, but we can use a sufficiently large number of galaxy images, available from wide-field survey, to detect the correlated shear signals existing in-between different galaxy images. Weak lensing is a unique method of measuring the total matter distribution including dark matter, free of galaxy bias uncertainty, and allows a direct comparison of the measurement with theory that is in most case about the statistical properties of the dark matter distribution. The theoretical predictions are obtained using N-body simulations (e.g. Springel et al. 2006) and/or analytical approaches (e.g. Bernardeau et al. 2002;Cooray & Sheth 2002). The cosmological weak lensing signal has been measured by several groups (e.g. Hamana et al. 2003;Hoekstra & Jain 2008;Schrabback et al. 2010, and also see references therein), and we are waiting for measurements with much higher statistical precision from upcoming surveys.
Most previous works on weak lensing, in theoretical and observational studies, have focused on the shear two-point correlation function or equivalently its Fourier transform, the power spectrum, as the statistics to quantify the lensing field. Although these statistics contain the full information when the field is a Gaussian random field as in the cosmic microwave background (CMB) field (Komatsu et al. 2011), it is not the case for the lensing field because non-linear clustering in structure formation causes a coupling between different Fourier modes, and the mass density field at redshifts relevant for lensing surveys is not Gaussian. In fact, various studies have shown that the information content carried by the power spectrum might be saturated at multipole scales of a few 10 3 [see Hamilton et al. (2006); Takahashi et al. (2009);Neyrinck et al. (2009) for the 3D mass density field, Sato et al. (2009) and Seo et al. (2011) for the lensing field, and Lee & Pen (2008) for the result from the actual data]. In particular, Sato et al. (2009) used 1000 ray-tracing simulation realizations for a Λ cold dark matter (ΛCDM) model to study the power spectrum covariance and the information content of the power spectrum. They found that the information content is reduced by a factor of 2 at multipoles l ≃ 10 3 compared to the Gaussian case for a survey with typical source redshift of zs ≃ 1. Further, they showed that large-scale mass density fluctuations of scales outside the simulation area contribute significantly to non-Gaussian terms of the covariance. They developed a formalism to describe the new non-Gaussian contribution by the number fluctuations of massive haloes based on halo bias theory, which we hereafter call the halo sample variance (HSV; also see Hu & Kravtsov 2003;Takada & Bridle 2007).
Some fundamental questions remain unresolved: how important and useful are the non-Gaussian signals in the lensing field for cosmology? Which statistical method to extract the non-Gaussian signals is most useful? Can we recover the Gaussian information content, which should have existed in the linear field or the primordial field, by combining the power spectrum and the non-Gaussian signals? For weak lensing, there is additional expectation that the non-Gaussian signals will be useful for cosmology, because the skewness, for example, has been shown complementary to the power spectrum in its dependence on cosmological parameters (Bernardeau et al. 1997;Jain & Seljak 1997;Hui 1999;Jain et al. 2000;White & Hu 2000;Hamana & Mellier 2001;Van Waerbeke et al. 2001;Cooray & Hu 2001b;Takada & Jain 2002, 2004Dodelson & Zhang 2005;Kilbinger & Schneider 2005;Semboloni et al. 2008;Bergé et al. 2010;Munshi et al. 2011;Pires et al. 2012). The attempt to measure the non-Gaussian signals from actual data was also made by several groups (Bernardeau et al. 2002;Zhang et al. 2003;Jarvis et al. 2004), and the first significant detection was recently reported by Semboloni et al. (2011), showing an improvement in cosmological parameters compared to the two-point statistics alone.
In this paper, we study the lensing bispectrum, which contains the lowest-order non-Gaussianity of the weak lensing field and is a natural extension of the power spectrum. We consider all triangle configurations available from a given range of multipoles and their full covariance matrix including the non-Gaussian contributions up to six-point correlation functions as well as the HSV term, while only the Gaussian errors have been assumed in most previous work (Takada & Jain 2004;Martin et al. 2012). We use the 1000 simulation realizations to study the usefulness and complementarity of the lensing bispectrum compared to the power spectrum, and also develop an analytical formula to describe the bispectrum covariance for a given cosmology. In particular, we will show that the HSV gives a significant contribution to the bispectrum covariance at l a few 10 2 , and that the bispectrum does carry additional information to the power spectrum even in the presence of these significant correlations. Thus we will give a quantitative answer to the fundamental questions above. This paper is organized as follows. After briefly reviewing the lensing power spectrum and bispectrum in Section 2, we develop a formulation to describe the bispectrum covariance in Section 3. In Section 4, we show the main results: we study the bispectrum covariance using both simulations and analytical model predictions. We quantify the information content of the lensing bispectrum by including contributions from all triangle configurations. Section 5 is devoted to discussion and conclusion.

PRELIMINARIES: LENSING POWER SPECTRUM AND BISPECTRUM
In the context of cosmological gravitational lensing, the convergence field is expressed as the weighted projection of the three-dimensional density fluctuation field between source and observer (see Bartelmann & Schneider 2001;Schneider 2006, for a thorough review) where θ is the angular position on the sky, χ is the comoving distance, and χH is the distance to the Hubble horizon. We assume a flat geometry throughout this paper, and the radial distance χ is equivalent to the comoving angular diameter distance. The comoving distance χ(a) from an observer at a = 1 to a source at a is expressed in terms of the Hubble expansion rate H(a) as χ(a) = 1 a da ′ /[H(a ′ )a ′2 ]. For source galaxies at a single redshift, the lensing efficiency function WGL(χ) is defined as where χs is the distance to the source galaxies. See equation (4) in Takada & Jain (2009) for the lensing efficiency functions for tomographic redshift bins. Under the flat-sky approximation, the Fourier transform of the lensing field is defined as When the sky coverage of a survey is finite, we need to use the discrete Fourier decomposition, rather than the infinite-range Fourier decomposition (also see appendix in Takada & Bridle 2007, for details).
In this paper, we study the bispectrum of lensing field and the covariance of the bispectrum. The n-point power spectra relevant for the bispectrum covariance are defined in terms of the ensemble averages of the convergence fields in Fourier space as κ l 1κ l 2κ l 3 ≡ (2π) 2 B(l1, l2, l3)δ 2 D (l1 + l2 + l3), κ l 1κ l 2κ l 3κ l 4 c ≡ (2π) 2 T (l1, l2, l3, l4)δ 2 D (l1 + l2 + l3), κ l 1κ l 2 . . .κ ln c ≡ (2π) 2 Pn(l1, l2, . . . , ln)δ 2 D (l1 + l2 + · · · + ln) if n 5, where δ 2 D (l) is the Dirac delta function; P , B and T are the lensing power spectrum, bispectrum and trispectrum, respectively; Pn is the n-point power spectrum. For the bispectrum covariance, we need to include up to the six-point power spectra P6. The higher-order correlation function than the bispectrum is the connected part of the n-point function, which characterizes the non-Gaussianity of the lensing field and cannot be expressed in terms of products of the power spectrum (or any other lower-order correlation functions). Exchange symmetry of wavevectors, li ↔ lj, reflects that the n-point correlation functions are invariant under permutations of the arguments. The delta functions come from the parallel translation invariance for a statistically homogeneous field as is the convergence field.
The lensing power spectrum and bispectrum can be given as the weighted line-of-sight projection of the three-dimensional power spectrum and bispectrum of the underlying mass distribution. Employing Limber's approximation (Limber 1954) and the flat-sky approximation, the lensing power spectrum and bispectrum are expressed as where P δ and B δ are the power spectrum and bispectrum of the mass distribution at each redshift χ(= χ(z)). Thus once the n-point spectra of the mass density field are given for a given cosmological model, we can compute the n-point spectra of the lensing field. The above equations also mean that statistical properties of the lensing field arise from those of the mass density field, since the prefactors such as the WGL(χ) are pure geometrical quantities, not statistical variables. The power spectrum measurement for an actual survey is affected by the intrinsic shape noise due to a finite sampling of source galaxy shapes: where σǫ is the rms of intrinsic ellipticities per component, andng is the mean number density of source galaxies per unit steradian. In the following, we will often omit the notation obs to refer P obs (l) for notational simplicity. Throughout this paper we assume that the orientation of intrinsic galaxy shape is random and the shapes of different galaxies are uncorrelated; the shape noise is a Gaussian random field. The bispectrum is a measure of the non-Gaussianity, so is not affected by the shape noise.

LENSING COVARIANCE
The covariances of the lensing spectrum and bispectrum describe a measurement accuracy of the spectra for a given survey. There are several sources of the measurement errors: the shot noise arising due to a finite sampling of galaxy shapes and the sample variance arising due to a finite survey area. If the lensing field is Gaussian, the different Fourier modes with l = l ′ are independent, and therefore the sample variance is determined by the number of independent Fourier modes for a given multipole bin l that are available from the survey, yielding a simple formula of the sample variance contribution (e.g. Knox 1995). However, this is not the case for the lensing field, because the lensing field is highly non-Gaussian at scales of interest (Takada & Jain 2004, 2009Sato et al. 2009), and the different Fourier modes correlate with each other. In the following, we discuss theory for the lensing covariance matrices.

Power spectrum covariance
The power spectrum covariance has been well studied by previous works Cooray & Hu 2001a;Takada & Bridle 2007;Takada & Jain 2009;Sato et al. 2009). In particular, Sato et al. (2009) derived an expression of the power spectrum covariance including a new contribution from the mass density fluctuations at larger scales than the survey area, and showed that the formula well reproduces their ray-tracing simulation results. According to this work, the power spectrum covariance is given as where δ K l i l j is the Kronecker delta, δ K l i l j = 1 if li = lj within the bin width and otherwise δ K l i l j = 0; Ωs is the survey area in units of steradian; A(li) is the area of the above integration in Fourier space, given as As(li) ≡ |l|∈li d 2 l, where the integration range is confined to the wavevectors satisfying the condition li − ∆l/2 |l| li + ∆l/2 (∆l is the bin width around the ith bin, li); the third term CovHSV is the new contribution which we call the halo sample variance (HSV) contribution (see below). The quantity Npairs(li) is the number of independent pairs of two vectors l and −l in Fourier space, where the vector l has the length li within the bin width and "independent" here means different pairs discriminated by the fundamental Fourier mode of a given survey, l f ≃ 2π/Θs (Θs is the angular scale of the survey area). At the limit li ≫ l f , A(li) ≃ 2πli∆l and the number of independent Fourier modes is given as where f sky is the sky fraction defined as f sky ≡ Ωs/4π. See Takada & Bridle (2007) for a pedagogical derivation of the power spectrum covariance based on the discrete Fourier decomposition formulation (except for the third term Cov P P HSV ). In equation (11), we ignored effects of non-trivial survey geometry for simplicity.
The first and second terms on the r.h.s. of equation (11) are the standard covariance terms studied in most previous works. The first term describes the Gaussian covariance term that vanishes when li = lj , i.e. no correlation between different multipole bins. The second term gives a non-Gaussian term arising from the lensing trispectrum (four-point correlation function), which describes the mode coupling between different multipole bins. Both terms scale with survey area as 1/Ωs; the amplitudes decrease with increasing the survey area. It should also be noted that the Gaussian term depends on the multipole bin width, while the non-Gaussian terms do not; a larger bin width relatively reduces the Gaussian term contribution at the multipole bin.
The third term of equation (11) arises from the mode coupling of the Fourier mode of our interest, li, with large-scale modes of scales comparable with or even outside the survey region. Such large-scale modes cannot be seen by an observer, but affect the power spectrum estimation. If the entire survey region happens to be in an overdense/underdense region, which is caused by the large-scale mass density fluctuations, the number of massive haloes becomes larger/smaller than the ensemble average according to the halo bias theory (Mo & White 1996;Sheth et al. 2001). Thus the number of massive haloes found in a finite survey volume correlates with the mass density fluctuations of scales comparable with or larger than the survey field (Hu & Kravtsov 2003). To be more explicit, the number fluctuations of haloes in a mass M and in the redshift slice centred at z are given as where d 2 V /dzdΩ is the comoving volume per unit redshift interval and per unit solid angle, d 2 V /dzdΩ = χ 2 for a flat universe; dn/dM is the ensemble-averaged mass function of haloes in the mass range [M, M + dM ]; b(M ) is the halo bias parameter; δm(Θs) is the mass density fluctuation averaged within the survey volume in the redshift slice which has area Ωs and the redshift width ∆z. The lensing power spectrum amplitudes at small angular scales are sensitive to the number of massive haloes in the survey region and then correlates with the number fluctuations, which results in the HSV. Note that the ensemble average of the power spectrum is not affected by the large-scale mode due to the fact δ m(Θs) = 0. At the limit of li, lj ≫ 1, the HSV contribution is given as whereW (l) is the Fourier transform of the survey window function; P L m is the linear mass power spectrum;κ l (χ) is the Fourier transform of the convergence field for which we assume a Navarro-Frenk-White (NFW) halo (Navarro et al. 1997) (see equation 28 in Oguri & Takada 2011 or section 3.2 in Takada & Jain 2003b for the expression ofκ l ). For notational simplicity, we omit to denote the redshift dependence of dn/dM , b(M ) andκ l . In this paper, we simply consider the window function given byW (x) = 2J1(x)/x, which corresponds to a circle-shaped survey geometry with a radius of Θs. For the halo model ingredients [dn/dM, b(M ) and the NFW profile], we will throughout this paper employ the same models as used in Takada & Jain (2009). Roughly speaking, the HSV term can be expressed as Cov P P HSV ∼ P 1h (li)P 1h (lj )σ 2 m (Ωs), where P 1h is the one-halo term of the lensing power spectrum and σm(Ωs) is the rms of the projected linear mass density fluctuations smoothed with the angular scale of the survey area. As implied from the above equation, the HSV affects the band powers of different multipoles in the same way, and does not change the shape of the power spectrum. This HSV term cannot be realized as long as the discrete Fourier decomposition is used for deriving the power spectrum covariance, because the large-scale modes outside the survey region cannot be described by Fourier modes confined inside the survey region.
Another important feature is that the HSV contribution depends on the survey area via the integration of the linear mass power spectrum, kdkP L m (k) W (kχΘs) 2 unlike the other terms which scale as 1/Ωs. For a power-law linear power spectrum, P L m (k) ∝ k n , the HSV term is found to scale as Cov P P HSV ∝ 1/(Ωs) 1+n/2 . Hence, for n < 0, which is indeed the case at k keq (keq is the wavenumber of the matter-radiation equality), the HSV amplitude decreases more slowly with increasing area coverage than other terms. On the other hand, when n > 0 or k keq, the HSV term decreases more quickly. Thus the HSV term depends on the survey area in a non-trivial way, and we will study the relative importance of the HSV term for different survey areas.
There are even other sources of non-Gaussian errors arising from a correlation of the lensing field in the weakly non-linear regime with the mass density fluctuations of scales comparable with or larger than the survey area. This contribution can be formulated based on the perturbation theory of structure formation, which is valid for the mass density field in the weakly non-linear regime. Rimes & Hamilton (2005) first studied this effect for the 3D mass power spectrum, and named this new contribution the beat-coupling mode (also see Hamilton et al. 2006;Sefusatti et al. 2006;Takahashi et al. 2009). Furthermore, de Putter et al. (2012 recently pointed out that the large-scale density fluctuations cause an apparent modulation in the mean density estimated from a finite survey region, and add an additional negative contribution to the covariance. This term was shown to have a similar amplitude to the beat-coupling mode. These large-scale mode contribution, which is relevant for the weakly non-linear regime, differs from the HSV effect, and Takada & Jain (2009) showed that the contribution to the lensing power spectrum covariance is negligible compared to the non-Gaussian errors arising from the four-point function (the second term in equation 11) at multipoles of l a few 10 2 . Hence, in this paper, we ignore the non-Gaussian errors arising from the mode-coupling in the weakly non-linear regime for simplicity.

Bispectrum covariance
As we discussed above, the lensing bispectrum is given as a function of triangle configurations. An estimator of the lensing bispectrum from the finite-area lensing survey can be found, by extending the method developed in Takada & Bridle (2007) (see Appendix A for the details), aŝ B(l1, l2, l3) = 1 ΩsNtrip(l1, l2, l3) q iκ q 1κ q 2κ q 3 ∆q 123 (l1, l2, l3), where q 123 ≡ q 1 + q 2 + q 3 and the summation runs over all the pixels of q 1 , q 2 and q 3 . The function ∆q 123 denotes the selection function which is unity if each vector has a target length of li − ∆li/2 qi li + ∆li/2 (i = 1, 2, 3) and the three vectors form the triangle configuration in Fourier space, q 123 = 0; otherwise ∆q 123 (l1, l2, l3) = 0. The prefactor 1/Ωs is from our definition of the discrete Fourier decomposition (see Takada & Bridle 2007). The quantity Ntrip is the number of the triplets in Fourier space that form a given triangle configuration specified by three side lengths l1, l2, l3 with the bin widths. This is calculated from the selection function as Ntrip = q i ∆q 123 (l1, l2, l3). As described in Joachimi et al. (2009) and Appendix A, for the limit of large multipole bins, l1, l2, l3 ≫ l f , we can analytically estimate Ntrip as where ∆li is the bin width of the ith side length. For the small-angle scales li ≫ 1 (flat-sky approximation limit) we are interested in, this equation gives a good approximation to the Wigner-3j symbols which appears in the bispectrum covariance derived under the full-sky approach (see equation 16 in Takada & Jain 2004). The bispectrum covariance can be similarly defined as The bispectrum covariance arises from the six-point correlation function of the lensing field. Figure 1. Illustration of the different terms of the bispectrum covariance. The triangle configurations for two bispectra in the covariance matrix, B(l 1 , l 2 , l 3 ) and B(l ′ 1 , l ′ 2 , l ′ 3 ), are specified by sets of the three vectors (l 1 , l 2 , l 3 ) or (l ′ 1 , l ′ 2 , l ′ 3 ), denoted by the solid and dashed lines, respectively. The three vectors satisfy the triangle conditions l 1 + l 2 + l 3 = 0 and l ′ 1 + l ′ 2 + l ′ 3 = 0. (i) The Gaussian part of the bispectrum covariance, which arises only if the two triangle configurations have the same shape (within the coarseness of the bin widths). The two vectors with the same length but in the opposite direction such as l 1 and l ′ 1 yield the power spectrum after the ensemble average. Hence the Gaussian term amplitude is proportional to P (l 1 )P (l 2 )P (l 3 ) and contributes to the diagonal terms of the covariance matrix. (ii) Non-Gaussian part of the bispectrum covariance which arises if one side of the two triangles is in the same length and direction; here l 1 = l ′ 1 is shown as an example. The amplitude of this covariance term is proportional to B(l 1 , l 2 , l 3 )B(l ′ 1 , l ′ 2 , l ′ 3 ). (iii) Non-Gaussian part that arises if one side length of the two triangles is the same but in the opposite direction, l 1 = −l ′ 1 . The amplitude is O(P T ); the figure shows a case of P (l 1 )T (l 1 , l 3 , l ′ 2 , l ′ 3 ). The surrounding four vectors (l 2 , l 3 , l ′ 2 , l ′ 3 ) form a quadrangular configuration satisfying the condition l 2 + l 3 + l ′ 2 + l ′ 3 = 0, which gives the trispectrum contribution. (iv) Non-Gaussian part which arises for generic triangle configurations and therefore contributes to all diagonal and off-diagonal terms of the covariance matrix. The amplitude is proportional to the connected part of the six-point correlation function, P 6 . As indicated, the 6 vectors of (l 1 , l 2 , l 3 , l ′ 1 , l ′ 2 , l ′ 3 ) arise from the two triangles that form the six-point configuration in Fourier space, although the two vertices of the 6 points are collapsed to one point due to the triangle conditions. We present the detailed derivation of the bispectrum covariance in Appendix A based on the discrete Fourier decomposition formulation. Here we just give the expression of the bispectrum covariance, which has three contributions, the Gaussian and non-Gaussian errors and the HSV contribution: where the notation 'NG' stands for the 'non-Gaussian' error contribution, and P6 denotes the connected part of the six-point correlation function. Fig. 1 shows a diagram picture of these covariance terms from the first to the fourth lines on the r.h.s. of equation (18). When further including the intrinsic shape noise contribution, we just replace the power spectra in the terms of the above equation (O(P 3 ) and O(P T ) terms) with the power spectrum including the shot noise contribution (equation 10).
The terms of the first line on the r.h.s. are the Gaussian covariance terms, which contribute only to the diagonal terms of the bispectrum covariances. The combination of the Kronecker deltas δ K is non-vanishing only if the two triangle configurations are in the same "shape" within the bin widths. In particular, if triangle configurations have symmetry such as isosceles or equilateral triangles, the combination of Kronecker deltas (the terms in the square bracket on the first line) yield a factor of 2 or 6 for isosceles and equilateral triangles, respectively. The factors account for the fact that different triangles transformed by parity and permutation transformations (li ↔ l ′ j ) are not independent for a statistically homogeneous and isotropic field. For a general triangle configuration li = lj , the factor becomes unity. The prefactor Ntrip(l1, l2, l3) is given by equation (16). The survey area dependence of the HSV contribution for the power spectrum and the bispectrum, relative to that of the Gaussian error, which scales as Cov Gauss ∝ 1/Ωs. Here we consider the scale of l = 2000 for P (l) and Beq(l), employ the flat sky approximation, and assume the survey geometry given by Ωs = πΘ 2 s .
The terms from the second to the fourth lines are the non-Gaussian error contributions, which arise from the higher-order correlation functions of the lensing field. The coefficient of each term such as 2π/(Ωsl1∆l1) is from the number of independent configurations in Fourier space that form a given configuration of six wavevectors (l1, l2, l3, l ′ 1 , l ′ 2 , l ′ 3 ) in Fig. 1 (also see Appendix A for the mathematical derivation). The angular integration in the fourth term including P6 is over the angle ψ between the vectors l1 and l ′ 1 in order to include contributions over all the possible six-point configurations in Fourier space. Note that the terms in the first, second and third lines on the r.h.s. depend on the multipole bin widths such as ∆l1, while the term including P6 and the HSV do not depend on the bin widths.
As in the power spectrum covariance (equation 11), the number fluctuations of massive haloes due to the large-scale mass fluctuations affect the bispectrum estimated from a finite area survey. The HSV contribution to the bispectrum covariance is given as This contribution has not been realized in previous works. By comparing with ray-tracing simulations, we will show below that the HSV contribution is dominant over other covariance terms at l 1000, in the non-linear regime, and adding the HSV term to the model predictions significantly improves agreement with the simulation results. Again note that the HSV term scales with the survey area in a non-trivial way via the linear mass power spectrum (see the discussion below equation 14), while the other terms scale with survey area as 1/Ωs. The left-hand and middle panels of Fig. 2 show which haloes of mass and redshift range contribute to the HSV term of the bispectrum covariance for an equilateral triangle configuration of a given side length. For the triangle with l = 1000, haloes with masses 10 14 M⊙ and at redshift z 0.4 give a dominant contribution to the HSV effect. These haloes are relatively easy to be identified from a concentration of galaxies on the sky, X-ray or the Sunyaev-Zel'dovich effect. In other words, identifying such massive haloes in a survey region and comparing the number with the expected number for a fiducial cosmological model will help understand the HSV effect for the weak lensing observables in the survey region (also see Takada & Bridle 2007, for the similar discussion). The right-hand panel shows how the HSV terms scale with the survey area, in comparison with the Gaussian covariance term for the power spectrum covariance and the bispectrum covariance of equilateral triangles. Here we consider the particular multipole bin l = 2000, and assume the flat-sky approximation and a circular survey geometry, Ωs = πΘ 2 s , for simplicity. The plot shows that the HSV term decreases more slowly than other terms as Ωs increases up to a few 100 deg 2 .

Cross-covariance between power spectrum and bispectrum
Since the power spectrum and bispectrum are not independent for the non-Gaussian field, we need to account for the cross-covariance between the two observables in order not to double-count the information content. We can derive the cross-covariance (see Appendix B) similarly to the power spectrum and the bispectrum covariance as where P5 is the five-point correlation function, and ψ is the angle between the vectors l and l1 as in equation (18). The HSV term is given as We will use these equations when computing the total information content for a combined measurement of the lensing power spectrum and bispectrum for a given survey.

Halo sample variance contribution to the n-point correlation function measurement
As we have seen for the HSV contributions to the lensing spectrum covariances (equations 14, 19 and 21), the equations have similar forms at the small-angle limit, where the one-halo term is dominated in the halo model picture. Extending these findings, we can find the HSV contribution to the covariance matrix between any nand n ′ -point correlation functions in Fourier space: Roughly speaking, the amplitude of the HSV term simply scales as CovHSV ∼ P 1h n (li)P 1h n ′ (l ′ j )σ 2 m (Θs). Thus any n-point correlation functions at small angle scales can be affected by the large-scale mass fluctuations of scales comparable with or outside the survey area.
We should also emphasize that the HSV contribution affects measurements of any two-or three-dimensional correlation functions from a finite area survey, and can be very important if one is interested in the small-scale signals which are sensitive to the abundance of haloes in the finite survey region (e.g. see Shaw et al. 2009;Zhang & Sheth 2007, for a similar discussion on the SZ power spectrum measurement).

Halo model predictions for the lensing covariances
As we have described up to the preceding section, the power spectrum and bispectrum covariance calculations require to compute the four-, five-and six-point correlation functions in addition to the power spectrum and bispectrum. For the power spectrum and bispectrum, some theoretical models are proposed by comparing with simulations, e.g. Smith et al. (2003) and Valageas & Nishimichi (2011a) for the power spectrum and Scoccimarro & Frieman (1999) and Valageas & Nishimichi (2011b) for the bispectrum. The higher-order functions are, however, fairly uncertain because there are fewer studies to compare the model predictions with simulations (see e.g. Takada & Jain 2002, for an attempt to compute the kurtosis which is the collapsed four-point function), partly because the higher-order correlations require a substantial amount of computational costs. Instead of pursuing a reliable model for the higher-order correlation functions, in this paper we employ the halo model approach to compute the higher-order functions (Peacock & Smith 2000;Seljak 2000;Ma & Fry 2000;Scoccimarro et al. 2001;Cooray & Sheth 2002) in which the correlations of the mass distribution are expressed as two separate contributions: correlations of dark matter particles within the same halo and correlations between particles in different haloes. We have found that, up to the four-point correlation functions, the halo model predictions are accurate at 10-30 per cent level in the amplitude compared to N-body simulations (in particular for lensing fields; Takada & Jain 2002, 2003a. Since our purpose of this paper is to assess the importance of non-Gaussian error contributions to the bispectrum covariance matrices, we consider that the halo model approach is adequate enough. We know that most of the lensing information comes from small angle scales in the non-linear clustering regime to which the onehalo term, the correlation arising from the same halo, provides a dominant contribution. In addition, the non-Gaussian errors are important only at the small scales, as can be explicitly found from fig. 5 in Sato et al. (2009). For these reasons, we include only the one-halo terms to compute the non-Gaussian error contributions to the lensing covariances, which significantly simplifies the computation. Although the n-point correlation function depends on n wavevectors such as Pn(l1, l2, . . . , ln), the one-halo term does not depend on any angle between the vectors, but rather depends only on the length of each vector; P 1h n (l1, l2, . . . , ln) = P 1h n (l1, l2, . . . , ln), reflecting spherical mass distribution around halo in a statistical average sense. To be more explicit, assuming the Limber's approximation, the one-halo term of the n-point correlation function can be computed as follows (see around equation 30 in Takada & Jain 2003a): The one-halo term is computed by a two-dimensional integration, as the Fourier transform of the lensing field due to an NFW halo,κM (l), can be computed analytically for a given halo with mass M and at redshift z. As can be found from equations (11), (18) and (20), some of the non-Gaussian terms can be further simplified; e.g. one of the non-Gaussian terms in the power spectrum covariance (equation 11) can be simplified as 1 Ωs |l|∈li where we have assumed that the lensing trispectrum does not largely change within the multipole bin ∆l around the bins li and lj. Thus the above approach allows a faster computation of the power spectrum and bispectrum covariances.

Ray-tracing simulations
To study the lensing covariance matrices, we use 1000 realizations of ray-tracing simulations for a ΛCDM model in Sato et al. (2009). Although the simulations were done for various source redshifts ranging from zs = 0.6 to 3, we use the outputs of zs = 1 in this paper. In brief, each realization has an area of 5 × 5 deg 2 (Ωs = 0.0076 sr) in square shaped geometry. The ΛCDM model adopted is characterized by cosmological parameters: the matter density parameter Ωm = 0.238, the baryon density parameter Ω b = 0.042, the initial spectral index ns = 0.958, the amplitude of the density fluctuations σ8 = 0.76 and the Hubble constant of H0 = 100h km s −1 Mpc −1 with h = 0.732. The linear matter power spectrum used to set the initial conditions of N-body simulations is computed by the public code CAMB (Lewis et al. 2000). It was shown that the ray-tracing simulations are reliable to within a 5 per cent accuracy up to multipole l ≃ 6000 for the power spectrum and up to l ≃ 4000 for the bispectrum. See Sato et al. (2009) and Valageas et al. (2012) for more details. 9 As can be found from fig. 1 in Sato et al. (2009), the ray-tracing simulations were done in a light cone of area 5×5 deg 2 , viewed from an observer position (z = 0). The projected mass density fields in intermediate-redshift slices were generated from N-body simulations which have a larger simulation box than the volume covered by the light cone. Hence the lensing fields have contributions from the mass density field of scales outside the ray-tracing simulation area, although, exactly speaking, the modes outside the N-body simulation box were not included (see Appendix C for the effect). Thus the ray-tracing simulations allow us to study the HSV effect on the covariance matrices.

Measuring power spectrum, bispectrum and the covariance matrices from simulations
In each simulation realization, the lensing convergence field, κ(θ), is given on 2048 × 2048 grids. We used the FFT method to compute the Fourier transformed field,κ(l). The fundamental mode of the discrete Fourier decomposition is l f = 72(= 2π/5 • ) and the Nyquist frequency is ∼ 70000, which is large enough compared to the resolution limit of the N-body simulations. We use multipole bins that are logarithmically spaced by ∆ ln l = 0.3 (∆ log l ≃ 0.13), which significantly reduces the number of triangle configurations (the number of different bispectra) we need to consider, compared to ∆l = 1 as in the CMB case. We consider 16 multipole bins in total; the first bin is in the range l = [72, 97.2], and the 16th bin is in the range l = [6481.2, 8748.8], so lmin = 72 and lmax = 8748.8.
The power spectrum of the ith multipole bin li, P (li), is estimated from each realization by computing an azimuthal average of the estimator,κ lκ−l , where |l| resides in the target bin. We then averaged the estimated power spectra in 1000 realizations to estimate the ensemble-averaged power spectrum (corresponding to the power spectrum for the area of 1000 × 25 = 25000 deg 2 ). Then we computed the scatters among the power spectra of 1000 realizations in order to estimate the covariance matrix of the power spectra. Hence the covariance matrix is for an area of 25 deg 2 .
The bispectrum is given as a function of triangle configuration. We use three side lengths (l1, l2, l3) to parametrize triangle configuration, where the triangle conditions are given as |lj − l k | li lj + l k . Although the multipole bin used has a logarithmically spaced bin width, we impose the triangle conditions on the central values of the multipole bins. In addition we impose the condition l1 l2 l3 so that every triangle configuration is counted once. For the 16 multipole bins above, we have 204 triangle configurations in total.
For a triangle configuration that is specified by the side lengths (l1, l2, l3) (with the bin widths), we can estimate the bispectrum from each ray-tracing simulation by averaging the estimator ] over all the triplets (l ′ 1 , l ′ 2 , l ′ 3 ) which satisfy the triangle conditions; the length of each vector is in the triangle bin such as l1,min l ′ 1 l1,max. Note that, as long as the triplets of (l ′ 1 , l ′ 2 , l ′ 3 ) in Fourier space have the same side lengths l1, l2, l3 within the bin widths, all the triangles transformed by the rotation and parity transform are equivalent to yield the same bispectrum due to the rotation and parity invariance for a statistically homogeneous and isotropic field. In our simulations, the Fourier transformed convergence field,κ(l), is given on 2048 × 2048 grids in Fourier space, where the grids are linearly spaced by the fundamental mode, l f = 2π/Θs = 72. To have an efficient computation over 1000 realizations, we first built the table of three vectors (grids), and is assigned to one of the triangle configurations binned by three side lengths (l1, l2, l3). Then we used the same table of triplets for the 1000 realizations to compute the average and scatters of the estimated bispectra as a function of triangle configurations.
To be more precise, we built the table of three vectors (grids) in the way illustrated in Fig. 3. First, we choose the first vector l ′ 1 from one of the grids in the upper half of Fourier space, i.e. 0 arg(l ′ 1 ) < π, by imposing the condition that the length l ′ 1 is in the range of the multipole bin, lmin l ′ 1 < lmax. Then we survey for the second vector l ′ 2 from the region where the length is lmin l ′ 2 < lmax and the ] over the triangles (see text for details). First, the first vector l ′ 1 is chosen from the annulus (shaded region), which is in the upper half plane and has the radius in the range [l 1,min , l 1,max ] (a bin of l 1 ). Then, the second vector l ′ 2 is chosen from the annulus which has position angle in the range arg(l ′ 1 ) arg(l ′ 2 ) π and has radius in the range [l 2,min , l 2,max ]. For the given pair of l ′ 1 and l ′ 2 , the third vector l ′ 3 is determined by the triangle condition l ′ lmax, the triplet of (l ′ 1 , l ′ 2 , l ′ 3 ) is accepted and otherwise discarded. The angle ϕ 12 is used for discussion in Appendix A. position angle satisfies arg(l ′ 1 ) arg(l ′ 2 ) < π (more precisely, l ′ 1 l ′ 2 if and only if arg(l ′ 1 ) = arg(l ′ 2 )). For a given pairs of l ′ 1 and l ′ 2 , we choose the third vector l ′ 3 from the condition l ′ 3 = −l ′ 1 − l ′ 2 and then accept the triplet of (l ′ 1 , l ′ 2 , l ′ 3 ) (otherwise discard it) if the length l ′ 3 is in the range satisfies the condition lmin l ′ 3 < lmax. Then we assign each set of three vectors, (l ′ 1 , l ′ 2 , l ′ 3 ), to one of the triangle configuration bins labelled by (l1, l2, l3) by sorting l ′ 1 , l ′ 2 and l ′ 3 in the ascending order. For any of the sets of three vectors chosen in this way, two of the three vectors are in the upper half plane of Fourier space, while the remaining one is in the lower half plane. Hence, we miss triangles with configurations for which two vectors are in the lower plane and the other in the upper. However, we can recover these triangles by just flipping the signs of all three vectors and obtain the same value of bispectrum, due to the symmetryκ l =κ * −l , which comes from the real condition of the lensing field. Thus, we do not need to count the latter cases, but just twice the number of actually counted triangles to obtain Ntrip, the number of independent triplets. For some of the following results, we will use the measured Ntrip when computing the Gaussian error contributions to the bispectrum covariance (the first term in equation 18). Although equation (16) gives a good approximation to Ntrip for the limit of li ≫ 1, the Ntrip directly estimated from the simulation properly takes into account the effect of discrete grids in the lensing map.
Using the table of three vectors (l ′ 1 , l ′ 2 , l ′ 3 ) obtained in the process above, we estimate the bispectrum by averaging Re[κ l ′ 1κ l ′ 2κ l ′

3
] from each realization. Note that the value ofκ l ′ i is taken from the field at the grid that has two coordinate components (l ′ ix , l ′ iy ), not from the field at the grid denoted by the arrow of the vector l ′ i in Fig. 3. Although we take the real part ofκ l ′ 1κ l ′ 2κ l ′ 3 for the average, this is not essential, because the estimated bispectrum satisfies the real condition to a very good approximation after the average over many triangles.
The dimensions of the resulting covariance matrices are: 16 × 16, 204 × 204 and 16 × 204 (or 204 × 16) for the power spectrum covariance, the bispectrum covariance and the cross-covariance, respectively. The 1000 realizations are enough to compute the 204×204 elements of the covariance matrix (see appendix of Takahashi et al. 2011), although a larger number of the realizations are ideally needed for a more accurate estimate of the covariance matrix. The situation will be worse in a case where more triangle configurations are considered, e.g. as in the case of lensing tomography where different redshift slices are further needed to include (Takada & Jain 2004). Hence an analytical approach to compute the covariance matrices is useful.
In fact, Hartlap et al. (2007) pointed out that, by assuming a multivariate Gaussian distribution for a statistical variable, the number of realizations used to estimate the covariance is important. They showed that the inverse covariance matrix can be biased if the covariance matrix is estimated from a finite number of the realizations. In our case, we use 1000 realizations to estimate the bispectrum covariance matrix for 204 triangle configurations (the dimension of the bispectrum covariance is 204 × 204), which may result in an overestimate of 10 per cent for the signal-to-noise ratio (S/N) for the bispectrum. However, we will show below the simulation results without any correction because the bias is not large comparing to the accuracy of our concern or we do not know whether the bispectrum estimators obey the multi-variate Gaussian distribution. For the three-halo term, we used the tree-level perturbation theory prediction. Note that the bispectra for l 4000 may be affected by the resolution limit of the ray-tracing simulations. 1000. An improvement of the halo model accuracy may be available by refining the halo model calculation, as performed in Valageas et al. (2012). However, we do not pursue this possibility in this paper.

Comparison of the simulation results and the halo model predictions
In Figs 5-7, we study the covariance matrices of the lensing bispectra for some representative triangle configurations of the 204 triangles. Note that, in these results, the two triangle configurations of the covariance have the same length(s) for at least one side length (e.g. l1 = l ′ 1 ). Hence the covariance terms O(BB) and O(P T ) in Fig. 1 are not vanishing for the covariances shown in these figures. First, in Fig. 5, we study the diagonal terms of the bispectrum covariance matrix for equilateral triangles, Cov[Beq(l), Beq(l)], as a function of the side length l. The points are the simulation results. The jagged feature at small l bins is due to the effect of discrete pixels in the lensing maps. The different curves are the halo model predictions for the covariance matrix (equations 18 and 19). The dotted curve shows the Gaussian error contribution that scales with P (l) 3 , where we use the number of triangles directly computed from the simulated lensing map. The long dot-dashed, short dot-dashed, long dashed and short dashed curves are the different non-Gaussian terms. For these calculations, we use the P (l) and Beq(l) directly estimated from the 1000 ray-tracing simulations, while we use the halo model in Section 3.5 to compute the the higher-order functions. The solid curve is the total power of the covariance matrix, the sum of the terms above. The figure clearly shows that the non-Gaussian errors become significant at l a few 100, and the model predictions are in good agreement with the simulation results if we include the HSV term. The HSV term dominates the other terms at l 1000. These findings are similar to the results in the power spectrum covariance (see figs 5-7 in Sato et al. 2009). Fig. 6 shows the covariance matrices for isosceles triangle configurations (left panel) and the off-diagonal components between different triangle configurations, equilateral and isosceles triangles (right panel). For this plot, we do not use the condition l1 l2 l3 for presentation purposes, but all the triangles shown here are indeed included in the covariance matrix elements we will use in the following analysis. Even for these more general triangle configurations, the model predictions including the HSV effect are in good agreement with the simulations. One may find that the covariance amplitudes or some non-Gaussian terms peak at some particular value of l3 in the x-axis. This happens when the isosceles triangles have higher symmetry, equilateral shape; the number of independent triangles is smaller for higher-symmetry triangles, leading to the greater covariance amplitudes.
In Fig. 7, we study the cross-covariance between the power spectrum and the bispectrum of equilateral triangles, Cov[P (l), Beq(l)], as a function of the multipole bin l. There is no Gaussian error contribution because the cross-covariance arises from the five-point functions. The figure again shows that the model predictions including the HSV term well reproduces the simulation results. The covariance amplitude is shown in the unit of l 6 Cov[Beq(l), Beq(l)], since l 3 Beq(l) gives the contribution to the skewness κ 3 . The points are the simulation results estimated from the scatters of 1000 simulations. The other curves are the halo model predictions, which are computed based on the method described in Section 3.2. The dotted, long dot-dashed, short dot-dashed, and long dashed curves are the contributions that are proportional to P 3 , B 2 , P T and P 6 , respectively (also see Fig. 1). For thick curves, we used the power spectra P (l), the bispectra B(l) and the number of triangles directly estimated from the simulations. For comparison, the thin curves show the results if we use the halo model for P (l) and B(l) as well as equation (16) for the number of triangles. The short dashed curve is the HSV contribution, which dominates over other terms at multipole bins, l 2000. The solid curves are the total contribution, the sum of all the terms. Figure 6. Similar to the previous plot, but for different triangle configurations. Left-hand panels: the covariance matrix for isosceles triangle configurations, B(l 1 , l 2 , l 3 ), with l 1 = l 2 . The different panels are for different side lengths l 1 and l 2 ; within each panel, the covariance matrix is shown as a function of l 3 . Note that the covariance matrix is shown in the multipole range where the triangle conditions |l j − l k | l i l j + l k are satisfied, but the condition l 1 l 2 l 3 is not imposed. Right-hand panels: the covariance matrix elements between the bispectra of isosceles and equilateral triangles. The Gaussian terms of P 3 denoted by the vertical dotted lines appear at a particular value of l 3 , where the two triangles become the same, i.e. equilateral triangles with l 1 = l 2 = l 3 .

Figure 7.
Cross-covariance between the power spectrum and the bispectrum of equilateral triangles, Cov[P (l), Beq(l)], as a function of l. The covariance amplitude is plotted in the units of l 5 Cov[P, Beq], because l 2 P and l 3 B contribute to κ 2 and κ 3 , respectively. Note that there is no Gaussian error contribution to the cross-covariance, because it arises from the five-point correlation functions. The different curves are the model predictions that are computed based on the method in Section 3.3.
In Fig. 8, we compare the halo model predictions with the simulation results for all the matrix elements of the power spectrum covariance, the bispectrum covariance and the cross-covariance, in one figure. To do this, we use the correlation coefficients of the covariance matrices defined as where X and Y are the power spectrum or the bispectrum, and the subscript i or j in X or Y denote the ith multipole bin or the ith triangle configuration; Xi = P (li), Xi = B(li), and so on. The diagonal components rii = 1 by definition. If rij = 0, it means no correlation between the spectra Xi and Yj, while the the higher values of rij mean stronger correlations. For illustration purposes, we use the following indices of the 204 triangles so that the different triangles are indexed in increasing order of l3: (1, 1, 2), (1, 2, 2), (2, 2, 2), where we have used the 16 logarithmically spaced multipole bins of l. Note that, for a given l3-bin, the other multipole bins (l1, l2) are listed in increasing order of l1 (l1 l2 for each triangle index). With this triangle index, the higher-index triangle configuration corresponds to the triangles having higher multipoles or larger side lengths. The figure shows that the halo model well reproduces the two-dimensional features of the covariance matrices seen from the simulations. 10 Most of the off-diagonal terms of the bispectrum covariance consist of only the term O(P6) and the HSV term, because in general the shapes of two triangles are different from each other in contrary to the cases of Figs 5 and 6. The correlation coefficients become greater at higher multipoles, almost rij ≃ 1. For comparison, Fig. 9 shows the results without the HSV term in the halo model predictions (the simulation results are the same to the previous figure), where the discrepancy is clear.

Information content of the lensing bispectrum
As we have studied, the non-linear structure formation induces non-Gaussian errors in the weak lensing field, provoking significant correlations between the power spectra of different multipoles and the bispectra of different triangle configurations as well as significant crosscorrelations between the power spectra and the bispectra. Then a more fundamental, important question arises. How much additional information do the lensing bispectra carry to the lensing power spectrum? Can joint measurements of the power spectra and bispectra recover the information content of the Gaussian field, which the primordial density field of structure formation should have had as in the CMB case? In this section, we address these questions. A useful quantity to quantify the impact of the non-Gaussian errors is the expected signal-to-noise ratio (S/N) for a measurement of the lensing power spectra and bispectra in a given survey that is characterized by its area and shot noise parameters. The S/N is sometimes called the information content (Tegmark et al. 1997) (also see Takada

S N
where the summation indices i, j run over multipole bin indices up to a given maximum multipole lmax, and [C P ] −1 is the inverse of the power spectrum covariance matrix. The inverse of S/N is equivalent to a precision of measuring the logarithmic amplitude of the power spectrum up to a given maximum multipole lmax, assuming that the shape of the power spectrum is perfectly known. The S/N is independent of the multipole bin width as long as the bin width is fine enough to capture the shape of the lensing power spectrum (on the other hand, the relative strength of the non-Gaussian errors to the Gaussian errors depends on the width). Similarly, the S/N for the bispectrum measurement or the information content about the bispectrum amplitude is defined as where the summation indices i, j run over triangle configurations, and we include all the triangle configurations whose side lengths are smaller than a given maximum multipole lmax.
We also consider the S/N for a combined measurement of the lensing power spectra and bispectra up to a given lmax. In the presence of the non-Gaussian errors, the total S/N is not simply a sum of the two estimates of S/N of the power spectra and the bispectra due to the cross-covariance. To study this, we first define the data vector for the joint measurement as D = P1, P2, · · · , Pn b , B1, B2, · · · , Bi triang,max .
The covariance matrix for the data vector D is given as where the C P B is the cross-covariance between the power spectrum and the bispectrum. Then, the total (S/N)P +B is similarly defined as Fig. 10 shows the expected S/N for measurements of the power spectra and the bispectra for a survey area of 25 square degrees (i.e. the area of the ray-tracing simulation), as a function of the maximum multipole lmax up to which the power spectrum and/or bispectrum information are included. The minimum multipole is fixed to lmin = 72. We do not include the shot noise contamination to the error covariance matrices, so the results solely correspond to the cosmological information contents. The circle, triangle and square symbols are . The minimum multipole is set to l min = 72. We do not include the shape noise contamination here -it is shown in the next figure. The circle, triangle and square symbols are the simulation results for P , B and P + B measurements, respectively, computed from the 1000 realizations. The thick short-dashed, long-dashed and solid curves are the corresponding halo model predictions. The corresponding thin curves are the results without the HSV contributions. For comparison, the dotted curve shows the S/N for the power spectrum for the Gaussian field, which the primordial density field should have contained. Note that the simulation results for B and P + B could be overestimated by about 10 per cent due to a finite number of the simulation realizations used to estimate the covariance matrices (Hartlap et al. 2007). the simulation results for the S/N of the power spectra, the bispectra and the joint measurements, respectively, which are computed using the 1000 realizations. The thick/thin short-dashed, long-dashed and solid curves are the halo model predictions with/without the HSV terms. First of all, the lensing bispectra add new information content to the power spectrum measurement. To be more quantitative, adding the bispectrum measurement increases the S/N by about 50 per cent for lmax ≃ 10 3 compared to the power spectrum measurement alone. Note that the lmax of a few thousands is the typical maximum multipole for upcoming weak lens surveys. This improvement is equivalent to about 2.3 larger survey area for the power spectrum measurement alone; that is, the same data sets can be used to obtain the additional information, if the bispectrum measurement is combined with the power spectrum measurement. Secondly, the halo model predictions are in nice agreement with the simulation results. Note that the total S/N for the joint measurement (P + B) is close to the linear sum of the S/N values ((S/N)P and (S/N)B), not the sum of their squared values (S/N) 2 , due to the significant cross-covariance between P and B (see Appendix C in Takada & Bridle 2007, for the similar discussion). If ignoring the cross-covariance, adding the bispectrum measurement does not much improve the S/N (only by 5 per cent or so). Hence it is important to take into account the correlation between the two measurements.
Next, let us compare the result above with the case of a Gaussian random field, which is shown in the dashed curve in Fig. 10. The S/N for a Gaussian field is equivalent to the number of independent Fourier modes up to a lmax for a given survey area. The figure clearly shows that the joint measurement of the power spectrum and the bispectrum does not recover the full information content of the Gaussian field. This implies that the higher-order statistics beyond the bispectrum are also important to recover the full information content. One may argue that the initial memory of the field cannot be recovered due to the non-linear structure formation. However, we would like to note that, if ignoring the HSV contribution to the covariance, adding the bispectrum can recover about 75 per cent of the Gaussian information, as shown by the thin curves. Hence the loss of the information contents is mostly due to the the HSV contribution. As discussed in Section 3.4, the HSV alters the overall amplitude but preserves the shape of the lensing spectra. Hence the HSV may give the worst case degradation of the amplitude parameter, but may not cause any serious degradation of parameters that are sensitive to the shapes of the lensing spectra. A genuine impact of the HSV on cosmological parameters needs to be further studied and this is our future work.
In Fig. 11 and Table 1, we show the S/N expected for the upcoming wide-field weak lensing surveys, the Subaru Hyper Suprime-Cam (HyperSC) survey and the DES, which are characterized by the survey area, the mean source redshift and the mean number density of source galaxies of Ωs = 1500 sq. degrees,zs = 1 andng = 20 arcmin −2 for the HyperSC survey, while Ωs = 5000 deg 2 ,zs = 0.7 and ng = 10 arcmin −2 for the DES, respectively. Here we employ the halo model to compute the S/N and assume a circular survey geometry for simplicity. The figure and table show that these surveys promise a significant detection of the lensing bispectrum; (S/N) ≃ 26 or 29 for the HyperSC or the DES, respectively, when assuming lmax ≃ 2000 and including the shot noise effect. It also means that the theoretical prediction of the lensing bispectrum needs to be as accurate as a few per cent for the upcoming surveys. We find that the bispectrum adds   Cumulative S/N of the power spectrum (P ), the bispectrum (B) and the joint measurement (P + B) expected for the upcoming weak lensing surveys, the Subaru Hyper Suprime-Cam survey (HyperSC) and the DES, as in Fig. 11. Here we consider lmax ≃ 1000 and 2000 for the maximum multipole (more exactly, the bins of lmax = 1245 and 2268). The column denoted by 'w/o' or 'with' in the row 'shot noise' gives the S/N values with and without the intrinsic shape noise contribution to the covariances. The percentage in the parenthesis shows the improvement of S/N for the joint measurement (P + B) compared to the power spectrum alone (P ).
new information, increasing the total S/N by about 20 -30 per cent compared to the power spectrum alone, which is equivalent to a factor of 1.4 -1.7 larger survey area. Fig. 11 also shows that the HSV is significant for these surveys. The thin curves in the upper plots of the left-panel are the S/N computed by scaling the values for 25 deg 2 in Fig. 10 assuming S/N ∝ Ω 1/2 s . Since the covariance terms except for the HSV term scale as 1/Ωs, the differences between the thick and thin curves are due to the HSV term which depends on the survey area via the shape of the matter power spectrum convolved with the survey window function (see the discussion below equation 14). The S/N values are smaller than the naively-scaled results, which means that the HSV decreases more slowly than 1/Ωs. The upper plot of the right panel shows the similar plot, but for the DES results with scaled values for the HyperSC in the left panel. The S/N are smaller than the scaled HyperSC results, because the typical source redshift ofzs = 0.7 for the DES is lower than that of the HyperSC ofzs = 1 and the DES is more sensitive to the non-linear density fluctuations. The HSV has significant influence on the S/N even in the presence of the shot noise as shown in the lower panels of the figure. Although the shot noise leads to a saturation of the S/N at large multipoles, note that the systematic effects such as the highly non-linear clustering effect and/or the baryonic effect become more significant at these high multipoles (White 2004 Fig. 11, but for (S/N) 2 instead of S/N (the points denote the central value of each multipole bin). The short-dashed and dotted curves show the differential contribution to the (S/N) 2 at each multipole bin, with and without the HSV effect, respectively. The solid curve shows how the (S/N) 2 value is recovered by adding the PCA eigenmodes. Note that the PCA results are shown as a function of the order of the PCA eigenmodes, where the PCA eigenmodes are ranked in decreasing order of the differential contribution to the (S/N) 2 (see equation 33). Right-hand panel: the projection matrix |S ai | for the first eight eigenmodes, where the index a denotes the ath eigenmode.

Principal component analysis of the lensing covariance matrices
A principal component analysis (PCA) of the power spectrum and bispectrum covariance matrices is useful to quantify how the different power spectra and/or bispectra are correlated with each other and how many independent modes or triangles contribute to most of the information contents (Takada & Jain 2009) (also see Scoccimarro 2000, for the the 3D bispectrum case). Since the covariance matrix is symmetric by definition, it can always be decomposed as where X is the lensing power spectrum or bispectrum, λ X a is the ath eigenvalue or principal component, (S X ) −1 = (S X ) T , k S X ik S X jk = δ K ij and S X is normalized so as to satisfy j (S X ij ) 2 = 1. The matrix S X ai is considered as the projection matrix as it describes how the power in the ith multipole bin or triangle configuration is projected onto the ath eigenmode. Using this representation, the inverse of the covariance matrix is given by [C X ] −1 ij = S X ia (1/λ X a ) 2 S X ja . Hence, the S/N values (equations 27 or 28) can be rewritten as where Xi is either P (li), B(li) or the joint (P + B). Thus, since (1/λ X a ) i S X ai Xi can be considered as the S/N for the ath eigenmode, the above equation expresses the total S/N (for a given lmax) as a sum of contributions from the independent eigenmodes. We can then re-order the eigenmodes such that the lower eigenmode has the higher contribution to the S/N. Fig. 12 shows the PCA results for the power spectrum measurement for the Subaru HyperSC-type survey including the shot noise effect as in Fig. 11. This figure can be compared with Fig. 5 in Takada & Jain (2009), where the HSV effect was not included. The short-dashed and dotted curves show each contribution of the ith multipole bin to the total (S/N) 2 with and without the HSV, respectively. It should be stressed that the power spectrum of l ∼ 10 3 has a local minimum, due to the significant HSV contribution.  (26); the larger indices correspond to triangles with the larger side length of l 3 . The vertical shaded regions denote the triangles with the same l 3 bin, as indicated by the top label of 'l 3 '-bin values. Upper panel: the short-dashed curve shows the differential (S/N) 2 of each triangle configuration, while the long dashed curve shows the cumulative (S/N) 2 as in the left-lower panel of Fig. 11, but for the finer multipole bins (stepped by each triangle configuration). The solid curve shows the cumulative (S/N) 2 as one adds the different PCA eigenmodes. About 70 eigenmodes out of 204 triangles can recover 90 per cent of the total (S/N) 2 of the bispectrum measurement over 72 l 8748.8. Lower panel: the solid curve shows the number of independent triangles available from Fourier modes for the HyperSC-type survey (1500 deg 2 ), for each triangle configuration in the x-axis.
The long-dashed curve shows the total (S/N) 2 as a function of the maximum multipole lmax, up to lmax = 8748.8 (the 16-th multipole bin). On the other hand, the solid curve shows how the cumulative (S/N) 2 value increases by adding the ath eigenmode (equation 33) 11 . The figure shows that, among the 16 multipole bins, including about seven eigenmodes (about half of the multipole bins) can recover about 90 per cent of the total (S/N) 2 . The other eigenmodes are relatively less important due to strong correlations between the different multipole bins. The right panel shows the projection matrix |Sia| for each multipole bin, showing how the neighboring multipole bins are correlated with each other. The projection matrix for each eigenmode peaks at some multipole bin, but has tails to different multipole bins. In particular, the eigenmodes around a few 10 3 have longer tails, reflecting significant correlations between different multipole bins due to the non-Gaussian errors. Fig. 13 shows the PCA analysis for the bispectrum covariance matrix, where we include the 204 triangles over 72 l 8748.8. The different triangle configurations contribute to the total (S/N) 2 in a complex way. We use the triangle index, given by equation (26), and the higher-index triangle corresponds to triangle with larger side length of 'l3', under the condition l1 l2 l3. However, for a given l3 bin, either of the two side lengths l1 or l2 can be much smaller than l3, and such a triangle configuration generally has a smaller contribution to (S/N) 2 . As a result, the differential (S/N) 2 curve has jagged features. The sold curve shows the cumulative (S/N) 2 by adding the new PCA eigenmodes, and manifests that about 70 eigenmodes, one third of the 204 triangles, carry 90 per cent of the total (S/N) 2 . Fig. 14 shows which triangle configurations contribute to the first 6 PCA eigenfunctions, which therefore have the 6 highest contributions to the (S/N) 2 . More exactly, each panel shows the projection matrix |Sia| for the ath eigenmode (a = 1, 2, . . . , 6). In most panels, the projection matrix |Sia| has a single or a few peaks, reflecting that the higher-S/N eigenmode arises mostly from a single or a few different triangles. For example, the highest (S/N) 2 eigenfunction arises from isosceles triangles, which have side lengths of a few 10 3 (11-and 14-th Figure 14. The curve in each panel shows which triangle configurations contribute to each PCA eigenmode for the bispectrum measurement in Fig. 13. We consider the first 6 PCA eigenmodes. They are labeled by the index a from the highest contribution to the (S/N) 2 and the absolute values of |S ia | are plotted, where S ia is the projection matrix of the ath eigenfunction onto the ith triangle configuration. The triangle configurations shown on the right side of each panel are the triangle configurations with S ia = 0, where the triangle sizes are plotted on the linear scale and the thickness of the lines corresponds to the amplitude of |S ia |. The configurations with the largest contribution are labelled in the side length bins (l 1 , l 2 , l 3 ) at top. Note that we are using logarithmically-spacing multipole bins in this paper, so about 3/4 triangles available in the Fourier space are isosceles triangles. multipole bins). However, note that we have used rather wide, logarithmically-spaced multipole bins, and most of triangle configurations available from Fourier space are close to isosceles triangle configurations. Hence, an exact shape of triangles giving a large contribution to the total (S/N) 2 slightly changes with the bin width. If we take a finer multipole bin, the triangles with the highest (S/N) 2 may differ from isosceles triangles.

CONCLUSION AND DISCUSSION
We have studied the covariance matrices of the lensing power spectrum and the bispectrum by using both 1000 ray-tracing simulations and the analytical halo model. We have found that there are significant non-Gaussian error contributions to the lensing covariance matrices; the power spectra or the bispectra at higher multipoles than a few hundreds are highly correlated with each other. In particular, we have shown that the mass density fluctuation at scales comparable with or outside a survey region causes significant non-Gaussian error contributions, which we call the HSV and has not been fully studied in the literature. With the HSV contributions, the halo model predictions reproduce the covariance matrices measured from the simulations (Figs 4-8).
Then, we have addressed how much information the lensing bispectrum adds to the power spectrum by including all the triangle configurations available up to a maximum multipole lmax. Adding the bispectrum measurement improves the S/N by about 20-50 per cent compared to the power spectrum alone measurement, at lmax of a few 10 3 (Fig. 10). We have also studied the prospects for upcoming weak lensing surveys, including the shot noise contamination due to intrinsic galaxy shapes ( Fig. 11 and Table 1). The improvement in S/N is equivalent to about 1.4-2.3 larger survey area for the power spectrum measurement alone. Hence, our results show that the same imaging data can be used to improve the constraining power of cosmological parameters, if we combine the power spectrum and bispectrum measurements. We have also found that the HSV effect is significant, leading to a large degradation in the S/N. By using a PCA of the covariance matrices, we have shown that about 1/3 eigenmodes of all the triangle configurations over the range of 10 2 l 10 4 (70 eigenmodes compared to the 204 triangles for the multipole binning we assumed) carry most of the total information for the bispectrum measurement (Fig. 13). Thus, our results give a quantitative answer to the longstanding question of how the bispectrum can be useful and complementary to the power spectrum by fully taking into account non-Gaussian errors and all triangle configurations.
Future surveys such as the Subaru Hyper Suprime-Cam survey or the DES allow for a significant detection of the lensing bispectra ( Fig. 11 and Table 1). However, there are some practical issues for making the bispectrum measurement feasible for future surveys. First, we must take into account effects of the complex survey geometry and masked regions such as saturated bright stars. This can be done by extending the method in Hikage et al. (2011) for the power spectrum measurement to the bispectrum, but has not yet been fully addressed in the literature. It is also important to explore a method of cleanly decomposing three-point correlations of E/B modes in the presence of the survey window function. An alternative approach is to measure the real-space three-point correlation functions, rather than the bispectrum, which does not suffer from the geometrical problem. However, there are even more significant correlations between the three-point correlations of different triangles, and the estimate of the covariance matrices involves multi-dimensional integrations of the bispectrum covariance if one uses the halo model approach. Although the alternative to the theoretical approach may be to use a sufficient number of simulations, it will not be feasible to construct the covariance matrices for many cosmological models. Hence, which of the Fourier-or real-space is more useful/tractable for the three-point correlation measurements is still an open issue.
Although we have focused on the S/N as a measure of the information content, this is just one measure to quantify the complementarity of the lensing bispectrum. Specifically, the S/N quantifies the expected precision of the amplitude, assuming that the shapes of the power spectrum of the bispectrum are known. Hence, the S/N is not necessarily the best measure to quantify the power of the bispectrum for constraining cosmological parameters, especially parameters that are sensitive to the shape of the bispectrum (also see Takada & Jain 2009, for the similar discussion). For instance, Takada & Jain (2004) showed that the bispectrum can significantly improve the accuracy of parameter estimate from the power spectrum measurement alone, by efficiently breaking parameter degeneracies, while the S/N of the bispectrum itself is smaller than that of the power spectrum by a factor of few up to lmax of a few 10 3 (see Fig. 5). Although they ignored the non-Gaussian error contributions to the covariances, the relative amplitude of the S/N between the power spectrum and the bispectrum is similar to our case with the full non-Gaussian terms, because the non-Gaussian terms degrade both the power spectrum and the bispectrum. Thus we can expect a similar improvement in cosmological parameters when including the bispectrum information, even with the inclusion of non-Gaussian errors. Exploring the genuine usefulness of the lensing bispectrum for cosmology is very interesting -the full forecasts for various cosmological parameters will be presented in future work.
In particular, to do such parameter forecasts, it would be much more interesting and useful if including tomographic information of the lensing field that is available from photometric redshift of source galaxies. It has been found that the lensing tomography can significantly improve the constraining power by recovering the radial-direction information of the lensing field (e.g. Takada & Jain 2004). However, for the lensing bispectrum tomography, we need to further include all the triangle configurations in different multipole and redshift bins in order to estimate the full potential. For example, if we consider three redshift bin tomography, the total number of triangles is 3 3 × 204 = 5508 for the same multipole binning as used in this paper. Thus, even 1000 realizations are not enough to reliably estimate the covariance matrices of the lensing bispectra. Hence we believe that the analytical formula we developed in this paper are essential to address these issues.
Perhaps as important as the improvement in statistical precision, the bispectrum might enable a self-calibration of systematic errors inherent in the lensing tomography measurements such as photometric redshift errors and imperfect galaxy shape measurement. Lensing bispectra depend on systematic errors in a different way from the power spectrum, but the two spectra share the same large-scale structure, hence the same cosmological parameters that describe the underlying true cosmology. Thus combining the power spectrum and bispectrum measurements can be used to self-calibrate systematic errors and improve cosmological constraints (Huterer et al. 2006). Again, it is important to fully take into account the non-Gaussian errors in order to quantify how well self-calibration works for upcoming lensing surveys.
Planned multi-colour imaging surveys can also be used to measure other cosmological probes such as the abundance of galaxy clusters and the correlation functions of galaxies. Recently, Oguri & Takada (2011) proposed a new method of using the halo-shear correlations or the stacked lensing signals around massive haloes as a function of the cluster redshifts, and showed that the stacked lensing tomography allows us to constrain cosmological measurements to a high precision. In halo (cluster)-shear correlations, the signals arise from the largescale structure at the cluster redshift. As we have shown, massive clusters are key observables for understanding the HSV effect on the lensing power spectra. It would be interesting to pursue how the lensing power spectrum and bispectrum measurements can be combined with the cluster observables to improve cosmology parameter estimate by calibrating the HSV effect. The formulation we have developed in this paper can be straightforwardly extended to the halo-shear bispectra such as the halo-shear-shear and halo-halo-shear three-point correlation functions.
Finally, there are several effects that we have ignored in this paper and need to be studied. One is that we only considered a simple circular-shaped survey geometry to compute the HSV term. For a general geometry of the survey, the Fourier transformed survey window function becomes anisotropic and causes additional apparent correlations between the different spectra. Since the large-scale density fluctuations which contribute to the HSV are in the linear regime, it may be straightforward to take into account the effect by extending our formulation. Secondly, we have used the flat-sky approximation, which is not valid for very wide surveys such as the LSST. For a full-sky survey as the ultimate case, the HSV effect is caused by the horizon-scale perturbations with wavenumbers comparable to the matter-radiation equality wavenumber keq. In this case, the HSV terms decrease more rapidly with increasing survey area than other covariance terms. Hence the HSV effect may not be that significant for such an all-sky survey. This is worth studying and the formulation and the method we developed in this paper can be used for such studies.

A2 Bispectrum estimator and the covariance matrix
The bispectrum is given as a function of triangle configurations. We use a parametrization of three side lengths (l1, l2, l3) to specify a triangle configuration. The three parameters are enough, because all the triangles transformed by rotation, permutation and parity transformations in Fourier space are equivalent to yielding the same bispectrum in an ensemble average sense for a homogeneous and isotropic field as expected for the lensing field. Hence, within the framework of the discrete Fourier decomposition we described above, we can define an estimator of the bispectrum: where we have introduced the abbreviated notation q 123 ≡ q 1 + q 2 + q 3 and the summation runs over all the grids of q 1 , q 2 and q 3 in Fourier space. The function ∆q 123 (l1, l2, l3) denotes the selection function defined so that ∆q 123 = 1 if each vector has the target length, li − ∆li/2 qi li + ∆li/2 (i = 1, 2, 3) as well as the three vectors satisfy the triangle condition, q 123 = 0; otherwise ∆q 123 = 0. The quantity Ntrip is the number of triplets of grids (vectors) which form the triangle configuration of (l1, l2, l3). We average the estimator ofκq 1κ q 2κ q 3 over all the triangles that have the side lengths of (l1, l2, l3) within the bin widths. For the limits of li ≫ l f , the number of independent triplets of a target triangle configuration, Ntrip, can be estimated as where "independent triplets" means different combinations of triplets that form the target triangle configuration, and ∆li is the bin width around each multipole bin, ϕ12 is the angle between the two side lengths q1 and q2 in Fig. 3, and ∆ϕ12 is the angle extent allowed by the multipole bins as we explain below. First, 2πl1∆l1/(2π/Θs) 2 is the number of grids in the annulus, which has the radius of l1 with width of ∆l1, as in the coefficient in the power spectrum covariance (see equation 12). Then, for a grid of the vector q 1 , we can find other two vertices (two grids), denoted by the vectors q 2 and q 3 in Fig. 3, so that the two vectors have the target lengths, q2 = l2 and q3 = l3 within the bin widths. Given the bin widths, the number of grids, which are covered by variations in the two vectors q 2 and q 3 , can be estimated as (l2∆ϕ12 × ∆l2)/(2π/Θs) 2 , where ∆ϕ12 is the variation in the angle ϕ12 allowed by the variations in the two vectors q 1 and q 2 . Finally, the prefactor 2 in equation (A5) accounts for counter-part triplet of grids that are transformed by parity transformation ϕ12 → −ϕ12. As can be found from Fig. 3, the angle ϕ12 is given as where the two vectors have the target lengths within the bin widths: li − ∆li/2 qi li + ∆li/2 (i = 1, 2). If varying the side length q3 by ∆l3, we can find the angle variation as ∆ϕ12 ≃ (sin ϕ12) −1 l3∆l3 l1l2 = 2l3∆l3 2l 2 1 l 2 2 + 2l 2 1 l 2 3 + 2l 2 2 l 2 3 − l 4 1 − l 4 2 − l 4 3 .
where we have dropped terms including the Kronecker delta such as δ K q 12 , because we are not interested in such triangles with q 3 = 0 under the triangle conditions q 123 = 0. As we explained in Fig. 1, the bispectrum covariance has different contributions that arise from the terms proportional to P 3 , O(BB), O(P T ) and P6, respectively. These terms are further simplified as shown below.

A2.1 Gaussian error contribution to the bispectrum covariance
First, we consider the term in equation (A11) proportional to the power spectra cubed P 3 . We call this term the Gaussian error contribution, although the bispectrum itself is a measure of non-Gaussianity in the lensing field. Many previous works only considered this term when studying the lensing bispectrum, e.g. for parameter forecasts, mainly for simplicity. The Gaussian term is simplified as
= Ωs Ntrip P (l1)P (l2)P (l3) δ K where we use the facts that the power spectrum does not largely change within the bin width and that the terms including the Kronecker delta are non-vanishing only if the two triangles of the two bispectra have the same shape. The terms including the Kronecker delta give 1, 2 or 6 for general, isosceles and equilateral triangles, respectively, which correspond to the factor ∆ in equation (28) of Takada & Jain (2004). The Gaussian error covariance term contributes only to the diagonal terms of the bispectrum covariance matrix. Recalling the fact Ntrip ∝ 1/Ω 2 s (equation 16), this contribution scales with survey area as CovGauss ∝ 1/Ωs.

A2.2 Non-Gaussian error contributions to the bispectrum covariance
Other terms in equation (A11) are non-Gaussian error contributions to the bispectrum covariance that arise from the higher-order functions.
In the following, we derive further simplified expressions for each term.

A2.2.1 Non-Gaussian terms of O(B 2 )
We consider the terms that are proportional to the bispectra squared, the terms of O(B 2 ). Here we consider the first term of O(B 2 ) terms in equation (A11) as an example as 1 Ntrip Inserting this equation into Eq.(B2) and using the similar calculation we have used, we can find the flat-sky formula for the cross-covariance as in equation (20) but without the HSV term.

APPENDIX C: DEPENDENCE OF HSV ON SIMULATION BOX SIZE
As we have studied, the large-scale mass fluctuations of scales comparable with or larger than a survey area cause a significant contribution to the non-Gaussian errors of the lensing field. However, the ray-tracing simulations we have used are generated from a finite volume of N-body simulations, and do not include contributions from the modes of scales beyond the N-body simulation box. Hence the simulation results we showed in the main text might underestimate the genuine effect of the HSV, which we estimate in this appendix. As described in detail in Sato et al. (2009), the simulated convergence maps are generated from N-body simulations with box size of 240h −1 Mpc up to the redshift of zs = 1. Thus the simulation maps do not include Fourier modes of k < k box ≡ 2π/240 ≃ 0.026 hMpc −1 . Using our halo model approach, we can estimate the HSV effect seen in the ray-tracing simulations by inserting k box into the lower bound of the integral of the HSV terms in equations (14), (19) and (21) instead of zero. The thick curves in Fig. C1 show the resulting S/N, which differs from the simulation results by 10-20 per cent. Comparing the thick and thin curves manifests that the box size effect is larger at higher multipoles, because the HSV contribution at higher multipoles arise preferentially from lensing structures at higher redshifts (Fig. 2) and therefore the finite box effect is more important at higher redshift (as the ratio of the simulation box to the ray-tracing simulation area becomes smaller). A possible reason of the discrepancy between the halo model and the ray-tracing simulation results is that the halo model calculation ignores some non-Gaussian contributions; for instance, we included only the one-halo terms to the four-, five-and six-point correlation functions in the covariance calculation. Hence, the halo model tends to underestimate the non-Gaussian error amplitudes, and in turn to overestimate S/N. However, we would also like to note that the ray-tracing simulations may not be accurate enough at higher multipoles beyond l ≃ 6000 (see Sato et al. 2009). Most importantly, in this paper, we were able to develop the model to describe the non-Gaussian covariance to a 10-20 per cent accuracy, and therefore we leave an issue on the 10-20 per cent discrepancy for future work.