Measuring Limb Darkening of Stars in high magnification Microlensing Events by the Finite Element Method

TL;DR

This paper introduces a finite element method to measure stellar limb darkening from high-magnification microlensing events, allowing model-independent analysis and application to real observational data.

Contribution

It develops and applies a FEM-based inversion technique to determine stellar limb darkening profiles without assuming specific brightness models.

Findings

FEM accurately recovers limb darkening profiles from simulated data.

Application to real events yields profiles consistent with linear and square-root models.

Method demonstrates stability and model-independence in analyzing microlensing light curves.

Abstract

The finite-size effect in gravitational microlensing provides a possibility to measure the limb darkening of distant stars. We use the Finite Element Method (FEM) as an inversion tool for discretization and inversion of the magnification-limb darkening integral equation. This method makes no explicit assumption about the shape of the brightness profile more than the flatness of the profile near the center of the stellar disk. From the simulation, we investigate the accuracy and stability of this method and we use regularization techniques to stabilize it. Finally, we apply this method to the single lens, high magnification transit events of OGLE-2004-BLG-254 (SAAO-I), MOA-2007-BLG-233/OGLE-2007-BLG-302 (OGLE-I, MOA-R), MOA-2010-BLG-436 (MOA-R), MOA-2011-BLG-93 (Canopus-V), MOA-2011-BLG-300/OGLE-2011-BLG-0990 (Pico-I) and MOA-2011-BLG-325/OGLE-2011-BLG-1101 (LT-I) in which light curves have been observed with a high cadence near the peak \citep{Ch}. The recovered intensity profile of stars from our analysis for five light curves are consistent with the linear limb darkening and two events with the square-root profiles. The advantage of FEM is to extract limb darkening of stars without any assumption about the limb darkening model.

Figures

Tables

Table 1. Table 1 : First column represents the name of microlensing event and the abbreviation name of the observatory that data are taken for our analysis. Here we adapt a single filter light curve. The next four columns are the primary parameters we used in our analysis obtained by (Choi et al., 2012 ) . Next two columns are the magnitude of baseline and the blending parameter. The other columns are the average of uncertainties in derivation of ℓ ​ ( σ ℓ ) ℓ subscript 𝜎 ℓ \ell\;(\sigma_{\ell}) and number of data with ℓ < 1 ℓ 1 \ell<1 .

Event	$t_0$(HJD-2450000)	$u_0$	$t_E$ (days)	${\rho }_{\ast }$	$m_{BL}$	b	average of ${\sigma }_{\ell }$	No. of data
($Observator{y}_{passband}$ )								()
OGLE-2004-BLG-254	$3166.823$	$0.0111$	$12.84$	$0.0418$	16.33	1.039	0.01	56
($SAA{O}_I$)	$\pm 0.001$	$\pm 0.0004$	$\pm 0.09$	$\pm 0.0004$	$\pm 0.014$	$\pm 0.032$
MOA-2007-BLG-233/	$4289.269$	$0.0060$	$15.90$	$0.0364$	16.31	1.021	0.004	9
OGLE-2007-BLG-302	$\pm 0.001$	$\pm 0.0002$	$\pm 0.05$	$\pm 0.0001$	$\pm 0.013$	$\pm 0.033$
($OGL{E}_I$)
MOA-2007-BLG-233/	$4289.269$	$0.0060$	$15.90$	$0.0364$	16.31	1.021	0.004	143
OGLE-2007-BLG-302	$\pm 0.001$	$\pm 0.0002$	$\pm 0.05$	$\pm 0.0001$	$\pm 0.013$	$\pm 0.033$
($MO{A}_R$)
MOA-2010-BLG-436	$5395.791$	$0.0002$	$12.78$	$0.0041$	16.96	0.026	0.092	16
($MO{A}_R$)	$\pm 0.001$	$\pm 0.0002$	$\pm 1.08$	$\pm 0.0003$	$\pm 0.038$	$\pm 0.003$
MOA-2011-BLG-093	$5678.555$	$0.0292$	$14.97$	$0.0538$	16.7	1.8	0.005	206
($Canopu{s}_I$)	$\pm 0.001$	$\pm 0.0002$	$\pm 0.05$	$\pm 0.0002$	$\pm 0.01$	$\pm 0.052$
MOA-2011-BLG-300/	$5758.691$	$0.0151$	$6.70$	$0.0199$	18.49	0.983	0.023	229
OGLE-2011-BLG-0990	$\pm 0.001$	$\pm 0.0004$	$\pm 0.07$	$\pm 0.0003$	$\pm 0.05$	$\pm 0.3$
( $Pic{o}_I$)
MOA-2011-BLG-325/	$5823.574$	$0.0485$	$29.06$	$0.0979$	15.18	1.038	0.06	19
OGLE-2011-BLG-1101	$\pm 0.002$	$\pm 0.0005$	$\pm 0.11$	$\pm 0.0006$	$\pm 0.008$	$\pm 0.014$
($L{T}_I$)

Table 2. Table 2 : The first column represents the name of microlensing event and the abbreviation name of the observatory that data are taken for our analysis. Here we adapt a single filter light curve. The second column is the source type and the third column is the best fitted value of LDC from previous studies, the fourth column is the best value of LDC to the regularised FEM solutions that are shown in Figures ( 16 - 22 ), the fifth column is the χ 2 / N d ​ o ​ f superscript 𝜒 2 subscript 𝑁 𝑑 𝑜 𝑓 \chi^{2}/N_{dof} from fitting the LLD profile to the recovered profile. The sixth, seventh columns are the best values of square root LDCs (see equation 36 ) and the eighth column is the corresponding χ 2 / N d ​ o ​ f superscript 𝜒 2 subscript 𝑁 𝑑 𝑜 𝑓 \chi^{2}/N_{dof} .

(1) Event name	(2) Source type	(3) Best fitted $u_{LLD}$	(4) $u_{LLD}$	(5) ${\chi }^2/N_{dof}$	(6) $c$	(7) $d$	(8) ${\chi }^2/N_{dof}$
(Observatory-passband)	( log g, $T_{eff})$	of previous studies		(LLD)			(square-root)
	Choi, et al. 2012
OGLE-2004-BLG-254	KIII	$0.55\pm 0.06$	$0.575$	$29.51/56$	$0.262$	$0.48$	$8.14/55$
		Choi et al. (2012)	$\pm 0.028$		$\pm 0.03$	$\pm 0.04$
(SAAO-I)	(2.0, 4750 K)	${0.45}_{-0.06}^{+0.03}$
		Cassan et al. (2006)
