A universal approach to the calculation of the transit light curves

We have developed a universal approach to compute accurately the brightness of eclipsing binary systems during the transit of a planet in front of the stellar disk. This approach is uniform for all values of the system parameters and applicable to most limb-darkening laws used in astrophysics. In the cases of linear and quadratic limb-darkening laws we obtained analytical expressions for the light curve and its derivatives in terms of elementary functions, elliptic integrals and piecewise-defined function of one variable. In the cases of logarithmic and square root laws of limb darkening the flux and its derivatives were expressed in terms of integrals which can be efficiently computed using Gaussian quadrature formula, taking into account singularities of the integrand.


INTRODUCTION
Recently several authors have developed algorithms for the calculation of transit light curves, see, e.g., Mandel&Agol (2002), Pal (2008), Pal (2012). However, the problem of calculation of the light curves is still relevant, because the existing algorithms are not applicable to all values of the system parameters for some limb-darkening laws. Besides, they do not allow sufficiently accurate calculations of the light curves for some limb-darkening laws. In addition, calculations of derivatives of the light-curve as a function of system parameters is important, because they can be used to solve the inverse problem of interpretation of the light curve.
The paper Mandel&Agol (2002) contains an analytical expression of the light curve by elliptic integrals, for the cases of linear and quadratic limbdarkening laws. In doing so, 13 variants of relations between the parameters are considered. For other limb-darkening laws (law of square-root and its power)only an approximate method of light-curve calculation at the radius of the planet more than 10 times smaller than the radius of the star is being used. In this case, the accuracy is 2% of the depth of the ⋆ E-mail: marat@sai.msu.ru † E-mail: ngostev@sai.msu.ru eclipse. In the paper Pal (2012)), directly, there is only an expression of the light curve in the linear and quadratic limb-darkening law, and the derivatives of the light curve are calculated by difference methods (This work contains no direct analytical expressions for the derivatives), that is less favorable in terms of time and accuracy of the computation. In addition, none of the above works the logarithmic law of darkeningbare not considered, which is the most preferred for early-type stars (Klinglesmith&Sobieski (1970) and Van Hamme (1993)).
The approach presented in this paper allows us to calculate a light curve and with almost machine accuracy for any values of the parameters (including near singularities). Binary system parameters are the radii of the components, and the distance between the centers of the components in the projection on the picture plane. In general, the algorithm is uniform for all values of the system parameters, which significantly facilitates its implementation. We obtained analytical expressions for the transit light curve of the eclipsing binary system and for its derivatives in the cases of the linear and quadratic limb-darkening laws. These quantities are expressed in terms of piecewise-defined function of one variable and incomplete elliptic integrals, which can be computed with effective methods proposed by Carlson (1994). In the cases of the logarithmic limb-darkening law and the square-root limb-darkening law the light function is expressed through integrals that can be efficiently computed using Gaussian quadrature formula. In this respect, the computation time of the light curve is not much more than the computation time by analytical expressions.

MODEL DESCRIPTION
We considered the model of the eclipse of a spherically symmetric star with thin atmosphere by another spherical opaque component (the other spherical star or a spherical planet). Fig. 1 shows the geometry of the stellar disks in eclipse.
The brightness at the point of the disk of the eclipsed star with polar coordinate ρ is given by: Here J(0) is the brightness at the center of this stellar disk, where functions f k are such that f k (1) = 0, are defined by the law of limb-darkening in question, and Λ k are the coefficients of limb-darkening. In this paper we consider the following frequently used limb-darkening laws: Linear limb-darkening law, for which f (µ) = = Λ l f l (µ) = Λ l (1 − µ); Square law of limb-darkening, which is characterized by the presence of the term Λqfq(µ) = Λq(1 − µ) 2 in the expression for f ; Logarithmic limb-darkening law, which is characterized by the presence of the term ΛLfL(µ) = = −ΛLµ ln µ in the expression for f ; Square root limb-darkening law, which is characterized by the presence of the term ΛQfQ(µ) = = ΛQ(1− √ µ) in the expression for f . Также полученные далее результаты можно очевидным образом обобщить на случай закона потемнения к краю, характеризуемого членом µ l в выражении для яркости, где l -положительное целое число.

