Free Access
Issue
A&A
Volume 540, April 2012
Article Number A34
Number of page(s) 12
Section Astronomical instrumentation
DOI https://doi.org/10.1051/0004-6361/201117129
Published online 20 March 2012

© ESO, 2012

1. Introduction

Weak gravitational lensing provides a method of directly mapping the distribution of dark matter in the universe (Mellier 1999; Bartelmann & Schneider 2001; Van Waerbeke et al. 2001; Mellier 2002; Refregier 2003; Munshi et al. 2008; Pires et al. 2010). This method is based on the weak distortions that lensing induces in the images of background galaxies as it travels toward us through intervening structures. The induced shear is small, typically around an order of magnitude less than the ellipticity seen in background galaxies, therefore the weak lensing effect cannot be measured on a single galaxy. To measure the shear field, it is necessary to measure ellipticities of many background galaxies and construct a statistical estimate of their systematic alignment. This measured shear field is directly related to the distribution of the (dark) matter and thus to the geometry and the dynamics of the universe. As a consequence, weak gravitational lensing offers unique possibilities for probing the statistical properties of dark matter and dark energy in the universe.

The mapping of the dark matter distribution has become a central topic in weak lensing since the very first reconstructed 2D mass maps demonstrated that you can see the dark components of the universe. Recently, this field has seen great success in reconstructing the deepest and largest 2D dark matter distribution in the COSMOS field (Massey et al. 2007). From this distribution, the nature of dark matter can be better understood, so better constraints can be set on dark energy (Schrabback et al. 2010). This method is now widely used but, the amplitude of the weak lensing signal is so weak that its detection relies on the accuracy of the techniques used to analyze the data. Each step in the analysis has required the development of advanced techniques dedicated to these applications.

The reconstruction of the dark matter mass map from shear measurements is a difficult inverse problem because of observational effects such as noise and the complex geometry of the field. Classical methods based on Kaiser & Squires (1993) are not optimal. Indeed they require a convolution of the observed shear with a kernel to be performed over the entire field, usually expressed in the Fourier domain, and the FFT imposes a condition of periodicity that is not true and introduces a mixing between E and B modes. An alternative interesting method was proposed in Seitz & Schneider (1996), which consists in convolving the shear field with a local kernel that depends on the shape of the domain. It was shown to be optimal in Lombardi & Bertin (1998), and it was reformulated later in Seitz & Schneider (2001).

This approach is relatively expensive in computational resources. Much effort has been made in past years to properly derive a cosmic shear two-point correlation on a finite interval (Kilbinger et al. 2006; Schneider & Kilbinger 2007; Fu & Kilbinger 2010). In Fu & Kilbinger (2010), a new statistic is proposed, maximizing the signal-to-noise ratio and a figure of merit based on the Fisher matrix of the cosmological parameters Ωm and σ8. However, this approach does not allow reconstructing a map, and cannot be used for higher order statistics studies such as bispectrum analysis (Van Waerbeke et al. 2001; Pires et al. 2009b; Munshi et al. 2011), peak counting (Marian et al. 2011), Minkowski Functionals (Kratochvil et al. 2011), or wavelet peak counting (Pires et al. 2009a).

In this paper, we study a new approach for weak lensing mass inversion in the ideal case of noise free data based on a Helmholtz decomposition in the wavelet domain. We first show that there is an explicit relation between the weak lensing E/B mode decomposition and the Helmholtz decomposition which decomposes a vector field into two components, one curl-free and the second divergence-free. In contrast to the method presented in Seitz & Schneider (1996), our method does not rely on a convolution with a kernel, and no intermediate field u derived from the shear field is needed. The Helmholtz decomposition is performed in the wavelet space, and the Fourier transform is never used. This gives us more flexibility to handle borders or to add some constraints in our solution. We introduced a new wavelet construction based on border-wavelets, and we present two kinds of wavelet Helmholtz decomposition algorithms. The first one imposes a global zero B modes constraint and the second a local zero B mode constraint only on the borders.

