© ESO 2019

1. Introduction

Weak gravitational lensing is recognised as a powerful probe of the large-scale structure of the Universe. Its reach, however, will always be limited by the accuracy to which terrestrial and astrophysical contaminating signals can be controlled. Known sources of astrophysical systematics include the intrinsic alignment of neighbouring galaxies (see Joachimi et al. 2015, and references therein) and the impact of baryon feedback when modelling the non-linear matter power spectrum (Semboloni et al. 2011) as well as the more subtle effect of the clustering of background “source” galaxies (Schneider et al. 2002a). Known sources of terrestrial systematics arise from residual distortions resulting from uncertainty in the point-spread function (PSF) model (Hoekstra 2004), biases in the adopted source redshift distributions (Hildebrandt et al. 2012), object selection bias (Hirata & Seljak 2003), shear calibration bias (Heymans et al. 2006) and detector-level effects (Massey et al. 2014; Antilogus et al. 2014). As weak lensing surveys have grown in size, the list of known sources of error has also grown, with accompanying mitigation strategies (see Mandelbaum 2018). This progress is impressive, but there will always be the possibility that hitherto unknown sources are contaminating the cosmic shear signals that we observe. In this paper we therefore explore the sensitivity of the “COSEBIs” weak lensing statistic to blindly uncover a range of different contaminating signals.

Complete Orthogonal Sets of E/B-Integrals, “COSEBIs”, were defined by Schneider et al. (2010). They provide a complete set of filter functions which cleanly separate a measured cosmic shear signal into its curl-free (E-mode) and divergence-free (B-mode) distortion patterns over a finite angular range. Weak lensing can only produce E-modes¹, and as such any detected B-modes in the measured cosmic shear signal will have a non-lensing origin. The most popular statistic used in current cosmic shear analyses are the shear two-point correlation functions, ξ_±, (2PCFs; Jee et al. 2016; Joudaki et al. 2017a; Hildebrandt et al. 2017; Troxel et al. 2018a). As these direct measurements of the cosmic shear signal mix E and B modes, other methods are required in order to extract and identify any contaminating non-lensing signal through its B-mode distortion pattern.

A range of different statistics exist to filter E/B-modes in 2PCFs, for example, aperture mass statistics (Schneider et al. 2002b), ξ_E/B (Crittenden et al. 2002) and ring statistics (Schneider & Kilbinger 2007). The aperture mass statistics and ξ_E/B rely on knowing the 2PCFs at either very small angular separations (where the galaxy images blend) or very large angular scales, (beyond the surveyed area). Using these statistics therefore results in biased estimates of E/B-modes. For the aperture mass statistic, the E/B-mode leakage is ∼10% for a typical case (Kilbinger et al. 2013). Both ring and aperture mass statistics suffer from a loss of information due to their filtering method.

Alternatives to real-space estimators decompose the cosmic shear signal into its E and B-mode convergence power spectrum. Quadratic estimators can be used (Köhlinger et al. 2016), but this method is sensitive to the modelling of the noise and is also challenging to use to estimate the power at large Fourier modes due to its computational speed (Köhlinger et al. 2017). Faster methods estimate “pseudo” power spectra where in an ideal case the E/B- power spectra can be easily separated. Unfortunately the presence of masks mixes Fourier modes, and hence E/B-modes, making this method sensitive to the modelling of the mask (Asgari et al. 2018; Hikage et al. 2019). Power spectra can also be estimated from 2PCFs, if the 2PCFs are known over all scales. In practice this is not feasible, hence the integrals over 2PCFs are truncated, which can produce biases in the estimates (van Uitert et al. 2018). Alternatives to band-power spectrum estimation from 2PCFs have also been suggested (Becker & Rozo 2016), which attempt to minimize the information leakage from the out-of-range angular scales. Another approach to power spectrum inference uses hierarchical Bayesian modelling (Alsing et al. 2016, 2017). Although this method is not sensitive to masking effects, it is highly computationally expensive as it relays on estimating posterior density distributions for all combinations of power spectra, which in turn can produce inaccuracies in the analysis.

In this paper we adopt the COSEBIs statistic as it is the only method that can cleanly, without loss of information, separate E and B-modes over a finite angular range from realistic lensing survey data. They are also efficient as a small number of COSEBIs modes (∼5 per tomographic redshift bin) can essentially capture the full cosmological information (Asgari et al. 2012). With data compression, using linear combinations of the tomographic COSEBIs modes that are most sensitive to the parameters to be estimated, the total number of data points can also be significantly decreased (Asgari & Schneider 2015). This compression then makes the method less sensitive to the accuracy to which the covariance matrix of the data can be estimated from numerical simulations.

COSEBIs have been used to analyse the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS; Kilbinger et al. 2013; Asgari et al. 2017), finding significant B-mode signals in the tomographic analysis that were not detected by a range of other systematic analyses (Heymans et al. 2012). The COSEBIs statistic is therefore more sensitive and stringent in detecting B-mode distortions. It is not immediately apparent, however, how a COSEBIs B-mode detection can be used in order to uncover the origin of the observed non-lensing distortions. In contrast, the ξ_B and aperture mass statistics are rather intuitive. For example. a peak in the measured B-mode at the angular scale of the CCD chip can be readily associated with an issue on the chip-level. It is also unclear how detected COSEBIs B-modes impact the cosmological parameters from the measured E-modes. For example, is a significant high-order COSEBIs B-mode detection an issue, when all the cosmological information is contained in the first five COSEBIs E-modes?

By using a range of different simulated systematic errors and analysing three public weak lensing surveys, this paper explores how B-mode statistics can be used to diagnose data-related systematic errors as follows. We describe the COSEBIs, ξ_E/B and compressed COSEBIs (CCOSEBIs) statistics as well as their covariance matrices in Sect. 2. In Sect. 3, we introduce the three public weak lensing surveys that we analyse; the Kilo-Degree Survey (KiDS-450, Hildebrandt et al. 2017), the science verification data from the Dark Energy Survey (DES-SV, Dark Energy Survey Collaboration 2016) and CFHTLenS (Heymans et al. 2013), presenting a full B-mode analysis of these surveys in Sect. 4. We then use mock weak lensing surveys to explore how the COSEBIs and ξ_E/B statistics respond to a range of different observationally motivated systematics, introduced in Sect. 5, with results presented in Sect. 6. We compare the results for the mocks and real data in Sect. 7 and conclude in Sect. 8. In Appendix A we discuss the biases that exist in published 2PCF analyses that arise from the angular binning of the 2PCFs. We also show how these biases can be mitigated. Appendix B determines the σ₈ − Ω_m degeneracy direction for a CCOSEBIs analysis of KiDS data. We discuss how to optimise B-mode null-tests using differing selections of the data vector in Appendix C and present supplementary material for the tomographic data analysis in Appendix D.

2. Methods

The most familiar two-point statistics used in cosmic shear analysis are the shear two-point correlation functions, ξ_±, which correlate γ_t/×, the tangential and cross components of shear, of two galaxies separated by an angle θ in the sky. They are defined as

$$\begin{array}{cc}\hfill {\xi }_{\pm }(\theta )& =⟨{\gamma }_{\mathrm{t}}{\gamma }_{\mathrm{t}}⟩(\theta )\pm ⟨{\gamma }_{\times }{\gamma }_{\times }⟩(\theta ).\hfill \end{array}$$(1)

In practice, galaxy ellipticities, ϵ, are measured with differing accuracies, accounted for using weights, w. In this case, an unbiased estimator for ξ_± is given by

$$\begin{array}{c}\hfill {\widehat{\xi }}_{\pm }(\theta )=\frac{{\sum }_{\mathit{ab}}{w}_a{w}_b[{ϵ}_{\mathrm{t}}({\mathit{x}}_a){ϵ}_{\mathrm{t}}({\mathit{x}}_b)\pm {ϵ}_{\times }({\mathit{x}}_a){ϵ}_{\times }({\mathit{x}}_b)]}{{\sum }_{\mathit{ab}}{w}_a{w}_b},\end{array}$$(2)

where the sum goes over all galaxy pairs in an angular bin labelled as θ (see Appendix A for binning choices). w_a is the weight associated with the measured ellipticity at x_a and ϵ_t/× are the tangential and cross components of the measured ellipticity (Schneider et al. 2002b). Here the ellipticity is defined such that its expectation value is equal to the reduced shear, in absence of systematics (Schramm & Kayser 1995; Seitz & Schneider 1997). If the ellipticity measurements require a multiplicative correction, m (see for example Miller et al. 2013), then the correlation functions may be calibrated by dividing them with the following correction,

$$\begin{array}{c}\hfill 1+K(\theta )=\frac{{\sum }_{\mathit{ab}}{w}_a{w}_b(1+m_a)(1+m_b)}{{\sum }_{\mathit{ab}}{w}_a{w}_b}·\end{array}$$(3)

Theoretically the 2PCFs can be calculated through their relation to the shear power spectrum, P_γ,

$$\begin{array}{cc}& {\xi }_{+}(\theta )={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{J}_0(\ell \theta )P_{\gamma }(\ell ),}\hfill \\ \hfill & {\xi }_{-}(\theta )={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{J}_4(\ell \theta )P_{\gamma }(\ell ),}\hfill \end{array}$$(4)

where ℓ is the Fourier conjugate of θ and J₀ and J₄ are the ordinary Bessel functions of zeroth and fourth order (see Kaiser 1992, for example). The shear power spectrum is in turn related to the three-dimensional matter power spectrum. This relation can be simplified assuming a flat-sky and Limber approximation (Kaiser 1998), although, these approximations start to fail for small ℓ-modes (corresponding to large scales). Various approximations and corrections are investigated in Kilbinger et al. (2017). Their “hybrid” case which is used in most of the recent cosmic shear data analysis (see also Loverde & Afshordi 2008), can be written for redshift bins i and j as follows,

$$\begin{array}{c}\hfill P_{\gamma }^{\mathit{ij}}(\ell )=\frac{9H_0^4{\mathrm{\Omega }}_{\mathrm{m}}^2}{4c^4}{\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi \frac{g^i(\chi )g^j(\chi )}{a^2}{P}_{\delta }(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},\chi ),}\end{array}$$(5)

where H₀ is the Hubble constant, Ω_m is the matter density parameter, c is the speed of light in vacuum, a is the scale factor normalized to one at the present, P_δ is the 3D matter power spectrum and χ is the comoving radial coordinate. The geometric factor for redshift bin i, gⁱ(χ), is given by

$$\begin{array}{c}\hfill g^i(\chi )={\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{h}}}\mathrm{d}{\chi }^{\prime }{p}_{\chi }^i({\chi }^{\prime })\frac{f_{\mathrm{K}}({\chi }^{\prime }-\chi )}{f_{\mathrm{K}}({\chi }^{\prime })},}\end{array}$$(6)

where χ_h is the comoving horizon scale, $p_{\chi }^i(\chi )$ is the probability density of sources in comoving distance for redshift bin i and f_K(χ) is the comoving angular diameter distance, which is equal to χ for a Universe with flat spatial geometry.

The correlation functions calculated using Eq. (4) need to be binned in θ before they can be compared to the measurement. As we usually compress the data by binning ξ_± into broad θ-bins, we should apply the same binning to the theory, to take the functional form of the 2PCFs over the angular bin into account. Additionally, the number of galaxy pairs is roughly proportional to their angular separation. Therefore in the binned data, 2PCFs values for larger θ contribute a larger weight to the mean signal in the bin. In Appendix A we calculate the biases introduced by binning 2PCFs data, showing that using a point estimate for the expected values of ξ_± can produce biases of up to ±10% for the angular range and binning adopted in Hildebrandt et al. (2017) and Troxel et al. (2018b).

In practice we need to modify the relation between the 2PCFs and shear power spectrum in Eq. (4), to accommodate any B-mode power spectra that may exist in the data,

$$\begin{array}{cc}\hfill {\xi }_{+}(\theta )& ={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{J}_0(\theta \ell )[P_{\mathrm{E}}(\ell )+P_B(\ell )],}\hfill \\ \hfill {\xi }_{-}(\theta )& ={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{J}_4(\theta \ell )[P_{\mathrm{E}}(\ell )-P_B(\ell )],}\hfill \end{array}$$(7)

where²P_γ(ℓ)=P_E(ℓ) and P_B(ℓ) is the B-mode power spectrum (Schneider et al. 2002a). In the following subsections we introduce three methods that separate E/B-modes in cosmic shear data: ξ_E/B, COSEBIs and compressed COSEBIs. COSEBIs, being the most robust two-point statistic method, will be used as our primary E/B-mode separation method.

2.1. E/B-mode 2PCFs

The correlation functions, ξ_±, can be separated into E/B-modes following Crittenden et al. (2002) and Schneider et al. (2002a), where

$$\begin{array}{c}\hfill {\xi }_{\mathrm{E}}=\frac{{\xi }_{+}(\theta )+\xi \prime (\theta )}{2}\mathrm{and}{\xi }_B=\frac{{\xi }_{+}(\theta )-\xi \prime (\theta )}{2},\end{array}$$(8)

with

$$\begin{array}{c}\hfill \xi \prime (\theta )={\xi }_{-}(\theta )+4{\displaystyle {\int }_{\theta }^{\infty }\frac{\mathrm{d}\vartheta }{\vartheta }{\xi }_{-}(\vartheta )-12{\theta }^2{\int }_{\theta }^{\infty }\frac{\mathrm{d}\vartheta }{{\vartheta }^3}{\xi }_{-}(\vartheta ).}\end{array}$$(9)

The above definition makes ξ_E/B pure E/B-modes and hence we can write

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{E}}(\theta )& ={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{J}_0(\theta \ell )P_{\mathrm{E}}(\ell ),}\hfill \\ \hfill {\xi }_B(\theta )& ={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{J}_0(\theta \ell )P_B(\ell ).}\hfill \end{array}$$(10)

From Eqs. (7)–(10) we can immediately see that for a B-mode free case ξ_E = ξ₊(θ).

In Schneider et al. (2002a), ξ_E/B(θ) is denoted ξ_E + /B+(θ) as they also provide an alternative definition for E/B two point correlation functions, ξ_{E − /B−}(θ), in terms of integrals over ξ₊(ϑ). In that case the integrals are taken from ϑ = 0 up to ϑ = θ, instead. Although in both cases the integral is taken over a range of angular separations that are not observable, it is preferable to use Eq. (8) since, at least for a B-mode free case, ξ₋(θ)/θ is very small for large θ (ξ₋(θ)∝θ⁻³ at large scales). In this case we can truncate the integrals in Eq. (9) without needing to extrapolate to infinitely large ϑ. However, we may lose some B-mode information by this truncation, as there is no guarantee that the B-mode signal is negligible for large angular scales. One way to extend the integral to large angular scales that are not available in the data is to use the theoretical value of ξ₋(θ) for these angular ranges. In this paper we use measurements over an angular range of [0.5′,300′] and a theoretical ξ₋(θ) from θ = 300′ out to θ = 1000′. We find that the inclusion of the theoretical extension of the integral has less than 5% effect on the largest angular bin (used in KiDS-450) centred at 50′ and drops to subpercent level for θ ≲ 20′.

2.2. COSEBIs

COSEBIs (Complete Orthogonal Sets of E/B-Integrals) modes live neither in Fourier nor real space. The filter functions for COSEBIs form sets of basis functions which transform 2PCFs and shear power spectra to the COSEBIs modes. The two sets of COSEBIs basis functions are the Lin- and Log-COSEBIs filters, which are written in terms of polynomials in ϑ and ln(ϑ) in real space, respectively. In this analysis we use the Log-COSEBIs, as they require fewer modes compared to the Lin-COSEBIs to capture essentially all the cosmological information (see Schneider et al. 2010 for a single redshift bin and Asgari et al. 2012 for the tomographic case).

The COSEBIs can be written in terms of the 2PCFs as

