Fast computation of quadrupole and hexadecapole approximations in microlensing with a single point-source evaluation

TL;DR

This paper introduces a faster method for computing quadrupole and hexadecapole approximations in gravitational microlensing, significantly reducing computation time and enabling more efficient analysis of large datasets.

Contribution

It presents an optimized computational approach that accelerates quadrupole and hexadecapole calculations, making microlensing modeling more efficient and accessible.

Findings

Quadrupole approximation is about twice as slow as point-source evaluation.

Hexadecapole approximation is approximately four times faster than previous methods.

Open-source Python routines are provided for practical use.

Abstract

The exoplanet detection rate from gravitational microlensing has grown significantly in recent years thanks to a great enhancement of resources and improved observational strategy. Current observatories include ground-based wide-field and/or robotic world-wide networks of telescopes, as well as space-based observatories such as satellites Spitzer or Kepler/K2. This results in a large quantity of data to be processed and analyzed, which is a challenge for modeling codes because of the complexity of the parameter space to be explored, and the intensive computations required to evaluate the models. In this work, I present a method that allows to compute the quadrupole and hexadecapole approximation of the finite-source magnification with more efficiency that previously available codes, with routines about x6 and x4 faster respectively. The quadrupole takes just about twice the time of a point-source evaluation, which advocates for generalizing its use to large portion of the light curves. The corresponding routines are available as open-source python codes.

Figures

Tables

Table 1. Table 1: Coefficients c n − p , n k superscript subscript 𝑐 𝑛 𝑝 𝑛 𝑘 c_{n-p,n}^{k} (with k 𝑘 k referring to column W k subscript 𝑊 𝑘 W_{k} ) as defined in ( 16 ), for orders 2 ≤ p ≤ 5 2 𝑝 5 2\leq p\leq 5 . Since the formal expression of Q i , j subscript 𝑄 𝑖 𝑗 Q_{i,j} is the same as Q j , i subscript 𝑄 𝑗 𝑖 Q_{j,i} (but Q i , j ≠ Q j , i subscript 𝑄 𝑖 𝑗 subscript 𝑄 𝑗 𝑖 Q_{i,j}\neq Q_{j,i} ), interchanging all indexes ( i , j ) 𝑖 𝑗 (i,j) appearing in the expressions of Q p − n , n subscript 𝑄 𝑝 𝑛 𝑛 Q_{p-n,n} gives the expression of Q n , p − n subscript 𝑄 𝑛 𝑝 𝑛 Q_{n,p-n} , so it is enough to provide c n − p , n k superscript subscript 𝑐 𝑛 𝑝 𝑛 𝑘 c_{n-p,n}^{k} for n ≤ ⌊ p / 2 ⌋ 𝑛 𝑝 2 n\leq\lfloor p/2\rfloor .

	$W_3$	$W_4$	$W_5$	$W_6$
$Q_{2,0}$	$a_{1,0}^2$	—	—	—
$Q_{1,1}$	$a_{1,0}{a}_{0,1}$	—	—	—
$Q_{3,0}$	$3a_{1,0}{a}_{2,0}$	$a_{1,0}^3$	—	—
$Q_{2,1}$	$2a_{1,0}{a}_{1,1}+a_{0,1}{a}_{2,0}$	$a_{1,0}^2{a}_{0,1}$	—	—
$Q_{4,0}$	$4a_{1,0}{a}_{3,0}+3a_{2,0}^2$	$6a_{1,0}^2{a}_{2,0}$	$a_{1,0}^4$	—
$Q_{3,1}$	$3a_{1,0}{a}_{2,1}+a_{0,1}{a}_{3,0}+3a_{1,1}{a}_{2,0}$	$3a_{1,0}^2{a}_{1,1}+3a_{1,0}{a}_{0,1}{a}_{2,0}$	$a_{1,0}^3{a}_{0,1}$	—
$Q_{2,2}$	$2a_{1,0}{a}_{1,2}+2a_{0,1}{a}_{2,1}+a_{2,0}{a}_{0,2}+2a_{1,1}^2$	$a_{1,0}^2{a}_{0,2}+4a_{1,0}{a}_{0,1}{a}_{1,1}+a_{0,1}^2{a}_{2,0}$	$a_{1,0}^2{a}_{0,1}^2$	—
$Q_{5,0}$	$10a_{2,0}{a}_{3,0}+5a_{1,0}{a}_{4,0}$	$15a_{1,0}{a}_{2,0}^2+10a_{1,0}^2{a}_{3,0}$	$10a_{1,0}^3{a}_{2,0}$	$a_{1,0}^5$
$Q_{4,1}$	$6a_{2,0}{a}_{2,1}+4a_{1,1}{a}_{3,0}+4a_{1,0}{a}_{3,1}+a_{0,1}{a}_{4,0}$	$12a_{1,0}{a}_{1,1}{a}_{2,0}+3a_{0,1}{a}_{2,0}^2+6a_{1,0}^2{a}_{2,1}+4a_{0,1}{a}_{1,0}{a}_{3,0}$	$4a_{1,0}^3{a}_{1,1}+6a_{0,1}{a}_{1,0}^2{a}_{2,0}$	$a_{0,1}{a}_{1,0}^4$
$Q_{3,2}$	$3a_{1,2}{a}_{2,0}+6a_{1,1}{a}_{2,1}+3a_{1,0}{a}_{2,2}+a_{0,2}{a}_{3,0}+2a_{0,1}{a}_{3,1}$	$6a_{1,0}{a}_{1,1}^2+3a_{1,0}^2{a}_{1,2}+3a_{0,2}{a}_{1,0}{a}_{2,0}+6a_{0,1}{a}_{1,1}{a}_{2,0}+6a_{0,1}{a}_{1,0}{a}_{2,1}+a_{0,1}^2{a}_{3,0}$	$a_{0,2}{a}_{1,0}^3+6a_{0,1}{a}_{1,0}^2{a}_{1,1}+3a_{0,1}^2{a}_{1,0}{a}_{2,0}$	$a_{0,1}^2{a}_{1,0}^3$

