On the Role of Caustic in Solar Gravitational Lens Imaging

We consider scattering of electromagnetic waves from a distant point source by the gravitational field of the sun, taking the field oblateness due to the quadrupole moment of the sun into account. Effects of the field oblateness can play an important role in the high resolution solar gravitational lens imaging in the sub-micrometer wavelength range of the electromagnetic spectrum.


Introduction
The idea of using the sun as a powerful telescope goes back to Eshleman [1]: The gravitational field of the sun acts as a spherical lens to magnify the intensity of radiation from distant objects along a semi-infinite focal line with the nearest point of observations being about 550AU (for a general introduction, see e.g. [3], [4], [5]). For example, the intensity from a distant point source of electromagnetic (EM) radiation at λ = 1µm wavelength can be pre-magnified by the sun gravitational lens up to µ ≈ 10 11 times. Depending on observation device, the resolvable angle between two point sources at this wavelength could be as small as 10 −10 arcsec.
Recently, properties of the solar gravitational lens attracted attention both due to discovery of numerous exo-(and possibly earth-like) planets and the success of the Voyager-1 spacecraft, presently operating at about 140AU. Possibilities of mega-pixel imaging of such planets from the focal line of solar gravitational lens are now being discussed.
In the present work we consider effects of oblateness of the gravitational field and that of rotation of the sun on the image formation and the diffraction pattern of the lens. Although the quadrupole moment of the sun is very small, effects of oblateness, nevertheless, turn out to be important: The focal line caustic unfolds and can have several hundred meters in cross section at the distances up to several thousand AU from the sun. Moreover, for the wavelengths in a micrometer range, the diffraction pattern of the point monochromatic source changes significantly and the maximum of the amplification of the EM energy flux radiated by such a source can decrease up to several orders of magnitude, depending on direction of observation and the distance between the sun and an observer.
We stress that the last statement concerns the maximal magnification of density of flux radiated by a point source, and not the intensity magnification for a realistic extended object. Since the "focal blur" of the gravitational lens is comparable with the size of the whole image 1 , there will be no deflection of initially parallel light rays coming from an infinitely distant point source.
A trajectory of light in the gravitational field of the sun can be found using post-Newtonian approximations for the null-geodesics of the post-Minkowskii metric element (see e.g. [12], [13], [16], [14], [15]) where Φ is the scalar (Newtonian) gravitational potential and A is the gravitomagnetic vector potential.
In the asymptotically Cartesian heliocentric coordinate system where parallel beams are incoming from z = −∞, the post-Newtonian deflection angle α, which is the difference between the incoming and outgoing beam direction vectors, equals the two-dimensional gradient of the potential Ψ where Z is the z-coordinate of an observer. The condition Z √ x 2 + y 2 r g , where r g ≈ 3 × 10 3 m is the gravitational radius of the sun, is also imposed on eq. (1). In this limit the two dimensional gradient of Ψ is independent of Z: The potential Ψ is a sum of the x, y-dependent (Z-independent) and the Z-dependent (x, y-independent) terms (see below), so it is essentially a two-dimensional potential.
Here, one can apply the thin lens approximation which leads to the following picture (see Figure  1): a light ray is incoming from z = −∞ and hitting the z = 0 "lens plane" at (x, y). At this plane the ray is deflected by the angle given by the Z → ∞ limit in eq. (1). Then, it follows that the outgoing ray intersects the observer plane z = Z > 0 at the point whose position (X, Y ) is determined by the extremum of the Fermat potential S ("lens equation"): where This equation is a manifestation of the Fermat principle for the beam delay time S/c. For a compact lens, a combined contribution to the two-dimensional potential Ψ in eq. (1) from the dipole terms of Φ and A can be cancelled by a translation. Since the gravitomagnetic field of the sun, produced by its rotation, is a dipole field, without loss of generality we can set A = 0 (for more details see eg [16]).
The exterior Newtonian potential of the sun can be approximated by that of the quadrupole where n is a unit vector in the direction of the polar axis of the sun, R 0 ≈ 7 × 10 8 m is the sun radius and I 2 ≈ 2 × 10 −7 is its dimensionless quadrupole moment . Without loss of generality we select the coordinate system where n = (0, sin β, cos β) with β being the angle between the polar axis of the sun and the incoming from z = −∞ beams (see Figure 1). Figure 1: Diagram of the geometrical optics problem. Section of the caustic surface by the observer z = Z plane (X, Y -plane) is schematically shown on the right, while the corresponding critical line in the lens z = 0 plane (x, y-plane) is schematically shown on the left (see eqs. (13,14)).
The term cT (Z) = r g log(−4Zz source /r 2 g ) can be dropped without loss of generality: It does not affect the geometrical optical values 4 since it is independent of x, y.
In order not to carry numerous constants through the computations, we re-scale both the lens plane and the observer plane lengths with the scaling length parameter 5 The new dimensionless polar coordinates (ρ, φ) in the lens plane and the dimensionless Cartesian coordinates (ξ, η) in the observer plane then read In these coordinates Ψ = 2r g ψ, ψ = log(ρ) − ρ 2 cos 2φ, 3 For details of this simple computation one can also refer to e.g. [15]. 4 The above term can be also dropped in the wave optical computations of the field intensity, which is the square of absolute value of the complex field amplitude, since it contributes only to the common phase factor to the amplitude (see next section). 5 We would like to stress that b is not an impact parameter. The latter will be introduced further.
In the case of the sun ≤ about 10 −7 . The small parameter is maximal when β = π/2, i.e. when the source is placed in the sun equatorial plane. It decreases as the source is displaced towards the sun polar axis, on which it vanishes (at β = 0): The light radiated by a source from the polar axis is deflected as if the sun were spherically symmetric. So, we will refer to the both situations of I 2 = 0 and β = 0 as the "spherically symmetric", "degenerate" or the "monopole" case. The parameter also decreases when the observer plane moves away from the sun (i.e. as Z increases). It is worth mentioning that one can also account for the light refraction in the solar plasma by adding corresponding correction term to the potential ψ. However, this contribution can be discarded for the sub-micrometer diapason of wavelengths. Evaluation of this contribution is given in the Appendix 2.
It follows from (2,3,8,9) that coordinates (ρ, φ) of the images of the point (ξ, η) are solutions of the lens equation, which has the following form in the complex notations We recall that the solution of this equation gives the "impact parameter" bρ(ξ, η; ) and the corresponding polar angle φ(ξ, η; ) in the lens plane for the ray(s) arriving to the observer plane at X = bξ, Y = bη. Both scaling factor b and a small parameter depend on distance Z between the planes. In the geometrical optics, the inverse intensity magnification equals the ratio of the corresponding surface elements on the observer and the lens planes (the Jacobian of transformation (11) or the Hessian of the Fermat potential S) The magnification µ diverges at the critical line (see Figure 1) or, according to (11), at the astroid (tetracuspid) caustic 6 in the observer plane Here we recall the standard procedure of solving the lens equation (11). When the deviation of the observer from the z-axis is much smaller than b, i.e. when |ξ| 1 and |η| 1, we have ρ = 1 + δ, δ 1.
The solutions φ of (15) are the "Einstein ring" coordinates of images of the point (ξ, η). It is not difficult to see that those are angles between the ξ-axis and the tangents to the astroid (14) drawn from the point (ξ, η) (see Figure 2). When the point (ξ, η) is inside the astroid, equation (15) has four solutions. images of a small (centered at the z-axis) disc source viewed by an observer at the z-axis ("Einstein Cross"). Fig.D: Spans of the strong and weak images of the same small disc source viewed by an observer near the cusp: Maximal angular span of the "strong limb" is approximately proportional to the cubic root of the ratio between the size of the heliocentric projection of the source to the observer plane and the size of the astroid (while the span of the "weak limb", as well as limbs of Fig.C, is linearly proportional to the above ratio). Note that image of the disc source does not form a full ring if the apparent size of the disc is smaller than d astroid /2Z.
Otherwise it has two solutions. The magnification of the jth image µ j is inversely proportional to the distance from the observer to the corresponding tangency point on astroid (see also Appendix 1).
In the = 0 case the lens equation has two solutions when ξ 2 + η 2 = 0 or infinite number of solutions forming a unit circle when ξ = η = 0. Since the number of solutions does not exceed four in the non-degenerate case, the image of a source of small (wrt astroid) apparent size never forms the ring 7 if = 0 (see Figure 2).
Returning to the = 0 case we note that the caustic (14) degenerates to the focal line ξ = 0, η = 0 and the deflected beams converge towards the z-axis at the "Einstein" angles 8 The "on-axis" observer "sees" the whole critical line. In other words, the rays are coming from the 6 According to the recent data the dimensionless octopole moment of the sun ∼ 10 −9 . Since is at most ≈ 10 −7 , the caustic cross-section can be considered as a pure astroid for our purposes. 7 Obviously, this does not imply that images of bigger sources, such as e.g. exo-planets of interest, do not form rings. 8 For a beam grazing the edge of the sun α E ≈ 0.85 · 10 −5 rad ≈ 1.75 arcsec.
circle of radius b in the lens plane towards an observer at X = Y = 0 in the z = Z plane. Therefore, b > R 0 , where R 0 is the sun radius. The distance from the sun to the closest focal point is determined by the condition b = R 0 and equals Z min = R 2 0 /(2r g ) ≈ 550 AU. In the wave optics, the maximal spatial resolution of the spherical lens in the neighborhood of the focal line is restricted by the radius of diffraction (radius of the Airy disc), which is of order λ/α E , where λ is the light wavelength [8], [9], [10], [11]: The circularly symmetric diffraction pattern of a point source oscillates in the radial direction and its intensity reaches maximum at the z-axis. The spatial scale of the oscillations is of order of the diffraction radius. For λ = 1µm, at the position of the closest observation Z ≈ 550AU, this radius is about decimeters.
On the other hand, from eqs. (7,8,10,14) it follows that the non-spherical model produces an astroid caustic of the diameter 9 which reaches up to d max ≈ 5.6 × 10 2 meters when β → π/2 and Z → Z min . Thus, effects of oblateness clearly lead to significant changes of the diffraction pattern of the point source when it moves from the sun polar axis to the equatorial plane. Indeed, the maximum of magnification is now reached in a neighborhood of the astroid, where the geometrical optics magnification diverges, thousands of diffraction radiuses away from the z-axis, which is now nonsingular. We note that, as follows from (17), the size of the astroid is proportional to Z −1/2 , i.e. the size varies slowly with the distance from the sun. For example, the maximal astroid diameter is about 400 meters at 1000AU, while sizes of the heliocentric projection of possible objects of observation are about several kilometers across.
Before going to the detailed wave optics computations in the next section, it is worthwhile to mention some heuristic arguments explaining significant difference in the maximal EM energy flux amplification in the β = 0 and β = π/2 cases for point source and small wavelengths: In the = 0 case (spherical lens) the lens plane image produced by a small distant source consists either of one single ring or two opposite arc-shaped "limbs" (one limb inside and another outside the critical curve ρ = ρ c = 1). Whether the image is a ring or limbs, as well as the size of limbs, depends on the source size and the observer position X, Y in the z = Z plane (the further away the observer is from the z-axis, the smaller are the limbs).
In difference from the symmetric case, in the case when the lens caustic is the astroid (14) the image of a small source 10 never forms the whole ring (see Figure 2). Such an image consists of two to four disjoint small limbs, some of them being weak and some strong, depending on the observer position wrt caustic. When the source size goes to zero, the limbs become a set of (two to four) points. Considering the principal Fresnel zones 11 around these points one can explain decrease of maximal magnification in the wave optics.
In more details (see Figure 3): the above zones have form of "limbs", whose dimensions depend on the size of caustic ∼ b, on the diffraction radius ∼ λ/α E and on a position of the observer: It is easy to see that the thickness (i.e. radial dimension) of the limb is approximately the same in symmetric and non-symmetric cases when is small. arrows. An observer is denoted by the cross. The critical line and caustic are schematically drawn by dashed lines. In the symmetric case (left) the maximal zone corresponds to the on-axis position of an observer. The zone is an annulus that contracts towards the circle as λ → 0. In the = 0, λ r g case (right) the maximal zones correspond to positions near cusps. In this case a zones contract both radially and tangentially towards point images as λ → 0. The angular span of the maximal zone in the = 0, λ r g case is proportional to (λ/ r g ) 1/4 , since e.g. near the cusp corresponding to φ = 0: φ = δφ, ρ = ρ c (0) + δρ and the variation of the Fermat potential for an observer at the cusp equals δS ≈ 2r g δρ 2 + r g (3δρδφ 2 − 3δρ 2 δφ 2 + δφ 4 /2). The radial dimension of zones is proportional to (λZ) 1/2 . However, in the spherically symmetric case even a point source can produce a circular geometrical optics image, and the maximal zone spans the whole circle even if λ → 0 (thickness of maximal zone goes to zero while its angular span always equals 2π). In contrast to the symmetric case, in the = 0 case the geometrical optics image of a point is always a point and therefore the maximal principal Fresnel zone contracts to a point as λ → 0. This leads to decrease of the maximal EM energy flux amplification in comparison with the symmetric case 12 (since the thicknesses of zones have practically the same dependence on λ when is small).

