Fast Computational Convolution Methods For Extended Source Effects In Microlensing Lightcurves

TL;DR

This paper introduces fast convolution-based algorithms for calculating microlensing light curves involving extended sources, enabling real-time analysis of large datasets to study stellar profiles and limb darkening.

Contribution

It develops a general formalism and analytic solutions for efficient computation of microlensing light curves with extended sources, including novel profiles like parabolic.

Findings

Analytic solutions achieve 0.5-3% accuracy.

Algorithms are significantly faster than existing routines.

Applicable to large-scale microlensing data analysis.

Abstract

Extended source effects can be seen in gravitational lensing events when sources cross critical lines. Those events probe the stellar intensity profile and could be used to measure limb darkening coefficients to test stellar model predictions. A data base of accurately measured stellar profiles will be needed to correctly subtract the stellar flux in planetary transient events. The amount of data that is being and will be produced in current and future microlensing surveys, from the space and the ground, requires algorithms that can quickly compute light curves for different source-lens configurations. Based on the convolution method we describe a general formalism to compute those curves for single lenses. We develop approximations in terms of quadratures of elliptic integrals that we integrate by solving the associated first order differential equations. We construct analytic solutions for a limb darkening and, for the first time, for a parabolic profile that are accurate at the $\sim 1-3\%$ and $0.5\%$ level, respectively. These solutions can be computed orders of magnitude faster than other integration routines. They can be implemented in pipelines processing large data sets to extract stellar parameters in real time.

Tables

Table 1. Table 1 Coefficients of the series expansion of the limb darkening profile out to 7th order.

	$f_0$	$e_0$	$f_1$	$e_1$	$f_2$	$e_2$
$a_0$	692583091200	1385166182400	400313026713600	400313026713600	144513002643609600	144513002643609600
$a_1$	1339448033263	3894814720117	-104737154144542	586396203521047	-41601097496530466	129881415373196381
$a_2$	312948159897	243184124641	118605243978153	130620128454965	25747524261728378	302061597305296826
$a_3$	72606018032	41938764464	-12149180921700	2732363448097	16403635906248033	17159057255113677
$a_4$	502361257	252468377	19577824287	233199164080	-760401005566596	-717449096541348
$a_5$	4109895	1997995	5187655332	840367513	683702230503	569981448915
$a_6$	79511	41279	7438983	4461547	29226547953	87184849284
$a_7$	835	527	111191	66815	120111823	82551251
$a_8$	308	308	963	655	1392911	931319
$a_9$			308	308	1107	799
$a_{10}$					308	308

Equations

$$\zeta =z+\sum _{i=1}^n\frac{m_i}{{\bar{z}}_i-\bar{z}},$$

$$z_{\mathrm{E}}=\sqrt{\frac{4GM{D}_{ls}{D}_l}{c^2{D}_s}}\mathrm{and}{\zeta }_{\mathrm{E}}=\frac{D_s}{D_l}{z}_{\mathrm{E}},$$

$$\frac{r_s}{{\zeta }_{\mathrm{E}}}\simeq 0.29\left(\frac{R}{400R_{\odot }}\right){\left[\left(\frac{D_l}{5\mathrm{Kpc}}\right)\left(\frac{5\mathrm{Kpc}}{D_{ls}}\right)\left(\frac{10\mathrm{Kpc}}{D_s}\right)\right]}^{1/2},$$

$${\mu }_j=\frac{1}{\det J}={{\left(1-\frac{\partial \zeta }{\partial \bar{z}}\overline{\frac{\partial \zeta }{\partial \bar{z}}}\right)}^{-1}|}_{z=z_j}.$$

$${\mu }_{\mathrm{t}}(\xi ,\eta )=\sum _{j=1}^k\left|{\mu }_j(\xi ,\eta )\right|=\sum _{j=1}^k\left|\mu (x_j(\xi ,\eta ),y_j(\xi ,\eta ))\right|.$$

$${\mu }_{\mathrm{ext}}({\xi }_0,{\eta }_0)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mu (\xi ,\eta )s({\xi }_0-\xi ,{\eta }_0-\eta )d\xi d\eta ={\int }_{-r_s}^{r_s}{\int }_{-r_s}^{r_s}\mu (\xi -{\xi }_0,\eta -{\eta }_0)s(\xi ,\eta )d\xi d\eta ,$$

$$V_s={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }s(\xi ,\eta )d\xi d\eta =1,$$

$$V_{\mu >1}={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }({\mu }_{\mathrm{ext}}(\xi ,\eta )-1)d\xi d\eta =V_s{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }({\mu }_{\mathrm{point}}(\xi ,\eta )-1)d\xi d\eta ,$$

$$\mu (r_0)=\frac{1}{2}\left(\frac{r_0}{\sqrt{r_0^2+4}}+\frac{\sqrt{r_0^2+4}}{r_0}\right).$$

$$V_{\mu >1}=2\pi {\int }_0^{\infty }(\mu (r)-1)rdr=\pi {\left[r\sqrt{4+r^2}-r^2\right]}_0^{\infty }=2\pi$$

$$V_{\mu >1}=2\pi (m_1+m_2)$$

$$c_{\mu }=\sum _i{\mu }_i^{+}-\sum _j{\mu }_j^{-}=const.$$

