A&A 485, 363-376 (2008)
DOI: 10.1051/0004-6361:20078631

## Weak lensing goes bananas: what flexion really measures

P. Schneider1 - X. Er1,2

1 - Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
2 - Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

Received 7 September 2007 / Accepted 27 February 2008

Abstract
In weak gravitational lensing, the image distortion caused by shear measures the projected tidal gravitational field of the deflecting mass distribution. To lowest order, the shear is proportional to the mean image ellipticity. If the image sizes are not small compared to the scale over which the shear varies, higher-order distortions occur, called flexion. For ordinary weak lensing, the observable quantity is not the shear, but the reduced shear, owing to the mass-sheet degeneracy. Likewise, the flexion itself is unobservable. Instead, higher-order image distortions measure the reduced flexion, i.e., derivatives of the reduced shear. We derive the corresponding lens equation in terms of the reduced flexion and calculate the resulting relation between brightness moments of source and image. Assuming an isotropic distribution of source orientations, estimates for the reduced shear and flexion are obtained and then tested with simulations. In particular, the presence of flexion affects the determination of the reduced shear. The results of these simulations yield the amount of bias of the estimators as a function of the shear and flexion. We point out and quantify a fundamental limitation of the flexion formalism in terms of the product of reduced flexion and source size. If this product increases above the derived threshold, multiple images of the source are formed locally, and the formalism breaks down. Finally, we show how a general (reduced) flexion field can be decomposed into its four components. Two of them are due to a shear field, carrying an E- and B-mode in general. The other two components do not correspond to a shear field, and they can also be split up into corresponding E- and B-modes.

Key words: gravitational lensing - galaxies: evolution - galaxies: statistics - cosmology: diffuse radiation

### 1 Introduction

Weak gravitational lensing provides a powerful tool for studying the mass distribution of clusters of galaxies, as well as the large scale structure in the Universe (see Mellier 1999; Bartelmann & Schneider 2001; Refregier 2003; Schneider 2006; Munshi et al. 2006, for reviews on weak lensing). It has led to constraints on cosmological parameters, such as those characterizing structure formation and the mass density of the Universe.

In weak lensing, one employs the fact that the image ellipticity of a distant source is modified by the tidal gravitational field of the intervening matter distribution. Based on the assumption that the orientation of distant sources is random, the ellipticity of each image yields an unbiased estimate of the line-of-sight integrated tidal field, usually called shear in lensing. The shear thus carries information about the properties of the mass distribution. Formally, the shear is described in terms of a first-order expansion of the lens equation, i.e., the locally linearized lens equation. This yields a valid description of the mapping from the image to the source sphere, as long as the images are small compared to the length scale on which the shear varies. However, this linear approximation breaks down for larger sources, or does so in regions of the lens plane where the shear varies rapidly. The most visible failure of the linearized lens equation is the occurrence of giant arcs, which in most cases correspond actually to multiple images of a background source. To model them, the full lens equation needs to be studied. However, there is an intermediate regime where the linearized lens equation breaks down, although (locally) no multiple images are formed - the arclets regime. Arclets are fairly strongly distorted images of background sources (Fort et al. 1988; Fort & Mellier 1994), though they do not correspond to multiple images.

Arclets are the most natural application for flexion. Flexion has been introduced by Goldberg & Bacon (2005) and Bacon et al. (2006), and describes the lowest-order deviation of the lens mapping from its linear expansion; it has also been termed sextupole lensing'' and been treated by Irwin & Shmakova (2005, 2006, and references therein). It corresponds to the derivative of the shear; in combination with a strong shear, this can deform round images into arclets, giving rise to images which resemble the shape of a banana. In their original paper, Goldberg & Bacon (2005) considered only a single component of flexion which, however, only provides an incomplete description of shear derivatives. In Bacon et al. (2006), the need for a second flexion component was recognized.

In the first part of this paper, we present the general theory of flexion; in contrast to earlier work, we explicitly consider the quantities that can be actually observed, by accounting for the mass-sheet degeneracy (Falco et al. 1985; Gorenstein et al. 1988). That is, a change of the surface mass density of the form leaves the shape of all observed images invariant. In usual weak lensing, this is accounted for by recognizing that not the shear can be obtained from observations, but only the reduced shear (Schneider & Seitz 1995). The difference of shear and reduced shear is typically small, in particular in applications of cosmic shear, since along most lines-of-sight, the value of is very much smaller than unity. In applications of flexion, however, we expect that the surface mass density no longer is very small; for instance, arclets occur in the inner parts of clusters where . Therefore, the difference between shear and reduced shear can no longer be neglected. Gradients of the shear are not directly observable; only derivatives of the reduced shear are, and thus we define the (reduced) flexion in terms of derivatives of g. In Sect. 2.1 we briefly recall the irreducible tensor components which are defined in term of their behavior under rotations of the coordinate system. It turns out that a complex notation for these tensor components is very useful. In Sect. 2.2 we expand the lens equation to second order, before deriving the corresponding lens equation (and relation for the local Jacobian) which is invariant under mass-sheet transformations. The second-order term in this lens equation is fully described by our reduced flexion components G1 and G3.

As is known from usual weak lensing studies, a measured shear is not necessarily accounted for by an (equivalent) surface mass density. Since the shear is a two-component quantity, it has one degree of freedom more than the field. Therefore, shear fields are decomposed into E- and B-modes (Crittenden et al. 2002; Schneider et al. 2002), where the former are due to a  field, whereas the latter describes the remaining (curl'') part. A similar situation occurs in flexion, which has four components. Therefore, in Sect. 3 we consider the decomposition of a general flexion field into contributions due to the gradient of the shear and those not related to the shear field. The former one can then be further subdivided into flexion resulting from an E- and B-mode shear field. We carry out this decomposition for the flexion as well as for the reduced flexion.