This paper is structured as follows. In Sect. 2, we review the basis of the weak lensing formalism. In Sect. 3, we introduce of the Helmholtz decomposition to the weak lensing E/B modes decomposition, and the discretization is discussed in Sect. 4. Then, in Sect. 5, the Helmholtz decomposition is extended to the wavelet domain and, we present new constructions of divergence-free and curl-free wavelets that satisfy suitable boundary conditions. And finally, in Sect. 6, we do some numerical experiments for weak lensing and conclude in Sect. 7. In Appendix A, we provide the mathematical expressions for the vector wavelets and in Appendix B, we present the algorithm to perform the E mode convergence reconstruction in a bounded domain using a wavelet Helmholtz decomposition.

2. Weak lensing formalism

2.1. The weak lensing equations

In weak lensing surveys, the observable shear field γi can be written in terms of the lensing gravitational potential ψ(θi) (Bartelmann & Schneider 2001): (1)where the partial derivatives i are with respect to θi. Notice the difference in sign with respect to (Seitz & Schneider 1996; Lombardi & Bertin 1998).

The shear γi(θ) with i = 1,2 is derived from the shapes of galaxies at positions θ in the image. The γ1 component corresponds to a deformation in the horizontal/vertical direction and the γ2 component to a deformation in the diagonal direction (see Fig. 1).

The convergence κ(θ) can also be expressed in terms of the lensing potential ψ, (2)and is related to the mass density Σ(θ) projected along the line of sight by (3)where the critical mass density Σcrit is given by (4)where G is Newton’s constant, c the speed of light, and DOS, DOL, and DLS are respectively the angular-diameter distances between the observer and the galaxies, the observer and the lens, and the lens and the galaxies.

2.2. Difference between the sheared galaxies and the shear field γ

The deformation of a circle by a shear field provides an ellipse and corresponds to the local linear application of the matrix (5)where κ is the convergence, φ the angle formed by the larger axis of the ellipse and the (Ox) axis, and |γ| the amplitude of the deformation.

thumbnail Fig. 1

Galaxy shapes and corresponding gravitational shears γ for a circular galaxy.

Open with DEXTER

In the weak lensing approximation, we only consider the shear γ = |γ|(cos   2φ,sin   2φ) which concentrates the information on the angle of deformation and the amplitude of deformation and defines a 2D vector field: γ:Ω → R2, see Fig. 1. Although the shear is often represented by sheared galaxies which is a spin-2 representation, as one can see in the first row of Fig. 1, the shear field γ itself is a usual vector field, as seen in the second row of Fig. 1, so we can compute its divergence and its curl: (6)and apply the Helmholtz decomposition to it.

2.3. The mass inversion

The weak lensing mass inversion problem consists in reconstructing the projected (normalized) mass distribution κ(θ) from the measured shear field γi(θ). We invert Eq. (1) to find the lensing potential ψ and then apply formula Eq. (2) following the classical methods based on the pioneering work of Kaiser & Squires (1993). In short, this corresponds to (Kaiser & Squires 1993; Starck et al. 2006): (7)This method usually is called a Fourier transform and has been generalized to spherical harmonics in Pichon et al. (2010) to deal with full-sky maps. To recover κ from both γ1 and γ2 there is a degeneracy when k1 = k2 = 0. Therefore, the mean value of κ cannot be recovered from the shear maps. This is known as the mass-sheet degeneracy.

The most important drawback of this method is that it requires an observation of the shear over the entire sky. As a result, if the observed shear field has a finite size or a complex geometry, then the method can produce artifacts on the reconstructed convergence distribution near the boundaries of the observed field. The error caused by the artificial periodic boundary conditions can not be decreased by taking mirror conditions which remove the discontinuity.

There are local inversions that reduce the unwanted boundary effects (Seitz & Schneider 1996). But regardless of the local formula used, the reconstructed field has more noise than when obtained with a global inversion, and the inversion is unstable in the presence of uncorrelated noise.

The convergence κ can be computed directly in direct space (without the Fourier transform) thanks to the kernel integration (Seitz & Schneider 1996): (8)where κ0 stands for the mean value of κ. The kernel K depends on the geometry of the domain Ω. For Ω = R2, it is given by (9)In Seitz & Schneider (1996), the authors propose to combine the derivatives of (γi) (10)and then to apply the Helmholtz decomposition u = ∇κ(E) + ∇ × κ(B), in order to reconstruct the convergence κ = κ(E).