GENERAL INTEGRAL FORMULA FOR THE FLUX
The decrease of the flux of the binary system due to eclipse is: (2) where L F is the unobscured flux of the binary system, L is the obscured flux of the binary system, i.e. the light-curve value, S(D) is the area of overlapping disks, R is radius-vector of the point on the stellar disk.
To calculate the integral (2) we introduce the functions: and Then The relation (5) is obtained naturally by noting that A z = Re arccos z, Q x = Re √ z for complex number z with Im z = 0 and for the functions of complex argument arccos and √ · that are analytic continuations of the inverse cosine and square root of a real argument. Analyticity region is such that −π < arg z π for each z.
In the polar coordinate system the region of integration S(D) is given by: cos ϕ .
In case of integration (2) with respect to coordinate ϕ, for the values of ρ, which satisfy 1 variable ϕ takes the values such that For the values of ρ, for which Here R * is the radius of the eclipsed star, Ro is the radius of the eclipsing component, D is the distance between the centers of the disks of the components, ρ, Ψ are respectively the polar radius and the polar angle of a point on the disk of the eclipsed star. The origin is located at the center of the eclipsed star, the polar angle is measured counterclockwise from the radius-vector connecting centers of the star and the transiting component.
holds for every value of ϕ. In this case, integration with respect to ϕ runs from −π to π.
For the values of ρ, for which the last inequality (6) is not satisfied for any values of ϕ. Formally, at these values of ρ both integration limits by ϕ are set equal to zero. Next, using the notation (3) and introducing the function integral in (2) can be rewritten as: where r = Ro R * , δ = D R * , and Note that ∆L(δ, r) is the value of the decrease of the flux of the binary system when radius and brightness at the center of eclipsed star equals unity, radius of the second (eclipsing) component equals r and the distance between centers of disks equals δ. In view of (1) we can express the decrease of the flux as linear combination with limb-darkening coefficients: Unobscured flux L f of the star is when R is the radius of the star, and (10) When both components of the binary system are the stars, unobscured flux L F of the binary system is the sum of L f for each star. For binary star and planet L F equals L f for star.
Let g be a function that g(ρ) is a separate linear term in the expression for I( √ ρ) (given by (1)), We consider the integral of the general form, which is a contribution to ∆L(δ, r) caused by the term g(ρ) in the expression for I( √ ρ): We note that for δ > 0, r > 0 Using integration by parts, we obtain where (5) is used for differentiating Ψ.
By differentiating in (11) integrand with respect to δ and r we similarly find an expression for the corresponding partial derivative ∆Lg: The contribution to the L f caused by the term g(ρ) in the expression for I( √ ρ):
A universal approach to the calculation of the transit light curves 7 and ∂∆L3(δ, r) ∂r = 4r For term with square-root limb-darkening coefficients in (9): The formulas obtained for the square root, it is easy to generalize to the case limb-darkening law contained in the expression for the brightness the term µ l , where l is an odd positive number, putting in (15), (16) add (17) (1 − x) 1+l/4 . For an even l of non-multiple 4, light curve and its derivative can be expressed by elliptical integrals similarly to as the formulas (22)-(24) were obtained. If l is divisible by 4 light curve and its derivative can be expressed by elementary functions similarly to as (26)-(28) were obtained for quadratic limb darkening.

NUMERICAL CALCULATION OF INTEGRALS
Thus the calculation of the brightness for the logarithmic and square-root limb-darkening law is reduced to the calculation of integrals P for P L r where vi(x) are some trigonometric polynomials, n > 0, k > 0. We denote the maximum degree of these trigonometric polynomials as τ . K(y) = √ y ln y or K(y) = 4 √ y γ , respectively, has a logarithmic or fractional power singularity at t = 0.
In the case of calculating of P Q r P Q δ , P Q r , we put γ = 1. For the limb darkening of the general form, which is characterized by the presence of the term µ l in the expression for brightness (odd l) it is enough to put γ = l mod 4.
By applying Gaussian quadrature formula, we can find the numerical value of the integrals with high precision, producing a relatively small number of elementary computations (the amount of computation of the integrand is proportional to required number of significant digits). But at the same time, an integrable function must satisfy certain conditions. In particular, this can be achieved if the higher derivatives of the integrand (or its non-singular component) is uniformly bounded on the section of integration. To reduce the computation of the integrals (45) and (46) to computation of the integrals that satisfy the above conditions, we divide the interval of integration (0, Xi − Xi+1 1 max{τ, 2} for all i < M and k . (48) In the case of (45) we also require that the following inequality: n + cos Xi n + cos Xi+1 2 for all i < M and k .
(47)-(50) can be used as recurrent relation, allowing us to construct the sequence Xi. Thus,P For fixed values of n and k the two last integrals can be represented in general form: where U (x) = V (k, x) n + sin 2 x for (51) and U (x) = V (k, x) for (52). By linear substitution of variable of integration in (53), we turn to the integration from zero to unity: (53) In this form, Pi can be computed by applying the Gaussian quadrature formula: Here ω(t) > 0 ∀t ∈ (0, 1), nodes ti are the roots of the polinomial HN (t), where {Hi} is the system of orhtogonal polynomials with weight ω in the interval (0, 1): w l can be found as the solution of the system of N linear algebraic equations, which can be obtained if (54) and replace the approximate equality with exact equality. N can be adjusted so as to ensure the required accuracy of calculation of P0i(n, k), Pdi(k) and can be the same for all values of i, n, k. N is of the same order of magnitude as the number of significant digits in the result, and this allows to calculate the integral with the required accuracy in a reasonable time. So, after calculation of the roots of polynomials x l and weights w l (this may take a while), we can re-use them for computing P0i(n, k), Pdi(k) for all i, n, k.
In the case of i > 0 or of k > 1 we put in (54): Note that in this case Hi(t) ≡ Pi(2t − 1) where Pi are Legendre polynomials.
Let S1 = N l=1 w l h(t l ).
Calculations show that in all cases of the applications of the Gauss quadrature accuracy of 19 significant decimal digits (corresponding to 80-bit machine numbers) can be achieved by choosing the N to be 14. Value sets of points ti and weights wi corresponding to each of the considered forms of the function ω, can be downloaded from the Internet, along with other materials (see Conclusion).