# Analytical Methods in Gravitational Microlensing

Authors: V. I. Zhdanov; A. N. Alexandrov; E. V. Fedorova; V. M. Sliusar
Journal: ISRN Astronomy and Astrophysics
 (2012)

---

## Abstract

We discuss analytical results dealing with photometric and astrometric gravitational microlensing. The first two sections concern approximation methods that allow us to get solutions of the general lens equation near fold caustics and cusp points up to any prescribed accuracy. Two methods of finding approximate solutions near the fold are worked out. The results are applied to derive new corrections to total amplifications of critical source images. Analytic expressions are obtained in case of the Gaussian, power-law, and limb-darkening extended source models; here we present the first nonzero corrections to the well-known linear caustic approximation. Possibilities to distinguish different source models in observations are discussed on the basis of statistical simulations of microlensed light curves. In the next section, we discuss astrometric microlensing effects in various cases of extended sources and extended lenses, including a simple model of weak statistical microlensing by extended dark matter clumps. Random walks of a distant source image microlensed by stochastic masses are estimated. We note that the bulk motion of foreground stars induces a small apparent rotation of the extragalactic reference frame. Compact analytical relations describing the statistics of such motions are presented.


---

## Body

## 1. Introduction

Gravitational lensing theory contains a number of nice analytical results and beautiful theorems, which give a description of the observed effects on the level of quality and form a basis for quantitative calculations [1, 2]. This theory is based on the general relativity that describes the propagation of electromagnetic radiation in a curved space-time. Two effects are most important: the Einsteinian light bending in the gravitational field and the gravitational time delay. These effects are described by simple formulas that can be found in any textbook on relativistic gravity. However, their use in a particular gravitational lens system (GLS) can be rather complicated, and there is a considerable need for improvement of numerical methods and analytical studies. The present article deals mainly with some of the topics of the latter subject. We believe that there exists a considerable number of interesting problems concerning GLS that can be solved either without using extensive computations or where such computations can be reduced to a minimum.One of the most spectacular effects in observations of GLS is a considerable brightness enhancement in one of the source images due to the microlensing [1–3]; this is often referred to as a high amplification event (HAE). In an extragalactic system, this is typically associated with the source crossing of the GLS caustic. It is well known that the only stable caustics of the two-dimensional mapping are folds and cusps [1, 3]. The crossing of the fold caustic is most probable (see, e.g., [4] in case of Q2237 + 0305), though the role of the cusps is also significant [1].Below we discuss approximation methods that allow us to get solutions of the lens equation near caustics with any prescribed accuracy. We consider both the cases of a fold caustic and a cusp point using a general procedure of a power reduction. A special attention is paid to the case of the fold. Application of this procedure allows us to find an analytic form of the solution which is used further. In this view we discuss two methods of approximate solution of the lens equation [5, 6]. Both methods agree with each other in a common domain of validity; they are used to obtain first orders of expansion of the approximate solutions. The results are applied to derive the amplification of a point source near the fold and then to obtain the amplifications of extended sources (the Gaussian, power-law and limb-darkening source models). In order to obtain nonzero corrections to the total amplification of two critical images, it was necessary to take into consideration additional higher-order terms in the expansion of the lens mapping in comparison with earlier works (see, e.g., [7]). In Section 3.5, we discuss the possibility to distinguish different source models on the basis of computer simulations of microlensed light curves.The following sections deal with astrometric microlensing. Here, we develop some of our earlier results [8–14]. We show that the trajectory of the image centroid depends strongly on the source size and on an external gravitational field (external shear); a microlensing by an extended lenses is also discussed including a simple “toy” model of microlensing by extended dark matter clumps. Then, we consider a statistics of weak astrometric microlensing of a distant source by a large number of foreground point masses in case of a small optical depth; this case corresponds to a weak microlensing by the Milky Way objects. Variable gravitational field of moving foreground masses induces small motions of remote sources’ images. We present compact analytical relations describing a statistics of these motions; it is different in case of a continuous (e.g., dark) matter and discrete masses (stars). This effect leads to a small apparent rotation of an extragalactic reference frame.


## 2. Analytic Theory of Caustic Crossing Events

Investigation of the lens mapping near singular points [3] is an important part of the HAE studies. The lens equation near a fold may be expanded in powers of local coordinates; in the lowest order of this expansion the caustic is represented by a straight line; so this approximation is often referred to as “the linear caustic approximation.” In this approximation, the point source flux amplification depends on the distance to the caustic and contains two parameters [1]. In most cases, the linear caustic approximation is sufficient to treat the observed light curves over the range of HAEs at modern accuracy of the flux measurement. The need for a modification of this formula, for example, by taking into account the caustic curvature, is nevertheless being discussed for a long time [15–17]. One may hope that future improvement of the photometric accuracy will make it possible to obtain additional parameters of the lens mapping, which are connected to the mass distribution in a lensing galaxy. At the same time, consideration of the “postlinear” terms is sometimes sensible in order to explain even the presently available observational data [5, 6]. Note also that the corrections to the amplification on a macrolensing level were the subject of investigations in connection with the problem of “anomalous flux ratios” [7].


### 2.1. Initial Equations and Notations

In this subsection, we recall some general notions of the gravitational lensing that may be found for example, in book [1]. We use the normalized lens equation in the form:
(1)y=x-∇Φ(x),Φ(x) is the lens mapping potential; this equation relates every point x=(x1,x2) of the image plane to the point y=(y1,y2) of the source plane. In the general case, there are several solutions X(i)(y) of the lens equation (1) that represent images of one point source at y; we denote solution number with the index in parentheses.In case of no continuous matter on the line of sight, the potential must be a harmonic functionΔΦ=0. Below, we will assume that this condition is fulfilled in the considered neighborhood of the critical point. We note, however, if we suppose that during HAE the continuous matter density is constant, the lens equation may be represented in the same form by a suitable renormalization of the variables.The amplification of a separate image of a point source is|J|-1, where J(x)≡|D(y)/D(x)| is the Jacobian of the lens mapping. The critical curves of the lens mapping (1) are defined by equation J(x)=0; they are mapped onto the caustics on the source plane. The stable critical points of two-dimensional mapping may be folds and cusps only, the folds being more probable in HAE. When a point source approaches the fold caustic from its convex side, two of its images approach the critical curve and their amplification tends to infinity; they disappear after the source crossing of the caustic. These two images are called critical. In case of the cusp point, there can be either three or one bright images near the critical curve depending on location of the source inside or outside of the caustic.


### 2.2. Methods of Approximate Solution of the Lens Equation Near the Fold

Below, we outline general procedures that allow to find approximate solutions of the lens equation near folds and cusps, the main attention is being paid to the case of the folds. The standard consideration of the caustic crossing events deals with the Taylor expansion of the potential near some pointpcr of the critical curve in the image plane. Let this point be the coordinate origin and we suppose that (1) maps pcr onto the coordinate origin on the source plane. We assume that the lens equation gives the coordinates of a point source as analytic functions of its image coordinates measured from the critical point pcr. The problem is to construct the inverse dependence of the image position as functions of the source coordinates. Since the Jacobian is zero at the critical point, these functions cannot be analytic.We consider two different methods. The first one deals with analytical expansions in powers of a small parameter, however, it results in nonanalytical functions of coordinates leading to nonintegrable terms in the amplification. The second method does not lead to such problems though it uses a somewhat more complicated representation of the solutions of the lens equation (containing roots of analytical functions). Both methods agree with each other in a common domain of validity; moreover, we use the second method to justify some expressions in the amplification formulas in terms of distributions to validate applications to the extended source models.


#### 2.2.1. Expansion of the Lens Mapping Near the Fold

We consider the folds in case of a harmonic potential, that is in case of no continuous matter on the line of sight. Nevertheless, we note that in a more general case we obtain analogous results [18].For the harmonic potential near the fold, one can write (see, e.g., [1])
(2)y1=2x1+a(x12-x22)+2bx1x2+c(x13-3x1x22)-d(x23-3x2x12)+gx24+⋯,y2=b(x12-x22)-2ax1x2+d(x13-3x1x22)+c(x23-3x2x12)+fx24+⋯,a,b,c,d,g,f are the expansion coefficients. If the y2 axis is directed toward the convexity of the caustic, then b<0 (at the fold points b≠0).First, we will use a regular procedure [6] in order to construct solutions of (2) with a desired accuracy. This procedure is useful to study the light curve of the point source, which has a trajectory crossing the fold caustic under some nonzero angle. For y2>0, we substitute
(3)yi→t2y~i;x1→t2x~1,x2→tx~2,
where i=1,2 and t may be considered as a (small) parameter describing a vicinity to the caustic. In fact for small yi and for an appropriate choice of t, we can work with y~i=O(1),x~i=O(1). This is a formal substitution that makes easier operations with different orders of the expansion. On the other hand, for small yi, the value t=1 also may be considered as “small”; in this case, yi=y~i. Therefore, in the other treatment, after performing calculations, we can put t=1 and thus return to the initial variables.Substitution of (3) into (2) yields
(4)y~1=2x~1-ax~22+tF(t,x~1,x~2),y~2=-bx~22+tG(t,x~1,x~2),
where F,G are analytical functions of all their variables. However, it should be emphasized that in fact we need finite orders of t for some approximation order; therefore, we always deal with polynomials of a finite order in powers of t,y~1,y~2,x~1,x~2.


#### 2.2.2. Direct Expansion of a Solution in Powers oft

System (4) can be considered for fixed y~i and variable t, then the solution of (4) forms a curve which is useful to study a local behavior of critical image trajectories [19]. In this case, t=0 corresponds to a crossing of a caustic by a point source and t2 may be considered as a time counted from the moment when two critical images appear. The results of [19] show that the solutions of (2) in terms of x~i can be represented as the expansions in powers of t. This can be seen directly for y2>0,b<0 if we rewrite (4) in the form ready for iterations:
(5)x~1=12[y~1-aby~2-tF1(t,x~1,x~2)],x~2=ε{-y~2b}1/2[1-ty~2G(t,x~1,x~2)+⋯],F1≡F-(a/b)G, where either ε=1 or ε=-1 is corresponding to different critical images.Therefore, from the very beginning, we can look for the solutions in the form:(6)x~1=x~10+x~11t+⋯,x~2=x~20+x~21t+⋯.From zeroth iteration, we have(7)x~10=12[y~1-aby~2],x~20=ε{+y~2|b|}1/2.
The first iteration yields
(8)x~1=x~10-t2F1(0,x~10,x~20),x~2=x~20[1-ty~2G(0,x~10,x~20)+⋯],
whence we obtain x~11,x~21 and so on.Evidently, for every step of the iterative procedure, we obtain some polynomials oft, but the dependence on y~2>0 is not analytical one. Moreover, the higher orders of approximation involve higher orders of (y~2)-1 leading to formally nonintegrable terms in the amplification factor (see below). This is a well-known phenomenon when we represent a nonanalytical expression (e.g., square root) as formal asymptotic expansion in powers of some parameter  (e.g., expansion in powers of t: x+t=x1/2+(t/2x1/2)-(t2/8x3/2)+⋯).


#### 2.2.3. Alternative Representation of Solutions

Now, we aim to obtain the other representation of the solutionx=(x1,x2) as an explicit algebraic function of analytical expansions in powers of coordinates y=(y1,y2) of the source plane. For this purpose we use a power reduction which allows to obtain a simple approximate form of the lens equation on every step of the procedure.The first equation of system (5) for sufficiently small t can be solved iteratively with respect to x~1 yielding a solution in the form:
(9)x~1=12(y~1-aby~2)+tF2(t,y~1,y~2,x~2).
In fact, we need finite orders of approximations, so only finite number of the iterations are involved. Then, after appropriate truncation, F2(t,y~1,y~2,x~2) becomes a finite-order polynomial of all the arguments. Substitution of (9) into the second equation of system (4) leads to one equation with respect to one variable x~2:
(10)y~2=-bx~22+tG1(t,y~1,y~2,x~2)
with (again after some truncation) a polynomial G1(t,x~2) of t,x~2.Now, we use the power reduction procedure so as to eliminate the terms containingx~2m with m>2 and to obtain a quadratic equation with respect to x~2. We put z=x~22 and write the terms containing x~2 as
(11)(x~2)2k=zk,(x~2)2k+1=x~2zk,k=0,1,2,….
Then, in the r.h.s. of (10), we can write (b<0 !)
(12)G1=b·(G1(1)(z)+x~2H1(1)(z)),G1(1)(z)=∑k=0n‍gkzk,H1(1)(z)=∑k=0n‍hkzk,
where to be specific we have truncated the order of the expansion by some degree n, coefficients gk,hk are polynomial functions of t,y~1,y~2. Thus, (10) can be written as
(13)z=z0+tG1(1)(z)+tx~2H1(1)(z),z0=-y~2/b. Note that the right-hand side is a linear function of x~2.Now, we solve (13) iteratively with respect to z, and every step of the iteration procedure being combined with reduction (11) so as to exclude powers of x~2 larger than 1. In the order ~O(t), we have
(14)z=z0+tG1(1)(z0)+tx~2H1(1)(z0),
In the next order, we must substitute (14) into the right-hand side of (13). This will yield the higher orders of x~2; these orders must be reduced making use of substitution (11) and the result of previous iterations. Thus, we again have an equation analogous to (14) with a linear function of x~2 in r.h.s. At the end, we obtain
(15)z=z0+tA(t,y~1,y~2)+tx~2B(t,y~1,y~2),
with analytical functions A,B of their arguments.At last, we return fromz to x~2 to obtain a quadratic equation:
(16)x~22-tx~2B(t,y~1,y~2)-z0-tA(t,y~1,y~2)=0,
yielding a solution:
(17)x~2=t2B(t,y)+ε[z0+tA(t,y)+14t2B2(t,y)]1/2.
The above approximation procedure can be fulfilled to yield approximate solutions with any degree of accuracy required. Note that in reality the calculations up to the order ~O(t2) involve a very limited number of coefficients gk,hk after truncation of higher order terms in t on every step of approximation. After derivation of x~2, we find x~1 owing to (9).From these considerations, it is easy to see that analytical structure of the solutions of the lens Equations (2) near the fold is
(18)x~1=p+tεrw,x~2=ts+εw;ε=±1,
where p,r,w,s are analytical functions of t,y~1,y~2. This allows us to seek for the solution in this form from the very beginning by substitution of (18) into the initial equations. This will be considered in more detail in Section 2.4.


### 2.3. Approximation Procedure Near the Cusp

The coefficients of Taylor expansion in (2) depend on the caustic point near which the expansion is performed. If this point is a cusp, then b=0 and the above approximation methods of Sections 2.2.2 and 2.2.3 are no longer valid. In this case, the lens equation is
(19)y1=2x1+a(x12-x22)+c(x13-3x1x22)-d(x23-3x2x12)+⋯,y2=-2ax1x2+d(x13-3x1x22)+c(x23-3x2x12)+fx24+⋯.The lower order terms correspond to a familiar form of the lens equation near the cusp as described in [1, 21].In order to perform a formal expansion, we write, instead of (3),
(20)y1→t2y~1,y2→t3y~2,x1→t2x~1,x2→tx~2.
Then the substitution into (19) yields, up to the first-order terms in t,
(21)y~1=2x~1-ax~22+t(-dx~23+⋯),y~2=-2ax~1x~2+cx~23+t(fx~24-3dx~1x~22+⋯).A detailed analysis of properties of this lens mapping in zeroth approximation (with respect tot) was carried out in [21, 22]. First-order corrections in a general case when the matter density in a vicinity of the line of sight is not zero were found in [23].The solution can be sought by using the expansion in powers oft. Here, we outline an alternative method using the power reduction like one described in Section 2.2.3, which allows determining an analytical structure of the solutions. First one must express x~1 by means of the other variables; this can be done iteratively using the first equation of system (21). Practically, we need a finite number of iterations. Then, we apply a power reduction to eliminate terms x~2m with m>3 so as to obtain a cubic equation for x~2. We put
(22)x~23k=zk,x~23k+1=zkx~2,x~23k+2=zkx~22.
On account of the second equation of system (21), after substitutions, we get an equation for x~2 of the form:
(23)z=K0(y)+K1(y)x~2+t{L0(t,y,z)+L1(t,y,z)x~2+L2(t,y,z)x~22}.
Using an iterative procedure combined with (22), this enables us to obtain z as follows:
(24)z=K0(y)+K1(y)x~2+t{L~0(t,y)+L~1(t,y)x~2+L~2(t,y)x~22},
where K0,K1,L~0,L~1,L~2 are polynomial functions of t and y, the order of the polynomials depending upon the order of the approximation required. Thus, we come to the cubic equation with respect to x~2:
(25)x~23-tL~2x~22-(K1+tL~1)x~2-(K0+tL~0)=0.
The roots can be obtained via the Cardano-Tartaglia formulas or the Viete trigonometric solution. Depending on sign of the discriminant of this equation, there are either three real roots (inside the cusp), or one real root (outside the cusp). The analytic properties of the solution are completely defined by the coefficients of (25), in particular, for t→0, they are defined by K0 and K1. It is important that for c≠a2 these coefficients are finite order polynomials of t,y~1,y~2. This structure will be preserved in any order of approximation leading to (25).


### 2.4. Amplification Near the Fold up to the Terms ~O(t2)

#### 2.4.1.  Solutions of the Lens Equation

Further, we deal with the results concerning the folds that follow from Sections2.2.2 and 2.2.3.The result of the expansions near the fold:(26)x~1=x~10+x~11t+x~12t2,x~2=x~20+x~21t+x~22t2,
as applied to the lens equation (2) in zeroth order is given by (7):
(27)x~10=12[y~1-aby~2],x~20=ε[y~2|b|]1/2.The next orders are as follows:(28)x~11=-εy~2|b|(ac-aR2+bd)y~2+bR2y~12b2,x~21=12(a2-cb2y~2-aby~1),
where R2=a2+b2. The solutions up to this accuracy level have been obtained earlier [7, 19]. The contributions of this order are cancelled in calculations of the total amplification of two critical images. Therefore, to get a nontrivial correction to zero order amplification, we need the higher order approximations.The second-order terms contain an expression, which is singular iny~2:
(29)x~12=14b4(3a5+5a3b2+2ab4-2b3d-2b2g+3ac214+b4-6a3c+3bcd-8a2bd+2abf)y~22+12b3(2a2c-b2c+3abd-2a2R2-b2R2)y~1y~2+aR24b2y~12,x~22=εy~2|b|[R28b18b3(10a2c-5c2-5a2R2+10abd-4bf)y~2-34b2(ac+bd-aR2)y~1-R28by~12y~2].Now, we proceed to second approach to construct approximate solutions of the lens equation in a vicinity of folds, which is free from such singularities. We look for two critical solutions (18) x~1=p+tεrw,x~2=ts+εw corresponding to different signs of ε=±1. After substitution into the lens equation, we separate the terms containing integer and half-integer powers of w, For example,
(30)∑n,m‍an,m(p+εtrw)n(ts+εw)m=A0+εw1/2A1,
where A0,A1 are analytical in t,p,r,s,w. This yields a system of four independent equations for the variables p,r,w,s; this system can be reduced to a form convenient for iteration procedure [5]. It is important to note that at every iteration step we obtain an approximate solution in the form of finite order polynomials of t,y~1,y~2, as it was stated at the end of Section 2.2.3. The solution of this system up to the terms ~t2 is [5]
(31)p=x~10+t2x~12,r=-R2y~12b+y~22b2[aR2-(ac+bd)],s=-a2by~1+a2-c2b2y~2,w=-y~2+t2Zb,Z=-R24by~12+32b2[a(a2-c)+b(ab-d)]y~1y~2-y~224b3[5a2(R2-2c)+5c2-10abd+4bf].


#### 2.4.2. Total Amplification of Critical Images of a Point Source Near the Fold

The solutions of the lens equation are then used to derive the Jacobians of the lens mapping (for both images near the critical curve). The value ofJ-1 yields the amplification of individual images. As we pointed out above, we need the total amplification of two critical images (the sum of two amplifications of the separate critical images). In the second-order approximation (using the expansion up to the terms t2), this is
(32)Kcr=12Θ(y2)|b|y2[1+Py2+Qy1-κ4y12y2],
where the constants P,Q, are expressed via the Taylor expansion coefficients from (1), and
(33)κ=a2+b22|b|;Θ(y2) is the Heaviside step function. Note that κ is the caustic curvature at the origin which enters explicitly into the amplification formula. Parameters P and Q are independent; explicit formulae for them may be found in [5, 6]. However, this is not needed when we use (32) for fitting the observational data, because these constants are whatever considered as free fitting parameters.Formula (32) yields an effective approximation for the point source amplification near the coordinate origin provided that y2>0, and y2/y12 is not too small (see the term containing κ). For a fixed source position, this can be satisfied always by an appropriate choice of the coordinate origin, so that the source will be situated almost on a normal to the tangent to the caustic.If the source is on the caustic tangent or in the region between the caustic and the tangent, then formula (32) does not represent a good approximation to the point source amplification. Nevertheless, in case of an extended source, we will show that result (32) can be used to obtain approximations to the amplification of this source even as it intersects the caustic. However, to do this, we need to redefine correctly the convolution of (32) with a brightness distribution.