Equations

$$\zeta =z-\sum _{l=1}^L\frac{q_l}{(\overline{z-s_l})}=z-\overline{W_1},$$

$$W_k\equiv (-1{)}^{(k-1)}(k-1)!\sum _{l=1}^L\frac{q_l}{{\left(z-s_l\right)}^k}.$$

$$J=\frac{\partial \zeta }{\partial z}\frac{\partial \overline{\zeta }}{\partial \overline{z}}-\frac{\partial \zeta }{\partial \overline{z}}\frac{\partial \overline{\zeta }}{\partial z}=1-{\left|\frac{\partial \zeta }{\partial \overline{z}}\right|}^2=1-{\left|W_2\right|}^2.$$

$${\mu }_0=\frac{1}{1-{\left|W_2\right|}^2},$$

$$z\left(\xi ,\eta \right)=z\left({\xi }_0+\delta \xi ,{\eta }_0+\delta \eta \right)\simeq z_0+\sum _{p\ge 1}\frac{1}{p!}{\left[\frac{\partial z}{\partial \xi }\delta \xi +\frac{\partial z}{\partial \eta }\delta \eta \right]}^{[p]},$$

$$a_{p-n,n}=\frac{{\partial }^{p}z}{\partial {\xi }^{p-n}\partial {\eta }^n},$$

$$\frac{\text{d}{W}_k}{\text{d}z}=W_{k+1}.$$

$$\frac{\partial \zeta }{\partial \alpha }=\frac{\partial z}{\partial \alpha }-\overline{W_2}\frac{\partial \overline{z}}{\partial \alpha },$$

$$\frac{\partial z}{\partial \alpha }\left(1-W_2\overline{W_2}\right)=\frac{\partial \zeta }{\partial \alpha }+\overline{W_2}\frac{\partial \overline{\zeta }}{\partial \alpha }.$$

$$\frac{\partial \zeta }{\partial \xi }=1,\frac{\partial \zeta }{\partial \eta }=i,\forall p>1,\frac{{\partial }^p\zeta }{\partial {\xi }^{p-n}\partial {\eta }^n}=0,$$

$$a_{1,0}=$$

$$a_{0,1}=$$

$$\frac{{\partial }^2\zeta }{\partial \alpha \partial \beta }=\frac{{\partial }^{2}z}{\partial \alpha \partial \beta }-\frac{\partial \overline{W_2}}{\partial \beta }\frac{\partial \overline{z}}{\partial \alpha }-\overline{W_2}\frac{{\partial }^2\overline{z}}{\partial \alpha \partial \beta }=0.$$

$$\frac{\partial \overline{W_2}}{\partial \beta }=\frac{d\overline{W_2}}{d\overline{z}}\frac{\partial \overline{z}}{\partial \beta }=\overline{W_3}\frac{\partial \overline{z}}{\partial \beta },$$