MOA-2007-BLG-233/	GIII	$0.53\pm 0.04$	$0.537$	$16.00/8$	$0.698$	$-0.25$	$7.07/7$
OGLE-2007-BLG-302	(2.5, 5000 K)	Choi et al. (2012)	$\pm 0.03$		$\pm 0.03$	$\pm 0.05$
( OGLE-I )
MOA-2007-BLG-233/	GIII	$0.56\pm 0.02$	$0.56$	$98.38/143$	$0.662$	$-0.25$	$78.95/142$
OGLE-2007-BLG-302	(3.0, 5500 K)	Choi et al. (2012)	$\pm 0.05$		$\pm 0.04$	$\pm 0.06$
(MOA-R )
MOA-2010-BLG-436	…	$0.52\pm 0.1$	$0.51$	$62.38/16$	$0.025$	$0.783$	$24.75/15$
(MOA-R)		Choi et al. (2012)	$\pm 0.08$		$\pm 0.03$	$\pm 0.04$
MOA-2011-BLG-093	GIII	$0.51\pm 0.03$	$0.50$	$1423.62/198$	$-0.046$	$0.876$	$409.76/197$
( Canopus-I)	(2.5, 5000 K)	Choi et al. (2012)	$\pm 0.03$		$\pm 0.03$	$\pm 0.04$
MOA-2011-BLG-300/	…	$0.56\pm 0.04$	$0.39$	$35.93/226$	$0.228$	$0.456$	$1.575/225$
OGLE-2011-BLG-0990		Choi et al. (2012)	$\pm 0.1$		$\pm 0.14$	$\pm 0.23$
(Pico-I)
MOA-2011-BLG-325/	KIII	…	$0.94$	$33.32/19$	$0.160$	$1.217$	$10.35/18$
OGLE-2011-BLG-1101	(2.0, 4250 K)		$\pm 0.07$		$\pm 0.051$	$\pm 0.08$
(LT-I)

Equations

$${\theta }_E=\sqrt{\frac{4GM}{c^2}\frac{D_{ls}}{D_l{D}_S}}.$$

$$A(u)=\frac{u^2+2}{u\sqrt{u^2+4}},u=(u_0^2+\frac{(t-t_0{)}^2}{t_E^2}{)}^{\frac{1}{2}}$$

$$A(\ell )$$

$$\ell =\sqrt{p^2+\frac{(t-t_0{)}^2}{t_{\ast }^2}},$$

$$A(\ell ,r)$$

$$A(\ell ,r,\varphi )$$

$$A(\ell ,r)$$

$$\text{w}here$$

$$\kappa$$

$$K(\kappa )$$

$$\Pi ({\alpha }^2,\kappa )$$

$$A(\ell ,\ell +\delta )$$

$$g(t)={\int }_{x_b}^{x_f}K(x,t)f(x)dx,$$

$$g(t_i)=\sum _{j=1}^M{\int }_{x_j}^{x_{j+1}}K(x,t_i)f(x)dx,i=\mathrm{1..}N$$

$$\bar{x}=2\frac{x-x_j^1}{{\Delta }_j}-1,{\Delta }_j=x_j^{n_{d2}}-x_j^1.$$

$$f_j(\bar{x})=\sum _{j_d=1}^{n_{d2}}{U}_j^{j_d}(\bar{x})f_j^{j_d},j=\mathrm{1..}M,$$

$$f_j(\bar{x})$$

$${\mathbf{U}}_j$$

$$\tilde{f}(x)=\sum _{j=1}^M{H}_j(x)f_j(\bar{x}),$$

$$\tilde{g}(t)$$

$$N$$

$${\mathbf{V}}_i(\bar{t})$$

$$\sum _{k=1}^N{H}_k(t)g_k(\bar{t})=\sum _{j=1}^M{\int }_{x_b}^{x_f}{H}_j(x)K(x,t)f_j(\bar{x})dx.$$

$$H_i(t){\mathbf{V}}_i(\bar{t}).{\mathbf{g}}_i=$$

$${\mathbf{g}}_i={\mathbf{G}}_{nn}^{-1}.\sum _{j=1}^M\frac{{\Delta }_j}{2}{\int }_{-1}^1{\int }_{-1}^{1}K(x,t){\mathbf{V}}_i(\bar{t})\otimes {\mathbf{U}}_j(\bar{x})d\bar{x}\bar{t}.{\mathbf{f}}_j,$$

$${\mathbf{G}}_{nn}={\int }_{-1}^1{\mathbf{V}}_n(\bar{t})\otimes {\mathbf{V}}_n(\bar{t})d\bar{t}.$$

$${\mathbf{g}}_i$$

$${\mathbf{E}}_{ij}$$

$$\mathbf{g}=\mathbf{A}.\mathbf{f},\mathbf{A}=\sum _{j=1}^M\sum _{i=1}^N{\hat{\mathbf{E}}}_{\mathbf{ij}},\mathbf{g}=\sum _{i=1}^N{\hat{\mathbf{g}}}_{\mathbf{i}},$$

$$A(\ell )=\sum _i^N{H}_i{U}_i(\bar{\ell }).A_i,\bar{\ell }=2\frac{\ell -{\ell }_i}{{\ell }_{i+1}-{\ell }_i}-1$$

$$I(r)=\sum _j^M{H}_j{V}_j(\bar{r}).I_j,\bar{r}=2\frac{r-r_i}{r_{i+1}-r_i}-1$$

$${\mathbf{A}}_i$$

$$\text{w}here:$$

$${\mathbf{E}}_{ij}$$

$${\mathbf{G}}_{ii}$$

$$\mathbf{A}=\mathbf{M}.\mathbf{I},\mathbf{M}=\sum _{j=1}^M\sum _{i=1}^N{\hat{\mathbf{E}}}_{\mathbf{ij}},\mathbf{A}=\sum _{j=1}^M{\mathbf{A}}_i,$$

Full text

††pagerange: Measuring Limb Darkening of Stars in high magnification Microlensing Events by the Finite Element Method–References††pubyear: 2016

1 Introduction

The limb of stellar disks is dimmer and redder than their center, this effect is known as Limb Darkening (LD) effect. LD happens because photons coming from the center of the stellar disk originate deeper in the photosphere than photons from the limb where the temperature is lower (Gray, 1992). The result is that the light from the center of the stellar disk is more intense and its temperature is higher. Direct study of stellar LD is possible by interferometric photometry of nearby giant stars in one or several bands (Burns et al., 1997; Perrin et al., 2004; Aufdenberg et al., 2006; Wittkowski et al., 2006; Montargès et al., 2014). The other method is studying special events such as eclipsing binaries (Popper, 1984; Southworth et al., 2005, 2015)  or occulting systems (Richichi & Lisi, 1990).

The other new method is the finite-size effect in the gravitational microlensing events (Albrow et al., 2001; An et al., 2002; Yoo et al., 2004; Choi et al., 2012; Rahvar, 2015) which is the subject of our study and can be used as a probe to scan the intensity profile of distant stars. Gravitational microlensing happens when a massive astronomical object inside the Milky Galaxy intervenes as a background star and bends its light toward the observer. Since the observer, lens and source are moving inside the Galaxy, the angular position of the lens compared to the source changes by time and as the angular separation gets closer, it results in the increase of the magnification of the source star (Paczyński, 1986). The time-scale of magnification scales with the square root of the lens mass and can take from hours to almost one month.