Diffraction Optics and Gravitational PSF
The geometrical theory of diffraction provides algorithms for finding near-caustic intensity from its geometric optics asymptotics [13], [17]. For caustics of the type (14) the near-field is expressed through the Airy function for caustic folds (see e.g. [13], [12], [19]), and through the Pearcey integral near 12 An estimate of the maximal magnification, based on evaluation of dimensions of the Fresnel zone corresponding to the cusp of the astroid (14) was presented in [1] using an analogy with focusing by atmosphere of an oblate planet [2]. There, the angular span ("horizontal dimension" in terminology of [2]) of the maximal Fresnel zone was estimated by using the following assumption (Section 6 of [2]): The variation, along the half-span of the zone, of distance between the cusp and the involute of the astroid equals the "vertical" dimension (i.e. thickness) of the zone. From this it follows that the angular span of the zone is proportional to (λ/ 2 r g ) 1/8 as λ → 0, while from direct evaluation of the variation of the Fermat potential S it follows that the span ∼ (λ/ r g ) 1/4 . caustic cusps (see e.g. [17], [20]). In cases of focal lines the near-field is expressed through Bessel functions [17], [19].
Application of these algorithms to our problem is presented in the Appendix 1. Below, we derive the gravitational point spread function as a single one-dimensional integral having all the above mentioned limits.
The change of polarisation angles of light in weak gravitomagnetic fields is of a post-post Newtonian order and can be neglected (see eg [14], [18] and references therein). The deflection angles are also small and space is asymptotically flat, so we can apply the scalar Huygens-Fresnel principle (Fresnel-Kirchhoff diffraction formula): In the thin lens, short wavelength approximation (Z R 0 r g , λ r g , Z min /r g ≈ 3 × 10 10 ) the diffraction field in vicinity of the z axis is a sum of contributions by spherical waves propagating from the lens plane z = 0 with the phase delays corresponding to the sum of the gravitational and geometric delays. This total time delay equals S/c, where S is given by (3). Then, the complex amplitude of the electromagnetic (EM) field at the observer position equals (up to the phase factor e ikZ ) The intensity magnification, i.e. amplification of the EM energy flux at the point X, Y of the observer plane, equals the square of the absolute value of the amplitude 13 This function of X, Y (µ also depends on Z, and k) is a point spread function (PSF) of a gravitational lens only, so we call it gravitational PSF or GPSF (the PSF of a combination of the gravitational lens and a telescope is discussed in Section 5). For a detailed derivation and justification of the Fresnel-Kirchhoff diffraction integral (18) see e.g. [12] or [19]. In (18), the geometrical optics magnification (12) is recovered as lim k→∞ |u| 2 . When kr g 1 and the point on the observer X, Y -plane is far away from caustic, the integral (18) can be expressed in a simple manner through the geometrical-optics data as (see e.g. [19], [12]) Here the sum is taken over the number of images in the lens plane, µ j is the geometrical optics magnification of the jth image, S j = S(X, Y ; x j , y j ) is the extremal value of the Fermat potential (3) for the jth image and n j = 0, 1, 2 corresponds to x j , y j being the minimum, saddle and the maximum point respectively. The above approximation breaks down in the neighborhood of caustic, the case we are mainly interested in. So, one has to either evaluate (18) exactly or to apply suitable asymptotic methods. Up to a common (X, Y, Z-dependent) phase factor where q is the dimensionless wavenumber For λ = 10 −6 m, q ≈ 3.7 × 10 10 .
The main purpose of this section is the direct numerical evaluation of the 2d integral (20). Before presenting the numerical results we would like to make several remarks: Remark 1: It is worthy to note that in the = 0 case, the exact 2d integration in (20) is possible: where J 0 is the zero-order Bessel function. After integration in ρ one can express µ in terms of the confluent hypergeometric function (see e.g. [19]) In the short-wavelength limit q 1 and when the argument iq(ξ 2 + η 2 )/2 of 1 F 1 is small, i.e.
the hypergeometric function 1 F 1 degenerates to the zero-order Bessel function (see e.g. [19], [11]) The maximum µ = µ 0 of the GPSF is reached at the focal line ξ = η = 0 and equals Remark 2: Condition (23) is, in fact, a condition of validity of the stationary phase integration at ρ = 1 in the last integral in (22). Indeed, the stationary phase approximation can be applied in (22) when the width of the stationary phase region δρ ∼ 1/ √ q is much more smaller than the scale of oscillations of the Bessel function δρ ∼ 1/(q √ ξ 2 + η 2 ), which leads to (23). Remark 3: Condition (23) will be encountered in the next section, when a similar type of the stationary phase integration will be performed for the general case = 0: As follows from (20), in the general situation where the function of three variables F (·, ·, ·) is defined as follows 14 Similarly to the = 0 case, the stationary phase integration can be performed at ρ = 1 when all the arguments of F in (25) are much smaller than √ q, i.e. when 14 F (χ, κ, ν) degenerates to J 0 in two special cases: 1) χ = 0 (see Remark 1) and 2) κ = ν = 0 (see eq. (54)).
and (23) holds. This will be demonstrated in detail in the next section.
Consider now the general case, i.e. the one when (23) is not necessarily true. Below we will proceed with the main subject of the present section, estimating the 2d integral (20) numerically without any assumptions.
Since q 1 (e.g. q ≈ 3.7 × 10 10 for λ = 10 −6 m) one can reduce (20) to a one-dimensional integral in φ by the stationary-phase integration 15 in ρ at fixed φ: The stationary phase integration in ρ produces a relative O(1/q) error, which is negligible for wavelengths of interest.
First, one should find the "stationary phase line" ρ = ρ st (ξ, η; φ) a such that 16 Then, up to a common phase factor From (27) we obtain 17 Therefore the GPSF (intensity magnification) expresses through the one-dimensional integral in φ: where two functions of φ: τ = τ (ξ, η; φ) and ρ st = ρ st (ξ, η; φ) are given in (29). The above integral can be particularly easy taken when = 0, ξ = η = 0, giving the value of µ 0 , obtained earlier (24). In the general case, the integral (31) should be taken numerically. Now, we present results of the direct numerical computation of the GPSF (30,31): Figure 4 shows a ratio of the maximal intensity µ at = 0 to that of the symmetric case µ 0 = πq for different values of the parameter q . The computations presented at Figure 4 are performed at fixed q ≈ 3.7 × 10 10 (which corresponds to λ ≈ 1µm) for It is not difficult to see that solutions of the lens equation (11) lie on this line. When the observer is far away from the caustic, the stationary phase integration in φ can be also performed around these points together with the above integration in ρ. Such a double stationary phase integration results in (19). This is not the case when the observer is in a neighborhood of the caustic since the second tangential derivative of V vanishes on the critical line and one has to apply other methods for computing the integral in φ. 17 Eq.(27) has two solutions: ρ 1 (φ) = ρ st (φ) and ρ 2 (φ) = −ρ st (φ + π). A positive solution has to be chosen, since the integration in ρ in (20) is performed for ρ > 0. Note, however, that permutation of the solutions only reverses the sign of integral (28) since the both solutions parametrize differently the same curve. The former image (q = 10) is 3× zoomed in wrt to the latter. Right Figure: The log − log plot of µ /µ 0 as a function of q. Numerical evaluation is shown by the point plot. High q asymptotics (32) is shown by the dashed line. Low q approximation (see eq.(54)) is shown by gray solid line. The point corresponding to observations in the equatorial plane of the sun at Z = 550 AU and λ = 10 −6 m is marked with a cross. At this point q ≈ 3.7 × 10 3 and µ /µ 0 ≈ 4 × 10 −3 . Computations are performed for q = 3.7 × 10 10 that corresponds to λ ≈ 10 −6 m. Small parameter varies from ≈ 2.7 × 10 −12 to ≈ 2.0 × 10 −7 . Error of numerical integration δµ /µ does not exceed ≈ 1 percent. different values of . However, the numerical results (as well as the analysis in the next section) show that µ /µ 0 approximately depends only on q when (26) holds.
As seen from Figure 4 one naturally recovers the unit ratio for q → 0 (More precisely µ /µ 0 ≈ J 2 0 ( q) when q <≈ 1.4. See eq. (54) and the end of this section). On the other hand when q 1. In this case four equal maxima of the intensity magnification are symmetrically situated on the X, Y axes: two at the X-axis and another two at Y -axis in the caustic interior close to four cusps (X, Y ) = (±4b , 0), (X, Y ) = (0, ±4b ). When q is large, the distance between the cusp and the neighboring maximum is proportional to b/ √ q, while value of µ at the cusp is proportional to µ (i.e. The GPSF oscillates and the amplitude of oscillations grows as one moves towards the neighborhood of the caustic folds/cusps. The amplitude falls off similarly to the = 0 case and the pattern becomes more and more radially symmetric, as one moves far away from caustic in the exterior direction (which is in agreement with (19)).
As q becomes smaller the four global maxima approach (discontinuously, see below) the center X = Y = 0. The diffraction pattern becomes circularly symmetric at q = 0.
It is useful to introduce the parameter χ = q: As follows from (32) In the case of the sun λ 0 ≈ 3.7 × 10 −3 m, i.e. the effects due to quadrupole moments of the sun can be noticeable already at the far-infrared part of the EM spectrum. For λ = 10 −6 m, the parameter χ can be as large as ≈ 3.7 × 10 3 . For these wavelengths the maximum magnification of the energy flux from the point source can decrease up to several orders of magnitude when the source goes from the polar axis of the sun β = 0 towards its equatorial plane β = π/2. Concluding this section we would like to mention that (as can be seen from Figure (4)) µ /µ 0 is a non-smooth function of χ: its derivative wrt χ jumps at certain points (e.g. the first jump occurs at χ ≈ 1.4, the second at χ ≈ 2.5 etc). A jump takes place every time when some of the local maxima of µ(X, Y ) become the global ones: For example, the extremum at the center 54)) is the global maximum when χ <≈ 1.4. At χ ≈ 1.4 this maximum becomes smaller than four maxima at the distance ≈ 0.54λ/α E from the center. As χ increases, positions of new global maxima change continuously until the next jump at χ ≈ 2.5 etc.
In the next section we give an analytic explanation of the numerical results obtained.

