Issue 
A&A
Volume 540, April 2012



Article Number  A34  
Number of page(s)  12  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/201117129  
Published online  20 March 2012 
Wavelet Helmholtz decomposition for weak lensing mass map reconstruction
^{1} Laboratoire de Mécanique Modélisation et Procédés Propres – UMR6181 CNRS, IMT La Jetée, Technopôle de ChâteauGombert, 38, rue Frédéric JoliotCurie, 13451 Marseille Cedex 20, France
email: erwan.deriaz@l3m.univmrs.fr
^{2} Laboratoire AIM, UMR CEACNRSParis 7, Irfu, SEDISAP, Service d’Astrophysique, CEA Saclay, 91191 GifsurYvette Cedex, France
Received: 22 April 2011
Accepted: 6 January 2012
To derive the convergence field from the gravitational shear γ of the background galaxy images, the classical methods require a convolution of the shear to be performed over the entire sky, usually expressed by the fast Fourier transform (FFT). However, it is not optimal for an imperfect geometry survey. Furthermore, FFT implicitly uses periodic conditions that introduce errors into the reconstruction. A method has been proposed that relies on computation of an intermediate field u that combines the derivatives of γ and on convolution with a Green kernel. In this paper, we study the wavelet Helmholtz decomposition as a new approach to reconstructing the dark matter mass map. We show that a link exists between the Helmholtz decomposition and the electric and magnetic component separation. We introduce a new wavelet construction that has a property that gives us more flexibility in handling the border problem, and we propose a new method of reconstructing the dark matter mass map in the wavelet space. A set of experiments based on noisefree images illustrates that this Wavelet Helmholtz decomposition reconstructs the borders better than all other existing methods.
Key words: methods: data analysis / gravitational lensing: weak
© 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 twopoint 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 signaltonoise 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 curlfree and the second divergencefree. 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 borderwavelets, 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 divergencefree and curlfree 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 D_{OS}, D_{OL}, and D_{LS} are respectively the angulardiameter 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.
Fig. 1 Galaxy shapes and corresponding gravitational shears γ for a circular galaxy. 
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: γ:Ω → R^{2}, see Fig. 1. Although the shear is often represented by sheared galaxies which is a spin2 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 fullsky maps. To recover κ from both γ_{1} and γ_{2} there is a degeneracy when k_{1} = k_{2} = 0. Therefore, the mean value of κ cannot be recovered from the shear maps. This is known as the masssheet 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 Ω = R^{2}, 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 borderwavelet, 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 R^{2} or on the torus T^{2}, i.e. γ ∈ (L^{2}(R^{2}))^{2} or γ ∈ (L^{2}(T^{2}))^{2}, then Γ_{E} and Γ_{B} define orthogonal spaces for the usual scalar product in L^{2}. 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 biharmonic function ψ such that Δ^{2}ψ = 0 produces modes that are elements both of Γ_{E} and Γ_{B} (Bunn et al. 2003). Indeed, the borders cause modemixing 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 (L^{2}(R^{2}))^{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 curlfree component and the transverse component to the divergencefree 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_{κ} = −A_{P} + A_{Q}. 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} = ∇p_{0} = ∇ × ω_{0} with Δp_{0} = Δω_{0} = 0. Alternatively we can write γ_{harm} = ∇p_{0} + ∇ × ω_{0} where p_{0} and ω_{0} are harmonic functions.
Considering the Helmholtz decomposition of γ = Pγ + Qγ, the new term γ_{harm} can be included either in the divergencefree part or in the curlfree 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 γ ∈ (L^{2}( [0,1] ^{2})^{2} is split into a divergencefree part Pγ and a curlfree 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., p_{0} 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 ∇p_{0} and ∇ × ω_{0}.
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. 
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, p_{0} and ω_{0} can be estimated accurately by minimizing the following quantity: (33)where ∇p and ∇ × ω can be obtained from the relation (30) and ∇p_{0} and ∇ × ω_{0} can be obtained by decomposing p_{0} 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 divergencefree and curlfree 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 Bzero border constraint.
5. Construction of wavelets for Helmholtz decomposition
In this section, we first review briefly the construction of divergencefree and curlfree wavelets on an unbounded domain. In Appendix A, we present in detail the bases that we have selected. Then, we exhibit new divergencefree and curlfree wavelets on the square that satisfy suitable boundary conditions.
5.1. Unbounded domain
Divergencefree and curlfree wavelets
The compactly supported divergencefree wavelets in (L^{2}(R^{n}))^{n} originally developed by LemariéRieusset (1992) are based on the existence of two different 1D multiresolution analyses (MRA) of L^{2}(R) related by differentiation and integration. Given a 1D multiresolution analysis (ϕ_{1},ψ_{1}) with ϕ_{1} ∈ C^{1 + ϵ} for ϵ > 0, there is another 1D MRA (ϕ_{0},ψ_{0}) such that (40)Different MRA divergencefree wavelets can be constructed from different families of 1D MRA’s. The construction of MRA curlfree wavelets follows the same approach. Figure 2 shows an example of a family of 1D MRA functions and Fig. 3 shows the corresponding 2D divergencefree wavelets. These wavelets form a basis of the space of the divergencefree functions on R^{2}. Nevertheless, they are not suited to implementing the Helmholtz decomposition.
Fig. 3 Example of a 2D divergencefree wavelet basis designed from the Lemarié & Rieusset (1992) scheme. 
Fig. 4 Examples of hyperbolic 2D divergencefree and curlfree wavelets. 
For this purpose, we prefer to use the hyperbolic divergencefree and curlfree 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 divergencefree and curlfree wavelets.
The Helmholtz decomposition in the wavelet domain uses both the divergencefree and curlfree wavelets. As the divergencefree and curlfree wavelets form biorthogonal bases in their respective spaces, the coefficients and corresponding to the divergencefree and curlfree 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 nonperiodic boundary conditions by the construction of divergencefree and curlfree wavelets on the interval. Since divergencefree and curlfree 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.
Fig. 5 Example of divergencefree wavelets with a free (left) and a sliding boundary condition (right). 
An example of spline divergencefree wavelets issued from this construction is shown in Fig. 5 for free boundary (left) and for sliding conditions (right). The boundary curlfree wavelets can be constructed similarly. We applied the resulting boundary divergencefree and curlfree 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 Nbody 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 256^{3} 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.
Fig. 6 Left, simulated κ^{S} image and its corresponding shear field γ^{S}. Right, extracted subfield γ^{D}. The extraction breaks the periodicity of the data. 
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 Bmode border constraint, and iv) the wavelet Helmholtz (50 iterations) with a zero Bmode 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 Bmode map for the Fourier method and the wavelet Helmholtz with a zero B modes border constraint.
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). 
Fig. 8 Bmode reconstructed maps. 
Fig. 9 Maps of the reconstruction error on κ_{E} for three different methods when B modes contaminate the data. 
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 Bmode 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 zeroborder 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.
Fig. 10 Original B mode and reconstructed error maps for the B modes map κ_{B}. 
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 Bmode border constraint. From this experiment, we can conclude that the wavelet Helmholtz method with zeroborder Bmode constraint is the best tradeoff 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 divergencefree and curlfree 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 FFTbased 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 divergencefree and curlfree 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: Divergencefree wavelets on a square domain
Expression of the divergencefree and curlfree wavelets
Fig. A.1 Quadratic spline scaling functions adapted to the left boundary of the interval. 
Fig. A.2 Linear spline scaling functions adapted to the interval. 
With the 1D wavelets of Eq. (40), the 2D isotropic divergencefree 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: .