$$\begin{array}{c}\hfill E_n^{(ij)}=\frac{1}{2}{\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\mathrm{d}\vartheta \vartheta [T_{+n}(\vartheta ){\xi }_{+}^{(ij)}(\vartheta )+T_{-n}(\vartheta ){\xi }_{-}^{(ij)}(\vartheta )],}\end{array}$$(11)

$$\begin{array}{c}\hfill B_n^{(ij)}=\frac{1}{2}{\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\mathrm{d}\vartheta \vartheta [T_{+n}(\vartheta ){\xi }_{+}^{(ij)}(\vartheta )-T_{-n}(\vartheta ){\xi }_{-}^{(ij)}(\vartheta )],}\end{array}$$(12)

where $E_n^{(ij)}$ and $B_n^{(ij)}$ are the E and B-mode COSEBIs for redshift bins i and j, and n, a natural number, is the order of the COSEBIs modes. T_±n(ϑ) are the COSEBIs filter functions, (given in Eqs. (28)–(37) in Schneider et al. 2010). These are oscillatory functions with n + 1 roots in their range of support, as shown in Fig. 1. Therefore, the COSEBIs modes with larger n values have more oscillations in their range of support and can pick up features in the 2PCFs that appear as smaller-scale variations, compared to the modes with small n and few oscillations. These small n-modes are more sensitive to the larger-scale variations in the input ξ_± or the overall behaviour of these functions.

Fig. 1. Log-COSEBIs filter functions, T±n(θ). These filter functions convert ξ± to COSEBIs E and B modes through Eqs. (11) and (12). We show four example n-modes for each filter for the angular separation range of [0.5′,100′]. By definition T±n(θ) are equal to zero outside of the range of their support.

The E/B-COSEBIs can also be expressed as a function of the convergence power spectra,

$$\begin{array}{cc}\hfill E_n^{(ij)}& ={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{P}_{\mathrm{E}}^{(ij)}(\ell )W_n(\ell ),}\hfill \\ \hfill B_n^{(ij)}& ={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{P}_{\mathrm{B}}^{(ij)}(\ell )W_n(\ell ),}\hfill \end{array}$$(13)

where $P_{E(B)}^{(ij)}$ are the E(B)-mode convergence power spectra and the W_n(ℓ) are the Hankel transforms of T_±n(ϑ),

$$\begin{array}{cc}\hfill W_n(\ell )& ={\displaystyle {\int }_{{\vartheta }_{\min }}^{{\vartheta }_{\max }}\mathrm{d}\vartheta \vartheta T_{+n}(\vartheta )J_0(\ell \vartheta ),}\hfill \\ \hfill & ={\displaystyle {\int }_{{\vartheta }_{\min }}^{{\vartheta }_{\max }}\mathrm{d}\vartheta \vartheta T_{-n}(\vartheta )J_4(\ell \vartheta ).}\hfill \end{array}$$(14)

Figure 2 shows the W_n(ℓ) functions corresponding to the T_±n(θ) filters shown in Fig. 1. The first peak in W_n(ℓ) is set by the value of ϑ_max and n. As can be seen, the higher order W_n pick up more power from larger ℓ. We use Eq. (13) to calculate the theoretical value of the E-mode COSEBIs as theories, in general, give their predictions in terms of the power spectrum. However, in practice the shear 2PCFs are more straightforward to measure from data, hence Eqs. (11) and (12) are used to calculate the E/B-mode COSEBIs from data and simulations. To evaluate these integrals in the angular range of [0.5′,100′] we use 4 × 10⁵ linear angular bins (see Asgari et al. 2017, for a discussion on optimising the number of bins for this type of analysis).

Fig. 2. Log-COSEBIs weight functions, Wn(ℓ), normalized to their maximum value. These weight functions convert E and B shear power spectra to COSEBIs modes through Eq. (13). Four example n-modes are shown for the angular range of [0.5′,100′].

2.3. CCOSEBIs

We use the data compression method of Asgari & Schneider (2015) to explore the effect of systematics on cosmological parameter estimation, as this method is informed by the sensitivity of the data to the parameters. This method, which can be applied to any statistic, reduces the number of data points, which is important to minimise errors when estimating covariance matrices from simulations (see Hartlap et al. 2007, for example).

To compress COSEBIs we need to have an estimate for their inverse covariance matrix (see Sect. 2.4), as well as their first and second-order derivatives with respect to the parameters to be measured. We then linearly combine the COSEBIs modes using the sensitivity of each mode to the given parameter(s) as their coefficient. For the first-order compressed E-COSEBIs we have,

$$\begin{array}{c}\hfill E_{\mu }^{\mathrm{c}}={\displaystyle \sum _{n,m=1}^{n_{\max }}\frac{\partial E_m}{\partial \mu }{({\mathbb{C}}^{-1})}_{\mathit{mn}}{E}_n,}\end{array}$$(15)

where μ is a cosmological parameter, ℂ⁻¹ is the inverse covariance matrix of E_n and n_max is the number of COSEBIs modes considered in the compression. This first order compression is equivalent to a Karhunen–Loeve compression where the covariance matrix is known (see Tegmark et al. 1997, for example), but using the first order compression alone can result in a loss of information when the covariance matrix estimate is inaccurate. We therefore follow Asgari & Schneider (2015) by adding the following second-order compressed quantities to the data,

$$\begin{array}{c}\hfill E_{\mu \nu }^{\mathrm{c}}={\displaystyle \sum _{n,m=1}^{n_{\max }}\frac{{\partial }^2{E}_m}{\partial \mu \partial \nu }{({\mathbb{C}}^{-1})}_{\mathit{mn}}{E}_n,}\end{array}$$(16)

where ν is a second cosmological parameter and second order derivatives of E_n are taken. In short, we can write both first and second-order CCOSEBIs as the following matrix equation,

$$\begin{array}{c}\hfill {\mathit{E}}^{\mathrm{c}}=\mathbf{\Gamma }\mathit{E}\mathrm{and}{\mathit{B}}^c=\mathbf{\Gamma }\mathit{B},\end{array}$$(17)

where the elements of the compression matrix, Γ, are formed from combinations of the derivatives of E_n with respect to the parameters and their inverse covariance matrix.

2.4. Covariance matrix

To quantify the significance of the B-modes measured from the data or in the simulations, we need to know the covariance matrix of the data vector. Aside from currently negligible physical effects that can produce B-modes (discussed in Sect. 1) and in the absence of systematics, we expect any observed B-modes to be consistent with random noise arising from galaxy shape-noise. Therefore, to calculate the B-mode covariance, we assume that they are only due to noise and find the covariance matrix for each of the above statistics.

Assuming a field of galaxies with ellipticities randomly picked from a Gaussian with zero mean and ${\sigma }_{ϵ}^2$ variance, we can write the covariance of the 2PCF as

$$\begin{array}{c}\hfill ⟨{\xi }_{\pm }^{ij,\mathrm{n}oise}(\theta ){\xi }_{\pm }^{kl,\mathrm{n}oise}(\vartheta )⟩=\frac{{\sigma }_{ϵ}^4}{2N_{\mathrm{pair}}^{\mathit{ij}}(\theta )}{\delta }_{\theta \vartheta }[{\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}+{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}}],\end{array}$$(18)

where $N_{\mathrm{pair}}^{\mathit{ij}}(\theta )$ is the number of galaxy pairs, in redshift bin pair ij, within an angular separation bin with the label θ. The Kronecker symbols, δ_ij and δ_θϑ, are equal to unity if their arguments are equal and are otherwise zero (for example see Eq. (34) in Joachimi et al. 2008).

An approximation for N_pair(θ) can be determined by calculating the number of pairs in an infinite field, scaled by the true finite field area, A, where

$$\begin{array}{c}\hfill N_{\mathrm{pair}}^{\mathit{ij}}{(\theta )}^{\mathrm{approx}}=2\pi A\theta \mathrm{\Delta }\theta {\overline{n}}_{\mathrm{gal}}^i{\overline{n}}_{\mathrm{gal}}^j,\end{array}$$(19)

and ${\overline{n}}_{\mathrm{gal}}^i$ is the mean number density of galaxies in redshift bin i. This approximation fails, however, as it does not account for intricate small-scale survey geometry, source clustering or any variable depth effects. Furthermore, as we get closer to the field size, it does not account for the pairs of galaxies which are lost due to the discontinuities in the observed field (see for example Joachimi et al. 2008). As the significance of any measured B-modes is determined entirely by the shot-noise, we therefore choose to use a direct measurement of N_pair(θ) from the data. We follow the method of Schneider et al. (2002b), who determine the full covariance matrix for 2PCFs for a weighted ellipticity field, to find the shape-noise-only term of the covariance matrix, with the number of galaxy pairs given as

$$\begin{array}{c}\hfill N_{\mathrm{pair}}(\theta )=\frac{{({\sum }_{\mathit{ab}}{w}_a{w}_b)}^2}{{\sum }_{\mathit{ab}}{w}_a^2{w}_b^2}·\end{array}$$(20)

Here the sums are over galaxies in the given angular separation bin. Determining N_pair from Eq. (20) instead of the approximation in Eq. (19), enlarges the covariance at large scales where there are fewer pairs of galaxies due to geometry effects. On small scales where variable depth and source clustering become important, the covariance is decreased.

Inserting Eq. (18) into the following expression for the COSEBIs covariance (Schneider et al. 2010)

$$\begin{array}{cc}\hfill C_{\mathit{mn}}^{ij,kl}& =\frac{1}{4}{\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\mathrm{d}\theta \theta {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\mathrm{d}\theta \prime \theta \prime }\hfill \\ \hfill & \times {\displaystyle \sum _{\mu \nu =+,-}{T}_{\mu m}(\theta )T_{\nu n}(\theta \prime )C_{\mu \nu }^{ij,kl}(\theta ,\theta \prime ),}\hfill \end{array}$$(21)

where C_±±(θ, θ′) is the covariance of ξ_±, we find the B-mode covariance for COSEBIs,

$$\begin{array}{cc}\hfill C_{\mathit{mn}}^{ij,kl}& =\frac{{\sigma }_{ϵ}^4}{8}{\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\mathrm{d}\theta \theta {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\frac{\mathrm{d}\theta \prime \theta \prime }{n_{\mathrm{pair}}^{\mathit{ij}}(\theta \prime )}[{\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}+{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}}]}\hfill \\ \hfill & \times {\delta }_{\mathrm{D}}(\theta -\theta \prime )[T_{+m}(\theta )T_{+n}(\theta \prime )+T_{-m}(\theta )T_{-n}(\theta \prime )],\hfill \end{array}$$(22)

where $n_{\mathrm{pair}}^{\mathit{ij}}(\theta )\mathrm{d}\theta =N_{\mathrm{pair}}^{\mathit{ij}}(\theta )$, δ_D is the Dirac delta function and we have used δ_θθ′ = δ_D(θ − θ′) Δθ to remove the Kronecker symbol. Taking the inner integral in Eq. (22) results in,

$$\begin{array}{cc}\hfill C_{\mathit{mn}}^{ij,kl}& =\frac{{\sigma }_{ϵ}^4}{8}[{\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}+{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}}]{\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\frac{\mathrm{d}\theta {\theta }^2}{n_{\mathrm{pair}}^{\mathit{ij}}(\theta )}}\hfill \\ \hfill & \times [T_{+m}(\theta )T_{+n}(\theta )+T_{-m}(\theta )T_{-n}(\theta )].\hfill \end{array}$$(23)

We calculate the COSEBIs B-mode covariance using trapezoidal integration with fine θ-bins and verified that these equations accurately predict the noise-only covariance, by analysing 1000 shape-noise-only mock simulations. We find that for a 100 deg² field the noise term for the COSEBIs covariance is underestimated by 30% if we use $N_{\mathrm{pair}}^{\mathrm{approx}}$ from Eq. (19), while using Eq. (20) recovers the measured covariance from the simulation.

The corresponding covariance for CCOSEBIs is simply equal to the COSEBIs covariance sandwiched between two compression matrices,

$$\begin{array}{c}\hfill {\mathbb{C}}^{\mathrm{c}}=\mathbf{\Gamma }\mathbb{C}{\mathbf{\Gamma }}^{\mathrm{t}},\end{array}$$(24)

where ^t denotes a transposed matrix³.

The covariance matrix of ξ_B can also be calculated from Eq. (18),

$$\begin{array}{c}\hfill ⟨{\xi }_B^{\mathit{ij}}(\theta ){\xi }_B^{\mathit{kl}}(\vartheta )⟩=\frac{{\sigma }_{ϵ}^4}{4N_{\mathrm{pair}}^{\mathit{ij}}(\theta )}{\delta }_{\theta \vartheta }[{\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}+{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}}].\end{array}$$(25)

Note that the only difference between Eqs. (25) and (18) is a factor of 2, which arises from the fact that ξ_± depends on both E/B-modes and their associated noise, while ξ_E/B only depends on a single component, as can be seen in Eqs. (10) and (7). As a result, ξ_E/B is only sensitive to the noise components that resemble E/B-modes. The power spectrum of the noise can be equally divided into an E-mode and a B-mode component, and as such the noise covariance for ξ_E/B is half the amplitude of the corresponding covariance for 2PCFs.

In addition to B-modes, we show E-mode measurements for the data with error bars calculated assuming Gaussian covariances. We choose not to include the non-Gaussian and super sample terms in the error calculation which primarily affect the off-diagonal terms of the covariance matrix⁴. As we do not analyse the E-modes in a quantitative way in this study, and use the E-mode covariances solely for plotting purposes, our chosen Gaussian treatment of the covariance is sufficient. We can write the Gaussian covariance for the E-modes in terms of three contributors,

$$\begin{array}{c}\hfill \mathbb{C}=\mathrm{Cosmic}\mathrm{variance}+\mathrm{Mixed}+\mathrm{Noise},\end{array}$$(26)

where the Mixed term depends on both cosmology and noise. The Noise term here is estimated in the same manner as the B-modes covariance, (Eqs. (23)–(25)), taking all the survey effects into account. For the other two contributions, however, we assume a simple survey geometry and follow Eqs. (53) and (54) in Joachimi et al. (2008) for the covariance of power spectra and correlation functions⁵, respectively. The Gaussian mixed and cosmic variance terms for COSEBIs covariance are given in Eq. (11) in Asgari et al. (2012) for the tomographic case.

3. Data

We use three sets of cosmic shear catalogues that are in the public domain, KiDS-450, DES-SV and CFHTLenS. Our focus in this paper is the analysis of their B-mode signal, but we also compare the corresponding measured E-mode signals to theoretical predictions, based on the published best fitting cosmological parameters from each survey, as given in Table 1. This allows for the level of B-modes to be assessed, relative to the E-modes, but we leave a full E-mode cosmological parameter analysis to a future paper.

The theoretical predictions are calculated using COSMOSIS (Zuntz et al. 2015)⁶ with linear matter power spectra calculated with CAMB (Lewis et al. 2000; Howlett et al. 2012)⁷. Takahashi et al. (2012) is used to model the non-linear evolution of the matter power spectrum. A Limber approximation is employed to estimate the lensing power spectrum as described in Sect. 2. For the intrinsic alignment of galaxies we adopt the non-linear model from Bridle & King (2007)⁸, which is equivalent to the models used in the analysis of all three surveys. The 2PCFs are measured from the data and the simulations using ATHENA⁹ (Kilbinger et al. 2014).

3.1. CFHTLenS

Heymans et al. (2012) present the CFHTLenS, a completed survey with 154 square degrees of observed data in 5 photometric bands. The public data products that we analyse here are processed by THELI (Erben et al. 2013), with galaxy ellipticities measured using lensfit (Miller et al. 2013) and photometric redshifts determined using the Bayesian photometric redshift code BPZ (Benítez 2000; Hildebrandt et al. 2012).