Diffraction Optics: Limiting Cases
For analytic description of the above numerical result we expand the argument of exponential in (31) in τ and (we recall that τ = X b cos φ + Y b sin φ, see eq. (29)) where, to the second order in τ and and to the first order in τ , We are interested in the short wavelength limit q 1 of (34). The U 2 -term in (34, 35) can be neglected if q |U 2 | π. This condition holds when The above conditions have been already encountered in (23), (26).
Since |τ | 2/ √ q and q 1, the τ -term in (36) can be dropped. Therefore, provided (37) holds, up to the constant phase factor Rewritten in terms of X and Y , the condition imposed on τ in (37) where R is the distance between the observer and the z-axis. For λ ∼ 10 −6 m (i.e. for q ∼ 10 10 − 10 11 ), the radius R v is about 10 kilometers, while the maximum possible radius of the caustic is about 300 meters, so the above condition of validity of (38) clearly holds in the region of the interest 18 . Therefore, to get the amplitude of EM field one can integrate over the circle ρ = 1 in the lens plane provided (39) holds. This happens due to the fact that the "optical path" S is extremal 19 on the ray trajectories (see eq. (2)) and small deformations of the integration contour do not significantly change contribution from the "monopole part" S 0 of S = S 0 + S 1 when (39) holds. Since the contour deformations are of order of , the error in the quadrupole contribution S 1 is of order of 2 , which is also negligible. In other words, the width of the stationary phase integration region ("thickness of Fresnel zone", see Figure 3) significantly exceeds the deviation of the integration contour ρ = ρ st (ξ, η, φ) from the unit circle when (39) is true.
Apart from the situation when (38, 39) overlaps with the approximation (19), analytical study of (38) can be performed for the three asymptotic cases: (1) The "degenerate" case χ = q 1, i.e the case of observations in directions that are close to the sun polar axis |β| λ λ 0 Z Z min (see eq.(33)). At Z ∼ 1000AU and λ ∼ 10 −6 m this corresponds to β's that are smaller than a fraction of a degree. These directions cover less than 0.01% of the celestial sphere.
(2) The "strongly non-degenerate" case χ = q 1, or equivalently λ/r g . In this case the scale of the diffraction pattern is much more smaller than the transverse caustic size, which takes place for |β| λ λ 0 Z Z min . At Z ∼ 1000AU and λ ∼ 10 −6 m, this corresponds to the directions with β's bigger than few degrees. 18 The ratio of the radius R v to the caustic radius 4 b equals 1/(2 √ q), which leads to the condition 1/ √ q. Since for the sun is at most about 10 −7 , this condition clearly holds in our case. 19 The second tangential derivative of S also vanishes on the critical line.