In the following, we propose to use the Helmholtz decomposition directly on (γi) to perform the E/B mode decomposition, and this can be seen as a kind of shortcut in comparison with the Seitz & Schneider method. Furthermore, we also propose to perform the Helmholtz decomposition in the wavelet space using a new kind of border-wavelet, which allows us to reduce border effects during the reconstruction.

2.4. The E and B mode decomposition

Just as a vector field can be decomposed into a gradient or electric (E) component, and a curl or magnetic (B) component, the shear field γi(θ) can be decomposed into two components, which for convenience we also call (E) and (B). The decomposition of the shear field into each of these components can be easily performed by noticing that a pure E mode can be transformed into a pure B mode by a clockwise rotation of the sheared galaxies by 45°: γ1 →  − γ2, γ2 → γ1 (see Fig. 1).

Because the weak lensing arises from a scalar potential (the Newtonian potential Φ), it can be shown that weak lensing only produces E modes. As a result, the least square estimator of the E modes convergence field is simply (11)The E mode deformations γ(E) which satisfy Eq. (1) can be equivalently defined by (12)On the other hand, residual systematics arising from imperfect correction of the instrumental PSF or telescope aberrations, generally generates both E and B modes. The presence of B modes can thus be used to test for any residual systematic effects in current weak lensing surveys. For this purpose, the following estimator for the B mode convergence field can be formed (Starck et al. 2006), (13)and it is consistency with zero will confirm the absence of systematics.

The B mode deformations γ(B) derive from potentials ψB as (14)and form the set (15)If we consider γ defined on the infinite plan R2 or on the torus T2, i.e. γ ∈ (L2(R2))2 or γ ∈ (L2(T2))2, then ΓE and ΓB define orthogonal spaces for the usual scalar product in L2. This is no longer the case for domains with boundaries. In this case, ΓE and ΓB have nonzero common elements.

Solutions Eqs. (11) and (13) are therefore no longer valid for bounded domains. The gravitational shear that derives from the bi-harmonic function ψ such that Δ2ψ = 0 produces modes that are elements both of ΓE and ΓB (Bunn et al. 2003). Indeed, the borders cause mode-mixing or, equivalently, a leakage of E modes into B modes. Therefore, even a pure gravitational shear will produce both E and B modes in a finite field. Consequently, it will be more difficult to confirm the absence of residual systematics with the presence of border effects. More details can be found in Crittenden et al. (2002); Schneider et al. (2002).

In the next section, we propose to explore the use of the Helmholtz decomposition to reduce boundary effects in E and B mode decomposition.

3. Introduction of the Helmholtz decomposition into the weak lensing formalism

In this section, we introduce the relation between the Helmholtz decomposition and the weak lensing E/B mode decomposition. The Helmholtz decomposition in two orthogonal components only stands for unbounded domains. Then the unicity of the κE reconstruction on bounded domains can only be obtained under some constraints.

3.1. The weak lensing E/B mode decomposition

On (L2(R2))2, we can express the relations Eqs. (1) & (14) thanks to the Fourier transform: (16)It allows us to derive the convergence κE deriving from the E modes (the actual convergence) and the convergence κB deriving from the B modes (the systematic errors) similarly to Eqs. (11) and (13): (17)where hat symbols denote the Fourier transforms.

3.2. Relation between Helmholtz and E/B modes decompositions on an unbounded domain

The Helmholtz decomposition theorem says that every vector field γ, defined everywhere in space can be decomposed into a rotational part γdiv   0 and an irrotational part γcurl   0. The Helmholtz decomposition consists, for a vector field γ, in writing (18)where γdiv   0 = ∇ × ψ and γcurl   0 = ∇p with ψ and p differentiable functions.

This decomposition is unique for domains without boundaries. In the Fourier domain, it consists in decomposing the vector field in the longitudinal direction, i.e. parallel to k and in the transverse direction, i.e. perpendicular to k. The longitudinal component corresponds to the curl-free component and the transverse component to the divergence-free component.

Using the following notations: Pγ = γdiv   0 and Qγ = γcurl   0, we obtain (19)with: From Eqs. (17), (20), and (21), we have Aκ =  −AP + AQ. Then (22)This provides the following E and B modes convergence:

3.3. Relation between Helmholtz and E/B mode decompositions on a bounded domain