In Sect. 4 we then define brightness moments of sources and images and derive the transformation laws between them. This approach is very similar to the HOLICs approach developed by Okura et al. (2007a,b) and later considered by Goldberg & Leonard (2007), except that we explicitly write all relations in terms of the reduced shear and the reduced flexion. Generalizing the usual assumption that the expectation value of the source ellipticity is zero - due to the phase averaging over source orientations - to the expectation values of all source shape parameters which are not invariant under coordinate rotations (as appropriate for a statistically isotropic Universe), we obtain in Sect. 5 estimates for the reduced shear and reduced flexion in terms of the brightness moments of the images. In Sect. 6 we perform a number of numerical experiments to test the validity of our approach and the accuracy of the estimators derived. In particular, we point out that there is a fundamental limit where the theory of flexion has to break down - the second-order lens equation is non-linear and will in general have critical curves, leading to multiple images of the source (or parts of it). If the source is cut by a caustic, different parts of it will have different numbers of images, and the assumption of random source orientation (which underlies all weak lens applications) will break down - the caustic introduces a preferred orientation into the source plane. In Appendix B we provide a full classification of the critical curves of the second-order lens equation and use these results in order to obtain the maximum source size (for given values of the reduced flexion) for which the flexion concept still makes sense. We discuss our results in Sect. 7.

### 2 Complex lensing notation

Like in many other instances in weak lensing, flexion is best described by using complex notation, which we shall briefly introduce next and which will be used for vectors and tensor components throughout this paper.

#### 2.1 Irreducible tensor components

For a two-dimensional vector , we define the complex number . Under rotations of the coordinate system by an angle , x gets multiplied by the phase factor . For a tensor of second rank, whose Cartesian components are Qij, we define the complex numbers and . A rotation of the coordinate systems by an angle multiplies Q2 by the phase factor , whereas Q0 remains unchanged. This is most easily seen by considering that the prototype of a second rank tensor is Qij=xi yj, where and are vectors; the foregoing statements are then obtained by noting that the complex numbers xy and x* y are multiplied by and 1, respectively, under coordinate rotations. According to this transformation behavior, we shall loosely speak about Q0 as a spin-0 quantity, whereas x and Q2 are spin-1 and spin-2 quantities, respectively.

We shall be dealing only with totally symmetric tensors. If Qijis symmetric, then

 (1)

If Tijk is a symmetric third-rank tensor, we define its spin-3 and spin-1 components as

 T3 := T1 := (2)

Furthermore, if Fijkl denotes a symmetric fourth-rank tensor, we decompose it into its spin-4, spin-2 and spin-0 components, respectively,

 F4 := F2 := F0 := F1111+2F1122+F2222 . (3)

Apart from notational simplicity, the complex lensing notation provides a check for the validity of equations. In a valid equation, each term has to have the same spin. The product of a spin-m and a spin-n quantity has spin m+n. The complex conjugate of a spin-nquantity has spin -n.

#### 2.2 Second-order expansion of the local lens equation

In weak lensing, the lens equation is linearized locally by writing the relative source coordinate in terms of the image position as , where is the deflection potential, indices separated by a comma denote partial derivatives with respect to , and summation over repeated indices is implied. Note that the form of this equation implies that the origin of the lens plane, , is mapped onto the origin of the source plane. The surface mass density and the complex shear at the origin are given in terms of the deflection potential, , , being spin-0 and spin-2 fields, respectively. In our complex notation, the locally linearized lens equation reads

 (4)

We next generalize this result to a second-order local expansion of the lens equation, which in Cartesian coordinates reads /2. The third-order derivatives of are related to the gradient of and . To write these derivatives also in complex form, we define the differential operators

 (5)

The differential operator turns a spin-n field into a spin-(n+1) field, whereas reduces the spin by one unit. One finds, for example,

 = = = (6)

and we recognize the combinations of third derivatives of which form the spin-1 and spin-3 combinations defined in (2). The final relation in (6) is the relation between first derivatives of  and  found by Kaiser (1995), here expressed in compact form. It expresses the fact that the third-order derivatives of the deflection potential can be summarized in the spin-3 field and the spin-1 field , where we introduced the usual notation for the two flexion quantities. The second-order lens equation in our complex notation then reads

 (7)

Since this is no longer a linear equation, a source at may have more than one image. In fact, up to four images of a source can be obtained, as can be seen for the special case of and by placing the source at . In this case, if we set , then one solution is , and the other three are , with , and . Of course, the origin for the occurrence of these solutions lies in the fact that is a spin-3 quantity. We shall later need the Jacobian determinant of this lens equation, which is

 = = (8)

where the first expression is just the first-order Taylor expansion of the Jacobian around the origin, and in the second step we made use of the relation . We point out that (8) is not the full Jacobian of the lens Eq. (7), but only its first-order expansion; the full Jacobian contains quadratic terms in . We will return to this important issue further below.

#### 2.3 Accounting for the mass-sheet degeneracy

The observables of a gravitational lens system are unchanged if the surface mass density is transformed as (Gorenstein et al. 1988). In the case of weak lensing, the shape of images is unchanged under this transformation (Schneider & Seitz 1995). Because of this mass-sheet degeneracy, the shear is not an observable in weak lensing, but only the reduced shear . In fact, since we expect that the most promising applications of flexion will come from situations where is not much smaller than unity, the distinction between shear and reduced shear is likely to be more important for flexion than for the usual weak lensing applications. Hence, at best we can expect from higher-order shape measurements to obtain an estimate for the reduced shear and its derivatives. For this reason, we shall rewrite the foregoing expressions in terms of the reduced shear.

The mass-sheet transformation is equivalent to an isotropic scaling of the source plane coordinates. Hence, we divide (7) by to obtain

 (9)