2.
Let us now pass to the strongly non-degenerate case q 1. In this case the size of caustic ∼ b greatly exceeds the diffraction radius ∼ λ/α E . First we consider the asymptotic of the GPSF in the cusp neighborhoods, where it reaches the maxima. For convenience we choose the cusp at X = 4 b, Y = 0, (i.e. at ξ = 4 , η = 0).
Since q 1, theṼ 1 term in (42) can be neglected when Then, provided the above conditions hold, from (41, 42) it follows that where Pe(x, y) is the Pearcey integral It follows from (30) that, in terms of unscaled deviationX = X − 4 b,Ỹ = Y from the cusp, the near-cusp GPSF equals The domain of validity (43) of the above asymptotics rewrites in terms of the unscaled deviations from the cusp as follows We recall that the above conditions are in agreement with (39), since, as has been mentioned before, the "radius of validity" of (38) greatly exceeds maximal possible size of caustic when 1/ √ q. The absolute value of the Pearcey integral Pe(ξ,η) reaches maximum atξ ≈ −2.02,η = 0 which is inside the domain of validity (43) of eq. (45). The maximum of GPSF in the q 1 limit then equals Therefore the ratio of the maximum µ to that of the spherically symmetric case µ 0 equals The distance between point of maximum and the cusp equals ≈ 2.02 b/ q/2. This confirms the numerical results of the Section 3.
In difference from the symmetric case, where the circular invariant diffraction pattern has radial oscillations, the near-cusp pattern in the q 1 case has a complicated two-dimensional lattice-like structure (see Figures 4,5). The latter transforms towards locally one dimensional structure of smaller intensity (52) as one moves along the caustic away from the cusp.
We now evaluate the GPSF in the vicinity of caustic folds (regular points of caustic) far from the cusps 20 .
This confirms an obvious fact that in a neighborhood of folds the intensity is smaller than the maximum µ in the cusp region. According to (45, 52) the GPSF oscillates and amplitude of oscillations falls as one moves away from the caustic and the "far field" GPSF can be approximated by (19).
3. Concluding this section we mention the case when the observer is on the z-axis, i.e. when X = Y = 0: Here, the integral (38) is expressed in terms of the zero-order Bessel function and µ axis = πqJ 2 0 ( q), which is in agreement with numerical results (see Figure 4). When q 1, i.e. when the size of the caustic exceeds greatly the diffraction radius 21 which is in agreement with (19).