According to the original paper by Einstein (1936), he investigated the gravitational lensing of a background star by another star and he stated that ”there is no great chance of observing this phenomenon”. This phenomenon has also investigated by pioneer astronomers such as F. Link. For more details refer to Valls‐Gabaud (2006). Due to instrumental progress, nowadays thousands of microlensing events are observed towards the center of Galaxy by OGLE111http://ogle.astrouw.edu.pl/ and MOA222http://www.phys.canterbury.ac.nz/moa/ surveys and other follow-up groups as $\mu$-Fun333http://www.astronomy.ohio-state.edu/ microfun/, MindStep444http://www.mindstep-science.org/, Planet555http://www.planet-legacy.org/. The gravitational microlensing has broad astrophysics applications such as investigating dark compact objects so-called MACHOs 666Massive Astrophysical Compact Halo Objects in the Milky Way halo (Paczyński, 1986). The two observational groups of EROS and MACHO after a decade monitoring of Magellanic Clouds for the microlensing events concluded that MACHOs don’t have a significant contribution in the dark matter contribution of the Galactic halo  (Alcock et al., 2000; Afonso et al., 2003).

The other application of microlensing is using this method for detecting extrasolar planets (Gaudi et al., 2008; Gaudi, 2012; Tsapras, 2018) and even detecting signals from Extraterrestrial intelligent life (Rahvar, 2016). Also, the microlensing can be used for studying the stellar spots on the source star by polarimetry (Simmons, Willis & Newsam, 1995; Sajadian, 2015) and time variation of center of light of the source star by astrometry (Walker, 1995; Sajadian & Rahvar, 2015). Studying the structure of  Milky Way through the combination of photometry and astrometry observations with GAIA is another important application of gravitational microlensing (Rahvar & Ghassemi, 2005; Moniez et al., 2017).

In this work, our aim is the application of gravitational microlensing for studying the LD of the source stars during the lensing. In microlensing events with the minimum impact parameter comparable with the size of the source star, the source star cannot be taken as a point-like object. The result of this effect, so-called finite-size effect (Schneider & Weiss, 1986; Schneider & Wagoner, 1987; Witt & Mao, 1994) is the deviation of the light curve from a point-like source around the peak. The other feature of this effect is that when the lens crosses over the source star, the main contribution of light is received from the location of the lens on the source star. This effect turns the microlensing effect to an astronomical scanner that can probe the surface of the source stars. This kind of source scanning also happens in the binary lenses where the lensing system produces caustic lines. The observation of these events with high cadence allows us to probe the detailed structure of the source star such as LD  (Valls-Gabaud, 1995; Witt, 1995; Valls-Gabaud, 1998; Fields et al., 2003; Cassan et al., 2006) and stellar spots (Heyrovský & Sasselov, 2000; Hendry, Bryce & Valls-Gabaud, 2002).

Here, we study the single-lens finite-size effect in the microlensing events to recover the LD of the source star, using the Finite Element Method (FEM) (Zienkiewicz, Taylor & Zhu, 2005). This method is an inversion tool to numerically solve the magnification-LD equation. There are other inverse methods that have been used for recovering limb darkening of source stars such as discretization (Gaudi & Gould, 1999; Bogdanov & Cherepashchuk, 1996), Backus-Gilbert method (Gray & Coleman, 2000) and PCA 777Principal Component Analysis inversion (Heyrovský, 2003). Heyrovský (2003) presents a detailed review about numerical methods for recovering LD.  The magnification-LD equation is a Fredholm integral equation of the first-kind (Wazwaz, 2011) and the result of solving this equation is recovering the limb-darkening profile data from the light curve data around the peak.

In section (2), we briefly introduce gravitational microlensing and finite-size effect. In section (3) we explain the FEM approach and apply it to a generic Fredholm integral equation of the first-kind. In section (4) we apply FEM to the magnification-LD equation by suitable adjustment of stellar disk mesh and adequate numerical integration technique. To examine systematic errors we recover the brightness profile from non-noisy simulated light curves; then to examine the accuracy and stability of our inversion method, we add different ranges of Gaussian noise to flux values and estimate the variation of LD profile induced by the flux noise. To achieve stability and better accuracy we introduce a data selection algorithm and use regularization techniques on FEM solutions. Finally, in this section, we apply our method to the single lens, high magnification data ofOGLE-2004-BLG-254 (SAAO-I888South African Astronomical Observatory, South Africa, I  passband), MOA-2007-BLG-233/OGLE-2007-BLG-302 (OGLE-I999Las CampanasObservatory, Chile, MOA-R101010Mt. John Observatory,NewZealand), MOA-2010-BLG-436 (MOA-R), MOA-2011-BLG-93 (Canopus-V111111Canopus Hill Observatory), MOA-2011-BLG-300/OGLE-2011-BLG-0990 (Pico-I121212Observatorio do Pico dos Dias, Brazil) and MOA-2011-BLG-325/OGLE-2011-BLG-1101 (LT-I131313Liverpool Telescope, La Palma, Spain) events to extract directly the LD profile of the source star. The conclusion is given in section (5).

2 Gravitational Microlensing and Finite Size Effect

When the light ray of a star (source) passes closely enough to another astronomical object (lens) it bends due to the gravitational field of the lens (Einstein, 1936). This effect causes secondary images from the source (Eddington, 1920) or a ring image in the case that we have perfect alignment of the source, lens, and observer (Chwolson, 1924). The angular radius of this ring is called angular Einstein radius,

[TABLE]

If the separation of resultant images is of the order of micro arcsec (microlensing events) the images are not resolvable but the source star will be magnified (Paczyński, 1996). During a microlensing event, the apparent brightness of the source star will rise and finally drops to the baseline. The time-dependent magnification of a point-like source by a single lens is as follows:

[TABLE]

in which $u$ is the angular separation of the lens and source in units of $\theta_{E}$, $u_{0}$ is the minimum impact parameter, $t_{0}$ is the time of maximum magnification and $t_{E}$ is the Einstein timescale, corresponds to the time-scale that source move an angular distance of $\theta_{E}$ .

For the case that the minimum impact parameter (i.e. $u_{0}$) becomes comparable with the angular size of source star (normalized to the Einstein angle), the lens amplifies different parts of the source star with different weights and it makes possible to study the surface brightness profile and size of the source star (Witt, 1995). In this case, the magnification is obtained from the convolution of equation (1) and surface brightness of the source. In the observations, the high magnification events alerted well before the peak so that a network of follow-up telescopes can perform high cadence observation (Alcock et al., 1997; Choi et al., 2012). From the measurement of the light curve around the peak, one can calculate the angular size of the source, $\theta_{\star}$ in units of $\theta_{E}$ (i.e. $\rho_{\ast}=\theta_{\star}/\theta_{E}$). Knowing the type of the source star and the distance of the source star from the observer which is mainly located at the Galactic Bulge, one can measure the Einstein angle of the lens.