$$a_{2,0}=$$

$$a_{1,1}=$$

$$a_{0,2}=$$

$$a_{p-n,n}={\mu }_0\left(\overline{Q_{p-n,n}}+\overline{W_2}{Q}_{p-n,n}\right),$$

$$Q_{p-n,n}=\sum _{k=3}^{p+1}{c}_{n-p,n}^k{W}_k,$$

$$R_p^k=\frac{\partial R_{p-1}^k}{\partial {\alpha }_p}+\left(R_{p-1}^{k-1}+\delta (k,3)\frac{{\partial }^{p-1}z}{\partial {\alpha }_{p-1}\dots \partial {\alpha }_1}\right)\frac{\partial z}{\partial {\alpha }_p},$$

$$Z(\rho,\theta)\equiv z\big{(}\xi_{0}+\rho\cos\theta,\eta_{0}+\rho\sin\theta\big{)}-z_{0}\,.$$

$$S=\text{Im}\left[{\oint }_{\mathcal{C}}\frac{\overline{Z}\text{d}Z}{2}\right]=\frac{1}{2}\text{Im}\left[{\int }_0^{2\pi }\overline{Z}\frac{\partial Z}{\partial \theta }\text{d}\theta \right],$$

$$S(\rho )=\sum _{p\ge 1}\frac{\pi {\rho }^{p+1}}{(p-1)!}{\mu }_{p-1}$$

$${\mu }_0=$$

$${\mu }_2=$$

$$\displaystyle\displaystyle+\;a_{0,1}\;(\overline{a_{1,2}}+\overline{a_{3,0}})\Big{]}\,,$$

$${\mu }_4=$$

$$+4a_{1,1}\overline{a_{4,0}}+6(a_{0,3}\overline{a_{1,2}}+a_{1,2}\overline{a_{2,1}}+a_{2,1}\overline{a_{3,0}})$$

$$+2a_{0,3}\overline{a_{3,0}}+4\overline{a_{1,1}}{a}_{0,4}+4a_{0,2}(\overline{a_{1,3}}+\overline{a_{3,1}})$$

$$\displaystyle\displaystyle+a_{0,1}\;(\overline{a_{1,4}}+2\;\overline{a_{3,2}}+\overline{a_{5,0}})\Big{]}\,.$$

$$I(r)=1-\Gamma \left(1-\frac{3}{2}\sqrt{1-\frac{r^2}{{\rho }^2}}\right),0\le r\le \rho ,$$

$${\iint }_{\mathcal{S}}I(r)\text{d}s={\int }_0^{2\pi }\text{d}\varphi {\int }_0^{\rho }I(r)r\text{d}r=\pi {\rho }^2.$$

$$\mu \equiv \frac{{\int }_0^{\rho }I(r)\text{d}S}{{\iint }_{\mathcal{S}}I(r)\text{d}s}=\frac{1}{{\rho }^2}\sum _{p\ge 1}\frac{(p+1){\mu }_{p-1}}{(p-1)!}{\int }_0^{\rho }I(r)r^p\text{d}r,$$

$$\mu ={\mu }_0+\frac{{\mu }_2}{2!}\left(1-\frac{1}{5}\Gamma \right){\rho }^2+\frac{{\mu }_4}{4!}\left(1-\frac{11}{35}\Gamma \right){\rho }^4+\mathcal{O}\left({\rho }^6\right).$$

$$A=\sum\limits_{j}\big{|}\mu(z_{j})\big{|}=\sum_{j}\epsilon_{j}\,\mu(z_{j})\,,$$

Full text

††pubyear: 2015††pagerange: Fast computation of quadrupole and hexadecapole approximations in microlensing with a single point-source evaluation–Fast computation of quadrupole and hexadecapole approximations in microlensing with a single point-source evaluation

1 Introduction

Since the visionary work of Mao & Paczynski (1991), Galactic gravitational microlensing has led to the discovery of dozens of exoplanets and brown dwarfs,111http://exoplanet.eu/catalog/ and revealed an unexpected population of cold, low-mass exoplanets located beyond the snow-line of their stars. Statistical studies have more recently settled that exoplanets in the Milky Way are the rule rather than the exception (Cassan et al., 2012), thereby opening exceptional prospects to discover exoplanets in a variety of systems and configurations. A recent highlight of exoplanets microlensing search is the characterization of the mass and distance from Earth of planetary microlenses through space parallax measurements: such observations are performed simultaneously with ground-based observatories and from space, using Spitzer (e.g. Udalski et al., 2015; Street et al., 2016) and Kepler/K2 (campaign C9, 2016 April 7 through July 1, Henderson et al., 2016), with a strong involvement of the international microlensing community.