$$\sum _i{\mu }_{i,\mathrm{ext}}^{+}-\sum _j{\mu }_{j,\mathrm{ext}}^{-}={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{c}_{\mu }s(\xi ,\eta )d\xi d\eta =c_{\mu }$$

$$\sum _i{A}_i^{+}-\sum _j{A}_j^{-}=A_s.$$

$$\frac{I(r)}{I(0)}=1-u_1-u_2+u_1\sqrt{1-\frac{r^2}{r_s^2}}+u_2\left(1-\frac{r^2}{r_s^2}\right).$$

$$s_{\mathrm{star}}(r)=\frac{6(1-u_1-u_2)s_{\mathrm{disk}}(r)+4u_1{s}_{\mathrm{limb}}(r)+3u_2{s}_{\mathrm{para}}(r)}{6-2u_1-3u_2},$$

$$s_{\mathrm{disk}}(r)$$

$$s_{\mathrm{limb}}(r)$$

$$s_{\mathrm{para}}(r)$$

$${\mu }_{\mathrm{ext}}(r_s,r_0)={\int }_0^{2\pi }{\int }_0^{r_s}s(r,\phi )\mu (|\mathbf{r}-{\mathbf{r}}_0|)rdrd\phi ,$$

$${\mu }_{\mathrm{ext}}(r_s,r_0)={\int }_0^{r_s}k\sqrt{\frac{r}{r_0}}[(r-r_0{)}^2\Pi (n,k)+2K(k)]s(r)dr,$$

$$n=\frac{4r_{0}r}{(r_0+r{)}^2},k=\sqrt{\frac{4n}{4+(r_0-r{)}^2}},$$

$$\frac{d\mu (r,r_0)}{dr}=k\sqrt{\frac{r}{r_0}}[(r-r_0{)}^2\Pi (n,k)+2K(k)]s(r),$$

$${\mu }_{\mathrm{disk}}(r,r_0)=a_1(r)K(k)+a_2(r)E(k)+a_3(r)\Pi (n,k)$$

$$a_1(r_s)=\frac{k(r_s)}{\sqrt{r_0{r}_s}}\frac{(r_s^2-r_0^2)p(r_s)}{\pi r_s^2},a_2(r_s)=\frac{\sqrt{r_0{r}_s}}{k(r_s)}\frac{q(r_s)}{\pi r_s^2}\mathrm{and}{a}_3(r_s)=\frac{k(r_s)}{\sqrt{r_0{r}_s}}\frac{(r_0-r_s{)}^2(r_s^2+1)}{2\pi r_s^2}$$

$$p(r)=\frac{1}{8}(8+r_0^2-r^2)\mathrm{and}q(r)=2.$$

$$\underset{r_s\to 0}{\lim }{\mu }_{\mathrm{ext}}(r_s,r_0)=\frac{2+r_0^2}{r_0\sqrt{4+r_0^2}}.$$

$$\underset{r_0\to 0}{\lim }{\mu }_{\mathrm{ext}}(r_s,r_0)=\frac{\pi }{2}(a_1(r_s)+a_2(r_s)+a_3(r_s)).$$

$$m_{{\lambda }_1}-m_{{\lambda }_2}=-2.5\log \left(\frac{{\mu }_{{\lambda }_1}}{{\mu }_{{\lambda }_2}}\right),$$

$$m_{\mathrm{disk}}-m_{\mathrm{para}}$$

$$m_{\mathrm{disk}}-m_{\mathrm{limb}}$$

$$m_{\mathrm{disk}}-m_{\mathrm{point}}$$

Full text

1. Introduction

Limb darkening is the variation of the specific intensity of a star from its center to its limb. The main probes of the effect in distant stars have been the light curves of eclipsing-binary systems (see Zola et al. 2015). Optical interferometry has allowed direct imaging of those few stars whose angular diameter could be resolved (Hestroffer, 2001; Aufdenberg, Ludwig & Kervella, 2005). Exoplanetary transits are potentially very useful to determine stellar profiles although empirical values derived from data show poor agreement with atmospheric models (Howarth, 2011). Equally relevant is that these transits provide accurate measurements of exoplanet properties. To study the atmosphere of transiting planets accurate limb-darkening stellar models are needed to correctly subtract the stellar flux (, Morello et al.2017). Therefore, different and independent data of the limb-darkening effect are required to correctly determine the properties of planetary atmospheres and to test stellar models. Microlensing has become a powerful observational tool with different astrophysical applications since the observation of the first event (Mao, 2012). In particular, gravitational lensing is sensitive to the effect of limb darkening (Witt, 1995). Extended source effects can be seen when the source crosses a caustic; the parameters of the star intensity profile can be measured from these events and can be compared with the predictions of stellar models. They were first seen in high magnification single-lens events (Alcock et al., 1997; Jiang et al., 2004; Fouqué et al., 2010; , Zub, Cassan, Heyrovský et al.2011). These events, although rare, are important as they provide a way to measure the lens mass (An et al., 2002; Jiang et al., 2004). Albrow, Beaulieu, Caldwell et al. (1999) were the first to determine the limb-darkening coefficients of a K giant star using microlensing, the first ever for a star in the Galactic bulge. Afonso et al. (2000) measured the limb-darkening coefficients for a metal-poor A star, in five bands from I to V, in a binary lens event. Sources crossing caustic lines in the lens plane are more common in multiple lens systems and extended effects could be more easily seen (Afonso et al., 2000; Albrow, An, Beaulieu et al., 2001; Choi et al., 2012; Chung et al., 2017; , Han, Jung & Udalski2018). Satellite observations like Spitzer offer a different line of sight than ground telescopes. In Chung et al. (2017), the center of the source passed very close to the projected position of the lens as seen by Spitzer and the size of the star could be measured from a single point, although no extended source effects were seen from Earth. Like in transient events, the profile coefficients of stars measured by gravitational lensing are not always in agreement with theoretical expectations (Fouqué et al., 2010). The understanding of radiative transfer models needs to be improved to fit microlensing observations (Cassan et al., 2006), demonstrating the need for additional data.

The effect of the source extension on the magnification pattern can be readily computed by convolving the source profile with the gravitational amplification of a point source since the magnification of different sources are additive for a given lens. Early calculations were centered on constant profiles (Gould, 1994; Witt & Mao, 1994). Convolution with an extended source was usually applied in the vicinity of a single caustic (e.g. Schneider & Weiss 1987) to compute the shape of the light curve. Deconvolution has been used to recover the profile of the lensed source (, Grieger, Kayser & Schramm1991). Although convolution was also applied to compute gravitational magnifications (Witt, Mao & Schechter, 1995), no account to the numerical method was given. In the near future we can expect a significant increment in data and fast methods to obtain light-curves will be needed. Facilities like the future Wide Field Infrared Survey Telescope (WFIRST) satellite, currently under construction, will devote a significant fraction of its observing time to the Galactic bulge, with a detection of $\sim 50,000$ microlensing events being expected (Spergel et. al, 2015; Gaudi, Akeson, Anderson et al., 2019). The ESA Euclid satellite is also capable of high cadence observations on a wide field of view, and in two colors, and it will also be a useful instrument for microlensing. Being at a different L2 location than WFIRST, events that are simultaneously observed by both satellites, and possibly also from the ground, will be measured from different lines of sight, facilitating the detection of extended source effects and allowing the lens mass to be estimated directly from the data (Penny, Bachelet, Johnson et al., 2019). The photometric error of the WFIRST W149 filter would be $\sim 10^{-2}$ mag per exposure for a AB$\sim 21$mag star (Penny,Gaudi, Kerins et al., 2019), comparable or better than observations from the ground. Processing those forthcoming data sets would require algorithms that can obtain light curves of lensing events with an accuracy of $\sim 0.01$ mag or better. This is most efficiently done using analytic expressions instead of light curves by numerical methods.

In this article we show that convolution offers a novel technique to compute very efficiently extended source effects for single point-mass lenses. It is an exact and computationally less expensive method than integration of the image contour described in Witt & Mao (1994) and Gould & Gaucherel (1997). Our approximations are more general than those used in (Yoo, DePoy, Gal-Yam et al.2004) and Chung et al. (2017), are faster than currently available methods and can be easily implemented in data analysis pipelines for the forthcoming data sets. Briefly, in §2 we summarize the basic formalism and we present the convolution method used to derive analytic solutions; we also introduce the models of stellar limb-darkening that will be considered in the article. In §3 we present our estimates for single lens events, in terms of simple quadratures, that can be applied to any source profile. In §4 we compute the exact analytic expressions for a constant intensity profile and in §5 we derive approximations for a parabolic and limb-darkening profiles that are accurate at the level of $\leq 0.005$mag and ${{}_{<}\atop{}^{\sim}}0.02-0.03$mag, respectively. These analytic formulae can not be generalized to binary or multiple lens system and an analysis of these cases will be given elsewhere. Finally in §6 we summarize our results and present our conclusions.

2. Basic Results and Source Profiles

A full account of the theory of lensing can be found in (Schneider, Ehlers & Falco1992). The lens equation for point masses is more easily solved if expressed in terms of complex quantities (Witt, 1990). Let $\zeta=\xi+\rm i\eta$ denote the (complex) position of a point source in the $(\xi,\eta)$-source plane, $z=x+{\rm i}y$ its image in the $(x,y)$-lens plane and $z_{i}=x_{i}+{\rm i}y_{i}$ the positions of a field of $n$ point mass stars of mass $m_{i}=M_{i}/M$ also in the lens plane. For this configuration, the lens equation is given by

[TABLE]

where $\bar{z}$ is the complex conjugate of $z$. The quantities $z$ and $\zeta$ are in units of the Einstein radius of the deflector/lens plane and source plane, respectively, and stellar masses $m_{i}$ can be expressed in units of the total mass $M$ so that $\sum_{i}m_{i}=1$. The normalized (Einstein) units are given by

[TABLE]

where $D_{l}$, $D_{s}$, $D_{ls}$ are the distances to the lens, the source and the distance between the lens and the source, respectively. The solutions of eq. (1) are the image positions $z$ for each source position $\zeta$.

In this article we will express the stellar radii $r_{s}$ in units of the Einstein radius of the lens-source system. A star in the bulge of physical radius $R$ will have a radius $r_{s}$ in the source plane

[TABLE]

with $R_{\odot}$ the solar radius. The fiducial value $R\simeq 400R_{\odot}$ is characteristic of C-rich giant stars (, Paladini, van Belle, Aringer et al.2011, van Belle, Paladini, Aringer et al.2013) so stellar radii of bulge stars will typically be in the range $r_{s}\simeq 10^{-3}\zeta_{\rm E}-0.3\zeta_{\rm E}$.

The effect of gravitational lensing on a background source is to magnify its intensity. A field of $n$ point-mass lenses in the deflector plane will produce $k$ images of a given point source located at position $\zeta$ in the source plane. The magnification of one image located at $z_{j}$ is given by

[TABLE]

The total magnification of the background source is the sum of the absolute magnification of each image,

[TABLE]

Images and total magnification can be computed for each source location $(\xi,\eta)$ so the relative motion of the lens and source will produce a magnification pattern over the whole source plane area.

From eq. (5) it follows that the magnification of an arbitrary number of sources by the same set of lenses is additive. Once the magnification of a point-like star, $\mu(\xi,\eta)$, is known the magnification pattern for an extended star can be derived by convolving the point-like pattern with the extended source profile. If the center of a background star of radius $r_{s}$ is located at $(\xi_{0},\eta_{0})$ then its magnification is

[TABLE]

where $s(\xi,\eta)$ is the source profile. If the source is not lensed, then the magnification pattern will be $\mu(\xi,\eta)=1$ everywhere and the magnification will remain unchanged for any source profile. Then, it follows the normalization condition

[TABLE]

valid for any profile $s$.

2.1. The Conservation of Magnification

Fubini theorem guarantees that the magnification is conserved, irrespectively of the source profile. The theorem states that the volume enclosed by the convolved function is equal to the product of the single volume enclosed by the two functions provided the integral exists. Although the volume enclosed by the magnification is infinite, the volume enclosed by $\mu(\xi,\eta)-1$ is finite. The integral in eq.(6) remains valid if we replace $\mu(\xi,\eta)$ by $\mu(\xi,\eta)-1$ since the total magnification is $\mu(\xi,\eta)\geq 1$ everywhere (cf. Schneider 1984). We can write

[TABLE]

where the identity follows from Fubini theorem and is valid for a field of $n$ point-like lenses. From the normalization condition given in eq. (7) it follows that the enclosed volume above $\mu>1$ does not change for arbitrary source profiles and is equal to that of a point-like source, and then, finite. Then eq. (8) provides a self-consistency check to verify the accuracy of the estimated amplification $\mu_{\rm ext}$ for any source profile.

We can solve eq. (8) for simple cases. Let $r_{0}=(\xi_{0}^{2}+\eta_{0}^{2})^{1/2}$ denote the separation of the center of source to the projected position of the lens, in the source plane. For a single point mass lens, the magnification is (Refsdal, 1964)

[TABLE]

Introducing polar coordinates in eq. (8) by defining $r=(\xi^{2}+\eta^{2})^{1/2}$ we obtain

[TABLE]

For more complicated lens configurations, eq. (8) can not be solved analytically. For instance, the light curve of a close binary is very different from that of single lens located at the same distance. The differences are largest when the distance from the center of the source to the lens, $r_{0}$, is of the same order of magnitude or smaller than the projected separation of the binary lens $d_{\rm lens}$. On the contrary, when $r_{0}\gg d_{\rm lens}$ the magnification pattern is that of a single lens. We solved eq. (8) for several binary configurations and carry out the integration from the center of mass of the binary to a distance $r_{\max}$ ten times the maximum of the binary lens separation and the radius of the source star: $r_{\max}=10\,\max(r_{s},d_{\rm lens})$. Our numerical estimates showed that, as expected, the enclosed volume of the magnification $V_{\mu>1}$ remains constant and depends only on the total mass of the system. We obtained

[TABLE]

where $m_{1}$ and $m_{2}$ are the masses of the binaries in arbitrary units. If we express the mass of each component in units of the total mass of the system, then $V_{\mu>1}=2\pi$ again. We conjecture that for a general lensed point mass system the volume of the magnification pattern will always depend on the total mass of the system $V_{\mu>1}=2\pi\sum_{i}m_{i}$.

2.2. Signed Magnification

A direct consequence of the convolution integral for lens models (eq. 6) is that in certain areas where the sum of signed magnification for point sources is constant the relation is the same for extended sources. Let us assume that we have a caustic network with an area

[TABLE]

where $\mu_{i}^{+}$ is the magnification of the image of positive parity and $\mu_{i}^{-}$ the magnification of the image of negative parity. By applying eq.(6), for an extended source we can write

[TABLE]

since we may exchange the integral and sum of each magnification. In particular for a single point mass lens we have $c_{\mu}=1$ everywhere and for a binary lens we have $c_{\mu}=1$ inside the caustic(s) (see Witt & Mao 1995, Witt & Mao 2000, Hunter & Evans 2001). Let $A_{i}^{+}$ be the area of the lensed image of positive parity and $A_{j}^{-}$ the area of the lensed image of negative parity. For the single point mass and for the binary lens if the area of the source $A_{s}$ is located inside the caustic, it follows that

[TABLE]

When computing light curves for extended sources by a contour method (cf. Gould & Gaucherel 1997), this equation provides a useful check on the accuracy of the numerical integration over the extended sources. Alternatively, it can be used to reduce the integration on the images of positive parity since this would be sufficient to obtain the total magnification.

2.3. Source Profiles

The intensity of a star can be very well modeled with a general limb-darkening profile of the form (e.g. Allen 1973, Claret, Díaz-Cordovés & Giménez1995)

[TABLE]

The coefficients $u_{1}$ and $u_{2}$ depend on the observed (wavelength) band of the star and may differ quite strongly (for instance, see Fouqué et al. 2010). Those quantities are not independent; assuming that the intensity profile (eq. 15) is everywhere positive and monotonically decreasing from the center to the limb, these coefficients verify $u_{1}>0$, $u_{1}+u_{2}<1$ and $2u_{2}+u_{1}>0$ (Kipping, 2013). The source profile normalized according to eq. (7) is

[TABLE]

with

[TABLE]

where each of these three profiles obeys individually the normalization condition of eq. (7). In the following sections we will use these profiles to illustrate our analysis.

3. The single point mass lens

For a single point mass lens, the magnification depends only on the distance between the source and the projected position of the lens (eq. 9). If the latter is at the origin of coordinates and $\bm{r}_{0}=(\xi_{0},\eta_{0})$, $\bm{r}=(\xi,\eta)$ denote the coordinates of the center of the source and of an arbitrary point on its surface, respectively, then eq. (6) can be written in polar coordinates as

[TABLE]

being $|\bm{r}-\bm{r}_{0}|=(r^{2}+r_{0}^{2}-2rr_{0}\cos\varphi)^{1/2}$ and $r_{s}$ the radius of the star. Computing the double integral on a grid gives profiles with rigging whose amplitud scales as $N_{p}^{-1/2}$, being $N_{p}$ the number of points on the grid. The rigging affects only those trajectories where the lens moves across the surface of the star and to improve the accuracy requires to increase the number of points, increasing the computational time.

Smoother light curves are obtained when the stellar intensity profiles are spherically symmetric. In this case, (see Appendix A and also Witt 1995, Heyrovský 2003) the angular part can be integrated to obtain the amplification in terms of a single quadrature of elliptical functions

[TABLE]

with

[TABLE]

In eq. (21) $K(k)$, $\Pi(n,k)$ denote the complete elliptic integral of the first and third kind. The integral has a singularity at $k=1$ or $r=r_{0}$ so that care needs to be taken to obtain accurate results. Gauss-Legendre integration over $10^{3}$ points introduces rigging at the 2-3% level. We used the Press et al. (2002) routines for improper integrals *qtrap.f * and *midpnt.f * for more accurate results. To reduce the amplitude of the rigging below the 0.1% level we run subdivisions up to JMAX=20 with an accuracy of EPS=$10^{-5}$. This high precision slows down the code compared with other integration routines.

In Fig. 1 we represent the magnification of a point source and of the three extended source profiles given in eqs. (17-19), computed with routines for improper integrals. Depending on the impact parameter and stellar radius the magnification could be larger or smaller than that of a point source. In all cases, the parabolic profile was the closest to the point source and the constant profile was the furthest. For instance, in plots (a,e,f) curves at the peak correspond, from top to bottom, to the amplification of a point source (black line) and extended sources with constant (blue), limb-darkening (red) and parabolic (green) profiles. In (b,c,d), the color code is the same but the ordering of the magnification curves is the reverse. In the x-axis, $v_{s}$ is the relative velocity of the source and the lens, and $t$ is the time of observation. The distance $r_{0}$ from the center of the source to the projected position of the lens is given by $r_{0}=[(v_{s}t)^{2}+b^{2}]^{1/2}$, being $b$ the impact parameter. Distances and stellar radii are given in units of the Einstein radius on the source plane, $\zeta_{\rm E}$. Plots (a-c) correspond an impact parameter $b=10^{-3}$ and stellar radii $r_{s}=10^{-3},10^{-2},0.1$; plots (d-f) correspond to impact parameters $b=(0.01,0.1,1)$ and stellar radius $r_{s}=0.3$. The largest differences between magnification of point and extended sources with any of the three profiles occur when the lens is within the projected surface of the star. Once the lens is outside the magnification decreases, extended source effects are diluted and the amplification coincides with that of a point source.

4. Exact Analytic solutions.

Faster numerical methods require the full integration of eq. (20). Exact solutions for particular cases can be found by solving the differential equation associated to eq. (21)

[TABLE]

In Appendix B we show that this differential equation can be solved exactly for the constant intensity profile $s(r)=s_{\rm disk}(r)=1/\pi r_{s}^{2}$. The solution is given in terms of elliptic functions and has the form

[TABLE]

with

[TABLE]

where

[TABLE]

The coefficients of eqs. (25, 26) were first derived by Witt & Mao (1994). Elliptic integrals can be evaluated very efficiently (for instance, with the routines given in Press et al. 2002) and this solution can be computed very quickly since no integration is required.

Fig 1 demonstrates that the gravitational magnification of all profiles converge when the source is far from the lens. Eventually, the parabolic and limb darkening profiles match smoothly the disc profile given by eqs. (24-26). We will assume that a general solution of eq. (23) can be represented by the same functional expression but with different coefficients or possibly functions $(a_{1},a_{2},a_{3})$, to be determined. Introducing this ansatz in eq. (24) gives a set of three coupled differential equations for the coefficients $a_{i}(r)$ that needs to be solved for each specific profile. To fix the constants of integration, the solution has to converge to some limiting cases. For instance, the magnification reduces to that of a point source when $r_{s}\rightarrow 0$

[TABLE]

The differential equations for the coefficients given in Appendix B are rather complicated. For the limb and parabolic case we found solutions only at two specific configurations, when star and lens are perfectly aligned, $r_{0}=0$, and when the lens is at the edge of the star, $r_{0}=r_{s}$. When $r_{0}=0$ eq. (21) becomes

[TABLE]

If $r_{0}=0$ or $r_{0}=r_{s}$ the elliptic integral simplify considerably (see eqs. D7, D8) and we derived simple series expansions that converge rather quickly when $r_{s}\ll 1$, the most common case. These expressions are given in Appendix B.

The results at these two locations, although limited, are still informative. Since the difference between profiles is largest when source and lens are perfectly aligned, we can derive an upper bound on the contribution of each profile to the overall magnification. Second, when the lens is on the edge of the star extended source effects start to be noticeable. If only a few (or just one) data point are available it is useful to have an accurate estimate of the difference between profiles at physically relevant locations. The differences between two source profiles are given by

[TABLE]

provided that they have the same intensity at wavelengths $\lambda_{1}$ and $\lambda_{2}$. Otherwise, an additional factor needs to be included which result in a constant magnitude between the same source at different wavelengths. In the limit $r_{0}=0$, $\mu\propto r_{s}^{-1}$ and, by taking ratios, the stellar radius cancels out. The maximum possible magnitude difference between the three extended profiles is independent of the radius. From eqs. (B23, B31, B35) we have

[TABLE]

The magnitude difference with respect to a point source can not be given since when perfectly aligned with the lens, the amplification diverges.

In Fig. 2 we represent $m_{\rm disk}-m_{\rm limb}$ (solid blue line) and $m_{\rm disk}-m_{\rm para}$ (dashed red line) for different stellar radii. The vertical dot-dashed lines correspond to $r_{0}=r_{s}$. When $b\approx 0$ and $r_{s}\gg b$, then $\Delta m$ are very close to the limits given in eqs. (30). The magnitude difference decreases with increasing distance and reaches a minimum at $r_{0}=r_{s}$. The series expansion given in Appendix B can be used to estimate the differences between profiles. From eqs. (B28, B32,B36) we have

[TABLE]

also in agreement with the results shown in Fig 2. Notice that at a distance $r_{0}=r_{s}$ the magnification for the limb darkening and parabolic profiles is smaller than that of a constant disk since for the latter a larger fraction of the stellar flux is closer to the lens.

When the effect of an extended source is measured from very few or just one single point (as in Chung et al. 2017) one can not expect to measure the relative contributions of different profiles with significant accuracy. The estimates given in eqs. (30, 4), together with the inequalities given in §2.3 could still be used to test if the stellar profile deviates from a constant even though the data could not determine the coefficients $u_{1}$ and $u_{2}$ of eq. (16).

5. Approximate analytic solutions.

The differential equations obtained by the technique described in section §4 and Appendix B are very difficult to solve for a generic profile. Approximate analytic solutions can be obtained by a similar, albeit more simplified, technique. Fig. 1 demonstrates that when the distance from the center of a source to the projected position of a lens is larger than a few stellar radii all light curves match the constant stellar profile. Since the disk solution is analytic, we are only required to construct approximate and accurate solutions in the range where the magnifications deviate, approximately $r_{0}\leq 3r_{s}$. When the stellar profiles depend only on the radial distance $r$, we can compute the amplification of extended sources by Taylor expanding around the origin the magnification given in eq. (9)

[TABLE]

This expansion converges for $|r|<2$ (in units of the Einstein radius) and can not be applied to larger distances. The integration of eq. (20) can be carried out term by term giving

[TABLE]

where $c_{n}$ are the coefficients in eq. (32). If we express the impact parameter in units of the source radius $r_{0}=ur_{s}$ and set the integration variable to $t=r/r_{s}$ then

[TABLE]

with

[TABLE]

and $\hat{s}(t)=s(r_{s}t)\pi r_{s}^{2}$. In Appendix C we demonstrate that the integrations in eq. (35) can be expressed in terms of elliptic integrals and pairs of functions $e_{n}(t),f_{n}(t)$ and take the form

[TABLE]

Our treatment provides analytical approximations for the disk, parabolic and limb profiles in terms of elliptic integrals of the first and second kind. The solutions given in Appendix C have to be evaluated at $t=1$, at the radius of the star. The dominant contribution is $n=0$ and higher orders contribute more with increasing source mass and distance to the lens. For completeness we present the results up to $n=2$.

• •

For the disk profile $\hat{s}_{\rm disk}(t)=1$ we have

[TABLE]

• •

For the parabolic profile $\hat{s}_{\rm para}(t)=2(1-t^{2})$ we obtain

[TABLE]

• •

For the limb-darkening profile $\hat{s}_{\rm limb}(t)=(3/2)(1-t^{2})^{1/2}$ the solution can not be computed in closed form. We also need to expand the profile in Taylor series and find a solution at each order. The technique quickly produces cumbersome expressions. Increasing the order in the Taylor expansion improves the accuracy at $r_{0}\simeq b$, but the resulting series diverge faster. We found that the results up to the 7th order expansion provided the most accurate approximations in the range $r_{0}=(0,2.5r_{s})$ after which the solution can be matched to the exact disk profile with minimal error. In this case,

[TABLE]

The coefficients $(a_{0},\dots,a_{10})$ are given in Table 1. A routine to compute these functions is available upon request.

In Fig. 3 we plot the magnitude difference $\Delta m=m_{\rm numerical}-m_{\rm analytic}$ between the analytic approximations of eqs. (• ‣ 5) and eqs. (• ‣ 5) with respect to the numerical integration of eq. (21). The dashed (red) and solid (blue) lines correspond to limb and parabolic profiles. The limb analytic profile is computed out to $r_{0}=2.5r_{s}$ and beyond that distance it is matched to the exact constant disc expression of eqs. (24-26). In (a) the small wiggle of amplitude $\sim 10^{-3}$mag is due to instabilities on the numerical solution. Stellar radii and the impact parameter are given in units of the Einstein radius.

The zero order term in the Taylor expansion is always the dominant. For sources with radii $r_{s}\leq 0.01$, the first order contributes with 0.1% to the total magnification and the second order is negligible. For a star of $r_{s}=0.3$, the first and second order term contribution is about 30% and a 3% at large distances. The differences between the exact and approximated parabolic profiles (dashed red line) is always negligible and only for $r_{s}=0.3$ and at distances $v_{s}t\geq 4r_{s}$ the error approaches $0.01$mag. The amplitude is smaller than the numerical result indicating that terms higher than the 2nd order in the expansion of eq. (33) are important. However, at those distances the exact disk solution is very accurate and, like in the limb darkening case, is simpler to switch to it instead of computing terms like $I_{3}$ or higher. In Fig. 3c the thin dot-dashed red line shows the analytic solution matched with the exact disc solution. We chose to match both solutions at $r_{0}=3r_{s}$, which reduces the error in magnitude between the approximate and exact solutions to less than $5\times 10^{-3}$mag.

The largest magnitude difference in the limb darkening case occurs at $r_{0}\approx r_{s}$. The error was always smaller than $0.02$ mag except in the narrow range $(0.27,0.30)$ where it was less than $0.03$ mag. Comparing with the differences between profiles presented in Fig 2 those errors are 5-30% of the difference between disk, parabolic and limb profiles, proving that our approximations are sensitive to the different contributions. In the case of the limb profile, we overestimate the correct amplification and, consequently, $u_{1}$ in eq. (15) will be underestimated. Once the range of parameter space that fits a light curve has been identified, the bias can be corrected by fitting numerically integrated profiles.

To summarize, our approximations give the effect of an extended source in a microlensing event with errors $<3\%$ in the range $r_{s}{{}_{<}\atop{}^{\sim}}0.3$ and are free from instabilities that affect numerical integration. Computing light curves by solving eq. (20) on a circular grid of $\sim 7,800$ points, eq. (21) using the Gauss-Legendre on a linear grid of $10^{3}$ points, the same integral using the qtrap.f and midpnt.f routines was $1.1\times 10^{3}$, $430$ and $11\times 10^{3}$ times slower than using the analytic approximations presented here. Our results generalize those of (Yoo, DePoy, Gal-Yam et al.2004) and the technique can be extended to other profiles with even powers of the stellar radius.

6. Conclusions

We have developed an analytical formalism to compute the light curves of extended sources gravitationally amplified by single lenses. For the limb darkening profile we have obtained more accurate analytic approximations than those currently being used in the literature. We present even more accurate approximations for a parabolic profile. To our knowledge, these are the first approximations in the literature valid for this profile. To generate our approximations we first showed that eq. (20) can be reduced to a single quadrature involving elliptic integrals (eq. 21). We integrated this expression exactly for a constant stellar intensity profile by solving its associated first order differential equation, given by eq. (23), using the ansatz of eq. (24). Applying this formalism to the limb-darkening and parabolic profiles we obtained accurate expressions for two specific configurations: when the lens and source were perfectly aligned and when the lens was on the boundary of the source. These solutions were expressed in terms of the stellar radius $r_{s}$ and they converge very quickly for of $r_{s}\ll r_{s}$.

Next, we constructed approximations by Taylor expanding the magnification out to second order. These approximations were valid for lens-source distances $r_{0}\sim 3r_{s}$. For the limb darkening profile, our approximations were accurate to the 1% level, except for a small interval when the lens is at the edge of the stellar disk, that was only better than a 2-3%. The net effect is that the analytic solution will underestimate the $u_{1}$ coefficient in eq. (16). This shortcoming is not very relevant since once extended source effects have been identified, one can fit the data to light curves computed with more accurate numerical methods but centered around the parameter space selected by the analytic approximations. For the parabolic profile, this solution was almost identical to the numerical integration for a wide range of parameter space and only for very large stars ($r_{s}\geq 0.3$) and at large distances ($r_{0}\geq 3.0$) the error approached 1% (see Fig.3). For those configurations, the parabolic profile has already converged to the constant disc solution and the error is negligible, smaller than the expected photometric errors of forthcoming satellites.

In the next decade, WFIRST, Euclid and telescopes from the ground will greatly increase the measured number of microlensing events. Although their primary aim will be to search for exoplanets and characterize their properties, they will also probe stellar intensity profiles. Measuring limb darkening coefficients at different wavelengths will help to improve the numerical models of stellar atmospheres, a necessary step to study planetary atmospheres in transient events and correctly subtract the stellar flux. The data analyses will be helped by faster algorithms using analytically generated light curves.

Acknowledgments

We thank Shude Mao for comments on the manuscript. F. A.-B. acknowledges financial support from Grants No. PGC2018-096038-B-I00 (MINECO/FEDER) and No. SA083P17 from the Junta de Castilla y León.

Appendix A A. Angular integration of the magnification equation.

In eq. (20) the integration of the angular variable $\varphi$ is

[TABLE]

where $R^{2}=r^{2}+r_{0}^{2}$. The change of variables $t=\cos\varphi$ yields to

[TABLE]

where $a=(4+r^{2}+r_{0}^{2})/2rr_{0}$ and $b=(r^{2}+r_{0}^{2})/2rr_{0}$. Using eqs. (253.00,253.02) of Byrd & Friedmann (1971) we finally obtain

[TABLE]

with $k^{2}=[2(a-b)]/[(a-1)(b+1)]$. This expression gives eq. (21).

Appendix B B. Differential Equations of the Coefficients

We will integrate eq. (21) by assuming that the solution of the associated differential equation (eq. 23) has the functional form

[TABLE]

Using the properties of the complete elliptic functions given in Appendix D, we obtain the following three differential equations for the coefficients $a_{i}(r)$

[TABLE]

The function $k=k(r)$ is defined in eq. (22), $k{\prime}=dk/dr$ and

[TABLE]

Eq. (B4) is not coupled to the other two equations and is the simplest to solve. The solution of the homogenous equation, $a_{3,H}(r)$ can be found by integration

[TABLE]

where $c_{3}$ is a constant. A solution of the inhomogeneous equation can be found by variation of the constant. We write

[TABLE]

By inserting this expression in eq. (B4) we obtain $c_{3}^{\prime}(r)=4rs(r)$ whose solution is

[TABLE]

The solutions for the three different source profiles of eqs. (17-19) are presented below.

The coupled differential eqs. (B2,B3) are more difficult to solve. We decoupled the two equations by differentiating again and the solutions of the homogeneous equations are given in terms of elliptic integrals. It is easy to verify that two independent solutions of the homogeneous eqs. (B2) and (B3) are: (1) $a_{11}=E(k(r))$, $a_{21}=-K(k(r))$ and (2) $a_{12}=(E(\hat{k})-K(\hat{k}))$, $a_{22}=K(\hat{k})$ with $\hat{k}=\sqrt{1-k^{2}}$ (cf. Whittaker & Watson 1915). The general homogenous solution is

[TABLE]

The solution of the inhomogeneous equation can be obtained by variation of parameters $c_{1}$ and $c_{2}$. We write

[TABLE]

We insert this ansatz into eqs. (B2,B3) to obtain

[TABLE]

This equation can be readily solved for $c_{1}^{\prime}$ and $c_{2}^{\prime}$ and then integrated. Finally we can write

[TABLE]

since by the Legendre relation (see eq.(D6)) we have $a_{11}a_{22}-a_{12}a_{21}=\pi/2$. Then, from eq. (B) we obtain

[TABLE]

Introducing this results into eq. (B1) we obtain

[TABLE]

Since the integrals in eqs. (B) are rather involved we will solve them by separating variables $r$ and $r_{0}$ in a series expansion. We define

[TABLE]

with

[TABLE]

By differentiating eq. (B15) and after some algebra we obtain two polynomial equations of the form

[TABLE]

with

[TABLE]

Since the solution of the homogenous differential equation has the form

[TABLE]

one obtains first an integral for $p_{0}(r)$ and $q_{0}(r)$ by variation of the constants $c_{1}$ and $c_{2}$. The higher order $p_{2n}(r)$ and $q_{2n}(r)$ can be then obtained by recursive integration and by using the homogenous solution. Further if $h_{1}(r)$ and $h_{2}(r)$ have even power then the solutions $p(r)$ and $q(r)$ of the inhomogeneous equations must have even power as well as is the case of the three source profiles discussed here.

Solving the system of eqs. (B2-B4) is rather complicated. Only for a constant disk profile it is possible to obtain a simple exact solution, given below. For other profiles, simple solutions exist for specific configurations, when the approximations given in eqs. (D7,D8) are valid. The general (and cumbersome) case can be solved in terms of power series expansions (eq. B16) but we shall not discuss it here. We defer to Appendix C the derivation of accurate approximations for the limb and parabolic profiles.

B.1. Disk of Constant Surface Brightness

The profile of a disk of constant surface brightness is simply $s(r)=1/(\pi r_{s}^{2})$ and

[TABLE]

where the constant $C_{3}$ needs to be determined by adequate limiting conditions. For $a_{3}(r)$ we obtain

[TABLE]

To fix the constant of integration, we take the the limit $r_{0}\approx 0$. In this limit, the r.h.s of eq. (B14)) can be transformed into a quadrature over elementary functions. When $r_{s}\ll 4$, $k\approx\sqrt{n}\approx 0$ and we can apply the approximation given in eq. (D8). In the limit $r_{0}=0$ the elliptic integrals are $K(0)=\Pi(0,0)=E(0)=\pi/2$ and eq. (21) reduces to