Curlfree wavelets: (43)with the scale parameter j = (j_{1},j_{2}) ∈ Z^{2} and the position parameter k = (k_{1},k_{2}) ∈ Z^{2}. The wavelet is localized at (2^{ − j1}k_{1},2^{ − j2}k_{2}).
One can notice that the divergencefree wavelets correspond to the rotational of the wavelet basis and the curlfree wavelets correspond simply to its gradient.
MultiResolution 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 divergencefree and curlfree 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 MultiResolution 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)
Divergencefree wavelets on a bounded domain
From the differentiation relation introduced above (56)we form the boundary divergencefree wavelets: (57)For these boundary divergencefree 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 H_{div 0}(R^{2}) is the space of divergencefree vector functions and H_{curl 0}(R^{2}) is the space of curlfree 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 divergencefree and curlfree wavelet bases, and, denoting P_{div 0} as the non orthogonal wavelet projector on H_{div 0} – see (Deriaz & Perrier 2006) for its definition – we have (60)where and are complementary projectors of P_{div 0} and Q_{curl 0}, respectively.
Unfortunately, the projectors on the divergencefree wavelet and curlfree 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 divergencefree and curlfree parts , we have implemented the following iterative algorithm:
Set u^{0} = v^{0} = 0 and the maximum number of iterations to I_{max}.
Set i = 0, while γ − u^{i} − v^{i}^{2} > ϵ and i < I_{max} then iterate:
Set
Set
Set
Set
Set i = i + 1.
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
Set u^{0} = v^{0} = 0 and the maximum number of iterations I_{max}.
Set i = 0, while γ − u^{i} − v^{i}^{2} > ϵ and i < I_{max} then iterate:
Set
Set
Set
Set
Set
Set
Set
Set
Set i = i + 1.
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

Set u^{0} = v^{0} = 0 and the maximum number of iterations to I_{max}.