The recent upgrades of ground-based telescopes, including robotic ones, have dramatically increased the amount of photometric data that need to be processed, with thousands of data points for which the models have to be computed. In fact, modelling is currently the most difficult task in microlensing and in most cases the bottleneck of detections delivery. Hence, the improvement of both the strategy of the exploration of the parameter space and the efficiency of the computations are of prime interest, in their mathematical and numerical aspects.

Improved strategies to search the parameter space first include the exploitation of features in the light curves to limit the region to be explored. Albrow et al. (1999) have introduced a model to fit individual caustic crossings independently from the whole light curve. This strategy has been extended to binary-lens caustic-crossing events through the definition of a specific parametrization (dates of caustic entry and exit and corresponding positions of the source centre on the caustics) that allows to limit the search to light curves producing caustic magnification peaks at the dates seen in the data (Cassan, 2008; Cassan et al., 2010). To avoid unnecessary calculations of light curves, Penny (2014) has developed in a similar manner the concept of caustic regions of influence that are defined as empirical analytic expressions limiting the parameter space to regions where most low-mass planetary signals lie. A better understanding of the link between the caustic topography and the resulting light curves is a key ingredient to limit the region in the parameter space to explore. A detailed study of this aspect has been conducted by Liebig et al. (2015) in the binary lens case. In that work, the authors established a classification of all possible peaks in the light curves into four types only, and arranged the corresponding possible light-curve morphologies into 73 categories. Hence, inspecting the characteristics of the observed light curves (number and shape of peaks for example) naturally provides clever initial guesses for the subsequent minimization algorithms (in particular Bayesian algorithms, e.g. Kains et al., 2012).

A second direction of improvement is to design more efficient mathematical methods and numerical codes to perform the calculations of the probed microlensing models. Skowron & Gould (2012) have proposed a new algorithm to solve complex polynomial equations, which can solve the lens equation faster than classical roots finder codes. Nevertheless when finite-source calculations are required, the computation time significantly increases. This has triggered the development of a number of methods or approximations to overcome the problem. While magnification maps obtained by inverse ray-shooting are well suited for a large number of lenses (Wambsganss, 1997), they are in general much too slow to be computed in real time, even for triple lenses. Pre-computed magnification maps are however useful for statistical studies where a large number of simulated light curves needs to be computed (e.g. Kubas et al., 2008; Cassan et al., 2012). A refined image-centred ray-shooting algorithm has been developed by Bennett (2010) to more specifically address the calculation of high-magnification planetary models, in which the source images are highly elongated. Image contouring provides an interesting alternative to ray-shooting (Gould & Gaucherel, 1997; Dong et al., 2006; Dominik, 2007), because it is far less demanding in computation time. Bozza (2010) has significantly improved the contouring method by including an error control algorithm which optimizes the sampling of the contour of the images.

When finite-source effects are noticeable but are weak enough (e.g. caustic or cusp passages without caustic crossing), multipole approximations of the finite-source magnification (Pejcha & Heyrovský, 2009) have proven to be of great help because the computation time is several order of magnitudes below that of the exact finite-source magnification. A simple implementation of the quadrupole and hexadecapole approximations has been proposed by Gould (2008), using respectively 9 and 13 (point source) resolutions of the polynomial lens equation to evaluate numerically the corresponding coefficients of the expansion.

In this work, I present an improved implementation of the quadrupole and hexadecapole approximations, based on the construction of the image contours through a Taylor expansion around the individual images of the source centre, and which makes use of a single resolution of the lens equation. In section 2, I present the method that yields the multipole coefficients of the expansion, and in section 3, I discuss its implementation and numerical efficiency.

2 Multipole expansion

The main steps of the method are as follows. For a given position of the source centre at affix $\zeta_{0}=\xi_{0}+i\eta_{0}\in\mathbb{C}$ (Witt, 1990) in the source plane, I first expand the image position $z\in\mathbb{C}$ around the exact position of the image $z_{0}$ of the source in the lens plane. I then use this expansion to perform the integration of the area of the image through the Green-Riemann formula, from which I obtain the magnification as a series of powers of $\rho^{2}$.