#### 2.4.3. Dark Matter on the Line of Sight

It is well known that the nonbaryonic dark matter (DM) dominates in galactic masses, though exact small-scale distribution of DM is a subject of studies. Simulations show that a complicated subhalo structure is possible [24]. In this connection, it is interesting to study how a variable DM density can affect the HAE characteristics.Formulae (27), (32) have been obtained under an assumption that there is no continuous matter on the line of sight k=0. It is easy to take into account an effect of a nonzero constant density k=const≠0 by means of a rescaling of variables and change of coefficients in  (2) [5].Fork≠0, the lens potential in (1) satisfies equation ΔΦ=2k(x). The Taylor expansion analogous to (2) after substitution (3) takes on the form in the second approximation:
(34)y~1=2(1-k0)x~1-a2x~22+t(2b1x~1x~2-dx~23)+t2(a1x~12-3c1x~1x~22+gx~24),y~2=-b2x~22+t(-2a2x~1x~2+c2x~23)+t2(b1x~12-3dx~1x~22+fx~24).
Here, k0=k(0) is the matter density at the origin, and the coefficients of expansion in (34) are expressed by means of derivatives of the lens potential Φ. If k is constant, then a1=a2=a, b1=b2=b, c1=c2=c. Thus, system (34) contains four additional parameters as compared to the previously studied case k(x)≡0. The solutions of system (34) have been studied in [18]; they are analogous to (27)–(29) and they have the same analytic structure. The formula for the amplification also preserves its form (32); the difference is only due to a change of explicit expressions for P,Q and κ in terms of new Taylor coefficients, which anyway cannot be determined from observations for realistic models of mass distributions. We, however, must make a reservation that here we suppose that the Taylor expansion in the lens equation is possible and it is effective (which presupposes that the continuous matter is sufficiently smooth); therefore, very inhomogeneous case of small objects with size ≪RE is not involved in our consideration.


### 2.5. Resume of Section2

The main results of this section deal with solutions of the lens equation in the caustic region, namely, near the folds and cusps. We propose a general procedure of the power reduction that provides an approximate solution with a prescribed accuracy. Most attention is paid to the case of the fold caustic. We outlined two methods that enable us to obtain the critical solutions of the gravitational lens equation near a fold with any desired accuracy.In order to obtain nontrivial corrections toKcr near the fold obtained in the linear caustic approximation, the higher orders of the expansion of the lens equation must be taken into account as compared to works [7, 19]. The modified formula for Kcr contains 3 extra parameters in addition to those of the linear caustic approximation. We point out that any presence of a continuous (cf. dark) matter on the on the line of sight (with the same typical scales of Taylor expansion) does not change analytical structure and the number of fitting parameters in the formula for Kcr.


## 2.1. Initial Equations and Notations

In this subsection, we recall some general notions of the gravitational lensing that may be found for example, in book [1]. We use the normalized lens equation in the form:
(1)y=x-∇Φ(x),Φ(x) is the lens mapping potential; this equation relates every point x=(x1,x2) of the image plane to the point y=(y1,y2) of the source plane. In the general case, there are several solutions X(i)(y) of the lens equation (1) that represent images of one point source at y; we denote solution number with the index in parentheses.In case of no continuous matter on the line of sight, the potential must be a harmonic functionΔΦ=0. Below, we will assume that this condition is fulfilled in the considered neighborhood of the critical point. We note, however, if we suppose that during HAE the continuous matter density is constant, the lens equation may be represented in the same form by a suitable renormalization of the variables.The amplification of a separate image of a point source is|J|-1, where J(x)≡|D(y)/D(x)| is the Jacobian of the lens mapping. The critical curves of the lens mapping (1) are defined by equation J(x)=0; they are mapped onto the caustics on the source plane. The stable critical points of two-dimensional mapping may be folds and cusps only, the folds being more probable in HAE. When a point source approaches the fold caustic from its convex side, two of its images approach the critical curve and their amplification tends to infinity; they disappear after the source crossing of the caustic. These two images are called critical. In case of the cusp point, there can be either three or one bright images near the critical curve depending on location of the source inside or outside of the caustic.


## 2.2. Methods of Approximate Solution of the Lens Equation Near the Fold

Below, we outline general procedures that allow to find approximate solutions of the lens equation near folds and cusps, the main attention is being paid to the case of the folds. The standard consideration of the caustic crossing events deals with the Taylor expansion of the potential near some pointpcr of the critical curve in the image plane. Let this point be the coordinate origin and we suppose that (1) maps pcr onto the coordinate origin on the source plane. We assume that the lens equation gives the coordinates of a point source as analytic functions of its image coordinates measured from the critical point pcr. The problem is to construct the inverse dependence of the image position as functions of the source coordinates. Since the Jacobian is zero at the critical point, these functions cannot be analytic.We consider two different methods. The first one deals with analytical expansions in powers of a small parameter, however, it results in nonanalytical functions of coordinates leading to nonintegrable terms in the amplification. The second method does not lead to such problems though it uses a somewhat more complicated representation of the solutions of the lens equation (containing roots of analytical functions). Both methods agree with each other in a common domain of validity; moreover, we use the second method to justify some expressions in the amplification formulas in terms of distributions to validate applications to the extended source models.


### 2.2.1. Expansion of the Lens Mapping Near the Fold

We consider the folds in case of a harmonic potential, that is in case of no continuous matter on the line of sight. Nevertheless, we note that in a more general case we obtain analogous results [18].For the harmonic potential near the fold, one can write (see, e.g., [1])
(2)y1=2x1+a(x12-x22)+2bx1x2+c(x13-3x1x22)-d(x23-3x2x12)+gx24+⋯,y2=b(x12-x22)-2ax1x2+d(x13-3x1x22)+c(x23-3x2x12)+fx24+⋯,a,b,c,d,g,f are the expansion coefficients. If the y2 axis is directed toward the convexity of the caustic, then b<0 (at the fold points b≠0).First, we will use a regular procedure [6] in order to construct solutions of (2) with a desired accuracy. This procedure is useful to study the light curve of the point source, which has a trajectory crossing the fold caustic under some nonzero angle. For y2>0, we substitute
(3)yi→t2y~i;x1→t2x~1,x2→tx~2,
where i=1,2 and t may be considered as a (small) parameter describing a vicinity to the caustic. In fact for small yi and for an appropriate choice of t, we can work with y~i=O(1),x~i=O(1). This is a formal substitution that makes easier operations with different orders of the expansion. On the other hand, for small yi, the value t=1 also may be considered as “small”; in this case, yi=y~i. Therefore, in the other treatment, after performing calculations, we can put t=1 and thus return to the initial variables.Substitution of (3) into (2) yields
(4)y~1=2x~1-ax~22+tF(t,x~1,x~2),y~2=-bx~22+tG(t,x~1,x~2),
where F,G are analytical functions of all their variables. However, it should be emphasized that in fact we need finite orders of t for some approximation order; therefore, we always deal with polynomials of a finite order in powers of t,y~1,y~2,x~1,x~2.


### 2.2.2. Direct Expansion of a Solution in Powers oft

System (4) can be considered for fixed y~i and variable t, then the solution of (4) forms a curve which is useful to study a local behavior of critical image trajectories [19]. In this case, t=0 corresponds to a crossing of a caustic by a point source and t2 may be considered as a time counted from the moment when two critical images appear. The results of [19] show that the solutions of (2) in terms of x~i can be represented as the expansions in powers of t. This can be seen directly for y2>0,b<0 if we rewrite (4) in the form ready for iterations:
(5)x~1=12[y~1-aby~2-tF1(t,x~1,x~2)],x~2=ε{-y~2b}1/2[1-ty~2G(t,x~1,x~2)+⋯],F1≡F-(a/b)G, where either ε=1 or ε=-1 is corresponding to different critical images.Therefore, from the very beginning, we can look for the solutions in the form:(6)x~1=x~10+x~11t+⋯,x~2=x~20+x~21t+⋯.From zeroth iteration, we have(7)x~10=12[y~1-aby~2],x~20=ε{+y~2|b|}1/2.
The first iteration yields
(8)x~1=x~10-t2F1(0,x~10,x~20),x~2=x~20[1-ty~2G(0,x~10,x~20)+⋯],
whence we obtain x~11,x~21 and so on.Evidently, for every step of the iterative procedure, we obtain some polynomials oft, but the dependence on y~2>0 is not analytical one. Moreover, the higher orders of approximation involve higher orders of (y~2)-1 leading to formally nonintegrable terms in the amplification factor (see below). This is a well-known phenomenon when we represent a nonanalytical expression (e.g., square root) as formal asymptotic expansion in powers of some parameter  (e.g., expansion in powers of t: x+t=x1/2+(t/2x1/2)-(t2/8x3/2)+⋯).


### 2.2.3. Alternative Representation of Solutions

Now, we aim to obtain the other representation of the solutionx=(x1,x2) as an explicit algebraic function of analytical expansions in powers of coordinates y=(y1,y2) of the source plane. For this purpose we use a power reduction which allows to obtain a simple approximate form of the lens equation on every step of the procedure.The first equation of system (5) for sufficiently small t can be solved iteratively with respect to x~1 yielding a solution in the form:
(9)x~1=12(y~1-aby~2)+tF2(t,y~1,y~2,x~2).
In fact, we need finite orders of approximations, so only finite number of the iterations are involved. Then, after appropriate truncation, F2(t,y~1,y~2,x~2) becomes a finite-order polynomial of all the arguments. Substitution of (9) into the second equation of system (4) leads to one equation with respect to one variable x~2:
(10)y~2=-bx~22+tG1(t,y~1,y~2,x~2)
with (again after some truncation) a polynomial G1(t,x~2) of t,x~2.Now, we use the power reduction procedure so as to eliminate the terms containingx~2m with m>2 and to obtain a quadratic equation with respect to x~2. We put z=x~22 and write the terms containing x~2 as
(11)(x~2)2k=zk,(x~2)2k+1=x~2zk,k=0,1,2,….
Then, in the r.h.s. of (10), we can write (b<0 !)
(12)G1=b·(G1(1)(z)+x~2H1(1)(z)),G1(1)(z)=∑k=0n‍gkzk,H1(1)(z)=∑k=0n‍hkzk,
where to be specific we have truncated the order of the expansion by some degree n, coefficients gk,hk are polynomial functions of t,y~1,y~2. Thus, (10) can be written as
(13)z=z0+tG1(1)(z)+tx~2H1(1)(z),z0=-y~2/b. Note that the right-hand side is a linear function of x~2.Now, we solve (13) iteratively with respect to z, and every step of the iteration procedure being combined with reduction (11) so as to exclude powers of x~2 larger than 1. In the order ~O(t), we have
(14)z=z0+tG1(1)(z0)+tx~2H1(1)(z0),
In the next order, we must substitute (14) into the right-hand side of (13). This will yield the higher orders of x~2; these orders must be reduced making use of substitution (11) and the result of previous iterations. Thus, we again have an equation analogous to (14) with a linear function of x~2 in r.h.s. At the end, we obtain
(15)z=z0+tA(t,y~1,y~2)+tx~2B(t,y~1,y~2),
with analytical functions A,B of their arguments.At last, we return fromz to x~2 to obtain a quadratic equation:
(16)x~22-tx~2B(t,y~1,y~2)-z0-tA(t,y~1,y~2)=0,
yielding a solution:
(17)x~2=t2B(t,y)+ε[z0+tA(t,y)+14t2B2(t,y)]1/2.
The above approximation procedure can be fulfilled to yield approximate solutions with any degree of accuracy required. Note that in reality the calculations up to the order ~O(t2) involve a very limited number of coefficients gk,hk after truncation of higher order terms in t on every step of approximation. After derivation of x~2, we find x~1 owing to (9).From these considerations, it is easy to see that analytical structure of the solutions of the lens Equations (2) near the fold is
(18)x~1=p+tεrw,x~2=ts+εw;ε=±1,
where p,r,w,s are analytical functions of t,y~1,y~2. This allows us to seek for the solution in this form from the very beginning by substitution of (18) into the initial equations. This will be considered in more detail in Section 2.4.


## 2.2.1. Expansion of the Lens Mapping Near the Fold

We consider the folds in case of a harmonic potential, that is in case of no continuous matter on the line of sight. Nevertheless, we note that in a more general case we obtain analogous results [18].For the harmonic potential near the fold, one can write (see, e.g., [1])
(2)y1=2x1+a(x12-x22)+2bx1x2+c(x13-3x1x22)-d(x23-3x2x12)+gx24+⋯,y2=b(x12-x22)-2ax1x2+d(x13-3x1x22)+c(x23-3x2x12)+fx24+⋯,a,b,c,d,g,f are the expansion coefficients. If the y2 axis is directed toward the convexity of the caustic, then b<0 (at the fold points b≠0).First, we will use a regular procedure [6] in order to construct solutions of (2) with a desired accuracy. This procedure is useful to study the light curve of the point source, which has a trajectory crossing the fold caustic under some nonzero angle. For y2>0, we substitute
(3)yi→t2y~i;x1→t2x~1,x2→tx~2,
where i=1,2 and t may be considered as a (small) parameter describing a vicinity to the caustic. In fact for small yi and for an appropriate choice of t, we can work with y~i=O(1),x~i=O(1). This is a formal substitution that makes easier operations with different orders of the expansion. On the other hand, for small yi, the value t=1 also may be considered as “small”; in this case, yi=y~i. Therefore, in the other treatment, after performing calculations, we can put t=1 and thus return to the initial variables.Substitution of (3) into (2) yields
(4)y~1=2x~1-ax~22+tF(t,x~1,x~2),y~2=-bx~22+tG(t,x~1,x~2),
where F,G are analytical functions of all their variables. However, it should be emphasized that in fact we need finite orders of t for some approximation order; therefore, we always deal with polynomials of a finite order in powers of t,y~1,y~2,x~1,x~2.


## 2.2.2. Direct Expansion of a Solution in Powers oft

System (4) can be considered for fixed y~i and variable t, then the solution of (4) forms a curve which is useful to study a local behavior of critical image trajectories [19]. In this case, t=0 corresponds to a crossing of a caustic by a point source and t2 may be considered as a time counted from the moment when two critical images appear. The results of [19] show that the solutions of (2) in terms of x~i can be represented as the expansions in powers of t. This can be seen directly for y2>0,b<0 if we rewrite (4) in the form ready for iterations:
(5)x~1=12[y~1-aby~2-tF1(t,x~1,x~2)],x~2=ε{-y~2b}1/2[1-ty~2G(t,x~1,x~2)+⋯],F1≡F-(a/b)G, where either ε=1 or ε=-1 is corresponding to different critical images.Therefore, from the very beginning, we can look for the solutions in the form:(6)x~1=x~10+x~11t+⋯,x~2=x~20+x~21t+⋯.From zeroth iteration, we have(7)x~10=12[y~1-aby~2],x~20=ε{+y~2|b|}1/2.
The first iteration yields
(8)x~1=x~10-t2F1(0,x~10,x~20),x~2=x~20[1-ty~2G(0,x~10,x~20)+⋯],
whence we obtain x~11,x~21 and so on.Evidently, for every step of the iterative procedure, we obtain some polynomials oft, but the dependence on y~2>0 is not analytical one. Moreover, the higher orders of approximation involve higher orders of (y~2)-1 leading to formally nonintegrable terms in the amplification factor (see below). This is a well-known phenomenon when we represent a nonanalytical expression (e.g., square root) as formal asymptotic expansion in powers of some parameter  (e.g., expansion in powers of t: x+t=x1/2+(t/2x1/2)-(t2/8x3/2)+⋯).


## 2.2.3. Alternative Representation of Solutions

Now, we aim to obtain the other representation of the solutionx=(x1,x2) as an explicit algebraic function of analytical expansions in powers of coordinates y=(y1,y2) of the source plane. For this purpose we use a power reduction which allows to obtain a simple approximate form of the lens equation on every step of the procedure.The first equation of system (5) for sufficiently small t can be solved iteratively with respect to x~1 yielding a solution in the form:
(9)x~1=12(y~1-aby~2)+tF2(t,y~1,y~2,x~2).
In fact, we need finite orders of approximations, so only finite number of the iterations are involved. Then, after appropriate truncation, F2(t,y~1,y~2,x~2) becomes a finite-order polynomial of all the arguments. Substitution of (9) into the second equation of system (4) leads to one equation with respect to one variable x~2:
(10)y~2=-bx~22+tG1(t,y~1,y~2,x~2)
with (again after some truncation) a polynomial G1(t,x~2) of t,x~2.Now, we use the power reduction procedure so as to eliminate the terms containingx~2m with m>2 and to obtain a quadratic equation with respect to x~2. We put z=x~22 and write the terms containing x~2 as
(11)(x~2)2k=zk,(x~2)2k+1=x~2zk,k=0,1,2,….
Then, in the r.h.s. of (10), we can write (b<0 !)
(12)G1=b·(G1(1)(z)+x~2H1(1)(z)),G1(1)(z)=∑k=0n‍gkzk,H1(1)(z)=∑k=0n‍hkzk,
where to be specific we have truncated the order of the expansion by some degree n, coefficients gk,hk are polynomial functions of t,y~1,y~2. Thus, (10) can be written as
(13)z=z0+tG1(1)(z)+tx~2H1(1)(z),z0=-y~2/b. Note that the right-hand side is a linear function of x~2.Now, we solve (13) iteratively with respect to z, and every step of the iteration procedure being combined with reduction (11) so as to exclude powers of x~2 larger than 1. In the order ~O(t), we have
(14)z=z0+tG1(1)(z0)+tx~2H1(1)(z0),
In the next order, we must substitute (14) into the right-hand side of (13). This will yield the higher orders of x~2; these orders must be reduced making use of substitution (11) and the result of previous iterations. Thus, we again have an equation analogous to (14) with a linear function of x~2 in r.h.s. At the end, we obtain
(15)z=z0+tA(t,y~1,y~2)+tx~2B(t,y~1,y~2),
with analytical functions A,B of their arguments.At last, we return fromz to x~2 to obtain a quadratic equation:
(16)x~22-tx~2B(t,y~1,y~2)-z0-tA(t,y~1,y~2)=0,
yielding a solution:
(17)x~2=t2B(t,y)+ε[z0+tA(t,y)+14t2B2(t,y)]1/2.
The above approximation procedure can be fulfilled to yield approximate solutions with any degree of accuracy required. Note that in reality the calculations up to the order ~O(t2) involve a very limited number of coefficients gk,hk after truncation of higher order terms in t on every step of approximation. After derivation of x~2, we find x~1 owing to (9).From these considerations, it is easy to see that analytical structure of the solutions of the lens Equations (2) near the fold is
(18)x~1=p+tεrw,x~2=ts+εw;ε=±1,
where p,r,w,s are analytical functions of t,y~1,y~2. This allows us to seek for the solution in this form from the very beginning by substitution of (18) into the initial equations. This will be considered in more detail in Section 2.4.


## 2.3. Approximation Procedure Near the Cusp

The coefficients of Taylor expansion in (2) depend on the caustic point near which the expansion is performed. If this point is a cusp, then b=0 and the above approximation methods of Sections 2.2.2 and 2.2.3 are no longer valid. In this case, the lens equation is
(19)y1=2x1+a(x12-x22)+c(x13-3x1x22)-d(x23-3x2x12)+⋯,y2=-2ax1x2+d(x13-3x1x22)+c(x23-3x2x12)+fx24+⋯.The lower order terms correspond to a familiar form of the lens equation near the cusp as described in [1, 21].In order to perform a formal expansion, we write, instead of (3),
(20)y1→t2y~1,y2→t3y~2,x1→t2x~1,x2→tx~2.
Then the substitution into (19) yields, up to the first-order terms in t,
(21)y~1=2x~1-ax~22+t(-dx~23+⋯),y~2=-2ax~1x~2+cx~23+t(fx~24-3dx~1x~22+⋯).A detailed analysis of properties of this lens mapping in zeroth approximation (with respect tot) was carried out in [21, 22]. First-order corrections in a general case when the matter density in a vicinity of the line of sight is not zero were found in [23].The solution can be sought by using the expansion in powers oft. Here, we outline an alternative method using the power reduction like one described in Section 2.2.3, which allows determining an analytical structure of the solutions. First one must express x~1 by means of the other variables; this can be done iteratively using the first equation of system (21). Practically, we need a finite number of iterations. Then, we apply a power reduction to eliminate terms x~2m with m>3 so as to obtain a cubic equation for x~2. We put
(22)x~23k=zk,x~23k+1=zkx~2,x~23k+2=zkx~22.
On account of the second equation of system (21), after substitutions, we get an equation for x~2 of the form:
(23)z=K0(y)+K1(y)x~2+t{L0(t,y,z)+L1(t,y,z)x~2+L2(t,y,z)x~22}.
Using an iterative procedure combined with (22), this enables us to obtain z as follows:
(24)z=K0(y)+K1(y)x~2+t{L~0(t,y)+L~1(t,y)x~2+L~2(t,y)x~22},
where K0,K1,L~0,L~1,L~2 are polynomial functions of t and y, the order of the polynomials depending upon the order of the approximation required. Thus, we come to the cubic equation with respect to x~2:
(25)x~23-tL~2x~22-(K1+tL~1)x~2-(K0+tL~0)=0.
The roots can be obtained via the Cardano-Tartaglia formulas or the Viete trigonometric solution. Depending on sign of the discriminant of this equation, there are either three real roots (inside the cusp), or one real root (outside the cusp). The analytic properties of the solution are completely defined by the coefficients of (25), in particular, for t→0, they are defined by K0 and K1. It is important that for c≠a2 these coefficients are finite order polynomials of t,y~1,y~2. This structure will be preserved in any order of approximation leading to (25).