The second channel for the finite-size effect observation is during the caustic crossing of the binary lens where the LD also can be measured with this method (Albrow et al., 1999; Zub et al., 2011).In this paper, our aim is to study the single-lens very high magnification events where due to a large number of microlensing events, the number of these types of events has been increased in recent years. Such events have been reported by Choi et al. (2012) where standard linear Limb Darkening Coefficient (LDC) of source stars have been obtained from the fitting model with the light curve.

Throughout this paper we adapt the normalized impact parameter to the angular size of star by $\ell={u}/{\rho_{\ast}}$, normalized minimum impact parameter to the size of star by $\;p={u_{0}}/{\rho_{\ast}}$ and the transit time scale of lens crossing over the star by $\;t_{\ast}=\rho_{\ast}t_{E}$. Here, $\ell=1$ corresponds to when the lens enters or leaves the source disk. Now we can calculate the magnification of a source with circular symmetric LD profile (i.e. $I(r)$) from a single lens  as follows:

[TABLE]

where $r$ is within the range of $[0,1]$ and $\ell$ depends on the location of lens with respect to the centre of source star and is given by

[TABLE]

$F_{0}=2\pi\int_{0}^{1}\,I(r)\,r\mathrm{d}r$ and $A(\ell,r)$ is the angle-integrated amplification:

[TABLE]

where $\epsilon={1}/{\rho_{\ast}}$. Equation (4) can be written in terms of elliptic integrals as follows (Gray, 2000) :

[TABLE]

$K$ and $\Pi$ are the first and the third-kind elliptic integrals, as follows:

[TABLE]

$A(\ell,r)$ has a logarithmic divergence in $r=\ell$ as shown in Figure (1), meaning thereby that at $r\to\ell$, parts of source star close to this point are amplified much more than the other parts. This property turns microlensing to a natural surface scanner.

The angle-integrated amplification can also be approximated near this divergency as follows (Heyrovský, 2003):

[TABLE]

3 Finite Element Method in Fredholm Integral Equations of the first-kind

The equation (2) as the magnification-LD integral is a Fredholm integral equation of the first-kind. We deal with such equations in several other situations in the astronomy (Craig & Brown, 1986). As the history of application of this method in astronomy, that was used in the galactic dynamics for modeling perturbed stellar systems (Jalali, 2010) and  constructing smooth distribution functions of stellar systems (Jalali & Tremaine, 2011). Let us take the integral as follows

[TABLE]

in which $f(x)$ is unknown and $K(x,t)$ and $g(t)$ are known. Let us first introduce the Product Integration Method (PIM) which is simpler than FEM but similar in some aspects, in PIM one takes $N$ data points : $(t_{i},g(t_{i}))$ and divides the $x$-domain into $M$ parts ($M\geqslant N$)and  writes discrete version of equation(9) as follows:

[TABLE]

Then by choosing a $f(x)$ to be piecewise constant or piecewise linear over each part, a simple algebraic formula can be derived for each data point.Then one can gather all above $N$ equations into a $N\times M$ set of algebraic equations.If one takes $M>N$ additional constraints such as monotonically and positiveness can be met. See Craig & Brown (1986) for more details.The main difference between FEM and PIM is that in FEM we assure the continuity of the solution and we can estimate data as a continuous piecewise polynomial as well. Below we explain this in more detail.

To solve equation (2) by use of FEM, we approximate $g(t)$ as a continuous piecewise $n$-degree polynomial and we find $f(x)$ as a continuous piecewise $m$-degree polynomial. To do this, first we divide $x$-domain into $M$ elements, each element contains $n_{d2}=m+1$ nodes on its boundaries or in its interior. Then we give two numbers to each node, one indicates the location of the node in the element (local numbering) and the other one is the location in the whole domain (global numbering), the total number of nodes $N_{2}$ is $N_{2}=M(n_{d2}-1)+1$.  See  Figure (2) as an example of a one-dimensional FEM-mesh and its numbering method where the size of elements essentially are not equal.

We can use certain local coordinates within each element, so that ranges of all elements become the same: $x_{i}^{1}\leq x\leq x_{i}^{n_{d2}}\to\,-1\leq\bar{x}\leq 1$. These local coordinates are defined as bellow:

[TABLE]

Within each element we write $f(x)$ as a linear combination of $n_{d2}$  functions which are $m$-degree polynomials. We choose these polynomials such that coefficients of the linear combination become nodal values ($f_{j}^{j_{d}}$, in which $j$ is the element number and $j_{d}$ is the local node number.). Therefore each polynomial should take the value of 1 in one node and value of 0 in other nodes, these basis functions are called shape functions in the FEM context (See Figure (3)).

Hence we can write the approximation of $f(x)$ within $j$th element as follows:

[TABLE]

where $U_{j}^{j_{d}}$s as seen in Figure (4) are the shape functions. We can write equation (12) in a compact form as a vector inner product (throughout this article we denote vectors and matrices by $\mathbf{bold}$ characters, dot product by ., transpose by T and $\otimes$ stands for the dyadic product.):

[TABLE]

Now we can write piecewise approximation of $f(x)$ in whole $x$-domain by summing over $f_{j}(\bar{x})$ of all elements:

[TABLE]

where $H_{j}(x)$ is the top hat function where it is zero every where and one within $j$th element.

Let us go back to equation (9), if $g(t)$ is known in $N_{1}$ points in $t$-domain, we can write its approximation as a piecewise $n$-degree polynomial by the same procedure:

[TABLE]

By substituting (16) and (15) in (9), in the FEM formalism this equation can be written as:

[TABLE]

If we multiply both sides of the above equation by $H_{i}(t)$ we get an equation for the $i$th element:

[TABLE]

The above equation shows the relation between $i$th element of $t$-domain and all elements of $x$-domain, hence we can proceed to derive a relation between nodal values of both sides. To do so we use Galerkin projection technique (Zienkiewicz, Taylor & Zhu, 2005): we operate (19) with $dt\mathbf{U}_{i}\otimes$.Then we integrate the result over the $t$-domain and transform to the local coordinate systems of $t$ and $x$-domains, we get:

[TABLE]

We can write the right hand side of (3) in terms of a summation over $M$, $n_{d1}\times n_{d2}$ matrices ($\mathbf{E}_{ij}$):

[TABLE]

Next step is assembling all $N$ set of equations (22) into one large set of algebraic equations, considering equality of nodal values in common nodes of neighbour elements ($g_{i}^{1}=g_{i-1}^{n_{d1}},\;f_{j}^{1}=f_{j-1}^{n_{d2}}$), see Figure (5)).