We will now express the coefficients in the lens Eq. (9) in terms of the derivatives of the reduced shear,

 (10)

The expression for in terms of the reduced shear and its derivatives has been derived by Kaiser (1995); in our notation it reads

 (11)

The expression for the derivative of in terms of the reduced shear can be easily obtained from differentiating the definition ,

 (12)

The derivatives G1,3 of the reduced shear are those quantities we can hope to observe; to distinguish them from and , one might call G1,3 the reduced flexion.

The Jacobian determinant of the mapping between the image position and the rescaled source position then becomes

 (13)

is a spin-1 quantity. Again, (13) is valid only to linear order in . Note that a similar equation for the determinant was obtained in Okura et al. (2007a; their Eq. (A1)), but they consider only the case of ; this has also consequences for the relations between source and image brightness moments, to be derived further below.

### 3 Compatibility relations

#### 3.1 Compatibility relations for flexion

Flexion has a total of four components, namely the real and imaginary parts of and . A measurement of flexion will thus yield four components, and we might ask whether these components are independent. We recall a similar situation in shear measurements. The shear has two components; on the other hand, the shear is defined as second partial derivatives of the deflection potential, which is a single scalar field. Therefore, the two shear components cannot be mutually independent if they are due to a gravitational lensing signal. Of course, the measured shear is not guaranteed to satisfy the condition that the two shear components can be derived from a single scalar deflection potential, since observational noise or intrinsic alignments of galaxies may affect the measured shear field. Therefore, one has introduced the notion of E- and B-modes in shear measurements (Crittenden et al. 2002). The E-mode shear is the one that can be written in terms of a deflection potential, whereas the B-mode shear cannot.

Formally, the E- and B-mode decomposition can be written in terms of a complex deflection potential and a complex surface mass density (Schneider et al. 2002). Each component of satisfies its own Poisson equation, , . Making use of this decomposition, the shear becomes

 = = (14)

The distinction between E- and B-mode shear can be obtained by considering second partial derivatives of the shear components. Taking the derivative of (14), one obtains
 = = (15)

which can be expressed in more compact form as

 (16)

A further derivative yields for the components
 = = (17)

However, it is easier to consider directly the complex derivative of , from which we obtain

 (18)

Thus, if the shear field is a pure E-mode field, is real. An imaginary part of is due to a B-mode field. This then yields the local distinction between E- and B-mode shear.

Since the flexion has four components, whereas the lens can be described by a single scalar field, we expect that there are three constraint relations a flexion field has to satisfy if it is due to a lensing potential. In fact, even if we leave the shear field arbitrary (that is, even if we allow it to be composed of E- and B-modes), then we expect two constraint equations, since the flexion field has two components more than the shear field. These constraint equations are easy to obtain. First, if the flexion field is due to a shear field, then we have

 (19)

where we defined the spin-2 quantity . It may describe contributions to the flexion which are not caused by a shear field, such as due to noise, intrinsic source alignments or higher-order terms (such as lens-lens coupling) in the propagation equation for light bundles. As a spin-2 field, a non-zero can be decomposed into its E- and B-modes. If , then the spin-3 flexion is completely determined by the spin-1 flexion up to an additive constant, as can be best seen in Fourier space, for which (19) yields . Second, if the flexion field is solely caused by a gravitational lens effect, i.e., by a pure E-mode shear field, then is real, i.e.,

 (20)

Thus, flexion from a pure E-mode shear field is characterized by the three constraint equations and , where the former is a two-component equation.

#### 3.2 The axially-symmetric case

To illustrate these compatibility relations, we consider the case of an axi-symmetric flexion field. For that, we introduce polar coordinates ; hence, in this subsection only, is the radial coordinate, not a complex number. The gradient operators then become

We first assume that the flexion derives from a shear field, which in the axi-symmetric case takes the form . In the case of a pure E-mode shear, is real, whereas for pure B-modes, is imaginary. The two flexions then read
 = = (21)

A further differentiation then yields the result that

so that the function defined in (19) vanishes, .

 Figure 1: The four different flexion fields discussed in the text. The upper left (right) panel shows the flexion corresponding to an axially-symmetric E-mode (B-mode) shear field, where arrows indicate the spin-1 flexion and the skeletons the spin-3 flexion component. In the lower left (right) panel, the flexion fields are displayed which are not due to a shear field, but a non-zero E-mode (B-mode) field. Open with DEXTER

If flexion does not derive from a shear field, then ; for example, this is the case if , which we shall take in the following. Owing to their spin properties, we can write

which then leads to the differential equations

with the solutions

 (22)

where F0 and G0 are constants of integration. We further can distinguish between the cases of being an E-mode field, in which case is real, or a B-mode field, when is imaginary, or a superposition of both.

As an explicit example, we consider the isothermal case. For a singular isothermal sphere with Einstein radius , one then has

 (23)

A further differentiation then shows that

again confirming that . The corresponding case for a B-mode shear field is obtained by multiplying all expressions in (23) by .

To obtain a similar example for the case that flexion is not derived from a shear field, we choose first a pure E-mode spin-2 field for ,

By appropriately choosing the integration constants in (22), the flexions then become

 (24)

Thus, the flexion fields are very similar to those given in (23), except that the relative signs are different. An analogous case for a pure B-mode is obtained by multiplying the foregoing expressions by . A graphical illustration of the four different cases is provided in Fig. 1.

#### 3.3 Compatibility relations for reduced flexion

Turning now to the reduced flexion, the compatibility equations can be obtained as follows. First, if the flexion is due to a shear field, we have

 (25)