## 2.4. Amplification Near the Fold up to the Terms ~O(t2)

### 2.4.1.  Solutions of the Lens Equation

Further, we deal with the results concerning the folds that follow from Sections2.2.2 and 2.2.3.The result of the expansions near the fold:(26)x~1=x~10+x~11t+x~12t2,x~2=x~20+x~21t+x~22t2,
as applied to the lens equation (2) in zeroth order is given by (7):
(27)x~10=12[y~1-aby~2],x~20=ε[y~2|b|]1/2.The next orders are as follows:(28)x~11=-εy~2|b|(ac-aR2+bd)y~2+bR2y~12b2,x~21=12(a2-cb2y~2-aby~1),
where R2=a2+b2. The solutions up to this accuracy level have been obtained earlier [7, 19]. The contributions of this order are cancelled in calculations of the total amplification of two critical images. Therefore, to get a nontrivial correction to zero order amplification, we need the higher order approximations.The second-order terms contain an expression, which is singular iny~2:
(29)x~12=14b4(3a5+5a3b2+2ab4-2b3d-2b2g+3ac214+b4-6a3c+3bcd-8a2bd+2abf)y~22+12b3(2a2c-b2c+3abd-2a2R2-b2R2)y~1y~2+aR24b2y~12,x~22=εy~2|b|[R28b18b3(10a2c-5c2-5a2R2+10abd-4bf)y~2-34b2(ac+bd-aR2)y~1-R28by~12y~2].Now, we proceed to second approach to construct approximate solutions of the lens equation in a vicinity of folds, which is free from such singularities. We look for two critical solutions (18) x~1=p+tεrw,x~2=ts+εw corresponding to different signs of ε=±1. After substitution into the lens equation, we separate the terms containing integer and half-integer powers of w, For example,
(30)∑n,m‍an,m(p+εtrw)n(ts+εw)m=A0+εw1/2A1,
where A0,A1 are analytical in t,p,r,s,w. This yields a system of four independent equations for the variables p,r,w,s; this system can be reduced to a form convenient for iteration procedure [5]. It is important to note that at every iteration step we obtain an approximate solution in the form of finite order polynomials of t,y~1,y~2, as it was stated at the end of Section 2.2.3. The solution of this system up to the terms ~t2 is [5]
(31)p=x~10+t2x~12,r=-R2y~12b+y~22b2[aR2-(ac+bd)],s=-a2by~1+a2-c2b2y~2,w=-y~2+t2Zb,Z=-R24by~12+32b2[a(a2-c)+b(ab-d)]y~1y~2-y~224b3[5a2(R2-2c)+5c2-10abd+4bf].


### 2.4.2. Total Amplification of Critical Images of a Point Source Near the Fold

The solutions of the lens equation are then used to derive the Jacobians of the lens mapping (for both images near the critical curve). The value ofJ-1 yields the amplification of individual images. As we pointed out above, we need the total amplification of two critical images (the sum of two amplifications of the separate critical images). In the second-order approximation (using the expansion up to the terms t2), this is
(32)Kcr=12Θ(y2)|b|y2[1+Py2+Qy1-κ4y12y2],
where the constants P,Q, are expressed via the Taylor expansion coefficients from (1), and
(33)κ=a2+b22|b|;Θ(y2) is the Heaviside step function. Note that κ is the caustic curvature at the origin which enters explicitly into the amplification formula. Parameters P and Q are independent; explicit formulae for them may be found in [5, 6]. However, this is not needed when we use (32) for fitting the observational data, because these constants are whatever considered as free fitting parameters.Formula (32) yields an effective approximation for the point source amplification near the coordinate origin provided that y2>0, and y2/y12 is not too small (see the term containing κ). For a fixed source position, this can be satisfied always by an appropriate choice of the coordinate origin, so that the source will be situated almost on a normal to the tangent to the caustic.If the source is on the caustic tangent or in the region between the caustic and the tangent, then formula (32) does not represent a good approximation to the point source amplification. Nevertheless, in case of an extended source, we will show that result (32) can be used to obtain approximations to the amplification of this source even as it intersects the caustic. However, to do this, we need to redefine correctly the convolution of (32) with a brightness distribution.


### 2.4.3. Dark Matter on the Line of Sight

It is well known that the nonbaryonic dark matter (DM) dominates in galactic masses, though exact small-scale distribution of DM is a subject of studies. Simulations show that a complicated subhalo structure is possible [24]. In this connection, it is interesting to study how a variable DM density can affect the HAE characteristics.Formulae (27), (32) have been obtained under an assumption that there is no continuous matter on the line of sight k=0. It is easy to take into account an effect of a nonzero constant density k=const≠0 by means of a rescaling of variables and change of coefficients in  (2) [5].Fork≠0, the lens potential in (1) satisfies equation ΔΦ=2k(x). The Taylor expansion analogous to (2) after substitution (3) takes on the form in the second approximation:
(34)y~1=2(1-k0)x~1-a2x~22+t(2b1x~1x~2-dx~23)+t2(a1x~12-3c1x~1x~22+gx~24),y~2=-b2x~22+t(-2a2x~1x~2+c2x~23)+t2(b1x~12-3dx~1x~22+fx~24).
Here, k0=k(0) is the matter density at the origin, and the coefficients of expansion in (34) are expressed by means of derivatives of the lens potential Φ. If k is constant, then a1=a2=a, b1=b2=b, c1=c2=c. Thus, system (34) contains four additional parameters as compared to the previously studied case k(x)≡0. The solutions of system (34) have been studied in [18]; they are analogous to (27)–(29) and they have the same analytic structure. The formula for the amplification also preserves its form (32); the difference is only due to a change of explicit expressions for P,Q and κ in terms of new Taylor coefficients, which anyway cannot be determined from observations for realistic models of mass distributions. We, however, must make a reservation that here we suppose that the Taylor expansion in the lens equation is possible and it is effective (which presupposes that the continuous matter is sufficiently smooth); therefore, very inhomogeneous case of small objects with size ≪RE is not involved in our consideration.


## 2.4.1.  Solutions of the Lens Equation

Further, we deal with the results concerning the folds that follow from Sections2.2.2 and 2.2.3.The result of the expansions near the fold:(26)x~1=x~10+x~11t+x~12t2,x~2=x~20+x~21t+x~22t2,
as applied to the lens equation (2) in zeroth order is given by (7):
(27)x~10=12[y~1-aby~2],x~20=ε[y~2|b|]1/2.The next orders are as follows:(28)x~11=-εy~2|b|(ac-aR2+bd)y~2+bR2y~12b2,x~21=12(a2-cb2y~2-aby~1),
where R2=a2+b2. The solutions up to this accuracy level have been obtained earlier [7, 19]. The contributions of this order are cancelled in calculations of the total amplification of two critical images. Therefore, to get a nontrivial correction to zero order amplification, we need the higher order approximations.The second-order terms contain an expression, which is singular iny~2:
(29)x~12=14b4(3a5+5a3b2+2ab4-2b3d-2b2g+3ac214+b4-6a3c+3bcd-8a2bd+2abf)y~22+12b3(2a2c-b2c+3abd-2a2R2-b2R2)y~1y~2+aR24b2y~12,x~22=εy~2|b|[R28b18b3(10a2c-5c2-5a2R2+10abd-4bf)y~2-34b2(ac+bd-aR2)y~1-R28by~12y~2].Now, we proceed to second approach to construct approximate solutions of the lens equation in a vicinity of folds, which is free from such singularities. We look for two critical solutions (18) x~1=p+tεrw,x~2=ts+εw corresponding to different signs of ε=±1. After substitution into the lens equation, we separate the terms containing integer and half-integer powers of w, For example,
(30)∑n,m‍an,m(p+εtrw)n(ts+εw)m=A0+εw1/2A1,
where A0,A1 are analytical in t,p,r,s,w. This yields a system of four independent equations for the variables p,r,w,s; this system can be reduced to a form convenient for iteration procedure [5]. It is important to note that at every iteration step we obtain an approximate solution in the form of finite order polynomials of t,y~1,y~2, as it was stated at the end of Section 2.2.3. The solution of this system up to the terms ~t2 is [5]
(31)p=x~10+t2x~12,r=-R2y~12b+y~22b2[aR2-(ac+bd)],s=-a2by~1+a2-c2b2y~2,w=-y~2+t2Zb,Z=-R24by~12+32b2[a(a2-c)+b(ab-d)]y~1y~2-y~224b3[5a2(R2-2c)+5c2-10abd+4bf].


## 2.4.2. Total Amplification of Critical Images of a Point Source Near the Fold

The solutions of the lens equation are then used to derive the Jacobians of the lens mapping (for both images near the critical curve). The value ofJ-1 yields the amplification of individual images. As we pointed out above, we need the total amplification of two critical images (the sum of two amplifications of the separate critical images). In the second-order approximation (using the expansion up to the terms t2), this is
(32)Kcr=12Θ(y2)|b|y2[1+Py2+Qy1-κ4y12y2],
where the constants P,Q, are expressed via the Taylor expansion coefficients from (1), and
(33)κ=a2+b22|b|;Θ(y2) is the Heaviside step function. Note that κ is the caustic curvature at the origin which enters explicitly into the amplification formula. Parameters P and Q are independent; explicit formulae for them may be found in [5, 6]. However, this is not needed when we use (32) for fitting the observational data, because these constants are whatever considered as free fitting parameters.Formula (32) yields an effective approximation for the point source amplification near the coordinate origin provided that y2>0, and y2/y12 is not too small (see the term containing κ). For a fixed source position, this can be satisfied always by an appropriate choice of the coordinate origin, so that the source will be situated almost on a normal to the tangent to the caustic.If the source is on the caustic tangent or in the region between the caustic and the tangent, then formula (32) does not represent a good approximation to the point source amplification. Nevertheless, in case of an extended source, we will show that result (32) can be used to obtain approximations to the amplification of this source even as it intersects the caustic. However, to do this, we need to redefine correctly the convolution of (32) with a brightness distribution.


## 2.4.3. Dark Matter on the Line of Sight

It is well known that the nonbaryonic dark matter (DM) dominates in galactic masses, though exact small-scale distribution of DM is a subject of studies. Simulations show that a complicated subhalo structure is possible [24]. In this connection, it is interesting to study how a variable DM density can affect the HAE characteristics.Formulae (27), (32) have been obtained under an assumption that there is no continuous matter on the line of sight k=0. It is easy to take into account an effect of a nonzero constant density k=const≠0 by means of a rescaling of variables and change of coefficients in  (2) [5].Fork≠0, the lens potential in (1) satisfies equation ΔΦ=2k(x). The Taylor expansion analogous to (2) after substitution (3) takes on the form in the second approximation:
(34)y~1=2(1-k0)x~1-a2x~22+t(2b1x~1x~2-dx~23)+t2(a1x~12-3c1x~1x~22+gx~24),y~2=-b2x~22+t(-2a2x~1x~2+c2x~23)+t2(b1x~12-3dx~1x~22+fx~24).
Here, k0=k(0) is the matter density at the origin, and the coefficients of expansion in (34) are expressed by means of derivatives of the lens potential Φ. If k is constant, then a1=a2=a, b1=b2=b, c1=c2=c. Thus, system (34) contains four additional parameters as compared to the previously studied case k(x)≡0. The solutions of system (34) have been studied in [18]; they are analogous to (27)–(29) and they have the same analytic structure. The formula for the amplification also preserves its form (32); the difference is only due to a change of explicit expressions for P,Q and κ in terms of new Taylor coefficients, which anyway cannot be determined from observations for realistic models of mass distributions. We, however, must make a reservation that here we suppose that the Taylor expansion in the lens equation is possible and it is effective (which presupposes that the continuous matter is sufficiently smooth); therefore, very inhomogeneous case of small objects with size ≪RE is not involved in our consideration.


## 2.5. Resume of Section2

The main results of this section deal with solutions of the lens equation in the caustic region, namely, near the folds and cusps. We propose a general procedure of the power reduction that provides an approximate solution with a prescribed accuracy. Most attention is paid to the case of the fold caustic. We outlined two methods that enable us to obtain the critical solutions of the gravitational lens equation near a fold with any desired accuracy.In order to obtain nontrivial corrections toKcr near the fold obtained in the linear caustic approximation, the higher orders of the expansion of the lens equation must be taken into account as compared to works [7, 19]. The modified formula for Kcr contains 3 extra parameters in addition to those of the linear caustic approximation. We point out that any presence of a continuous (cf. dark) matter on the on the line of sight (with the same typical scales of Taylor expansion) does not change analytical structure and the number of fitting parameters in the formula for Kcr.


## 3. Amplification of Extended Sources

### 3.1.   Preliminary Comments

In this section, we proceed to applications of the previous results to some extended source models typically used in fittings of the observed light curves. Interest to HAE in extragalactic GLSs is due to possibilities to study a brightness distributions in sources. This is especially interesting in connection with investigation of central regions of distant quasars. The idea, first proposed by Grieger et al. [25], uses an approximate formula of flux amplification during HAE, which contains a few fitting parameters. This makes possible some estimations of certain GLS characteristics, in particular, the source size [25]. For example, in case of the well-known GLS Q2237+0305 (Einstein Cross), several HAEs was observed [26–28] and the estimates of the source size have been obtained within different source models [16, 29–34].A possibility do distinguish different source models is widely discussed elsewhere. Typical problems arising in determination of the source brightness distribution are as follows.First, it is impossible to get complete information about the brightness distribution from the light curve observations. Information from separate HAE only, without making recourse to the whole light curve, and so forth, is still more limited. Observations provide only a one-dimensional luminosity profile of the source (integrated along the caustic); then, without using additional information, we cannot determine even the source size because we do not know the value and direction of the source velocity, ellipticity, and orientation with respect to the caustic.Second, even for a circularly symmetric source and known normal velocity with respect to the fold caustic, determination of the luminosity profile is a kind of ill-posed mathematical problems: small variations of input data may lead to considerable changes of the solution. A standard way to mitigate this difficulty involves additional restrictions and/or using some simple “fiducial” models for brightness distribution. Some of these models are considered below in this paper. However, it should be clearly understood that the real picture of the central quasar region is more complicated that the simple brightness distributions of the following section; these models can be considered rather as reference ones. On the other hand, in view of the present-day accuracy of observations, it is difficult to distinguish even these simple source models on the basis of observational data. For example, the authors of [35] argue that the accretion disk can be modelled with any brightness profile (Gaussian, uniform, etc.), and this model will agree with the available data provided that an appropriate source size is chosen. On the other hand, a number of authors [4, 16, 33–39] discussed delicate questions concerning determination a fine quasar structure from HAE. For example, the authors of [36] wrote that the GLITP data [28] on Q2237+0305 admit only accretion disc models (see also [4, 39]). Obviously, the presence of an accretion disk in a central region of quasar is beyond any doubts, as well as the fact that the real appearance of the quasar core can be quite different from our simplified models. However, is it possible to prove the existence of the accretion disk in a concrete GLS a posteriori? This is an open question.We note that since the work by Kochanek [37], followed by a number of authors [4, 35, 38–42], statistical methods dealing with the complete light curves of the GLS images have been developed. This approach is very attractive because it allows to take into account the whole aggregate of observational data on image brightness variations yielding estimates of the microlens masses and source model parameters. However, this treatment involves a large number of realizations of the microlensing field and requires a considerable computer time for such simulation. On the other hand, we must remember that (i) the source structure manifests itself only in the HAEs; far from the caustics the source looks like the point one and all the information about its structure is being lost; (ii) in reality, we have only one light curve, not a statistical ensemble, so probabilistic estimates of GLS parameters obtained in statistical simulations can differ from the real values in a concrete GLS. If we restrict ourselves to the HAE neighborhood, then part of the information is lost, but instead we use the most general model of the microlensing field described by a small set of coefficients in the lens mapping.


### 3.2. From a Point Source to Extended One

After these reservations, we turn to the amplification of some extended source models.LetI(y) be a surface brightness distribution of an extended source. If the source center is located at the point Y=(Y1,Y2) in the source plane, then the total microlensed flux from the source is
(35)F(Y)=∬I(y(x)-Y)dx1dx2,
where x=(x1,x2); x→y(x) is the lens mapping. The result of using (35) obviously is equivalent to the result of the well-known ray-tracing method [1] (when the pixel sizes tend to zero).An equivalent representation of this formula is(36)F(Y)=∬K(y)I(y-Y)dy1dy2,
where the point source amplification K(y)=∑iKi(y) is the sum of amplifications of all the images.Near a caustic, one can approximateK(y)=K0+Kcr(y), where K0 is an amplification of all noncritical images that is supposed to be constant during HAE, and Kcr is the amplification of the critical images. Due to a relative motion of the lensing galaxy and the source (quasar), the flux is a function of time representing the lightcurve of some quasar image in GLS.Formula (32) for Kcr contains the nonintegrable term~Θ(y2)(y2)-3/2. Therefore, the question arises of how formula (32) can be used in situation when the extended source intersects a caustic and some part of the source is in the zone between the tangent and the caustic. In view of Sections 2.2.3 and 2.4, it is evident that the mentioned term is a result of the asymptotic expansion of the root y2+κy12t2/2+⋯ in the approximate solution (18). Direct usage of the solution in the form (18) for calculation of the Jacobians of the lens mapping and then for the derivation of amplifications does not lead to any divergences and any nonintegrable terms in Kcr do not arise without using this expansion. Nevertheless, it is convenient to have a representation of Kcr in the form of expansion in powers of small parameter. Such an expansion can be fulfilled correctly after substitution of Kcr into integral (36). On this way, starting from form (18), it is easy to show that to define Kcr correctly, one must replace the term Θ(y2)(y2)-3/2 in (32) by the distribution (generalized function) (y2)+-3/2 [43]. We recall that the distribution y+-3/2 of the variable y is defined by the expression:
(37)∫y+-3/2f(y)dy=2∫0∞y-1/2∂f(y)∂ydy
for any test function f(y).After this redefinition, we have(38)Kcr=Θ(y2)2|b|y2[1+Py2+Qy1]-κ8|b|y12(y2)+-3/2.
This formula can be used to correctly derive an approximate amplification of a sufficiently smooth extended source including the case where the source crosses the caustic.


### 3.3. The Extended Source Models

#### 3.3.1.  Gaussian and Power-Law Models

Below, we list most simple and commonly used brightness distributions of a source in GLS; without loss of generality, they are chosen to be normalized to 1:(39)∬I(y)dy1dy2=1.
To compare different models of the brightness distribution, we have to use the same parameter that characterizes the size of an object. The r.m.s. size Rrms is often used:
(40)Rrms2=∬y2I(y)dy1dy2.
However, for slowly decreasing brightness profile (e.g., I(y)~|y|-α,α≤4), the r.m.s. size diverges. In case of the circularly symmetric sources the half-brightness radius R1/2 is also widely used; it is defined by the relation:
(41)2π∫0R1/2I(r)rdr=π∫0∞I(r)rdr.In case of Gaussian source model,(42)IG(r)=1πR2exp[-(rR)2],
where R stands for a size parameter; Rrms=R, R1/2=Rln(2).Limb-darkening model (see, e.g., [44]):
(43)ILD(r)=q+1πR2Ξ(rR;q),
where
(44)Ξ(ξ;q)=Θ(1-ξ2)(1-ξ2)q,
and Rrms=R/q+2. Here, we assume q>0. The half-brightness radius is R1/2=R1-(1/2)1/(q+1). For fixed Rrms and q→∞, the brightness distribution (43) tends to the Gaussian one.The models (43) and (42) describe a class of compact sources with fast brightness decrease. On the contrary, the power-law models [16, 33] describe a slow decrease at large r:
(45)IPL(r)=p-1πR2[1+r2R2]-p,
where p>1 is the power index, and R is related to the r.m.s. radius Rrms as R2=(p-2)Rrms2. The model (45) may be considered as an alternative to (43). The half-brightness radius of the source for this model is R1/2=R21/(p-1)-1. Like (43), for fixed Rrms and p→∞, the brightness distribution (45) tends to the Gaussian one. For small p, we have a “long-range” distribution; Rrms diverges for p≤2.Linear combinations of different distributions (42), (43), (45) with different parameters yield rather a wide class of symmetric source models to fit any kind of data. On the other hand, (42) may be considered as a fiducial model to determine some parameters such as the source size and (43), (45) are useful in case when in addition we are interested in investigation of a brightness behavior at large r. More physical models are considered in the next subsection.


#### 3.3.2. Accretion Disk Models