To do this we use global numbering of nodes, we rewrite all $\mathbf{E}_{ij}$ matrices which have the $n_{d1}\times n_{d2}$ dimension in form of $N_{1}\times N_{2}$ matrices (i.e. $\mathbf{\hat{E}}_{ij}$) and $\mathbf{g}_{i}$ in form of a $N1\times 1$ vector (i.e. $\mathbf{\hat{g}}_{i}$). In other words, $E_{ij}^{i_{d}j_{d}}\to\hat{E}_{ij}^{IJ}$ and $g_{i}^{i_{d}}\to\hat{g}_{i}^{I}$ where $(i_{d}=1..n_{d1},j_{d}=1..n_{d2},I=1..N_{1},J=1..N_{2})$.

This process is shown in a schematic way in Figure (6). We get $N$ set of algebraic equations, we assemble them all to get a global set of equations which is the approximation of the original integral-equation (9):

[TABLE]

in which $\mathbf{f}$ is the vector of nodal values of $f(x)$ ordered by their global numbering and the $(N_{1}\times N_{2})$ matrix ($\mathbf{A}$) is called Global Stiffness Matrix (GSM) (see Figure (7) for a schematic description of assembling process).

Now if we get a well-posed GSM we can derive nodal values of $f(x)$ by solving linear algebraic equation (23) otherwise one could use regularization techniques described by Craig & Brown (1986). Then using (15), the piecewise approximation of $f(x)$ can be calculated in entire $x$-domain.

4 FINITE ELEMENT METHOD IN MICROLENSING

In this section, we apply FEM in the high magnification microlensing event, with the finite-size effect to recover the LD profile (the integral-equation 2), following the procedure described in the former section. Suppose we have $A(\ell)$ in $N_{1}$ nodes, so we can divide the $\ell$-space into $N$ elements and write the Continuous Piecewise Polynomial (CPP) approximation of $A(\ell)$ by using shape functions that we introduced before as follows:

[TABLE]

We impose the error in each magnification as $A_{i}=A(l_{i})\pm\delta A_{i}$ where $l_{i}={u(t_{i})}/{\rho_{\ast}},\,i=1..N_{1}$.We write $1$-degree CPP approximation of normalised LD profile (i.e. $I(r)/{F_{0}}$) by total $N_{2}$ nodes and $M=N_{2}-1$ elements in $r$-space ($r_{1},r_{2},...,r_{N_{2}}$) :

[TABLE]

From equation (22) using notation for gravitational microlensing, we rewrite these equations as follows:

[TABLE]

$\mathbf{E}_{ij}$s in equation (26) due to logarithmic divergency of $A(l_{i},r)$ at $r=l_{i}$ (equation 2), can not be calculated by simple numeric integration methods . To carry out $\bar{r}$ integrations we use Runge-Kutta adaptive step size method (Press et al., 1992). With this method we calculate the $\bar{r}$ integrations up to precision of $\mathcal{O}(10^{-7})$. We use Gaussian quadrature for $\bar{\ell}$ integration. By assembling $\mathbf{E}_{ij}$s matrices we get the global algebra set of equations which is the approximation of the original integral-equation (2):

[TABLE]

in which $\mathbf{I}$ is the vector of nodal values of $I(r)$ ordered by their global numbering. By solving this set of algebraic equation we can derive $\mathbf{I}$ and $I(r)$ using equation (15). The resultant light curve of such $I(r)$ will go through all data points near peak ($A_{i},i=1..N_{1}$).

To compare FEM with a simpler numerical inversion technique, we apply PIM (Equation 10 ) to the magnification-LD equation (Equation 2) as well. If we choose the result to be piecewise constant we get:

[TABLE]

As in FEM we use Runge-Kutta adaptive step size method to derive $D_{ij}$s up to precision of $\mathcal{O}(10^{-7})$.

4.1 Simulation Details: Classical FEM

In this section, we simulate microlensing light curves with the finite-size effect to study the quality of recovered LD profile ($I_{FEM}$) by FEM. We also compare it with the results from the PIM (i.e. Equation 28). In this simulation, we examine the effect of different parameters i.e. minimum impact parameter, data cadence and error bars in the light curve data on the quality of the recovered LD profile.We simulate the light curves with the parameters of ($p,t_{\ast},t_{0},\rho_{\ast}$) with taking a Linear Limb Darkening (LLD) profile and normalize it to the unit total flux (NLLD):

[TABLE]

Where $\Gamma$ is the linear Limb Darkening Coefficient (LDC) and it depends on surface gravity ($\log g$), effective temperature ($T_{eff}$) and metallicity ($Z$) of a star. The LLD profiles that is more commonly used is in the following form of

[TABLE]

There is a straightforward relation between $\Gamma$ and $u_{LLD}$ that is: $u_{LLD}={3\Gamma}/{(2+\Gamma)}$.Using stellar atmospheric models, one can derive LDCs for stars with different $\log g,T_{eff},Z$ (Claret, 2000).

We use equation (2) to calculate $A_{i}$ in $N_{1}$ different moments of $t_{i}$ starting from $t^{-}_{transit}$ to $t^{+}_{transit}$, the time when the lens enters and leaves the source disk at $t^{\pm}_{transit}=t_{0}\pm t_{\ast}\,\sqrt{1-p^{2}}$. We do not use outside the range of $t^{-}_{transit}<t<t^{+}_{transit}$ in our analysis as the LD-effect on the light curve is negligible.

For simplicity, we choose uniform cadence of $t_{i}$ ($t_{i+1}-t_{i}=\Delta t$) in this simulation. For each $t_{i}$ the projected angular distance between the lens and the centre of source disk is $l_{i}=\sqrt{p^{2}+\frac{(t_{i}-t_{0})^{2}}{t_{\ast}^{2}}}$, where $i=1..N_{1}$; this leads to producing elements in the $\ell-$domain with different sizes, smaller near centre ($\ell=p$) and larger near edge ($\ell=1$). Then we consider an error bar of $\sigma_{i}$ for each point of the light curve and the magnification from the theoretical light curve shifted by a Gaussian distribution with the width of $\sigma_{i}$, where $\sigma_{i}$s results from the uncorrelated magnification error bars of the light curves observed by OGLE and MOA (Choi et al., 2012).

In the next step we discretize the source of the microlensing event by dividing the stellar disk into $N_{2}$ annuli ($r_{j}$ where $j=1..N_{2}$) and choose annuli such that they cover the whole stellar disk. We note that in the FEM, $N_{2}$ might be larger than $N_{1}$ (i.e. $N_{2}\geq N_{1}$). Here, in the application of the FEM method, we correspond a map between each data point in the magnification space and the source space where for each element in this space we have at least one corresponding data. We have tested that having an annuli-element in the source space with no corresponding data in the magnification space results in a large numerical error.

Moreover, we add a constrain from the physics of the LD to the set of algebraic equation. For the intensity of star at the centre where we call it $I(r_{1})$ according to our convention, the radial derivate along $r$-coordinate is zero. This means that $I(r_{2})\simeq I(r_{1})$ where $r_{2}$ is the second nodal value in the $r$-space. We use the convention of $I(r_{i})=I_{i}$ as we introduced in FEM formalism. Using $I_{2}-I_{1}=0$ constrain and equation (27), we have $N_{1}+1$ equations with $N_{2}$ unknowns intensities for each annuli (i.e. $I_{1},I_{2},...,I_{N_{2}}$). If we take $N_{2}=N_{1}+1$, then we get a unique solution with using a linear algebraic equations solver. Here we use LU141414Lower triangle, Upper triangle. -decomposition technique (Press et al., 1992).