The 2PCFs cosmic shear analysis for CFHTLenS is presented in Kilbinger et al. (2013) and Heymans et al. (2013). As summarised in Kilbinger et al. (2017), however, several improvements have been recognised since these publications, in particular with respect to the calibration of the photometric redshifts (see for example Choi et al. 2016; Joudaki et al. 2017b) and the shear measurements (see Kuijken et al. 2015; Fenech Conti et al. 2017). The resulting uncertainty in these calibrations will impact the E-mode cosmological parameter constraints from this survey. As our focus is on a B-mode analysis however, which is independent of these calibration corrections, we choose to use the redshift distributions and calibration corrections adopted by Heymans et al. (2013) for this study.

We follow Heymans et al. (2013) by dividing the data into six photometric redshift bins: z₁ ∈ (0.2, 0.39], z₂ ∈ (0.39, 0.58], z₃ ∈ (0.58, 0.72], z₄ ∈ (0.72, 0.86], z₅ ∈ (0.86, 1.02] and z₆ ∈ (1.02, 1.3], also including a single bin case that uses the full range of z ∈ (0.2, 1.3]. In Asgari et al. (2017), we analysed CFHTLenS using COSEBIs to find a significant level of B-modes. We extend this analysis to explore higher modes in COSEBIs, in addition to ξ_E/B, and we use an exact noise covariance (Eq. (20)) in contrast to our earlier work which used Eq. (19).

3.2. KiDS-450

The Kilo-Degree Survey (KiDS) will collect 1350 square degrees and in combination with VIKING (VISTA Kilo-degree Infrared Galaxy survey) will present data in nine photometric bands (see Kuijken et al. 2015 and de Jong et al. 2017). We analyse the data products released for the first 450 square degrees (KiDS-450), that has been processed by THELI (Erben et al. 2013) and Astro-WISE (Begeman et al. 2013). Galaxy ellipticities are measured with lensfit (Miller et al. 2013) and calibrated using the image simulations described in Fenech Conti et al. (2017). The 4-band photometric redshifts are calibrated using external overlapping spectroscopic surveys (Hildebrandt et al. 2017) and galaxies are binned into tomographic bins using BPZ.

The KiDS-450 2PCFs cosmic shear analysis is shown in Hildebrandt et al. (2017) and Joudaki et al. (2017a); Joudaki et al. (2018), with complementary cosmic shear power spectrum analyses calculated using quadratic estimators in Köhlinger et al. (2017), and integrals over 2PCFs in van Uitert et al. (2018). All these analyses reported significant but low-level traces of B-modes in the data.

As in the KiDS-450 cosmic shear analyses we divide the data into four photometric redshift bins: z₁ ∈ (0.1, 0.3], z₂ ∈ (0.3, 0.5], z₃ ∈ (0.5, 0.7] and z₄ ∈ (0.7, 0.9], including a single bin case that uses the full range of z ∈ (0.1, 0.9].

3.3. DES-SV

The Dark Energy Survey Collaboration (2005) introduce the Dark Energy Survey (DES) project which will produce 5000 square degrees of gravitational lensing data in five bands. The science verification data also known as DES-SV¹⁰ is the public dataset that we analyse here. The galaxy ellipticities in DES-SV are measured using NGMIX (Jarvis et al. 2016) and photometric redshifts are determined using a machine learning-based pipeline, SKYNET (Bonnett et al. 2016).

Becker et al. (2016) present the primary cosmic shear analysis of the DES-SV data using 2PCFs along with cosmic shear power spectrum measurements (also see Troxel et al. 2018a, for the analysis of the first 1300 square degrees of DES data). Fourier space B-mode measurements detected no significant B-modes on scales ℓ <  2500.

We divide the data into three photometric redshift bins following Becker et al. (2016): z₁ ∈ (0.3, 0.55), z₂ ∈ (0.55, 0.83) and z₃ ∈ (0.83, 1.3) and also consider a single bin case that uses the full range of z ∈ (0.3, 1.3). In order to compare our measured E-mode signal to the published best-fitting cosmological parameters, listed in Table 1, we also take into consideration the best-fitting DES-SV shear calibration and photometric redshift biases in our predictions, which Abbott et al. (2016) include as nuisance parameters in their fit. For our single bin analysis of DES-SV data we adopt zero bias for the photometric redshift and the same value as the first tomographic bin for the shear calibration bias, which is similar to the average of the biases measured for the three bins (see Table D.1).

4. Results: survey E/B-modes

In this section we present the measured COSEBIs (Eq. (11) and Eq. (12)), CCOSEBIs (Eq. (17)) and ξ_E/B (Eq. (8)) for KiDS-450, DES-SV and CFHTLenS. In Fig. 3 we show the COSEBIs measurement for a single redshift bin encompassing the full range of redshifts adopted by each survey. For the COSEBIs statistics we need to choose an angular range and throughout this paper we show results for three sets of angular ranges: the full angular range: [0.5′,100′], large scales: [40′,100′] and small scales: [0.5′,40′]. These were chosen to span both the survey-adopted ξ₊(θ) angular ranges: KiDS-450 (0.5′< θ₊ <  72′), DES-SV (2′≲θ₊ ≲ 60′), CFHTLenS (1.5′< θ₊ <  35′), whilst also probing some of the larger angular scales used in the corresponding ξ₋(θ) analysis: KiDS-450 (8.6′< θ₋ <  300′), DES-SV (24.5′≲θ₋ ≲ 245.5′) and CFHTLenS (1.5′< θ₋ <  35′). The large scale cut for ξ₊(θ) is generally employed to avoid biasing the results, when a constant additive bias term (c-term) is present in the shear catalogues. The same large scale angular-cut is not applied to ξ₋(θ), since this statistic is not sensitive to a constant c-term. COSEBIs share this insensitivity with ξ₋(θ) and hence any measured COSEBIs B-modes that use scales beyond the maximum θ₊ range are not a result of a constant c-term.

Fig. 3. COSEBIs E-modes (left) and B-modes (right) for a single broad redshift bin. Results for DES-SV are shown with blue squares, KiDS-450 with black stars and CFHTLenS with magenta triangles. The angular ranges are shown for each row in the upper right corner. In addition, the significance of the B-modes is shown as p-values for each survey and angular range. E-mode predictions are calculated using the best fitting cosmological parameter values given in Table 1 for DES-SV (solid), KiDS-450 (dashed) and CFHTLenS (dotted). Note that COSEBIs modes are discrete and the theory values are connected to each other only as a visual aid. A zero-line is also shown for reference.

Each row in Fig. 3 corresponds to one angular range, as denoted in the right panels, with E-modes on the left and B-modes on the right. The different symbols show the results for DES-SV (squares), KiDS-450 (stars) and CFHTLenS (triangles). Overlaid are the theoretical predictions, given the published best-fitting survey cosmological parameters from Table 1. We show these E-mode predictions as curves for ease of comparison, even though the COSEBIs modes are discrete. As the COSEBIs modes are correlated to their neighbouring modes (see Asgari et al. 2012, 2017, for plots of the covariance matrices), we caution that the goodness-of-fit to the model should not be deduced by simply looking at the graphs, a practice commonly known as “χ-by-eye”. Any goodness-of-fit exercise must take into account the significant correlations between the points.

Throughout this paper we truncate the COSEBIs measurements at n = 20. In principle the COSEBIs can be estimated for an infinite number of modes. However as can be seen in the left hand panels of Fig. 3 the E-mode predictions are equivalent to zero for n ≳ 7. Therefore, we do not expect to gain any cosmological information from these modes. On the other hand the signal from systematic effects does not necessarily follow the same behaviour. We expect a significant signal at larger n-modes for certain systematics (see Sect. 6). As a result we choose n = 20 as our maximum n-mode, which encompasses modes that are important for both cosmological and systematic analyses. Future analyses may however want to extend their diagnostic B-mode analysis to even higher n-modes, depending on the signature of the systematic that they are searching for.

Focusing first on the E-mode measurements (left panels of Fig. 3) we expect to measure signal in the lower n-modes and none for the modes n ≳ 8, as seen in the theoretical predictions. This arises from the fact that both the 2PCFs and shear power spectrum are relatively smooth functions with a few features that are captured, almost entirely, by the first few COSEBIs modes. Any significant detection of high-order COSEBIs modes indicates high-frequency variations in the 2PCFs, which are unexpected in a ΛCDM cosmology and therefore indicative of systematics. We remind the reader that our E-mode errors, which include both sampling variance assuming a Gaussian shear field and shot noise, will be slightly underestimated as we have not included the sub-dominant super sample and non-Gaussian contributions to the sampling variance terms (see Sect. 2.4).

Turning to the B-mode measurements (right panels of Fig. 3), we determine the significance of the measured B-modes using “p-values”, for each dataset and angular range, listed in Table 2. The p-value is equal to the probability of randomly producing a B-mode that is equally or more significant than the measured B-mode signal, given the model that B-modes are equal to zero and their distribution is Gaussian (see Appendix C for the mathematical definition of p-value). This model is appropriate for B-modes generated from random noise. The degrees-of-freedom here is equal to the number of COSEBIs modes (20 modes in the single redshift bin case), as the model has no free parameters to be fitted. The p-values take into account the correlations between the COSEBIs modes. Our error analysis for the B-modes is accurate, taking into account the weighted number of galaxy pairs in each dataset. We consider the B-modes to be significant when the measured p-values are p <  0.01 (highlighted in bold), corresponding to greater than 2.3σ detection of B-modes. We find that the B-modes of KiDS-450 and CFHTLenS are consistent with zero, finding p >  0.1 in all cases. DES-SV, however, shows significant 2.8σB-modes with p = 0.0026, when the full angular range is considered. We discuss the complexity of linking B-mode features with E-mode features in Sect. 6.

In Table 2 we also list the significance of the measured COSEBIs B-modes for a tomographic analysis of the three angular ranges, using the survey-defined photometric redshift bins (see Sect. 3). The COSEBIs tomographic measurements for each survey, adopting the full angular range, are shown in Appendix D. For all angular ranges, we find no significant COSEBIs B-modes for KiDS-450. In contrast, for DES-SV data we find a 4.0σ detection of B-modes for the large-scale angular range that includes angular scales used in the DES-SV cosmic shear analysis. For the full angular range, including small-scale information that was excluded from the DES-SV cosmic shear analysis, the significance of the detection increases to 5.5σ. For CFHTLenS we find a significant B-mode detection for small scales, but not at large scales. This result is in contrast to Asgari et al. (2017) who found significant CFHTLenS B-modes for large, but not small scales. We do however recover this result if we limit our p-value analysis to the first 7 COSEBIs modes adopted by Asgari et al. (2017). This demonstrates that the p-values are sensitive to the choice of modes considered in the analysis, motivating the study of how different systematics impact different COSEBIs E/B-modes in Sect. 6.

In Fig. 4 we show the measured compressed COSEBIs, where the COSEBIs modes are combined to produce a set of E-mode CCOSEBIs that, in a systematic-free dataset, are only sensitive to cosmological parameters (Eq. (17)). We compress the B-mode COSEBIs using the same compression matrix. Cosmic shear is mainly sensitive to a combination of σ₈ and Ω_m (see for example Jain & Seljak 1997), hence we choose only these two parameters to form the CCOSEBIs modes¹¹. The CCOSEBIs modes are highly correlated as σ₈ and Ω_m are degenerate in cosmic shear data, and we hence caution the reader, again, against a “χ-by-eye” analysis.

Fig. 4. CCOSEBIs E and B-modes for non-tomographic (left) and tomographic (right) analyses. The E-modes are shown as empty symbols, with the B-modes shown as filled symbols, for DES-SV (blue squares), KiDS-450 (black stars) and CFHTLenS (magenta triangles). The analysis is conducted over three different angular ranges, denoted in the upper right corner of each panel. The CCOSEBIs mode is indicated on the horizontal-axis. E-mode predictions are calculated using the best fitting cosmological parameter values given in Table 1 for DES-SV (solid), KiDS-450 (dashed) and CFHTLenS (dotted). A zero-line is also shown for reference.

Figure 4 shows the results for a single-bin analysis (left) and a tomographic analysis (right) for the three sets of angular ranges indicated in the right panels. The E-modes are shown as open symbols and the B-modes as filled symbols for DES-SV, KiDS-450 and CFHTLenS. Overlaid is the theoretical expectation for the E-mode signal, shown as curves for visual aid even though the CCOSEBIs modes are discrete. The horizontal axis shows which parameter (for the first-order modes) or two parameters (for the second-order modes) the CCOSEBIs mode is sensitive to. We highlight that CCOSEBIs represent a significant data compression, particularly in the tomographic case where, for example, we compress the 3-bin 120 data-point DES-SV analysis, and the 6-bin 420 data point CFHTLenS analysis, down to the same 5 CCOSEBIs modes.

Comparing the measured E-modes with the level of B-modes in Fig. 4 we find that, aside from the largest angular range that also has the lowest signal-to-noise ratio, the E-modes are about an order of magnitude larger than the B-modes. In all panels we see that the KiDS-450 E-mode signal is lower than DES-SV and CFHTLenS, resulting from a smaller upper photometric redshift cut of z_phot <  0.9 in this dataset.

Table 3 shows the p-values for CCOSEBIs B-modes. We do not show p-values for the E-modes since for this analysis we have not included the super sample covariance term and our E-mode errors are therefore underestimated. Readers concerned by the apparent offset between the highly correlated [40′,100′] E-mode measurements and expectation values, should note that the similarly highly correlated B-modes, for the non-tomographic case, are all consistent with zero, even for the cases where they look inconsistent. The significance of the B-modes is different from the values shown in Table 2, where we have used the first 20 COSEBIs modes to measure the p-values. This apparent inconsistency is not unexpected, as the bulk of the CCOSEBIs signal comes from the first few COSEBIs modes, which contain the cosmological signal and different levels of systematics in comparison to the full set of 20 COSEBIs modes. A good example of this difference comes in the tomographic analysis of DES-SV where we find a significant ∼5.5σ non-zero B-mode signal for COSEBIs, but the CCOSEBI B-mode is not significant at 2.2σ. This shows that the first few COSEBIs modes have a smaller contribution to the total DES-SV B-mode signal compared to the higher order modes, which can also be seen in Fig. D.1. KiDS-450, however, shows an insignificant B-mode signal when we consider both high and low COSEBIs modes, in contrast to a 2.7σB-mode detection when only the low n-modes are used to construct the CCOSEBIs. As can be seen in Table 3, and the upper right panel of Fig. 3, the low-nB-modes for KiDS-450 data only become significant when the small angular scales are included.

If the origin of the B-modes detected in the COSEBIs analysis was known to impact the E and B modes equally, then the CCOSEBIs result would be the most relevant for cosmic shear studies. If the systematics impact the E and B modes differently, however, then the compressed CCOSEBIs result, focused on only low-n modes, could lead to a false null-test for the survey. It is therefore important to study how different systematics impact the full range of E and B COSEBIs, which we carry out in Sect. 6, and discuss this matter further in Sect. 7. In Appendix C we also discuss how analysis choices, for example in this case tomographic or non-tomographic, COSEBIs or CCOSEBIs, can optimise or dilute the power of a B-mode null-test.

Finally we turn to Fig. 5 which shows the measured ξ_E/B statistic across the full redshift range for each survey. Overlaid are the best fitting theory curves for each dataset derived from the published cosmological parameters in Table 1. The p-values corresponding to the zero B-mode model are low in all cases, as given in the legend of the figure, with all surveys showing a tendency for increasing B-mode power and decreasing E-mode power at large scales, which, we discuss further in Sect. 7. B-modes are detected at greater than ∼2.6σ for all surveys. For DES-SV the significance of the B-modes is particularly high at ∼9σ, but this reduces to 2.3σ, or p = 0.012, when we select the angular scales [4.2′,72′] which roughly correspond to the angular cuts applied to ξ₊ in the DES-SV cosmic shear analysis¹². In Appendix D we present the ξ_E/B tomographic analysis for each survey where we find a significant B-mode detection for DES-SV (p ∼ 4 × 10⁻¹⁹), but no detection of B-modes for KiDS or CFHTLenS (p ∼ 0.7).