The accretion disk (AD) of Shakura-Sunyaev [45] has a more complicated profile. This model gives the energy density of the radiation from accretion disk around a nonrotating black hole as a function of radius whence the (normalized) brightness distribution is
(46)IAD(r)=3Rθ(r-R)2πr3[1-Rr],
here R being the radius of the inner edge of the accretion disk. For this AD model, the half-brightness radius is R1/2=4R, Rrms=∞. This formula describes a total brightness integrated over all radiation frequencies. The maximum brightness is at r=rm=49R/36.For a blackbody radiation, the temperature scales asT~IAD1/4, whence the specific intensity as a function of radius ~λ-5(ehc/λT-1)-1 for wavelength λ (AD1):
(47)IAD1(r)=CAD1R2[exp(κ1ρ3/4(1-ρ-1/2)1/4)-1]-1θ(r-R),ρ=rR,
here CAD1 is a normalization factor. For the maximum disk temperature (at r=rm), the peak value of intensity as a function of λ corresponds to κ1=2.422 and CAD1=9.61. Further, we adopt this value of κ1 for which Rrms/R=6.34, R1/2/R=2.021. Below in Section 3.5 we use these parameters with R/RE=0.1, so R1/2/RE≈0.2. Intensity (47) can be easily rewritten for any wavelength taking into account that κ1~λ-1.Though intensities (46), (47) are quite different, their light curves in HAE can look very similar for an appropriate choice of parameters. The main feature that distinguishes the accretion disk models from the Gaussian one is the concavity of the light curve owing to the dark “hole” in the accretion disk center. We will return to this question in Section 3.5.


### 3.4. Amplifications Including the Postlinear Corrections

#### 3.4.1.  Amplification for the Gaussian Source

Formula (38) has been used [5] to derive the amplification of a Gaussian source with the brightness distribution (42), the limb-darkening source, and the power-law source (see the next subsections).Further, we use the dimensionless coordinatess=Y1/R,h=Y2/R of the source center and the functions:
(48)Ik(h)=∫0∞uk-(1/2)exp(-u2+2uh)du=12∑n=0∞‍Γ(1/4+(k+n)/2)n!(2h)n.
These functions can be expressed in terms of the confluent hypergeometric function1F1 or the parabolic cylinder function D:
(49)Ik(h)=2-((k/2)+(1/4))Γ(k+12)e(h2/2)D-(k+1/2)(-2·h).The substitution of (38) and (42) in (36) yields
(50)KG(s,h)=12π|b|R{κπ|b|R2Φ0(h)12π|b|R+R[PΦ1(h)-κ2Φ2(h)12π|b|R+QsΦ0(h)-κs2Φ2(h)κ2]12π|b|R}.Here,(51)Φ0(h)=I0(h)exp(-h2),Φ1(h)=I1(h)exp(-h2),Φ2(h)=[hI0(h)-I1(h)]exp(-h2).
Note that the main term of (50) which corresponds to the linear caustic approximation was first obtained in the paper [20].


#### 3.4.2. Amplification for the Limb-Darkening Source

Analogous considerations allow us to obtain formulas for the amplification of extended sources for the limb-darkening (43) and power-law (45) brightness profiles; the results are represented analytically in terms of hypergeometric function2F1 [46].Denote(52)Xk,q(h)=Γ(q+2)Γ(q+3/2)∫0∞‍yk-1/2Ξ(y-h;q+12)dy,k=1,2. We have
(53)Xk,q(h)=2q+(1/2)(1+h)q+k+1Γ(q+2)Γ(k+1/2)Γ(q+k+2)×2F1(-q-12,q+32;q+k+2;1+h2)
for -1<h<1 and
(54)Xk,q(h)=π(h+1)k-(1/2)2F1(q+32,12-k;2q+3;2h+1)
for h>1.Also, fork=-1, we define
(55)X-1,q(h)=4(q+1)(hX0,q-1-X1,q-1).Then in case of the model with limb darkening (43), the critical images disappear when the source lies on the outer side of the caustic (i.e., for h<-1). The amplification due to critical images takes on the form:
(56)KLD(s,h)=12π|b|R{κ8(q+2)X0,q(h)12π|b|R+R[κ8(q+2)PX1,q(h)-κ8(q+2)X-1,q+1(h)12π|b|R+R+QsX0,q(h)-κ4s2X-1,q(h)κ8(q+2)]}.


#### 3.4.3. Amplification for a Power-Law Source

The result for the amplification involves integrals:(57)Ψk,p(h)=Γ(p-(1/2))Γ(p-1)∫0∞yk-(1/2)dy(1+(y-h)2)p-1/2=Γ(p-(1/2))Γ(p-1)B(k+12,2p-k-32)×(1+h2)(k/2)+(3/4)-p×2F1(k+12,2p-k-32;p;12(1+h1+h2))
for k=0,1, B(x,y) being the Beta-function.We extend this tok=-1 having in mind the definition of (y)+-3/2, so that
(58)Ψ-1,p(h)=4(p-1)[hΨ0,p+1(h)-Ψ1,p+1(h)].
Now, the amplification due to critical images takes on the form:
(59)KPL(s,h)=12π|b|R{κ4Ψ0,p(h)12π|b|R+R[PΨ1,p(h)-κ8(p-2)Ψ-1,p-1(h)12π|b|R+R+QsΨ0,p(h)-κ4s2Ψ-1,p(h)]}.
The zeroth approximation to this formula has been derived in [16].


### 3.5. Simulations of Microlensing Light Curves with Different Source Models

#### 3.5.1.  Equations of Microlensing and Parameters of Simulations

Though we expect that a source structure reveals itself just during HAEs, it is instructive to compare whole light curves with different source models; this is fulfilled below. Nevertheless, we point out that below we mainly deal with an estimate forηm (see  (62) below), which is just a characteristic of HAE. We assume the lens equation:
(60)y=x-∑n=1NRE,n2x-xn|x-xn|2+an2,
which describes microlensing by N extended masses (clumps) in absence of any external field; the surface mass density corresponding to (60) is
(61)σ(x)=∑n=1Nf(an,Mn,x-xn),
where the surface mass density of one clump is f(a,M,x)~a2M(x2+a2)-2; here an,Mn stand for the size and mass of nth clump. In case of point masses an=0; however, in simulations of the point mass microlensing, we have chosen some small positive values of an≪RE,n so as not to have problems with singularities. In Section 4 we also consider the case of a finite-size extended masses.Here, we present results of straightforward calculations of a microlensed flux according to (35) to obtain the light curves for different realizations of the point microlenses positions. The parameters of numerical integration are presented in Table 1 along with the parameters of the microlensing field. The microlense positions were chosen in a random way with uniform distribution over the field. The length of trajectory has been taken sufficiently long so as to provide the caustic crossings. However, the trajectory in the simulations was situated far from the field boundaries, and the size of the field was chosen large enough to avoid boundary effects.Table 1
Parameters of simulation.
Parameter
Value
Number of pixels
1.23
×
1
0
6
Pixel size
0.01RE
Source trajectory length
2RE; 5 RE
Radius of field
70RE
Microlensing optical depth (σ)
0.3
Source speed (V)
1
Time discretization (δt)
0.01RE/VHere we compare the light curves for Gaussian, PL, LD, AD, and AD1 source models with brightness distributions from Section3.3 The simulations were performed for a set of 100 random realizations of microlensing field with optical depth σ=0.3 in order to have possibility to compare with some of the results of papers [47–50]. All source models had the same half-brightness radius, the light curves for all source models have been calculated for the same microlensing fields. All the microlenses were static, the total number of microlenses was 1470. The typical magnification pattern is shown on Figure 1; here the speed of source is V=1 so we can identify the source position as a function of time.Figure 1
Magnification pattern with the trajectory of the source in the source plane. The source moves from left to right along the straight line with uniform speed.The simulations were carried out with half-brightness radiusR1/2=0.21 for all models; the power-law index was p=3/2 for the “long range” PL model; note that AD model also corresponds to this class of the power-law asymptotic dependence (with p=3/2). For LD model, we have chosen q=1.To compare the light curves for different source models1 and 2, we used the relative value:
(62)ηm=2〈maxt{η(t)}〉,η(t)=|K1(t)-K2(t)|K1(t)+K2(t),
where Ki(t) is the amplification for ith model along the trajectory of source linear motion; 〈⋯〉 is an average over a number of realizations. Because we expect that maximum difference of brightness on lightcurves for different source models takes place during HAE, so ηm typically is a characteristic of HAE.From the “light curves” on Figure2 we observe a significant difference between the “compact” (LD and Gaussian) and the “long-range” models. The long-range character of the latter reveals itself even on considerable distances from the caustics, where we expect that brightness of all the sources must have the same behavior as that of a point source. The differences between these two groups of models are essentially larger than the differences within each group (e.g., between Gaussian and LD models). This conclusion is confirmed by results of the statistical considerations over 100 realizations shown in Table 2 for half-brightness radius R1/2=0.21 as the example.Table 2
Relative difference between the light curves in HAE.
ith model
jth model
η
m
AD
Gaussian
0
.
074
±
0
.
001
AD
AD1
0
.
085
±
0
.
002
AD
LD
0
.
091
±
0
.
002
AD
PL
0
.
038
±
0
.
001
Gaussian
AD1
0
.
073
±
0
.
002
Gaussian
LD
0
.
042
±
0
.
001
Gaussian
PL
0
.
073
±
0
.
001
AD1
LD
0
.
094
±
0
.
002
AD1
PL
0
.
052
±
0
.
001
LD
PL
0
.
090
±
0
.
002Figure 2
“Light curves”: amplification as a function of time for different source models that correspond to magnification pattern of Figure1.We note that the results of comparison may depend on a complexity of the caustics involved into our consideration. One can expect, for example, that in case of complicated caustic crossing (such as crossing of the fold caustics in a vicinity of cusp, intersection of dense aggregations of caustics, among others) the light curve difference for different source models may be more significant. Our previous simulations involvedall possible realizations of the microlensing field that may invoke the complicated caustic crossings. However, in case of a concrete GLS, we deal with a single light curve. One may ask: is it relevant to apply the results of statistical simulations to this single light curve, where we may have a kind of an observational selection. Therefore, it is necessary to check how this “complexity” affects the average value of η. In this connection, we considered some modification of our statistical procedure: we have chosen (by eye) such realizations rather simple fold caustic crossings involved. However, the numerical results for simple caustic crossing events, such as shown as the example on Figure 3, appeared to be nearly the same as that of the Table 2. Though we must note that here the statistics have been worse because of a smaller number of the “simple” caustic realizations. So this conclusion requires further verification.Figure 3
Magnification pattern with simple caustic crossing events.


#### 3.5.2. Gaussian Fittings of Accretion Disk and Limb Darkening Models

Table2 concerns a comparison of different models with the sameR1/2. However, in reality, we do not know what source we must fit and one must check whether we can replace one model with a different one with some other source parameters to get better fitting.We have fitted the limb-darkening and accretion disk model light curves with that of the Gaussian source of different radii. The half-brightness radius of LD and AD models was fixed:R1/2=0.2RE (this choice was determined by the AD1 model, see comments after formula (47)) and that of the Gaussian source was varied to get better fitting; here P=3/2 for PL model and q=1 for LD model. One example of the relative light curve differences η(T) according to (62) between LD, AD1 and the best fitted Gaussian models is shown in Figure 4. As we see from this figure, usage of an optimal size of the Gaussian model can diminish ηm as compared to the average estimate from Table 2 dealing with the same R1/2 for all the models.Figure 4
Relative differencesη(t) between Gaussian and LD (q=1), AD1 models for the best fitted curves corresponding to magnification pattern of Figure 3.To test this statement statistically, we carried out simulations with 20 magnification patterns having the “simple” caustic crossing. We found that fitting of AD light curves by the Gaussian model (with differentR1/2) yields ηm=0.07±0.006; this is roughly the same as in Table 2. Fittings of LD model light curves by the Gaussian model (with different R1/2) yielding ηm=0.026±0.002 (less then half of Table 2); in this case, an improvement is noticeable. Therefore, the Gaussian source cannot reproduce all the models, though for some compact models the fitting results are rather good. Then, in the caustic crossing events different models may be distinguished on the accuracy level of modern photometric observations.


### 3.6. Resume of Section3

Here the formulas of Section2 are applied to the Gaussian, power-law and limb darkening models of an extended source. The analytical formulas for these source models are obtained. We note that a fitting of the light curve of GLS Q2237+0305C [5, 6] showed that some of these higher-order corrections can be statistically significant even at the present accuracy level. This means that when we are looking for some fine effects in HAE due to the source size, a consistent treatment must also take into account possible effects of the above corrections.In case of the accretion disk models, we have not succeeded to find compact “workable” expressions forKcr, so this case has been treated numerically to compare possible observational differences of various source models.Now we remind the main assumptions and the class of models used, which is necessary for understanding the results of the Section3.5. The simulations use the GLS parameters similar to that of Q2237+030 [47, 51]. We consider the equal mass microlenses; we do not consider any mass distributions and/or populations of small (planetary) masses. Next, we consider most simple circular symmetric source models without effects of ellipticity, and so forth. The results of Section 3.5 can be summed up as follows. There are special classes of models: (i) accretion disk models having the dark region in the center and (ii) source models with a slow decrease of the (long-range) brightness distribution at large distances from the center. From our results, we see that (i) and (ii) show in HAEs a behavior which is different from that of the Gaussian and the limb-darkening models on a relative accuracy level of~4÷7%. There is then a possibility to distinguish such models on account of observations during HAEs. Outside HAE the surface brightness profile typically has little effect on microlensing (see [35, 52]).We note that our conclusions involve concrete source models with a special choice of parameters (source sizes, microlensing optical depth, etc.) and a considerable work is to be fulfilled to obtain more general results. Also, our results are of a statistical nature, so our average estimate ofηm may be different from analogous one in a separate light curve observation.Note that adding a sufficient number of microlenses having a small (e.g., planetary) masses (cf. [51, 53–55]) must produce an additional caustics that in its turn can give rise to fine features of light curves during HAE. One may expect that this can mimic features due to brightness distribution over the source. This question deserves further investigation.


## 3.1.   Preliminary Comments

In this section, we proceed to applications of the previous results to some extended source models typically used in fittings of the observed light curves. Interest to HAE in extragalactic GLSs is due to possibilities to study a brightness distributions in sources. This is especially interesting in connection with investigation of central regions of distant quasars. The idea, first proposed by Grieger et al. [25], uses an approximate formula of flux amplification during HAE, which contains a few fitting parameters. This makes possible some estimations of certain GLS characteristics, in particular, the source size [25]. For example, in case of the well-known GLS Q2237+0305 (Einstein Cross), several HAEs was observed [26–28] and the estimates of the source size have been obtained within different source models [16, 29–34].A possibility do distinguish different source models is widely discussed elsewhere. Typical problems arising in determination of the source brightness distribution are as follows.First, it is impossible to get complete information about the brightness distribution from the light curve observations. Information from separate HAE only, without making recourse to the whole light curve, and so forth, is still more limited. Observations provide only a one-dimensional luminosity profile of the source (integrated along the caustic); then, without using additional information, we cannot determine even the source size because we do not know the value and direction of the source velocity, ellipticity, and orientation with respect to the caustic.Second, even for a circularly symmetric source and known normal velocity with respect to the fold caustic, determination of the luminosity profile is a kind of ill-posed mathematical problems: small variations of input data may lead to considerable changes of the solution. A standard way to mitigate this difficulty involves additional restrictions and/or using some simple “fiducial” models for brightness distribution. Some of these models are considered below in this paper. However, it should be clearly understood that the real picture of the central quasar region is more complicated that the simple brightness distributions of the following section; these models can be considered rather as reference ones. On the other hand, in view of the present-day accuracy of observations, it is difficult to distinguish even these simple source models on the basis of observational data. For example, the authors of [35] argue that the accretion disk can be modelled with any brightness profile (Gaussian, uniform, etc.), and this model will agree with the available data provided that an appropriate source size is chosen. On the other hand, a number of authors [4, 16, 33–39] discussed delicate questions concerning determination a fine quasar structure from HAE. For example, the authors of [36] wrote that the GLITP data [28] on Q2237+0305 admit only accretion disc models (see also [4, 39]). Obviously, the presence of an accretion disk in a central region of quasar is beyond any doubts, as well as the fact that the real appearance of the quasar core can be quite different from our simplified models. However, is it possible to prove the existence of the accretion disk in a concrete GLS a posteriori? This is an open question.We note that since the work by Kochanek [37], followed by a number of authors [4, 35, 38–42], statistical methods dealing with the complete light curves of the GLS images have been developed. This approach is very attractive because it allows to take into account the whole aggregate of observational data on image brightness variations yielding estimates of the microlens masses and source model parameters. However, this treatment involves a large number of realizations of the microlensing field and requires a considerable computer time for such simulation. On the other hand, we must remember that (i) the source structure manifests itself only in the HAEs; far from the caustics the source looks like the point one and all the information about its structure is being lost; (ii) in reality, we have only one light curve, not a statistical ensemble, so probabilistic estimates of GLS parameters obtained in statistical simulations can differ from the real values in a concrete GLS. If we restrict ourselves to the HAE neighborhood, then part of the information is lost, but instead we use the most general model of the microlensing field described by a small set of coefficients in the lens mapping.


## 3.2. From a Point Source to Extended One

After these reservations, we turn to the amplification of some extended source models.LetI(y) be a surface brightness distribution of an extended source. If the source center is located at the point Y=(Y1,Y2) in the source plane, then the total microlensed flux from the source is
(35)F(Y)=∬I(y(x)-Y)dx1dx2,
where x=(x1,x2); x→y(x) is the lens mapping. The result of using (35) obviously is equivalent to the result of the well-known ray-tracing method [1] (when the pixel sizes tend to zero).An equivalent representation of this formula is(36)F(Y)=∬K(y)I(y-Y)dy1dy2,
where the point source amplification K(y)=∑iKi(y) is the sum of amplifications of all the images.Near a caustic, one can approximateK(y)=K0+Kcr(y), where K0 is an amplification of all noncritical images that is supposed to be constant during HAE, and Kcr is the amplification of the critical images. Due to a relative motion of the lensing galaxy and the source (quasar), the flux is a function of time representing the lightcurve of some quasar image in GLS.Formula (32) for Kcr contains the nonintegrable term~Θ(y2)(y2)-3/2. Therefore, the question arises of how formula (32) can be used in situation when the extended source intersects a caustic and some part of the source is in the zone between the tangent and the caustic. In view of Sections 2.2.3 and 2.4, it is evident that the mentioned term is a result of the asymptotic expansion of the root y2+κy12t2/2+⋯ in the approximate solution (18). Direct usage of the solution in the form (18) for calculation of the Jacobians of the lens mapping and then for the derivation of amplifications does not lead to any divergences and any nonintegrable terms in Kcr do not arise without using this expansion. Nevertheless, it is convenient to have a representation of Kcr in the form of expansion in powers of small parameter. Such an expansion can be fulfilled correctly after substitution of Kcr into integral (36). On this way, starting from form (18), it is easy to show that to define Kcr correctly, one must replace the term Θ(y2)(y2)-3/2 in (32) by the distribution (generalized function) (y2)+-3/2 [43]. We recall that the distribution y+-3/2 of the variable y is defined by the expression:
(37)∫y+-3/2f(y)dy=2∫0∞y-1/2∂f(y)∂ydy
for any test function f(y).After this redefinition, we have(38)Kcr=Θ(y2)2|b|y2[1+Py2+Qy1]-κ8|b|y12(y2)+-3/2.
This formula can be used to correctly derive an approximate amplification of a sufficiently smooth extended source including the case where the source crosses the caustic.


## 3.3. The Extended Source Models

### 3.3.1.  Gaussian and Power-Law Models

Below, we list most simple and commonly used brightness distributions of a source in GLS; without loss of generality, they are chosen to be normalized to 1:(39)∬I(y)dy1dy2=1.
To compare different models of the brightness distribution, we have to use the same parameter that characterizes the size of an object. The r.m.s. size Rrms is often used:
(40)Rrms2=∬y2I(y)dy1dy2.
However, for slowly decreasing brightness profile (e.g., I(y)~|y|-α,α≤4), the r.m.s. size diverges. In case of the circularly symmetric sources the half-brightness radius R1/2 is also widely used; it is defined by the relation:
(41)2π∫0R1/2I(r)rdr=π∫0∞I(r)rdr.In case of Gaussian source model,(42)IG(r)=1πR2exp[-(rR)2],
where R stands for a size parameter; Rrms=R, R1/2=Rln(2).Limb-darkening model (see, e.g., [44]):
(43)ILD(r)=q+1πR2Ξ(rR;q),
where
(44)Ξ(ξ;q)=Θ(1-ξ2)(1-ξ2)q,
and Rrms=R/q+2. Here, we assume q>0. The half-brightness radius is R1/2=R1-(1/2)1/(q+1). For fixed Rrms and q→∞, the brightness distribution (43) tends to the Gaussian one.The models (43) and (42) describe a class of compact sources with fast brightness decrease. On the contrary, the power-law models [16, 33] describe a slow decrease at large r:
(45)IPL(r)=p-1πR2[1+r2R2]-p,
where p>1 is the power index, and R is related to the r.m.s. radius Rrms as R2=(p-2)Rrms2. The model (45) may be considered as an alternative to (43). The half-brightness radius of the source for this model is R1/2=R21/(p-1)-1. Like (43), for fixed Rrms and p→∞, the brightness distribution (45) tends to the Gaussian one. For small p, we have a “long-range” distribution; Rrms diverges for p≤2.Linear combinations of different distributions (42), (43), (45) with different parameters yield rather a wide class of symmetric source models to fit any kind of data. On the other hand, (42) may be considered as a fiducial model to determine some parameters such as the source size and (43), (45) are useful in case when in addition we are interested in investigation of a brightness behavior at large r. More physical models are considered in the next subsection.