as follows from the definition (10) of the two flexion components in terms of the reduced shear. Again, if this equation is satisfied, G3 is completely determined by G1, up to an additive constant. Second, if the flexion is caused by a pure E-mode shear, i.e., if the shear is due to a real surface mass density, then we employ the quantity , which is real and invariant under mass-sheet transformations, up to an additive constant. Therefore, must be real. We find:

 (26)

so that a flexion coming from an E-mode shear field satisfies K2=K2*.

We point out that the foregoing relation suggests a natural way to use flexion for finite-field mass reconstructions in weak lensing. Seitz & Schneider (2001) formulated the finite-field mass reconstruction from measured reduced shear in terms of a von Neumann boundary value problem for , whose solution determines K up to an additive constant. The source'' for was determined by the reduced shear and its derivatives, and is given by (26). In Seitz & Schneider (2001), the derivatives of the reduced shear were obtained by finite differencing of g. If flexion is measured, one can replace the source'' for by a weighted sum of the differentiated reduced shear field and the combination (K2+K2*)/2 of the flexion field, with the weights chosen according to the estimated noise properties of both contributions.

### 4 Brightness moments of source and image

We consider an image of a source, and denote the brightness distribution of the source by . Since surface brightness is conserved by lensing, the brightness distribution of the image is . Since the scaling of the source plane is unobservable, we shall only work in the following in terms of the scaled source plane coordinates, and therefore drop the hat on , as well as on .

We define the origin of the image (or lens) plane as the center-of-light of the image under consideration, i.e. we require

 (27)

Let be a function of the source coordinate; we define the operator as

 = (28)

where here and in the following, we use the linear approximation for . In particular, setting F=1, one finds that

 = (29)

since first-order moments of the light distribution in the lens plane vanish, due to our choice (27) of the coordinate system. Here, S is the flux of the lensed image, so that , as usual, where is the Jacobian at the origin .

#### 4.1 Centroid shift

The origin of the coordinates in the source plane is the image of the origin in the lens plane as mapped with the lens equation. In particular, this does not coincide with the center-of-light of the source, which is given by , or

 = (30)

Expanding the integrand, we note that terms linear in vanish, due to (27). Defining the second-order brightness moments of the image in the form

 (31)

we obtain for the source centroid shift

 = (32)

We now write these equations in a more compact form; for this, we define the matrix by , where the T'' denotes the transpose of the matrix. Then,

 (33)

where the coefficients of are given by

 b1 = b3 = b4 = (34)

The centroid shift in the source plane is thus given by the product of the derivatives of the reduced shear (expressed by G1 and G3) and the area of the image, which is proportional to Q0 and Q2. Of course, since the reduced shear and its derivatives are not directly observable, the centroid shift in unobservable as well. To get an order-of-magnitude estimate of , we assume that the source has a linear angular size , consider the reduced shear to be of order unity, and let be the angular scale on which the reduced shear varies. Then,

 (35)

#### 4.2 Transformation of second-order brightness moments

Next we consider the second-order brightness moments of the source, defined as and . By defining the third-order brightness moments of the image through

 (36)

we obtain

 = (37)

 (38)

Note that is real. In a more compact notation, (37) reads

 (39)

where the matrix has coefficients

 a1 = a2 = a3 = a4 = (40)

#### 4.3 Transformation of third-order brightness moments

We now define the third-order brightness moments of the source, separated into a spin-3 and a spin-1 component,

 = (41)

where we used that and . Similarly, we obtain

 = (42)

Defining the fourth-order brightness moments of the image by

 F0 = F4 = (43)

where Fn is a spin-n quantity, we obtain for the third-order moments of the source:

 (44)

where the matrix is defined by its transpose . The elements of are

 = T3*-3 g* T1*+3 g*2 T1- g*3 T3 ; = -g T3* +(1+2 g g*)T1* - g*(2 +g g*)T1 +g*2 T3 ; = (45)

where the last two relations are obvious. The matrix is given explicitly in Appendix A; each of its elements consists of a sum of terms proportional to fourth-order brightness moments, Fn, and terms proportional to squares of second-order brightness moments. Okura et al. (2007a) and Goldberg & Leonard (2007) have derived expressions similar to (44), though using a number of simplifying assumptions (such as ) and (in the latter paper), not considering the reduced flexion.

We will now consider the order-of-magnitudes of the various terms appearing in (39) and (44). Assuming that the third-order moments of the sources are small, then the third-order moments of the image are given by the product of and . With and , we find that . To get an estimate of the size of the various terms in (39), we note that the first three terms on the right-hand side (those proportional to the Qn) are of order , whereas and . Hence, the last two terms are of equal magnitude in general, each of them being smaller than the first three terms by a factor . Only if the source is of the same order as the scale over which the reduced shear varies do the last two terms in (39) contribute. In (44), we have neglected the terms , since they are two powers of smaller than the terms written down.

### 5 Shear and flexion estimates

#### 5.1 Estimate of the reduced shear

We see that (44) is a linear equation for , which can thus be solved,

 (46)

Inserting this into (39) then yields

 (47)

We are thus left with a single complex equation for g, which contains the observable brightness moments of the image, as well as the unobservable brightness moments of the source. This equation can be used to estimate the reduced shear if we make assumptions concerning the properties of the source brightness moments. We assume that the sources are oriented randomly, which implies that all quantities with spin unequal zero have a vanishing expectation value. That is, we set , , to arrive at

 (48)

where we have indicated that the right-hand side depends on the reduced shear (in fact it does so in a very complex manner). However, since we have argued above that the terms on the left-hand are much larger than those on the right-hand side, an iterative solution of this equation is suggested. Assume the right-hand side is given, then we get the solutions

 (49)