While the method itself seems fairly clear, in practice it is not straightforward to obtain simple expressions of the coefficients of the expansion, which should ideally be easy to calculate and fast to compute numerically, as it is of main interest in this work. I find that by using a combination of a Taylor expansion with respect to the two coordinates in the source plane $(\xi,\eta)$ and using properties of complex numbers provide an elegant and powerful way to express the expansion, as I show below.

Let us consider a microlensing system composed of $L$ components with mass ratio $q_{l}=M_{l}/M$ with respect to the total lens mass $M$, and located at positions $s_{l}\in\mathbb{C}$ in the lens plane. The lens equation then reads222A shear can be added in the equation, with minor changes in the expansion as explained in footnote 4.

[TABLE]

where I have introduced the $W_{k}$ factors, $k\geq 1$, as

[TABLE]

For a point-source $\zeta$, the lens equation provides several images $j$ located at $z_{j}$ (3 or 5 for binary lenses, $n+1$ to $5n-5$ for $n$ lens components, Rhie, 2003), but for clarity I have omitted any explicit reference to $j$ in $z$, $W_{k}$ and other quantities introduced later. The signed point-source magnification $\mu_{0}$ of an image is the inverse of the determinant of the Jacobi matrix of transformation $(z,\overline{z})\mapsto(\zeta,\overline{\zeta})$ (Witt, 1990; Daněk & Heyrovský, 2015), given by $J\equiv\det\;\partial(\zeta,\overline{\zeta})/\partial(z,\overline{z})$, or

[TABLE]

Hence the point-source magnification reads

[TABLE]

the sign of which (in fact, the sign of $J$) gives the parity of the image. The critical curves correspond to infinite values of $\mu_{0}$, or $W_{2}=-e^{-i\phi}$, with $\phi$ a phase parameter ranging from [math] to $2\pi$ (Witt, 1990, see Appendix A for two interpretations of $\phi$).

I proceed now with the first step, the expansion of $z$ around $z_{0}$, the exact image of the source centre $\zeta_{0}$. At order $p\geq 1$, the Taylor expansion of $z$ as a function of coordinates $(\xi,\eta)\in\mathbb{R}^{2}$ in the source plane is

[TABLE]

in which $[p]$ refers to the symbolic binomial expansion of the terms inside the brackets333For example, $[p]=2$ gives $\displaystyle\frac{\partial^{2}z}{\partial\xi^{2}}\;\delta\xi^{2}+2\,\frac{\partial^{2}z}{\partial\xi\partial\eta}\;\delta\xi\delta\eta+\frac{\partial^{2}z}{\partial\eta^{2}}\;\delta\eta^{2}$., and where the derivatives are evaluated at $(\xi_{0},\eta_{0})$. In the following, I will use a more compact notation of the derivatives,

[TABLE]

for $0\leq n\leq p$. Until this point, I have used the linearity of the complex notation as a convenient way to write the Taylor expansion of $z=x+iy$. A naive approach would be to differentiate $x$ and $y$ with respect to $\xi$ and $\eta$, but this requires to separate real and imaginary parts of the lens equation (1), which results in cumbersome, numerically time-consuming expressions of the derivatives, and furthermore depends on the detailed form of the adopted lens equation. In the following, I will therefore use the complex formalism and exploit the property that $W_{k}$ is a single function of $z$, so that444The method can be extended to any expression $W_{k}$ depending on $z$ only, as for example adding a shear $\gamma$. Then, $\zeta=z+\gamma\overline{z}-\overline{W_{1}}$, $W_{1}^{\prime}=-\gamma z+W_{1}$, $W_{2}^{\prime}=-\gamma+W_{2}$ and for $k\geq 3$, $W_{k}^{\prime}=W_{k}$.

[TABLE]

I use this property to compute the derivatives (6), starting with $p=1$ for which I provide the explicit expressions, and use $p=2$ to explain the general method for any $p$.

From the lens equation, one has

[TABLE]

where the differentiation is made with respect to $\alpha\in\{\xi,\eta\}$. To eliminating derivatives of $\overline{z}$, I conjugate (8), isolate $\partial\overline{z}/\partial\alpha$ and introduce it back into (8), which leads to

[TABLE]

Considering that $1-W_{2}\overline{W_{2}}=1/\mu_{0}$ and that

[TABLE]

I obtain the two derivatives $a_{1,0}$ and $a_{0,1}$ ($p=1$) by successively choosing $\alpha=\xi$ and $\eta$,

[TABLE]

For $p=2$, I derive (8) a second time with respect to variables $(\alpha,\beta)\in\{\xi,\eta\}\times\{\xi,\eta\}$, which leads to three different combinations $a_{2-n,n}$ ($n=0$, 1 and 2) satifying