Now we start with simulation of light curve to examine the FEM. For the first step, we simulate data points of the light curves without taking into account the error bars (i.e. $A_{i}=A(l_{i}),\,\sigma_{i}=0$). The result for the reconstructed LD profile is limited by the errors of the numerical method due to discretisation and the roundoff errors. We take the following set of parameters for our numerical experiment ($p=0.1,\;t_{*}=0.1458$ days, $t_{0}=0,\;\rho_{*}=0.02,\;u_{LLD}=0.6$) and the cadence of $\Delta t=(t_{transit}-t_{0})/(N_{1}-1),\,N_{1}=15$. Here the theoretical light curve for this event is shown in Figure (8) and the results from reconstructing the intensity of the source are shown in Figure (8). We can see that the residuals are larger near the limb as the intensity-derivative of star near the limb is larger than the central part of the star and in order to improve the results, we need more sampling near the limb area. In order to compare the FEM with that of PIM, we plot the result in Figure (8). The relative errors and residuals from this method are two order of magnitude higher than the case of FEM.

In the next step we study the effect of error bars on the reconstructed intensity from the FEM. We simulate a microlensing event with the parameters given in the first part of this section and the uncorrelated error bars from real MOA and OGLE observations. The average value of the error bars in terms of the magnitude is around $0.005$. Then we use the Monte Carlo simulation and produce $1000$ realization from the same event where each event is different than the other in terms of the measure value of the magnification which is given according to the Gaussian error-bar. Fig.9 represents a light curve and associated error bar for each data point. We take the mean value of 1000 solutions as the nodal values and associate the error bars of $I_{FEM}$ from the variance

[TABLE]

The results is shown in Fig.9. The average profile of the reconstructed profile (i.e. $<I_{FEM}>$) is close to the simulation’s input profile. However, if we take randomly a reconstructed profile from FEM, there are dispersions around some of the points from the input LD profile. However to overcome this issue we can use regularization techniques (Craig & Brown, 1986) where we will describe it in the next section.

We examine the effect of photometric precision of data on recovered intensity profile, By Monte Carlo simulation, we produce $1000$ realisation for a set of events ($p=0.1,\;t_{*}=0.1458$ days, $t_{0}=0,\;\rho_{*}=0.02,\;u_{LLD}=0.6,N_{1}=10$)) with different set of $\sigma_{m}$. For each set of light curve we calculate $\sigma_{I_{FEM}(r_{i})}$ and in order to have an overall dispersion around the FEM result, we calculate the median of $\sigma_{I_{FEM}(r_{i})}$, denoted by median($\sigma_{I}$) for whole of stellar-disk. This parameter provides a relation between the dispersion of the final result of FEM to the photometric accuracy. The other relevant parameter for the quality of the FEM is the median of absolute residuals $|I_{model}(r_{i})-I_{FEM}(r_{i})|$, also denoted by median($R$). The results are shown in Figure (10) where by reducing the photometric errors, the precision of LD data becomes better and converging to the model however there is a limit which might be depend on the sampling rate, data coverage and numerical errors. The same procedure is done for PIM. The result is shown in Figure (10), we see that median($\sigma_{I}$) for FEM is one order of magnitude smaller compare to PIM and we have the same situation for median($R$) comparing FEM with the PIM.

Next we study the effect of impact parameter on the quality of the reconstructed LD profile. We fix data cadence ($N_{1}=10$) and $\sigma_{m,i}\simeq 0.005$ while changing impact parameter. Results are shown in Figure (11) where the reconstructed function of LD is in favor of the small impact parameters.

The other important parameter is the uniformity of the data in the light curve. In our simulations we found that two neighboring data points being closer compared to the average time steps of the data set results in a larger error of the corresponding LD nodal value. To study this effect, we produce data with uniform cadence ($t_{i+1}-t_{i}=\Delta t$) then we put two neighboring points close together ($t_{n}=t_{n-1}+\delta\Delta t,\,0<\delta\leq 1$) while other points remain in the same place. Then for each $\delta=(t_{n}-t_{n-1})/(t_{n+1}-t_{n})$ we produce 1000 light curves and we recover LD nodal values and calculate the variance of $I_{FEM}$ in $r_{n-1}=\ell(t_{n-2})$, $r_{n}=\ell(t_{n-1})$, $r_{n+1}=\ell(t_{n})$, $r_{n+2}=\ell(t_{n+1})$ to compare variances of nodal values close to this defect. Results are shown in Figure (12) where the variances of $I(r_{n+2})$ and $I(r_{n-1})$ remain similar to the uniform sampling rate ($\delta=1$) however the variances at the denser sampling rate becomes higher (i.e. $r_{n+1}$ and $r_{n}$). Most real light curves are not uniform, in these cases we suggest to select several almost uniform $(\delta\gtrsim 0.5)$ subsets of data and recover intensity profile for each subset ($\mathbf{I}_{i},i=1..N_{subset}$) and then take the average of all to find the final intensity profile.

To study the stability of solution in FEM in terms of the number of data points, we adapt a constant ${\sigma_{i}}$ and let $N_{1}$ to change from $4$ to $90$. We take a uniform cadence of $\Delta t=(t_{transit}-t_{0})/{(N_{1}-1)}$ to generate light curves. The results are shown in Figure (13) where increasing the number of data points results in reconstruction of poor LD profile. This effect results from numerical errors and this problem is well known in inverse problems (Craig & Brown, 1986). A larger data sets leads to larger Global Stiffness Matrix (GSM) and there is more chance that its rows and columns become nearly linear dependent and GSM becomes near singular.

On the other hand for smaller $N_{1}$, the average of variance decreases but absolute residual increases, which results in inaccurate intensity profiles. For the case of uniform data cadence and specified set of primary parameters, the optimum number of data during the transit is obtained to be $N_{1,optimum}=10$. This number depends on the uniformity of data points, error bars and coverage during the transit. We can stabilize the FEM solutions using regularization technique and we will discuss it in the next section.

In this section we find out that for reducing artificial dispersion of the recovered intensity profiles and to overcome instability problem we need to use regularisation technique. Additionally to reduce the dispersion that results from the denser parts of data points in the light curve ($\delta\ll 0.5$) we propose to select data subsets that are almost uniform $(\delta\gtrsim 0.5)$ and recover LD profile for each subset and take average over them to find the final LD profile. We find out that light curves with smaller impact parameter and higher quality photometry give more reliable intensity profiles.

4.2 Simulation Details: Regularized FEM

