Weak gravitational lensing

Gravitational lensing can lead to the formation of multiple images of the same background source, especially if it is in near perfect alignment towards the line-of-sight of a massive lens. Owing to this condition of alignment, such phenomenon are rare (about 1 in 1000 massive elliptical galaxies act as lenses). This is the strong lensing regime. In most cases however, the background galaxy will not be perfectly aligned to produce multiple images, but nevertheless its shape will get distorted due to the shear due to object at its location.

This distortion can be characterized using the Jacobian of the mapping from the source plane to the lens plane at the observed location of the image.

\begin{eqnarray} \mathcal{M}^{-1}(\boldsymbol{\theta}) &=& \begin{pmatrix} 1 - \kappa - \gamma_1 & -\gamma_{2}\\ -\gamma_2 & 1-\kappa + \gamma_1 \end{pmatrix} \\ &=& (1-\kappa)\begin{pmatrix} 1-\frac{\gamma_1}{(1-\kappa)} & -\frac{\gamma_2}{(1-\kappa)} \\ -\frac{\gamma_2}{(1-\kappa)} & 1+\frac{\gamma_1}{(1-\kappa)} \end{pmatrix} \end{eqnarray}

where the two components of the shear are given by \begin{eqnarray} \gamma_1 &=& \frac{1}{2} \left[ \psi_{,11} - \psi_{,22} \right] \\ \gamma_2 &=& \psi_{,12} \end{eqnarray} and that the larger eigenvector of the distortion is inclined at an angle \(\phi\), with \begin{eqnarray} \tan 2\phi = \frac{\gamma_2}{\gamma_1}\,. \end{eqnarray}

A circular source thus gets distorted in to an ellipse, where the angle of orientation is related the relative ratio of the two shear components, and the ellipticity is related to the reduced shear at the location of the source. \begin{eqnarray} e=\frac{a-b}{a+b} &=& \frac{\frac{1}{1-\kappa-|\gamma|} - \frac{1}{1-\kappa+|\gamma|}}{\frac{1}{1-\kappa-|\gamma|} + \frac{1}{1-\kappa+|\gamma|}} \\ &=& \frac{|\gamma|}{1-\kappa} \end{eqnarray} One can also consider another way to quantify the ellipticity \begin{eqnarray} \chi=\frac{a^2-b^2}{a^2+b^2} &=& \frac{\frac{1}{(1-\kappa-|\gamma|)^2} - \frac{1}{(1-\kappa+|\gamma|)^2}}{\frac{1}{(1-\kappa-|\gamma|)^2} + \frac{1}{(1-\kappa+|\gamma|)^2}} \\ \end{eqnarray} This shows \begin{eqnarray} \chi =\frac{2|\gamma|}{(1-\kappa)} \end{eqnarray}

The above figure shows how the ellipses are oriented given the values of these components. Pay special attention to cases along the two different axes.

Weak lensing as mass probe

Let us now consider the case where we have determined the estimate of shear around a galaxy cluster. How do we infer the surface mass density of the cluster from the shear?

Consider the following figure which shows how the presence of shear distorts the shapes of galaxies. The left hand panel shows a distribution of background galaxies which are intrinsically circular in shape. The right hand panel shows the impact of including an azimuthally symmetric mass distribution. We see that the background shapes are now distorted in a coherent manner, all tangentially oriented with respect to the center of the mass distribution. The galaxies nearer to the center have larger distortions while those further away have weaker distortions.

We can define the tangential shear by defining a coordinate system with one axis which is the line joining the center of the mass distribution to the position of the source galaxy. Let \(\theta\) be the angle made by this line with respect to the x-axis. \begin{eqnarray} \gamma_{\rm t} = - {\rm Re}({\boldsymbol \gamma }e^{-2i\theta}) \end{eqnarray}

Tangential shear average over circles

Consider first Gauss’s law and integrate the Laplacian of the lensing potential within a disk inside radius R.