[TABLE]

Using the rule (7), the derivative of $\overline{W_{2}}$ can be expressed as

[TABLE]

so that (12) involves derivatives of $z$ and $\overline{z}$, and constants depending on $z_{0}$ only. With the same procedure as for $p=1$ (conjugating and replacing), I get

[TABLE]

More generally for $p\geq 1$ and $0\leq n\leq p$, it appears that $a_{p-n,n}$ can be expressed as

[TABLE]

where the $Q_{p-n,n}$ coefficients can be iteratively calculated using the prescriptions given above. Factoring the $W_{k}$, I define

[TABLE]

where the coefficients $c_{n-p,n}^{k}$ are given in Table 1 up to order $p=5$ (we shall see later that expanding $z$ up to this order provides the hexadecapole term of the finite-source magnification). For $p\geq 3$, I find that $Q_{p-n,n}$ can be obtained with the following algorithm: let us introduce $p$ variables $(\alpha_{p},\dots,\alpha_{1})\in\{\xi,\eta\}\times\dots\times\{\xi,\eta\}$. For $3\leq k\leq p+1$, the general expression of $c_{n-p,n}^{k}$ is obtained from

[TABLE]

with $\delta(3,k)=1$ only if $k=3$, and in which $n$ of the $\alpha_{i}$ variables are chosen to be $\eta$, and the remaining $p-n$ others $\xi$.

The second step of the method consists of calculating the area of the image using the previous expansion of $z$, with the provisional assumption that the source is uniformly bright. Let $\rho$ be the source radius in Einstein units. For $\rho\ll 1$, the (circular) contour of the source can be parametrized as $\zeta=\zeta_{0}+\rho\cos\theta+i\,\rho\sin\theta$, i.e. by choosing $\delta\xi=\rho\cos\theta$ and $\delta\eta=\rho\sin\theta$ in the expansion of $z$ written in (5). As a matter of fact, $z$ is now a function of $(\rho,\theta)$, and to avoid confusions I introduce

[TABLE]

This expression hence involves powers of $\rho$ as well as powers of $\cos\theta$ and $\sin\theta$. For a given source size $\rho$, the (signed) area $S$ of the image can be performed through the Green-Riemann (or Stokes) formula,

[TABLE]

since for $Z=X+iY$ one has $-Y\textrm{d}X+X\textrm{d}Y=\textrm{Im}[\overline{Z}\,\textrm{d}Z]$. As the expression of $Z$ obtained in (18) is analytical, it is also the case for $\partial Z/\partial\theta$. The integral (19) involves terms like $(\cos^{a}\theta\sin^{b}\theta)$, and the integration is easily handled formally with a software like mathematica. It is interesting to remark that all terms with odd values of $a$ or $b$ cancel out, so after the integration, only terms with even powers of $\rho$ remain. Furthermore, it can be found by inspecting $\overline{Z}$ and $\partial Z/\partial\theta$ that pursuing the expansion to order $p$ adds terms factors of $\rho^{m}$, where $m\geq p+1$. In other words, the expansion of $Z$ is complete in $\rho^{m}$ for $p=m-1$.

Expanding $S$ in the form

[TABLE]

yields, up to order $p=5$, the following non-vanishing terms

[TABLE]

The first term $\mu_{0}$ is indeed the point-source magnification defined in (21), since from (11) one has $\overline{a_{1,0}}\;a_{0,1}=i\mu_{0}^{2}(1-|W_{2}|^{2}-\overline{W_{2}}+W_{2})$ of which imaginary part is $\mu_{0}^{2}(1-\left|W_{2}\right|^{2})=\mu_{0}$. Let us now consider a limb-darkened source with brightness profile

[TABLE]

where for all $\Gamma$, the surface integral of $I(r)$ over the source face $\mathcal{S}$ always equals $\pi\rho^{2}$,

[TABLE]

The uniformly bright source has $\Gamma=0$. The signed magnification is now given by the ratio

[TABLE]

where the enumerator of the first integral is already integrated over the angle, and in which I have changed variable $\textrm{d}S=\frac{\textrm{d}S}{\textrm{d}r}\,\textrm{d}r$ according to (20). The integral can be performed analytically for any value of $p$, and yields

[TABLE]

As expected, the monopole of the finite-source expansion ($p=1$) is the point-source magnification of the source, and there is no dipole ($p=2$). The quadrupole is obtained for $p=3$, and the hexadecapole for $p=5$.