Fig. 5. ξE and ξBE/B-modes for a single broad redshift bin. The E-modes are shown as empty symbols, with the B-modes shown as filled symbols, for DES-SV (blue squares), KiDS-450 (black stars) and CFHTLenS (magenta triangles). The DES-SV and CFHTLenS results are horizontally offset relative to KiDS-450 to aid visualisation. E-mode predictions for ξE are calculated using the best fitting cosmological parameter values given in Table 1 for DES-SV (solid), KiDS-450 (dashed) and CFHTLenS (dotted). A zero-line is also shown for reference. We detect significant B-modes in all cases as shown by the p-values, in the legend, which determine the probability of the data B-modes given a null B-mode model.

Given the required extrapolation of the data in order to calculate the ξ_E/B statistic (see Eq. (9)) we emphasize that these results are, by nature, a biased measurement of ξ_E/B, which may not represent the data accurately. For this statistic, the errors on ξ_B are uncorrelated (see Eq. (25)) but also biased as the integral truncation when estimating ξ_E/B also affects its noise properties, which we have not taken into account. We therefore do not place too strong an emphasis on the high significance of the measured B-modes, or the lack of E-mode power on large-scales for all surveys, particularly as these are the scales that are most impacted by the choices made when extrapolating the data. That said, if surveys continue to use 2PCFs as a standard cosmic shear statistic, then it is still relevant to measure ξ_E/B as it is the B-mode measurement that is most closely related to the 2PCFs.

5. Modelling systematics

In Sect. 4 we detected significant B-modes in the DES-SV data as well as in certain tomographic combinations of CFHTLenS and KiDS-450 data. In this section we introduce models for a range of data-related systematic effects that are appropriate for the datasets described in Sect. 3. We consider three models of systematics that affect the shear measurement of all galaxies independently of their redshift. In addition we model one photometric redshift-dependent systematic, demonstrating how catastrophic errors in photometric redshifts can lead to shape selection bias. We add these systematic models to mock data to explore their effect on the 2 point statistics introduced in Sect. 2. We are particularly interested in measuring the B-modes associated with each systematic (see Sect. 6). By comparing simulated results with and without these systematic effects, the B-mode signatures can then be used as a tool for diagnosing the source of the B-modes in the surveys analysed in Sect. 4.

In this analysis, we do not model masking effects, since all of the methods we use rely on measuring 2PCFs, which are insensitive to masking, provided the mask is uncorrelated with the shear field. If such correlations exist, all statistics will be affected by them. This is in contrast to methods that rely on Fourier transforms of the shear field, where masks can cause significant systematic effects (see for example Asgari et al. 2018).

5.1. Shear measurement errors

For the case of weak shear with |γ|≪1, we can model the observed ellipticity as

$$\begin{array}{c}\hfill {ϵ}^{\mathrm{obs}}=(1+m)[{ϵ}^{\mathrm{int}}+\gamma ]+\eta +\alpha {ϵ}^{\ast }+\beta \delta {ϵ}^{\ast }+c,\end{array}$$(27)

where ϵ^int is the intrinsic galaxy ellipticity, γ is the shear, η is random noise on the ellipticity measurement, ϵ^* is the PSF model ellipticity, $\delta {ϵ}^{*}={ϵ}^{*}-{ϵ}_{\text{true}}^{*}$ is the residual ellipticity between the model and true PSF, and c is an additive shear that is uncorrelated with the PSF. For all these quantities we use complex notation where, for example, γ = γ₁ + i γ₂. For the two PSF-dependent terms, αϵ^* quantifies the fraction of the model PSF ellipticity that leaks into the shape measurement, and β δϵ^* quantifies the fraction of the residual PSF arising from PSF modelling errors, that leaks into the shape measurement. The term m is a multiplicative shear bias that is traditionally calibrated using image simulations.

We simulate each of the systematic terms in Eq. (27) in isolation, in order to characterise their B-mode signature. One exception is the shear calibration correction, m, which we set to zero, as an isotropic shear bias cannot introduce a B-mode signal, only scale an E-mode signal.

5.1.1. Point spread function (PSF) leakage: αϵ^*

In order to mimic the effect of the PSF leakage on cosmic shear measurements we use PSF models from KiDS to make a realistic spatially varying PSF model spanning 100 square degrees. We construct this large-scale PSF pattern on a 1 arcmin² resolution grid, mapping the KiDS PSF measurements onto a 10 ° ×10° field by stitching together two 5 ° ×10° sections from the G12 and G15 regions in KiDS-450 data (see Hildebrandt et al. 2017, for details). This provides us with a model for ${ϵ}_i^{*}(x,y)$, where ${ϵ}_i^{*}(x,y)$ is shown in the left panel of Fig. 6. In KiDS the PSF is modelled with polynomials of third order within each pointing, where the lowest order is allowed to vary between CCDs to allow for discontinuities between CCD chips (see Kuijken et al. 2015, for more details). Similar modelling approaches are taken by CFHTLenS and DES-SV. The mean of the PSF ellipticity and its one sigma deviation is ${ϵ}_i^{*}=0.006\pm 0.016$ and its full range is covered by $-0.1<{ϵ}_i^{*}<0.1$ for both components. Figure 6 shows how the PSF pattern changes within and across each ∼square degree pointing. In areas where the KiDS data are masked and the PSF model unconstrained, we linearly interpolate the value of the PSF ellipticity to accommodate all galaxy positions in our unmasked mock data analysis.

Fig. 6. First ellipticity component of the spatially varying systematic effects, simulated over a 10 ° ×10° field. Here the effects are normalized to their maximum value for a better visual comparison. From the left, the first panel shows the point spread function pattern used to model PSF leakage ($-0.01<\alpha {ϵ}_1^{*}<0.01$). The second panel shows a regular pattern using the detector chip bias model from OmegaCam multiplied by a factor of 5 (0.001 <  c1 <  0.025). The third panel shows the random correlated noise PSF residual model with a smoothing length similar to the chip size ($-0.006<\beta \delta {ϵ}_1^{*}<0.006$), while the last panel shows the same systematic for a roughly pointing size smoothing length ($-0.003<\beta \delta {ϵ}_1^{*}<0.003$).

We choose to apply a 10% PSF leakage by setting α = 0.1. This level of leakage corresponds to the α measured in the poorer-seeing KiDS i-band data (see Amon et al. 2018). For the high-quality KiDS r-band data that are used for the main cosmic shear analysis, α was found to be consistent with zero (Hildebrandt et al. 2017).

5.1.2. Regular repeating additive pattern: c(x,y)

In the absence of PSF-related errors, the amplitude of any remaining additive bias that is uncorrelated with the PSF, c, can be estimated directly from the data. Since we expect ⟨ϵ^int + γ_i + η⟩=0 over a large area, $⟨{ϵ}_i^{\text{obs}}⟩=c$ when α = β = 0. Surveys typically correct for any significant measurement of c before any analysis, but this empirical correction usually takes an average over all galaxies and is therefore insensitive to small scale spatial variations c(x, y) (van Uitert & Schneider 2016). Systematic effects that are stable and associated with the camera or telescope would result in a repeating pattern in the survey which is built from a series of different pointings. To determine the impact of such a systematic, we model a spatially varying, but repeating additive term, which remains constant between pointings.

We use the data from Hoekstra et al. (in prep.), who present a detailed analysis of imaging from the KiDS OmegaCam camera. Hoekstra et al. (in prep.) report low-level but significant detector and electronic defects that introduce an additive shear contribution per CCD that is uncorrelated with the PSF and spans the range $0.0002<{ϵ}_1^{*}<0.005$ and $-0.0004<{ϵ}_2^{*}<0.00015$, shown in the second panel in Fig. 6. The white “chip 15” of OmegaCam shows the largest bias at the level of ϵ₁ = 0.005. For the purposes of this analysis we multiply the Hoekstra et al. (in prep.) detector-bias model by a factor of five to amplify its effect, as we find that the original level of this effect is too small to show any significant E/B-modes for the current datasets.

5.1.3. Random but correlated noise: β δϵ^*

Errors in PSF modelling, δϵ^*, can be systematic if the stars used in the modelling are unrepresentative of the PSF experienced by the galaxies (Antilogus et al. 2014; Guyonnet et al. 2015; Gruen et al. 2015). In this case the resulting systematic behaves similarly to the PSF leakage model outlined in Sect. 5.1.1, and we therefore do not consider this type of PSF modelling error.

Instead we consider the random errors in the PSF modelling that arise from noise in the PSF measurement. The impact of measurement noise on the PSF model increases as the number of stars available to characterise the model at each position in the field decreases. The PSF modelling strategy (see Sect. 5.1.1) means that any local random errors from the sparse PSF measurement will propagate as random but correlated errors across the PSF model for the full field of view.

We mimic the impact of random but correlated PSF residual errors by assigning a randomly generated number from a Gaussian distribution with zero mean and unit dispersion to each 5 × 5 arcsec pixel within a 10 ° ×10° field. We first verify that the uncorrelated version of this systematic does not produce any coherent signal, as expected from a random error. We then correlate the random PSF measurement noise over each pointing using a Gaussian filter convolution defined within the boundaries of the pointing. These convolved fields are then renormalized to produce an overall dispersion equal to 10% of the shear dispersion in the mock data, σ_RCN = 0.1σ_γ. We investigate two kernel sizes with a correlation length of roughly the CCD chip level (∼1.6 arcmin) and the pointing scale level (∼43 arcmin). The resulting systematic patterns are shown in the two right panels of Fig. 6, where the systematic ranges are $-0.006<\beta \delta {ϵ}_i^{*}<0.006$ (chip-level correlation) and $-0.003<\beta \delta {ϵ}_i^{*}<0.003$ (pointing-level correlation). For this systematic we chose both components of the contaminating ellipticity to be equal.

5.2. Photometric redshift selection bias

Cosmic shear surveys rely on photometric measurements to estimate the redshifts of galaxies. The photometric redshift (photo-z) of a galaxy can be estimated by comparing the magnitude of the galaxy in several colour-bands to template catalogues of galaxy spectral energy distributions (SEDs) or to spectroscopic training samples (see Salvato et al. 2019, and references therein). The most probable value for the redshift of each galaxy, given the measured photometric colours, z_phot, is then used to divide the galaxy sample into tomographic redshift bins. The true redshifts of these galaxies may not all lie within the boundaries of their appointed photometric redshift bins but provided the true underlying redshift distribution of the galaxies is known, this can be accounted for in the theoretical predictions of the cosmic shear signal (Eq. (6)). The dispersion in true redshift within these tomographic bins will however depend on the precision of each galaxy’s photo-z estimate, which in turn depends on the error on the measured flux of the galaxy in each photometric band. As a galaxy imaged with different noise realisations can therefore appear in different photometric redshift bins, in cases where the flux error is correlated with the shape or orientation of the galaxy, selecting a galaxy sample based on photometric redshifts could therefore lead to an ellipticity-redshift selection bias and hence a systematic error in a cosmic shear analysis. In this section we explore and introduce this new concept of photometric redshift selection bias as a systematic for lensing surveys.

Consider two identical elliptical galaxies of fixed size and flux. The first galaxy is convolved with an elliptical PSF aligned with its major axis. The second is convolved with an elliptical PSF aligned perpendicular to its major axis. In the resulting image our second galaxy will appear to cover a larger area than our first galaxy and with a lower surface brightness and lower significance. It will therefore have larger photometric errors compared to the first galaxy. Kaiser (2000) recognised that this effect implied that any cuts on observed significance would introduce a PSF-dependent selection bias in the ellipticity of the galaxies (see also Bernstein & Jarvis 2002). The introduction of tomographic photo-z selection in a cosmic shear analysis, which implicitly depends on the significance of each detection, can therefore also lead to an ellipticity-dependent selection bias.

In addition to this core effect, flux errors that are correlated with the relative orientation of the galaxy and PSF can also arise simply from the methodology used to measure the photometry in each band. DES-SV use SEXTRACTOR automated aperture magnitudes where the aperture is fixed by the galaxy shape in the detection image (Bonnett et al. 2016; Rykoff et al. 2016). Whilst this method ensures that the physical apertures are matched between the bands, it does not take into account the differing PSFs. Hildebrandt et al. (2012) show that this approach leads to an overall degradation in the photometric redshifts. For example, if the PSF in the i-band is perpendicular to the PSF in the detection r-band, the resulting i-band flux, assuming a fixed-detection aperture, will be underestimated. This approach therefore results in flux errors that are correlated with the relative orientation of the galaxy and PSF in each band. Hildebrandt et al. (2012) demonstrate the importance of homogenising the PSFs between bands before determining the matched-aperture photometry. Both CFHTLenS and KiDS Gaussianise the PSFs before measuring the photometry using the methodology proposed by Kuijken (2008) and Kuijken et al. (2015). These surveys should therefore be fairly immune to this additional error and we note that the DES photometry methodology has been significantly improved since the release of the DES-SV data that we analyse in this paper (Drlica-Wagner et al. 2018).

In this analysis we make the first step to examine photometric-redshift selection bias, by simulating a mock galaxy catalogue where we introduce an anti-correlation between the signal-to-noise ratio of the measured flux and the ellipticity of the galaxies relative to the local PSF ellipticity, $|ϵ-{ϵ}_{\text{x}}^{*}|$, in four bands x = u, g, r, i. We use the following linear relation for the anti-correlation,

$$\begin{array}{c}\hfill \frac{\mathrm{Flux}}{\mathrm{Flux}\mathrm{error}}=a_x|ϵ-{ϵ}_x^{\ast }|+b_x,\end{array}$$(28)

where the value for a_x and b_x are determined by fitting to KiDS-450 multi-band data (see Table D.2). Given that KiDS can only measure the noisy observed ellipticity, we recognise that the majority of the anti-correlation that we find in the KiDS-450 data, derive from taking the mean of the absolute value of the observed ellipticity where the measurement noise, η in Eq. (27), increases with decreasing signal-to-noise. Using Eq. (28) to apply a correlation between the signal-to-noise of a galaxy and its relative ellipticity to the local PSF therefore provides an upper limit for this effect. Future work will use multi-band image simulations to determine values for a_x and b_x where the true ellipticity is known. Our current approach is however sufficient for the purposes of examining the B-mode signature that is introduced by such an effect.

We produce mock ellipticity catalogues by randomly associating ellipticities to galaxies, using a fit to the observed KiDS-450 galaxy ellipticity distribution, such that ⟨ϵ_mock⟩=0. We simulate 15 fields of 100 deg² each with a total galaxy number density of 5.5 arcmin⁻². We choose a simple model of constant PSF per 1 deg² pointing taken randomly from a uniform distribution between [ − 0.1, 0.1] for each component of the PSF ellipticity¹³.

We associate noise-free multi-band fluxes to the mock galaxies using simulations similar to the ones presented in Sect. 3.1 of Hildebrandt et al. (2009) but adapted to KiDS. These simulations were created with the HYPERZ package (Bolzonella et al. 2000) and are based on SED templates from the library of Bruzual & Charlot (1993), i-band number counts from COSMOS (Capak et al. 2007), and redshift distributions from BPZ (Benítez 2000). These magnitude simulations contain half a million galaxies with magnitudes given in each of the four bands, selected to recover the KiDS redshift distributions given in Fig. 7.

Fig. 7. True redshift distribution of the mock galaxies, separated into photometric redshift bins. The zB selection is shown in the legend. The cyan histogram shows the true redshift distribution of the galaxies in the parent noise-free sample. In order to determine the photometric redshifts we introduce flux errors that mimic a KiDS-like survey and depend on the relative ellipticity of the mock galaxy to the mock PSF.