[TABLE]

For a constant surface brightness

[TABLE]

Comparing this series with the integration of the l.h.s. of eq. (B14) given in eqs. (B15-B18) fixes the integration constant to the value $C_{3}=2$. Inserting the polynomials $p(r)$ and $q(r)$ into eqs. (B17) triggers rather simple results for $p_{0}(r)$ and $q_{0}(r)$. Further integration yields $p_{2}(r)=1/8$ and $q_{2}(r)=0$ and the series terminates. The final result is

[TABLE]

It is straightforward to verify that $p(r)$ and $q(r)$ satisfy eqs. (B17) for a disk profile. Finally, the full solution is

[TABLE]

When $r_{s},r_{0}\ll 1$ (in units of the Einstein radius) these expressions can be further simplified. In this limit we also have $k\approx\sqrt{n}$ and eq. (D8) allow us to write

[TABLE]

and for the total magnification one obtains

[TABLE]

where $u=r_{0}/r_{s}$. This result coincides with those presented in eqs. (36, 37) with $n=0$. The derivation is described in Appendix C.

Extended source effects start to become important when the projected position of the lens intersects the source disk, i.e., when $r_{0}\simeq r_{s}$. Using the properties of the elliptic integrals given in eq. (D7) we can write

[TABLE]

Similar limits at $r_{0}\simeq 0$ and $r_{0}=r_{s}$ can be found for the other profiles as we shall see below.