PSF in Focal Plane of Telescope
The GPSF (30, 31) is, in fact, a point spread function for a "zero-aperture" telescope that can be used only for the intensity scan in the observer plane. One should use the PSF of a compound system of the gravitational lens and a telescope (for more details see e.g. [21]) if the diffraction resolution of a telescope is finer than the angular radius of the Einstein ring α E . In the Fraunhofer approximation for the telescope lens of the focal length F, the PSF in the telescope focal plane is expressed through the Fourier transform of the complex field amplitude at aperture where γ and ω are the observation and the point source angles correspondingly (see Figure 7), and R = (X, Y ) are coordinates in the aperture (observer) plane. The complex amplitude of the EM field at the aperture plane u( R − ωZ)e ik ω R is expressed through u( R) given by (18). This expression is obtained by application of the small rotation 22 to the spatial arguments of the complex field amplitude e ikZ u( R). We recall that the phase factor e ikZ was dropped in (18), so we restored it to get the complex amplitude at aperture. After getting the amplitude by application of the above rotation, this factor can be dropped again in (55). For simplicity, we will perform our computations modulo common phase factors. After substituting (18) to (55) one gets (up to a common phase factor) where r ⊥ = (x, y) and A is the aperture function, i.e.
The last term ik R 2 2Z in the last exponential of eq. (56) can be dropped when a √ λZ, where a is the radius of aperture 23 . Therefore where function A k is proportional to the Fourier transform of the aperture function For a circular aperture of radius a, A k can be expressed through the Bessel function of the first order As in the case of the GPSF 24 , the stationary phase integration in ρ (recall that r ⊥ = bρ, see (8)) can be performed in (57). Indeed, a variation of the argument of the Bessel function in (57, 58) across the stationary phase region is ∼ kaδr ⊥ /Z, where δr ⊥ ∼ √ λZ is the width of this region. Therefore, such a variation can be neglected when a √ λZ. The last condition holds for any realistic aperture.
23 For Z ∼ 1000AU and λ ∼ 10 −6 m, √ λZ ∼ 10 4 m, so this term can be dropped for any realistic aperture 24 As in the case of the GPSF, we drop the cT (Z)-term in (9), since the variation of the optical path over the aperture due to this term is at most 2r g aω/Z ≈ aα 2 E ω. Indeed, even for ω ∼ α E this variation is ∼ a × 10 −15 , not to mention ω ∼ apparent sizes of object to be observed.
When condition (39) holds, one can integrate over the unit circle ρ = 1 in the observer plane and where function of one variable h(·) is defined as and Γ, θ are the dimensionless polar coordinates in the focal plane γ = (α E Γ cos θ, α E Γ sin θ).
In eq. (59), X, Y denote the deviation of the observer from the heliocentric projection of the point source ((X, Y ) equals − R on Figure 7): If an aperture is much more smaller than the diffraction radius λ/α E , the PSF (59) is independent of θ and it is proportional to the GPSF (the focal plane image is rather the Airy spot than the limbs for such apertures). On the other hand, when an aperture is big enough, so that the telescope diffraction resolution is much finer than α E and when the object apparent size ∆ω max is much more smaller than the diffraction resolution of telescope (i.e. for aα E λ and a|∆ω max | λ) this function is "concentrated" (within the telescope diffraction limit) in the neighborhood of the "Einstein circle" Γ = 1 (as on e.g. Figures 8, 11).
Note that the ratio M/M 0 , where M 0 = max γ M| R=0, =0 and R = (X, Y ), is essentially a function of six arguments Thus, our results are scalable. For example, max γ, R M/M 0 for λ = 10 −6 m, a = 1m and Z = 550AU is the same as for λ = 0.5 × 10 −6 m, a ≈ 0.71m and Z = 1100AU. It is also important to note that, in difference from the GPSF µ, the maximum of the focal plane PSF M is not necessarily smaller for equatorial observations in comparison with polar observations, when the aperture is big enough (see Figure 8).
In more detail: Consider first the = 0 case. (59) shows (see Figure 9) that the focal plane PSF, as a function of γ and R (at fixed a, λ, Z and F), reaches its global maximum