On the square Ω = (0,1)2, the Helmholtz decomposition becomes much more complicated, see (Dautray & Lions 1984, for further details). Indeed, the unicity of the Helmholtz decomposition γ = Pγ + Qγ disappears. However, it can be regained by considering three components in the orthogonal decomposition: (25)with div   γdiv   0 = 0 on the square Ω and γdiv   0·n = 0 on the borders of the square Ω, curl   γcurl   0 = 0 on Ω and γcurl   0·τ = 0 on Ω, and div   γharm = curl   γharm = 0 on Ω, where n and τ are the normal and tangential unit vectors on Ω. This new term γharm derives from a harmonic function on Ω: γharm = ∇p0 = ∇ × ω0 with Δp0 = Δω0 = 0. Alternatively we can write γharm = ∇p0 + ∇ × ω0 where p0 and ω0 are harmonic functions.

Considering the Helmholtz decomposition of γ = Pγ + Qγ, the new term γharm can be included either in the divergence-free part or in the curl-free part of the decomposition. This decomposition is no longer unique, and it is the same for the E/B mode decomposition (24).

However, the unicity of the reconstruction of κE can be obtained under the condition κB = 0. Indeed, if γ ∈ (L2( [0,1] 2)2 is split into a divergence-free part Pγ and a curl-free part Qγ, satisfying (26)then, γ is an element of ΓE as defined in Eq. (12). In addition, the convergence κE is uniquely defined, modulo an additive constant. Reciprocally, if γ is an element of ΓE, then we have the decomposition (27)and (28)so The same is true for the B modes. Under the constraint κE = 0, γ is an element of ΓB as defined in Eq. (15), and κB is uniquely defined, modulo an additive constant.

This method offers a way to recover the E mode convergence without the border effects. But the constraint κB = 0 should be used carefully, because it assumes that there are no residual systematics in the data due to the observational or instrumental effects. We see in the following that we can relax this constraint, and impose that κB = 0 only on the border, and have a good reconstruction of E modes, even in the presence of B modes.

4. Discretization of the problem

4.1. Using a harmonic decomposition

To solve the weak lensing problem on a bounded domain Ω using the formula Eq. (26), a possible solution consists in writing γ as (29)with (30)where p = 0 and ω = 0 on the borders of the square δΩ and (31)i.e., p0 and ω0 are harmonic functions on Ω.

From the relations Eqs. (24) and (29), we can rewrite the E/B modes decomposition as (32)This decomposition is not unique, and we have to find how to split the harmonic term γharm into ∇p0 and ∇ × ω0.

thumbnail Fig. 2

Example of two1D MRA related by differentiation and integration. Scaling functions and associated wavelets for biorthogonal MRA with 2 (left) and 3 (right) zero moments.

Open with DEXTER

By assuming γ derives from the gravitational potential without residual systematics, we can apply to the solution the constraint that κB = 0 on all the domain Ω. Then, the solution to this problem becomes unique and the convergence is given by Eq. (32).

Then, p0 and ω0 can be estimated accurately by minimizing the following quantity: (33)where ∇p and ∇ × ω can be obtained from the relation (30) and ∇p0 and ∇ × ω0 can be obtained by decomposing p0 and ω0 in a basis of harmonic functions on Ω ( and ) and minimizing the quantity (33). A set of harmonic functions ψn can be built as (34)Or, alternatively, (35)where ℛ and ℐ stand for the real and imaginary parts, and L is a constant such that 2π   L is larger than the size of the domain.

4.2. Using a wavelet decomposition

To solve the weak lensing problem on a bounded domain Ω, another possible solution consists in looking for a decomposition of γ of the form: (36)By considering a zero B modes constraint as in the previous section, (37)the problem can be solved by minimizing the following quantity: (38)This minimization problem may involve two bases and of the divergence-free and curl-free functions on Ω without any boundary conditions: (39)This is the option that has been considered in the paper. The basis functions and will consist of wavelets. The construction of these wavelets will be presented in the next section and in Appendix A. Only a few of the boundary elements and do not satisfy the conditions γdiv   0·n = 0 and γcurl   0·τ = 0 on the border Ω as in the decomposition (30). Then it is possible to restrain the zero B mode constraint to these border wavelets, and we call this the B-zero border constraint.

5. Construction of wavelets for Helmholtz decomposition

In this section, we first review briefly the construction of divergence-free and curl-free wavelets on an unbounded domain. In Appendix A, we present in detail the bases that we have selected. Then, we exhibit new divergence-free and curl-free wavelets on the square that satisfy suitable boundary conditions.

5.1. Unbounded domain

Divergence-free and curl-free wavelets

The compactly supported divergence-free wavelets in (L2(Rn))n originally developed by Lemarié-Rieusset (1992) are based on the existence of two different 1D multiresolution analyses (MRA) of L2(R) related by differentiation and integration. Given a 1D multiresolution analysis (ϕ1,ψ1) with ϕ1 ∈ C1 + ϵ for ϵ > 0, there is another 1D MRA (ϕ0,ψ0) such that (40)Different MRA divergence-free wavelets can be constructed from different families of 1D MRA’s. The construction of MRA curl-free wavelets follows the same approach. Figure 2 shows an example of a family of 1D MRA functions and Fig. 3 shows the corresponding 2D divergence-free wavelets. These wavelets form a basis of the space of the divergence-free functions on R2. Nevertheless, they are not suited to implementing the Helmholtz decomposition.

thumbnail Fig. 3

Example of a 2D divergence-free wavelet basis designed from the Lemarié & Rieusset (1992) scheme.

Open with DEXTER

thumbnail Fig. 4

Examples of hyperbolic 2D divergence-free and curl-free wavelets.

Open with DEXTER

For this purpose, we prefer to use the hyperbolic divergence-free and curl-free wavelets proposed by Deriaz & Perrier (2009), which form better conditioned bases. They are based on multidimensional hyperbolic wavelets introduced by DeVore et al. (1998). The advantage of these wavelets is that there is a fast algorithm to compute the Helmholtz decomposition in terms of divergence-free and curl-free wavelets.

The Helmholtz decomposition in the wavelet domain uses both the divergence-free and curl-free wavelets. As the divergence-free and curl-free wavelets form biorthogonal bases in their respective spaces, the coefficients and corresponding to the divergence-free and curl-free parts (Eq. (39)) are, in practice, computed with an iterative algorithm (see Deriaz & Perrier 2006, for more details).

5.2. Bounded domain

The Helmholtz decomposition (39) may be generalized to bounded domains with non-periodic boundary conditions by the construction of divergence-free and curl-free wavelets on the interval. Since divergence-free and curl-free wavelets are constructed from standard, compactly supported, biorthogonal wavelet bases, they can incorporate boundary conditions. The following section introduces a new construction of wavelets well suited to the Helmholtz decomposition on bounded domains.

5.2.1. Wavelet transforms on the interval

We obtain wavelets on the square by taking tensor products of 1D wavelets on the interval. The wavelets on the interval we introduce here differ from the original constructions (Meyer 1990; Cohen et al. 1993) and rely on the lifting scheme (Daubechies & Sweldens 1998).

The properties of these MRA on the interval allow us to apply the wavelet Helmholtz algorithm efficiently on the square: the partition of the unity property facilitates the approximation of the data, and that only ϕ(1) and ϕ(N) differ from 0 at the boundaries permits us to impose boundary conditions on the wavelets.

thumbnail Fig. 5

Example of divergence-free wavelets with a free (left) and a sliding boundary condition (right).

Open with DEXTER

An example of spline divergence-free wavelets issued from this construction is shown in Fig. 5 for free boundary (left) and for sliding conditions (right). The boundary curl-free wavelets can be constructed similarly. We applied the resulting boundary divergence-free and curl-free wavelets to perform a Helmholtz decomposition on the bounded domain. In the next section, we present the results obtained with the algorithms that are derived from this construction and presented in Appendix B.

6. Numerical experiment for weak lensing

6.1. Simulated data

We have run realistic simulated convergence mass maps derived from N-body cosmological simulations using the RAMSES code (Teyssier 2002). The cosmological model is chosen to be in concordance with the ΛCDM model. We chose a model with the following parameters close to WMAP: Ωm = 0.3, σ8 = 0.9, ΩL = 0.7, h = 0.7. Each simulation has 2563 particles with a box size of 162 Mpc/h. The resulting convergence map covers a field of 2°    ×    2° with 350  ×  350 pixels and assumes a galaxy redshift of 1. The overdensities correspond to the halos of groups and clusters of galaxies.

6.2. Example with no B modes

In this experiment, we computed the shear map γS from the simulated map κS, and we have extracted a sub area γD from γS. The lefthand side of Fig. 6 shows the simulated κS mass map and its corresponding shear field. The righthand side of Fig. 6 shows the subfield γD which is used for the reconstruction. Then we reconstructed and from γD, and compared to κT, where κT is the corresponding original area of κS. This allows us to study how well the reconstructing methods perform with the border problem.

thumbnail Fig. 6

Left, simulated κS image and its corresponding shear field γS. Right, extracted subfield γD. The extraction breaks the periodicity of the data.

Open with DEXTER

Figure 7 shows the error maps relative to four different reconstruction methods: i) the standard FFT approach, ii) Seitz & Schneider method, iii) the wavelet Helmholtz (50 iterations) with a zero B-mode border constraint, and iv) the wavelet Helmholtz (50 iterations) with a zero B-mode constraint. The error maps were obtained by .