B.2. Limb Darkening Profile

For the limb darkening profile (eq. 18) we have

[TABLE]

with $C_{3}$ a constant to be determined. Then,

[TABLE]

From eq. (B22) we can compute the amplification when the lens and source are perfectly aligned,

[TABLE]

with $k_{2}=r_{s}/\sqrt{4+r_{s}^{2}}$. In this case, comparison with eq. (B30) gives the limit $C_{3}=0$, that physically corresponds to the source profile vanishing at the edge of the star.

Using the simplifying properties of elliptic integrals given in eq. (D7) for $r_{0}=r_{s}$ we obtain

[TABLE]

Since most lensing events will have $r_{s}\ll 1$ (see the estimate given in eq. (3)) this series converges very quickly.

B.3. Parabolic Source Profile

For the parabolic brightness profile of eq. (19) we have

[TABLE]

and

[TABLE]

The constant $C_{3}$ can be determined using the solution at $r_{0}=0$. In this case the amplification is

[TABLE]

that, like in the limb case, yields $C_{3}=0$. Using the simplifying properties of elliptic function given by eq. (D7) when $r_{0}=r_{s}$, we can derive the following series expansion

[TABLE]

Again, as $r_{s}\ll 1$ this series converges very quickly.

Appendix C C. Series Expansion of the Single Lens Magnification