For each galaxy, we assign an error on the flux in each band using Eq. (28). Noise is then added to the mock galaxy flux, sampling from a Gaussian distribution. This approach correlates high values of observed galaxy ellipticity with high flux errors, as expected from the ellipticity measurement noise in the data. In addition, the flux error may also depend on the relative-orientation of the galaxy and the PSF, in each band, as expected from some photometry measurement methods in addition to the Kaiser (2000) effect. As this is the first investigation into photometric redshift selection bias, we have not tried to separate these effects in our analysis. We also note that this method of assigning noise to our mock galaxy sample ignores the additional dependence of the signal-to-noise ratio on galaxy size and magnitude. Future work will need to investigate this in more detail, using multi-colour image simulations.

We use the Bayesian Photometric Redshifts software BPZ to estimate photo-z’s for each of our mock galaxies using a template fitting method (Benítez 2000; Benítez et al. 2004; Coe et al. 2006). The inputs are the noisy flux measurements and their associated errors. The output is the best fitting photometric redshift, z_B, which we use to then bin the galaxies into the four redshift bins that were used in the KiDS-450 cosmic shear analysis, z_i ≤ z_B <  z_i + 1, with z_i = {0.1, 0.3, 0.5, 0.7, 0.9} as well as a broad single bin encompassing the full redshift range of KiDS-450, 0.1 ≤ z_B <  0.9. DES-SV and CFHTLenS use a similar number of tomographic bins, spanning similar ranges in photometric redshifts. Note that the SED templates that we use in BPZ to estimate the redshift of the mock galaxies is independent from the ones used to make the mocks, which shows the robustness of this method to the choice of templates.

Figure 7 shows the true redshift distribution of the mock galaxies for each tomographic redshift bin in z_B. The distributions are broad due to the noise with extended high and low redshift tails which we label as catastrophic outliers in the distribution. The mean and median of each tomographic bin is similar to those in the KiDS-450 data, demonstrating that our method to assign noise to our mock galaxy sample is sufficient for this analysis.

6. Results: mock E/B-modes

In this section we present the two-point statistic signatures of the systematics introduced in Sect. 5.1 using the statistics explained in Sect. 2. We show the effect of the three redshift-independent shear systematics in Sect. 6.2 considering only a single redshift distribution for the galaxies. The effect of the remaining systematic, photometric redshift selection bias, is explored in Sect. 6.3 where redshift binning is applied. With these signatures identified, our goal is to use B-mode measurements as a diagnostic tool to uncover the origin of the systematic signals identified in DES-SV, CFHTLenS and KiDS in Sect. 4. We can also determine the impact of these systematics on the measured E-mode signals. Our approach is complementary to previous studies by Amara & Réfrégier (2008), Kitching et al. (2016), Taylor & Kitching (2018) who propagated cosmic shear systematics through to cosmological parameters in order to set requirements on their (in)significance.

6.1. SLICS shear simulations

The basis of our systematics analysis makes use of the ensemble of mock KiDS-450 catalogues constructed from the SLICS¹⁴ simulations suite described in Harnois-Déraps & van Waerbeke (2015) and Harnois-Déraps et al. (2018). Each SLICS line of sight corresponds to a 10 ° ×10° field that includes galaxy positions, shear and their true and photometric redshifts. On average the redshift distribution and number density of galaxies in these mocks correspond to the KiDS-450 data, which is not so dissimilar from the properties of both DES-SV and CFHTLenS.

In Fig. 8 we present the SLICS cosmic shear measurements, ξ_±, ξ_E/B and COSEBIs, averaged over 10 shape noise-free lines of sight (i.e. ϵ^int = 0). On the left we show ξ_± and ξ_E/B, for θ ∈ [0.5′,300′] in 50 logarithmic bins. The top left panel shows ξ₊ and ξ_E and the lower left panel shows ξ₋ and ξ_B. The right panels belong to the E-mode COSEBIs for a range of angular scales. The measurements from SLICS can be compared to the theoretical prediction (Eqs. (4), (10), (13)), shown as thick solid curves for ξ and pluses for COSEBIs. Here we adopt a flat ΛCDM model given by the input cosmology of the SLICS simulations in Table 1. We use a Bond & Efstathiou (1984) transfer function to estimate the linear matter power spectrum and the Smith et al. (2003) halofit model for the non-linear scales. This combination, although dated, was chosen as the resulting theoretical predictions fit the mocks better than the more modern alternatives (Harnois-Déraps & van Waerbeke 2015). The thin coloured lines in Fig. 8 show the measured values for each line-of-sight (LOS), which show a considerable scatter, especially for larger angular scales. Even with the inclusion of a larger number of LOS, however, we do not expect the theory to match the mean of the mocks perfectly, as the finite box-size of the N-body simulations, where L_box = 505 h⁻¹ Mpc, results in a loss of power on large scales (Harnois-Déraps et al. 2018).

Fig. 8. SLICS 2-point statistics, ξ± and ξE/B (left) and E-mode COSEBIs (right), averaged over 10 noise-free lines-of-sight, which serves as our fiducial “systematics-free” measurement. The mean result can be compared to the theoretical expectation (smooth solid curves). For ξ± and COSEBIs we also show the measurements for each individual line-of-sight with thin solid curves with matching colors between different panels. The upper left panel shows the measured ξ+ (magenta squares) and ξE (blue diamonds), with the lower left panel showing ξ− (green pluses) and ξB (black crosses). The expectation value for ξB is zero, shown with the dashed black line. The COSEBIs E-modes (right panels) are shown for the three angular ranges indicated in each row. The COSEBIs B-modes in SLICS are 4 orders of magnitude smaller than the E-modes and are therefore not shown. The measurements are shown as squares and their expected theory value as plus symbols. Note that COSEBIs modes are discrete and the points are only connected together as a visual aid.

For a B-mode free dataset, ξ₊ = ξ_E. In the upper left panel of Fig. 8, we see that this is not the case for SLICS as at large angular scales ξ_E is smaller than ξ₊. Looking at the lower panel we see that ξ_B is non-zero for the same angular ranges. This leakage from E to B in the ξ_E/B statistic is a result of using the theoretical ξ₋ at large scales for calculating the integrals in Eq. (9), which differs from the ξ₋ of SLICS due to its finite box-size. We find that COSEBIs do not suffer from either of these effects, with the COSEBIs B-mode signal in the mocks found to be ∼4 orders of magnitude smaller than the E-modes (not shown). The reason for the robustness of COSEBIs to the finite box bias comes from the weight functions that convert the shear power spectrum to these statistics. The low ℓ behaviour of the weight functions for COSEBIs have a leading-order term proportional to ℓ⁴ such that the function reaches zero at small ℓ-values in contrast to the ξ₊ kernel, J₀, which has power at small arguments (see Eq. (4)). At high ℓ the COSEBI weights also diminish rapidly in contrast to the 2PCFs which include some degree of power from all scales (see Fig. 4 in Kilbinger et al. 2017, for a comparison between the kernels corresponding to COSEBIs and 2PCFs).

6.2. The B-mode signature of shear measurement systematics

We add the shear measurement systematic effects, developed in Sect. 5.1, in turn to the SLICS simulations. We follow the standard approach of applying an empirical systematics correction to each mock by subtracting the average observed ellipticity from each line of sight before commencing our statistical analysis. In Fig. 9 we compare the resulting 2PCFs with the signal measured in the systematic-free fiducial data (see Fig. 8). The 2PCFs are ξ₊ (squares), ξ₋ (pluses), ξ_E (blue diamonds) and ξ_B (crosses). The left panels show the fractional difference of ξ_± and ξ_E to their fiducial values, calculated from systematic free mocks and shown with a “fid” superscript, for each systematic, while the right panels show the difference between ξ_B and its fiducial value as well as ξ_E − ξ_B (pluses) for each case. Here we use 50 logarithmic θ-bins, to show the angular dependence of each systematic in detail.

Fig. 9. Impact of shear measurement systematics on ξ± and ξE/B for four different types of shear measurement systematics; From the top down: PSF leakage, a repeating additive pattern, and random but correlated noise, correlated on chip and pointing scales (see Fig. 6). For ξ+ (magenta squares), ξ− (green pluses) and ξE (blue diamonds) we present, in the left panels, the fractional difference between the measured signal in the systematic-induced KiDS-like SLICS mocks and the fiducial systematic-free case. As ξB (black crosses) and the E/B difference ξE − ξB (red pluses) tends to zero, we present, in the right panels, the difference between these measurements and the fiducial case, multiplied by the angular distance in arcminutes and scaled by 104. The measured B-modes can be compared to the expected shape-noise error for KiDS-450 (shaded area).

Each row in Fig. 9 shows the impact of the systematic which from top-down cover PSF leakage, a repeating additive pattern, and random but correlated noise (RCN), similar to PSF residuals, correlated over chip-scales and then pointing scales. The grey regions in the right panels show the level of noise expected for KiDS-450 data, which is similar to the noise in the DES-SV and CFHTLenS analyses. As the simulations are free of shape-noise, the error associated with them is negligible compared to the expected errors from either of the three surveys, therefore we have excluded the simulation error from this figure (Fig. 9) and the next one (Fig. 10). We note that all these systematics also produce parity violating signals ξ_× = ⟨ϵ_tϵ_×⟩, which is expected to be zero for a shear only field. We find that their amplitude is about an order of magnitude smaller than the B-modes, however, and are therefore harder to detect in the data. As a result, we limit our systematics study to the effect of systematics on E/B-modes.

Fig. 10. Impact of shear measurement systematics on E-mode (left) and B-mode (right) COSEBIs for four different types of shear measurement systematics; PSF leakage (blue squares), a repeating additive pattern (black stars), and random but correlated noise on chip (green triangles) and pointing (magenta diamonds) scales (see Fig. 6). The analysis is conducted for three different angular ranges spanning [0.5′,40′] (upper panels), [0.5′,100′] (middle panel), and [40′,100′] (lower panels). We present the difference between the measured signal in the systematic-induced KiDS-like SLICS mocks and the fiducial systematic-free case scaled by 1010. The measured B-modes and the resulting change to the E-mode can be compared to the expected shape-noise error for KiDS-450 (shaded area).

One interesting result from Fig. 9 comes from the non-zero signal in the $({\xi }_E-{\xi }_E^{\mathrm{fid}})-({\xi }_B-{\xi }_B^{\mathrm{fid}})$ curves. If a systematic adds equal power to both observed E and B-modes, P_sys, from Eq. (7) and Eq. (8) we find that the observed ξ_B is equal to the excess signal in the observed ξ_E due to this systematic,

$$\begin{array}{c}\hfill {\xi }_E^{\mathrm{obs}}-{\xi }_E^{\mathrm{true}}={\xi }_B^{\mathrm{obs}},\end{array}$$(29)

where ${\xi }_{E/B}^{\mathrm{obs}}$ are the observed E/B-mode 2PCFs and ${\xi }_E^{\mathrm{true}}$ is the E-mode signal produced from a shear only field. Furthermore as ξ₋ is proportional to P_E − P_B, an equal P_sys contribution to the E and B-modes will cancel such that ${\xi }_{-}^{\text{obs}}-{\xi }_{-}^{\text{true}}=0$. In this case there is a clear route to correct the measured E-mode by the measured B-mode or to select “clean” angular scales for the E-mode analysis which are B-mode free (see for example Hildebrandt et al. 2017). For the sample of systematics that we have simulated, however, we see that this common assumption of the equal contribution of systematic power to the E and B modes is far from reality, especially for larger angular scales¹⁵. This implies that if a ξ_B signal is detected at any angular scale, its origin should be identified and mitigated at the catalogue or image level. Without understanding the origin it is unclear how that systematic will contaminate the ξ_± signal.

In Fig. 10 we present the COSEBIs analysis of the mocks. We show the relative effect of systematics on the E-mode (left) and B-mode (right) COSEBIs, as the difference between their fiducial values and those estimated from the systematic induced mocks including PSF leakage, αϵ^* (blue squares), repeating additive pattern, c(x, y) (black stars), and random but correlated noise (RCN), βδϵ^*, correlated on chip scales (green triangles) and pointing scales (magenta diamonds). All values are shown for the mean of the 10 SLICS lines-of-sight. The grey regions show the one sigma errors corresponding to a KiDS-450-like survey. Random but correlated noise at the chip level shows small deviations from the fiducial values in agreement with Fig. 9. Within a single pointing, this systematic has a similar form to the additive pattern, where we find similar low-n behaviour between these two systematics with similar peaks at n = 2 and a dip at n = 4, 5, albeit at different amplitudes. PSF leakage and the random but correlated noise at the pointing level are more significant, exhibiting a similar signal from the lower COSEBIs modes. The repeating additive pattern has the most chaotic effect on COSEBIs, in comparison to the other systematics that we have simulated. The erratic high frequency changes that can be seen in the 2PCFs in Fig. 9 are reflected in the significant power seen in the higher COSEBIs modes. As these systematics produce varying correlations for different angular scales, COSEBIs modes are affected by them in differing amplitudes, which also depend on the angular range they probe. Comparing Fig. 9 and Fig. 10 gives insight into the sensitivity of COSEBIs modes to correlation at various angular scales.

For all four systematics, we find that their effects on the COSEBIs are more prominent when the full angular range is used (middle panels). Comparing the right and left panels we see that all systematics affect both E and B-modes, but not equally. In general a significant B-mode signal translates to a significant contamination to the E-modes on the same scales. The repeating additive pattern forms a clear exception to this rule though. This draws us to the same conclusion as the ξ_B analysis, the origin of any COSEBIs B-mode signal should be traced back to its source and corrected for at the pixel-data product level where phase information is still available, since it is unclear how these systematics will impact the E-mode at another angular scale.

The characteristic patterns that we have identified should be used in any future approach to diagnose and correct shear-measurement systematics. As an example, if COSEBI B-modes are found to be oscillatory and extend to high-n, the survey should investigate additive biases that repeat on a fixed angular scale across the survey, for example detector-level effects. If the B-modes are localised at low-n with little high-n power, the survey should investigate the PSF modelling. With COSEBIs alone we cannot distinguish between PSF leakage, αϵ^* or correlated noise in the PSF model, βδϵ^*, but these two effects can be separated by measuring the correlation between galaxy shape and PSF ellipticity, which will be significant if α is non-zero (Bacon et al. 2003). If our ellipticity model in Eq. (27) is reasonable in its approach to add systematic terms linearly, we would expect the B-modes from each individual effect to also add linearly. When we see both significant power at n <  7 and high-n oscillatory power, as we do for DES-SV for example, we can conclude that there is a likely superposition of systematics from both the PSF modelling errors and repeating additive biases.

The simulation approach that we use here should be specialised to the survey in question in future work. This would allow for a more precise exploration of how survey-specific issues flow through to cosmological biases. In this analysis we have presented results that use KiDS to motivate the angular dependence of the systematics that we have simulated. In addition, we have also tested a variety of alternative schemes in the development of this work such that we are confident that the global behaviour of the B-mode signatures, presented in Fig. 10, are broadly representative of how these systematics would feature in any weak lensing survey. These tests included modelling different possible patterns and fixed amplitudes for each of the systematics, with the pointing size fixed to 1 deg² in all cases in accordance to both KiDS and CFHTLenS. DES pointings are however hexagons of width 2.2°.

6.3. The B-mode signature of photometric redshift selection bias

We simulate photometric redshift selection bias following the approach developed in Sect. 5.2. As this effect is subtle, we choose to analyse the correlations in a random intrinsic ellipticity field alone, in contrast to the analysis in Sect. 6.2. where we also include the correlated SLICS cosmic shear field. Starting with the traditional statistics in Fig. 11, we compare two cases. The top panels show the 2PCFs using the full galaxy sample which are consistent with zero by construction. The lower panels show the same 2PCF analysis, including a photometric redshift selection with 0.1 ≤ z_phot ≤ 0.9. Comparing these two results we immediately see that the photometric redshift selection has produced a significant signal in all but the ξ₋ statistic. As such we find similar levels of E/B-modes in the lower right panel where the difference between ξ_E − ξ_B ≈ 0. Comparing the bottom panel of Fig. 11 with the top panel of Fig. 9 we see that the location of the peak in ξ₊ is similar. This is because the PSF-leakage systematic (Fig. 9) is related to photometric redshift selection bias, although, here we have used a simplified PSF model without the small scale variations that affect the PSF-leakage systematic (Fig. 6).