1) The numerical integration in
when the telescope is placed at some non-zero distance R = R m from the caustic line. The maximum is reached at γ corresponding to positions of the geometrical optics images. Let us now find M m and R m analytically.
Without loss of generality we can set X = R ≥ 0, Y = 0, so that one of the geometrical optics images 25 would be at Γ = 1, θ = 0. Considering the case when both the aperture and R are much 25 When = 0, R √ λZ, a √ λZ and λ r g , the PSF is symmetric wrt central inversion θ → θ + π . bigger than the diffraction radius (i.e. a λ/α E and R λ/α E ) we can approximate value of F in eq. (59) at Γ = 1, θ = 0 as follows where We would like to compare the intensity at R = 0, θ = 0 with maximal intensity of the circularly symmetric ring seen by the on-axis observer. From (59) it follows that for such a ring (up to a common, θ-dependent phase factor) Therefore (see (59)), the maximal intensity of the circularly symmetric ring seen by the on-axis observer equals q π a Fα E 2 , a λ/α E . Figure 9: Ratio of the maximal focal plane intensity for R = 0 to that for R = 0 in the = 0 case is shown by the gray point plot. It is obtained by numerical integration in (59) for aα E /λ ≈ 10. The ratio (vertical axis) is plotted against Q = ka 2 α E /R (horizontal axis). The ratio given by eq.(61) is plotted by the black line. Note that, at fixed R, the number of the focal plane points where maximum is reached is either two or four. The two point maxima (i.e. maxima at θ = 0 and θ = π) correspond to values of Q at which the above two plots (approximately) coincide. The maximum max γ, R M| =0 = M m ≈ 1.35M 0 is a two-point maximum. It is seen when Q = Q m ≈ 5.7 (see Figure 8B).
Note that for small apertures F | =0,Γ=1,R=0 ≈ 1, and therefore in the asymptotic cases we have For an intermediate range of apertures one has to evaluate M 0 numerically. We will take M 0 as a reference value in what follows. Integral in eq. (60) can be expressed in terms of the Bessel functions 26 and we get This ratio can exceed unity and takes the biggest value ≈ 1.35 at Q = Q m ≈ 5.7 (see Figure 9 and Figure 8B). In other words, the maximum max γ,R M = M m ≈ 1.35M 0 is seen when the telescope is placed at the distance R = R m from the caustic line, where Note that once an aperture is much bigger than the diffraction radius, the similar condition for the distance R m λ/α E is satisfied automatically. The domain of validity of eq.(62) is also restricted by the condition (39), i.e. R m √ λZ or λ (a 4 α 2 E /Z) 1/3 = (2r g a 4 /Z 2 ) 1/3 . For a ∼ 1m and Z ∼ 1000AU this restriction reads as λ 10 −8 m. The fact that max γ, R M| =0 is not seen when the telescope is centered exactly at the caustic line is explained as follows: Although the telescope at R = R m takes less total energy flux, the size of the bright limb/spot in the focal plane is also smaller (compared to the full ring for telescope at R = 0), so the maximal flux density through the focal plane happens to be a bit bigger at R = R m than at R = 0. At big distances from the caustic line, the characteristic size of the focal plane images is defined only by the diffraction limit of the telescope and the ratio max γ M( = 0; R)/M 0 decays proportionally to 1/R (i.e. proportionally to Q) when R R m (i.e. Q Q m ) (see Figure 9). Now, let us return to the non-degenerate = 0 case.