Finally, the total magnification of the source $A(\zeta_{0})$ is the sum of the absolute values of the individual magnification factors $\mu(z_{j})$ of each of the images $j$. If $\epsilon_{j}$ denotes the parity of image $j$ ($\epsilon_{j}=1$ if $\mu>0$, $-1$ otherwise), one has

[TABLE]

so that after factorizing $\rho^{2}$ and $\rho^{4}$ amongst the different images, the total magnification $A$ has the same form555This is the same expansion as Gould (2008) but with a different choice of numerical factors. as (27) with $A_{q}$ instead of $\mu_{q}$, where $A_{q}$ is a combination of $\epsilon_{j}\mu_{q}$.

3 Application

An example of light curves obtained with the various finite-source approximations are displayed in Fig. 1. In this example, the lens is a binary with parameters $s=1.7$ (separation in Einstein units) and $q=0.2$ (lens mass ratio), and the source radius is $\rho=0.01$ in Einstein units (slightly larger than typical source sizes to better see the differences). The source crosses two caustics (entry at $t\sim-4$, exit at $t\sim 8$) and later approaches a cusp ($t\sim 25$). The exact finite-source magnification is displayed as the bold, dark grey curve, while the blue, red and green curves are respectively the monopole, quadrupole and hexadecapole approximations. The middle panel shows a zoom on the residuals (approximated minus exact magnifications) in linear scale, while the lower panel shows the absolute value of the residuals in logarithmic scale. The latter panel can be used to compare the precision of the approximations to the precision of the photometry expected from future space-based missions such as Wide-Field Infrared Survey Telescope, which is expected to be of order of a mmag.

It appears that the quadrupole expansion already provides a much better approximation than the monopole (point-source magnification), in particular, near cusp approaches. The hexadecapole appears essentially as an approximation that allows the source to approach the caustics slightly closer than the quadrupole before the approximation breaks down (see also the discussion in Gould, 2008).

Since one of the drivers of this work was to improve the numerical efficiency in calculating the quadrupole and hexadecapole approximations, I have tested a non-fully optimized python routine to estimate the potential gain for a binary lens. I find that the quadrupole (respectively hexadecapole) is about $\times\,2$ (respectively $\times\,5$) slower than point source. The quadrupole and hexadecapole implementations presented here are respectively $\times\,6$ and $\times\,4$ faster than the implementation of Gould (2008). It is likely that going to higher orders in the finite-source approximation will not help, though, not only because the gain in precision will be limited, but also because the additional calculations (such as $a_{p-n,n}$) grow substantially with increasing order $p$. It is also clear that the gain in time will increase with the number of lens components since solving the lens equation will take more time, while the additional calculations do not, as they depend only on the set of input complex numbers $W_{k}$.

Additionally, a further advantage of using the exact multipole expansion rather than a numerical estimation resides in the fact that the discontinuities of the magnification happen at the same positions in the light curve than with point-source, so there are no spurious wiggles of the magnification close to the caustics.

An open source code with a first implementation of the equations presented here is available for download in my GitHub repository666https://github.com/ArnaudCassan/microlensing/. The file multipoles.py includes two functions that return the quadrupole and hexadecapole approximations of the finite-source magnification: quadrupole(Wk,rho,Gamma) and hexadecapole(Wk,rho,Gamma), whose arguments are $W_{k}$ (a 2D complex array in which one of the dimension refers to the individual images), $\rho$ the source radius in Einstein units and $\Gamma$ the source linear limb-darkening coefficient. In practice (as shown in the example() function provided), for each position of the source, one needs to compute the images of the source centre, discard the virtual ones, and evaluate the corresponding $W_{k}$ (up to $k=4$ for the quadrupole and $k=6$ for the hexadecapole) that are inputs of the two functions. It is likely that these functions will gain in speed from a re-writing in cython or c++ with a more efficient use of complex number calculations. For information, further orders of $a_{p-n,n}$ for $p\geq 6$ can be displayed using function Q(p) in file Rkp.py, as an implementation of (17).

4 Conclusion

I have presented a method that allows to compute the quadrupole and hexadecapole approximations with more efficiency than previously available codes (respectively about six times and four times faster). It appears that the quadrupole approximation already provides a much more precise approximation than point-source, for a computing time only about two times slower. This advocates for using the quadrupole in place of point-source in most part of the light curve, except at baseline. The hexadecapole seems well suited to make the link between exact finite-source and quadrupole in limited regions of the light curves close to sharp magnification peaks. Open source codes of the algorithms presented here are available for download, and I welcome numerical optimization updates.