We calculated the relative percentage error on both E and B modes: and , and obtain ϵE = 28%, ϵE = 5.3%, ϵE = 1.2% and ϵE = 1.0% for the four methods respectively. For the B mode component, we obtained ϵB = 30.2% for FFT approach and ϵB = 1.3 % for the wavelet Helmholtz with a zero B modes border constraint. Figure 8 shows the reconstructed B-mode map for the Fourier method and the wavelet Helmholtz with a zero B modes border constraint.

thumbnail Fig. 7

Reconstruction error maps for four different methods: with the Fourier transform (top left), Seitz & Schneider (top right), with the constrained Helmholtz decomposition with the constraint only on the boundary (bottom left), and on all the domain (bottom right).

Open with DEXTER

thumbnail Fig. 8

B-mode reconstructed maps.

Open with DEXTER

thumbnail Fig. 9

Maps of the reconstruction error on κE for three different methods when B modes contaminate the data.

Open with DEXTER

This error is not due to the discontinuity on the boundary because we tested the mirroring technique, which consists in symmetrizing the data along the (Ox) and (Oy) axes in order to recover the continuity in exchange to a doubling of the size of the domain. It produced almost the same error in terms of norm as a straightforward Fourier transform.

From this experiment, we can conclude that the wavelet Helmholtz method allows us to dramatically reduce the error relative to the border effect, when no B modes contaminate our data.