\begin{eqnarray} \int_{0}^{\theta} d\vec\theta' \nabla. \nabla \psi &=& \theta \oint d\varphi \nabla \psi.\hat n \\ 2 \int_0^{\theta} d\vec\theta' \kappa(\vec\theta') &=& \theta \oint d\varphi \partial_{\theta} \psi \\ 2 \int_0^{\theta} d\varphi \vartheta d\vartheta \kappa(\vartheta, \varphi) &=& \theta \oint d\varphi \partial_{\theta} \psi \\ 2 \int_0^{\theta} 2 \pi \vartheta d\vartheta \langle\kappa(\vartheta)\rangle &=& \theta \oint d\varphi \partial_{\theta} \psi\,, \end{eqnarray} where we define \begin{eqnarray} \langle\kappa(\vartheta)\rangle &=& \frac{1}{2\pi} \int d\varphi \kappa(\vartheta, \varphi) \end{eqnarray} which denotes the scaled surface density \(\kappa\) averaged over all angles for a given \(\theta\).

Differentiating both sides with respect to \(\theta\), we obtain \begin{eqnarray} 4\pi \theta \langle\kappa(\theta)\rangle &=& \frac{\partial}{\partial \theta} \left[\theta \oint d\varphi \partial_\theta\psi\right] \\ &=& \frac{2}{\theta} \int_0^{\theta} d\vec\theta' \kappa(\vec\theta') + \theta \oint d\varphi \partial^2_{\theta} \psi\\ &=& \frac{2}{\theta} \int_0^{\theta} d\vec\theta' \kappa(\vec\theta') + \theta \oint d\varphi \left[\kappa(\theta) - \gamma_{\rm t}(\theta)\right] \\ &=& \frac{2 \pi \theta}{\pi \theta^2} \int_0^{\theta} d\vec\theta' \kappa(\vec\theta') + 2\pi \theta \left[\langle\kappa(\theta)\rangle - \langle\gamma_{\rm t}(\theta)\rangle\right] \\ &=& 2 \pi \theta \bar\kappa(<\theta)+ 2\pi \theta \left[\langle\kappa(\theta)\rangle - \langle\gamma_{\rm t}(\theta)\rangle\right]\\ \langle\gamma_{\rm t}\rangle &=& \bar\kappa(<\theta) - \langle\kappa(\theta)\rangle \end{eqnarray}

The third equality can be understood in the following way: consider a point on the x-axis \(\partial^2\psi/\partial\theta^2=\psi_{11}=\kappa+\gamma_{\rm 1}\). However in this case, \(\gamma_1=-\gamma_{\rm t}\). Thus for the point on the x-axis \(\partial^2\psi/\partial\theta^2=\kappa-\gamma_{\rm t}\). The double derivative can be evaluated at every point by performing a rotation, and it will always yield the value \(\kappa-\gamma_{\rm t}\).

The last equality is a very important result. It says that regardless of the absence of symmetry, when the tangential shear is averaged over an annulus surrounding any mass distribution it’s value is the difference between the average surface density inside that annulus and the surface density averaged over the annulus. Thus measurements of the averaged shear on annuli can be used to derive the surface density profile of any object.

Mass maps from shear estimates

Let us consider the definition for the lensing potential once again:

\begin{eqnarray} \psi(\vec\theta)&=& \left[ \frac{1}{\pi}\int d^2\vec \theta' \kappa(\vec \theta') \ln (|\vec\theta-\vec\theta'|)\right]\\ \end{eqnarray} and the definition of the shear components: \begin{eqnarray} \gamma_1 &=& \frac{1}{2} \left[ \psi_{,11} - \psi_{,22} \right] \\ \gamma_2 &=& \psi_{,12} \end{eqnarray}

One can consider the Fourier transform of the potential \begin{eqnarray} \tilde \psi({\vec\ell}) = \int d\vec\theta e^{-i\vec\ell.\vec\theta}\,\psi \end{eqnarray} The derivatives are simpler in Fourier space. Given that the lensing potential and the surface density obey Poisson equation \(\nabla^2\psi = 2 \kappa\), we have \begin{eqnarray} |\vec\ell|^2\tilde \psi({\vec\ell}) = -2\kappa(\vec\ell) \end{eqnarray} From the definition of double derivative we also obtain \begin{eqnarray} {\boldsymbol\gamma} &=& -\left(\frac{1}{2} \left[ \ell_{1}^2 - \ell_{2}^2 \right] + i \ell_1 \ell_2\right) \psi (\vec\ell) \\ &=& \frac{1}{|\vec\ell|^2}\left( \left[ \ell_{1}^2 - \ell_{2}^2 \right] + i 2 \ell_1 \ell_2\right) \kappa(\vec\ell)\\ \kappa(\vec\ell) &=& \frac{1}{|\vec\ell|^2} \left( \left[ \ell_{1}^2 - \ell_{2}^2 \right] - i 2 \ell_1 \ell_2\right) {\boldsymbol\gamma}(\vec\ell) \end{eqnarray} where in the last equality we have used \({\bf a}^{-1} = |a|^{-2}{\bf a}^{*}\).

This is another interesting and useful result. If we can measure the shear field in a region on a grid, take its Fourier transform using FFT and use the above equation, one can get the convergence field (surface density field). The corresponding result in real space is (Kaiser & Squires, 1993) \begin{eqnarray} \kappa = \frac{1}{\pi}\int d\vec\theta' {\cal D}_{\gamma\rightarrow\kappa}(\vec\theta-\vec\theta').{\boldsymbol\gamma}(\vec\theta') \end{eqnarray} where \begin{eqnarray} {\cal D}_{\gamma\rightarrow\kappa}(\vec\theta) = \frac{1}{|\vec\theta|^4}[\theta_2^2-\theta_1^2, 2\theta_1\theta_2] \end{eqnarray} Because \([\theta_1^2-\theta_2^2)^2 + (2\theta_1\theta_2)^2]=|\vec\theta|^2\), we can write this as \begin{eqnarray} {\cal D}_{\gamma\rightarrow\kappa}(\vec\theta) = \frac{1}{|\vec\theta|^2}[\cos 2 \varphi, \sin 2\varphi] \end{eqnarray} This kernel is similar to the ellipticity pattern induced by a point mass around its position. Thus there is an intuitive way to understand this equation in real space. It is a sort of a match filter, the value of \(\kappa\) will be large when the pattern matches the pattern induced by a point mass.

Another way to see some of the similarity between the expressions is as follows. Let us start from the equation for the deflection angle in terms of components \begin{eqnarray} \psi_{,i} = \alpha_i = \frac{1}{\pi}\int d\vec{\theta'} \kappa(\vec\theta') \frac{\theta_i-\theta_i'}{|\vec\theta-\vec\theta'|^2} \end{eqnarray} Differentiating once again with respect to \(\theta_j\), \begin{eqnarray} \psi_{,ij} = \frac{1}{\pi}\int d\vec{\theta'} \kappa(\vec\theta') \left[\frac{\delta^{\rm D}_{ij}}{|\vec\theta-\vec\theta'|^2} - 2\frac{(\theta_i-\theta'_i)(\theta_j-\theta'_j)}{|\vec\theta-\vec\theta'|^4}\right] \end{eqnarray} Using the definitions of the shear components we get \begin{eqnarray} {\boldsymbol \gamma} = \frac{1}{\pi}\int d\vec\theta' {\cal D}_{\kappa\rightarrow\gamma}(\vec\theta-\vec\theta') \kappa(\vec\theta') \end{eqnarray} where \begin{eqnarray} {\cal D}_{\kappa\rightarrow\gamma}(\vec\theta) = \frac{\theta_2^2 - \theta_1^2-2i\theta_1\theta_2}{|\vec \theta|^4} \end{eqnarray}

Shape measurements

Clearly once the shear of galaxies is measured, it allows us to create mass maps in the Universe and to study the matter distribution around galaxies. Let us next discuss how one can measure the shapes of galaxies and how one can obtain the value of the shear from them.

Let \(I(\vec\theta)\) be the surface brightness of an object. One can define the centroid of the object by computing the first moment of the object \begin{eqnarray} {\vec\theta}_{\rm cen} = \frac{\int d\vec\theta\,I(\vec \theta)\,q(I)\,\vec\theta}{\int d\vec\theta\,I(\vec\theta)\,q(I) } \end{eqnarray} where we have included a weight \(q(I)\). The second moments of the brightness distribution are given by \begin{eqnarray} Q_{ij} = \frac{\int d\vec\theta\,I(\vec \theta)\,q(I)\,(\theta_i-\theta_{\rm cen, i})(\theta_j-\theta_{\rm cen, j})}{\int d\vec\theta\,I(\vec\theta)\,q(I) } \end{eqnarray}

One can define two different types of ellipticities with this definition \begin{eqnarray} \chi = \frac{Q_{11}-Q_{22}+2iQ_{12}}{Q_{11}+Q_{22}}\,; \epsilon = \frac{Q_{11}-Q_{22}+2iQ_{12}}{Q_{11}+Q_{22}+[Q_{11}Q_{22}-Q_{12}^2]^{1/2}}\, \end{eqnarray} Consider an elliptical galaxy with a gaussian shape oriented in the positive x-axis, the two axes have sigmas equal to \(a\) and \(b\). In this case \(Q_{11}=a^2; Q_{22}=b^2; Q_{12}=0\). The first definition will result in \begin{eqnarray} \chi &=& \frac{a^2-b^2}{a^2+b^2}\,;\\ &=& \frac{1-|r|^2}{1+|r|^2}\\ \epsilon &=& \frac{a^2-b^2}{a^2+b^2+2ab}\,\\ &=& \frac{1-|r|}{1+|r|} \end{eqnarray} where r is the minor to major axis ratio.

Similar to the image ellipticities one can also infer the source ellipticities and figure out how the lensing transformation changes the ellipticities of galaxies. Locally the transformation is given by the inverse of magnification matrix. The moments of the source will be given by \begin{eqnarray} Q_{ij}^{\rm s} = \frac{\int d\vec\beta I^{\rm s}(\vec\beta) q[I^s] (\beta_{\rm cen,i}-\beta_i)(\beta_{\rm cen, j}-\beta_j)}{\int d\vec\beta I^{\rm s}(\vec\beta) q[I^s]}\,, \end{eqnarray} Using \(\beta_{\rm cen, i}-\beta_i=A_{ik} (\theta_{\rm cen, i} - \theta)\), and that the surface brightness is conserved, we obtain \begin{eqnarray} Q_{ij}^{\rm s} &=& \frac{\int d\vec\theta I(\vec\theta) q[I] A_{ik}(\theta_{\rm cen,i}-\beta_i)A_{jl}(\theta_{\rm cen, j}-\theta_j)}{\int d\vec\theta I(\vec\theta) q[I]}\,,\\ &=& A_{ik}A_{jl} Q_{kl} \end{eqnarray} Thus the lens local transformation can be used to relate the image ellipticities to the source ellipticities. The rule of transformation is given by \begin{eqnarray} \chi^{\rm s} &=& \frac{\chi-2g+g^2\chi^*}{1+|g|^2-2{\rm Re}(g\chi^*)}\\ \epsilon^{\rm s} &=& \frac{\epsilon-g}{1-g^*\epsilon} {\,\,\,\rm if\,\,|g|\le 1}\\ \epsilon^{\rm s} &=& \frac{1-g\epsilon^*}{\epsilon^*-g^*} {\,\,\,\rm if\,\,|g|\gt 1}\\ \end{eqnarray} If the source galaxies are randomly oriented with respect to each other the ensemble average \(\langle\chi^{\rm s}\rangle=\langle\epsilon^{\rm s}\rangle = 0\), which allows us to determine the value of g by performing an ensemble average of the measured ellipticities. \begin{eqnarray} \langle \epsilon \rangle &=& g {\,\,\,\rm if\,\,|g|\le 1}\\ &=& 1/g^* {\,\,\,\rm if\,\,|g|\gt 1} \end{eqnarray} The case related to the choice \(|g|\gt 1\) is related to negative parity images (elliptical sources do not care about parity), which can occur in the case of strong lensing. These form in a very small region near the center of the mass distribution where such effects may be more important. But however local measurements of shapes cannot inform us about which regime we are at. One can define distortion which is invariant under change between \(g\) and \(1/g^{*}\), \begin{eqnarray} \delta = \frac{2g}{1+|g|^2} \end{eqnarray} The averaging of \(\chi^{\rm s}\) is a little more involved, and Schneider and Seitz (1995) show that the average of \(\chi^{\rm s}\) equal to \(0\) implies \begin{eqnarray} \sum_{i} u_i\frac{\chi_i-\delta}{1+{\rm Re}(\delta \chi^{*}_i)} = 0 \end{eqnarray} where \(u_i\)’s represent the different weight functions for galaxies depending upon how far they are from the point that we want to estimate the shear. These weights could also be adjusted to give weights to the source galaxies depending upon their SNR, for example. This equation can be solved iteratively to obtain the value of \(\delta\) given the measurements of \(\chi_i\).

In the weak lensing regime, it can be shown that \begin{eqnarray} \gamma \approx g \approx \langle \epsilon \rangle \approx \frac{1}{2} \langle \chi \rangle \approx \frac{1}{2} \langle \delta \rangle \end{eqnarray}

These relations imply that even if we cannot use a single galaxy to perform the weak lensing measurement, when we average over an ensemble of objects with similar shear (say in a small region of space) we will be able to obtain an estimate of the shear. Once the shear estimate can be obtained we can use either the averaging on circles or the mass inversion technique mentioned above in order to obtain mass distributions in the Universe.

Practical issues

We will look at the practical steps taken in order to reduce images from real telescopes in the next chapter. Here we discuss some of the issues relevant for the shape measurements following techniques in Kaiser, Squires and Broadhurst 1995.

Weights in moment measurements

The weight function used to define the moments above that we chose were a function of the surface brightness. However, in real life this is much more complicated. Measurements of local surface brightness are difficult specially close to the noise background. The noise can enter the measured ellipticities in a non-linear fashion in this case. To avoid these issues typically the weight is chosen as a function of \(\vec\theta\) instead (say a circular gaussian), which avoids the noise in the image from entering the ellipticity measurement. However in this case the expression for how to go from the ellipticity measurements to shear no longer remain perfectly correct for such a weighting scheme.

Point spread function

The second important issue is that of the point spread function, primarily due to atmosphere for ground based data and that of the telescope optics for space based data. The image observed on the detector is convolution of the true image with the PSF. \begin{eqnarray} I^{\rm obs}(\vec\theta) = \int d\vec\vartheta I(\vec\vartheta) P(\vec\theta-\vec\vartheta) \end{eqnarray} where \(P\) represents the PSF and is normalized to unity and is centered at zero. A delta-function image gets converted into the same shape as PSF. In general the PSF smears the image of an object. An isotropic PSF will reduce the ellipticity by smoothing out the image, especially for small objects whose size is comparable to the PSF. In addition, any anisotropy of the PSF will leak into the measured ellipticity.

Ideal weighting

If the weights used for the determination of ellipticity of galaxies are \(q(I)=I\), in that case there is a straightforward relation between the PSF ellipticity, the image ellipticity and the observed ellipticity on the detector after the PSF convolution,

\begin{eqnarray} \chi^{\rm obs} = \frac{\chi+T\chi^{\rm PSF}}{1+T} \end{eqnarray} where \begin{eqnarray} T&=&\frac{P_{11}+P_{22}}{Q_{11}+Q_{22}}\\ \chi^{\rm PSF} &=& \frac{P_{11}-P_{22}+2iP_{12}}{P_{11}+P_{22}} \end{eqnarray} T represents the ratio of the PSF size to the galaxy size before convolution. If T is small the galaxy is well resolved (the PSF is much smaller than the galaxy), then the correction factor is small. However for getting more and more galaxies to perform the average, we need to use smaller galaxies and so we run into galaxies with sizes comparable to the PSF size.

In reality however we use weight functions which are not the same as the ideal one above, and so the relation between how the PSF changes the ellipticity needs to be calibrated using realistic image simulations.

Dealing with anisotropic PSF

The total PSF can be decomposed into an isotropic and an anisotropic part:

\begin{eqnarray} P(\vartheta) = \int d\vec\varphi q(\vec\varphi) P^{\rm iso}(\vec\vartheta-\vec\varphi) \end{eqnarray}

One defines \begin{eqnarray} I^{\rm iso}(\vec\theta) &=& \int d\vec\varphi I(\vec\varphi) P^{\rm iso}(\vec\theta-\vec\varphi)\\ I^{\rm iso, s}(\vec\theta) &=& \int d\vec\varphi I^{\rm s}(\vec\varphi) P^{\rm iso}(\vec\theta-\vec\varphi)\,\\ \end{eqnarray} and can use these definitions to compute moments (using the position dependent weights) and ellipticities \(\chi^{\rm iso}\) and \(\chi^{\rm iso, s}\). The isotropically smeared ellipticities of the unlensed source still average out to zero, and thus such an average can be used to figure out the shear at a particular location. Both these quantities cannot be observed but we can connect them to the observed ellipticity \(\chi^{\rm obs}\) via a two step process: \(\chi^{\rm obs}\rightarrow \chi^{\rm iso}\) and \(\chi^{\rm iso}\rightarrow \chi^{\rm iso, s}\).

Since the anisotropies and shears tend to be small one can compute the effect up to linear order as \begin{eqnarray} \chi^{\rm obs}_{i} = \chi^{\rm iso, s}_{i} + P^{\rm sm}_{ij}q_{j} + P^{\rm g}_{ij} g_{j} \end{eqnarray} The quantity \(P^{\rm sm}\) is called the smear polarizability and is the linear response of the observed ellipticity to the anisotropy of the PSF. The quantity P^{\rm g} is called the shear polarizability and is the linear response of the observed ellipticity to the shear. The reduced shear can then be obtained as \begin{eqnarray} g = (P^{\rm g})^{-1}(\chi^{\rm obs} - P^{\rm sm}q) \end{eqnarray} Although these quantities can be computed for each image given the PSF size and the brightness distribution of the observed image.

Adaptive moments for shapes

There are new advances in shape measurements which reduce issues with the KSB techniques. In particular using weight functions which are adaptive to the size of the galaxy. The ellipticity of a galaxy can be measured by asking for the transform that when applied in inverse results in a round galaxy. In this case the optimal weights turn out to be (see Bernstein and Jarvis 2002) \begin{eqnarray} w(\theta) = \frac{1}{\theta}\frac{d\ln I}{d\theta} \end{eqnarray} with the shape of the weight adapted to the shape of the galaxy under consideration. In principle one can change the weight to reflect the profile of each galaxy, but it turns out to be easier and more practical to use a Gaussian weight as it downweights the regions away from the center of the object (which can be a problem in crowded fields). Such shapes can be found by minimizing \begin{eqnarray} \chi^2 = \int d\vec\theta \left[ I(\vec\theta) - A \exp\left(-\frac{1}{2}(\vec\theta-\vec\theta_0)^{\rm T}M^{-1}(\vec\theta-\vec\theta_0)\right) \right]^2 \end{eqnarray} This allows one to fit for the shape moments of a galaxy \(M\) (3 independent numbers), in addition to the amplitude A and the position coordinates \(\vec\theta_0\) (2 more numbers). One can show that the shapes which minimize the above \(\chi^2\) are weighted moments with \begin{eqnarray} \theta_0 &=& \frac{\int d\vec\theta I(\vec\theta) w(\vec\theta) \theta}{\int d\vec\theta I(\vec\theta) w(\vec\theta)} \\ M_{ij} &=& 2\frac{\int d\vec\theta I(\vec\theta) w(\vec\theta) (\theta_i-\theta_{0, i})(\theta_j-\theta_{0, j})}{\int d\vec\theta I(\vec\theta) w(\vec\theta)} \end{eqnarray}

Bernstein and Jarvis 2002 also present an estimator for the distortion \(\delta\) for shapes measured using such adaptive moments when considering an ensemble of sources. They consider how the distribution of the ellipticities transform under the application of a distortion to a random distribution and how to use this change to estimate the distortion. They compute the optimal weights in the case of general distribution \(P(e)\), but for a Gaussian distribution for \(P(e)\) and in the case of small measurement noise, these weights are approximately given by \begin{eqnarray} w \approx \frac{1}{{e_{\rm rms}^2+\sigma_{\rm meas}^2}} \end{eqnarray} where \(e_{\rm rms}\) is the rms intrinsic distortion per component. Thus one can use the estimator \(\hat\delta\) \begin{eqnarray} \hat\delta = \frac{\sum w_i e_i}{\sum w_i} \end{eqnarray} The estimator needs to be corrected for the responsivity \begin{eqnarray} {\cal R} &=& \frac{\partial \hat\delta}{\partial \delta}\\ &\approx& 1 - e_{\rm rms}^2 \end{eqnarray} Thus the shear estimate is given by \begin{eqnarray} 2\gamma &=& \delta &=& \frac{1}{{\cal R}}\frac{\sum w_i e_i}{\sum w_i} \end{eqnarray}

PSF corrections

PSF convolution is straightforward when the weighting used to determine the moments are ideal, however as mentioned before such methods have to deal with problems of noise which enter the ellipticity in a non-linear manner. There are several ways proposed to deal with PSF corrections when weighted moments are used for ellipticity determination.

One can potentially measure the dilution of PSF using detailed image simulation, however the correction factors can be dependent upon how realistic the simulations are in terms of their unlensed galaxy properties. Perturbative methods like KSB described above have their own limitations when applied to realistic PSFs (Bernstein & Jarvis 2002). Another approach is to deconvolve the PSF by using stacks of galaxies in similar locations (Kuijken 1999). However this requires the galaxies to be similarly sheared and a large number of background galaxies in a region where the PSF does not vary substantially. Bernstein & Jarvis (2002) suggest the use Laguerre exponential functions as basis functions to decompose the image and PSFs into. Such expansions have convenient properties in order apply shears or for convolutions and thus can be used to correct for the effects of the PSF. However, these do require some truncation in the basis fucntion expansion. For this Bernstein & Jarvis (2002) suggest applying a rounding kernel to the PSF first and give an expression to correct for the PSF designed to be accurate in the case of Gaussian PSF and Gaussian objects and which works for well resolved galaxies.

Hirata and Seljak (2003) expand on this and suggest that non-Gaussian PSFs can be dealt with in the following manner:

Consider a PSF \(g(\vec{x})\) and approximate it with a Gaussian \(G(\vec x)\), and consider the residual function \(\epsilon(\vec x) = g-G\). Let \(I^{\rm pre}(x)\) be the original pre-seeing image, then the image \(I^{\rm post}(x) = g\otimes I^{\rm pre}(x) = G \otimes I^{\rm pre}(x) + \epsilon\otimes I^{\rm pre}(x)\). Thus one can determine \(G \otimes I^{\rm pre}(x) = I^{\rm post}(x) - \epsilon\otimes I^{\rm pre}(x)\). We do not really know what is \(I^{\rm pre}(x)\) but since this is convolved with a smaller function \(\epsilon\) we can use an approximate function \(\tilde I^{\rm pre}\) which is a deconvolved image assuming the PSF and the original image to be Gaussians. The Bernstein and Jarvis 2002 can now be made even more accurate. This technique called re-Gaussianization is used to then perform the PSF correction.

Tangential shear computation

The shape measurements for galaxies are carried out typically in sky coordinates. Suppose we want to use method one of averaging the tangential shear around an object (say galaxy G1) at coordinate with right ascension and declination \((\alpha_1, \delta_1)\). Let us suppose we have measured the shape of another galaxy (G2) at \((\alpha_2, \delta_2)\). We want to define the tangential shear with respect to the line joining G1 and the galaxy G2 whose shape we have measured.

Consider this Figure from Kiblinger et al. 2013.

Angles should always be measured with respect to great circles. We need to consider the angle between the great circle joining G1G2 and the great circle passing through G2, in order to define the tangential shear. The great circle passing through G1G2 has a normal given by \(\hat n_{G_1} \times \hat n_{G_2}/\sin(\vartheta)\), where \(\hat n\) are the unit vectors which point towards the corresponding galaxies. The unit vectors are given by \(\cos\delta \cos\alpha \hat i + \cos\delta \sin\alpha \hat j + \sin\delta \hat k\). The cosine of the angle between the two positions is given by

\begin{eqnarray} \cos\vartheta &=& \cos\delta_1\cos\alpha_1\cos\delta_2\cos\alpha_2 + \cos\delta_1 \sin\alpha_1 \cos\delta_2 \sin\alpha_2 + \sin\delta_1\sin\delta_2\\ &=& \cos\delta_1\cos\delta_2 \cos(\alpha_1-\alpha_2) + \sin\delta_1\sin\delta_2 \end{eqnarray}

The great circle passing through G2 and the north pole can be calculated using a similar dot product between the unit vector denoting the galaxy and the unit vector passing through the north pole. This product will give a vector of magnitude \(\sin(90-\delta_2)=\cos\delta_2\).

import sympy as s
alpha1, delta1, alpha2, delta2,theta = s.symbols("alpha1,delta1,alpha2,delta2 theta")

# Define unit vectors
nG1 = s.Matrix([s.cos(delta1)*s.cos(alpha1), s.sin(alpha1)*s.cos(delta1), s.sin(delta1)])
nG2 = s.Matrix([s.cos(delta2)*s.cos(alpha2), s.sin(alpha2)*s.cos(delta2), s.sin(delta2)])

# Define north pole
nNP = s.Matrix([0, 0, 1])

# G1 \cross G2/sin(theta)
L_G1G2=nG1.cross(nG2)/s.Abs(s.sin(theta))

# G2 \cross NP/cos(delta2)
L_GC = nG2.cross(nNP)/s.cos(delta2)

# Compute dot product for cos(beta2) = sin(phi2)
s.simplify(((L_GC.dot(L_G1G2))))

$\displaystyle \frac{- \sin{\left(\delta_{1} \right)} \cos{\left(\delta_{2} \right)} + \sin{\left(\delta_{2} \right)} \cos{\left(\delta_{1} \right)} \cos{\left(\alpha_{1} - \alpha_{2} \right)}}{\left|{\sin{\left(\theta \right)}}\right|}$

# Compute dot product G1 G2 to get \cos\theta
s.simplify(nG1.dot(nG2))

$\displaystyle \sin{\left(\delta_{1} \right)} \sin{\left(\delta_{2} \right)} + \cos{\left(\delta_{1} \right)} \cos{\left(\delta_{2} \right)} \cos{\left(\alpha_{1} - \alpha_{2} \right)}$

# Obtain sin(beta2) = cos(phi2)
x = L_G1G2.cross(L_GC)
s.simplify(x)

$\displaystyle \left[\begin{matrix}- \frac{\sin{\left(\alpha_{1} - \alpha_{2} \right)} \cos{\left(\alpha_{2} \right)} \cos{\left(\delta_{1} \right)} \cos{\left(\delta_{2} \right)}}{\sin{\left(\theta \right)}}\\- \frac{\sin{\left(\alpha_{2} \right)} \sin{\left(\alpha_{1} - \alpha_{2} \right)} \cos{\left(\delta_{1} \right)} \cos{\left(\delta_{2} \right)}}{\sin{\left(\theta \right)}}\\- \frac{\sin{\left(\delta_{2} \right)} \sin{\left(\alpha_{1} - \alpha_{2} \right)} \cos{\left(\delta_{1} \right)}}{\sin{\left(\theta \right)}}\end{matrix}\right]$

The above calculations show that \begin{eqnarray} \cos\beta_2 &=& \sin\phi_2 = \frac{-\sin\delta_1\cos\delta_2+\sin\delta_2\cos\delta_1\cos(\alpha_2-\alpha_1)}{|\sin\theta|}\\ \sin\beta_2 &=& \cos\phi_2 = \frac{ \sin(\alpha_2-\alpha_1)\cos\delta_1}{|\sin\theta|} \end{eqnarray}

One has to thus compute the tangential shear with respect to the line at an angle \(\phi_2\) given the shape of the galaxy. Therefore,

\begin{eqnarray} e_{\rm tan} &=& - {\rm Re}\left[(e_1 + i e_2)e^{-i2\phi_2}\right]\\ &=& -e_1\cos 2\phi_2 - e_2 \sin2\phi_2 \end{eqnarray}

These values of the tangential ellipticity can be averaged (in a weighted manner) over a large number of sources to yield measurements of the tangential shear profile.

\begin{eqnarray} \langle e_{\rm tan}\rangle = \frac{\sum w_i e_{\rm tan, i}}{w_i} \end{eqnarray} which can then be related to the tangential shear depending upon the definition of ellipticity used.

Measurement of the weak lensing signal

The weak lensing signal can be measured from an ensemble average of ellipticities. \begin{eqnarray} \langle \gamma_{\rm t} (\theta) \rangle &=& \bar\kappa(<\theta) - \langle\kappa\rangle(\theta)\\ &=& \frac{\bar\Sigma(<\theta)-\langle\Sigma\rangle(\theta)}{\Sigma_{\rm crit}} \\ &=& \frac{\Delta \Sigma}{\Sigma_{\rm crit}}\\ &=& \frac{1}{2{\cal R}}\frac{\sum w_i e_{\rm tan}}{\sum w_i} \end{eqnarray} In the last line we have assumed that the ellipticity measure used gives an estimate of distortion. The fact that the background source galaxies are not circular results in a statistical noise in the measurement of the average tangential shear from the shape of galaxies. This noise can be beaten down by averaging over a large number of source galaxies. A single object may not have enough number of source galaxies to beat down the shape noise. In that case, we can average over a large number of lensing objects. The source galaxies that we consider may also not lie at the same redshift, and so each source galaxy will have a different critical surface density. In this case, it is better to devise an estimator for \(\Delta\Sigma\) instead of \(\gamma_{\rm t}\). An individual galaxy gives an estimate of \(\gamma_{\rm t}\Sigma_{\rm crit}/(2{\cal R})\). The error on this individual estimate is \([e_{\rm rms}^2+e_{\rm meas}^2]^{1/2}\Sigma_{\rm crit}\). The estimator \begin{eqnarray} \tilde{\Delta\Sigma}\left[R_{\rm p}=D_{\rm d}\theta\right] = \frac{\sum w_{ls} \left(e_{\rm tan} \Sigma_{\rm crit}\right)}{(2{\cal R})\sum w_{\rm ls}} \end{eqnarray} would then result in a minimum variance estimate if \begin{eqnarray} w_{\rm ls} = \frac{1}{\Sigma_{\rm crit}^2(e_{\rm rms}^2+e_{\rm meas}^2)} \end{eqnarray}

Also it is common to measure the signal as a function of projected radius from the center of the lensing object. In this case one has to use the critical surface density also in comoving units. \begin{eqnarray} \Sigma_{\rm crit}^{\rm com} = \frac{1}{(1+z_{\rm d})^2}\frac{c^2}{4\pi G}\frac{D_{\rm s}}{D_{\rm d} D_{\rm ds}} \end{eqnarray} The distances here are angular diameter distances, because they were used in order to convert physical distances into angles (\(D_{\rm s}, D_{\rm ds}\) to convert deflection angle to scaled deflection angle, and \(D_{\rm d}\) to convert \(\theta\) to the physical impact parameter). Thus we will use \begin{eqnarray} \tilde{\Delta\Sigma}^{\rm com}\left[R_{\rm p}=(1+z_{\rm d})D_{\rm d}\theta\right] &=& \frac{\sum w_{ls} \left(e_{\rm tan} \Sigma^{\rm com}_{\rm crit}\right)}{(2{\cal R})\sum w_{\rm ls}}\\ \end{eqnarray}

Weak lensing magnification

Until now we have talked about the effects of weak lensing on the shapes of galaxies. However weak lensing also magnifies the flux coming from an object. At the same time, weak lensing conserves surface brightness, which means the solid angle corresponding to the images seen in the sky is in reality smaller. The surface density of sources in a given region with magnification \(\mu\) will thus \begin{eqnarray} n(>S, \vec\theta, z) = \frac{n(>S/\mu(\vec\theta, z))}{\mu(\vec\theta, z))} \end{eqnarray} The numerator implies in area of magnification greater than unity, we expect to see galaxies with lower fluxes, while the denominator represents the fact that the increase in observed effective area, implies the observed number of galaxies from a smaller area will get spread out over a larger area and therefore the observed surface density should be lower.

