{'fraction_non_alphanumeric': 0.04324270835910249, 'fraction_numerical': 0.0170447517502412, 'mean_word_length': 4.469354620484373, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Shape, shear and flexion II - Quantifying the flexion formalism for extended sources with the ray-bundle method

Flexion-based weak gravitational lensing analysis is proving to be a useful adjunct to traditional shear-based techniques. As flexion arises from gradients across an image, analytic and numerical techniques are required to investigate flexion predictions for extended image/source pairs. Using the Schwarzschild lens model, we demonstrate that the ray-bundle method for gravitational lensing can be used to accurately recover second flexion, and is consistent with recovery of zero first flexion. Using lens plane to source plane bundle propagation, we find that second flexion can be recovered with an error no worse than 1% for bundle radii smaller than {\Delta}{\theta} = 0.01 {\theta}_E and lens plane impact pararameters greater than {\theta}_E + {\Delta}{\theta}, where {\theta}_E is the angular Einstein radius. Using source plane to lens plane bundle propagation, we demonstrate the existence of a preferred flexion zone. For images at radii closer to the lens than the inner boundary of this zone, indicative of the true strong lensing regime, the flexion formalism should be used with caution (errors greater than 5% for extended image/source pairs). We also define a shear zone boundary, beyond which image shapes are essentially indistinguishable from ellipses (1% error in ellipticity). While suggestive that a traditional weak lensing analysis is satisfactory beyond this boundary, a potentially detectable non-zero flexion signal remains.


INTRODUCTION
Weak gravitational lensing provides one of the most direct probes of the matter distribution of the Universe as it is independent of both the dynamical state and the nature of the matter. Building on the pioneering attempts by Valdes et al. (1983) and Tyson et al. (1984) to measure coherent changes in the shapes of background galaxies due to a foreground lens population, weak lensing techniques have now come of age. In recent years, weak lensing has successfully been applied in the cases of galaxy-galaxy lensing (e.g. Brainerd, Blandford & Smail 1996;Hudson et al. 1998;Fischer et al. 2000;Smith et al. 2001;Guzik & Seljak 2002;Hoekstra, Yee & Gladders 2004;Sheldon et al. 2004;Heymans et al. 2006;Mandelbaum et al. 2006;Parker et al. 2007;Johnston et al. 2007; Research undertaken as part of the Commonwealth Cosmology Initiative (CCI: www.thecci.org), an international collaboration supported by the Australian Research Council † cfluke@swin.edu.au The conventional mathematical basis for weak lensing analysis assumes that a shear field causes an additional ellipticity to the shape of a background source, which can be calculated by measuring the moments of the images (Kaiser 1995). However, this approach does not account for higherorder shape distortions that occur when there are strong tidal fields across the image. Recently, an additional lensing effect called flexion has been investigated as an extension to shear-based measurements (Goldberg & Natarajan 2002;Goldberg & Bacon 2005;Bacon et al. 2006) 1 . Flexion has two components denoted first and second flexion, which are a shift of the image centroid and a representation of the "arciness" of the image respectively.
One main issue associated with shear-based gravitational lensing is that galaxies are intrinsically elliptical in shape (ellipses at some inclination to the line-of-sight seen in projection), so it is necessary to disentangle lensinduced shear from the intrinsic shape. Resolved galaxies, however, are not intrisincally flexed -although systems undergoing a merger, or galaxies with substantial asymmetric sub-structure such as a large starforming region, may be mis-intepreted as flexion signals. It has been suggested that flexion may provide a stronger constraint on dark matter Bacon, Amara & Read 2010;Hawken & Bridle 2009;Leonard & King 2010), galaxy cluster mass models (Leonard, King & Wilkins 2009) and also on delensing gravitational wave signals (Shapiro et al. 2010) than shear on its own, notwithstanding the challenges in measuring flexion (e.g. Okura, Umetsu & Futamase 2007; Massey et al. 2007;Irwin & Shmakova 2006;Schneider & Er 2008).
In a previous paper we presented analytic flexion results for a range of popular mass density profiles: Schwarzschild lens, singular isothermal sphere (SIS), Navarro-Frenk-White (NFW) profile and Sérsic-like profiles (Lasky & Fluke 2009; hereafter Paper I). Our analytic solutions present a flexion formalism where we consider a two-dimensional (2D) field in the lens plane, which allows for the treatement of extended sources. In this paper we extend our previous work by considering the following key question: over what range of image/source sizes is the flexion approach valid? That is, if flexion appears as a gradient of shear across an image, which is implictly assumed to be zero for traditional weak-lensing analysis, how well can we recover flexion for extended sources?
The approach we use is via the ray-bundle method introduced by Fluke, Webster & Mortlock (1999). Here, bundles of light rays with a known initial configuration are propagated through one or more lens planes, and the deflection of each light ray is determined using the gravitational lens equation, see equation (1) below. The initial and final shapes of the bundles can now be used to obtain numerical estimates for convergence, shear, magnification, first and second flexion as a function of source size. By starting with simple lens models, we use known analytic results to test the accuracy of our approach, and hence assess its applicability to cases where analytic flexion results do not exist (e.g. asymmetric mass profiles, such as from N -body, dark matter halo simulations).
The paper is set out as follows: in Section 2, we summarise the key results for the analytic flexion formalism. In Section 3, we describe how the ray-bundle method can be used to obtain flexion along individual lines-of-sight. We present the results of our numerical testing, and determine a range of image locations and bundle radii for which the flex-1 see also Irwin & Shmakova (2005, 2006); Irwin, Shmakova & Anderson (2007) who independently derived equivalent higher-order gravitational lensing effects which they call sextupole lensing. ion approximation is valid for a Schwarzschild lens model. We demonstrate the existence of a flexion zone, which defines a physical region where the flexion formalism is valid for extended sources. In Section 4 we present calculations of the size of the flexion zone as a function of lens mass, lens and source redshift. We present our conclusions and identify future directions for this work in Section 5.

THE FLEXION FORMALISM
In this section, we summarise relevant results of the flexion formalism -for full details, see Paper I and references therein.

The gravitational lens equation
The thin-lens gravitational lens equation relates the location of a background source and its image(s) due to an intervening mass distribution: For the simple case of a single lens at the origin of the lens plane, ξi is the image location in the image plane, ηi is the source location in the source plane, andαi is the deflection angle. Throughout this paper, indices i, j, k = 1, 2 signify the two-dimensional components, with summation assumed over repeated indices, and tensors are implicity functions of the image position, ξi, unless otherwise indicated. The thinlens approximation requires that the spatial extent of the lens distribution be much smaller than the angular diameter distances between the observer and lens, DL, observer and source, DS, and lens to source, DLS. A convenient scaling of the lens equation into angular coordinates uses the substitutions βi = ηi/DS, θi = ξi/DL and αi = D LS D Sα i, so we have the compact form: We introduce a further scaling for the lens equation in section 3.4.

The lens matrix
As a mapping between the lens and source planes, equation (2) is often rewritten in terms of its Jacobian matrix: where the matrix elements are identified with the convergence, κ, and two orthogonal components of shear, γ1 and γ2. The magnification along a given line-of-sight, µ, is related to the area distortion of the lens mapping: where the total shear is γ = γ 2 1 + γ 2 2 . By setting the location of the light ray as the origin of coordinates in the lens and source planes, one finds for small changes in the position of the light ray: The interpretation of this mapping is as follows: if the gradients of the shear and convergence do not change significantly across the image, then the first-order weak lensing effect is to produce an additional ellipticity in the shape of background sources. While this approach is satisfactory for point sources, it breaks down for extended sources. This is seen even for the conceptually simplest lens model, the point mass or Schwarzschild lens, which produces increasingly "arcy" images for extended sources at small impact parameters to the lens. In order to account for the effects of gradients in the shear and convergence, a second-order Taylor expansion of the gravitational field is required: Following the approach of Bacon et al. (2006), the tensor term, D ijk , is associated with two additional lensing distortion terms By further defining F = F1 + iF2 and G = G1 + iG2, we note that the components of first (F) and second (G) flexion can be written: The components of first and second flexion only require a subset of the tensor components, D ijk , and do not depend on knowledge of D112, D121, D212 or D221. The relevance of this is discussed in Section 3.3.

Circularly symmetric lens models
For a circularly symmetric mass profile, the deflection angle is where |ξ| = ξ 2 1 + ξ 2 2 , and the projected mass, M , is for surface density, Σ(|ξ|). This enables us to re-write equation (1) as with critical surface density Defining the function we can now express the shear, convergence and first and second flexion with their explicit dependence on physical coordinates and surface mass density (see Paper I for full details):

FLEXION WITH THE RAY-BUNDLE METHOD
In this section, we describe how the ray-bundle method can be used as a numerical means of estimating flexion along a given line-of-sight. We use the Schwarzschild lens model (hereafter, SL), as it provides us with simple analytic solutions for all of the relevant lensing properties. However, since the SL has zero first flexion everywhere except at the origin, we emphasise the recovery of second flexion with the RBM, whilst demonstrating that results remain consistent with zero first flexion. We consider both backwards (lens plane to source plane) and forwards (source plane to lens plane) ray-bundle propagation, as this separately allows us to constrain the appropriate bundle radius to use and to quantify the extent of the flexion zone.

Inverse ray-tracing
While knowledge of an image location uniquely defines the source location, in general, equation (1) is not easily invertible to give all image locations for a given source position. While "brute force" solution methods can be used (e.g. Paczyński 1986), the main alternative is to use inverse raytracing in its direct form, as was introduced by Kayser, Refsdal & Stabell (1986) and Schneider & Weiss (1986;1987), or in its hierachical tree-code form (Wambsganss 1990;1999).
Here, light rays are projected backwards from the observer, through the lens plane to the source plane, which is represented by a two-dimensional grid of source pixels. The deflection of each light ray is calculated with equation (1). While ray-tracing methods are extremely well-suited to studying statistical lensing effects (e.g. the creation of magnification maps for studying probabilities of high magnification events in quasar microlensing), they are less well-suited for studying lensing effects along a given line-of-sight. If it were possible to write down an explicit analytic form for the null geodesic equation of general relativity for an arbitrary mass (i.e. lens) distribution, the optical scalar equations (Sachs 1961;Dyer & Roeder 1974) could be used to Figure 1. Shapes in the image and source planes showing nonlinear lensing effects for extended image/source pairs. For the simplest case of a Schwarzschild lens, a circular image is mapped to a "teardrop" source (top row), and a circular source is mapped to an arc-like image (bottom row). measure the changing shape of a (small) bundle of light rays as it propagates from the source to the observer. Unfortunately, specific solutions only exist for a limited number of cases, such as light propogation in Swiss Cheese "inhomogeneous" cosmological models (Harper 1991), so an alternative approach is required.

The ray-bundle method
The ray-bundle method (RBM) was developed as an alternative to grid-based inverse ray-tracing for studying gravitational lensing along specified sightlines (Fluke, Webster & Mortlock 1999). An image shape is defined in the image plane (or first lens plane in the case of multiple lens planes), comprising a central light ray, which represents the "null geodesic", surrounded by a bundle of Nray light rays. The light rays in this bundle are traced backwards to the source plane, producing a distorted source shape. As such, RBM provides a numerical analogy to the optical scalar solution, while retaining properties of the inverse ray-tracing approach, such as the use of the thin lens equation for deflection calculations.
In earlier work (Fluke, Webster & Mortlock 2002), the RBM was used to obtain magnification probability distributions for dark matter-only cosmological models. However, the RBM's emphasis on bundle shape means it is ideally suited for studying both linear and higher-order gravitational lensing phenomena (viz. shear and flexion) in detail along any given line-of-sight.
Moreover, the RBM gives a great deal of flexibility in choosing the size and shape of the initial bundle (i.e. the image), so that we can constrain the length-scale over which the flexion approximation of equation (6) is valid. If the gradient of the shear and convergence across an image is negligible, i.e. we are in the regime of equation (5), which is the truly weak-field limit of gravitational lensing, then a circular source will appear as an elliptical image, and a circular image 2 would be indicative of an intrinsically elliptical source. When tidal gradients do occur across an image, most readily due to the extended nature of image/source pairs, the 2 although unlikely to occur naturally. relationship betwen image and source shapes is less obvious. For the SL, circular sources are mapped to "arc-like" images, and circular images are due to "teardrop" shaped sources. This correspondence is shown in Fig. 1.

Flexion with the ray-bundle method
By analogy with equation (5), small changes in the position of a light ray in the image (δθ1, δθ2) and source (δβ1, δβ2) planes are related through the following equations, which are explicit expansions of equation (6): By choosing the number of rays in each bundle, it is possible to obtain a complete solution for the unknown Aij and D ijk terms along an arbitrary line-of-sight. Each light ray in the bundle experiences the tidal field around the central light ray, and so each pair of image and source rays provides a unique solution to equations (24) and (25). Since there are five unknowns we wish to solve for, we need only use six light rays in total per ray-bundle: a central light ray plus Nray = 5 rays defining the circumference. Writing these last equations in matrix form for each of the five δθin and δβ ik bundle pairs (k = 1 . . . 5): The unknown column vectors of the Aij and D ijk terms are the solutions to a set of simultaneous equations. We note that there is a degeneracy in the co-effecients of the cross-terms, δθ1δθ2, which means we cannot solve uniquely for D112, D121, D212 and D221. Instead, these terms can be obtained from symmetries of D ijk , viz. D ijk = D ikj and D12i = D21i. However, as we highlighted earlier, these four elements of the D ijk tensor do not contribute to first or second flexion [see equations (8-11)], and are instead associated with the twist and turn operators identified by Bacon & Schäfer (2009).

Image-to-source plane ray-bundles
Our first aim is to examine the (physical) scale in the image plane over which the RBM can be used to recover flexion. Whereas we used angular coordinates for the lens equation and flexion terms in sections 2.2 and 3.3, it is now more convenient to use a scaling relative to the size of the angular (point-mass) Einstein radius: in radians, for a lens with mass, M ; G and c are the gravitational constant and speed of light respectively. In terms of scaled coordinates, xi = θi/θE and yi = βi/θE, the two-dimensional lens equation for the SL model has solutions: so that the Einstein ring has a radius of 1 unit, and we can specify the bundle radius, ∆x, as a fraction of the Einstein radius. This scaling changes the Taylor expansions: and the column vectors in equations (28) and (29) once we substitute δyi = δβi/θE and δxi = δθi/θE. Note that the dimensionless nature of shear, through the Aij terms, is maintained with this coordinate change: as expected. Next, we use the relationship D ijk = F ijk +G ijk (Bacon et al. 2006), which reduces to D ijk = G ijk for the SL model, as all the F ijk = 0. In terms of components, and we use equations (22) and (23), to calculate the analytic solutions for the two components of second flexion for the SL model: Combining these last two equations implies Fortuitously, the θ −1 E term in equations (38) and (39) effectively cancels the θE coefficient of the D ijk terms in the column vectors, equation (33), so we can obtain the scale independent second flexion with the RBM without requiring a specific value of θE.
Each ray-bundle is characterised by the radius, ∆x, twodimensional location of the central ray, xi, and the number of bundle-rays, Nray. As shown in Section 3.3, we choose Nray = 5 in order to solve for the unknown Aij and D ijk terms. For convenience in plotting and making comparisons with analytic results, we use polar coordinates, (r, φ), where r = x 2 1 + x 2 2 is the radial impact parameter, and φ = tan −1 (x2/x1) is the polar angle. The analytic and RBMrecovered solutions are sampled on a grid with dimension Nr × N φ in polar space, which means that the sampling of solution space is sparser with increasing r.
For each ∆x, we have to choose an appropriate range of r values for sampling. The minimum impact parameter is rmin = 1 + ∆x(1 + ε), where ε → 0 is a small numerical offset (we use ε = 0.01). The offset avoids light rays in the image bundle overlapping the Einstein radius, thus ensuring that the 5 × 5 square matrix of δx terms is invertible. While the lens equation describes a mapping from the image plane backwards to the source plane, not all RBM images are permissable. For example, it is not possible to have a regular polygonal image straddling the Einstein ring that corresponds to a single source shape. However, the RBM does work for images either completely outside the Einstein radius, or completely within it. 3 For more complex lens models, where there are non-degenerate caustics, the solution is less obvious.
The upper limit, rmax, for the impact parameter is chosen to be Max (10, 10∆x), with the latter limit only relevant for bundle radii greater than the Einstein radius. In practice, the majority of bundle radii we use were ∆x < 0.1 (see Fig. 6). Our choice of rmax = 10 was based on numerical tests, which included varying the number of radial samples, Nr. Beyond this impact parameter, differences between the analytic and RBM-recovered flexion values were consistent with numerical (i.e. precision) limits.
then build the matrices of simultaneous equations for each ray-bundle and use Gaussian elimination with back substitution to solve for the Aij and D ijk terms, and hence Gi and Fi.
To test the matrix inversion code, we used the SL solutions of equations (34-35) and (38-39), as the analytic solutions for each (r, φ) sample. Within the limits of numerical accuracy afforded by our implementation, we find that we correctly recover the input Fi = 0 and Gi, independent of Nr, N φ and ∆x. This demonstrates the accuracy and utility of our simultaneous equation-solving code, under the proviso that there is no degeneracy in the mapping of image ray-bundles to source ray-bundles. This assumption is appropriate in the "almost weak" regime where flexion acts, and by ensuring that no part of the source crosses the caustic point for the Schwarzschild lens.
Next, we use the scaled lens equation, equation (29), to deflect the individual light rays in each bundle, so that we have pairs of image and source bundle shapes, and build the matrices for each bundle. Inverting the matrices, we obtain RBM-estimates for second flexion and use recovered first flexion results as a consistency check on numerical effects (see below).
In Fig. 2, we plot an indicative second flexion surface, G1(r, φ)θE; second flexion values must be scaled by the angular Einstein radius (in radians) for comparison with a particular system. The spin-3 oscilliatory nature of second flexion is visible, with the amplitude increasing towards the origin. The surface is plotted over a limited range in impact parameter (1.0101 ≤ r ≤ 5.0), for an image bundle with radius ∆x = 0.01. The surface is sampled over a Nr × N φ = 100 × 100 polar grid.
We compare our RBM-based estimates of second flexion with the analytic values using the mean-square error, Ms:

Ms
Here, G(p, q) is the grid-sampled analytic surface (one of Gi or |G|), andĜ(p, q) is the grid-sampled, RBM estimate. We define a peak signal-to-noise ratio as: which provides a quantitative value for the equality of surfaces: PS → ∞ if the surfaces are identical, and Max(Ĝ) is the maximum value of the surfaceĜ(p, q). We find minimal dependence in the calculated PS on the gridded surface resolution, with a variation ∆PS = ±0.4 when PS > 20 for grids with dimensions 50 ≤ Nr, N φ ≤ 250. We plot PS versus ∆x for RBM-recovery of G1 in Fig.  3. Results for G2 and |G| are comparable with G1. For the components of second flexion, the peak signal-to-noise continues to grow as image size is reduced, and we approach the theoretical point-source/point-image case. This result demonstrates that the RBM can indeed be used to recover second flexion for the SL model, and gives us confidence that this technique can also be applied for other, more complex lens models in a backwards ray-tracing mode (i.e. when the lens equation cannot be inverted to give image locations as a function of a source location).
Since the SL model has zero first flexion, we are not able to demonstrate that the RBM is actually sensitive to the centroid shifts from a lens model with non-zero first flexion components. What we can show, however, is that RBM is consistent with recovering zero first flexion for the SL model. Using an approach similar to that outlined above for Gi and |G|, we determine estimates for F1 from the matrix solutions as shown in Table 1. The maximum (absolute) values of the recovered F1 values are based on bundle radii in the range 10 −4 ≤ ∆x ≤ 1 and polar grid resolution 50 ≤ Nr, N φ ≤ 500. Results for F2 and |F| are comparable.  A slight increase in the maximum recovered first flexion is seen at the lowest bundle radius: from ∼ 10 −7 for ∆x = 10 −4 to ∼ 10 −9 for ∆x = 10 −3 , suggesting that we have reached a numerical limit (in our implementation). We propose that ∆x = 10 −3 is thus an appropriate lower bundle radius to maximise accuracy in using the RBM for an arbitrary lens configuration for recovering flexion.
Due to our choice of Nray = 5, there is a slight dependence on the orientation of the bundle with respect to the lens. Suppose we create bundles that are evenly spaced in angle, φ, around the lens, so that the bundle centres are: and the individual bundle rays are located at twodimensional coordinates with ψ k = 2kπ/Nray for k = 1 . . . Nray, and ψ0 is a phase term that is fixed per bundle. If we set ψ0 = 0 for all bundles, we get the situation shown in the left-hand panel of Fig. 4, where the image (regular pentagons) and source bundles (irregular pentagons) are shown with respect to the Einstein radius (dotted line). While the bundles retain the same orientation in the image plane, independent of φ, the radial separation between the lens and the closest point of the bundle varies with φ. This means that the tidal field experienced by each bundle depends weakly on φ. For the case of image bundles that appear more circular (right-hand panel, Nray = 100) this effect is much less pronounced, and we have the expected rotational symmetry. The consequence of this is demonstrated in Fig. 5, where we plot |Ĝ(r, φ)| on a 100 × 100 polar grid for ∆x = 0.2 and ψ0 = 0 (top left), ∆x = 0.2 and uniform randomly selected ψ0 (top right), ∆x = 0.01 and uniform randomly selected ψ0 (bottom left) and analytic |G(r, φ)| (bottom right). A sinusoidal variation with φ is apparent when ψ = 0 for fixed r, and a scatter is introduced when ψ0 is chosen uniformly randomly between [0, 2π). The magnitude of the φ dependence is reduced by choosing a smaller ∆x.
To quantify this effect we determine the minimum impact parameter, rτ , at which G(rτ ) −Ĝ(rτ ) /G(r) ≤ τ % for a given image size, and G(rτ ) = Max [G(r, φ)]. In Fig. 6 we plot the maximum (thick line) and minimum (thin line) percentage deviation between |Ĝ(r)| and |G(r)| as a function of the image bundle radius. Samples were evaluated on a 250 × 250 polar grid. The horizontal line is a nominal error of 1%, indicating that a sufficiently high level of accuracy occurs for ∆x = 0.01 (vertical line). Moreover, we find that for ∆x < 0.01, the error in Gi is less than 1% for all impact parameters r > 1.01, which suggests that we can use the RBM approach for impact parameters extending almost to the Einstein radius. From this result, we infer that any variation we see in the RBM-recoveredĜi for an arbitrary lens configuration is evidence of a φ-dependent second flexion, which differs from the predictions of a circularly-symmetric lens models.

Source-to-image plane ray-bundles
For the SL we have the fortunate situation that there are analytic solutions for both θ → β (backwards) and β → θ (forwards) versions of the lens equation. This gives us more flexibility in defining the image shape in order to test the range over which the flexion formalism is appropriate. In the previous section, we started with a regular image, which is deflected to a teardrop source. Here, we use a regular polygonal source, and obtain the corresponding image shape. Propagating this modified (i.e. flexed) image backwards as a raybundle, we hope to recover the input source shape.
We note that an image can cover a small portion of the image plane, and hence the tidal fields can be weak over the extent of the image. For a strongly flexed ("banana- shaped") image, more of the image plane is covered by the image, and hence tidal effects are more prevalent. For an arbitrary lens model or configuration, where we do not have analytic solutions for β → θ, we must work in the "regular polygonal image maps to unknown source shape" regime, so it is critical to understand the limits of our technique. Fig. 7 shows six sample configurations. In each case, a circular source (A) is mapped to a flexed image (B) using the inverse solution, equation (30), of the lens equation for the SL. Next, equation (6) maps the differential image ray coordinates, δθi, back to the source plane to produce a new source profile (C). Finally, equation (5) is used to create an elliptical, shear-only image shape (D) from the circular source. The major and minor axis lengths of image D are comparable to the equivalent (distorted) axes for the flexed image (C). Source C and image D are presented as mirror images of their true positions for clarity in the figure. We have used Nray = 100 for the bundles in this figure in order to show the true shapes, hence we can neglect the φdependence (see Section 3.4). While the Scwharzschild lens model produces two images for each source position, we only consider the more highly magnified image outside of the Einstein radius. Significant variation between the source shapes, A and C, is an indication that we are in a regime where the second order Taylor series expression for flexion is not valid Figure 7. Sample plots showing the limitations of the flexion formalism for extended sources using the Schwarzschild lens model, with respect to the Einstein radius (dotted circles). In each case, a circular source (A) is mapped as a ray-bundle to a flexed image (B) using the lens equation (i.e. reverse RBM). This flexed image is mapped back to the source plane (C) using equation (6). The elliptical, shear-only version of the image (D), obtained by applying equation (5) to the circular source, is shown for comparison. Note that source C and image D are presented as mirror images of their actual positions for clarity in this figure. Analytic expressions for flexion of extended sources are valid in cases where source shapes A and C are comparable; shear-only weak lensing analysis is satisfactory when image shapes B and D are comparable. The level of discrepancy between true and recovered source shapes increases for larger image radii and small impact parameters. Each bundle comprises 100 rays. Source bundle radii are: ∆y = 0.5 (left column), 0.2 (top right and middle right), and 1.0 (bottom right). I and S values are noted for each scenario, as defined in equations (45) and (46).
for extended sources -this is the true strong lensing regime where images are arcs rather than arclets.
We examine the level of agreement between the Taylor expansion, equation (6), and the forwards RBM solution by comparing the relative locations of each image or source ray in the bundle. We introduce two quantities: where the two-dimensional vectors, IB,n and I ′ D,n , are the Nray light rays in the actual (B) and elliptical (D) images, and similarly S = log 10 σ 2 S = log 10   1 Nray where the two-dimensional vectors, SA,n and S ′ C,n , are the light rays in the initial (A) and recovered (C) source bundles.
We calculate I and S for source bundles with radii ∆y = 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1 and 2. In Fig. 8, we plot I and S as functions of the image bundle impact parameter, xc. We select I ≥ −4 and S ≥ −4 as indicative that the distortion is significant; this was confirmed by eye using plots similar to Fig. 7. For higher values of I and S, there were clear differences between sources A and C, and images B and D. We find xc 1 + 2.2∆y, for comparison of source bundle shapes A and C (S ≥ −4) and for comparison of image bundle shapes B and D (I ≥ −4). Next, we consider a more quantitative approach, based on comparing bundle ellipticities, in a manner comparable  (45) and (46), plotted as a function of log 10 (xc − 1). Individual lines are for different source radii (top to bottom) ∆y = 2.0, 1.0, 0.5, 0.2, 0.1, 0.05, 0.02 and 0.01. The horizontal dashed line in each panel is at I = S = −4, taken as the limit above which shape distortion is significant. The dots represent the empirical fits (left) x hi 1 + 3.6∆y for I and (right) x lo 1 + 2.2∆y for S. The noise in the lowest lines is due to the use of single precision floating point numbers.
to the standard analysis for examining weak lensing-induced shear. Defining quadrupole terms relative to the bundle centroid we consider a complex ellipticity (e.g. Schneider 2005) of the form: which has norm: We define error terms: for comparison between ellipticities, and determine the image plane impact parameter, xc, at which a source with radius, ∆y, first exceeds E = 1%, 5% and 10%. Since a circular source has |χ| = 0, we use elliptical sources with axis ratios b/a = 0.8, 0.9 and 0.99. We consider the two cases where the semi-major axis is aligned tangentially or radially to the Einstein radius, computing these limits and also the average (based on original data values) of these two orientations.
Plotting results as log 10 (xc − 1) versus log 10 ∆y, see Figure 9. Minimum impact parameter, xc, at which the error, E AC , is first below 5% as a function of the source radius, ∆y, for elliptical source shapes with axis ratios b/a = 0.8 (red, dotdashed line), 0.9 (blue, dashed line) and 0.99 (black, solid line). Lines are the least-squares fits in log-log space to the functional form xc = 1 + ǫ∆y n . Fig. 9, we see a relationship that is highly suggestive of a functional form: We perform a least-squares fit in log-log space to obtain the parameters ǫ and n. Results of these fits are presented in Table 2 -in all cases, the calculated Pearson coefficient is r > 0.994, indicative that equation (55) is an appropriate functional form. There is variation in the fitted parameters based on the chosen source axis ratio. This is not suprising, as the tidal gravitational field across the resultant image depends on the relative separations and orientation of individual image rays from the lens (c.f. with discussion on ori-  (56) and (57) respectively. The white region is the true strong lensing regime, while the light grey-shaded region represents the preferred zone for both weak lensing shear and flexion analysis. The dark greyshaded region, ∆y 0.5 and hence source sizes comparable to the Einstein radius, means the stronger tidal fields "flex" rather than "stretch" images. Table 2. Least-squares fitting parameters for the functional form xc = 1 + ǫ∆y n based on relative errors in source A and C ellipticities. Fits are made for sources with axis ratios b/a = 0.8, 0.9 and 0.99 with the semi-major axis aligned either tangentially or radially to the Einstein radius. Fits were performed in log-log space, but averaging is performed with original data values. In all cases, the Pearson coefficient is r ≥ 0.994. The inferred ǫ and n values based on the limit S ≥ −4 are shown for comparison. entation of image bundles in section 3.4). As a best estimate, we average over the three chosen axis ratios for EAC = 5%, to obtain the second-to-last row in Table 2. We propose the following interpretation: if we see an image that looks like B, and we use the flexion formalism to determine what the source would look like, we would be wrong (error of 5% in source ellipticity) unless: We refer to the boundary defined by this expression as the start of the "flexion zone" (blue line and both grey-shaded regions in Fig. 10). For image impact parameters closer to the lens than this limit (white region in Fig. 10), the second-order Taylor series approximation given by equation (6) is not sufficiently accurate when applied to an extended source, and we are in the true strong lensing regime. The radius of the flexion zone boundary increases for larger source sizes, ∆y, relative to the Einstein radius, which is expected as there will be greater tidal field variations across an image/source bundle. We perform a similar analysis for the relative ellipticity error between images B (RBM) and D (shear-only), although we now revert to using a circular source only. Parameters are presented in Table 3. We find that for xc 1 + 3.87∆y 1.15 , the lensed image shape and a shear-only intereptation are essentially the same (error 1% in ellipticity) implying the first-order Taylor expansion is adequate. We refer to the boundary defined by this expression as the start of the "shear zone", where image shapes are essentially indistinguishable from ellipses, and a traditional weak-lensing (i.e. shear and convergence only) analysis is satisfactory. However, our result does not preclude use of the flexion formalism at larger impact parameters, as a potentially measurable non-zero flexion remains.
There is a crossover between the boundaries when ∆y ∼ 0.5, indicated by the dashed line in Fig. 10. In some sense, the light grey-shaded region is the preferred region for weak lensing shear and flexion analysis, as it corresponds to both shear-only ellipticity errors < 1% and flexion-recovered source ellipticity errors < 5%. The dark grey-shaded region, ∆y 0.5 and hence source sizes comparable to the Einstein radius, means the stronger tidal fields "flex" rather than "stretch" images. The effect of this is demonstrated more clearly in Fig. 12, when we consider cosmologically-realistic scenarios.

APPLICATION
We approach our application of the results from the preceding section with the understanding that the Schwarzschild lens is not an ideal description of the extended mass distribution of, for example, a galaxy cluster lens. Analytic predictions for first and second flexion do depend on the density profile of the lens model. Indeed, from Fig. 2 in Paper I, it can be seen that the Schwarzschild lens has zero κ and |F|, while these values are non-zero for the extended mass profiles (ie. SIS, NFW and Sérsic profiles). However, Birkhoff's theorem allows us to consider all truncated lens models to be Schwarzschild-like outside of the truncation radius. For comparisons with observations, we need to convert our Einstein-radius scaled results from Section 3 back to angular units on the sky. With θi as the angular position in the lens plane and βi the angular position in the source plane, we have in the small angle limit (which is appropriate for sources and lenses at cosmological distances): where θE, the angular (point mass) Einstein radius, was defined in equation (28).
The desired criteria (lower dashed line) at 0.01 for a 1% flexion error, is only met for low lens redshifts and high mass lenses, typical of galaxy clusters, and for small source galaxy radii. A more conservative limit at 0.1 (upper dashed line), results in a 10% error in recovered flexion values for extended sources -see Fig. 6. While this may appear somewhat discouraging for flexion programs, we note that the required size criteria could be reached by considering isophotes of an image, corresponding to smaller source sizes. We do not discuss this further in the present work.
In Fig. 12, we plot the location of the inner boundaries of the flexion zone (dashed lines) and the shear zone (solid lines), for several typical scenarios: lens masses 10 13 M⊙ (top panel), 10 14 M⊙ (middle panel) and 10 15 M⊙ (bottom panel). In each panel (from top to bottom), the lens redshifts are z d = 0.1 (black), 0.2 (red), 0.5 (green) and 1.0 (blue), the source radius is ∆η = 10 kpc, and z d < zs ≤ 2. Note the crossover between these regions that occurs for M = 10 13 M⊙, and in most cases presented here, the shear zone actually starts closer to the lens than the flexion zone. Fig. 13 shows the effect of changing the source radius, z d = 0.2, with ∆η = 20 kpc (dashed lines) and 5 kpc (solid lines). Lens masses used were 10 13 M⊙ (black), 10 14 M⊙ (red) and 10 15 M⊙ (blue).

CONCLUDING REMARKS
While weak lensing with shear is now well-established both theoretically and observationally in the galaxy-galaxy lensing and cosmic shear cases, weak lensing via flexion is still in its infancy.
In Paper I, we considered analytic solutions for circularly-symmetric lens models. In this paper, we have demonstrated how the ray-bundle method can be used to recover the analytic second flexion results for the Schwarzschild lens model to high accuarcy, and is consistent with the recovery of zero first flexion. Indeed, we recover the Schwarzschild model second flexion solutions with errors no worse than 1% if bundle radii are ∆θ 0.01θE. In such circumstances, the second-order Taylor series expansion used by Bacon et al. (2006) is appropriate for extended sources.
Furthermore, we have identified the existence of a flexion zone in the image plane, which can be considered an optimal region for applying the analytic flexion formalism to extended sources.
The ray-bundle method now provides us with a valuable numerical tool for studying flexion due to complex lens distributions, such as asymmetric lens models or cosmological structures, where no such analytic solutions exist.

ACKNOWLEDGEMENTS
This research was supported under Australian Research Council's Discovery Projects funding scheme (project num-