6.3. Example with B modes

In this experiment, we have added B modes. The B-mode map consists in a Gaussian with a maximum equal to 10% of the maximum of κT from the mass map. Since the constrained wavelet Helmholtz decomposition assumes no B modes, or zero-border B modes, it is interesting to see what happens when we introduce some actual B modes into the model.

Figure 9 shows the reconstruction error maps with the three methods. We can see that the reconstruction error when using Fourier is not amplified in the presence of B modes, while the error increases when using the wavelet Helmholtz with a constraint on the B modes. This is expected since the model is not completely true anymore. We note, however, that even if the error increases, it is still much below the FFT error. In this case, wavelet Helmholtz with zero border B modes is more robust than wavelet Helmholtz with zero B modes.

thumbnail Fig. 10

Original B mode and reconstructed error maps for the B modes map κB.

Open with DEXTER

Figure 10 shows the original B mode map and the reconstructed error maps for both the Fourier transform and the wavelet Helmholtz decomposition with a zero B-mode border constraint. From this experiment, we can conclude that the wavelet Helmholtz method with zero-border B-mode constraint is the best trade-off when some residual B modes may contaminate the data.

7. Conclusion

We have sorted out a direct and original relationship between the Helmholtz decomposition and the E/B mode reconstruction. The Helmholtz decomposition can be performed without the FFT, by using divergence-free and curl-free wavelets. This wavelet Helmholtz decomposition is iterative, and we have shown that we have different strategies to handle the border, thanks to the versatility of wavelets. This means that, periodic borders are not a curse anymore. Indeed, we were able to design two kinds of Helmholtz decompositions, one imposing B modes as equal to 0, and another imposing B modes as zero only on the border. We have shown that both approaches outperform the FFT-based standard E/B mode reconstruction, which suffers from the serious drawback of imposing periodic border conditions that concentrate in the low frequencies. The smaller the field, the more important the bias introduced by this periodic border condition. In our experiments, using the Helmholtz wavelet decomposition reduces dramatically the error. We have shown that the wavelet Helmholtz decomposition outperforms also Seitz and Schneider method (Seitz & Schneider 1996). Our approach also offers the advantage of no needing to pixelate the shear map. Indeed, both divergence-free and curl-free wavelet transforms can be applied directly on the shear catalog without having to pixelate first on a uniform grid.