is the complex ellipticity of the image. Obviously, there are two solutions g for a given value of Y. This situation is similar to that of ordinary'' weak lensing, where this ambiguity also occurs: as shown by Schneider & Seitz (1995), from shape measurements of background galaxies, and cannot distinguish locally between an estimate g and 1/g*=g/|g|2. The same occurs here; we therefore assume that we pick one of the two solutions, say the one corresponding to the -'' sign; this then yields for small shear . It should be stressed that flexion impacts the determination of shear from the second-order brightness moments, due to its impact on higher-order brightness moments; hence, in general the determination of shear and flexion are coupled.

We start the iteration by setting Y0=0. This yields a first-order solution for the estimate of g,

 (50)

We then use the iteration equations

 (51)

This procedure converges quickly to one of the two solutions (g,G1,G3); the other solution is obtained by taking the +'' sign in the above equations.

Of course, our approach of setting yields a biased estimator for g; this is true even in the absence of flexion (e.g., Schneider & Seitz 1995). The reason is that, although the expectation value of vanishes, the resulting estimator for g is a non-linear function of and thus biased. The bias depends on the ellipticity distribution of the sources. It should be stressed, however, that a modified definition of image ellipticity exist such that its expectation value is an unbiased estimate of the reduced shear (Seitz & Schneider 1997).

#### 5.2 Estimates for the reduced flexion

The flexion estimator is given by (46). Since the matrix contains many terms, this is a fairly complicated equation in general. A simpler estimate is obtained if we assume that the reduced shear is small, , in which case the matrix simplifies considerably - see Appendix. Furthermore, if we assume that the brightness moments of spin are much smaller than the corresponding ones with spin 0, then we find the simple relations

 (52)

If we then set the , as would be true for the expectation value, then we obtain as estimates for the reduced flexion

 (53)

Thus, the flexion is then given by the third-order brightness moments of the image, divided by a quantity that just depends on the size of the image. Similar relations to (53) have been given in Goldberg & Leonard (2007), whereas Okura et al. (2007a) obtain a different expression for G1. We will check the accuracy of (53) in Sect. 6 below.

A more accurate estimate is obtained if we consider the reduced shear as well as the ratios of non-zero spin brightness moments to zero spin moments (such as |Q2/Q0| or |F2,4/F0|) to be of order , and then expand the flexion to first order in the (small) parameter to obtain

 (54)

### 6 Numerical tests of flexion estimators

In this section we describe some simulations that we have performed in order to test the behavior of the estimators given in the previous section.

#### 6.1 Description of the simulations

We model the sources as elliptical Gaussians, truncated at three times the scale radius'' chosen such that the area of a source is independent of its ellipticity. The ellipticity of the sources follows a Gaussian distribution, with a dispersion of of R=0.4 (i.e., we use the same ellipticity distribution as in Schneider & Seitz 1995). However, for reasons explained in the next section, we truncate the intrinsic ellipticity distribution at . For each source, we map a grid of pixels from the lens plane to the source plane using the lens equation to obtain the brightness distribution in the lens plane. From this distribution, the brightness moments of the image are measured. A shift in the lens plane coordinates is applied as to satisfy (27). We then apply the shear and flexion estimators described above to the resulting brightness moments Qn, Tn and Fn. The shear and flexion estimates are then averaged over the Gaussian ellipticity distribution of the sources, in particular over their random orientation.

It should be noted that flexion is a dimensional quantity . As can be checked explicitly from Sect. 4, the way flexion appears in the equations is always with one order higher in the source (or image) size than the other terms in the equations. As an example, we consider (44); the left-hand side and the first term on the right-hand side are , whereas the coefficients of the matrix . This then implies that the accuracy of the flexion estimates does not depend on the magnitude of the flexion and the source size individually, but only on the product . Therefore, the following results are quoted always in terms of this product.

 Figure 2: Constraints on the combination of source size and reduced flexion for the validity of the concept of flexion. Each curve shows the dividing line between a circular source of limiting isophote being cut by a caustic (above the curve) or not (below the curve); in the former case, the assumptions underlying the flexion concept break down. The different curves in each panel are for different values of g, chosen as g=0.4,0.2,0.1,0.05,0, as indicated by different line types. Without loss of generality, we choose g to be real and non-negative. The four panels differ in the phase of the reduced flexion, as indicated. E.g., in the upper left panel, the phases of G1, G3 are the same as that of g. Open with DEXTER

#### 6.2 Multiple images, and the breakdown of flexion

 Figure 3: Accuracy of the estimates for reduced shear and flexion. The left panel shows contour of constant fractional error of , and , on the estimate of the reduced shear g, as a function of , where we chose g=0.05 as input value, and assumed the phases of G1, G3 to be the same as that of g. The estimate was obtained by solving the iteration equations given in Sect. 5. The right panel shows the fractional error levels at 3, 5, and 10% for the reduced flexion, as quantified by (55), where the estimate was obtained again with the iterative procedure. In both cases, we assumed circular sources. Open with DEXTER

 Figure 4: Comparison of the reduced flexion estimators (53) with the full expression (46) and the input value. The horizontal and vertical axis show . For both panels, we take g=0.05, and G3=0 (G1=0) for the left (right) panel. The line indicates the input value, the plus symbols show the simplified reduced flexion estimate (53), and the crosses result from the full expression of reduced flexion (46). As can be seen from the left-hand panel, the full estimator for the reduced flexion yields a more biased result that the approximate expression (53); we have not found a reasonable explanation for this behavior. Open with DEXTER