### 3.3.2. Accretion Disk Models

The accretion disk (AD) of Shakura-Sunyaev [45] has a more complicated profile. This model gives the energy density of the radiation from accretion disk around a nonrotating black hole as a function of radius whence the (normalized) brightness distribution is
(46)IAD(r)=3Rθ(r-R)2πr3[1-Rr],
here R being the radius of the inner edge of the accretion disk. For this AD model, the half-brightness radius is R1/2=4R, Rrms=∞. This formula describes a total brightness integrated over all radiation frequencies. The maximum brightness is at r=rm=49R/36.For a blackbody radiation, the temperature scales asT~IAD1/4, whence the specific intensity as a function of radius ~λ-5(ehc/λT-1)-1 for wavelength λ (AD1):
(47)IAD1(r)=CAD1R2[exp(κ1ρ3/4(1-ρ-1/2)1/4)-1]-1θ(r-R),ρ=rR,
here CAD1 is a normalization factor. For the maximum disk temperature (at r=rm), the peak value of intensity as a function of λ corresponds to κ1=2.422 and CAD1=9.61. Further, we adopt this value of κ1 for which Rrms/R=6.34, R1/2/R=2.021. Below in Section 3.5 we use these parameters with R/RE=0.1, so R1/2/RE≈0.2. Intensity (47) can be easily rewritten for any wavelength taking into account that κ1~λ-1.Though intensities (46), (47) are quite different, their light curves in HAE can look very similar for an appropriate choice of parameters. The main feature that distinguishes the accretion disk models from the Gaussian one is the concavity of the light curve owing to the dark “hole” in the accretion disk center. We will return to this question in Section 3.5.


## 3.3.1.  Gaussian and Power-Law Models

Below, we list most simple and commonly used brightness distributions of a source in GLS; without loss of generality, they are chosen to be normalized to 1:(39)∬I(y)dy1dy2=1.
To compare different models of the brightness distribution, we have to use the same parameter that characterizes the size of an object. The r.m.s. size Rrms is often used:
(40)Rrms2=∬y2I(y)dy1dy2.
However, for slowly decreasing brightness profile (e.g., I(y)~|y|-α,α≤4), the r.m.s. size diverges. In case of the circularly symmetric sources the half-brightness radius R1/2 is also widely used; it is defined by the relation:
(41)2π∫0R1/2I(r)rdr=π∫0∞I(r)rdr.In case of Gaussian source model,(42)IG(r)=1πR2exp[-(rR)2],
where R stands for a size parameter; Rrms=R, R1/2=Rln(2).Limb-darkening model (see, e.g., [44]):
(43)ILD(r)=q+1πR2Ξ(rR;q),
where
(44)Ξ(ξ;q)=Θ(1-ξ2)(1-ξ2)q,
and Rrms=R/q+2. Here, we assume q>0. The half-brightness radius is R1/2=R1-(1/2)1/(q+1). For fixed Rrms and q→∞, the brightness distribution (43) tends to the Gaussian one.The models (43) and (42) describe a class of compact sources with fast brightness decrease. On the contrary, the power-law models [16, 33] describe a slow decrease at large r:
(45)IPL(r)=p-1πR2[1+r2R2]-p,
where p>1 is the power index, and R is related to the r.m.s. radius Rrms as R2=(p-2)Rrms2. The model (45) may be considered as an alternative to (43). The half-brightness radius of the source for this model is R1/2=R21/(p-1)-1. Like (43), for fixed Rrms and p→∞, the brightness distribution (45) tends to the Gaussian one. For small p, we have a “long-range” distribution; Rrms diverges for p≤2.Linear combinations of different distributions (42), (43), (45) with different parameters yield rather a wide class of symmetric source models to fit any kind of data. On the other hand, (42) may be considered as a fiducial model to determine some parameters such as the source size and (43), (45) are useful in case when in addition we are interested in investigation of a brightness behavior at large r. More physical models are considered in the next subsection.


## 3.3.2. Accretion Disk Models

The accretion disk (AD) of Shakura-Sunyaev [45] has a more complicated profile. This model gives the energy density of the radiation from accretion disk around a nonrotating black hole as a function of radius whence the (normalized) brightness distribution is
(46)IAD(r)=3Rθ(r-R)2πr3[1-Rr],
here R being the radius of the inner edge of the accretion disk. For this AD model, the half-brightness radius is R1/2=4R, Rrms=∞. This formula describes a total brightness integrated over all radiation frequencies. The maximum brightness is at r=rm=49R/36.For a blackbody radiation, the temperature scales asT~IAD1/4, whence the specific intensity as a function of radius ~λ-5(ehc/λT-1)-1 for wavelength λ (AD1):
(47)IAD1(r)=CAD1R2[exp(κ1ρ3/4(1-ρ-1/2)1/4)-1]-1θ(r-R),ρ=rR,
here CAD1 is a normalization factor. For the maximum disk temperature (at r=rm), the peak value of intensity as a function of λ corresponds to κ1=2.422 and CAD1=9.61. Further, we adopt this value of κ1 for which Rrms/R=6.34, R1/2/R=2.021. Below in Section 3.5 we use these parameters with R/RE=0.1, so R1/2/RE≈0.2. Intensity (47) can be easily rewritten for any wavelength taking into account that κ1~λ-1.Though intensities (46), (47) are quite different, their light curves in HAE can look very similar for an appropriate choice of parameters. The main feature that distinguishes the accretion disk models from the Gaussian one is the concavity of the light curve owing to the dark “hole” in the accretion disk center. We will return to this question in Section 3.5.


## 3.4. Amplifications Including the Postlinear Corrections

### 3.4.1.  Amplification for the Gaussian Source

Formula (38) has been used [5] to derive the amplification of a Gaussian source with the brightness distribution (42), the limb-darkening source, and the power-law source (see the next subsections).Further, we use the dimensionless coordinatess=Y1/R,h=Y2/R of the source center and the functions:
(48)Ik(h)=∫0∞uk-(1/2)exp(-u2+2uh)du=12∑n=0∞‍Γ(1/4+(k+n)/2)n!(2h)n.
These functions can be expressed in terms of the confluent hypergeometric function1F1 or the parabolic cylinder function D:
(49)Ik(h)=2-((k/2)+(1/4))Γ(k+12)e(h2/2)D-(k+1/2)(-2·h).The substitution of (38) and (42) in (36) yields
(50)KG(s,h)=12π|b|R{κπ|b|R2Φ0(h)12π|b|R+R[PΦ1(h)-κ2Φ2(h)12π|b|R+QsΦ0(h)-κs2Φ2(h)κ2]12π|b|R}.Here,(51)Φ0(h)=I0(h)exp(-h2),Φ1(h)=I1(h)exp(-h2),Φ2(h)=[hI0(h)-I1(h)]exp(-h2).
Note that the main term of (50) which corresponds to the linear caustic approximation was first obtained in the paper [20].


### 3.4.2. Amplification for the Limb-Darkening Source

Analogous considerations allow us to obtain formulas for the amplification of extended sources for the limb-darkening (43) and power-law (45) brightness profiles; the results are represented analytically in terms of hypergeometric function2F1 [46].Denote(52)Xk,q(h)=Γ(q+2)Γ(q+3/2)∫0∞‍yk-1/2Ξ(y-h;q+12)dy,k=1,2. We have
(53)Xk,q(h)=2q+(1/2)(1+h)q+k+1Γ(q+2)Γ(k+1/2)Γ(q+k+2)×2F1(-q-12,q+32;q+k+2;1+h2)
for -1<h<1 and
(54)Xk,q(h)=π(h+1)k-(1/2)2F1(q+32,12-k;2q+3;2h+1)
for h>1.Also, fork=-1, we define
(55)X-1,q(h)=4(q+1)(hX0,q-1-X1,q-1).Then in case of the model with limb darkening (43), the critical images disappear when the source lies on the outer side of the caustic (i.e., for h<-1). The amplification due to critical images takes on the form:
(56)KLD(s,h)=12π|b|R{κ8(q+2)X0,q(h)12π|b|R+R[κ8(q+2)PX1,q(h)-κ8(q+2)X-1,q+1(h)12π|b|R+R+QsX0,q(h)-κ4s2X-1,q(h)κ8(q+2)]}.


### 3.4.3. Amplification for a Power-Law Source

The result for the amplification involves integrals:(57)Ψk,p(h)=Γ(p-(1/2))Γ(p-1)∫0∞yk-(1/2)dy(1+(y-h)2)p-1/2=Γ(p-(1/2))Γ(p-1)B(k+12,2p-k-32)×(1+h2)(k/2)+(3/4)-p×2F1(k+12,2p-k-32;p;12(1+h1+h2))
for k=0,1, B(x,y) being the Beta-function.We extend this tok=-1 having in mind the definition of (y)+-3/2, so that
(58)Ψ-1,p(h)=4(p-1)[hΨ0,p+1(h)-Ψ1,p+1(h)].
Now, the amplification due to critical images takes on the form:
(59)KPL(s,h)=12π|b|R{κ4Ψ0,p(h)12π|b|R+R[PΨ1,p(h)-κ8(p-2)Ψ-1,p-1(h)12π|b|R+R+QsΨ0,p(h)-κ4s2Ψ-1,p(h)]}.
The zeroth approximation to this formula has been derived in [16].


## 3.4.1.  Amplification for the Gaussian Source

Formula (38) has been used [5] to derive the amplification of a Gaussian source with the brightness distribution (42), the limb-darkening source, and the power-law source (see the next subsections).Further, we use the dimensionless coordinatess=Y1/R,h=Y2/R of the source center and the functions:
(48)Ik(h)=∫0∞uk-(1/2)exp(-u2+2uh)du=12∑n=0∞‍Γ(1/4+(k+n)/2)n!(2h)n.
These functions can be expressed in terms of the confluent hypergeometric function1F1 or the parabolic cylinder function D:
(49)Ik(h)=2-((k/2)+(1/4))Γ(k+12)e(h2/2)D-(k+1/2)(-2·h).The substitution of (38) and (42) in (36) yields
(50)KG(s,h)=12π|b|R{κπ|b|R2Φ0(h)12π|b|R+R[PΦ1(h)-κ2Φ2(h)12π|b|R+QsΦ0(h)-κs2Φ2(h)κ2]12π|b|R}.Here,(51)Φ0(h)=I0(h)exp(-h2),Φ1(h)=I1(h)exp(-h2),Φ2(h)=[hI0(h)-I1(h)]exp(-h2).
Note that the main term of (50) which corresponds to the linear caustic approximation was first obtained in the paper [20].


## 3.4.2. Amplification for the Limb-Darkening Source

Analogous considerations allow us to obtain formulas for the amplification of extended sources for the limb-darkening (43) and power-law (45) brightness profiles; the results are represented analytically in terms of hypergeometric function2F1 [46].Denote(52)Xk,q(h)=Γ(q+2)Γ(q+3/2)∫0∞‍yk-1/2Ξ(y-h;q+12)dy,k=1,2. We have
(53)Xk,q(h)=2q+(1/2)(1+h)q+k+1Γ(q+2)Γ(k+1/2)Γ(q+k+2)×2F1(-q-12,q+32;q+k+2;1+h2)
for -1<h<1 and
(54)Xk,q(h)=π(h+1)k-(1/2)2F1(q+32,12-k;2q+3;2h+1)
for h>1.Also, fork=-1, we define
(55)X-1,q(h)=4(q+1)(hX0,q-1-X1,q-1).Then in case of the model with limb darkening (43), the critical images disappear when the source lies on the outer side of the caustic (i.e., for h<-1). The amplification due to critical images takes on the form:
(56)KLD(s,h)=12π|b|R{κ8(q+2)X0,q(h)12π|b|R+R[κ8(q+2)PX1,q(h)-κ8(q+2)X-1,q+1(h)12π|b|R+R+QsX0,q(h)-κ4s2X-1,q(h)κ8(q+2)]}.


## 3.4.3. Amplification for a Power-Law Source

The result for the amplification involves integrals:(57)Ψk,p(h)=Γ(p-(1/2))Γ(p-1)∫0∞yk-(1/2)dy(1+(y-h)2)p-1/2=Γ(p-(1/2))Γ(p-1)B(k+12,2p-k-32)×(1+h2)(k/2)+(3/4)-p×2F1(k+12,2p-k-32;p;12(1+h1+h2))
for k=0,1, B(x,y) being the Beta-function.We extend this tok=-1 having in mind the definition of (y)+-3/2, so that
(58)Ψ-1,p(h)=4(p-1)[hΨ0,p+1(h)-Ψ1,p+1(h)].
Now, the amplification due to critical images takes on the form:
(59)KPL(s,h)=12π|b|R{κ4Ψ0,p(h)12π|b|R+R[PΨ1,p(h)-κ8(p-2)Ψ-1,p-1(h)12π|b|R+R+QsΨ0,p(h)-κ4s2Ψ-1,p(h)]}.
The zeroth approximation to this formula has been derived in [16].


## 3.5. Simulations of Microlensing Light Curves with Different Source Models

### 3.5.1.  Equations of Microlensing and Parameters of Simulations

Though we expect that a source structure reveals itself just during HAEs, it is instructive to compare whole light curves with different source models; this is fulfilled below. Nevertheless, we point out that below we mainly deal with an estimate forηm (see  (62) below), which is just a characteristic of HAE. We assume the lens equation:
(60)y=x-∑n=1NRE,n2x-xn|x-xn|2+an2,
which describes microlensing by N extended masses (clumps) in absence of any external field; the surface mass density corresponding to (60) is
(61)σ(x)=∑n=1Nf(an,Mn,x-xn),
where the surface mass density of one clump is f(a,M,x)~a2M(x2+a2)-2; here an,Mn stand for the size and mass of nth clump. In case of point masses an=0; however, in simulations of the point mass microlensing, we have chosen some small positive values of an≪RE,n so as not to have problems with singularities. In Section 4 we also consider the case of a finite-size extended masses.Here, we present results of straightforward calculations of a microlensed flux according to (35) to obtain the light curves for different realizations of the point microlenses positions. The parameters of numerical integration are presented in Table 1 along with the parameters of the microlensing field. The microlense positions were chosen in a random way with uniform distribution over the field. The length of trajectory has been taken sufficiently long so as to provide the caustic crossings. However, the trajectory in the simulations was situated far from the field boundaries, and the size of the field was chosen large enough to avoid boundary effects.Table 1
Parameters of simulation.
Parameter
Value
Number of pixels
1.23
×
1
0
6
Pixel size
0.01RE
Source trajectory length
2RE; 5 RE
Radius of field
70RE
Microlensing optical depth (σ)
0.3
Source speed (V)
1
Time discretization (δt)
0.01RE/VHere we compare the light curves for Gaussian, PL, LD, AD, and AD1 source models with brightness distributions from Section3.3 The simulations were performed for a set of 100 random realizations of microlensing field with optical depth σ=0.3 in order to have possibility to compare with some of the results of papers [47–50]. All source models had the same half-brightness radius, the light curves for all source models have been calculated for the same microlensing fields. All the microlenses were static, the total number of microlenses was 1470. The typical magnification pattern is shown on Figure 1; here the speed of source is V=1 so we can identify the source position as a function of time.Figure 1
Magnification pattern with the trajectory of the source in the source plane. The source moves from left to right along the straight line with uniform speed.The simulations were carried out with half-brightness radiusR1/2=0.21 for all models; the power-law index was p=3/2 for the “long range” PL model; note that AD model also corresponds to this class of the power-law asymptotic dependence (with p=3/2). For LD model, we have chosen q=1.To compare the light curves for different source models1 and 2, we used the relative value:
(62)ηm=2〈maxt{η(t)}〉,η(t)=|K1(t)-K2(t)|K1(t)+K2(t),
where Ki(t) is the amplification for ith model along the trajectory of source linear motion; 〈⋯〉 is an average over a number of realizations. Because we expect that maximum difference of brightness on lightcurves for different source models takes place during HAE, so ηm typically is a characteristic of HAE.From the “light curves” on Figure2 we observe a significant difference between the “compact” (LD and Gaussian) and the “long-range” models. The long-range character of the latter reveals itself even on considerable distances from the caustics, where we expect that brightness of all the sources must have the same behavior as that of a point source. The differences between these two groups of models are essentially larger than the differences within each group (e.g., between Gaussian and LD models). This conclusion is confirmed by results of the statistical considerations over 100 realizations shown in Table 2 for half-brightness radius R1/2=0.21 as the example.Table 2
Relative difference between the light curves in HAE.
ith model
jth model
η
m
AD
Gaussian
0
.
074
±
0
.
001
AD
AD1
0
.
085
±
0
.
002
AD
LD
0
.
091
±
0
.
002
AD
PL
0
.
038
±
0
.
001
Gaussian
AD1
0
.
073
±
0
.
002
Gaussian
LD
0
.
042
±
0
.
001
Gaussian
PL
0
.
073
±
0
.
001
AD1
LD
0
.
094
±
0
.
002
AD1
PL
0
.
052
±
0
.
001
LD
PL
0
.
090
±
0
.
002Figure 2
“Light curves”: amplification as a function of time for different source models that correspond to magnification pattern of Figure1.We note that the results of comparison may depend on a complexity of the caustics involved into our consideration. One can expect, for example, that in case of complicated caustic crossing (such as crossing of the fold caustics in a vicinity of cusp, intersection of dense aggregations of caustics, among others) the light curve difference for different source models may be more significant. Our previous simulations involvedall possible realizations of the microlensing field that may invoke the complicated caustic crossings. However, in case of a concrete GLS, we deal with a single light curve. One may ask: is it relevant to apply the results of statistical simulations to this single light curve, where we may have a kind of an observational selection. Therefore, it is necessary to check how this “complexity” affects the average value of η. In this connection, we considered some modification of our statistical procedure: we have chosen (by eye) such realizations rather simple fold caustic crossings involved. However, the numerical results for simple caustic crossing events, such as shown as the example on Figure 3, appeared to be nearly the same as that of the Table 2. Though we must note that here the statistics have been worse because of a smaller number of the “simple” caustic realizations. So this conclusion requires further verification.Figure 3
Magnification pattern with simple caustic crossing events.


### 3.5.2. Gaussian Fittings of Accretion Disk and Limb Darkening Models

Table2 concerns a comparison of different models with the sameR1/2. However, in reality, we do not know what source we must fit and one must check whether we can replace one model with a different one with some other source parameters to get better fitting.We have fitted the limb-darkening and accretion disk model light curves with that of the Gaussian source of different radii. The half-brightness radius of LD and AD models was fixed:R1/2=0.2RE (this choice was determined by the AD1 model, see comments after formula (47)) and that of the Gaussian source was varied to get better fitting; here P=3/2 for PL model and q=1 for LD model. One example of the relative light curve differences η(T) according to (62) between LD, AD1 and the best fitted Gaussian models is shown in Figure 4. As we see from this figure, usage of an optimal size of the Gaussian model can diminish ηm as compared to the average estimate from Table 2 dealing with the same R1/2 for all the models.Figure 4
Relative differencesη(t) between Gaussian and LD (q=1), AD1 models for the best fitted curves corresponding to magnification pattern of Figure 3.To test this statement statistically, we carried out simulations with 20 magnification patterns having the “simple” caustic crossing. We found that fitting of AD light curves by the Gaussian model (with differentR1/2) yields ηm=0.07±0.006; this is roughly the same as in Table 2. Fittings of LD model light curves by the Gaussian model (with different R1/2) yielding ηm=0.026±0.002 (less then half of Table 2); in this case, an improvement is noticeable. Therefore, the Gaussian source cannot reproduce all the models, though for some compact models the fitting results are rather good. Then, in the caustic crossing events different models may be distinguished on the accuracy level of modern photometric observations.


## 3.5.1.  Equations of Microlensing and Parameters of Simulations