For power law source surface densities, \(n_0(>S) \propto S^{-\alpha}\). Thus we obtain

\begin{eqnarray} n(>S) &=& \frac{n_0(>S)}{\mu^{-\alpha}\mu} \\ \frac{n(>S)}{n(>S_0)}&=& \mu^{\alpha-1}\\ \end{eqnarray}

Given the two competing effects the magnification depends upon the value of the local slope of the number counts at the limiting magnitude of the survey. For values of \(\alpha>1\), the source counts are enhanced in regions of \(\mu>1\). Such magnification effects have also been seen in real observations, especially the dependence of this ratio on the flux limit has been used to verify that the behaviour seen is indeed due to magnification.

In the weak lensing regime, the magnification is given by \begin{eqnarray} \mu = \frac{1}{(1-\kappa)^2+\gamma^2} \approx 1 + 2\kappa \end{eqnarray} and therefore the enhancement/depletion is given by \begin{eqnarray} \frac{n(>S)}{n(>S_0)}&\approx& 1 + 2(\alpha-1)\kappa \end{eqnarray}

\begin{eqnarray} \gamma' &=& \gamma e^{-2i\theta}\\ \gamma'_1 + i \gamma'_2 &=& (\gamma_1 + i \gamma_2) (\cos2\theta - i\sin2\theta)\\ \gamma'_1 &=& \gamma_1 \cos2\theta + \gamma_2 \sin 2\theta \\ &=& |\gamma|(\cos2\phi \cos2\theta + \sin2 \phi \sin 2\theta )\\ &=& |\gamma| \cos(2(\phi-\theta)) \\ \gamma'_2 &=& -\gamma_1 \sin2\theta + \gamma_2 \cos 2\theta \\ &=& |\gamma|(-\cos2\phi \sin2\theta + \sin2 \phi \cos 2\theta )\\ &=& |\gamma| \sin(2(\phi-\theta)) \end{eqnarray}