Fig. 11. Photometric redshift selection bias: the upper panels show the correlations measured from a random intrinsic ellipticity field using the full galaxy sample. Lower panels: correlations measured from the same random intrinsic ellipticity field after a 0.1 ≤ zphot ≤ 0.9 photometric redshift selection has been applied. The 2PCFs shown on the left are ξ+ (squares) and ξ− (pluses). On the right panel, ξE (diamonds), ξB (crosses) and the difference between the two (pluses) are shown. The one-sigma error bars correspond to the level of ellipticity noise in the mock data. The signals are shown multiplied by θ in arcminutes and scaled by 104.

To see the full effect of the photometric redshift selection bias on the measured correlation function, we use COSEBIs to analyse four photometric redshift bins, corresponding to those chosen in the KiDS-450 analysis. Fig. 12 shows the results revealing significant signal in the different tomographic slices, with the strongest effect in the highest redshift bins. To quantify the significance of this effect we use a theoretical covariance to estimate a χ² value for COSEBIs relative to the null hypothesis of zero signal. We then calculate the p-value corresponding to that χ² finding vanishingly small values of p ∼ 10⁻²⁸ and p ∼ 10⁻¹⁵ for the E and B modes respectively, clearly confirming the significance of the bias.

Fig. 12. Impact of photometric redshift selection bias for a COSEBIs 4-bin tomographic analysis of a random intrinsic ellipticity field. The upper triangle shows the E-modes and the lower triangle shows the B-modes, where the error bars in both cases correspond to the level of ellipticity noise in the mocks. Each panel shows COSEBIs for a tomographic redshift bin pair, z − ij, corresponding to the correlation between photometric redshift bins i and j. As our photometric redshift mocks are devoid of any cosmological correlations, in the absence of any selection bias, we would expect both the E and B modes to be consistent with zero.

As with the analysis of the shear measurement systematics in Sect. 6.2 we see similarities (generally low-n) and differences (generally high-n) in the measured E and B-modes, again leading us to the conclusion that B-modes can be used as a diagnostic but cannot blindly be used to correct the E-modes.

6.4. Cosmological parameter inference

Although we can see the signature of each systematic in Figs. 9 to 12, it is not immediately clear how they would affect cosmological parameter inference. One could carry out a likelihood analysis to find any biases introduced by these systematics (see for example Amara & Réfrégier 2008), but our preference is to use compressed COSEBIs (CCOSEBIs, see Sect. 2.3) as a faster alternative approach. CCOSEBIs are formed of linear combinations of COSEBIs that are sensitive to cosmological parameters. If the systematics that we identify in the COSEBIs analyses are null in both the E and B mode CCOSEBIs case, then we can conclude that the systematics are unlikely to be detrimental to the cosmological inference.

In this analysis we focus on the CCOSEBIs that are sensitive to Σ₈ = σ₈(Ω_m/0.3)^α, as this is the combination of parameters that cosmic shear data are mostly sensitive to. We find that for a KiDS-450 redshift distribution, α = 0.65 best describes the COSEBIs degeneracy direction (see Appendix B). Because we are interested in Σ₈ we only consider the 5 first and second-order CCOSEBIs for σ₈ and Ω_m; $E_{{\sigma }_8}^c$, $E_{{\mathrm{\Omega }}_{\text{m}}}^c$, $E_{{\sigma }_8{\sigma }_8}^c$, $E_{{\sigma }_8{\mathrm{\Omega }}_{\text{m}}}^c$ and $E_{{\Omega }_{\text{m}}{\Omega }_{\text{m}}}^c$, where E^c is defined in Eq. (15). Although it is possible to construct a CCOSEBIs mode that is sensitive to Σ₈ directly, we choose to look at these 5 modes instead in order to also provide an internal consistency check. For each of the five modes we calculate a single compressed value, $E_{\mu (\nu )}^c$, that can be compared to its expectation value, given a set of cosmological parameters, noise covariance and source redshift distribution.

Figure 13 shows the measured E-mode (left panel) and B-mode (right panel) CCOSEBIs for the full angular range of [0.5′,100′]. The symbols correspond to the range of shear measurement systematics¹⁶ simulated using the SLICS cosmic shear simulations in Sect. 5.1. The lines connecting the E-mode points indicate which of the 5 CCOSEBIs modes are shown and also show their theoretical value. Each measured E-mode can be compared to the value of Σ₈ = σ₈(Ω_m/0.3)^0.65 that would be inferred from the measurement. In the absence of systematic errors, we would expect to find the inferred parameters to be consistent with each other and the input SLICS cosmology with Σ₈ = 0.808. We would also expect to find the B-mode signal consistent with zero, but looking at the fiducial “no-systematics” mocks (circles), we do recover a very small residual B-mode and a slightly high best-fit Σ₈ = 0.815. This result is expected, however, given the imperfect match between the two-point statistics measured from SLICS and the theoretical expectation shown in Fig. 8. Here no errors are associated to the E-mode CCOSEBIs, because the mock data used to produce them are free of shape-noise.

Fig. 13. Left panel: inferred values of Σ8 = σ8(Ωm/0.3)0.65 from an E-mode CCOSEBIs analysis of four mock cosmic shear surveys that suffer from PSF leakage (blue squares), a repeating additive pattern (black stars), and random but correlated noise on chip (green triangles) and pointing (magenta diamond) scales. The curves show the theoretical values of the 5 E-mode CCOSEBIs when varying Σ8, calculated using the KiDS-450 noise only covariance matrix and redshift distribution; Eσ8σ8 (dashed), EΩmΩm (dotted), Eσ8 (middle solid), EΩm (dot-dashed) and Eσ8Ωm (lower solid). The inferred cosmology for each mock systematic survey can be compared between the 5 different modes. The higher the recovered E-mode is relative to the fiducial “no-systematics case” (circles), the stronger the bias is on the inferred value of Σ8 and the more discrepant the inferred cosmology is between the different CCOSEBIs modes. The bias in Σ8 can be compared to the grey region which shows the one-sigma error for Σ8 from the KiDS-450 cosmic shear analysis, centred on the fiducial case. Right panel:B-mode CCOSEBIs from the four mock cosmic shear surveys. The measured B-mode signal, which does not depend on Σ8, can be compared to the shape noise on a KiDS-450-like survey (shown in grey).

We find that the introduction of the random, but correlated noise (RCN) increases the recovered Σ₈ value, but within the statistical tolerance of KiDS (shown as a grey bar) in the case of the chip-scale correlation. The PSF leakage and repeating additive pattern result in the largest bias in cosmological parameters with a ∼5% deviation from the true input cosmology. Applying this level of bias to either σ₈ or Ω_m can produce excess correlations of only up to 13% which is significantly less than the up to 40% biases seen in the two-point correlation function analysis in Fig. 9, from which we can conclude that the impact of these systematics on the data can be only weakly correlated with the impact of varying cosmological parameters.

We find that the stronger the bias in the recovered cosmology, the larger the inconsistency between the 5 CCOSEBIs modes, providing another important diagnostic tool. We also note that all the shear measurement systematics tested in this analysis serve to increase the inferred value of Σ₈. If these types of systematics were present in the weak lensing data, correcting for them would decrease the recovered Σ₈, exacerbating the current hints of cosmological parameter tension between weak lensing surveys and Planck (see for example Troxel et al. 2018b).

Comparing the power in the measured E-modes (left panel) and B-modes (right panel) reveals a close connection, where a larger bias in the E-modes corresponds to larger B-modes. We note that although the magnitude of the E/B-modes are connected, they can take opposite signs. For example, in Fig. 10 we see that for the large angular scale analysis with the repeating additive pattern systematic, the sign of the first four E-modes differs from the sign of the first four B-modes. Although we do not show the large-scale CCOSEBIs result, we can confirm that, this difference in sign is also reflected there, as the CCOSEBIs are sensitive to the first few modes that contain a large proportion of the cosmological information. This example demonstrates that although the large angular scales have the lowest signal to noise, they can and should be used as an investigative tool for hunting systematics that could also impact small angular scales.

7. Discussion

In this section we discuss how we can use the measured COSEBIs B-mode signatures from our systematics mocks in Sect. 6 to diagnose the origin of the B-modes recovered in the CFHTLenS, DES and KiDS surveys in Sect. 4. Firstly, we focus on the non-tomographic COSEBIs B-mode measurements in Fig. 3. One feature that stands out for all three surveys is the high n-mode oscillatory pattern in the full angular range, shown in the middle right panel. This oscillatory pattern is the signature of a repeating additive systematic, shown in Fig. 10, which we find to be the systematic that was most detrimental to cosmological parameter estimation in Fig. 13.

We find the level of B-modes for KiDS-450 and CFHTLenS for higher n-modes to be small and hence the repeating additive signature is not highly significant in these cases. The similarity between the B-modes in the data and this systematic signature does however warrant further exploration, particularly as we also see similarities in the E-modes for KiDS and the repeating pattern E-mode signature. Here the unexpected E-mode “secondary peak” seen at n ∼ 6 in the small angular scales of KiDS-450 E-modes (upper left panel of Fig. 3) is replicated at n ∼ 6 in the E-mode analysis of the same angular scales of the repeating additive bias mocks (upper left panel of Fig. 10). If a repeating additive systematic persists it would likely become significant in future releases of KiDS. It could also be responsible for the power seen in the low-n modes that lead to the significant detection of the KiDS CCOSEBIs B-modes.

For DES-SV we find a significant detection of B-modes, noting that in addition to the high-n oscillatory pattern, DES-SV presents significant additional signal for modes around n = 8 and n = 4. For an instrument-based repeating additive pattern, the resulting B-mode signature will depend on the dithering strategy and camera field-of-view. Both KiDS-450 and CFHTLenS have a field-of-view of ∼1 deg² with small dithers. DES-SV, however, has a hexagonal field-of-view, 2.2° across, and uses half-field dithers. This means the frequency of any repeating additive pattern will differ for DES-SV in comparison to the KiDS-like imaging strategy that we have simulated in our mocks. Looking at only the first few modes for DES-SV data, however, we find that the signal resembles the signature of both PSF leakage, and random but correlated noise on the pointing level. This result is consistent with the findings of Zuntz et al. (2018) who report and correct for a small but significant PSF residual in their analysis of the first year of DES observations. We therefore conclude that the B-mode signature recovered for DES-SV is likely a superposition of different shear measurement systematics.

By comparing the p-values in Table 2 we can see that for DES-SV the significance of the COSEBIs B-modes substantially increases when the data are separated into tomographic bins. This could be understood by considering the photometric redshift selection bias explained in Sect. 5.2. This systematic correlates the PSF ellipticity with the redshift estimation for a galaxy, and can produce significant B-modes when the data are binned into smaller photometric redshift bins. It is likely that all surveys will suffer from this systematic to some degree, but the level will depend on how the multi-band photometry is measured in each survey and how the PSF ellipticity varies in each optical band. We cannot directly compare our mock analysis with the B-modes in the DES-SV tomographic analysis, but our first-look at this effect certainly motivates further exploration with more detailed simulations that fully mimic the photometric redshift measurement in each survey.

Interestingly, comparing the DES-SV tomographic and non-tomographic p-values in Table 3, we find that for the analyses that include small-scale information, the significance of the B-modes, measured using the cosmological-parameter-sensitive CCOSEBIs, decreases when the data are separated into tomographic bins. This promising result means that if the systematic that was introduced when the DES-SV tomographic selection is applied adds equal power to the E and B modes, that systematic would not introduce modifications to the E-mode signal that would bias the inferred cosmological parameters. Unfortunately, however, the photometric redshift selection bias systematic was found to exhibit different E and B mode signals in Fig. 12. Passing the CCOSEBIs B-mode null-test therefore cannot validate the CCOSEBIs E-mode measurement. In addition, this CCOSEBIs B-mode result does not hold for the large angular scales, [40′,100′], where again we see a substantial increase in the measured B-mode when the data are separated into tomographic bins.

Our findings for DES-SV contradict Becker et al. (2016) who conclude that the B-modes in DES-SV are insignificant using two Fourier space methods. We argue that as power spectra pick up features of the data at different scales compared to ξ_±, they are not suitable statistics for verifying the absence of B-modes in a ξ_± cosmic shear analysis. In addition, power spectra measurements are binned in a range of Fourier modes, such that any high-frequency variations in Fourier space will average out. COSEBIs are sensitive to these variations and can therefore be used to diagnose the origin of the B-modes in the data.

At first sight our findings for KiDS also contradict Hildebrandt et al. (2017) who report a low-level but significant detection of ξ_B. This is in contrast to our tomographic ξ_B analysis which concludes that ξ_B is consistent with zero (see Fig. D.5). We find that the ξ_B statistic is sensitive to the choice of the maximum θ-scale measured from the data and the maximum θ-scale used for completing the integral to infinity using a theoretical prediction (in this analysis we use 1000′ instead of 3000′ used in Hildebrandt et al. 2017). We also find that ξ_B is sensitive to the method used to bin ξ_± as explained in Appendix A (see Eq. (A.10)). This sensitivity to data analysis choices provides another reason to archive the traditional ξ_E/B approach. In this paper we promote COSEBIs as the optimal statistic for both E and B mode measurements as it can be estimated accurately and free of any biases connected to binning and extrapolating the data. Analysing all 20 COSEBIs modes, we find no significant evidence for B-modes in KiDS. In our compressed CCOSEBIs analysis, however, we arrive at the same conclusion of both Hildebrandt et al. (2017) and van Uitert et al. (2018), that low-level but significant B-modes are present in KiDS-450. In our [0.5′,100′] CCOSEBIs analysis, we find a ∼2.7σ detection of a B-mode signal that is less than 10% of the amplitude of the E-mode. This difference between the significance of the COSEBIs and CCOSEBIs B-mode analysis might seem confusing or even contradictory. We therefore refer the reader to Appendix C where we explore how choices over the number of modes used in a null-test can dilute or optimise the detection of systematics.

8. Conclusions

Two-point shear correlation functions (2PCFs) have been the primary observables in cosmic shear analysis to date, but they are not immune to systematics. These statistics mix E and B-modes in the data, giving rise to a mixed lensing and non-lensing signal in the presence of systematic errors. In order to test for systematics most surveys turn to alternative statistics to separate E/B-modes, using ξ_E, B or power spectrum measurements. We argue that these alternative statistics are biased as they depend on infinite integrals over 2PCFs and are sensitive to binning choices. In addition, treating the E/B-mode decomposition with a statistic that has a different scale-dependence to the statistic used in the cosmological parameter inference, causes a disparity in the analysis. For future cosmic shear analyses, we therefore advocate the use of COSEBIs for both parameter inference and systematic analyses (see Sect. 2.2). COSEBIs cleanly and completely separate E/B-modes over a finite angular range, without loss of information. They have discrete modes and therefore are insensitive to binning choices. The first few modes of COSEBIs contain almost all of the cosmological information and as such a COSEBIs analysis is also an efficient approach to data compression. For a B-mode analysis, however, a larger number of modes need to be considered, as systematics can affect the E and B-mode at different scales.