Though we expect that a source structure reveals itself just during HAEs, it is instructive to compare whole light curves with different source models; this is fulfilled below. Nevertheless, we point out that below we mainly deal with an estimate forηm (see  (62) below), which is just a characteristic of HAE. We assume the lens equation:
(60)y=x-∑n=1NRE,n2x-xn|x-xn|2+an2,
which describes microlensing by N extended masses (clumps) in absence of any external field; the surface mass density corresponding to (60) is
(61)σ(x)=∑n=1Nf(an,Mn,x-xn),
where the surface mass density of one clump is f(a,M,x)~a2M(x2+a2)-2; here an,Mn stand for the size and mass of nth clump. In case of point masses an=0; however, in simulations of the point mass microlensing, we have chosen some small positive values of an≪RE,n so as not to have problems with singularities. In Section 4 we also consider the case of a finite-size extended masses.Here, we present results of straightforward calculations of a microlensed flux according to (35) to obtain the light curves for different realizations of the point microlenses positions. The parameters of numerical integration are presented in Table 1 along with the parameters of the microlensing field. The microlense positions were chosen in a random way with uniform distribution over the field. The length of trajectory has been taken sufficiently long so as to provide the caustic crossings. However, the trajectory in the simulations was situated far from the field boundaries, and the size of the field was chosen large enough to avoid boundary effects.Table 1
Parameters of simulation.
Parameter
Value
Number of pixels
1.23
×
1
0
6
Pixel size
0.01RE
Source trajectory length
2RE; 5 RE
Radius of field
70RE
Microlensing optical depth (σ)
0.3
Source speed (V)
1
Time discretization (δt)
0.01RE/VHere we compare the light curves for Gaussian, PL, LD, AD, and AD1 source models with brightness distributions from Section3.3 The simulations were performed for a set of 100 random realizations of microlensing field with optical depth σ=0.3 in order to have possibility to compare with some of the results of papers [47–50]. All source models had the same half-brightness radius, the light curves for all source models have been calculated for the same microlensing fields. All the microlenses were static, the total number of microlenses was 1470. The typical magnification pattern is shown on Figure 1; here the speed of source is V=1 so we can identify the source position as a function of time.Figure 1
Magnification pattern with the trajectory of the source in the source plane. The source moves from left to right along the straight line with uniform speed.The simulations were carried out with half-brightness radiusR1/2=0.21 for all models; the power-law index was p=3/2 for the “long range” PL model; note that AD model also corresponds to this class of the power-law asymptotic dependence (with p=3/2). For LD model, we have chosen q=1.To compare the light curves for different source models1 and 2, we used the relative value:
(62)ηm=2〈maxt{η(t)}〉,η(t)=|K1(t)-K2(t)|K1(t)+K2(t),
where Ki(t) is the amplification for ith model along the trajectory of source linear motion; 〈⋯〉 is an average over a number of realizations. Because we expect that maximum difference of brightness on lightcurves for different source models takes place during HAE, so ηm typically is a characteristic of HAE.From the “light curves” on Figure2 we observe a significant difference between the “compact” (LD and Gaussian) and the “long-range” models. The long-range character of the latter reveals itself even on considerable distances from the caustics, where we expect that brightness of all the sources must have the same behavior as that of a point source. The differences between these two groups of models are essentially larger than the differences within each group (e.g., between Gaussian and LD models). This conclusion is confirmed by results of the statistical considerations over 100 realizations shown in Table 2 for half-brightness radius R1/2=0.21 as the example.Table 2
Relative difference between the light curves in HAE.
ith model
jth model
η
m
AD
Gaussian
0
.
074
±
0
.
001
AD
AD1
0
.
085
±
0
.
002
AD
LD
0
.
091
±
0
.
002
AD
PL
0
.
038
±
0
.
001
Gaussian
AD1
0
.
073
±
0
.
002
Gaussian
LD
0
.
042
±
0
.
001
Gaussian
PL
0
.
073
±
0
.
001
AD1
LD
0
.
094
±
0
.
002
AD1
PL
0
.
052
±
0
.
001
LD
PL
0
.
090
±
0
.
002Figure 2
“Light curves”: amplification as a function of time for different source models that correspond to magnification pattern of Figure1.We note that the results of comparison may depend on a complexity of the caustics involved into our consideration. One can expect, for example, that in case of complicated caustic crossing (such as crossing of the fold caustics in a vicinity of cusp, intersection of dense aggregations of caustics, among others) the light curve difference for different source models may be more significant. Our previous simulations involvedall possible realizations of the microlensing field that may invoke the complicated caustic crossings. However, in case of a concrete GLS, we deal with a single light curve. One may ask: is it relevant to apply the results of statistical simulations to this single light curve, where we may have a kind of an observational selection. Therefore, it is necessary to check how this “complexity” affects the average value of η. In this connection, we considered some modification of our statistical procedure: we have chosen (by eye) such realizations rather simple fold caustic crossings involved. However, the numerical results for simple caustic crossing events, such as shown as the example on Figure 3, appeared to be nearly the same as that of the Table 2. Though we must note that here the statistics have been worse because of a smaller number of the “simple” caustic realizations. So this conclusion requires further verification.Figure 3
Magnification pattern with simple caustic crossing events.


## 3.5.2. Gaussian Fittings of Accretion Disk and Limb Darkening Models

Table2 concerns a comparison of different models with the sameR1/2. However, in reality, we do not know what source we must fit and one must check whether we can replace one model with a different one with some other source parameters to get better fitting.We have fitted the limb-darkening and accretion disk model light curves with that of the Gaussian source of different radii. The half-brightness radius of LD and AD models was fixed:R1/2=0.2RE (this choice was determined by the AD1 model, see comments after formula (47)) and that of the Gaussian source was varied to get better fitting; here P=3/2 for PL model and q=1 for LD model. One example of the relative light curve differences η(T) according to (62) between LD, AD1 and the best fitted Gaussian models is shown in Figure 4. As we see from this figure, usage of an optimal size of the Gaussian model can diminish ηm as compared to the average estimate from Table 2 dealing with the same R1/2 for all the models.Figure 4
Relative differencesη(t) between Gaussian and LD (q=1), AD1 models for the best fitted curves corresponding to magnification pattern of Figure 3.To test this statement statistically, we carried out simulations with 20 magnification patterns having the “simple” caustic crossing. We found that fitting of AD light curves by the Gaussian model (with differentR1/2) yields ηm=0.07±0.006; this is roughly the same as in Table 2. Fittings of LD model light curves by the Gaussian model (with different R1/2) yielding ηm=0.026±0.002 (less then half of Table 2); in this case, an improvement is noticeable. Therefore, the Gaussian source cannot reproduce all the models, though for some compact models the fitting results are rather good. Then, in the caustic crossing events different models may be distinguished on the accuracy level of modern photometric observations.


## 3.6. Resume of Section3

Here the formulas of Section2 are applied to the Gaussian, power-law and limb darkening models of an extended source. The analytical formulas for these source models are obtained. We note that a fitting of the light curve of GLS Q2237+0305C [5, 6] showed that some of these higher-order corrections can be statistically significant even at the present accuracy level. This means that when we are looking for some fine effects in HAE due to the source size, a consistent treatment must also take into account possible effects of the above corrections.In case of the accretion disk models, we have not succeeded to find compact “workable” expressions forKcr, so this case has been treated numerically to compare possible observational differences of various source models.Now we remind the main assumptions and the class of models used, which is necessary for understanding the results of the Section3.5. The simulations use the GLS parameters similar to that of Q2237+030 [47, 51]. We consider the equal mass microlenses; we do not consider any mass distributions and/or populations of small (planetary) masses. Next, we consider most simple circular symmetric source models without effects of ellipticity, and so forth. The results of Section 3.5 can be summed up as follows. There are special classes of models: (i) accretion disk models having the dark region in the center and (ii) source models with a slow decrease of the (long-range) brightness distribution at large distances from the center. From our results, we see that (i) and (ii) show in HAEs a behavior which is different from that of the Gaussian and the limb-darkening models on a relative accuracy level of~4÷7%. There is then a possibility to distinguish such models on account of observations during HAEs. Outside HAE the surface brightness profile typically has little effect on microlensing (see [35, 52]).We note that our conclusions involve concrete source models with a special choice of parameters (source sizes, microlensing optical depth, etc.) and a considerable work is to be fulfilled to obtain more general results. Also, our results are of a statistical nature, so our average estimate ofηm may be different from analogous one in a separate light curve observation.Note that adding a sufficient number of microlenses having a small (e.g., planetary) masses (cf. [51, 53–55]) must produce an additional caustics that in its turn can give rise to fine features of light curves during HAE. One may expect that this can mimic features due to brightness distribution over the source. This question deserves further investigation.


## 4. Astrometric Gravitational Microlensing

### 4.1. Preliminary Comments

Besides photometric observations of microlensed sources, additional information about the source and lens structure may come from astrometric effects. Gravitational fields of masses cause the bending of light rays; this effect has passed through numerous experimental tests within the solar system. This bending leads to brightness amplification of distant sources in GLSs. We say about the astrometric gravitational microlensing when we deal with shifts and motions of the distant source images due to the gravitational fields of foreground stellar-mass objects (stars, black holes, dark matter substructures, etc., either in the Milky Way or in the other galaxies).Possibilities to detect an image motion of the remote source caused by the gravitational field of the Milky Way stars have been discussed for a long time [56–60]. The problem gained a new interest in connection with investigations of massive compact halo objects [61–63], see also [64–71] and references therein; for the extended source and extended lens effects see [11, 72–74]. The extragalactic topics are less studied [12, 71, 75, 76]. Nevertheless, the microarcsecond accuracy level appears to be sufficient to resolve the shape of a source image trajectory in extragalactic GLS in microlensing events. The corresponding image shifts can achieve the level of~10μas.Astrometric measurements can provide a valuable information about mass distributions of microlenses, which can complement the photometric observations of GLS images [77, 78]. Observations of trajectories of the remote source images during the microlensing events, in addition to the photometric data, would give us a possibility to estimate the foreground field parameters and the optical depth of the continuous and discrete matter. This may be helpful in connection with the dark matter (DM) problem; in particular, this can be used to test the existence of the putative DM subhalos [79]. Though such positional effects are still beyond the accuracy of observations, the astrometric microlensing attracts a considerable attention during the last decade in connection with perspectives of spaceborn microarcsecond astrometry (the forthcoming Gaia mission, space-based VLBI missions, optical Space Interferometry Mission). It should be pointed out that, typically, the accuracy of astrometric positioning is essentially higher than the image resolution. The other positive feature of astrometric microlensing events is the slow dependence on the impact parameter (~p-1) and, hence, there is a larger probability of observing such events, once a required accuracy will be achieved.There is a number of candidates for astrometric microlensing in our galaxy that can be detected in the near future [80–82]. On the other hand, Klioner [83] points out that unpredictable microlensing noise can spoil the determination of positions and proper motions of the objects resulting from future astrometric missions on a submicroarcsecond level. The gravitational image motion (GIM) due to microlensing effects may be comparable to proper motion of quasars [8, 9, 84–86].


### 4.2. Extended Source and Extended Lens

Further, dealing with trajectories of the microlensed source images, we mean the trajectory of the source image centroid (SIC), that is, a weighted average over all elements of all the source images. The position of SIC can be defined as(63)R(Y)=1F(Y)∬xI(y(x)-Y)dx1dx2,x=(x1,x2), y(x) is the lens mapping, Y=(Y1,Y2) is the center of the source in the source plane; F(Y) is given by formula (35).Having in mind that most of the astrometric microlensing effects are still outside observational possibilities, for astrometric purposes, it is quite reasonable to confine ourselves to the most simple source brightness profiles described in Section3.3 We only note that for the microlensing of a star the limb darkening model is most appropriate, and structures of real extragalactic sources (e.g., quasars) are more complicated.In this subsection, we consider microlensing by one mass in presence of an external shear. This can be the case of microlensing of a remote source by a small mass (planet) in presence of a background gravitational field of a star or microlensing by a star in the field of a putative DM clump.First, we consider a weak microlensing system, that is, the linear distance from a microlensing mass to the line of sight is much larger than its Einstein radius. In this case, if a remote point source is microlensed by a foreground moving point mass, the trajectory of the source image is approximately the circle in the rest frame of the source [9, 66]. In case of strong lensing HAE (linear distance from the microlens to the line of sight is of the order of the Einstein radius), this circle will be distorted. In any case, the point source image is “repulsed” from the lensing point mass.The trajectory qualitatively changes when the projection of the point mass onto the source plane crosses the extended source. In this case, the SIC motion trajectory typically has a three-leaf form [11, 72]. This has been first shown in case of a source with uniform surface brightness [72] and in case of the Gaussian source [11]. The reason for such behavior is that when the lens is projected onto the internal part of the source, the “repulsion” of SIC is changed by its “attraction” to the lens due to larger amplification near it. For sufficiently small values of impact parameter P of the lens with respect to the extended source center, the SIC trajectory bends and passes through the source center forming a self-intersection.The external field effects, which are essential in extragalactic GLS, distort considerably the trajectory of SIC [12, 13, 87]. Below, we present some simple examples of microlensing of an extended Gaussian source in presence of an external shear (Chang-Refsdal lens [88]). The configuration of the figures is as follows: the left hand panels show the SIC trajectories; the right-hand ones show corresponding dependencies of the amplification K(t) as functions of time. The SIC trajectories are shown in the rest frame of the source; in the right upper corner the source center trajectory is shown (schematically) with respect to the caustic. We show in figures the impact parameter (P) and the shear (γ).Figure5 shows the differences between the point-like source and the extended Gaussian one (42) in presence of an external shear. To compare astrometrical and photometrical effects for different source models, we represent in (Figures 6 and 7) the SIC trajectories and light curves of the Gaussian source and accretion disk (47) in band V, both with the same half-brightness radius of 0.08RE. The last one, Figure 8 shows the SIC trajectories in case of finite-size microlensing masses. We note that the SIC trajectory depends on the source trajectory and upon the source model more distinctly than the light curve.Figure 5
Microlensing in the Chang-Refsdal lens for two different sources: point-like one and Gaussian one (42) with r.m.s. radius R=0.1RE (see (40); this corresponds to the half-brightness radius (41) R1/2≈0.08RE). The source center crosses two caustic folds.Figure 6
Microlensing in the Chang-Refsdal lens: the source center passes over the cusp of the caustic. Source models: Gauss (42) and AD1 (47).Figure 7
Microlensing in the Chang-Refsdal lens with large convergence. Source models: Gauss (42) and AD1 (47).Figure 8
Microlensing of the Gaussian source by the Zakharov and Sazhin lens [89] of different lens sizes RL. The source size is 0.1RE.The simplest situation is when the source track passes far from the caustic. In this case, the SIC trajectory is oval-like and light curve is continuous. However, when the source center track is close to the caustic, even if it does not cross it, the trajectory becomes much more complicated. In general, it can be self-intersected several times. When point source crosses the caustic, its SIC suffers a jump (Figure5). This is explained by the emergence or by disappearance of pairs of images of infinite brightness, leading to sudden jumps of the weighted average of positions of the source elements. Analytical expressions for such jumps can be found in [13]. For an extended source, corresponding discontinuities of the trajectory are smoothed out but may be quite noticeable. When the small source passes near the cusp or crosses it, the total images brightness increases dramatically, but no jump of image motion occurs because of the existence of one more image of infinite (in case of a point source) brightness which does not disappear. This distinguishes the high amplification events corresponding to the fold and cusp caustic crossings.The effect of high optical depth of continuous matter is demonstrated in Figure7 to show possible observational signals of continuous DM near the line of the sight.In papers [73, 89], astrometric and photometric effects of GLS with an extended microlens were considered. The microlensing effect caused by a compact body (e.g., a star) must be distinguished from the eclipse-like effect in case of a presence of quite dense but essentially continuous background (cf. “subhalos”) of DM. Such “subhalos” are supposed to be extended clumps of DM. The size of these clumps predicted by ΛCDM model is larger than typical Einstein radius (see, e.g., Diemand et al. [90]). Such an object can manifest itself in microlensing process acting together with another compact object. In this case, smooth-shaped or even quasielliptical caustics can exist instead of, or together with, cusped ones (Figure 7). When the source crosses the fold of the smooth caustic, two images of it disappear. As a result, after some peak of the image brightness, there may be swift diminution. Thus DM can cause the effect like “eclipse,” demonstrated on Figure 7. In an opposite case, when a subhalo have finite size smaller than the Einstein radius, finite-size lens models can be applied to describe the situation. As it was shown in [89], the finite-size opaque lens can cause “ears” on the trajectory when one of the images is obscured by the lens, and the convergence causes the trajectory disruption when the source crosses caustic. In the light curve, a finite lens size manifests itself in the magnification jumps when the image enters in the lens shadow. In case of an overcritical convergence, the magnification of the point source images becomes infinite when the source crosses the caustic, and the brightness decays swiftly when the source is inside the caustic.We show here also some results obtained within the other “toy-example” for the microlensing by a subhalo of the sizea modelled by lens equation (60) with the only extended lens (N=1). The results obtained within this model are shown on the Figure 9. As we can see here, the finite size of the subhalo causes an apparition of the additional quasielliptical caustic. Meanwhile, the SIC trajectories do not differ principally from the point-lens situation. However, here, we can see the dip on the lightcurve similar to that obtained in Zakharov and Sazhin lens model [89].Figure 9
Microlensing by the “subhalo” lens: Gauss (42) and AD1 (47).The effect of the finite-size lens manifest itself in smoother way if the source is extended (see Figure8). In this case, it changes the SIC trajectory and the light curve has no jumps or gaps. The most prominent effect on the light curve is in the recession of the dip between two caustic crossings (i.e., when source is inside the caustic).


### 4.3. Statistics of Weak Astrometric Microlensing

#### 4.3.1.   Correlation Tensor of Image Shifts