As we mentioned before, the lens Eq. (7) can give rise to multiple images. As can be seen from the example given after (7), if the flexion is sufficiently small, all but one of these multiple images will be located at a large distance from the origin, and the central image of an extended source will be isolated. In this case, this central, or primary, image (the shape of which we measure here) is not crossed by a critical curve, and thus the source is not crossed by a caustic. The multiple images at large distances from the origin then result from the low-order Taylor expansion of the lens equation, which most likely breaks down at these image positions anyway; hence, these additional images are of no relevance. If, however, the flexion becomes sufficiently large - or if the source is large enough - this is no longer the case, and the multiple images of an extended source will merge. If that happens, the whole method of determining shear and derivatives thereof from brightness moments will break down. This can be most easily seen by considering the caustic curve cutting the source. Different parts of the source will be mapped onto a different number of image points in the lens plane, and the caustic curve introduces a direction into the situation. Hence, the assumption of an isotropic orientation of sources can no longer be employed. Mathematically, this can be seen from (28); there, the transformation between source and image plane no longer is correct if multiple images do occur. More precisely, the transformation between source and image coordinates in the calculation of the brightness moments implicitly assumes that within the limiting isophote of the primary image, the lens equation is invertible. Owing to what was said above, the condition that the central image is isolated (so that locally no multiple images occur) can be expressed solely by the products . These products approximately measure the fractional change of the reduced shear across the image of a source.

In our simulations we can check whether a critical curve crosses our central image, just by controlling the sign of the Jacobian determinant (the true one, not the linear approximation Eq. (13)). If the source size becomes too large, some points in the image will have a negative Jacobian. In the Appendix B, we consider the critical curves and caustics of the lens Eq. (9), which allows us to determine the regions in flexion space where no local multiple imaging occurs. Some examples of this are plotted in Fig. 2. Each panel shows the dividing line between parameter pairs for a circular source of limiting isophotal radius ; below the curves, no local multiple images occur, whereas for parameter pairs above the lines, the flexion formalism using moments necessarily breaks down. The different lines in each panel correspond to different values of g. The occurrence of critical curves also is the reason why we truncated the intrinsic ellipticity distribution of the sources in the simulations, since in the limit of , keeping the source area fixed, there will be orientation angles for which the source will hit a caustic.

#### 6.3 Estimates of the accuracy

We now present some results of our numerical simulation regarding the accuracy with which the reduced shear and flexion can be obtained with our moment approach. For given input values of g, G1 and G3, we either measure the brightness moments for a single circular source, or average the results over an ellipticity distribution, as described above. It should be noted that we have to deal with a 5-dimensional parameter space, namely the 3 complex parameters g, G1 and G3, minus one overall phase that can be chosen, e.g., to make g real and positive. Thus, instead of sampling the parameter space comprehensively, we only give a few selected results.

We start by considering a circular source, and determine the effect of flexion on the determination of the reduced shear. The left-hand panel of Fig. 3 shows contours of constant fractional deviation , in the flexion parameter plane. Here it is assumed that the phase of both flexion components is the same as that of g (as would be the case in an axially-symmetric lens potential). Errors of order 5% occur already for , and the fractional error increases approximately linearly with the strength of flexion (or with the source size), although it does not scale equally with both flexion components. The reason for this effect has been mentioned before - flexion affects the transformation between source and image quadrupole moments, as can be seen in (37).

In Fig. 4, we show the expectation value of the reduced flexion components, as a function of the input flexion. The expectation value has been determined by averaging over an isotropic ensemble of elliptical sources, as described before. The left and right panel show the behavior of the expectation value of G1 and G3, respectively, where the other flexion component was set to zero. The dashed curve shows the identity, the plus symbols were obtained by using the approximate estimator (53), whereas the crosses show the expectation values as obtained by employing the full expression (46), where the corresponding value of gwas obtained by the iterative process described in Sect. 5. It is reassuring that the expectation value closely traces the input value, i.e., that the estimates have a fairly small bias. Furthermore, we see that the approximate estimator (53) performs remarkably well. It is seen that the estimates for G3 behave better than those for G1. This can also be seen from the right-hand panel of Fig. 3, where we plot contours of constant fractional error

 (55)

where is the deviation of the estimate of Gn from its input value. For simplicity, we have assumed a circular source. We see that the accuracy decreases much faster with increasing G1 than with increasing G3. The reason for that may be related to the fact that the estimator of G1 is more strongly affected by the non-linearity of the equations, as can also be seen in (54).

### 7 Conclusions and further work

In this paper, we have studied the effect of flexion in weak gravitational lensing. The main results are summarized as follows:

• Owing to the mass-sheet degeneracy, flexion itself cannot be determined, but only reduced flexion. We have therefore written the second-order lens equation (which contains the derivatives of the reduced shear, i.e., flexion) as well as the relations between the brightness moments of source and image strictly in terms of the reduced shear and the reduced flexion.
• We pointed out that a general flexion field can be decomposed into a pair of components which is due to a shear field, i.e., its derivatives, and a pair of components not related to shear. The former pair can be further separated into flexion due to an E- and B-mode shear, with only the E-mode flexion expected to arise from gravitational lensing. For the second pair of components, no physical interpretation is available; if they arise in measurements, they are most likely due to noise or intrinsic shape effects of sources. General relations to separate these components are given.
• We derived the relations between low-order brightness moments of source and image, taking into account that the presence of flexion leads to a centroid shift, and it also affects the relation between second-order brightness moments - and thus the estimate of the reduced shear. Hence, the presence of flexion has an impact on the shear measurements. Starting from these moment equations, we obtained approximate estimates for the reduced shear and flexion.
• We pointed out a limit where the flexion formalism ceases to be valid, namely when the product of source size and flexion is sufficiently large that parts of the source are multiply imaged locally, i.e., where a caustic cuts through the source. We quantified this with numerical simulations, and also provided a complete classification of the critical curves of the second-order lens equation employed in flexion studies.
• We performed a number of numerical experiments to study the bias of the reduced shear and flexion estimators. However, due to the high dimensionality of parameter space, no comprehensive study was presented here. We also pointed out that only the product of flexion and source size matters in the accuracy of estimates.
The possible occurrence of critical curves in highly distorted images may provide a serious obstacle to applications of flexion. Perhaps the most promising application of flexion measurements are those in regions where the shear field varies on small scales, i.e., close to galaxies (and thus can be used for galaxy-galaxy lensing) or in the inner regions of clusters. However, if one finds a strongly distorted image of a background galaxy as in the case of the arclet A5 in Abell 370 (Fort et al. 1988), how can one be sure that it is not due to a merged double image of the source? Using flexion for studying small-scale structure in mass distributions can therefore be affected by the occurrence of multiple imaging.