However, at this stage, our wavelet Helmholtz cannot be properly used yet on real data. Indeed, two main issues have not been addressed in this paper, the noise and the missing data. Both may be problematic: stability, convergence, and unicity of the solution must be studied in detail. This will be done in the future.

Appendix A: Divergence-free wavelets on a square domain

Expression of the divergence-free and curl-free wavelets

thumbnail Fig. A.1

Quadratic spline scaling functions adapted to the left boundary of the interval.

Open with DEXTER

thumbnail Fig. A.2

Linear spline scaling functions adapted to the interval.

Open with DEXTER

With the 1D wavelets of Eq. (40), the 2D isotropic divergence-free scaling function is defined as (Lemarié-Rieusset 1992) (41)and the details on each scale are obtained from the three wavelets

  • Horizontal wavelet: ;

  • Vertical wavelet: ;

  • Diagonal wavelet: .

The hyperbolic divergence-free vector wavelets are expressed as

  • Divergence-free wavelets: (42)

  • Curl-free wavelets: (43)with the scale parameter j = (j1,j2) ∈ Z2 and the position parameter k = (k1,k2) ∈ Z2. The wavelet is localized at (2 − j1k1,2 − j2k2).

The fast wavelet decomposition in these vector wavelet bases proceeds as follows. First, in 2D, we need to define two differentiable wavelet bases (Deriaz & Perrier 2006): (44)We decompose the field γ in these bases. Then, as the divergence-free and curl-free wavelets are defined by , and . The vector wavelet expansion is obtained by a change of basis (cf. Deriaz & Perrier 2006).

One can notice that the divergence-free wavelets correspond to the rotational of the wavelet basis and the curl-free wavelets correspond simply to its gradient.

Multi-Resolution Analyses on the interval linked by derivation

We extend this construction to bounded domains using special 1D wavelets on the interval. To obtain the MRAs on  [0,1] , we modify the filters of multiresolution analyses on R on the boundaries. The associated scaling functions (father wavelets) form a partition of unity in the sense that their sum is equal to one. For instance the quadratic spline filter, (45)is modified on the left boundary to (46)and (47)We thus we obtain the corresponding scaling functions plotted in Fig. A.1.

To form the divergence-free and curl-free wavelets, we need two MRAs linked by derivation. In the following, we present a method to integrate an MRA that was partially derived from Stevenson (2010).

We start from a Multi-Resolution Analysis (MRA) on the interval  [0,N] , verifying (48)We put (49)Then we define the MRA by (50)(51)for      i = 2,...,M, (52)Then the set forms a partition of the unity, (53)and satisfies . Thanks to this construction, to any wavelet defined in the MRA formed by we can associate a wavelet in the MRA formed by by the relation .

In the example of Sect. 5.2.1, the corresponding MRA are the linear splines on the interval. The filter is given by (54)and the left boundary filter by (55)

Divergence-free wavelets on a bounded domain

From the differentiation relation introduced above (56)we form the boundary divergence-free wavelets: (57)For these boundary divergence-free wavelets, the border conditions are given by , perpendicular to the border and by , tangential to the border.

Appendix B: Algorithms

In this appendix, we present a new algorithm based on the wavelet Helmholtz decomposition that improves the weak lensing E/B mode decomposition.

Wavelet Helmholtz decomposition algorithm

The Helmholtz decomposition may be viewed as the following orthogonal space splitting: (58)where Hdiv   0(R2) is the space of divergence-free vector functions and Hcurl   0(R2) is the space of curl-free vector functions. According to the wavelet decomposition (Deriaz & Perrier 2006), we have the following direct sums of MRA spaces: (59)These correspond to the decompositions in the divergence-free and curl-free wavelet bases, and, denoting Pdiv   0 as the non orthogonal wavelet projector on Hdiv   0 – see (Deriaz & Perrier 2006) for its definition – we have (60)where and are complementary projectors of Pdiv   0 and Qcurl   0, respectively.