Acknowledgements

AC acknowledges financial support from Universit Pierre et Marie Curie under grants mergence@Sorbonne Universit s 2016 and mergence-UPMC 2012. I would like to thank the anonymous referee for useful remarks on the text.

Appendix A Two interpretations of Witt’s $\phi$

In his original article, Witt (1990) introduced a phase parameter $\phi\in[0,2\pi]$ to compute the critical curves777The two $-$ signs reflect the different conventions between the formalisms.,

[TABLE]

In this appendix I propose two interpretations of $\phi$.

The first one is geometrical. Starting from the lens equation (1), one can write

[TABLE]

which, after conjugating the whole expression and replacing $d\overline{z}/d\phi$ in (30), yields

[TABLE]

Since $d\zeta/d\phi$ defines in the source plane a vector $T$ tangent to the caustic curve parametrized by $\phi$, the former equation means that $\phi$ is the oriented angle between the symmetrical vector of $T$ with respect to the horizontal axis and $T$. In other words, $\phi/2$ modulo $\pi$ is the geometric angle between the horizontal axis and the tangent to the caustic.

A second interpretation of $\phi$ is related to curves in the lens plane obtained with $\phi={\rm cst}$. From the definition of the Jacobian $J$ in (3) and since $J<1$, one has

[TABLE]

from which we write

[TABLE]

Solving this equation for a given $\phi\in[0,2\pi[$ and $J\in]-\infty,1[$ leads to four possible solutions $z$ in the lens plane. Varying $\phi$ for a given $J$ draws iso-magnification curves (as four distinct branches that connect), that I will refer to as $J$-curves. Let us study the curves obtained for $\phi$ constant. Differentiating both sides of equation (33) with respect to $\phi$ and $J$ yields

[TABLE]

and

[TABLE]

Since $\partial z/\partial J$ defines in the lens plane a vector $T_{J}$ tangent to the $J$-curves, $\partial z/\partial\phi$ hence defines a vector $T_{\phi}$ perpendicular to those, as long as $T_{\phi}\neq 0$ (or equivalently $T_{J}\neq 0$). It means that $\phi={\rm cst}$ lines can be interpreted as field lines (perpendicular to the $J$-curves) parametrized by $J$ and crossing the critical lines at $J=0$. I will call them $\phi$-lines.

An example is given in Fig. 2 for a binary lens with separation $s=1$ and mass ratio $q=0.25$. The $\phi$-lines are the green, radial lines crossing the (thin black) $J$-curves orthogonally (except at saddle points, marked as black $+$, see below). The thick black line is the critical curve. $\phi$-lines can be of finite length or semi-infinite. For $J\rightarrow-\infty$, the only possibility is that at least one of the $z\rightarrow s_{l}$, so that all $\phi$-lines start at one of the lens component positions. For $J\rightarrow 1$, one has $W_{2}\rightarrow 0$, which results in two possibilities for $z$: either $z\rightarrow\infty$ and the line goes to infinity, or $z$ converges to a finite value, which geometrically necessarily correspond to an extrema of mapping $J(z)$. Expanding (33) and setting $\sqrt{1-J}\rightarrow 0$, one finds

[TABLE]

Therefore for $L$ lens components, there are $2(L-1)$ extrema of $J(z)$. In the binary lens case, with the convention that the origin of the coordinate system is at the position of the more massive body $z_{1}=0$ and that the less massive body is located at $z_{2}=-s$, with $q=\mu_{2}/\mu_{1}<1$ their relative mass ratio, there are two extrema at

[TABLE]

which are marked as off-axis $\times$ in Fig. 2.

Finally, the boundary values of $\phi$ which separate regions of different kinds of $\phi$-lines (finite or semi-infinite, starting at one or another lens position) correspond to lines that necessarily cross the $J$-curves at saddle points as stated before, where$W_{3}=0$ (Daněk & Heyrovský, 2015) i.e.

[TABLE]

The boundary $\phi$-lines may start and/or end at one of the lens position, at a maxima of $J(z)$ or at infinity, but they all necessarily pass through a saddle point. Once the saddle points are found through (36), the corresponding values of $\phi$ are obtained from (29) and the corresponding boundary $\phi$-lines can be drawn. In the binary lens case, it appears that there exists six such boundary $\phi$-lines, displayed in Fig. 2 as thick magenta lines.