Similar to the situation in shear measurements, the moment approach for flexion as presented here must be modified in several ways to be applicable to real data. First, brightness moments must be weighted in order not to be dominated by the very noisy outer regions of the image. As is known from shear measurements, such a weighting affects the relation between source and image brightness moments. Secondly, one needs to account for the effects of a point-spread function. Both of these modifications were successfully achieved for second-order brightness moments by Kaiser et al. (1995; see also Luppino & Kaiser 1997). Goldberg & Leonard (2007) consider these effects in the context of flexion. It should be noted, though, that their consideration of the PSF effects is restricted to unweighted moments, for which these effects are given by a simple convolution. In the case of weighted brightness moments, however, the PSF effects are much more subtle. Okura et al. (2007b) indeed developed a PSF correction scheme similar to that of Kaiser et al. (2005), now accounting for higher-order brightness moments of the images and the PSF; their application of this method to synthetic data and to images of the cluster A1689 is encouraging.

But even disregarding these complications, the present paper only scratches the surface in investigating estimators for reduced flexion and their properties. As mentioned before, the second-order lens equation contains five essential parameters. The bias of an estimator for reduced shear and flexion will depend on these parameters, as well as on the intrinsic ellipticity (and higher-order moments) distribution of sources. One might ask whether it is possible to find an unbiased flexion estimator, such as was possible to construct for the reduced shear. Unfortunately, we have been unable to make analytic progress: even for a circular Gaussian source, the brightness moments of the image cannot be calculated analytically. Our ray-tracing algorithm with which we conducted our numerical simulations is almost certainly sub-optimal; a more advanced method should be developed to reduce the numerical efforts in calculating brightness moments. Beside the bias, it would be interesting to calculate the variance of the various estimators, or more precisely, their covariance.

It may turn out that measurements of flexion, and PSF corrections, are more conveniently done with shapelets, as was originally considered by Goldberg & Bacon (2005), Bacon et al. (2006) and Massey et al. (2007a). Even if this turns out to be the case (see Leonard et al. 2007, for an application of flexion measurements in the galaxy cluster A1689), the moment approach provides a more intuitive picture of the effects of flexion. In addition, the weak lensing community has profited substantially from the existence of several different methods to measure shear (see Heymans et al. 2006; Massey et al. 2007b, for the first results of a comprehensive Shear TEsting Programme, in which these various methods are studied and compared); therefore, the development of different techniques for measuring flexion will certainly be of interest once the flexion method is put to extensive use.

Acknowledgements
We thank Jan Hartlap and Ismael Tereno for useful comments on this paper. This work was supported by the Deutsche Forschungsgemeinschaft under the project SCHN 342/6-1 and the TR33 The Dark Universe''. X.E. was supported for this research through a stipend from the International Max-Planck Research School (IMPRS) for Radio and Infrared Astronomy at the University of Bonn.

### Appendix A: The matrix C

 Figure A.1: The critical curves ( left-hand panel) and caustics ( right-hand panel) of the lens Eq. (9) for the cases of hyperbolic critical curves, as described in Sect. B.2. The parameters chosen here are g=0.05, , . A circular source is mapped onto two images, as indicated. If the source size were increased, it would hit the caustic, the two images would merge, and the flexion concept would break down. The unit of the reduced flexion is the inverse of the unit in which coordinates are measured. Open with DEXTER

 Figure A.2: Same as Fig. A.1, but for the parabolic case, with parameters g=0.05, G1=-0.04, G3=0.112. Open with DEXTER

In this Appendix, we list the coefficients of the matrix which occurs in (44):

 4(1-gg*) C11 = -2g F2* +(9g g*-3) F0 +6 g*(1-2 g g*) F2 +g*2(5 g g*-3) F4 +6 g Q2* Q0 -12 g g* Q02 +(3 -9 g g*) Q2* Q2 + 6 g* (4 g g*-1) Q0 Q2 +3 g*2(1- 3 g g*) Q22 4(1-gg*) C12 = 5 g F4* -2(5+6 g g*) F2* +9 g*(3+g g*) F0 -2 g*2(12+g g*) F2 +7 g*3 F4 -9 g Q2*2 +6(3 +4 g g*) Q2* Q0 -12 g*(3+g g*) Q02 -3 g*(5+3 g g*) Q2* Q2 +6 g*2(8+g g*) Q0 Q2 -15 g*3 Q22 4(1-gg*) C13 = -7 F4* +26 g* F2* -36 g*2 F0 +22 g*3 F2 -5 g*4 F4 + 15 Q2*2 -54 g* Q2* Q0 +48 g*2 Q02 +24 g*2 Q2* Q2 -42 g*3 Q0 Q2 +9 g*4 Q22 4(1-gg*) C14 = -2 g* F4* +6 g*2 F2* - 6 g*3 F0 +2 g*4 F2 +6 g* Q2*2 -18 g*2 Q2* Q0 +12 g*3 Q02 +6 g*3 Q2* Q2 -6 g*4 Q0 Q2 4(1-gg*) C21 = 2 g2 F2* -6 g2 g* F0+[4 g g*(1+g g*)-2] F2 +2 g*(1-2 g g*) F4 -6 g2 Q2* Q0 +4 g (1+2 g g*) Q02 +6 g2 g* Q2* Q2 +[2-4 g g*(3+2 g g*)] Q0 Q2+2 g* (4 g g*-1) Q22 4(1-gg*) C22 = -5 g2 F4* +2 g (7+4 g g*) F2* -3[3+g g*(8+g g*)] F0+2 g*(8+5 g g*) F2 -7 g*2 F4 + 9 g2 Q2*2 -2 g (13+8 g g*) Q2* Q0 +4[3+g g*(8+g g*)] Q02 +[5+g g*(16+3 g g*)] Q2* Q2 -2 g*(16+11 g g*) Q0 Q2 +15 g*2 Q22 4(1-gg*) C23 = 7g F4* -2(4 +9 g g*) F2* +3 g*(7+5 g g*) F0 -2 g*2(9+2 g g*) F2+5 g*3 F4 -15 g Q2*2 +(16+38 g g*) Q2* Q0 - 4 g*(7 +5 g g*) Q02 - g*(13+11 g g*) Q2* Q2 +2 g*2(17+4 g g*) Q0 Q2 -9 g*3 Q22 4(1-gg*) C24 = (3 g g*-1) F4*-6 g g*2 F2* +3 g*2(1 + g g*) F0 - 4 g*2(2+g g*) Q02-3 g*2(1+g g*) Q2* Q2 +6 g*3 Q0 Q2.