Unfortunately, the projectors on the divergence-free wavelet and curl-free wavelet bases are biorthogonal projectors, so we need an iterative method to provide γdiv   0 and γcurl   0 and to write (61)Since the wavelet bases are conditioned well, it is possible to apply them recursively in order to construct sequences and , which converge to γdiv   0 and γcurl   0 in a straight way. Then, in order to decompose the data γ into its divergence-free and curl-free parts , we have implemented the following iterative algorithm:

  1. Set u0 = v0 = 0 and the maximum number of iterations to Imax.

  2. Set i = 0, while ||γ − ui − vi||2 > ϵ and i < Imax then iterate:

    1. Set

    2. Set

    3. Set

    4. Set

    5. Set i = i + 1.

The algorithm converges when the residual goes to 0 and has been proven numerically. We tested this algorithm with periodic boundary conditions and obtained almost the same result as with the Fourier transform. Another point is that this algorithm still converges on a bounded domain with free boundary condition wavelets (see the left panel of Fig. 5).

Weak lensing reconstruction on a bounded domain

In this section, we present the iterative algorithm that aims at adding the constraint κB = 0 in the wavelet Helmholtz decomposition algorithm to improve the reconstruction of the E mode convergence κE in a bounded domain: (62)The algorithm that we have implemented is

  1. Set u0 = v0 = 0 and the maximum number of iterations Imax.

  2. Set i = 0, while ||γ − ui − vi||2 > ϵ and i < Imax then iterate:

    1. Set

    2. Set

    3. Set

    4. Set

    5. Set

    6. Set

    7. Set

    8. Set

    9. Set i = i + 1.

This algorithm converges with a slow decay rate.

An alternative to this algorithm consists in restricting the B modes constraint only to the boundary wavelets . With such an approach, we are able to isolate some of the B modes that are not harmonic. The alternative algorithm will be

  1. Set u0 = v0 = 0 and the maximum number of iterations to Imax.

  2. Set i = 0, while ||γ − ui − vi||2 > ϵ and i < Imax then iterate:

    1. Set

    2. Set with:

    3. Compute

    4. Set

    5. Set

    6. Set

    7. Set with

    8. Compute

    9. Set

    10. Set .

Acknowledgments

The authors thank Peter Schneider for kindly providing his code, and Adrienne Leonard for her help in using it. This work was supported by the French National Agency for Research (ANR -08-EMER-009-01) and the European Research Council grant SparseAstro (ERC-228261).

References

All Figures

thumbnail Fig. 1

Galaxy shapes and corresponding gravitational shears γ for a circular galaxy.

Open with DEXTER
In the text
thumbnail Fig. 2

Example of two1D MRA related by differentiation and integration. Scaling functions and associated wavelets for biorthogonal MRA with 2 (left) and 3 (right) zero moments.

Open with DEXTER
In the text
thumbnail Fig. 3

Example of a 2D divergence-free wavelet basis designed from the Lemarié & Rieusset (1992) scheme.

Open with DEXTER
In the text
thumbnail Fig. 4

Examples of hyperbolic 2D divergence-free and curl-free wavelets.

Open with DEXTER
In the text
thumbnail Fig. 5

Example of divergence-free wavelets with a free (left) and a sliding boundary condition (right).

Open with DEXTER
In the text
thumbnail Fig. 6

Left, simulated κS image and its corresponding shear field γS. Right, extracted subfield γD. The extraction breaks the periodicity of the data.

Open with DEXTER
In the text
thumbnail Fig. 7

Reconstruction error maps for four different methods: with the Fourier transform (top left), Seitz & Schneider (top right), with the constrained Helmholtz decomposition with the constraint only on the boundary (bottom left), and on all the domain (bottom right).

Open with DEXTER
In the text
thumbnail Fig. 8

B-mode reconstructed maps.

Open with DEXTER
In the text
thumbnail Fig. 9

Maps of the reconstruction error on κE for three different methods when B modes contaminate the data.

Open with DEXTER
In the text
thumbnail Fig. 10

Original B mode and reconstructed error maps for the B modes map κB.

Open with DEXTER
In the text
thumbnail Fig. A.1

Quadratic spline scaling functions adapted to the left boundary of the interval.

Open with DEXTER
In the text
thumbnail Fig. A.2

Linear spline scaling functions adapted to the interval.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.