In a general case of a complicated microlensing system, we need some statistical measure of image motion; possible characteristic is proposed below.Letδx describes a deviation of the point source image from unperturbed position. In case of one microlens, this value typically varies as r-1, where r is a distance of the microlens from the line of sight. Because of this, we expect that in case of uniform distribution of many microlenses usual correlation function and 〈δx2〉 diverge, and so we prefer to deal with a correlation tensor like that used in the turbulence theory [91]:
(64)Dij(τ)≡〈(δxi(t)-δxi(t-τ))(δxj(t)-δxj(t-τ))〉.
This tensor is symmetric and it has a diagonal form in a system where the source moves along one of coordinate axes, for example, along y1-axis. So a determination of Dij(τ) from observations may help to constrain a direction of the source motion.We will illustrate calculation ofDij in case of weak astrometric microlensing by N clumps having a uniform distribution over the lens plane. The system is described by the lens equation (60). We suppose that all the clumps are independent, and xn has a uniform distribution over the lens plane. To be precise, we assume that the spatial distribution function for every xn is
(65){S-1,xn∈Ω,0,xn∉Ω,
where S is the area of domain Ω. Further considerations deal with the limit of Ω that will be blown up to cover the whole plane assuming
(66)S(Ω)N=s(a,M)=const,N→∞,
where we take into account that the surface per one clump s(a,M) may be different for clumps of different mass and size. Density f(a,M,x)~|x|-4 goes to zero sufficiently fast as x→∞, so the clumps have a finite size and integral M=∫‍d2xf(a,M,x)<∞ is convergent. To save the space, we do not consider a distribution of parameters a,M assuming that these parameters are the same. This means that in the general case the results below must be averaged over possible a,M.According to (60), in case of weak microlensing the shift of the point source image is
(67)δx=∑n=1NRE,n2y-xn|y-xn|2+an2.For a straight line motion of a source, we puty1=vt,y2=const. Then, we have a diagonal form of Dij with
(68)D11(τ)=RE4s∫d2x[x1,-x1,-2+x22+a2-x1,+x1,+2+x22+a2]2,D22(τ)=RE4s∫d2x[x2x1,-2+x22+a2-x2x1,+2+x22+a2]2,x1,±=x1±yτ,yτ≡|v|τ2.Calculation of the integrals yields(69)D11(τ)=πRE4sF11(yτa),D22(τ)=πRE4sF22(yτa),
where
(70)F11(ξ)=ηln1+η1-η,F22(ξ)=1ηln1+η1-η-2,η≡ξξ2+1.


#### 4.3.2. Statistical Effect of Gravitational Dragging of a Distant Source Image by Point Masses

In this section, we consider statistical effects in the motion of the GLS images due to gravitational field of a large number of moving stars. Here, we consider a rarefied collection of lensing masses, when the typical distances between the source and the microlens are much larger than the corresponding Einstein radius. For brevity, we call this “weak” microlensing; this must be distinguished from the well-known extragalactic weak lensing, that is, image distortions caused by the cosmic shear. Unlike the previous sections, here we rule out strong microlensing events that are accompanied by a considerable brightness amplification.In some situations it will be difficult to separate GIM from real proper motions, and this is important for accuracy of the fundamental reference frame based on extragalactic sources [83, 85, 86]. The gravitational field of point masses (which is essentially inhomogeneous) induces a stochastic motion of the source image in the reference frame of the source. On the other hand, smoothly distributed masses also induce a stochastic motion, but typical amplitude of the effect will be weaker. Indeed, we cannot take into account positions of all stars [83]; we only can work with their probability distribution leading to distribution of image motions of a source [8, 9]. The situation could be further complicated in presence of invisible DM clumps, if they exist on small scales.In case of the point masses, it is possible to estimate the probability distribution of image motionsu (i.e., angular velocities) of the distant source image microlensed by collection of a randomly moving point masses [8]:
(71)P(u)=u02π(u02+u2)3/2,
where u0 is a characteristic value depending upon the transverse velocity of point masses, their value, and their spatial density. It should be noted that the distribution (71) is obtained under the approximation of the weak microlensing [8], so that HAEs must be excluded.The distribution (71) has a slowly decreasing tail yielding an infinite estimate of velocity squared 〈u2〉. This tail is a formal consequence of pointness of the microlenses that makes possible very high GIMs. In fact this case is not covered by formula (71), which is valid for weak microlensing only. In case of smooth DM concentrations such velocities would be absent.If a bulk motion of microlensing masses is present, the average〈u〉 requires a special consideration. It has been pointed out in [8, 9] that stochastic GIM of an extragalactic source induced by stellar motions is accompanied with a systematic component (not a statistical average) which depends upon bulk velocity of microlensing stars, including stars that are far away from the line of sight. The typical velocity of this component remains nonzero even for a homogeneous distribution of microlensing masses which move in certain direction (all masses have the same collective velocity component). This needs some explanation because the average position of the unmoving source image (averaged over a very long time) must be unchanged and, evidently, the statistical average of image velocity must be zero. However, when dealing with the extragalactic reference frame, in reality, we do not deal with a very large observational time and/or with a sufficiently large number of sources to obtain a statistical average. In this case, because the strong microlensing events yielding the fast image motions are very rare, they in practice are not observable unless the number of simultaneously observed sources is of the order of millions or more. This is well known from theory and observations of MACHOs. Then the statistical average over set of all events is not a good characteristic of GIM; the most probable value of GIM (which is different from zero) would be a better measure.Otherwise, we can estimate probability for all necessary velocity intervals. Detailed approach on this way, however, involves consideration of strong lensing events with different extended source models. However, to have an order-of-magnitude estimate, we can avoid calculations of probability distributions, if we confine ourselves to the domain of weak microlensing events (“W”, for brevity) [9, 10]. Typical velocity which is a characteristic of such events is a statistical average over the domain W, with HAE being excluded. Note that these rare strong lensing events (S) essentially contribute to the usual statistical average velocity 〈u〉W+S (i.e., average over all the space of events S+W). If we exclude S-events (this is easy to do in real observations because of strong brightness amplification of a source), then we deal with “truncated” average velocity 〈u〉W over W events. This value is nonzero even in case of a homogeneous distribution of microlenses having a collective velocity. It is important to note that this value practically does not depend upon the exact definition of S-domain (i.e., upon a value of minimal impact parameter of microlensing masses with respect to the line of sight, which determines events “S”).Note that the “truncated” average would better correspond to realistic observations when the typical velocity is estimated as an arithmetic mean over all reference sources used in extragalactic International Celestial Reference Frame (ICRF). In a near future, the number of these sources will not exceed several hundreds. Therefore, in fact, we deal with most probable dispositions of lensing masses which typically do not include theS-events: it is well known that the probability of the such events is of the order of 10-6 per year in the Milky Way galaxy. Then, in ICRF considerations, it is reasonable to exclude such rare events that would hardly occur during the century.In case of microlensing of distant sources by Milky Way stars, we may neglect cosmological curvature. Evidently, this does not contradict to assumption that the radiation source is at the infinity. Let unperturbed light ray move from the infinitely distant radiation source in negative direction ofz-axis of Cartesian coordinates {x,y,z}, the observer being at the origin. Let position of microlensing point mass M be (r,z), where r=(x,y) is a two-dimensional vector in the transverse plane; that is, r=|r| is the impact distance of the unperturbed ray with respect to the mass. In the case of weak microlensing, r≫(mz)1/2, m=GM/c2. In this case,
(72)Ψ=mΨ~,Ψ~=-2rr2[1+z(z2+r2)-1/2],
where Ψ=(Ψ1,Ψ2) is a two-dimensional vector describing the source image angular shift [9, 85, 92].FurtherVp=(v,w) stands for velocity of the point microlens, w is the velocity component parallel to the line of sight, and v represents the transversal components. In virtue of (72), this leads to the source GIM that equals to mU (in radians per unit of time), where
(73)U=dΨ~dt=-2{wr-zv/2(z2+r2)3/21r4[vr2-2r(r·v)]U=dΨ~dt=-2·[1+zz2+r2+zr22(z2+r2)3/2]U=dΨ~dt=-2+wr-zv/2(z2+r2)3/2}.
Here, U is a function of the microlens position (r,z) and its velocity v.Because of the smallness of the effect, the action of all Milky Way stars will be taken into account linearly. This will be performed by integrating (73). However, this requires justification. In fact, we are dealing with extended sources; this is not taken into account by  (73). As we will see below, this may be important in our problem.Thus, we return to (73), which describes contribution of a single star. We assume that a star at the point (r,z) has velocity Vp(r,z). Then, we consider sum over all stars in the Milky Way, which must further be averaged with the mass density ρ(r,z). This enables us to pass to integration to yield an average GIM in domain W:
(74)〈Utot〉W=Gc2∫dz∫d2rρ(r,z)U(r,z,Vp(r,z)),
where we suppose that ρ vanishes outside a bounded domain. Also, we suppose that ρ(r) is sufficiently smooth so as we can use approximation of weak microlensing inside the domain where ρ(r)≠0.Taking into account the explicit form (73), it is easy to see that in a general case considerable contribution may be due to stars at large impact distances from the line of sight. Also, one can show that the singularities in (74) for small r are integrable and the integral is convergent. This allows us to avoid the question about exact value of the lower limit of r in the definition of domain W.It is interesting to compare GIM (74) with an average over all S+W events in case of a general continuous mass distribution ρ=ρ(r,z,t). This distribution must satisfy the continuity equation:
(75)∂ρ∂t+div(ρVp)=0.In this case, (72) should be integrated with the density ρ:
(76)Ψtot=Gc2∫dz∫d2rρ(r,z,t)Ψ~(r,z).Note that this formula forΨtot can be considered either as the average over all S+W events or simply as a shift in case of mass distribution ρ(r).The shift of remote source image in case of continuous mass distribution has been treated in [93]. These shifts appear to be of the order of 10-5÷10-6 radians, but this value cannot be observed from the solar system. The observable (at least in principle) value is a change of (76) with time:
(77)Utot*≡dΨtotdt=Gc2∫dz∫d2r∂ρ∂tΨ~(r,z).
Formula (77) is, in fact, an average image velocity over all S+W events:
(78)Utot*=〈Utot〉S+W.
One might think that this must be the same as (74), but it is not true. The reason is that operations of averaging and differentiation here do not commute. To compare (74) with (77), one must use the continuity equation (75) and avoid singular points at r=0 in Utot* [10].We haveUtot*=U(0)+U(1), where U(0) and U(1) are two-dimensional vectors, and the components of the first term (i=1,2):
(79)Ui(0)=Gc2∫dz∫d2rρ(Vp·∇Ψ~i),
are the same as the components of (74): U(0)=〈Utot〉W. The components of the second term can be written as
(80)Ui(1)=-Gc2∫0∞dz∫r>ε→0d2rdiv(VpρΨ~i).
After some calculation on account of the Ostrogradsky-Gauss theorem, we have [10]
(81)〈Utot〉S+W-〈Utot〉W≡U(1)=-4πGc2∫0∞dzρ(z)V⊥(z),
where ρ(z) stands for the density of microlensing masses along the line of sight, V⊥ is the transverse velocity of microlenses near the line of sight. Evidently, (81), differs from zero. For example, if an observer at rest is situated in the center of Milky Way (stationary mass distribution), we have Utot*=0 (see (77)), but GIM does not equal to zero.For an observer in the solar system this effect yields a nonzero GIM value leading to apparent rotation of extragalactic reference frame in the direction of the Milky Way rotation [10]. For sources near the Galactic plane, the effect amounts about ~10-8 arcseconds per year [10]. The effect will be even smaller for sources with higher Galactic latitudes. Therefore, in principle (if the corresponding accuracy be achieved!), the effect can be observed from the Solar system.


### 4.4. Resume of Section4

In this Section, we review main points of astrometric microlensing. We demonstrate effects due to a finite size of the source and extended lens as well; a special attention is paid to effects that may serve as signals of dark matter. Several examples of microlensing by one mass in presence of an external shear are given. This can be the case of microlensing of a remote source by a small mass (planet) in presence of a background gravitational field of a star or microlensing by a star in the field of a putative DM clump.Then, we proceed to some statistical subjects of astrometric microlensing. We introduced the correlation tensor as a statistical characteristic of a stochastic image motion and applied this to a “toy” model of weak microlensing by dark matter clumps.We have shown that gravitational field of foreground stars induces an image motion, which differs for a continuous and discrete mass density distribution. In the latter case, image performs random walks, as distinct from a regular motion in case of continuous matter, plus some additional motion due to the collective motion of the microlensing masses. This effect leads to apparent rotation of extragalactic reference frame in the direction of the Milky Way rotation. In principle, the effect can be observed, however, it is far beyond modern possibilities.


## 4.1. Preliminary Comments

Besides photometric observations of microlensed sources, additional information about the source and lens structure may come from astrometric effects. Gravitational fields of masses cause the bending of light rays; this effect has passed through numerous experimental tests within the solar system. This bending leads to brightness amplification of distant sources in GLSs. We say about the astrometric gravitational microlensing when we deal with shifts and motions of the distant source images due to the gravitational fields of foreground stellar-mass objects (stars, black holes, dark matter substructures, etc., either in the Milky Way or in the other galaxies).Possibilities to detect an image motion of the remote source caused by the gravitational field of the Milky Way stars have been discussed for a long time [56–60]. The problem gained a new interest in connection with investigations of massive compact halo objects [61–63], see also [64–71] and references therein; for the extended source and extended lens effects see [11, 72–74]. The extragalactic topics are less studied [12, 71, 75, 76]. Nevertheless, the microarcsecond accuracy level appears to be sufficient to resolve the shape of a source image trajectory in extragalactic GLS in microlensing events. The corresponding image shifts can achieve the level of~10μas.Astrometric measurements can provide a valuable information about mass distributions of microlenses, which can complement the photometric observations of GLS images [77, 78]. Observations of trajectories of the remote source images during the microlensing events, in addition to the photometric data, would give us a possibility to estimate the foreground field parameters and the optical depth of the continuous and discrete matter. This may be helpful in connection with the dark matter (DM) problem; in particular, this can be used to test the existence of the putative DM subhalos [79]. Though such positional effects are still beyond the accuracy of observations, the astrometric microlensing attracts a considerable attention during the last decade in connection with perspectives of spaceborn microarcsecond astrometry (the forthcoming Gaia mission, space-based VLBI missions, optical Space Interferometry Mission). It should be pointed out that, typically, the accuracy of astrometric positioning is essentially higher than the image resolution. The other positive feature of astrometric microlensing events is the slow dependence on the impact parameter (~p-1) and, hence, there is a larger probability of observing such events, once a required accuracy will be achieved.There is a number of candidates for astrometric microlensing in our galaxy that can be detected in the near future [80–82]. On the other hand, Klioner [83] points out that unpredictable microlensing noise can spoil the determination of positions and proper motions of the objects resulting from future astrometric missions on a submicroarcsecond level. The gravitational image motion (GIM) due to microlensing effects may be comparable to proper motion of quasars [8, 9, 84–86].


## 4.2. Extended Source and Extended Lens

Further, dealing with trajectories of the microlensed source images, we mean the trajectory of the source image centroid (SIC), that is, a weighted average over all elements of all the source images. The position of SIC can be defined as(63)R(Y)=1F(Y)∬xI(y(x)-Y)dx1dx2,x=(x1,x2), y(x) is the lens mapping, Y=(Y1,Y2) is the center of the source in the source plane; F(Y) is given by formula (35).Having in mind that most of the astrometric microlensing effects are still outside observational possibilities, for astrometric purposes, it is quite reasonable to confine ourselves to the most simple source brightness profiles described in Section3.3 We only note that for the microlensing of a star the limb darkening model is most appropriate, and structures of real extragalactic sources (e.g., quasars) are more complicated.In this subsection, we consider microlensing by one mass in presence of an external shear. This can be the case of microlensing of a remote source by a small mass (planet) in presence of a background gravitational field of a star or microlensing by a star in the field of a putative DM clump.First, we consider a weak microlensing system, that is, the linear distance from a microlensing mass to the line of sight is much larger than its Einstein radius. In this case, if a remote point source is microlensed by a foreground moving point mass, the trajectory of the source image is approximately the circle in the rest frame of the source [9, 66]. In case of strong lensing HAE (linear distance from the microlens to the line of sight is of the order of the Einstein radius), this circle will be distorted. In any case, the point source image is “repulsed” from the lensing point mass.The trajectory qualitatively changes when the projection of the point mass onto the source plane crosses the extended source. In this case, the SIC motion trajectory typically has a three-leaf form [11, 72]. This has been first shown in case of a source with uniform surface brightness [72] and in case of the Gaussian source [11]. The reason for such behavior is that when the lens is projected onto the internal part of the source, the “repulsion” of SIC is changed by its “attraction” to the lens due to larger amplification near it. For sufficiently small values of impact parameter P of the lens with respect to the extended source center, the SIC trajectory bends and passes through the source center forming a self-intersection.The external field effects, which are essential in extragalactic GLS, distort considerably the trajectory of SIC [12, 13, 87]. Below, we present some simple examples of microlensing of an extended Gaussian source in presence of an external shear (Chang-Refsdal lens [88]). The configuration of the figures is as follows: the left hand panels show the SIC trajectories; the right-hand ones show corresponding dependencies of the amplification K(t) as functions of time. The SIC trajectories are shown in the rest frame of the source; in the right upper corner the source center trajectory is shown (schematically) with respect to the caustic. We show in figures the impact parameter (P) and the shear (γ).Figure5 shows the differences between the point-like source and the extended Gaussian one (42) in presence of an external shear. To compare astrometrical and photometrical effects for different source models, we represent in (Figures 6 and 7) the SIC trajectories and light curves of the Gaussian source and accretion disk (47) in band V, both with the same half-brightness radius of 0.08RE. The last one, Figure 8 shows the SIC trajectories in case of finite-size microlensing masses. We note that the SIC trajectory depends on the source trajectory and upon the source model more distinctly than the light curve.Figure 5
Microlensing in the Chang-Refsdal lens for two different sources: point-like one and Gaussian one (42) with r.m.s. radius R=0.1RE (see (40); this corresponds to the half-brightness radius (41) R1/2≈0.08RE). The source center crosses two caustic folds.Figure 6
Microlensing in the Chang-Refsdal lens: the source center passes over the cusp of the caustic. Source models: Gauss (42) and AD1 (47).Figure 7
Microlensing in the Chang-Refsdal lens with large convergence. Source models: Gauss (42) and AD1 (47).Figure 8
Microlensing of the Gaussian source by the Zakharov and Sazhin lens [89] of different lens sizes RL. The source size is 0.1RE.The simplest situation is when the source track passes far from the caustic. In this case, the SIC trajectory is oval-like and light curve is continuous. However, when the source center track is close to the caustic, even if it does not cross it, the trajectory becomes much more complicated. In general, it can be self-intersected several times. When point source crosses the caustic, its SIC suffers a jump (Figure5). This is explained by the emergence or by disappearance of pairs of images of infinite brightness, leading to sudden jumps of the weighted average of positions of the source elements. Analytical expressions for such jumps can be found in [13]. For an extended source, corresponding discontinuities of the trajectory are smoothed out but may be quite noticeable. When the small source passes near the cusp or crosses it, the total images brightness increases dramatically, but no jump of image motion occurs because of the existence of one more image of infinite (in case of a point source) brightness which does not disappear. This distinguishes the high amplification events corresponding to the fold and cusp caustic crossings.The effect of high optical depth of continuous matter is demonstrated in Figure7 to show possible observational signals of continuous DM near the line of the sight.In papers [73, 89], astrometric and photometric effects of GLS with an extended microlens were considered. The microlensing effect caused by a compact body (e.g., a star) must be distinguished from the eclipse-like effect in case of a presence of quite dense but essentially continuous background (cf. “subhalos”) of DM. Such “subhalos” are supposed to be extended clumps of DM. The size of these clumps predicted by ΛCDM model is larger than typical Einstein radius (see, e.g., Diemand et al. [90]). Such an object can manifest itself in microlensing process acting together with another compact object. In this case, smooth-shaped or even quasielliptical caustics can exist instead of, or together with, cusped ones (Figure 7). When the source crosses the fold of the smooth caustic, two images of it disappear. As a result, after some peak of the image brightness, there may be swift diminution. Thus DM can cause the effect like “eclipse,” demonstrated on Figure 7. In an opposite case, when a subhalo have finite size smaller than the Einstein radius, finite-size lens models can be applied to describe the situation. As it was shown in [89], the finite-size opaque lens can cause “ears” on the trajectory when one of the images is obscured by the lens, and the convergence causes the trajectory disruption when the source crosses caustic. In the light curve, a finite lens size manifests itself in the magnification jumps when the image enters in the lens shadow. In case of an overcritical convergence, the magnification of the point source images becomes infinite when the source crosses the caustic, and the brightness decays swiftly when the source is inside the caustic.We show here also some results obtained within the other “toy-example” for the microlensing by a subhalo of the sizea modelled by lens equation (60) with the only extended lens (N=1). The results obtained within this model are shown on the Figure 9. As we can see here, the finite size of the subhalo causes an apparition of the additional quasielliptical caustic. Meanwhile, the SIC trajectories do not differ principally from the point-lens situation. However, here, we can see the dip on the lightcurve similar to that obtained in Zakharov and Sazhin lens model [89].Figure 9
Microlensing by the “subhalo” lens: Gauss (42) and AD1 (47).The effect of the finite-size lens manifest itself in smoother way if the source is extended (see Figure8). In this case, it changes the SIC trajectory and the light curve has no jumps or gaps. The most prominent effect on the light curve is in the recession of the dip between two caustic crossings (i.e., when source is inside the caustic).


## 4.3. Statistics of Weak Astrometric Microlensing

### 4.3.1.   Correlation Tensor of Image Shifts

In a general case of a complicated microlensing system, we need some statistical measure of image motion; possible characteristic is proposed below.Letδx describes a deviation of the point source image from unperturbed position. In case of one microlens, this value typically varies as r-1, where r is a distance of the microlens from the line of sight. Because of this, we expect that in case of uniform distribution of many microlenses usual correlation function and 〈δx2〉 diverge, and so we prefer to deal with a correlation tensor like that used in the turbulence theory [91]:
(64)Dij(τ)≡〈(δxi(t)-δxi(t-τ))(δxj(t)-δxj(t-τ))〉.
This tensor is symmetric and it has a diagonal form in a system where the source moves along one of coordinate axes, for example, along y1-axis. So a determination of Dij(τ) from observations may help to constrain a direction of the source motion.We will illustrate calculation ofDij in case of weak astrometric microlensing by N clumps having a uniform distribution over the lens plane. The system is described by the lens equation (60). We suppose that all the clumps are independent, and xn has a uniform distribution over the lens plane. To be precise, we assume that the spatial distribution function for every xn is
(65){S-1,xn∈Ω,0,xn∉Ω,
where S is the area of domain Ω. Further considerations deal with the limit of Ω that will be blown up to cover the whole plane assuming
(66)S(Ω)N=s(a,M)=const,N→∞,
where we take into account that the surface per one clump s(a,M) may be different for clumps of different mass and size. Density f(a,M,x)~|x|-4 goes to zero sufficiently fast as x→∞, so the clumps have a finite size and integral M=∫‍d2xf(a,M,x)<∞ is convergent. To save the space, we do not consider a distribution of parameters a,M assuming that these parameters are the same. This means that in the general case the results below must be averaged over possible a,M.According to (60), in case of weak microlensing the shift of the point source image is
(67)δx=∑n=1NRE,n2y-xn|y-xn|2+an2.For a straight line motion of a source, we puty1=vt,y2=const. Then, we have a diagonal form of Dij with
(68)D11(τ)=RE4s∫d2x[x1,-x1,-2+x22+a2-x1,+x1,+2+x22+a2]2,D22(τ)=RE4s∫d2x[x2x1,-2+x22+a2-x2x1,+2+x22+a2]2,x1,±=x1±yτ,yτ≡|v|τ2.Calculation of the integrals yields(69)D11(τ)=πRE4sF11(yτa),D22(τ)=πRE4sF22(yτa),
where
(70)F11(ξ)=ηln1+η1-η,F22(ξ)=1ηln1+η1-η-2,η≡ξξ2+1.


### 4.3.2. Statistical Effect of Gravitational Dragging of a Distant Source Image by Point Masses

In this section, we consider statistical effects in the motion of the GLS images due to gravitational field of a large number of moving stars. Here, we consider a rarefied collection of lensing masses, when the typical distances between the source and the microlens are much larger than the corresponding Einstein radius. For brevity, we call this “weak” microlensing; this must be distinguished from the well-known extragalactic weak lensing, that is, image distortions caused by the cosmic shear. Unlike the previous sections, here we rule out strong microlensing events that are accompanied by a considerable brightness amplification.In some situations it will be difficult to separate GIM from real proper motions, and this is important for accuracy of the fundamental reference frame based on extragalactic sources [83, 85, 86]. The gravitational field of point masses (which is essentially inhomogeneous) induces a stochastic motion of the source image in the reference frame of the source. On the other hand, smoothly distributed masses also induce a stochastic motion, but typical amplitude of the effect will be weaker. Indeed, we cannot take into account positions of all stars [83]; we only can work with their probability distribution leading to distribution of image motions of a source [8, 9]. The situation could be further complicated in presence of invisible DM clumps, if they exist on small scales.In case of the point masses, it is possible to estimate the probability distribution of image motionsu (i.e., angular velocities) of the distant source image microlensed by collection of a randomly moving point masses [8]:
(71)P(u)=u02π(u02+u2)3/2,
where u0 is a characteristic value depending upon the transverse velocity of point masses, their value, and their spatial density. It should be noted that the distribution (71) is obtained under the approximation of the weak microlensing [8], so that HAEs must be excluded.The distribution (71) has a slowly decreasing tail yielding an infinite estimate of velocity squared 〈u2〉. This tail is a formal consequence of pointness of the microlenses that makes possible very high GIMs. In fact this case is not covered by formula (71), which is valid for weak microlensing only. In case of smooth DM concentrations such velocities would be absent.If a bulk motion of microlensing masses is present, the average〈u〉 requires a special consideration. It has been pointed out in [8, 9] that stochastic GIM of an extragalactic source induced by stellar motions is accompanied with a systematic component (not a statistical average) which depends upon bulk velocity of microlensing stars, including stars that are far away from the line of sight. The typical velocity of this component remains nonzero even for a homogeneous distribution of microlensing masses which move in certain direction (all masses have the same collective velocity component). This needs some explanation because the average position of the unmoving source image (averaged over a very long time) must be unchanged and, evidently, the statistical average of image velocity must be zero. However, when dealing with the extragalactic reference frame, in reality, we do not deal with a very large observational time and/or with a sufficiently large number of sources to obtain a statistical average. In this case, because the strong microlensing events yielding the fast image motions are very rare, they in practice are not observable unless the number of simultaneously observed sources is of the order of millions or more. This is well known from theory and observations of MACHOs. Then the statistical average over set of all events is not a good characteristic of GIM; the most probable value of GIM (which is different from zero) would be a better measure.Otherwise, we can estimate probability for all necessary velocity intervals. Detailed approach on this way, however, involves consideration of strong lensing events with different extended source models. However, to have an order-of-magnitude estimate, we can avoid calculations of probability distributions, if we confine ourselves to the domain of weak microlensing events (“W”, for brevity) [9, 10]. Typical velocity which is a characteristic of such events is a statistical average over the domain W, with HAE being excluded. Note that these rare strong lensing events (S) essentially contribute to the usual statistical average velocity 〈u〉W+S (i.e., average over all the space of events S+W). If we exclude S-events (this is easy to do in real observations because of strong brightness amplification of a source), then we deal with “truncated” average velocity 〈u〉W over W events. This value is nonzero even in case of a homogeneous distribution of microlenses having a collective velocity. It is important to note that this value practically does not depend upon the exact definition of S-domain (i.e., upon a value of minimal impact parameter of microlensing masses with respect to the line of sight, which determines events “S”).Note that the “truncated” average would better correspond to realistic observations when the typical velocity is estimated as an arithmetic mean over all reference sources used in extragalactic International Celestial Reference Frame (ICRF). In a near future, the number of these sources will not exceed several hundreds. Therefore, in fact, we deal with most probable dispositions of lensing masses which typically do not include theS-events: it is well known that the probability of the such events is of the order of 10-6 per year in the Milky Way galaxy. Then, in ICRF considerations, it is reasonable to exclude such rare events that would hardly occur during the century.In case of microlensing of distant sources by Milky Way stars, we may neglect cosmological curvature. Evidently, this does not contradict to assumption that the radiation source is at the infinity. Let unperturbed light ray move from the infinitely distant radiation source in negative direction ofz-axis of Cartesian coordinates {x,y,z}, the observer being at the origin. Let position of microlensing point mass M be (r,z), where r=(x,y) is a two-dimensional vector in the transverse plane; that is, r=|r| is the impact distance of the unperturbed ray with respect to the mass. In the case of weak microlensing, r≫(mz)1/2, m=GM/c2. In this case,
(72)Ψ=mΨ~,Ψ~=-2rr2[1+z(z2+r2)-1/2],
where Ψ=(Ψ1,Ψ2) is a two-dimensional vector describing the source image angular shift [9, 85, 92].FurtherVp=(v,w) stands for velocity of the point microlens, w is the velocity component parallel to the line of sight, and v represents the transversal components. In virtue of (72), this leads to the source GIM that equals to mU (in radians per unit of time), where
(73)U=dΨ~dt=-2{wr-zv/2(z2+r2)3/21r4[vr2-2r(r·v)]U=dΨ~dt=-2·[1+zz2+r2+zr22(z2+r2)3/2]U=dΨ~dt=-2+wr-zv/2(z2+r2)3/2}.
Here, U is a function of the microlens position (r,z) and its velocity v.Because of the smallness of the effect, the action of all Milky Way stars will be taken into account linearly. This will be performed by integrating (73). However, this requires justification. In fact, we are dealing with extended sources; this is not taken into account by  (73). As we will see below, this may be important in our problem.Thus, we return to (73), which describes contribution of a single star. We assume that a star at the point (r,z) has velocity Vp(r,z). Then, we consider sum over all stars in the Milky Way, which must further be averaged with the mass density ρ(r,z). This enables us to pass to integration to yield an average GIM in domain W:
(74)〈Utot〉W=Gc2∫dz∫d2rρ(r,z)U(r,z,Vp(r,z)),
where we suppose that ρ vanishes outside a bounded domain. Also, we suppose that ρ(r) is sufficiently smooth so as we can use approximation of weak microlensing inside the domain where ρ(r)≠0.Taking into account the explicit form (73), it is easy to see that in a general case considerable contribution may be due to stars at large impact distances from the line of sight. Also, one can show that the singularities in (74) for small r are integrable and the integral is convergent. This allows us to avoid the question about exact value of the lower limit of r in the definition of domain W.It is interesting to compare GIM (74) with an average over all S+W events in case of a general continuous mass distribution ρ=ρ(r,z,t). This distribution must satisfy the continuity equation:
(75)∂ρ∂t+div(ρVp)=0.In this case, (72) should be integrated with the density ρ:
(76)Ψtot=Gc2∫dz∫d2rρ(r,z,t)Ψ~(r,z).Note that this formula forΨtot can be considered either as the average over all S+W events or simply as a shift in case of mass distribution ρ(r).The shift of remote source image in case of continuous mass distribution has been treated in [93]. These shifts appear to be of the order of 10-5÷10-6 radians, but this value cannot be observed from the solar system. The observable (at least in principle) value is a change of (76) with time:
(77)Utot*≡dΨtotdt=Gc2∫dz∫d2r∂ρ∂tΨ~(r,z).
Formula (77) is, in fact, an average image velocity over all S+W events:
(78)Utot*=〈Utot〉S+W.
One might think that this must be the same as (74), but it is not true. The reason is that operations of averaging and differentiation here do not commute. To compare (74) with (77), one must use the continuity equation (75) and avoid singular points at r=0 in Utot* [10].We haveUtot*=U(0)+U(1), where U(0) and U(1) are two-dimensional vectors, and the components of the first term (i=1,2):
(79)Ui(0)=Gc2∫dz∫d2rρ(Vp·∇Ψ~i),
are the same as the components of (74): U(0)=〈Utot〉W. The components of the second term can be written as
(80)Ui(1)=-Gc2∫0∞dz∫r>ε→0d2rdiv(VpρΨ~i).
After some calculation on account of the Ostrogradsky-Gauss theorem, we have [10]
(81)〈Utot〉S+W-〈Utot〉W≡U(1)=-4πGc2∫0∞dzρ(z)V⊥(z),
where ρ(z) stands for the density of microlensing masses along the line of sight, V⊥ is the transverse velocity of microlenses near the line of sight. Evidently, (81), differs from zero. For example, if an observer at rest is situated in the center of Milky Way (stationary mass distribution), we have Utot*=0 (see (77)), but GIM does not equal to zero.For an observer in the solar system this effect yields a nonzero GIM value leading to apparent rotation of extragalactic reference frame in the direction of the Milky Way rotation [10]. For sources near the Galactic plane, the effect amounts about ~10-8 arcseconds per year [10]. The effect will be even smaller for sources with higher Galactic latitudes. Therefore, in principle (if the corresponding accuracy be achieved!), the effect can be observed from the Solar system.


## 4.3.1.   Correlation Tensor of Image Shifts

In a general case of a complicated microlensing system, we need some statistical measure of image motion; possible characteristic is proposed below.Letδx describes a deviation of the point source image from unperturbed position. In case of one microlens, this value typically varies as r-1, where r is a distance of the microlens from the line of sight. Because of this, we expect that in case of uniform distribution of many microlenses usual correlation function and 〈δx2〉 diverge, and so we prefer to deal with a correlation tensor like that used in the turbulence theory [91]:
(64)Dij(τ)≡〈(δxi(t)-δxi(t-τ))(δxj(t)-δxj(t-τ))〉.
This tensor is symmetric and it has a diagonal form in a system where the source moves along one of coordinate axes, for example, along y1-axis. So a determination of Dij(τ) from observations may help to constrain a direction of the source motion.We will illustrate calculation ofDij in case of weak astrometric microlensing by N clumps having a uniform distribution over the lens plane. The system is described by the lens equation (60). We suppose that all the clumps are independent, and xn has a uniform distribution over the lens plane. To be precise, we assume that the spatial distribution function for every xn is
(65){S-1,xn∈Ω,0,xn∉Ω,
where S is the area of domain Ω. Further considerations deal with the limit of Ω that will be blown up to cover the whole plane assuming
(66)S(Ω)N=s(a,M)=const,N→∞,
where we take into account that the surface per one clump s(a,M) may be different for clumps of different mass and size. Density f(a,M,x)~|x|-4 goes to zero sufficiently fast as x→∞, so the clumps have a finite size and integral M=∫‍d2xf(a,M,x)<∞ is convergent. To save the space, we do not consider a distribution of parameters a,M assuming that these parameters are the same. This means that in the general case the results below must be averaged over possible a,M.According to (60), in case of weak microlensing the shift of the point source image is
(67)δx=∑n=1NRE,n2y-xn|y-xn|2+an2.For a straight line motion of a source, we puty1=vt,y2=const. Then, we have a diagonal form of Dij with
(68)D11(τ)=RE4s∫d2x[x1,-x1,-2+x22+a2-x1,+x1,+2+x22+a2]2,D22(τ)=RE4s∫d2x[x2x1,-2+x22+a2-x2x1,+2+x22+a2]2,x1,±=x1±yτ,yτ≡|v|τ2.Calculation of the integrals yields(69)D11(τ)=πRE4sF11(yτa),D22(τ)=πRE4sF22(yτa),
where
(70)F11(ξ)=ηln1+η1-η,F22(ξ)=1ηln1+η1-η-2,η≡ξξ2+1.


## 4.3.2. Statistical Effect of Gravitational Dragging of a Distant Source Image by Point Masses

In this section, we consider statistical effects in the motion of the GLS images due to gravitational field of a large number of moving stars. Here, we consider a rarefied collection of lensing masses, when the typical distances between the source and the microlens are much larger than the corresponding Einstein radius. For brevity, we call this “weak” microlensing; this must be distinguished from the well-known extragalactic weak lensing, that is, image distortions caused by the cosmic shear. Unlike the previous sections, here we rule out strong microlensing events that are accompanied by a considerable brightness amplification.In some situations it will be difficult to separate GIM from real proper motions, and this is important for accuracy of the fundamental reference frame based on extragalactic sources [83, 85, 86]. The gravitational field of point masses (which is essentially inhomogeneous) induces a stochastic motion of the source image in the reference frame of the source. On the other hand, smoothly distributed masses also induce a stochastic motion, but typical amplitude of the effect will be weaker. Indeed, we cannot take into account positions of all stars [83]; we only can work with their probability distribution leading to distribution of image motions of a source [8, 9]. The situation could be further complicated in presence of invisible DM clumps, if they exist on small scales.In case of the point masses, it is possible to estimate the probability distribution of image motionsu (i.e., angular velocities) of the distant source image microlensed by collection of a randomly moving point masses [8]:
(71)P(u)=u02π(u02+u2)3/2,
where u0 is a characteristic value depending upon the transverse velocity of point masses, their value, and their spatial density. It should be noted that the distribution (71) is obtained under the approximation of the weak microlensing [8], so that HAEs must be excluded.The distribution (71) has a slowly decreasing tail yielding an infinite estimate of velocity squared 〈u2〉. This tail is a formal consequence of pointness of the microlenses that makes possible very high GIMs. In fact this case is not covered by formula (71), which is valid for weak microlensing only. In case of smooth DM concentrations such velocities would be absent.If a bulk motion of microlensing masses is present, the average〈u〉 requires a special consideration. It has been pointed out in [8, 9] that stochastic GIM of an extragalactic source induced by stellar motions is accompanied with a systematic component (not a statistical average) which depends upon bulk velocity of microlensing stars, including stars that are far away from the line of sight. The typical velocity of this component remains nonzero even for a homogeneous distribution of microlensing masses which move in certain direction (all masses have the same collective velocity component). This needs some explanation because the average position of the unmoving source image (averaged over a very long time) must be unchanged and, evidently, the statistical average of image velocity must be zero. However, when dealing with the extragalactic reference frame, in reality, we do not deal with a very large observational time and/or with a sufficiently large number of sources to obtain a statistical average. In this case, because the strong microlensing events yielding the fast image motions are very rare, they in practice are not observable unless the number of simultaneously observed sources is of the order of millions or more. This is well known from theory and observations of MACHOs. Then the statistical average over set of all events is not a good characteristic of GIM; the most probable value of GIM (which is different from zero) would be a better measure.Otherwise, we can estimate probability for all necessary velocity intervals. Detailed approach on this way, however, involves consideration of strong lensing events with different extended source models. However, to have an order-of-magnitude estimate, we can avoid calculations of probability distributions, if we confine ourselves to the domain of weak microlensing events (“W”, for brevity) [9, 10]. Typical velocity which is a characteristic of such events is a statistical average over the domain W, with HAE being excluded. Note that these rare strong lensing events (S) essentially contribute to the usual statistical average velocity 〈u〉W+S (i.e., average over all the space of events S+W). If we exclude S-events (this is easy to do in real observations because of strong brightness amplification of a source), then we deal with “truncated” average velocity 〈u〉W over W events. This value is nonzero even in case of a homogeneous distribution of microlenses having a collective velocity. It is important to note that this value practically does not depend upon the exact definition of S-domain (i.e., upon a value of minimal impact parameter of microlensing masses with respect to the line of sight, which determines events “S”).Note that the “truncated” average would better correspond to realistic observations when the typical velocity is estimated as an arithmetic mean over all reference sources used in extragalactic International Celestial Reference Frame (ICRF). In a near future, the number of these sources will not exceed several hundreds. Therefore, in fact, we deal with most probable dispositions of lensing masses which typically do not include theS-events: it is well known that the probability of the such events is of the order of 10-6 per year in the Milky Way galaxy. Then, in ICRF considerations, it is reasonable to exclude such rare events that would hardly occur during the century.In case of microlensing of distant sources by Milky Way stars, we may neglect cosmological curvature. Evidently, this does not contradict to assumption that the radiation source is at the infinity. Let unperturbed light ray move from the infinitely distant radiation source in negative direction ofz-axis of Cartesian coordinates {x,y,z}, the observer being at the origin. Let position of microlensing point mass M be (r,z), where r=(x,y) is a two-dimensional vector in the transverse plane; that is, r=|r| is the impact distance of the unperturbed ray with respect to the mass. In the case of weak microlensing, r≫(mz)1/2, m=GM/c2. In this case,
(72)Ψ=mΨ~,Ψ~=-2rr2[1+z(z2+r2)-1/2],
where Ψ=(Ψ1,Ψ2) is a two-dimensional vector describing the source image angular shift [9, 85, 92].FurtherVp=(v,w) stands for velocity of the point microlens, w is the velocity component parallel to the line of sight, and v represents the transversal components. In virtue of (72), this leads to the source GIM that equals to mU (in radians per unit of time), where
(73)U=dΨ~dt=-2{wr-zv/2(z2+r2)3/21r4[vr2-2r(r·v)]U=dΨ~dt=-2·[1+zz2+r2+zr22(z2+r2)3/2]U=dΨ~dt=-2+wr-zv/2(z2+r2)3/2}.
Here, U is a function of the microlens position (r,z) and its velocity v.Because of the smallness of the effect, the action of all Milky Way stars will be taken into account linearly. This will be performed by integrating (73). However, this requires justification. In fact, we are dealing with extended sources; this is not taken into account by  (73). As we will see below, this may be important in our problem.Thus, we return to (73), which describes contribution of a single star. We assume that a star at the point (r,z) has velocity Vp(r,z). Then, we consider sum over all stars in the Milky Way, which must further be averaged with the mass density ρ(r,z). This enables us to pass to integration to yield an average GIM in domain W:
(74)〈Utot〉W=Gc2∫dz∫d2rρ(r,z)U(r,z,Vp(r,z)),
where we suppose that ρ vanishes outside a bounded domain. Also, we suppose that ρ(r) is sufficiently smooth so as we can use approximation of weak microlensing inside the domain where ρ(r)≠0.Taking into account the explicit form (73), it is easy to see that in a general case considerable contribution may be due to stars at large impact distances from the line of sight. Also, one can show that the singularities in (74) for small r are integrable and the integral is convergent. This allows us to avoid the question about exact value of the lower limit of r in the definition of domain W.It is interesting to compare GIM (74) with an average over all S+W events in case of a general continuous mass distribution ρ=ρ(r,z,t). This distribution must satisfy the continuity equation:
(75)∂ρ∂t+div(ρVp)=0.In this case, (72) should be integrated with the density ρ:
(76)Ψtot=Gc2∫dz∫d2rρ(r,z,t)Ψ~(r,z).Note that this formula forΨtot can be considered either as the average over all S+W events or simply as a shift in case of mass distribution ρ(r).The shift of remote source image in case of continuous mass distribution has been treated in [93]. These shifts appear to be of the order of 10-5÷10-6 radians, but this value cannot be observed from the solar system. The observable (at least in principle) value is a change of (76) with time:
(77)Utot*≡dΨtotdt=Gc2∫dz∫d2r∂ρ∂tΨ~(r,z).
Formula (77) is, in fact, an average image velocity over all S+W events:
(78)Utot*=〈Utot〉S+W.
One might think that this must be the same as (74), but it is not true. The reason is that operations of averaging and differentiation here do not commute. To compare (74) with (77), one must use the continuity equation (75) and avoid singular points at r=0 in Utot* [10].We haveUtot*=U(0)+U(1), where U(0) and U(1) are two-dimensional vectors, and the components of the first term (i=1,2):
(79)Ui(0)=Gc2∫dz∫d2rρ(Vp·∇Ψ~i),
are the same as the components of (74): U(0)=〈Utot〉W. The components of the second term can be written as
(80)Ui(1)=-Gc2∫0∞dz∫r>ε→0d2rdiv(VpρΨ~i).
After some calculation on account of the Ostrogradsky-Gauss theorem, we have [10]
(81)〈Utot〉S+W-〈Utot〉W≡U(1)=-4πGc2∫0∞dzρ(z)V⊥(z),
where ρ(z) stands for the density of microlensing masses along the line of sight, V⊥ is the transverse velocity of microlenses near the line of sight. Evidently, (81), differs from zero. For example, if an observer at rest is situated in the center of Milky Way (stationary mass distribution), we have Utot*=0 (see (77)), but GIM does not equal to zero.For an observer in the solar system this effect yields a nonzero GIM value leading to apparent rotation of extragalactic reference frame in the direction of the Milky Way rotation [10]. For sources near the Galactic plane, the effect amounts about ~10-8 arcseconds per year [10]. The effect will be even smaller for sources with higher Galactic latitudes. Therefore, in principle (if the corresponding accuracy be achieved!), the effect can be observed from the Solar system.


## 4.4. Resume of Section4

In this Section, we review main points of astrometric microlensing. We demonstrate effects due to a finite size of the source and extended lens as well; a special attention is paid to effects that may serve as signals of dark matter. Several examples of microlensing by one mass in presence of an external shear are given. This can be the case of microlensing of a remote source by a small mass (planet) in presence of a background gravitational field of a star or microlensing by a star in the field of a putative DM clump.Then, we proceed to some statistical subjects of astrometric microlensing. We introduced the correlation tensor as a statistical characteristic of a stochastic image motion and applied this to a “toy” model of weak microlensing by dark matter clumps.We have shown that gravitational field of foreground stars induces an image motion, which differs for a continuous and discrete mass density distribution. In the latter case, image performs random walks, as distinct from a regular motion in case of continuous matter, plus some additional motion due to the collective motion of the microlensing masses. This effect leads to apparent rotation of extragalactic reference frame in the direction of the Milky Way rotation. In principle, the effect can be observed, however, it is far beyond modern possibilities.



---

*Source: 906951-2012-09-25.xml*