A universal approach to the calculation of the transit light curves
aa r X i v : . [ a s t r o - ph . S R ] A p r Mon. Not. R. Astron. Soc. , 1– ?? (2012) Printed 5 сентября 2018 г. (MN L A TEX style file v2.2)
A universal approach to the calculation of the transitlight curves
M.K. Abubekerov ⋆ and N.Yu. Gostev † Lomonosov Moscow State University, Sternberg Astronomical Institute, Russia
Accepted 2013 April 2. Received 2013 March 22; in original form 2012 November 20
ABSTRACT
We have developed a universal approach to compute accurately thebrightness of eclipsing binary systems during the transit of a planet infront of the stellar disk. This approach is uniform for all values of thesystem parameters and applicable to most limb-darkening laws used inastrophysics. In the cases of linear and quadratic limb-darkening lawswe obtained analytical expressions for the light curve and its derivativesin terms of elementary functions, elliptic integrals and piecewise-definedfunction of one variable. In the cases of logarithmic and square root lawsof limb darkening the flux and its derivatives were expressed in termsof integrals which can be efficiently computed using Gaussian quadratureformula, taking into account singularities of the integrand.
Key words: stars, binary systems, eclipse.
Recently several authors have developedalgorithms 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 curvesis still relevant, because the existing algorithms arenot applicable to all values of the system parametersfor some limb-darkening laws. Besides, they do notallow sufficiently accurate calculations of the lightcurves for some limb-darkening laws. In addition,calculations of derivatives of the light-curve as afunction of system parameters is important, becausethey can be used to solve the inverse problem ofinterpretation of the light curve.The paper Mandel&Agol (2002) contains ananalytical expression of the light curve by ellipticintegrals, for the cases of linear and quadratic limb-darkening laws. In doing so, 13 variants of relationsbetween the parameters are considered. For otherlimb-darkening laws (law of square-root and itspower)only an approximate method of light-curvecalculation at the radius of the planet more than 10times smaller than the radius of the star is being used.In this case, the accuracy is 2% of the depth of the ⋆ E-mail: [email protected] † E-mail: [email protected] eclipse. In the paper Pal (2012)), directly, there isonly an expression of the light curve in the linear andquadratic limb-darkening law, and the derivatives ofthe light curve are calculated by difference methods(This work contains no direct analytical expressionsfor the derivatives), that is less favorable in terms oftime and accuracy of the computation. In addition,none of the above works the logarithmic law ofdarkeningbare not considered, which is the mostpreferred for early-type stars (Klinglesmith&Sobieski(1970) and Van Hamme (1993)).The approach presented in this paper allows usto calculate a light curve and with almost machineaccuracy for any values of the parameters (includingnear singularities). Binary system parameters are theradii of the components, and the distance betweenthe centers of the components in the projectionon the picture plane. In general, the algorithm isuniform for all values of the system parameters,which significantly facilitates its implementation. Weobtained analytical expressions for the transit lightcurve of the eclipsing binary system and for itsderivatives in the cases of the linear and quadraticlimb-darkening laws. These quantities are expressed interms of piecewise-defined function of one variable andincomplete elliptic integrals, which can be computedwith effective methods proposed by Carlson (1994). Inthe cases of the logarithmic limb-darkening law and c (cid:13) M. K. Abubekerov and N. Yu. Gostev the square-root limb-darkening law the light functionis expressed through integrals that can be efficientlycomputed using Gaussian quadrature formula. In thisrespect, the computation time of the light curve is notmuch more than the computation time by analyticalexpressions.
We considered the model of the eclipse of aspherically symmetric star with thin atmosphereby another spherical opaque component (the otherspherical star or a spherical planet).Fig. 1 shows the geometry of the stellar disks ineclipse.The brightness at the point of the disk of theeclipsed star with polar coordinate ρ is given by: J ( ρ ) = J (0) I (cid:18) ρR ∗ (cid:19) . Here J (0) is the brightness at the center of this stellardisk, I ( r ) = (1 − f ( µ ( r ))) , (1) µ ( r ) = p − r f ( µ ) = X k Λ k f k ( µ ) , where functions f k are such that f k (1) = 0 , are definedby the law of limb-darkening in question, and Λ k arethe coefficients of limb-darkening.In this paper we consider the following frequentlyused limb-darkening laws:Linear limb-darkening law, for which f ( µ ) == Λ l f l ( µ ) = Λ l (1 − µ ) ;Square law of limb-darkening, which ischaracterized by the presence of the term Λ q f q ( µ ) = Λ q (1 − µ ) in the expression for f ;Logarithmic limb-darkening law, which ischaracterized by the presence of the term Λ L f L ( µ ) == − Λ L µ ln µ in the expression for f ;Square root limb-darkening law, which ischaracterized by the presence of the term Λ Q f Q ( µ ) == Λ Q (1 −√ µ ) in the expression for f . Также получен-ные далее результаты можно очевидным образомобобщить на случай закона потемнения к краю, ха-рактеризуемого членом p µ l в выражении для яр-кости, где l – положительное целое число. The decrease of the flux of the binary system dueto eclipse is: L F − L ( D, R ∗ , R o ) = ∆ L ( D, R ∗ , R o ) = Z Z S ( D ) J ( | R | ) d R , (2)where L F is the unobscured flux of the binary system, L is the obscured flux of the binary system, i.e. thelight-curve value, S ( D ) is the area of overlappingdisks, R is radius-vector of the point on the stellardisk.To calculate the integral (2) we introduce thefunctions: A x ≡ π, x < − x, − x , x > . (3)and Q x ≡ ( √ x, x > , x < . (4)Then d A xd x = Q (cid:18) − x (cid:19) . (5)The relation (5) is obtained naturally by notingthat A z = Re arccos z , Q x = Re √ z for complexnumber z with Im z = 0 and for the functions ofcomplex argument arccos and √· that are analyticcontinuations of the inverse cosine and square rootof a real argument. Analyticity region is such that − π < arg z π for each z .In the polar coordinate system the region ofintegration S ( D ) is given by: S ( D ) = ρ < R ∗ − π < ϕ πρ + D − R o ρD cos ϕ . (6)In case of integration (2) with respect tocoordinate ϕ , for the values of ρ , which satisfy (cid:12)(cid:12)(cid:12)(cid:12) ρ + D − R o ρD (cid:12)(cid:12)(cid:12)(cid:12) variable ϕ takes the values suchthat ρ + D − R o ρD cos ϕ ⇔ | ϕ | arccos (cid:18) ρ + D − R o ρD (cid:19) . Hence, the integration over ϕ isfrom − arccos (cid:18) ρ + D − R o ρD (cid:19) to arccos (cid:18) ρ + D − R o ρD (cid:19) .For the values of ρ , for which ρ + D − R o ρD < − inequality ρ + D − R o ρD < cos ϕ c (cid:13) , 1– ?? universal approach to the calculation of the transit light curves d R* R o D Figure 1.
The model of eclipsing binary system. The projection on the picture plane. Here the smaller component is a staror an exoplanet. The geometry of stellar disks in eclipse. Here R ∗ is the radius of the eclipsed star, R o is the radius of theeclipsing component, D is the distance between the centers of the disks of the components, ρ , Ψ are respectively the polarradius and the polar angle of a point on the disk of the eclipsed star. The origin is located at the center of the eclipsedstar, the polar angle is measured counterclockwise from the radius-vector connecting centers of the star and the transitingcomponent. holds for every value of ϕ . In this case, integrationwith respect to ϕ runs from − π to π .For the values of ρ , for which ρ + D − R o ρD > ,the last inequality (6) is not satisfied for any valuesof ϕ . Formally, at these values of ρ both integrationlimits by ϕ are set equal to zero.Next, using the notation (3) and introducing thefunction Ψ( D, x, y ) ≡ A (cid:18) x + D − y x D (cid:19) , integral in (2) can be rewritten as: ∆ L ( D, R ∗ , R o ) = R ∗ Z ρdρ Ψ( D, ρ, R o ) Z − Ψ( D, ρ, R o ) dϕJ ( ρ ) == 2 J (0) R ∗ Z ρ Ψ( D, ρ, R o ) I (cid:18) ρR ∗ (cid:19) dρ == J (0) R ∗ ∆ L (cid:18) DR ∗ , , R o R ∗ (cid:19) = J (0) R ∗ ∆L ( δ, r ) , (7)where r = R o R ∗ , δ = DR ∗ , and ∆L ( δ, r ) = 2 Z ρ Ψ( δ, ρ, r ) I ( ρ ) dρ == Z Ψ( δ, √ ρ, r ) I ( √ ρ ) dρ . (8)Note that ∆L ( δ, r ) is the value of the decrease of theflux of the binary system when radius and brightness c (cid:13) , 1– ?? M. K. Abubekerov and N. Yu. Gostev at the center of eclipsed star equals unity, radius ofthe second (eclipsing) component equals r and thedistance between centers of disks equals δ . In viewof (1) we can express the decrease of the flux as linearcombination with limb-darkening coefficients: ∆L ( δ, r ) = ∆L ( δ, r ) + Λ l ∆L l ( δ, r )++ Λ q ∆L q ( δ, r ) + Λ L ∆L L ( δ, r ) + Λ Q ∆L Q ( δ, r ) . (9)Unobscured flux L f of the star is L f = J (0) R π L f when R is the radius of the star, and L f = Z I ( √ ρ ) dρ == L f + Λ l L fl + Λ q L fq + Λ L L fL + Λ Q L fQ . (10)When both components of the binary system are thestars, unobscured flux L F of the binary system is thesum 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 linearterm in the expression for I ( √ ρ ) (given by (1)), g ( − is one of its primitives: g ( ρ ) = dg ( − ( ρ ) dρ . Weconsider the integral of the general form, which is acontribution to ∆L ( δ, r ) caused by the term g ( ρ ) inthe expression for I ( √ ρ ) : ∆L g ( δ, r ) = Z Ψ( δ, √ ρ, r ) g ( ρ ) dρ . (11)We note that for δ > , r > ρ → Ψ( δ, √ ρ, r ) = π Θ( r − δ ) , where Θ( t ) ≡ , t > , t = 00 , t < , Using integration by parts, we obtain ∆L g ( δ, r ) = Ψ( δ, , r ) g ( − (1) − π Θ( r − δ ) g ( − (0) −− Z ( δ − r − ρ ) g ( − ( ρ )2 ρ ×× Q (cid:0) ρ − ( δ − r ) (cid:1) (cid:0) ( δ + r ) − ρ (cid:1) ! dρ , (12)where (5) is used for differentiating Ψ .For non-negative r and δ the integrand in (12) isnon-zero only if ( δ − r ) < ρ < ( δ + r ) . Thereforethe integral in (12) is zero if | δ − r | > . And if | δ − r | < , integrating can be performed over theinterval (cid:0) ( δ − r ) , min (cid:0) ( δ + r ) , (cid:1)(cid:1) . In this interval arccos (cid:18) δ + r − ρ δ r (cid:19) is a monotone function of ρ , so we can perform change of variable in integration inthe following way: x = arccos (cid:18) δ + r − ρ δ r (cid:19) . (13)Then ρ = δ + r − rδ cos x . (14)Integration with respect to x will be perfomedover the interval , arccos δ + r − min (cid:0) ( δ + r ) , (cid:1) δ r !! . Taking into account the fact that ddρ arccos (cid:18) δ + r − ρ δ r (cid:19) == 1 q(cid:0) ρ − ( δ − r ) (cid:1) (cid:0) ( δ + r ) − ρ (cid:1) , and arccos δ + r − min (cid:0) ( δ + r ) , (cid:1) δ r ! == A (cid:18) δ + r − δ r (cid:19) = Ψ( δ, r, , we obtain: ∆L g ( δ, r ) = Ψ( δ, , r ) g ( − (1) − π Θ( r − δ ) g ( − (0)++Ψ( δ, r, Z ( r − rδ cos x )) g ( − ( δ + r − rδ cos x ) δ + r − rδ cos x dx (15)Expression (15) is obtained assuming | δ − r | < .If | δ − r | > value of Ψ( δ, r,
1) = 0 and integral in(15) vanishes. As noted above, when | δ − r | > , theintegral in (12) is zero because the integrand vanishes.Thus, the expression (15) is valid for all positive valuesof δ and r .By differentiating in (11) integrand with respectto δ and r we similarly find an expression for thecorresponding partial derivative ∆L g : ∂ ∆L g ( δ, r ) ∂δ = − r Ψ( δ, r, Z cos x g ( δ + r − rδ cos x ) dx (16) ∂ ∆L g ( δ, r ) ∂r = 2 r Ψ( δ, r, Z g ( δ + r − rδ cos x ) dx . (17)The contribution to the L f caused by the term g ( ρ ) in the expression for I ( √ ρ ) : L fg = π ( g ( − (1) − g ( − (0)) . (18) c (cid:13) , 1– ?? universal approach to the calculation of the transit light curves The expression for the decrease of the flux dueto eclipse of the stellar disk with uniform brightness(with zero coefficients of limb darkening) can beobtained if we put in (15), (16) and (17) g ( x ) = 1 , g − ( x ) = x . Then: L f = π , ∆L ( δ, r ) = Ψ( δ, , r ) + Ψ( δ, r, r − Q ( δ, r ) , (19) ∂ ∆L ( δ, r ) ∂δ = − Q ( δ, r ) δ , (20)and ∂ ∆L ( δ, r ) ∂r = 2Ψ( δ, r, r . (21)Here Q ( δ, r ) ≡ Q (cid:0)(cid:0) − ( δ − r ) (cid:1) (cid:0) ( δ + r ) − (cid:1)(cid:1) . Putting g ( x ) = µ ( √ x ) = √ − x, g ( − ( x ) == − (1 − x ) , we get: L f = 2 π , ∆L ( δ, r ) = 2 π r − δ )++ Q (cid:18) − ( r − δ ) (cid:19) (cid:20) δ + r )3( δ − r ) ˆΠ −− (cid:0) δ − r ) + (1 − ( r − δ ) )(( r + δ ) − (cid:1) ˆ F (cid:21) ++ 29 Q (cid:0) − ( r − δ ) (cid:1) (7 r + δ −
4) ˆ
E . (22)Here ˆΠ ≡ Π (cid:18) − δr ( r − δ ) ; Ψ( δ, r, (cid:12)(cid:12)(cid:12)(cid:12) δr − ( r − δ ) (cid:19) ˆ F ≡ F (cid:18) Ψ( δ, r, (cid:12)(cid:12)(cid:12)(cid:12) δr − ( r − δ ) (cid:19) ˆ E ≡ E (cid:18) Ψ( δ, r, (cid:12)(cid:12)(cid:12)(cid:12) δr − ( r − δ ) (cid:19) , where F, E and Π are incomplete elliptic integrals ofthe first, second and third kind: F ( φ | m ) ≡ φ Z dθ p − m sin ( θ ) , E ( φ | m ) ≡ φ Z q − m sin ( θ ) , Π( n ; φ | m ) ≡ φ Z dθ (1 − n sin ( θ )) p − m sin ( θ ) . The efficient algorithms for their calculations weresuggested by Carlson (1994). When | δ − r | → or | δ − r | → , the limit of the term containingthe factor ˆΠ in (22) is equal to zero. Note that asimilar expression was obtained by Pal (2012) forthe integral (primitive) of the appropriately chosenvector field along the limb of the eclipsed component.However, application of this expression for calculationof the flux of the system requires further accountof its singularities. Expressions (22) and (19) givedirect algorithm for calculating of the brightness,and the possible singularities are taken into accountautomatically by piecewise smooth functions of onevariable A and Q ∂ ∆L ( δ, r ) ∂δ == − r Ψ( δ, r, Z cos( x ) r − ( δ − r ) − δr sin (cid:16) x (cid:17) dx == − δ Q (cid:0) − ( r − δ ) (cid:1) h(cid:0) ( r + δ ) − (cid:1) ˆ F ++ (cid:0) − δ − r (cid:1) ˆ E i . (23) ∂ ∆L ( δ, r ) ∂r == 2 r Ψ( δ, r, Z r − ( δ − r ) − δr sin (cid:16) x (cid:17) dx == 4 r Q (cid:0) − ( r − δ ) (cid:1) ˆ E . (24)For term with linear limb-darkening coefficients in (9): ∆L l ( δ, r ) = ∆L ( δ, r ) − ∆L ( δ, r ) , (25) L fl = L f − L f = − π . Assuming g ( x ) = x, g ( − ( x ) = x / we get: L f = π , ∆L ( δ, r ) = 12 Ψ( δ, , r )++ r (cid:0) δ + r (cid:1) Ψ( δ, r, −− (cid:0) δ + 5 r + 1 (cid:1) Q ( δ, r ) . (26)The partial derivatives ∆L : c (cid:13) , 1– ?? M. K. Abubekerov and N. Yu. Gostev ∂ ∆L ( δ, r ) ∂δ = 2 δr Ψ( δ, r, − δ + r + 12 δ Q ( δ, r ) , (27) ∂ ∆L ( δ, r ) ∂r = 2 r (cid:0) δ + r (cid:1) Ψ( δ, r, − rQ ( δ, r ) , (28) ∆L q ( δ, r ) = 2 ∆L ( δ, r ) − ∆L ( δ, r ) − ∆L ( δ, r ) , (29) L fq = 2 L f − L f − L f = π . Further, we note that Ψ( δ, r, π − A Q (cid:18) − ( δ − r ) δr (cid:19) . Assuming g ( x ) = √ − x ln(1 − x ) , g ( − ( x ) == (1 − x ) / (4 / − / − x )) we obtain for thelogarithmic limb-darkening law: L fL = − π , ∆L L ( δ, r ) = ∆L ( δ, r ) (cid:18) ln(4 δr ) − (cid:19) −− π δr )Θ( r − δ ) −− √ δr (cid:20) r P L (cid:18) ( δ − r ) δr , − ( δ − r ) δr (cid:19) −− δr P L (cid:18) ( δ − r ) δr , − ( δ − r ) δr (cid:19)(cid:21) (30)where P L ( n, k ) == π − A ( Q k ) Z (cid:0) k − sin x (cid:1) ln (cid:0) k − sin x (cid:1) n + sin x dx , (31) P L ( n, k ) == π − A ( Q k ) Z cos 2 x (cid:0) k − sin x (cid:1) ln (cid:0) k − sin x (cid:1) n + sin x dx (32)The partial derivatives of ∆L L : ∂ ∆L L ( δ, r ) ∂δ = ln(4 δr ) ∂ ∆L ( δ, r ) ∂δ −− r √ rδ P Lδ (cid:18) − ( δ − r ) δr (cid:19) (33) P Lδ ( k ) == π − A ( Q k ) Z cos 2 x p k − sin x ln (cid:0) k − sin x (cid:1) dx (34) ∂ ∆L L ( δ, r ) ∂r = ln(4 δr ) ∂ ∆L ( δ, r ) ∂r ++ 8 r √ rδ P Lr (cid:18) − ( δ − r ) δr (cid:19) (35) P Lr ( k ) == π − A ( Q k ) Z p k − sin x ln (cid:0) k − sin x (cid:1) dx (36)Assuming g ( x ) = √ − x, g ( − ( x ) = −
45 (1 − x ) / , weobtain the following expression for the case of squareroot limb-darkening law: L f = 4 π , ∆L ( δ, r ) = 4 π r − δ ) −− √ δr (cid:20) r P Q (cid:18) ( δ − r ) δr , − ( δ − r ) δr (cid:19) −− δr P Q (cid:18) ( δ − r ) δr , − ( δ − r ) δr (cid:19)(cid:21) , (37)where P Q ( n, k ) = π − A ( Q k ) Z (cid:0) k − sin x (cid:1) n + sin x dx , (38) P Q ( n, k ) = π − A ( Q k ) Z cos 2 x (cid:0) k − sin x (cid:1) n + sin x dx . (39)The partial derivatives of ∆L Q : ∂ ∆L ( δ, r ) ∂δ = − r √ rδ P Qδ (cid:18) − ( δ − r ) δr (cid:19) (40)where P Qδ ( k ) = π − A ( Q k ) Z cos 2 x (cid:0) k − sin x (cid:1) dx (41) c (cid:13) , 1– ?? universal approach to the calculation of the transit light curves and ∂ ∆L ( δ, r ) ∂r = 4 r √ rδ P Qr (cid:18) − ( δ − r ) δr (cid:19) (42)where P Qr ( k ) = π − A ( Q k ) Z (cid:0) k − sin x (cid:1) dx . (43)For term with square-root limb-darkening coefficientsin (9): ∆L Q ( δ, r ) = ∆L ( δ, r ) − ∆L ( δ, r ) , (44) L fQ = L f − L f = − π . The formulas obtained for the square root, itis easy to generalize to the case limb-darkening lawcontained in the expression for the brightness the term p µ l , where l is an odd positive number, putting in(15), (16) add (17) g ( x ) = p (1 − x ) l , g ( − ( x ) == −
44 + l (1 − x ) l/ . For an even l of non-multiple4, light curve and its derivative can be expressed byelliptical integrals similarly to as the formulas (22)–(24) were obtained. If l is divisible by 4 light curveand its derivative can be expressed by elementaryfunctions similarly to as (26)–(28) were obtained forquadratic limb darkening. Thus the calculation of the brightness for thelogarithmic and square-root limb-darkening law isreduced to the calculation of integrals P L , P L , P Lr P Lδ , P Lr , P L , P Q , P Qr P Qδ , P Qr (depending on parameters).These integrals can be represented in general form: ˜ P ( n, k ) = π − A ( Q k ) Z V ( k, x ) K (cid:0) k − sin x (cid:1) n + sin x dx , (45)for P L , P L , P Q , P Q or ¯ P ( k ) = π − A ( Q k ) Z V ( k, x ) K (cid:0) k − sin x (cid:1) dx , (46)for P Lr P Lδ , P Lr , P Qr P Qδ , P Qr . Here V ( k, x ) == s P i =1 u i ( k ) v i ( x ) where v i ( x ) are some trigonometricpolynomials, n > , k > . We denote the maximumdegree of these trigonometric polynomials as τ . K ( y ) = √ y ln y or K ( y ) = √ y γ , respectively, has a logarithmic or fractional power singularity at t = 0 .In the case of calculating of P Qr P Qδ , P Qr , we put γ = 1 .For the limb darkening of the general form, which ischaracterized by the presence of the term p µ l in theexpression for brightness (odd l ) it is enough to put γ = l mod 4 .By applying Gaussian quadrature formula, wecan find the numerical value of the integrals withhigh precision, producing a relatively small number ofelementary computations (the amount of computationof the integrand is proportional to required number ofsignificant digits). But at the same time, an integrablefunction must satisfy certain conditions. In particular,this can be achieved if the higher derivatives of theintegrand (or its non-singular component) is uniformlybounded on the section of integration. To reducethe computation of the integrals (45) and (46) tocomputation of the integrals that satisfy the aboveconditions, we divide the interval of integration (0 , π −A ( Q k )) sequence X > X > .... > X M , such that X = π − A ( Q k )) , X M = 0 and k − sin X i +1 k − sin X i for i > and k = 1 , (47) X i − X i +1 { τ, } for all i < M and k . (48)In the case of (45) we also require that thefollowing inequality: n + cos X i n + cos X i +1 for all i < M and k . (49)If k > , inequality from (47) also holds for i = 0 .If k < , then k − sin X ( X − X ) sin(2 X ) and X > X / . (50)(47)–(50) can be used as recurrent relation,allowing us to construct the sequence X i .Thus, ˜ P ( n, k ) = M − X i =0 ˜ P i ( n, k ) , ¯ P ( k ) = M − X i =0 ¯ P i ( k ) . where ˜ P i ( n, k ) = X i Z X i +1 V ( k, x ) K (cid:0) k − sin x (cid:1) n + sin x dx , (51)and ¯ P i ( k ) = X i Z X i +1 V ( k, x ) K (cid:0) k − sin x (cid:1) dx , (52)For fixed values of n and k the two last integrals can c (cid:13) , 1– ?? M. K. Abubekerov and N. Yu. Gostev be represented in general form: P i = X i Z X i +1 U ( x ) K (cid:0) k − sin x (cid:1) dx , where U ( x ) = V ( k, x ) n + sin x for (51) and U ( x ) = V ( k, x ) for (52).By linear substitution of variable of integration x ( t ) = X i +1 + t ( X i − X i +1 ) in (53), we turn to the integration from zero to unity: P i = ( X i − X i +1 ) Z U ( x ( t )) K (cid:0) k − sin ( x ( t )) (cid:1) dt . (53)In this form, P i can be computed by applying theGaussian quadrature formula: Z h ( t ) ω ( t ) dt ≈ N X l =1 w l h ( t l ) . (54)Here ω ( t ) > ∀ t ∈ (0 , , nodes t i are the roots ofthe polinomial H N ( t ) , where { H i } is the system oforhtogonal polynomials with weight ω in the interval (0 , : Z H l ( t ) H j ( t ) ω ( t ) dt = 0 for l = j .w l can be found as the solution of the system of N linear algebraic equations, which can be obtained ifwe put h ( t ) ≡ , h ( t ) ≡ t, . . . , h ( t ) ≡ t N in (54) andreplace the approximate equality with exact equality. N can be adjusted so as to ensure the requiredaccuracy of calculation of P i ( n, k ) , P di ( k ) and can bethe same for all values of i, n, k . N is of the same orderof magnitude as the number of significant digits in theresult, and this allows to calculate the integral withthe required accuracy in a reasonable time. So, aftercalculation of the roots of polynomials x l and weights w l (this may take a while), we can re-use them forcomputing P i ( n, k ) , P di ( k ) for all i, n, k .In the case of i > or of k > weput in (54): ω ( t ) = 1 ∀ t ∈ (0 , , h ( t ) == ( X i − X i +1 ) U ( x ( t )) K (cid:0) k − sin x ( t ) (cid:1) . Then P i ≈ N X l =1 w l h ( t l ) . Note that in this case H i ( t ) ≡ P i (2 t − where P i areLegendre polynomials.In the case of i > , k < and logarithmiclimb-darkening law ( K ( y ) = √ y ln y ) we represent theintegrand from (53) in the form: U ( x ( t )) √ − t s k − sin ( x ( t ))1 − t ×× (cid:20) ln (cid:18) k − sin ( x ( t ))1 − t (cid:19) + ln(1 − t ) (cid:21) Next, we put in (54): ω ( t ) = √ − t ∀ t ∈ (0 , , h ( t ) = U ( x ( t )) s k − sin ( x ( t ))1 − t ln (cid:18) k − sin ( x ( t ))1 − t (cid:19) . Note that here H i ( t ) ≡ P ( , i (2 t − where P ( , i are Jacobi polynomials. Let S = N P l =1 w l h ( t l ) .Next, we put in (54): ω ( t ) == −√ − t ln(1 − t ) ∀ t ∈ (0 , , h ( t ) = − U ( x ( t )) s k − sin ( x ( t ))1 − t . The polynomials corresponding to this value of ω can be obtained through the standard procedure oforthogonalization. Let S = N P l =1 w l h ( t l ) . Then P ≈ ( S + S )( X − X ) .In the case of i > , k = 1 and logarithmiclimb-darkening law ( K ( y ) = √ y ln y ) we represent theintegrand from (53) in the form: U ( x ( t )) cos( x ( t )) (cid:20) ln (cid:18) cos( x ( t ))1 − t (cid:19) + ln(1 − t ) (cid:21) Next, we put in (54): ω ( t ) = 1 ∀ t ∈ (0 , , h ( t ) = 2 U ( x ( t )) cos( x ( t )) ln (cid:18) cos( x ( t ))1 − t (cid:19) . Let S = N P l =1 w l h ( t l ) .Next, we put in (54): ω ( t ) == − ln(1 − t ) ∀ t ∈ (0 , , h ( t ) = − U ( x ( t )) cos( x ( t ) . Let S = N P l =1 w l h ( t l ) . Then P ≈ ( S + S )( X − X ) .In the case of i > , k < and square-root limb-darkening law ( K ( y ) = √ y γ ), we put in (54): ω ( t ) == (1 − t ) γ ∀ t ∈ (0 , , h ( t ) = U ( x ( t )) (cid:18) k − sin x ( t )1 − t (cid:19) γ . Note that here H i ( t ) ≡ P ( γ , i (2 t − where P ( γ , i areJacobi polynomials. Then P ≈ ( X − X ) N P l =1 w l h ( t l ) .In the case of i > , k = 1 and square-root limb-darkening law, we put in (54): ω ( t ) = √ − t ∀ t ∈∈ (0 , , h ( t ) = U ( x ( t )) (cid:18) cos( x ( t ))1 − t (cid:19) γ . Then P ≈ ( X − X ) N P l =1 w l h ( t l ) .Calculations show that in all cases of theapplications of the Gauss quadrature accuracy of19 significant decimal digits (corresponding to 80-bitmachine numbers) can be achieved by choosing the N to be 14. Value sets of points t i and weights w i c (cid:13) , 1– ?? universal approach to the calculation of the transit light curves corresponding to each of the considered forms of thefunction ω , can be downloaded from the Internet,along with other materials (see Conclusion). We have derived the expression for the calculationof the eclipsing binary flux and its derivatives.We considered the linear limb-darkening law, thequadratic limb-darkening law, the logarithmic limb-darkening law, and the square root limb-darkeninglaw. In general, the decrease of the flux is givenby the expression (9). In (19)-(21), ∆L correspondsto uniform brightness and its derivatives, it isexpressed in terms of easily computed piecewise-defined functions of one variable A (3) and Q (4). ∆L l corresponds to the linear limb-darkening law,given as linear combination of ∆L and ∆L (44),where ∆L with its derivatives is expressed in termsof incomplete elliptic integrals in (22)-(24). ∆L q corresponds to the quadratic limb-darkening law,given as linear combination of ∆L , ∆L and ∆L (29), where ∆L with its derivatives is given in(26)-(28). ∆L L corresponds to the logarithmic limb-darkening and ∆L Q corresponds to the square-rootlimb-darkening expressed by two- and one-parametricintegrals. Further, we described how these integralscan be found numerically by multiple application ofthe Gaussian quadrature formula. It is importantthat the nodes for this formula can be found onceand re-used for the calculations for different valuesof parameters. Also, the general integral form (15)–(17) of the flux component allows us to extend thisapproach to other limb-darkening laws.The algorithm described above was tested byAbubekerov et al. (2010, 2011); Gostev (2011) forthe interpretation of the high-precision polychromelight-curves of the binary system with exoplanetsHD 209458 Brown et al. (2001) , HD 189733Pont et al. (2007) and monochrome light-curves ofKepler-5b, Kepler-6b, Kepler-7b Koch et al. (2010);Dunham et al. (2010); Latham et al. (2010).The algorithm is implemented in ANSI Cin the form of the functions for computationof the individual component ∆L ( δ, r ) and itsderivatives. This implementation is available fromhttp://lnfm1.sai.msu.su/ ∼ ngostev/algorithm.htmlThis work was supported by the Presidentof Russian Federation grant MK-893.2012.2, RFBRgrant 12-02-31466. ACKNOWLEDGMENTS
We thank Professor Anatoly Cherepashchuk forsome helpful suggestions and useful comments thatimproved the presentation of the paper.
Список литературы
Abubekerov M.K., Gostev N.Yu., CherepashchukA.M., 2010, Astron. Rep., 54, 1105Abubekerov M.K., Gostev N.Yu., CherepashchukA.M., 2011, Astron. Rep., 55, 1051Brown T.M., Charbonneau D., Gilliland R.L. et al.,2001, ApJ., 552, 699Carlson B.C., 1994, preprint (arXiv:math.CA/9409227 v1)Dunham E.W., Borucki W.J., Koch D.G. et al., 2010,ApJ., 713, L136Gostev N.Yu., 2011, Astron. Rep., 55, 649Klinglesmith D.A., Sobieski S., 1970, AJ, 75, 175Koch D.G., Borucki W.J., Rowe J.F. et al., 2010,ApJ., 713, 131Latham D.W., Borucki W.J., Koch D.G. et al., 2010,ApJ, 713, L140Mandel K., Agol E., 2002, ApJ., 580, L171Pal A., 2008, MNRAS, 390, 281Pal A., 2012, MNRAS, 420, 1630Pont F., Gilliland R.L., Moutou C. et al., 2007 ,A&A,476, 1347Van Hamme W., 1993, AJ., 106, 2096 c (cid:13) , 1–, 1–