The integration of the angular part of eq. (35) can be expressed in terms of elliptic integrals. If we denote $k_{1}=2\sqrt{ut}/(u+t)$, the first three terms in the series expansion become

[TABLE]

Similar combinations of complete elliptical integrals of the first, $K(k_{1})$, and second kind, $E(k_{1})$, with higher order polynomials are obtained for $n>2$. We will compute approximations up to $n=2$ since above this order the contribution to the light curve is at most a 1% (see Fig. 3).

The indefinite integrals can be evaluated assuming the result has the form

[TABLE]

By differentiating eqs. (C1, C2) the coefficients verify the following differential equations

[TABLE]

where the prime denotes the derivative with respect to the variable $t$, and for $n=0,1,2$ we have

[TABLE]

Then, solving eqs. (C3, C4) for $n=0,1,2$ is equivalent to integrate eqs. (C1). To simplify the differential equations of the coefficients, we define the new functions $e_{n},f_{n}$ as $a_{n}(t,u)=(t-u)f_{n}(t)$ and $b_{n}(t,u)=(t+u)e_{n}(t)$. Using the identities

[TABLE]

we obtain

[TABLE]

These equations can be combined to obtain a second order differential equation for $e_{n}$

[TABLE]

The r. h. s. of the differential equation $Q_{n}(t)$ results always into a polynomial of even power as long as the source profile $\hat{s}(t)$ has even power. We have