In this paper we analysed the E and B-mode signals in three public cosmic shear surveys, CFHTLenS, DES-SV and KiDS-450. We compared the ξ_E, B statistic with COSEBIs and CCOSEBIs, using p-values to quantify the level of B-modes in the data. To determine COSEBIs filter functions we need to first define an angular range of interest. For this study we chose three sets of angular separation ranges: small separations, [0.5′,40′], large separations, [40′,100′], and the overall separation range, [0.5′,100′]. We measure COSEBIs up to mode n = 20. We considered two cases for each survey; one using the same redshift bins as used in each survey’s primary cosmic shear analysis, and another combining those bins into a single redshift bin. We see that for DES-SV data the tomographic cases show significant B-modes at a level between 4σ and 5.5σ. For the non-tomographic DES-SV analyses, B-modes are detected at the level of 2.8σ. For KiDS-450 and CFHTLenS, we find no significant detection of B-modes for the majority of our analyses. There is however some exceptions in each case. The CCOSEBIs analysis of the small separations (non-tomographic case only) and the analysis of the full angular range show B-modes at up to 2.7σ for KiDS-450. The tomographic COSEBIs analysis over the small angular range [0.5′,40′] detects a B-mode signal at 2.8σ for CFHTLenS.

In order to diagnose the origin of the B-modes detected in each survey, we modelled several non-astrophysical systematic effects relevant to current data in order to determine their E/B mode signature and assess their impact on cosmological parameter inference. We modelled four shear measurement systematics. PSF-leakage, was modelled using the mosaic PSF pattern from KiDS-450 assuming a 10% leakage with α = 0.1. An instrument-based additive bias term resulting in a repeating pattern from pointing to pointing. Here we used the low-level CCD bias of OmegaCam (Hoekstra et al., in prep.), multiplied by a factor of five to amplify and model this effect. To model biases arising from random PSF modelling errors, we correlated low levels of random noise using two kernel sizes, corresponding to roughly KiDS CCD and pointing scales. In addition to these shear measurement systematics, we have introduced a new effect by modelling the impact of photometric redshift selection bias that arises from the correlation between the relative orientation of PSF ellipticity and galaxy ellipticity, and the measured signal-to-noise of the galaxy.

All of the systematics simulated were detected in our B-mode analysis. The PSF-leakage and random but correlated noise systematics introduced low n-mode COSEBIs signal. This was in contrast to the repeating additive bias which introduced high frequency variations in the shear field which are picked up as oscillatory behaviour in the high n-mode COSEBIs measurements. Photometric redshift selection bias also resulted in high n-mode power in the high photometric redshift bins. Comparing the B-mode signatures recovered by our mocks to the B-modes measured in each survey we conclude that DES-SV is likely to suffer from a combination of all the systematics that we have simulated. The significant increase in DES-SV B-modes when the tomographic redshift selection is applied is particularly striking, motivating future work to enhance the realism of the first-look photometric redshift simulations that we have analysed in this paper. KiDS-450 and CFHTLenS show oscillatory behaviour in the high n-mode indicating a repeating additive bias in the data, although this result is not significant.

The simulated systematics produce E-modes that would bias cosmological parameter inference. For the analysed shear measurement systematics we found that Σ₈ = σ₈(Ω_m/0.3)^0.65 is biased high in all cases. As a result, we conclude that these types of systematics, if present, cannot explain the mild tension between some current cosmic shear and Planck results and could in principle exacerbate the tension as they bias Σ₈ to even higher values (see Fig. 13). It is interesting to note that the DES-SV cosmological parameter constraints on S₈ = σ₈(Ω_m/0.3)^0.5 are higher than those from KiDS-450, CFHTLenS and the first year DES results, which include a number of improvements over the DES-SV analysis. Given the significant DES-SV B-modes detected in our analysis, the direction of the difference in S₈ between the surveys is expected. The published cosmological parameter constraints from all three surveys are, however, in good agreement.

For the analysis of KiDS-450, we find an interesting case where the survey formally passes the COSEBIs B-mode analysis, but a flag is raised with a 2.7σB-mode detection in the compressed CCOSEBIs analysis. Here the COSEBIs B-modes that are insignificant overall are weighted in such a way that the resulting CCOSEBIs signal becomes significant. We therefore recommend measuring both CCOSEBIs and COSEBIs B-modes in future analyses, ensuring that both are consistent with zero. The CCOSEBIs B-modes will robustly identify systematics that will lead to a bias in the cosmological parameter inference, if the systematic impacts the E/B-modes in the same way. In contrast the COSEBIs B-modes detect systematics that can affect E and B-modes differently. For this reason we have shown that it is not sufficient to correct data by simply subtracting the B-modes from the E-modes. As we have seen from our repeating additive bias systematic, where the E and B mode behaviour is very different, it will be crucial to look at both a COSEBIs and CCOSEBIs analysis and return to correct the input catalogues if a significant B-mode is detected in either case.

Acknowledgments

We thank Joe Zuntz for his help with DES-SV data and COSMOSIS. We also thank Benjamin Joachimi, Alex Hall, Tilman Troester, Michael Troxel and Angus Wright for useful discussions. CH, MA, AA and JHD acknowledge support from the European Research Council under grant number 647112. This work was carried out in part at the Aspen Center for Physics, which is supported by the National Science Foundation grant PHY-1607611, where CH, HH and KK were also supported by a grant from the Simons Foundation. PS is supported by the Deutsche Forschungsgemeinschaft in the framework of the TR33 “The Dark Universe”. KK acknowledges support by the Alexander von Humboldt Foundation. HH is supported by Emmy Noether (Hi 1495/2-1) and Heisenberg grants (Hi 1495/5-1) of the Deutsche Forschungsgemeinschaft as well as an ERC Consolidator Grant (No. 770935). LM acknowledges support from STFC grant ST/N000919/1. AC acknowledges support from NASA grant 15-WFIRST15-0008. Computations for the N-body simulations were performed in part on the Orcinus supercomputer at the WestGrid HPC consortium (www.westgrid.ca), in part on the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund – Research Excellence; and the University of Toronto. The KiDS-450 results in this paper are based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017 and 177.A-3018, and on data products produced by Target/OmegaCEN, INAF-OACN, INAF-OAPD and the KiDS production team, on behalf of the KiDS consortium. The CFHTLenS results in this paper are based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/IRFU, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. CFHTLenS data processing was made possible thanks to significant computing support from the NSERC Research Tools and Instruments grant program. This DES-SV results in this paper are based on public archival data from the Dark Energy Survey (DES). Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the OzDES Membership Consortium, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. Author contributions: All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (MA, CH) followed by two alphabetical groups. The first alphabetical group includes those who are key contributors to both the scientific analysis and the data products. The second group covers those who have either made a significant contribution to the data products, or to the scientific analysis.

References

Appendix A: Binning bias with ξ_±

When using the traditional ξ_± statistic an issue arises from the choices that can be made when binning the data. ξ_±(θ) is usually binned in broad θ-bins as a form of data compression. As we expect the number of galaxy pairs to roughly scale with θ, see Eq. (19), the sampling of ξ_±(θ) within the bin is non-uniform. If we bin these functions into broad angular bins, their value will therefore be biased towards larger θ scales in each bin. This is not an issue provided the theory is treated in the same way, but this is not the standard approach that is taken, as it is computationally more expensive. Troxel et al. (2018b) compare the differences in KiDS-450 cosmological parameter inference if one takes the logarithmic mid-point of the bin or the weighted mean value of θ in each bin and evaluate the theoretical 2PCFs at each θ value. They argue that the latter approach is correct, supported by Krause et al. (2017) who conclude that this approach is sufficiently accurate for the first year DES analysis. Here we provide more detail on the question of binning bias, quantifying how inexact each treatment of the theory is, where we find up to 10% biases in both approaches. A full integration of the theory within the bin is the correct approach to this problem. If future surveys wish to use an approximation, however, we demonstrate that using the linear mid-point of the bin provides the closest match to the binned data, with less than 2.5% bias.

A.1. Binning theory

Consider making measurements of a function f(x) from noisy data, with samples drawn non-uniformly in x. We denote the sampled data points by f_data(x) and the distribution of measured x by 𝒟(x). Given that the sampled data points are noisy, we want to combine them to find an estimate for the function for a given binning in x. One way to bin the data is to write it as

$$\begin{array}{c}\hfill {\widehat{f}}_{\mathrm{binned}}(x_{+\mathrm{b}})=\frac{1}{N_{\mathrm{bin}}}{\displaystyle \sum _x{f}_{\mathrm{data}}(x)\mathrm{\Delta }(x-x_{\mathrm{b}}),}\end{array}$$(A.1)

where N_bin is the number of data points in the given bin and Δ(x − x_b) is the binning function defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }(x-x_{\mathrm{b}})=\{\begin{array}{cc}1\hfill & \mathrm{if}x\mathrm{is}\mathrm{in}\mathrm{bin}{x}_{\mathrm{b}},\hfill \\ 0\hfill & \mathrm{if}x\mathrm{is}\mathrm{not}\mathrm{in}\mathrm{bin}{x}_{\mathrm{b}}.\hfill \end{array}\end{array}$$(A.2)

This estimate for the binned function corresponds to a weighted binning, with more weight given to the values where there are more sampled points. The expectation value of ${\widehat{f}}_{\mathrm{binned}}(x_{\mathrm{b}})$ is in general not equal to f(x_b),

$$\begin{array}{c}\hfill ⟨{\widehat{f}}_{\mathrm{binned}}(x_{\mathrm{b}})⟩=\frac{{\int }_{x_{\min }(x_{\mathrm{b}})}^{x_{\max }(x_{\mathrm{b}})}\mathrm{d}xf(x)\mathrm{D}(x)}{{\int }_{x_{\min }(x_{\mathrm{b}})}^{x_{\max }(x_{\mathrm{b}})}\mathrm{d}x\mathcal{D}(x)}\ne f(x_{\mathrm{b}}),\end{array}$$(A.3)

where x_min(x_b) and x_max(x_b) are the edges of bin x_b. Even if we define x_b as the weighted mean of the x-values in the bin, as advocated by Troxel et al. (2018b),

$$\begin{array}{c}\hfill x_{\mathrm{b}}=\frac{{\int }_{x_{\min }(x_{\mathrm{b}})}^{x_{\max }(x_{\mathrm{b}})}\mathrm{d}xx\mathcal{D}(x)}{{\int }_{x_{\min }(x_{\mathrm{b}})}^{x_{\max }(x_{\mathrm{b}})}\mathrm{d}x\mathcal{D}(x)},\end{array}$$(A.4)

we would only recover the true value of the binned data, if f(x) is either a constant or has a linear relation to x.

Let’s now take very fine x-bins, such that the variation in the sampling of f(x) is negligible, i.e. 𝒟(x) ≈ constant within each bin. The expectation value of the finely binned function $⟨{\widehat{f}}_{\mathrm{f}}(x_{\mathrm{f}})⟩$, where the subscript f represents “fine” and x_f is the mid-point of the fine bin, is given by

$$\begin{array}{cc}\hfill ⟨{\widehat{f}}_{\mathrm{f}}(x_{\mathrm{f}})⟩& =\frac{{\int }_{x_{\min }(x_{\mathrm{f}})}^{x_{\max }(x_{\mathrm{f}})}\mathrm{d}xf(x)\mathcal{D}(x)}{{\int }_{x_{\min }(x_{\mathrm{f}})}^{x_{\max }(x_{\mathrm{f}})}\mathrm{d}x\mathcal{D}(x)}\hfill \\ \hfill & \cong \frac{1}{2\delta x}{\displaystyle {\int }_{x_{\mathrm{f}}-\delta x}^{x_{\mathrm{f}}+\delta x}\mathrm{d}xf(x)\cong f(x_{\mathrm{f}}),}\hfill \end{array}$$(A.5)

where 2 δx is the width of the fine bin. If we first measure finely binned ${\widehat{f}}_{\mathrm{f}}(x_{\mathrm{f}})$ from the data, then we have the flexibility to re-bin the measurements as desired,

$$\begin{array}{c}\hfill {\widehat{f}}_{\mathrm{w}}(x_{\mathrm{b}})=\frac{{\sum }_x{}_{\mathrm{f}}w(x_{\mathrm{f}}){\widehat{f}}_{\mathrm{f}}(x_{\mathrm{f}})\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})}{{{\sum }_x}_{\mathrm{f}}w(x_{\mathrm{f}})\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})},\end{array}$$(A.6)

where w(x_f) is a weight function assigned to each fine bin. If we chose to set w(x_f) to 𝒟(x_f) we would recover the weighted binning defined in Eq. (A.1).

If we choose w(x_f) such that it does not vary between different realisations of the data, the expectation value of ${\widehat{f}}_{\mathrm{w}}(x_{\mathrm{b}})$ is given by

$$\begin{array}{c}\hfill ⟨{\widehat{f}}_{\mathrm{w}}(x_{\mathrm{b}})⟩=\frac{{{\sum }_x}_{\mathrm{f}}w(x_{\mathrm{f}})f(x_{\mathrm{f}})\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})}{{{\sum }_x}_{\mathrm{f}}w(x_{\mathrm{f}})\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})},\end{array}$$(A.7)

and the covariance of the binned data for bins x_b and y_b is given by

$$\begin{array}{cc}\hfill {\mathbb{C}}_{\mathrm{w}}(x_{\mathrm{b}},y_{\mathrm{b}})& =⟨{\widehat{f}}_{\mathrm{w}}(x_{\mathrm{b}}){\widehat{f}}_{\mathrm{w}}(y_{\mathrm{b}})⟩-⟨{\widehat{f}}_{\mathrm{w}}(x_{\mathrm{b}})⟩⟨{\widehat{f}}_{\mathrm{w}}(y_{\mathrm{b}})⟩\hfill \\ \hfill & =\frac{{{\sum }_x}_{\mathrm{f}}{{\sum }_y}_{\mathrm{f}}w(x_{\mathrm{f}})w(y_{\mathrm{f}}){\mathbb{C}}_{\mathrm{f}}(x_{\mathrm{f}},y_{\mathrm{f}})\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})\mathrm{\Delta }(y_{\mathrm{f}}-y_{\mathrm{b}})}{{{\sum }_x}_{\mathrm{f}}{{\sum }_y}_{\mathrm{f}}w(x_{\mathrm{f}})w(y_{\mathrm{f}})\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})\mathrm{\Delta }(y_{\mathrm{f}}-y_{\mathrm{b}})},\hfill \end{array}$$(A.8)

where ℂ_f(x_f, y_f) is the covariance of the finely binned measurements, f_f(x_f) and f_f(y_f). If we assume no cross-correlation between the bins, which is the case for a shape-noise only covariance, then Eq. (A.8) simplifies and the variance of f_w(x_b) can be written as,

$$\begin{array}{cc}& {\sigma }_{\mathrm{w}}^2(x_{\mathrm{b}})=⟨{\widehat{f}}_{\mathrm{w}}^2(x_{\mathrm{b}})⟩-{⟨{\widehat{f}}_{\mathrm{w}}(x_{\mathrm{b}})⟩}^2\hfill \\ \hfill & =\frac{{{\sum }_x}_{\mathrm{f}}{w}^2(x_{\mathrm{f}}){\sigma }_{\mathrm{f}}^2(x_{\mathrm{f}})\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})}{{{\sum }_x}_{\mathrm{f}}{w}^2(x_{\mathrm{f}})\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})}·\hfill \end{array}$$(A.9)

From this equation we can see that the variance of the binned data is also not equal to the variance of the function at x_b, ${\sigma }_{\text{w}}^2(x_{\text{b}})\ne {\sigma }^2(x_{\text{b}})$, which complicates the calculation of covariance matrices for binned data.

To simplify covariance calculation we can set the weights in Eq. (A.6) equal to unity and obtain an unweighted rebinned estimate,

$$\begin{array}{c}\hfill {\widehat{f}}_{\mathrm{unweighted}}(x_{\mathrm{b}})=\frac{{{\sum }_x}_{\mathrm{f}}\widehat{f}(x_{\mathrm{f}})\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})}{{{\sum }_x}_{\mathrm{f}}\mathrm{\Delta }(x_{\mathrm{f}}-x_{\mathrm{b}})}.\end{array}$$(A.10)