2)
Here, analysis of behaviour of the PSF is more involved, so we present only some general comments and preliminary (mostly numerical) results.
Behaviour of the PSF can be easily described in both small and big aperture limits. When the aperture is much smaller than the scale of the diffraction pattern (i.e. when α E is much more smaller than diffraction limit of telescope) the focal plane PSF is essentially a product of the GPSF and the PSF of the telescope lens. Therefore On the other hand when aperture is much more bigger than the diameter of astroid the maxima of the PSF should be approximately the same in the degenerate and non-degenerate cases. Numerical results as well as analysis below show that the maxima can become approximately the same already at apertures that are much more smaller than d astroid (see Figure 10). To get an analytic estimate for the corresponding range of apertures and absolute maximum of M/M 0 (or that of M/M m ) we apply the approach similar to that of the symmetric case. In more detail: Numerical computations show that when the aperture is big enough, max γ, R M is reached at four symmetrically situated points (two at the X-axis and two at the Y -axis) and at γ corresponding to a position of a brightest (at given X or Y ) geometrical optical image. Therefore, we place an observer at the X-axis at the point X = d astroid /2 +R, Y = 0, so that the distance from the cusp is much smaller than the diameter of astroid |R| d astroid . Value of M at the point Γ = 1, θ = 0 of the focal plane (i.e. at the position of the brightest geometrical optical image) can be obtained in a manner similar to that of the symmetric case (see eq. (60)). Then, up to a common phase factor, we get Similarly to the degenerate case, the maximum is observed whenR ≈ ±R m , i.e. at X ≈ d astroid /2 ± R m , Y = 0, which is confirmed by the numerical computations. Note that, due to approximation where . Note the oscillatory behavior of curves in an interior of the astroid (i.e. for 2X/d astroid < 1) due to the diffraction pattern (see Figure 5).
the ϕ 4 term is neglected in (63), we got eight equal maxima (four inside and four outside the astroid) instead of four: Although for a a 0 the difference between the maxima inside and outside the astroid is negligible (see e.g. curve D of Fig 10), only four or them are the global ones.
In the intermediate range of apertures a ∼ a 0 , a λ/α E the analysis of the PSF becomes non-trivial. However, it follows from (63) that for a λ/α E , the ratio maxR M| Γ=1,θ=0 /M 0 is approximately a function of a/a 0 only. Therefore, it is sufficient to perform a set of computations for a fixed q 1 and different p's (i.e. different apertures) to approximately get the absolute maximum of M/M 0 (i.e. that for all values of parameters) in the strongly non-degenerate case.
Numerical computations performed for the strongly non-degenerate case show that, once the aperture exceeds some value ∼ a 0 , the maximum of "non-symmetric" PSF becomes even bigger than that of the degenerate case, i.e. the maximum of M exceeds M m when the aperture is big enough. However, the difference in these maxima cannot exceed few (∼ 10) percents. The results of numerical computations are shown on Figure 10. Examples of the corresponding focal plane images are shown on Figure  8. Note that a 0 ∼ 1m for equatorial observations at Z ∼ 1000AU and λ ∼ 10 −6 m. Therefore, for a telescope of a modest aperture (∼ 1m), the maximum of the PSF (PSF as a function of γ and R) changes only by at most ≈ 10 percents when the source is moved from the equatorial plane to the polar axis.
Finally, let us estimate the PSF when a telescope whose diffraction resolution is much finer than α E is placed far away from the caustic. Here our computations facilitate due to the fact that one can perform the two-dimensional stationary phase integration in (57). Then, similarly to the case of the GPSF (cf. (19)) , one gets where r j = (x j , y j ) are coordinates of the jth image in the lens plane and S j = S( R; r j ), µ j , n j are same as in eq. (19). Computation of the PSF then reduces to the summation where P( ϑ) = |A k ( ϑ)| 2 is the PSF of the telescope lens. When the angular separation between the geometrical optics images of a point source is much more bigger than the diffraction limit of a telescope, the double sum in the last equation can be dropped.
Then the focal plane PSF equals the following convolution In other words, for a telescope of reasonably big aperture (aα E λ) placed away from caustic, the gravitational lensing can be described by the geometrical optics, while the wave effects due to finite aperture have to be taken into account. Related examples of the focal plane images are shown on Figure  11.