The other eight elements follow trivially from the foregoing ones, since the second half of the matrix is just the complex conjugate one of the first half, i.e., C44=C11*, C34=C21* etc., or in general, Cij=C5-i,5-j*.

 Figure A.3: Same as Fig. A.1, but for the elliptical case, with parameters . Open with DEXTER

### Appendix B: Critical curves and caustics

In this Appendix we consider the critical curves of the lens Eq. (9), with the goal of finding the pricipal limits of the applicability of the flexion formalism - which necessarily breaks down if parts of the source are multiply imaged. For this, we need to derive the full Jacobian, which can most easily be obtained from considering and as independent variables, and then use , , which can be inverted to yield , . With these relations, one finds that . Carrying out these derivatives, the Jacobian becomes

 (B.1)

with
 (B.2)

Note that A is a spin-2 quantity, whereas B is a real scalar, i.e., has spin-0. In the generic case, the critical curves ( ) are conical sections, which may be degenerate, though. We will now perform a complete classification of cases that can occur, as well as to derive the critical curve(s) in parametric form; the caustics are then obtained by inserting the parametric form of the critical curves into the second-order lens equation. As we shall see, the type of conical section is determined, amongst other parameters, by the discriminant

 (B.3)

#### B.1 Zero discriminant

We start with the case that , which implies B2=4 A A*, or . The case A=0=B either implies that , in which case also so that no critical curves occur, or that , for which in general. In this case, the critical curve is a straight line, satisfying . As can be seen by inspection, it reads

 (B.4)

If , the phase of A is defined. Since it is a spin-2 quantity, we write . Furthermore, we introduce the rotation . Then the equation for the critical curve reads, after dividing (B.1) by |A|,

 (B.5)

and the sign on the left-hand side of the equation depends on the sign of B, where we used . The parametric form of the critical curve, which takes the form of a parabola, can then be written as
 (B.6)

where the first (second) equation applies for B>0 (B<0). Note that the parabola degenerates into a straight line if is real (for B>0) or purely imaginary (for B<0).

#### B.2 Non-zero discriminant

If , we can perform a translation to eliminate the linear term in . Hence we define and choose such that terms linear in vanish. We then obtain for and for the critical curve condition

 (B.7)

with

 C = = (B.8)

where the second step was obtained by inserting the expression for in terms of the 's, and in the final one we defined Vas the expression in the parenthesis.

As the first case, we consider A=0 and (the case A=0=Bwas treated above), which implies that and . The equation for the critical curve then reduces to . Furthermore, , and C=-1. Thus, the critical curve is a circle of radius and center , or , .

We now consider the case ; then the phase of A is defined, as used before. Introducing a rotation by defining , the equation for the critical curve becomes

 (B.9)

The presence and topology of critical curves now depends on the signs of and C. We first consider the case C=0; then, if , no critical curves occur, except for the isolated point x=0. If , the critical curves are two straight lines, as can be obtained from (B.7): inserting the ansatz , one obtains . Thus, the critical curves are parametrized as
 (B.10)

For the case of , the consideration of (B.9) yields the result that for , the critical curves consist of two hyperbolae. From (B.8) we see that negative implies C>0. Also note that implies that 2|A|-B>0, 2|A|+B>0. The critical curves then read
 (B.11)

For the other case, , we find from (B.8) that C<0. If , we then see from (B.9) that no critical curves exist. If , which in particular implies B<0, the critical curve is an ellipse parametrized as
 (B.12)

This concludes the classification of critical curves of the lens Eq. (9). The caustics are obtained by inserting the parametrized form of the critical curves into the lens equation. In order to see whether a critical curve cuts through the primary image of a circular source of outer isophotal radius , we calculate the minimum value of along the caustics. If , the image is not cut by a critical curve. For an elliptical critical curve, the maximum source size allowed is ; these values are plotted in Fig. 2. In the cases where two critical curves exist (e.g., two straight lines or hyperbolae), the situation is slightly more complicated. Consider, e.g., the case of two straight critical curves. Only those sections of them that are closer to the origin are relevant for this consideration, since if the primary image of the source is not cut by these closer sections of critical curves, it will still be an isolated image; the caustics coming from the outer sections of the critical curves correspond to multiply imagedsource sections of secondary images. Accounting for this complication, the maximum sources size have been obtained, as plotted in Fig. 2.