In this case the expectation value of the estimator is,

$$\begin{array}{c}\hfill ⟨{\widehat{f}}_{\mathrm{unweighted}}(x_{\mathrm{b}})⟩=\frac{1}{2\mathrm{\Delta }x}{{\displaystyle {\int }_x}}_{\min }{(x_{\mathrm{b}})}^{x_{\max }(x_{\mathrm{b}})}\mathrm{d}xf(x),\end{array}$$(A.11)

where 2 Δx is the width of the bin. If the relative variation of the sampled points within a broad bin is large, then this estimator may not be optimal and can produce larger errors compared to the estimator in Eq. (A.6).

Equation (A.7) is useful for predicting the theoretical value of the binned function, especially when the sampling frequency of the data points, 𝒟(x), is derived from the data itself. In the case of cosmic shear ξ_±(θ) the sampling of the data points roughly scales with θ, however, survey geometry and masking effects together with variations in the depth of the images complicates the analytical estimation for the distribution of data in angular scale. Hence we suggest measuring 𝒟(θ) from the data and use Eq. (A.7) to predict the binned ξ_± values.

A.2. Application to cosmic shear

To demonstrate the level of bias introduced by partial treatment of the theory in a ξ_± cosmic shear analysis we use a theoretical prediction for ξ_± as our function, f(x), assuming a single KiDS-450 redshift bin. To sample ξ_±(θ) in a non-uniform way, we randomly pick θ values from a 𝒟(θ) = θ/arcmin × 2000 distribution in the angular range of [0.5′,300′]. We then add a constant Gaussian random noise with σ = 0.01 to each sampled point to produce the noisy sampled data points, ξ_±data(θ) and then bin ξ_±data(θ) into 1000 fine logarithmic bins to produce ${\widehat{\xi }}_{\pm \mathrm{f}}({\theta }_{\mathrm{f}})$ and 9 broad logarithmic bins to get ${\widehat{\xi }}_{\pm \mathrm{binned}}({\theta }_{\mathrm{b}})$ (see Eq. (A.1)).

Fig. A.1. 2PCF binning bias introduced for a range of analysis choices, shown as the ratio between the measured ${\widehat{\xi }}_{\pm }$ and their proposed theoretical value as a function of angular scale, θ. The legends in this figure are shared between the two panels. For the weighted broad bin estimator, ${\widehat{\xi }}_{\pm \mathrm{binned}}({\theta }_{\mathrm{b}})$, the bias is calculated assuming θb is given by the logarithmic mid-point of the bin, θmid,  log (green squares), the weighted mean of the bin, θw (black pluses, under the blue diamonds), the geometric mean or linear mid-point of the bin, θmid, lin (dark grey triangles), or an area-weighted bin centre, θarea (blue diamonds). These estimators can be compared to the fine binning case, ${\widehat{\xi }}_{\pm \mathrm{fine}}({\theta }_{\mathrm{mid},\mathrm{fine}})$, where the theory is estimated at the weighted mean of each bin (light grey circles), and the exact case (red crosses) where the theoretical value is calculated as a weighted integral over the signal within the bin (Eq. (A.7)). All points are plotted with errorbars, but in the case of broad binning the errors are too small to be visible.

In Fig. A.1 we show the binning bias introduced for a range of cases as a ratio between the measured and the proposed theoretical value of ξ₊ (top panel) and ξ₋ (bottom panel), as a function of angular scale, θ. The noisy finely binned data, ${\widehat{\xi }}_{\pm \mathrm{f}}(\theta )$ (light grey circles), shows no significant bias relative to its expectation value (see Eq. (A.5)). As was shown in Eqs. (A.3) and (A.7) the expectation value of the broad binned, ${\widehat{\xi }}_{\pm \mathrm{binned}}({\theta }_{\mathrm{b}})$, should be calculated using an integral over ξ_± with the appropriate weights. The red crosses in the figure correspond to this theoretical prediction which is unbiased as expected. The remaining curves show the biases introduced when broad binning is applied to the measurements and the theory is evaluated at a single point in the bin denoted as θ_b. The green squares assume that θ_b is given by the logarithmic mid-point of the bin (as used in Hildebrandt et al. 2017), for the black pluses θ_b is the weighted mean of the bin (as used in Heymans et al. 2013; Troxel et al. 2018a,b), the grey triangles use the geometric mean or linear mid-point of the bin (not used to date) and finally blue diamonds assume that θ_b is the area-weighted bin centre (as advocated by Krause et al. 2017), where

$$\begin{array}{c}\hfill {\theta }_{\mathrm{area}}=\frac{{\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\mathrm{d}\theta \theta 2\pi \theta }{{\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\mathrm{d}\theta 2\pi \theta }=\frac{2({\theta }_{\max }^3-{\theta }_{\min }^3)}{3({\theta }_{\max }^2-{\theta }_{\min }^2)}·\end{array}$$(A.12)

Here θ_min and θ_max are the minimum and maximum values of the bin.

We find that the weighted mid-point and the area weighted values are similarly biased, boosting the signal at the ∼3% level at 10 arcmin, rising to ∼10% bias at large scales for ξ₊. Taking the logarithmic mid-point of the bin has the opposite effect, decreasing the signal at ∼7% level at 10′ and ∼10% bias at large scales for ξ₊. That the biases work in the opposite sense here increases the inferred impact of binning bias when comparing the two KiDS analyses in Troxel et al. (2018b).

In all cases we see that the choice of binning affects ξ₊ more than ξ₋, since ξ₊ has more curvature than ξ₋. We note that these biases will be smaller for narrower angular bins and as such their effect will not be as significant for the first year DES analysis (Troxel et al. 2018a) which uses the weighted mean for θ_b with roughly twice as many bins in the same angular range as shown here.

If future surveys conclude that it is too computationally expensive to calculate the impact of binning theoretically, especially in the case of the covariance matrix, our proposed solution is to use the linear mid point of the θ-bin in the binned ξ_± analysis. We find that this approximation presents the weakest bias with at most 2.5% bias at large and small scales and below percent level bias between 0.5′< θ <  300′. Another alternative is to move to a COSEBIs analysis. As COSEBIs are discrete they are not subject to any of the binning biases presented in this Appendix.

Appendix B: σ₈ − Ω_m degeneracy

Fig. B.1. Degeneracy direction of σ8 and Ωm for a CCOSEBIs analysis of the KiDS-like data. The colours in the image show the value of the CCOSEBIs Eσ8 mode, in comparison to dashed lines of constant Σ8 = σ8(Ωm/0.3)α with α = 0.65. The repeating color scheme was chosen to capture the variations in the values of Eσ8. The lower left corner has the smallest value of Eσ8 which gradually increases, perpendicular to the dashed curves, towards the upper right corner.

Cosmic shear is most sensitive to a combination of σ₈ and Ω_m (Jain & Seljak 1997), where the degeneracy can be written as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_8={\sigma }_8{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{{\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{fid}}}\right)}^{\alpha }.\end{array}$$(B.1)

Here ${\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{fid}}$ is arbitrary but is usually taken to be 0.3. In the majority of cosmic shear analyses α has been taken to be α = 0.5, even though the optimal value of α will depend on the statistic used, the redshift distributions and the angular ranges used in the analysis. As an example, Hildebrandt et al. (2017) present joint Σ₈ − Ω_m constraints with α = 0.5. The tilt seen in their Fig. 6 of these constraints demonstrates that α = 0.5 does not best represent the degeneracy direction of Ω_m and σ₈ for the KiDS-450 2PCF tomographic analysis.

In Fig. B.1 we show the value of the CCOSEBIs mode E_σ8 (see colour bar) for a range of σ₈ and Ω_m values assuming a KiDS-like survey. The degeneracy shown in E_σ8 can be compared to the dashed lines of constant Σ₈ = σ₈(Ω_m/0.3)^α where α = 0.65. We have carried out this test for all the CCOSEBIs modes in our analysis; E_Ωm,E_σ8, σ8, E_σ8, Ωm and E_Ωm, Ωm to confirm that α = 0.65 is an optimal choice for our CCOSEBIs analysis.

Appendix C: Optimising the COSEBIs B-mode null-test

All null tests are subject to the choices we make in our data analysis. As an example, if we limit our B-mode analysis of CFHTLenS to the first 7 COSEBIs modes, following Asgari et al. (2017), we conclude there are no significant small scale B-modes in CFHTLenS. In contrast, our analysis of the first 20 COSEBIs modes, in Sect. 4, finds a significant B-mode detection for CFHTLenS on the same scales. In this Appendix we explore the question of how many COSEBIs modes should be used to determine the overall significance of the B-modes in a dataset.

Fig. C.1. Model comparison using 10 000 random samples of the PSF-leakage systematic model, MPSF (solid) defined in Sect. 5.1.1. The data points show the mean of the random samples, with the errors reflecting the noise on a single realisation. The data samples are analysed to determine how often the input model MPSF can be distinguished from the no-systematics zero B-mode model, M0 (dashed).

As an illustrative example, we take two parameter-free models for COSEBIs B-modes, shown in Fig. C.1: M₀ where B_n = 0 for all n, and M_PSF where B_n corresponds to the measured PSF-leakage systematic defined in Sect. 5.1.1. The difference between these two models is captured by the first few modes, with almost zero power for n ≳ 10. We create 10 000 random samples of B_n for the full angular range of [0.5′,100′] given the model M_PSF and the KiDS noise-only covariance for the non-tomographic case. Figure C.1 shows the mean of these samples (red squares) with errors corresponding to a single sample as well as the input model (blue curve).

Fig. C.2. χ2 distribution of the mock data given the true PSF leakage model, MPSF, (blue histogram) or given the null model, M0, (orange histogram). Left panel: analysis of the n ≤ 20 COSEBIs modes. In the right panel only the n ≤ 5 modes are considered.

Fig. C.3. Distribution of p-values for the 10 000 data samples, showing the probability that the MPSF model (blue histogram) or the M0 no systematics model (orange histogram) is the true underlying model, given each data sample. Left panel: p-values from an analysis of the n ≤ 20 COSEBIs modes. In the right panel only the n ≤ 5 modes are considered.

We can determine which of the two models best represents the data using a Bayesian evidence analysis. If we give the same weight to both models then the ratio of the Bayesian evidences for these models is given by the Bayes factor,

$$\begin{array}{c}\hfill \mathrm{Bayes}\mathrm{Factor}=\frac{P(D|{\mathsf{M}}_1)}{P(D|{\mathsf{M}}_2)}=\frac{\int P(D|{\mathsf{M}}_1,{\mathrm{\Phi }}_1)P({\mathrm{\Phi }}_1|{\mathsf{M}}_1)\mathrm{d}{\mathrm{\Phi }}_1}{\int P(D|{\mathsf{M}}_2,{\mathrm{\Phi }}_2)P({\mathrm{\Phi }}_2|{\mathsf{M}}_2)\mathrm{d}{\mathrm{\Phi }}_2},\end{array}$$(C.1)

where D is the data, M_i is model i and Φ_i represents the set of parameters for model i. For the simplified case of parameter free models that we consider here, Eq. (C.1) simplifies to,

$$\begin{array}{c}\hfill \mathrm{Bayes}\mathrm{Factor}=\frac{P(D|{\mathsf{M}}_0)}{P(D|{\mathsf{M}}_{\mathrm{PSF}})}·\end{array}$$(C.2)

The resulting Bayes factor will however depend on the number of n-modes that are included in the analysis. The Bayesian evidence can only be used when an alternative model exists, but in the case of null tests, such as a B-mode test, the only available model is the null hypothesis and therefore we need to use classical methods to identify the significance of the B-modes. Here we use χ² and p-values to test the null hypothesis. The p-value for the χ² is defined as the probability of calculating a χ² value larger than the measured one, ${\chi }_{\text{m}}^2$, given the model M,

$$\begin{array}{c}\hfill p-\mathrm{value}=\mathrm{Pr}({\chi }^2>{\chi }_{\mathrm{m}}^2|\mathsf{M})={\displaystyle {\int }_{{\chi }_{\mathrm{m}}^2}^{\infty }\mathrm{d}4{\chi }^2}\mathrm{Pr}({\chi }^2|\mathsf{M}).\end{array}$$(C.3)

Figure C.2 shows the distribution of the measured χ² across our 10 000 random samples when the data are fit using the input M_PSF model (blue histogram) and the M₀ no systematics model (orange histogram). In the left panel we take the null-test case where all modes up to n = 20 are included in the analysis (All-modes). In the right panel, only the first 5 modes (n <  6) are analysed. As M_PSF is the correct model, we naturally find better fits to the data, i.e. lower χ² values, for this model. The difference between the two distributions for the χ² values is however enhanced when the modes analysed are limited to the range where the two models differ significantly. This means that the power of the null test is optimised over this reduced, n ≤ 5, range.

Figure C.3 shows the distribution of p-values for the χ² values shown in Fig. C.2. If the model used to fit the data is the true underlying model, any particular p-value is as likely to be measured as the other. If the model is not representative of the data, however, then one is more likely to obtain smaller p-values from the sample. As expected with M_PSF as the correct model, we find a uniform distribution of p-values and a skewed distribution for the M₀ model. When all 20 COSEBIs modes are included this p-value distribution is less skewed compared to when we only include the n ≤ 5 modes. By adding more data points to the analysis, we have diluted the systematic signal of the PSF leakage, making this null-test less effective.

Based on this analysis, we must recognise that finding that the B-modes pass a null-test using a large data vector does not ensure that analysing a smaller dataset will give the same result. A good example of this is KiDS-450 passing the 20-mode COSEBIs null-test, but failing the CCOSEBIs null-test which is most sensitive to the n ≤ 5 modes. In contrast DES-SV and CFHTLenS fail the 20-mode COSEBIs null-test, even though they pass the CCOSEBIs null-test. Their B-modes therefore appear when adding in more data points to the analysis. As our example shows how increasing the size of your data set serves to reduce the stringency of the null-test, we can therefore conclude that the significant DES-SV and CFHTLenS B-mode, seen with COSEBIs and not with CCOSEBIs, is present in the high-n data that is not included in the CCOSEBIs analysis. If we had only performed a COSEBIs null-test, we would have missed the presence of a systematic signal in KiDS. If we had only performed a low-n CCOSEBIs null-test, we would have missed the presence of a systematic signal in DES-SV and CFHTLenS.

To illustrate our discussion of null-tests we have used COSEBIs, but the concept holds for any statistic or null-test. If a systematic produces a feature at a particular scale, but is otherwise identical to the standard model, by adding data from other scales we will dilute the power of the statistical test to distinguish between the two cases. As null B-mode tests are generally performed independently of alternative models, it is not clear which data points should be added to the null-test analysis. We therefore propose that future null-tests are performed with the B-mode signatures shown in Sect. 6 in mind. In this way one can optimise the modes over which to carry out a model comparison.

Appendix D: Supplementary data and figures

Figures D.1–D.3 show the tomographic COSEBIs measurements, using the angular range of [0.5′,100′], for DES-SV, KiDS-450 and CFHTLenS respectively. In each figure, the upper panels present the E-modes, the lower panels present the B-modes, and the significance of the B-modes are indicated with a p-value shown in the upper left corner. The p-values for the other two angular ranges analysed are given in Table 2. The predicted E-modes, given the best-fitting cosmology parameters listed in Table 1, are shown as curves.

Figures D.4–D.6 show ξ_E/B for the tomographic cases for DES-SV, KiDS-450 and CFHTLenS respectively. We show p-values for the significance of the B-modes in each figure, but caution the reader that due to binning and the truncated integrals discussed in Sect. 4, this method is not robust. However, as ξ_B data points are uncorrelated, they can help with identifying the source of the systematic even though it was seen in Sect. 6 that systematics do not always affect the same angular ranges for E and B-modes. The prediction for ξ_E, given the best-fitting cosmology parameters listed in Table 1, is shown as curves.