Discussion and Conclusions
In this section we would like to discuss implication of effects related to the quadrupole moment of the sun on prospective observations. For this we first summarize the main results of the present work: The transversal size of the astroid caustic (due to the quadrupole moment of the sun) could reach several hundred meters at the distances ranging from that of the closest observation (550 AU) and up to several thousand of AU. This size is comparable with sizes of heliocentric projections of possible objects of observation, which are about several kilometers across. At such distances, the diffraction pattern of a monochromatic point source transforms significantly (in the region of interest, at e.g. sub-micrometer wavelengths) when the direction of observation is changed from the one along the sun polar axis to that in the sun equatorial plane. The maximum of the gravitational point spread function (GPSF) can differ up to about two-three orders of magnitude. In the strongly non-degenerate case the maximum of GPSF is reached in a neighborhood of the cusp of a caustic. The GPSF can be expressed in terms of the Pearcey integral in this neighborhood.
On the other hand, behaviour of the PSF of a compound system of the gravitational lens and a telescope depends on the telescope's aperture. If the aperture is small, the focal plane PSF and GPSF are essentially the same. For big apertures (e.g. 2m aperture) the absolute maximum of the focal plane PSF can be even bigger in the non-symmetric case. Although these maxima do not differ very much, the formation of images can be significantly different in the symmetric and non-symmetric case. For instance, in contrast to the symmetric case, in the strongly non-degenerate case an image of a point source never forms a "bright" ring, but rather consists of small limbs/spots. The focal plane image is not generally centrally symmetric (which can be an advantage 27 ), number of images of a point source can be different etc.
We recall that the GPSF/focal plane PSF are magnifications of a monochromatic point source and not those of a realistic extended object. Therefore, in prospective missions, one does not expect to directly observe diffraction patterns. However, one needs the PSF for de-convolution of realistic images.
In more detail: the energy flux at the point R = (X, Y ) of the observer plane radiated by an extended, totally (spatially and temporally) incoherent source, equals the convolution of GPSF and the surface brightness of source Here I s ( R ; q)dq is a non-magnified energy flux (times the filtering function) in the interval q, q + dq of the spectrum radiated by a surface element on the source, and R = (X , Y ) stands for coordinates of the heliocentric projection of this element to the observer plane (see Figure 7).
The role of caustic in de-convolution process can be already seen in the geometrical optics: Suppose that one got the intensity (65) in the non-symmetric case and then performed de-convolution of I( R) as if µ were magnification of the monopole. The finest spatial resolution of such a de-convolution will be of order of the astroid diameter d astroid , since both non-symmetric and monopole µ( R) have the same asymptotic behavior at R d astroid and they start to differ significantly at R ∼ d astroid . As has been mentioned before, the typical size of the heliocentric projection of an exo-planet to the observer plane is about several kilometers at 1000AU and d astroid is about 10 percents of that size. Therefore, the area of a minimal pixel will be of order of 1 percent of that of the whole de-convoluted image, i.e. changes of µ play an important role starting from about hecto-pixel level of imaging, not to mention the mega-pixel imaging currently discussed in the literature [10], [11]. At the latter level of resolution one has to take the diffraction pattern of PSF into account.
Similar arguments can be applied to the de-convolution of the focal plane images. It is most likely that this type of de-convolution, rather than that of intensity scan (65) will be used in a possible mission where a sequence of images of the Einstein rings along a path of a spacecraft across the "high intensity" region will be taken. The intensity of a ring point in this set, in a sense, encodes "projection data" taken along a "section" of the source surface. Therefore, a kind of "tomographic" reconstruction algorithms should be developed for the de-convolution related with the focal plane PSF.
One might consider a hypothetic possibility of obtaining a relatively high resolution image during a single arbitrary passage in a neighborhood or/and through the heliocentric projection of an exo-planet. Indeed, at Z ∼ 1000 − 2000AU, a modest ∼ 1 meter telescope equipped with coronagraphs resolves the Einstein ring with about ∼ 20 circumferential elements. Therefore, taking about ∼ 10 2 samples of rings along e.g. a straight path crossing a high-intensity region could, in principle, result in a kilo-pixel image. Note that at this level of resolution the size and structure of caustic plays an important role.
Completing this section, we would like to mention the influence of the higher multi-pole moments of the sun on the PSF. Computation of the PSF accounting for higher moments is a straightforward generalisation of the case considered in this work: one should add the higher-order harmonic terms to the potential ψ in (9) where n are complex harmonic moments. Then one could perform the geometric optics analysis and the stationary phase integration with an account of a new potential.
Note, that the corrections to the PSF due to the fluctuations of plasma density in the solar atmosphere might be even more important than corrections accounting for higher harmonic moments. of the caustic surface equals to (see e.g. [13], [17]) where D is the distance to the caustic from its convex side. Here, R s stands for the radius of curvature of section of the caustic surface by the plane P containing a light ray that is tangent to caustic at Q. The plane P also contains the vector normal to the caustic at Q (see Figure 12). The pre-factor U = U (Q) in (66) is determined by matching the geometrical optics value of magnification (12) in vicinity of the caustic (taking into account the multiplicity of images) with the following asymptotics of (66) at D (R s /k 2 ) 1/3 To find all the above values, we use the expansion in the proximity of the critical line (13) ρ = ρ c (θ) + ∆ρ.
Vector ∆ r is tangent to the caustic at Q. Therefore, for small ∆ r, the distance from the caustic to the point (X, Y, Z) equals where R a is the radius of curvature of the astroid (14) R a = 6b sin(2θ).
The point (X, Y ) has two ("strong") pre-images in the close vicinity of the critical line 28 with ∆ρ = ± 1 2b 2R a D.
28 and up to four pre-images in total, depending on the observer position relatively to the caustic.
On the other hand, it follows from (68) and (12) that away from the caustic µ = 1 4∆ρ and with account of image multiplicity and signs we get Comparing the above equation with (67) we get since the plane P intersects the z-axis under the angle α E + O(α E ) (see Figure 12). Plugging the above values of U and R s into (66), with the help of (7, 16), we get the expected final expression (52) for the near-fold GPSF.
Let us now consider the pattern in regions near the turning points θ = 0, π/2, π, 3π/2, where the intensity reaches its maximum. Without loss of generality we take the cusp at θ = 0.
According to the theory of the uniform caustic expansions (see eg [17], [20]), in the local coordinates of the section planex,ỹ (see Figure 12), where equation of the caustic has the approximate form Here, the pre-factor W is determined by matching the geometrical optics value of magnification (12) in the cusp neighborhood with the corresponding asymptotics of (71). It is convenient to setỹ = 0 and use the asymptotics of (71) forx a 6k µ(x, 0) → πW x a 6k .
The section plane is parallel to the y-axis and intersects the z-axis under the angle α E + O(α E ) (see Figure 12). Therefore,X = α Ex ,Ỹ =ỹ, where (X,Ỹ ) is the deviation from the cusp in the observer z = Z plane. Taking (73) into account, from (70) and (14) we get a = 12 b/α 3 E . From (12) it follows that forỸ = 0 andX > 0, near the cusp µ → b 2X . Then with the help of (72, 73) we obtain Substituting the above values in (71) and taking (73) into account we get the expected expression (45) for the near-cusp GPSF.

Appendix 2
The solar atmosphere introduces corrections to the gravitational deflection picture. In the first approximation, the correction α pl to the total deflection angle α → α + α pl due to refraction in the solar plasma is the cylindrically symmetric vector field (for a review see e.g. [4] and references therein) where Λ ≈ 50 meters, A ≈ 1.1, B ≈ 2.28 × 10 2 , C ≈ 2.952 × 10 3 . The atmosphere bends the rays outwards while the gravity bends the rays inwards. The above correction corresponds to the circularly symmetric contribution ψ pl = ψ pl (ρ) to the dimensionless potential ψ → ψ + ψ pl (for definitions of ψ and ρ see eqs. (9,8)) is given by (7). The effect of the solar plasma is small for the sub-micrometer wavelengths. Indeed, for λ = 10 −6 m, pl ≈ 10 −10 .
It is easy to see that due to the circular symmetry of ψ pl and smallness of pl the atmospheric refraction can be taken into account in this approximation by the formal re-normalization of the gravitational radius r g →r g (Z) in all our previous results where α pl (Z) is the α pl evaluated at r ⊥ = b(Z). Note that the correction factor α pl /α E ≈ 5 × 10 −9 for λ = 10 −6 m and Z ≈ 1000AU, i.e. effect of refraction in sub-micrometer range of the EM spectrum is extremely small in this approximation. Higher order (non-symmetric) corrections to ψ pl are proportional to the product of pl and small deformation parameters such as the oblateness coefficient of the columnar density of the solar atmosphere etc. In principle, the contribution to ψ pl corresponding to this oblateness can be accounted for through an effective (Z-dependent) correction of the quadrupole moment of the sun. This and higher moment corrections can be also discarded in the context of the present work due to their smallness for wavelengths of interest. However, the question of deformations of caustic/PSF due to fluctuations in the solar atmosphere is worth studying in the context of the high-resolution or/and longer wavelength imaging.