In the previous section we studied the effects of photometric error bars and increasing the number of data points on the quality of recovered LD profiles. We found out that photometric error bars resulted in artificial dispersions in the recovered intensity profiles and increasing the number of data points increase the dispersions. We can overcome both of these problems using regularisation techniques. The general idea is that instead of solving the equation (27) that returns answers with minimum residual $R^{2}=\sum_{i=1}^{i=N_{1}}(A_{data,i}-A_{I_{FEM},i})^{2}$, we minimize an objective function that combines $R^{2}$ with other physical assumptions such as minimizing dispersion of the intensity profile in the nodes (in another word, minimizing the norm of the second derivative of the solution). The objective function is defined as follows (Craig & Brown, 1986):

[TABLE]

in which $\lambda$ is the smoothing parameter. To minimize this objective function we write this equation in the matrix form and differentiate it with respect to $I_{\lambda}$. The result is as follows:

[TABLE]

where $\mathbf{Q}$ is the smoothing matrix. The answer with $\lambda=0$ is the classical solution of Equation (27), as shown in Figure (9) and the larger $\lambda$ yields smoother solutions. At large $\lambda$s the second term is more prominent and the resultant profiles approach to a straight line since the second derivative of a straight line is zero. So further smoothing of the solution provides a straight line for the intensity profile which is not our desire. We change the variable of the intensity profile from $r$ to $\mu=\sqrt{1-r^{2}}$ in the second term of equation (31) which results in the extension of limb area and it is suitable for the recovering the limb of the source star. The second finite derivative of intensity profile in terms of the new variable, $\mu$ is

[TABLE]

[TABLE]

in which $\Delta_{i}=\mu(r_{i+1})-\mu(r_{i})$. We note that $\Delta_{i}$s are not equal. According to equation (32) there exists a solution of $\mathbf{I_{\lambda}}$ for each $\lambda$. To recover LD profile for each light curve we need to find the optimum value of $\lambda$ using simulations. We generate simulated light curve data with a similar cadence and photometric error bars as in the real data. Also we implement our known limb darkening profile ($I_{input}$) in the light curve.

In the real light curves, we have non-uniform data cadence and some of data points are very close to each other ($\delta<0.001$). From our investigation in the previous section, to avoid errors from the non-uniform sampling in the recovered profiles, we need to generate almost uniform data subsets. We recover LD profiles for each of the subsets and take the average over all recovered profiles. In order to select almost uniform subsets of data we produce histogram of data with number of $N_{hist}$ bins for each light curve. Then, we randomly select one data from each of the bins of the histogram to produce an almost uniform ($median(\delta)\gtrsim 0.5$) subset of data. We generate an ensemble of data subsets to use all the initial data in the observed light curve. If the original light curve satisty the condition $median(\delta)\gtrsim 0.5$ we use all data at once and skip bininng procedure. The final result is the average of all the profiles (i.e. $I_{FEM,\lambda}=1/N_{s}\sum_{i=1}^{N_{s}}\mathbf{I_{\lambda,i}})$.

Here, we set $N_{hist}=100$, also since $\lambda$ is a free parameter, we let this parameter change with a power-law function as $\lambda_{i}=10^{i}$ (instead of a linear function). We choose the parameter of $i$ in the range of $i\in[-5,3]$ with $160$ steps. For each $\lambda_{i}$ we calculate dispersion ($D_{\lambda_{i}}$) and squared weighted residual ($R_{\lambda_{i}}$) of the corresponding $I_{FEM,\lambda_{i}}$:

[TABLE]

in which $I_{FEM,\lambda_{i},j}$ is the $j$th node of $I_{FEM,\lambda_{i}}$ and

[TABLE]

in which $A(I_{FEM,\lambda_{i}})_{j}$ is the recovered amplification with the parameter of $\lambda_{i}$. A profile ($I_{FEM,\lambda_{i}})$ with $R_{\lambda_{i}}\simeq 1$ would provide the reasonable recovered profile. On the other hand, a profile with $R_{\lambda_{i}}\eqsim 0$ is an oscillatory solution and it generates a light curve that goes through all data points.

To find the optimum value of $\lambda$ we use a combined condition using both $D_{\lambda_{i}}$ and $R_{\lambda_{i}}$ functions. First we find minimum dispersion $min(D_{\lambda_{i}})$ and check $R_{\lambda_{i}}$ to satisfy the following condition of

[TABLE]

where $\alpha>1$ and this parameter corresponds to a neighborhood around the minimum value of $min(D_{\lambda_{i}})$. In order to specify the $\alpha$ value, we use a test with simulating light curves where data points have the same cadence and error bars as in the real light curve while functions of $I_{input}$s are different. We call these as the control simulations and take one of the input profiles to be a square root function:

[TABLE]

and the other one to be

[TABLE]

where $f_{1}$ and $f_{2}$ are chosen in such a way to normalize the overall flux of source star to be one.

These two function behave differently in $\mu$-space. In the first function from equation (36), $d^{2}I(\mu)/d\mu^{2}<0$ while the second function in equation (37) has reflection point with $d^{2}I(\mu)/d\mu^{2}=0$. From the first profile in equation (36) we obtain the best solution and from the second profile we exclude the the recovered profiles with the overly smoothed solutions; such a solution is very close to a straight line in the $\mu$ space regardless of the input profile. Figures (14) and (15) show the result of this procedure for a microlensing event with two different number of data points in the light curve and a high cadence near limb. We use all data at once and skip the binning procedure. Here is the parameters of the light curve: $\sigma_{A}/{A}=0.005$ and the lensing parameters as $\rho_{\ast}=0.02,p=0.005,t_{E}=7.29\text{ days},t_{0}=0,c=0.3,d=0.6$.

4.3 Results

We apply FEM to a sample of microlensing light curves with finite-source effect to obtain the limb darkening of the source stars. This sample is discovered by survey groups of MOA and/or OGLE and alerted to the follow-up collaborations of PLANET, $\mu$FUN, RoboNet, MiNDSTEp. These events have been analyzed by Choi et al. (2012) where they assumed a linear standard limb darkening function and found parameter of the profile from fitting to the light curves in a specific filter. Here, we choose $7$ light curves from six events that have a good data coverage near the peak. We note that the limb darkening is a wavelength dependent phenomenon and data of each observatory with different filters should be analyzed separately. The selected light curves are listed in Table (1) with the best values of the lensing parameters in Choi et al. (2012). Also from the standard $\chi^{2}$-fitting, we obtain the blending (i.e. $b$) and the baseline of the source stars (i.e. $m_{BL}$). We used high quality data as in Choi et al. (2012) with $t<3t_{E}$, $\sigma_{m}<0.1$ and $residual<3\sigma_{m}$.

Figures (16-22) represents the light curves of $7$ events in $\ell$-space (see equation (3)) together with their control simulations in the left-down panels and corresponding reconstructed limb darkening of the source stars in r-space at the middle panels and in the $\mu$-space at the right panels. In these figures, the grey region corresponds to the area that there is no observed data. In the white region we have data observed by the telescopes. The light grey area around the recovered profile is the standard deviation that is obtained from ten simulated light curves with the cadence and error bars as the real data. The list of events and corresponding analysis are as follows:

OGLE-2004-BLG-254: The OGLE-2004-BLG-254 is a bulge event discovered by the OGLE survey and alerted for the follow-up telescopes of PLANET and $\mu$FUN collaborations. We use SAAO I-band of PLANET collaboration which has good data coverage and better photometric quality in our analysis. This event was analysed for the first time by Cassan et al. (2006) and later by Choi et al. (2012). They determined the source type to be a KIII star, based on its location in the CMD151515Color-Magnitude Diagram. They both obtained linear LDC (i.e. $u_{LLD}$ in equation.(30)) with $\chi^{2}$ fiting. Cassan et al. (2006) used Gould (1994) approximation for a single lens event with an extended source of uniform intensity. They derived $u_{LLD}=0.45^{+0.03}_{-0.06}$ for the I-band observation of SAAO, the source size $\rho_{\ast}=0.04\pm 0.0002$ and lens parameters ($t_{0}=3166.8194\pm 0.0002,u_{0}=0.0046\pm 0.0008,t_{E}=13.23\pm 0.05$). Choi et al. (2012) also used inverse-ray shooting technique to compute light curve. They derived $u_{LLD}=0.55\pm 0.06$ in I-band observation of SAAO, the source size $\rho_{\ast}=0.0418\pm 0.0004$ and lens parameters ($t_{0}=3166.823\pm 0.001,u_{0}=0.0111\pm 0.0004,t_{E}=12.84\pm 0.09$). In our analysis we adapt these lens parameters for our light curve. Since we have good data coverage in this light curve, we choose smaller bins around the limb to recover the limb darkening with better accuracy. Figure (16) represents the recovered profile intensity profile for this event.

MOA-2007-BLG-233/OGLE-2007-BLG-302: This event was discovered and alerted by both MOA and OGLE surveys, the follow-up observations carried out by $\mu$FUN, PLANET and MiNDSTEp collaborations. The observational data suitable for our method (as discussed before) are Las Campanas Observatory’s I-passband (OGLE-I) and Mt. John Observatory’s R-passband (MOA-R). This event was analysed by Choi et al. (2012). The source star is a GIII star and using a best fit parametric method they derived $u_{LLD}=0.56\pm 0.02$ for the MOA-R data and $u_{LLD}=0.53\pm 0.04$ for the OGLE-I data. The recovered LD profile by regularised FEM is shown in Figure (17) in OGLE I-band Figure (18) in MOA R-band. We note that for this event OGLE has almost uniform coverage of the light curve, we skip the binning procedure.

MOA-2010-BLG-436: This event is observed by MOA survey in R-band. It was analysed by Choi et al. (2012). The type of the source star was not determined due to bad quality of data in V-band however the LDC was determined to be $u_{LLD}=0.52\pm 0.10$ from MOA-R data. For this event we have uniform data coverage and we skip the binning procedure. Figure (19) shows the LD from the regularised FEM.

MOA-2011-BLG-093: This event was discovered and alerted by both MOA and OGLE surveys. The follow-up observations were carried out by $\mu$FUN, PLANET, RoboNet and MiNDSTEp collaborations. The observational data suitable for our method is Canopus Hill Observatory’s I-passband (Canopus-I). This event was analysed by Choi et al. (2012) and they determined the source star as GIII and the best LDC as $u_{LLD}=0.51\pm 0.03$. We had a good data coverage for this event around the limb and we choose smaller binning for this part. Figure (20) shows the LD from the regularised FEM.

MOA-2011-BLG-300/OGLE-2011-BLG-0990: This event was discovered and alerted by both MOA and OGLE surveys. The follow-up observations carried out by $\mu$FUN, PLANET collaborations. The observational data suitable for our method is taken from Observatorio do Pico dos Dias, I-passband (Pico-I). This event was analysed by Choi et al. (2012). The type of the source star was not determined, but they derived the LDC using a best fit parametric method to be $u_{LLD}=0.56\pm 0.04$. Figure (21) shows the LD from the regularised FEM.

MOA-2011-BLG-325/OGLE-2011-BLG-1101: This event was discovered and alerted by both MOA and OGLE surveys. The follow-up observations carried out by $\mu$FUN, PLANET, RoboNet and MiNDSTEp collaborations. We use the Liverpool Telescope’s I-passband observation (LT-I) has better data coverage during transit over the source star. Figure (22) shows the LD profile from the regularised FEM.

Finally, in order to check the consistency of our recovered profiles with standard LD profiles, we compare our recovered profiles (from regularized FEM) with the linear limb darkening (LLD) (see equation 30) and the square root LD (see equation 36) functions. Table (2) represents the result of this comparison where we derive the parameters of these two LD profiles. We also compare our results with the LLD parameters that have been obtained from direct fitting to the light curves (Choi et al., 2012; Cassan et al., 2006). The result of our analysis from Table (2) shows that events MOA-2010-BLG-436 and MOA-2011-BLG-093 are more consistent with the square-root profile than the linear limb darkening profile. The rest of events are consistent with the linear LD profile.

5 Conclusion

In the finite-source effect during the microlensing events, the lens can transit over the source star for the small impact parameters. The intensity profile of the source star affect on the light curve around the peak on the microlensing events. Studying the light curve around the peak enables us to investigate the limb darkening (LD) profile of the source star.

In this paper we used the classical Finite Element Method (FEM ) as an inversion tool to recover the LD equation from the observed microlensing light curve. We found out that this method has the two main problems of (i) instability with increasing the number of data points, (ii) a disperse recovered solution. In the second step we used the regularized FEM where we minimize an objective function that combines the residual of observed data with the reconstructed light curve as well as minimizing dispersion of the intensity profile in the nodes, see equation (31). This method could resolve the problems we faced in the classical FEM.

Finally we applied our method to single lens transit microlensing events and select data points from the light curve that has good coverage around the peak. We applied the regularized FEM to the following microlensing events of OGLE-2004-BLG-254 (SAAO-I band), MOA-2007-BLG-233/OGLE-2007-BLG-302 (OGLE-I, MOA-R), MOA-2010-BLG-436 (MOA-R), MOA-2011-BLG-93 (Canopus) I-band, MOA-2011-BLG-300/OGLE-2011-BLG-0990 (Pico) I-band and MOA-2011-BLG-325/OGLE-2011-BLG-1101 (LT) I-Band and compare our model-independent results with that of other works where a simple model was assumed for LD of the star. We note that the advantage of this method would be the reconstruction of LD of source star without pre-assumption about the intensity profile of the source star.

Acknowledgments

We thank C. Han and J. -Y. Choi for providing us the microlensing data. We also thank M.A. Jalali and M. Dominik for useful discussion and sharing their insights with us. Also, we thank anonymous referee to his/her useful comments improving this work.