[TABLE]

Now we are looking for solutions of $e(t)$ if $Q(t)$ is simply given by $Q(t)=t^{2m}$. This differential equation has trivial solutions of a form of a polynomial $e(t)=P_{2m,u}(t)=\sum_{k=0}^{m}p_{2k}t^{2k}$ of the same power. Inserting the polynomial into eq. (C9) gives

[TABLE]

It is now easy to verify that $P_{2m,u}(t)$ obeys the recursion relation

[TABLE]

with

[TABLE]

and so on. In general we can now construct the solution $e_{n}(t)$ for any given $Q_{n}(t)$. The solution has the form

[TABLE]

In summary, $f_{n}(t)$ and $e_{n}(t)$ are polynomials of degree $2n$ in $t$ and $u$ for a disk and of degree $2n+2$ for a parabolic profile, respectively. If we expand the limb darkening profile until degree 10 the resulting $f_{0}$ and $e_{0}$ have powers of 10 in $t$ and $u$ respectively. Each higher order increase the power by 2.

Since we have constructed now a solution $e_{n}(t)$ and $f_{n}(t)$ we are able to write down for each source profile a set of solutions. For example, for $I_{0}$ we can write

[TABLE]

Let us remark that the lower limit of the integral always vanishes so we need to pay attention only to the upper limit. For each source profile

[TABLE]

Evaluating these expressions at the surface of the star, $t=1$, we obtain the results given in the main text. We note that the integral has the form

[TABLE]

with

[TABLE]

Note that in the main text we have defined $e_{n}(u)\equiv e_{n}(t=1)$ and $f_{n}(u)\equiv f_{n}(t=1)$.

Appendix D D. Functional Relations of Complete Elliptic Integrals

The following properties of the first, second and third complete elliptical integrals (see (Byrd & Friedmann, 1971))

[TABLE]

where used to derive eqs. (B2), (B3) and (B4). Also, the following relations

[TABLE]

and the Legendre relation

[TABLE]

with $\hat{k}=\sqrt{1-k^{2}}$, were also used. In the limiting case with $r\rightarrow r_{0}$ we have the following relations

[TABLE]

where we used the definition of $n$ and $k$ of eq.(22). For small values $r,r_{0}\ll 1$ we have $k\approx\sqrt{n}$ so that we can make use of the relation

[TABLE]