Set i = 0, while γ − u^{i} − v^{i}^{2} > ϵ and i < I_{max} then iterate:
Set
Set with:
Compute
Set
Set
Set
Set with
Compute
Set
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 08EMER00901) and the European Research Council grant SparseAstro (ERC228261).
References
 Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Bunn, E., Zaldarriaga, M., Tegmark, M., & de OliveiraCosta, A. 2003, Phys. Rev. D, 67, 023501 [NASA ADS] [CrossRef] [Google Scholar]
 Cohen, A., Daubechies, I., & Vial, P. 1993, Appl. Comput. Harmon. Anal., 1, 54 [CrossRef] [MathSciNet] [Google Scholar]
 Crittenden, R. G., Natarajan, P., Pen, U.L., & Theuns, T. 2002, ApJ, 568, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Daubechies, I., & Sweldens, W. 1998, J. Fourier Anal. Appl., 4, 247 [Google Scholar]
 Dautray, R., & Lions, J.L. 1984, Analyse mathématique et calcul numérique pour les sciences et les techniques, Vol. 5 (Masson) [Google Scholar]
 Deriaz, E., & Perrier, V. 2006, J. Turbulence, 7, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Deriaz, E., & Perrier, V. 2009, Appl. Comput. Harmon. Anal., 26, 249 [CrossRef] [Google Scholar]
 DeVore, R. A., Konyagin, S. V., & Temlyakov, V. N. 1998, Const. Appr., 14, 1 [CrossRef] [MathSciNet] [Google Scholar]
 Fu, L., & Kilbinger, M. 2010, MNRAS, 401, 1264 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, N., & Squires, G. 1993, ApJ, 404, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Kilbinger, M., Schneider, P., & Eifler, T. 2006, A&A, 457, 15 [Google Scholar]
 Kratochvil, J. M., Lim, E. A., Wang, S., et al. 2011 [arXiv:1109.6334] [Google Scholar]
 LemariéRieusset, P. 1992, Revista Matemática Iberoamericana, 8, 221 [CrossRef] [Google Scholar]
 Lombardi, M., & Bertin, G. 1998, A&A, 335, 1 [NASA ADS] [Google Scholar]
 Marian, L., Smith, R. E., Hilbert, S., & Schneider, P. 2011, MNRAS, submitted [arXiv:1110.4635] [Google Scholar]
 Massey, R., Rhodes, J., Ellis, R., et al. 2007, Nature, 445, 286 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Mellier, Y. 1999, ARA&A, 37, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Mellier, Y. 2002, Space Sci. Rev., 100, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Meyer, Y. 1990, Ondelettes et opérateurs, I. Ondelettes, II. Opérateurs de CalderónZygmund, III. Opérateurs multilinéaires (Hermann) [Google Scholar]
 Munshi, D., Valageas, P., van Waerbeke, L., & Heavens, A. 2008, Phys. Rep., 462, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Munshi, D., Heavens, A., & Coles, P. 2011, MNRAS, 411, 2161 [NASA ADS] [CrossRef] [Google Scholar]
 Pichon, C., Thiébaut, E., Prunet, S., et al. 2010, MNRAS, 401, 705 [Google Scholar]
 Pires, S., Starck, J., Amara, A., Réfrégier, A., & Teyssier, R. 2009a, A&A, 505, 969 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pires, S., Starck, J., Amara, A., et al. 2009b, MNRAS, 395, 1265 [NASA ADS] [CrossRef] [Google Scholar]
 Pires, S., Starck, J., & Refregier, A. 2010, IEEE Sig. Proc. Mag., 27, 76 [Google Scholar]
 Refregier, A. 2003, ARA&A, 41, 645 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, P., & Kilbinger, M. 2007, A&A, 462, 841 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., van Waerbeke, L., Kilbinger, M., & Mellier, Y. 2002, A&A, 396, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schrabback, T., Hartlap, J., Joachimi, B., et al. 2010, A&A, 516, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seitz, S., & Schneider, P. 1996, A&A, 305, 383 [NASA ADS] [Google Scholar]
 Seitz, S., & Schneider, P. 2001, A&A, 374, 740 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Starck, J., Pires, S., & Réfrégier, A. 2006, A&A, 451, 1139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stevenson, R. 2010, Appl. Comput. Harmon. Anal., in press [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van Waerbeke, L., Mellier, Y., Radovich, M., et al. 2001, A&A, 374, 757 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Figures
Fig. 1 Galaxy shapes and corresponding gravitational shears γ for a circular galaxy. 

In the text 
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. 

In the text 
Fig. 3 Example of a 2D divergencefree wavelet basis designed from the Lemarié & Rieusset (1992) scheme. 

In the text 
Fig. 4 Examples of hyperbolic 2D divergencefree and curlfree wavelets. 

In the text 
Fig. 5 Example of divergencefree wavelets with a free (left) and a sliding boundary condition (right). 

In the text 
Fig. 6 Left, simulated κ^{S} image and its corresponding shear field γ^{S}. Right, extracted subfield γ^{D}. The extraction breaks the periodicity of the data. 

In the text 
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). 

In the text 
Fig. 8 Bmode reconstructed maps. 

In the text 
Fig. 9 Maps of the reconstruction error on κ_{E} for three different methods when B modes contaminate the data. 

In the text 
Fig. 10 Original B mode and reconstructed error maps for the B modes map κ_{B}. 

In the text 
Fig. A.1 Quadratic spline scaling functions adapted to the left boundary of the interval. 

In the text 
Fig. A.2 Linear spline scaling functions adapted to the interval. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.