Free Access
Issue
A&A
Volume 584, December 2015
Article Number A8
Number of page(s) 13
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201425244
Published online 13 November 2015

© ESO, 2015

1. Introduction

There are two groups of astrophysical tasks generally solved by analysing observations of eclipsing binaries (EBs) and stars with transiting extrasolar planets. We can derive the physics of the present state of a system and its components, namely the dimensions and geometry of the system, the outer characteristics of both eclipsing bodies (e.g. their radii, shapes, masses, temperatures, limb darkening, gravity darkening/brightening, gravitational lensing, albedos, spottiness, pulsation), and the parameters of possible streams and disks. The required information is extracted by means of various, more or less sophisticated physical models of double systems (e.g. Wilson & Devinney 1971; Kallrath & Milone 1998; Hadrava 2004; Bradstreet 2005; Prša & Zwitter 2005; Prša et al. 2011; Pribulla 2012) applied to numerous and precise data of all kinds obtained by contemporary observational methods and approaches developed for this purpose.

Equally important are the description and classification of light curves (LCs) or studies of the evolution of the systems on the time scale of decades. These usually very tiny changes in light variation parameters (typically the period) may provide information on, for example, the rate of mass exchange between interacting components (e.g. Zhu et al. 2010, 2012; Mikulášek et al. 2012a), the presence and characteristics of possible invisible bodies (stars, planets, e.g. Van Hamme & Wilson 2007; Qian et al. 2005) in the system, or the degree of the mass concentration in stellar interiors (internal structure constants) if we study apsidal motion in some eccentric EBs (e.g. Kopal 1978; Claret & Gimenénez 1993; Zasche & Wolf 2013).

In period analyses of EBs we do not need all the information about the physics of the system, but only good templates of phase curves of all the data we want to analyse (Mikulášek et al. 2012a). These template phase curves (typically LCs in various photometric bands) can be obtained in principle by means of the sophisticated versions of physical EB models applied to the complete set of observational data (see in Van Hamme & Wilson 2007; Wilson & Van Hamme 2014) or to their best parts (e.g. Zasche et al. 2014; Zasche 2015). There is also the possibility of using the best observed LCs themselves (Pribulla et al. 2012). We offer another alternative using relatively simple phenomenological (mathematical) modelling of observed data variations.

Our aim is to establish a general model of LCs of eclipsing systems (ES; both eclipsing binaries and stars with transiting planets) that could fit the LCs with an accuracy of 1% of their amplitudes or better and that can be applied to the majority of observed ESs. Such a model could be used for most LC description tasks and detailed period analysis, for which other sources of phase information could also be used, especially individual eclipse timings and radial-velocity (RV) curves.

The remainder of this paper is organized as follows. In Sect. 2 we describe the general properties of the phase variation of periodically variable objects, Sect. 3 specifies the general properties of ES LC models, and Sect. 4 presents the phenomenological model of ES LCs. Section 5 is devoted to model parameters and to the search for them during analysis, in Sect. 6 we compare the results of the modelling XY Boo observations processed by means of phenomenological and physical methods, and in Sect. 7 we summarize and discuss our results of the phenomenological modelling of the ESs.

2. Periodically variable objects

Eclipsing binaries and stars with transiting exoplanets periodically change their light mainly due to the regular mutual eclipses of components (transits and occultations) and proximity effects (tidally induced ellipticity of components and reflection). The period of these ES variations corresponds to the observed orbital velocity of the system; the shape of LCs, which is dominated by eclipses and proximity effects, is thus more or less constant (for details see Sect. 3). That is why ESs are ranked into a group of periodically variable objects.

Most of the variations in periodic variables are more or less cyclic with an observed instantaneous period P(t), which is usually strictly constant or slightly variable. The period itself and its development over time cannot be observed directly, but both can be derived through analysing time series of light changes or extremum timings. For that purpose, it was useful to introduce a monotonically rising phase function ϑ(t), as a sum of the epoch E(t) and the phase in the common usage ϕ and its inversion function t(ϑ) (Mikulášek et al. 2008).

Functions ϑ(t) are tied with the instantaneous observed period P(t),P(ϑ) by the following differential equations (see Kalimeris et al. 1994; Mikulášek et al. 2008, 2012b): dϑ(t)dt=1P(t),ϑ(M0)=0;dt(ϑ)dϑ=P(ϑ);t(0)=M0,\begin{eqnarray} \label{phasefundef} \frac{\mathrm{d}\vartheta(t)}{\mathrm{d}t}=\frac{1}{P(t)},\ \vartheta(M_0)=0;\quad \frac{\mathrm{d}t(\vartheta)}{\mathrm{d}\vartheta}=P(\vartheta);\ t(0)=M_0, \end{eqnarray}(1)where M0 is the origin of counting of epochs (for ESs the time of the basic primary minimum).

Using t(ϑ) we can predict the zeroth phase time Θ(E) for the epoch E (primary minimum timing according to the ephemeris Eq. (1)), Θ(E) = Θ(ϑ = E).

2.1. Period models. Phase and time shifts

The linear period model supposes that the period P(t) of the variable object is constant. The corresponding linear phase function ϑ1(t,M0,P0), and its inversion t1(ϑ,M0,P0) are then given by ϑ1(t)=(tM0)/P0,t1(ϑ)=M0+P0ϑ.\begin{eqnarray} \label{linmod} \vartheta_1(t)=(t-M_0)/P_0,\quad t_1(\vartheta)=M_0+P_0\,\vartheta. \end{eqnarray}(2)A possible tiny modulation ΔP(t) of the basic period P0 causes a detectable phase function shift Δϑ(t) in LCs and shifts Δt(ϑ) = ΔΘ(E) in LC extrema timing: P(t)=P0+ΔP(t);ϑ(t)=ϑ1+Δϕ(t);t(ϑ)=t1+Δt(ϑ).\begin{eqnarray} P(t)=P_0+ \Delta P(t);\ \ \vartheta(t)=\vartheta_1+\Delta \varphi (t);\ \ t(\vartheta)=t_1+\Delta t(\vartheta).\label{shifts} \end{eqnarray}(3)Combining definitions in Eq. (3) with the fact that functions ϑ(t) and t(ϑ) are mutually inverse, we obtain the following relations \begin{eqnarray} &&\Delta t = -P_0\,\Delta\vartheta; \quad \Delta P(\vartheta)=\frac{\mathrm d \Delta t(\vartheta)}{\mathrm d \vartheta}=\frac{\mathrm d \hzav{t(\vartheta)-t_1(\vartheta)}}{\mathrm d \vartheta};\label{relat}\\ &&\Delta \vartheta(t)\doteq\int_{M_0}^t \hzav{\frac{\Delta P(\tau)}{P_0^2}-\frac{\zav{\Delta P(\tau)}^2}{P_0^3}}\,\mathrm{d}\tau \doteq \int_{M_0}^t \frac{\Delta P(\tau)}{P_0^2}\,\mathrm{d}\tau. \label{deltafi} \end{eqnarray}The last approximation in the Eq. (5) is valid for all known EBs, including SV Cen with a record-breaking decrease in its orbital period by /P0 = −2.36(5) × 10-5 yr-1 (Drechsel et al. 1982).

The time development of phase function shifts Δϑ(t) can be derived by analysing LCs, whilst ΔΘ(E) can be revealed using standard O–C diagrams constructed by means of extrema timings.

2.2. Basic period models. Finding O–C shifts

If the time development of the period P(t) is continuous and smooth, we can express it in the form of the Taylor polynomials with the centre at t = M0, ϑ1 = (tM0) /P0P(t)=P0+P00ϑ1+P02¨P0ϑ122!+···+P0kdkP0dtkϑ1kk!···\begin{eqnarray} \textstyle P(t)=P_0+P_0\dot{P}_0\,\vartheta_1 +P_0^2\ddot{P_0}\frac{\vartheta_1^2}{2!}+\cdots+ P_0^k \frac{\mathrm d^k\!P_0}{\mathrm dt^k}\frac{\vartheta_1^k}{k!}\cdots \label{perlaurin} \end{eqnarray}(6)Using Eqs. (2)–(4), and a simplified version of Eq. (5) we obtain ϑtM0+P0ϑ+P00ϑ22!+P02¨P0ϑ33!+···+P0kdkP0dtkϑk+1(k+1)!···Θ\begin{eqnarray} \vartheta&&\textstyle\simeq\vartheta_1 -\dot{P_0}\,\frac{\vartheta_1^2}{2!}- P_0\ddot{P}_0\,\frac{\vartheta_1^3}{3!}-\cdots - P_0^{k-1} \frac{\mathrm d^k\!P_0}{\mathrm dt^k}\frac{\vartheta_1^{k+1}}{(k+1)!}\cdots \label{thetalaurin}\\ \hspace*{-18mm}&& t\textstyle\simeq M_0\!+\!P_0\,\vartheta+P_0\dot{P}_0\frac{\vartheta^2}{2!}\!+\! P_0^2\ddot{P}_0\frac{\vartheta^3}{3!}\!+\!\cdots\!+\!P_0^k \frac{\mathrm d^k\!P_0}{\mathrm dt^k}\frac{\vartheta^{k+1}}{(k+1)!}\cdots \nonumber \\ \mathit\Theta&&\textstyle\simeq M_0\!+\!P_0\,E+P_0\dot{P}_0\frac{E^2}{2!}\!+\! P_0^2\ddot{P}_0\frac{E^3}{3!}\!+\!\cdots\!+\!P_0^k \frac{\mathrm d^k\!P_0}{\mathrm dt^k}\frac{E^{k+1}}{(k+1)!}\cdots \label{Thetalaurin} \end{eqnarray}Similarly, we can establish other arbitrarily complex period models of ϑ(t) determined by a set of free parameters, including cyclic period modulation of the light time effect (LiTE) caused by another body in the system (Mikulášek et al. 2011c; Liška et al. 2015) or the apsidal motion.

The real shape of the phase curve can also be approximated using so-called O–C time shifts of the observed phase curves versus the predicted LC derived by the period model with fixed parameters (typically P0,M0) expressed in time units (usually in days). We divide the whole time interval covered by observations into nOC appropriate time intervals (typically nights or seasons). The phase function ϑ(t,r) during a certain rth time interval with the (O–C)r time shift is then given by the formula (Mikulášek et al. 2011b) ϑ(t,r,ϑc,{OC})=ϑc(t)s=1nOCηrs(OC)sP,\begin{eqnarray} \label{oc} \vartheta(t,r,\vartheta_{\mathrm c},\{{\rm O{-}C}\})=\vartheta_{\mathrm c}(t) - \sum_{s\,=\,1}^{n_{\mathrm{OC}}} \,\eta_{\rm rs}\,\frac{({\rm O{-}C})_s}{P}, \end{eqnarray}(9)where ϑc(ti) is a predicted phase function calculated by an appropriate period model at the time ti, O–C is a set of all nOC values of the found (O–C)r time shifts versus this model. The symbol ηrs represents a discrete function (a table or a matrix), which for each individual ith observation from our data assigns either 1, if its order number of the interval ri = s, or 0 if ris.

The dependence of O–C values on the epoch serves as a common O–C diagram, the basic tool for the period analysis. Using the found individual O–C values and their uncertainties, we can calculate a set of so-called virtual minima timings (Mikulášek et al. 2011b,a, or Sect. 6.4.2 in this paper). Virtual minima timings can be combined with others derived, for instance, by other techniques (e.g. Brát et al. 2012; Mandel & Agol 2002; Kwee & van Woerden 1956; Mikulášek et al. 2006, 2014; Zasche et al. 2014; Zasche 2015).

3. General properties of an eclipsing system light curve model

3.1. Instrumental term of an observed light curve

The observed LC (or its segment) of a chosen ES in a particular colour of an effective wavelength λeff is defined by a time series { tri,yri } obtained during an observational interval r (observing night, part of it, or a season). The rth subset of observational data can generally be modelled by the function Yr(t,λeff), expressed in magnitudes1: Yr(t,λeff)=Y0r(t,λeff)+F(ϑ,λeff),\begin{eqnarray} Y_r(t,\lef)=Y_{0r}(t,\lef)+F(\vartheta,\lef), \label{general} \end{eqnarray}(10)where Y0r(t,λeff) is an additive term removing instrumental trends and unifying the observed magnitudes or magnitude differences { yi }, whilst F(ϑ,λeff) is an intrinsic LC model function of the phase function ϑ (see Sect. 2.1) and the effective wavelength λeff, free of instrumental and observational shifts and trends. The function F(ϑ,λeff) is normalized so that its mean value without eclipses is equal to zero. The function Y0r(t,λeff) can be approximated by a linear combination of gr dimensionless functions of time t, Ξj (t), with magnitude-like coefficients m0rj(λeff): Y0r(t,λeff)j=0grm0rj(λeff)Ξj(t),\begin{eqnarray} Y_{0r}(t,\lef)\simeq \sum_{j\,=\,0}^{g_r} \,m_{0rj}\,(\lef)\ \mathit{\Xi}_{j}\,(t),\label{instrument} \end{eqnarray}(11)where Ξj (t) can be, for example, normalized polynomials Ξrj(t)=[(ttr¯)/σ(tr)]j\hbox{$\mathit{\Xi}_{rj}\,(t)=[(t-\bar{t_r})/\sigma(t_r)]^j$} (σ(tr) is the weighted variance of observational times in the segment r) or Legendre polynomials or special quasi orthogonal functions combining polynomials with harmonic functions (see in Mikulášek & Gráf 2005). The set of coefficients { m0rj (λeff) } is found together with parameters describing the model function F(ϑ,λeff).

3.2. Bases and limitations of the phenomenological model of eclipsing system light curves

Light curves of ESs are nearly periodic functions, and it would be natural to express them in the form of Fourier series (see e.g. Rucinski 1973; Kallrath & Milone 1998; Selam 2004; Nedoroščík et al. 2012; Andronov 2012). This concept proves its worth in many types of extrinsic periodically variable stars, especially in the case of rotating variables with photometric spots on their surfaces (North 1984; Mikulášek et al. 2007) or non-eclipsing (e.g. purely elliptical) double stars, where we manage with harmonic polynomials of a low degree (e.g. Kallrath & Milone 1998).

However, it is generally known that the presence of eclipses in LCs asks for the use of harmonic polynomials of a rather high degree if a proper fit of the observed LC is required. Their usage is badly influenced by departures from the ideal equidistant distribution of observations according to the orbital phase. Even small phase gaps are then filled with unreal LC artefacts.

It seems that it is better to use properly selected phenomenological models of LCs of ESs described by a few parameters. Several more or less successful attempts on how to model these LCs have been proposed, for instance, by Tsesevich (1971) and Kholopov (1981) and see the reviews and references in Andronov (2012) and Chrastina et al. (2013), Table 2, and Figs. 1, 9.

Theory is able to explain the observed periodic light variations of an ES in general and in detail as the result of alternating mutual eclipses of components of the system, non-isotropic radiation of orbiting components, caused by their close proximity, and by unevenly distributed photometric spots on rotating components (e.g. Bradstreet 2005; Hilditch 2001; Pribulla et al. 2012). The periodicity of light changes caused by mutual eclipses and proximity effects, as well as variations connected with the rotation of synchronously rotating spotted components, is dictated by their observed orbital motion. These changes are periodic with an instantaneous period2P(t), the shapes of LCs remain more or less constant for decades.

For purely geometrical reasons, the prevailing majority of EBs ranks among relatively close, hence tidally interacting systems, where the processes of synchronizing the components’ rotation and orbit circularization are strong and effective (Zahn 1992; Goldman & Mazeh 1991, and references therein). That is why the rotations of components are usually synchronous, and more than 80% of their orbits are fairly circular (see CALEB, Bradstreet 2005). In a further introductory text, we will concentrate mainly on systems with more or less circular orbits3 and constant LCs4.

The intrinsic one-colour LC function (expressed in magnitudes) of an ES is a periodic function F(ϑ,λeff) that can be approximated as the sum of three more or less independent terms (see also in Andronov 2012): F(ϑ,λeff)=Fe(ϑ,λeff)+Fp(ϑ,λeff)+Fc(ϑ,λeff),\begin{eqnarray} \label{EBbas} F(\vartheta,\lef)=F_{\mathrm e}(\vartheta,\lef)+F_{\mathrm p}(\vartheta,\lef)+F_{\mathrm c}(\vartheta,\lef), \end{eqnarray}(12)where Fe(ϑ,λeff) describes a LC of eclipses (Fe ≡ 0 outside of eclipses), whilst Fp(ϑ,λeff) and Fc(ϑ,λeff) express contributions from proximity and O’Connell (1979) effects without eclipses (Fp=Fc=0\hbox{$\overline{F_{\mathrm p}}=\overline{F_{\mathrm c}}=0$}). The mathematical models are formulated and discussed in Sect. 4.1.2.

4. Phenomenological model of eclipsing system light curves

4.1. Model of one-colour light curves

4.1.1. Eclipses

The essential feature of all ES LCs are two nearly symmetrical depressions caused by mutual eclipses of synchronously rotating stellar or planetary components. The profiles of both minima are complex functions determined primarily by the geometry of the system and the relative brightness of components in a given spectral region centred at the effective wavelength λeff. The contribution of eclipses Fe(ϑ,λeff) to an ES light curve can be approximated by a sum of two special periodic functions of phase function ϑ. In the case of circular orbits, eclipses are exactly symmetrical around their centres at phases ϕ01 and ϕ02. If we put the origin of the phase function M0 at the time of the primary minimum, then ϕ01 = 0,ϕ02 = 0.5.

The model function was selected so that it describes as aptly as possible those parts of LCs that are in the vicinity of their inflex points, where their slopes are maximum. The functions are parameterized by their widths D1,D2, eclipse LC kurtosis coefficients Γ1,Γ2, dimensionless correcting factors C1,C2, and central depths A1(λeff),A2(λeff): Fe(ϑ,λeff)=k=1neAk(1+Ckϕk2Dk2)1{1exp[1cosh(ϕkDk)]}Γk,ϕk=ϑ0.5(k1)round[ϑ0.5(k1)],\begin{eqnarray} \label{modelecl} &&\displaystyle F_{\!\mathrm{e}}(\vartheta,\!\lef)=\! \sum_{k\,=\,1}^{n_{\mathrm{e}}}\!A_{k}\zav{1\!+\!C_{\!k}\!\frac{\varphi_k^2}{D_k^2}} \left\{1\!-\!\left\{1\!-\!\exp\!\hzav{1\!-\!\cosh\!\zav{\frac{\varphi_k}{D_k}}} \right\}^{\mathit{\Gamma\!_k}}\right\},\nonumber\\ &&\varphi_k=\vartheta-0.5\,(k-1) -\mathrm{round}\hzav{\vartheta-0.5\,(k-1)}, \end{eqnarray}(13)where the summation is over the number of eclipses during one cycle, ne: ne = 2 or ne = 1 (the common situation for exoplanet transits). Each eclipse in a given colour is thus described by only four parameters – its depth A, width D, kurtosis Γ, and the correcting parameter C.

thumbnail Fig. 1

Examples of elementary model LC Ψ = Fe/Ak for C = 0; see Eq. (13) allowed by our phenomenological EB models. Parameter d describes the duration of the eclipse, and Γ is a parameter expressing the kurtosis of the LC.

In the case of EBs with two minima in a cycle (ne = 2), we need eight parameters, but sometimes the number of needed parameters can be smaller. Inspecting the parameters D,Γ, and C for both eclipses of many EBs, we have concluded that they are as a rule nearly the same: especially D1D2,Γ1Γ2, and C1C2. We therefore usually need only five monochromatic parameters (A1,A2,D,Γ,C). The parameter C is mostly comparable to its uncertainty, so we can neglect it entirely. Then we need just four parameters! On the other hand, in EBs with totalities we see that the bottoms of their occultations are flat whilst transits are convex. It can be described by introducing of different parameters C1,C2 (see the case of EK Com in Table 1, Fig. 2).

The LCs of the exoplanet transits (ne = 1) need only four parameters (A,D,Γ,C – see Fig. 3); in cases of very precise measurements, we add another dimensionless parameter K (Mikulášek et al. 2015): Fe=A(1+Cϕ2D2+Kϕ4D4){1{1exp[1cosh(ϕD)]}Γ}.\begin{eqnarray} \label{modeltrs} \displaystyle F_{\!\mathrm{e}}=A\zav{1\!+\!C\frac{\varphi^2}{D^2}\!+ \!K\frac{\varphi^4}{D^4}} \left\{1\!-\!\left\{1\!-\!\exp\hzav{1\!-\!\cosh\zav{\frac{\varphi}{D}}} \right\}^{\mathit{\Gamma}}\right\}. \end{eqnarray}(14)Testing several dozens of LCs of various types of ESs, we found that the standard deviation of the fit is typically well below one per cent. The only minor inconvenience is the existence of a spike (a jump in derivatives) in the mid-eclipses for LCs with Γ< 1 (see Fig. 1).

thumbnail Fig. 2

Fit of the V curve of the overcontact spotted EB EK Com (the secondary minimum is a transit). The LC of the close binary is affected by O’Connell and proximity effects. The synthetic model curve is depicted by a solid line, and the phenomenological LC fit is depicted by the dashed line. For more information see Sect. 4.3 and Table 1.

thumbnail Fig. 3

Simultaneous fit to 15 transits of extrasolar planet TrES-3b corrected for trends (see Eq. (10)). The fit parameters according to Eq. (14): A = 0.02725(19) mag, D = 0.0117(7), C = −0.17(7), and Γ = 1.58(12), the parameter K was not introduced. The black line is the fit, circles are the normal points, grey circles indicate individual measurements with the area proportional to their weights.

4.1.2. Proximity effects. O’Connell effect

Light variations of EBs caused by eclipses are usually modified by the asphericity of the components, effects of gravity darkening/brightening, and mutual irradiation. All these proximity effects are the manifestation of the interaction between the components acting, namely, in close systems. Contrary to eclipses, the proximity effects modify LCs permanently in each phase.

Light curves of some EBs are influenced by the O’Connell effect that results in the asymmetry of some LCs of close EBs, manifesting itself as the difference in light maxima between eclipses. The standard explanation for this is the presence of one or more cool or hot spots on the surface of one of the synchronously rotating components or by asymmetrically distributed circumstellar material in the system (e.g. Wilsey & Beaky 2009; Pribulla et al. 2012, and citations therein). The amount and the sign of the O’Connell effect vary with time (Beaky & Koju 2012). The effect is also wavelength dependent – in blue it is usually stronger, but this is not a rule (Pribulla et al. 2012).

The contribution of proximity effects Fp(ϑ) should be an even function symmetric with the phases 0.0 and 0.5, consequently they can be satisfactorily expressed as a linear combination of np elementary cosine functions cos(2 πϑ), cos(4 πϑ), cos(6 πϑ), etc. The even terms are the consequence of the ellipticity of tidally interacting components, whilst the odd terms result from the differences between the near and far sides of components. As a rule we can limit ourselves only to the first two or three terms in the Fp (Russell & Merrill 1952; Kallrath & Milone 1998). The O’Connell effect contribution Fc(ϑ) can be modelled well by a simple sinusoid (Davidge & Milone 1984; Wilsey & Beaky 2009): Fp=k=ne+1np+neAkcos[2π(kne)ϑ],Fc=k=np+ne+1nc+np+neAksin(2πϑ),\begin{eqnarray} \label{proxcon} F\!_{\rm p} =\!\sum_{k\,=\,n_{\mathrm e}+1}^{n_{\mathrm p}+ n_{\mathrm e}}\! A_k\cos\hzav{2\pi(k\!-n_{\mathrm e})\vartheta},\ F\!_{\rm c}=\!\sum_{k\,=\,n_{\mathrm p}+n_{\mathrm e}+1}^{n_{\mathrm c}+ n_{\mathrm p}+n_{\mathrm e}}\! A_k\sin(2\pi\vartheta), \end{eqnarray}(15)where np is the number of terms in Fp(ϑ): np = 0,1,2,3,··· , nc = 0, if the O’Connell asymmetry is not present5, otherwise nc = 1.

thumbnail Fig. 4

V LC of a famous EB consisting of a very hot nucleus of a planetary nebula and solar type star. For the fit of the LC strongly affected by proximity effects, we need only 8 parameters phenomenological. For details see Sect. 4.3 and Table 1.

thumbnail Fig. 5

Fit of the primary minimum of the detached EB AR Aur in V and B, represented by normal points (circles), by physical (grey solid lines), and phenomenological (dashed lines) LCs (see Table 1).

The V LC of the close EB EK Com (see Fig. 2) with apparent O’Connell effect is described by nine parameters: A1,2,3,4,5,D,C1,2,  ne = 2, np = 2,nc = 1. The uncommon V LC of EB V477 Lyr (see Fig. 4) is determined by eight parameters: A1,2,3,4,5,D,C,Γ, ne = 2, np = 3,nc = 0.

4.2. The model of multicolour light curves

The parameters of the model functions defined above, especially the amplitudes Ak, and parameters C1,2, D1,2 and Γ1,2 are generally functions of the wavelength λ.

In principle, we can use the one-colour models formulated in Sect. 4.1 separately, assuming that all of the parameters are wavelength dependent. Fortunately, it follows from our experience with modelling of LCs of hundreds of real systems and their physical models that it is not necessary to take response curves of different photometric passbands into account; we manage with their effective wavelengths alone. It enables the association of photometric colours with different transparency widths and equal or close effective wavelengths. (Typically we are allowed to combine measurements done in V and y.)

On top of that, the dependencies of model parameters on the effective wavelength λeff are typically smooth and mostly monotonic, so we can approximate them by low-order polynomials of the dimensionless parameter Λ; Λ = λ0/λeff−1, where λ0 is an arbitrarily selected central wavelength of the data set (Mikulášek et al. 2015, see also Fig. 7): Ak=j=1gAkakjΛj1,Cl=j=1gClcljΛj1,Dl=j=1gDldljΛj1,Γl=j=1gΓlγljΛj1,\begin{eqnarray} \label{multi} &&\textstyle A_{k}=\sum_{j\,=\,1}^{g_{ {A\,k}}}a_{kj}\,\mathit{\Lambda}^{j-1},\quad C_l=\sum_{j\,=\,1}^{g_{C\,l}}c_{lj}\,\mathit{\Lambda}^{j-1},\\ &&\textstyle D_l=\sum_{j\,=\,1}^{g_{D\,l}}d_{lj}\,\mathit{\Lambda}^{j-1},\quad \mathit{\Gamma}_{l}=\sum_{j\,=\,1}^{g_{\mathit{\Gamma} l}}\gamma_{lj}\,\mathit{\Lambda}^{j-1}, \nonumber \end{eqnarray}(16)where gAk,gCl,gDl,gΓl, l = 1,2 or l = 1 are the numbers of degrees of freedom of the corresponding parameters of the model. The standard set of the one-colour LC model parameters of EBs (see Sect. 4.1): { Ak,Cl,Dl,Γl } can be considered as the special case of the multicolour decomposition Eq. (16) for λ0 = λeff = 0, j = 1: { ak1,cl1,dl1,γl1 }.

The set of relations in Eq. (16) enables the calculations of all parameters { Ak(Λ),Cl(Λ),Dl(Λ)l(Λ) } needed for calculating the model of a LC in any photometric band characterized by the parameter Λ. Figure 6 shows the fit of the extreme ten-colour photometry (350–980 nm) of an exoplanet transit (Knutson et al. 2007) – here we need nine parameters, namely A = a1 + a2Λ, D = d1 + d2Λ, C = c1 + c2Λ + c3Λ2, Γ = γ1 + γ2Λ, where Λ = 550 /λeff−1 (Mikulášek et al. 2015).

thumbnail Fig. 6

Ten-colour transit LCs for HD 209458b taken from Knutson et al. (2007). Effective wavelengths of individual colours are in nm. For the description of all 10 LCs we only need 9 parameters (courtesy of Mikulášek et al. 2015).

The fit of BVRcIc proper LCs of AV Del (see Fig. 8) needs only nine parameters, namely D1 = D2 = D, Γ1 = Γ2 = Γ, C1 = C2 = C, a11,a12,a21,a22,a41, and a42, ne=2, np = 2,nc = 0, A3 ≡ 0. See also Table 1.

Table 1

Parameters describing LCs of several EBs.

4.3. Brief description of selected eclipsing systems

TrES-3b is an extrasolar planet orbiting the star GSC  03089-00929 with a period of 31 hours. It belongs to the hot Jupiters that are undergoing orbital decay due to tidal effects. For the LC inspection (Fig. 3), the R photometry of 15 transits containing 2820 individual measurements in total (courtesy Vaňko et al. 2013) was used. The parameters of the LC fit of TrES-3b are in the legend of Fig. 3.

HD 209458 was the first star found to have a transiting planet (Charbonneau et al. 2000; Henry et al. 2000), and it remains the second brightest star known to have a transiting planetary companion. Knutson et al. (2007) obtained 1066 spectra over four distinct transits with the STIS spectrometer on HST allowing us to synthesize LCs in ten spectrophotometric bandpasses in 290–1030 nm (see also Sect. 4.2, Fig. 6 and Mikulášek et al. 2015).

Table 1 contains the parameters of the phenomenological fit of LCs and some other information on the following selected EBs:

  • AR Aurigae, a prototype of a detached EB(O’Connell 1979).

  • EK Comae, an overcontact, spotted EB with a short orbital period (Samec et al. 1996).

  • AV Delphinis, a “cool Algol” consisting of a F type primary on the main sequence and a K subgiant filling its Roche lobe (Mader et al. 2005).

  • V477 Lyrae, an unusual, detached EB consisting of a very hot and luminous nucleus of the planetary nebula as a primary component and a solar type star as secondary (Pollacco & Bell 1994).

Table 1 shows that the fit of LCs of the above-mentioned stars by our models, quantified by the ratio ρ is nearly the same or better than in the case of the fit of BM3 (CALEB) physical model (Bradstreet 2005).

5. Phenomenological model solution

5.1. Finding model parameters and their uncertainties

The procedure for finding model parameters is based on the simultaneous mathematical processing of all relevant photometric data consisting of individual photometric observations, including barycentric Julian date of the ith measurement ti, the measured magnitude or magnitude difference corrected for possible trend(s) during nights or seasons yi, and the estimate of its uncertainty σi. Furthermore, we should know the effective wavelengths of the photometric filter used, λeffi or Λi, and submission of an individual observation to one of the observational subsets ri (see Sect. 3.1).

thumbnail Fig. 7

Dependence of the amplitude A of transit of HD 209458b (see Fig. 6) on Λ = 550 /λeff−1 proves the validity of the approximation Eq. (16) for that case.

thumbnail Fig. 8

Fit of the BVRcIc LCs of a semidetached short periodic “cool Algol” AV Del. Dark dots are observations, and grey lines indicate the 13-parameter phenomenological fit. For details see Table 1 and Sect. 4.3.

For simplicity we assume that the shapes of LCs are constant, and the variability of an object is described by the unique model function Eq. (10), consisting of the instrumental term Yr0(t,λeff) (Eq. (11)) and the intrinsic phenomenological ES model LC function F(ϑ,λeff), specified in Sect. 4. The phase function ϑ(t,P0,M0,0, or { OC }, etc.) is a function of time and some free parameters of a variety of period models offered in Sect. 2. The result of the solution – the full set of g free parameters of the complete model including the estimate of the parameter uncertainties – was evaluated simultaneously using the non-linear least square method by minimizing the quantity χ2 by the standard technique (using tried Newton-Raphson method of the non-linear equation solution) described well in Press et al. (2002), Hayashi (2000), Hartkopf et al. (1989), and Mikulášek et al. (2011a), among others. With a good initial estimate of the parameter vectors, the iterations converge fairly quickly.

All estimates of uncertainties of model parameters were computed using formulae that take that our models fit phase curves of EBs with uneven accuracy. Since the models are not orthogonal, uncertainties of the functions of model coefficients (typically the fits of LCs or minima times) should be computed by the general law of uncertainty propagation assuming also correlations amongst individual coefficients (see e.g. Bevington & Robinson 2003; Mikulášek & Zejda 2013). It is advisable to orthogonalize the models at least in the ephemeris parameters, in accordance with what we did in Mikulášek (2007) and Mikulášek et al. (2008).

5.2. Selection of an optimal model of light curves

The models of LCs should be tailored to the studied object(s), available data, and the purpose of fitting the data. Especially the number of used free parameters should be as few as possible, but without any serious influence on the accuracy and reliability of the results. It follows from our experience that it pays off to adhere to the following general principles:

  • 1.

    It is always advantageous to process all available datasimultaneously and not divided into parts. It is also valid fordetermining mid-eclipse times of individual eclipses where weshould use Eq. (9).

  • 2.

    We have to pay attention to the right weighting of entered data because it is required by the used χ2 regression. If individual uncertainties of original data are not known in advance, we have to estimate them iteratively from the scatter of residuals { Δyi } for appropriately created data subsets.

  • 3.

    Fixing the LC model parameters, if they are known from previous analyses, is also advised. In addition, neglecting all model terms whose amplitudes are less than or comparable to their uncertainties is recommended.

  • 4.

    The number of free parameters can also be reduced by using simpler model functions of eclipses listed in Table 2. Suggested model functions are compared in Fig. 9.

  • 5.

    The complexity of the selected phenomenological model should correspond to the purpose of the modelling. Several tasks (e.g. the basic classification of LCs) allow for the use of simple, unified models with basic parametric outfit.

Unfortunately, the effort of using the optimum phenomenological models considerably encumbers automation of the computational process. The diversity of real LCs of particular ESs requires that they be solved individually.

6. Period analysis of XY Bootis in 1955–2009

The phenomenological modelling of periodic variable stars described above has been developed first and foremost for the sake of the simultaneous period analysis of all available data containing phase information. In the case of ESs, such data are times of light minima, LCs or their segments, and RV curves. The majority of studies of ES period changes were based on analysing O–C diagrams constructed using the times of light minima, where these timings were determined individually directly by observers who did not take LCs of the star observed before into account. Several other studies were based only on analysing LCs obtained for several years. Measurements of radial velocities were used, namely, for solving the geometry of the system.

The simultaneous analysis of data of various kind is possible using an extended version of the contemporary sophisticated codes for the physical solution of ESs (see e.g. Hadrava 2004; Van Hamme & Wilson 2007; Wilson & Van Hamme 2014). The concept of EB period analyses without using O–C diagrams and the physical solution of the system were outlined, among others, in Mikulášek et al. (2011b, 2012a, 2013a), and studies of period changes of individual EBs (Zhu et al. 2010, 2012; Mikulášek et al. 2013b). The results of the preliminary report on the ephemeris of AR Aur – a star with the light time effect (Mikulášek et al. 2011c) – were referred to and extensively discussed by Wilson & Van Hamme (2014), who used AR Aur as an example of systems with a third companion.

Wilson & Van Hamme (2014) also studied the W UMa type EB XY Bootis as a prototype of a close double star with a nearly steady period change. The study clearly proves the advantage of simultaneous processing of all relevant data. That is why we want to compare the results of this sophisticated study describing the period change of XY Boo during the time interval 1955–2009 with our results based on strictly the same observational data.

Table 2

List of alternatives for several phenomenological model functions of eclipses.

6.1. XY Bootis

The eclipsing binary XY Bootis (=BD+ 20°2874 = HIP 67431; Vmax = 10.3 mag; Sp. F5V) was discovered as a variable star by Hoffmeister (1935). Tsesevich (1950, 1954) observed the star visually and classified it as an EB of W UMa-type.

The first photoelectric observations were obtained by Hinderer (1960). Wood (1965) reanalysed Hinderer’s data and found the true period of the star P=0.d37054\hbox{$P = 0\fd37054$}. Binnendijk (1971) observed six minima of XY Boo in BV, improved the ephemeris (M0=2440389.7321,P=0.d37055)\hbox{$(M_0\,=\,2440389.7321,\ P\,=\,0\fd37055)$}, and revealed a remarkable increase in the period. Winkler (1977) also observed XY Booin 1976 in BV and derived three other times of minima (see Table 6). Awadalla & Yamasaki (1984) gave two new times of light minima and confirmed a period increase. The history of the XY Boo investigation from its discovery to 1998 is described more extensively in Molík & Wolf (1998). McLean & Hilditch (1983) spectroscopically measured the RV of both components and estimated the mass ratio q = 0.16 ± 0.04.

The first detailed period study was done by Molík & Wolf (1998), who used 43 moments of minimum spanning the interval 1944–98 and calculated the quadratic ephemeris. The found record-breaking period increase = 1.67(5) × 10-9 = 5.3 s per century was explained by mass transfer of 1.34 × 10-7M yr-1 from the secondary to the primary component. The results were confirmed and improved by Yang et al. (2005), who found = 1.711(6) × 10-9 = 5.4 s per century on the basis of standard O–C analysis of 54 minima times.

In their calculations of the quadratic ephemeris of XY Boo, Wilson & Van Hamme (2014) combined phase information hidden in BV LCs obtained by Binnendijk (1971), Winkler (1977), and Awadalla & Yamasaki (1984), further RV curves of McLean & Hilditch (1983), all times of minima listed in Yang et al. (2005), and other 33 minima timings collected in Table 1 of their article covering the interval 2005–09. Using all of those data, they found the mean period increase = 1.6348(8) × 10-9 = 5.159(24) s per century. They also noticed some oscillation from JD 2 448 000 to 2 455 000 with no indication of periodicity. We have now collected many other observations that prove the complexity of the XY Boo period variations. Nevertheless, for comparing the effectiveness of our method, we used exclusively those data used by Wilson & Van Hamme (2014) in the following short study.

6.2. Phenomenological model of XY Boo variability

thumbnail Fig. 9

Difference in mag between the physical model primary minimum of AR Aur and various alternatives of phenomenological basic function listed in Table 2.

Similar to Wilson & Van Hamme (2014), we assumed that the instantaneous period P(t) of XY Boo is lengthening with the constant rate (t) = const. Then the phase function ϑ(t) and the prediction of primary minimum times Θ(E) can be approximated (according to Eqs. (7), (8)) by simple relations: ϑ=ϑ1ϑ122;Θmin=M0+P0E+12P0E2,\begin{eqnarray} \vartheta=\vartheta_1-\dot{P}\,\frac{\vartheta_1^2}{2};\quad \label{xytheta} \mathit{\Theta}_{\mathrm{min}}=M_0+P_0\,E+\textstyle{\frac{1}{2}}P_0\,\dot{P}\,E^2, \end{eqnarray}(17)where ϑ1 = (tM0) /P0, M0 is the JD timing of the basic primary minimum – the origin of counting of epochs E, P0 = P(t = M0) is the instantaneous period at the time t = M0. Integer doubling of epoch E is done (even for times of primary minima, odd for times of secondary minima)6.

Light curves of XY Boo with almost equally deep minima (see Fig. 10) agree with W UMa classification. We found that all of them can be well-fitted by the following phenomenological model with nine free parameters: a11,a12,a21,a22,A3,A4,D1,D2, and Γ, FLC(ϑ)=k=12(ak1+ak2Λ)1{1exp[1cosh(ϕkDk)]}Γϕk=[ϑ(k1)/2]round[ϑ(k1)/2],Λ=550λeff1.\begin{eqnarray} &&\displaystyle F_{\mathrm{LC}}(\vartheta)=\sum_{k\,=\,1}^2 (a_{k1}+a_{k2}\mathit{\Lambda})\left\{1\!-\!\left\{1\!-\!\exp\hzav{1\!- \!\cosh\zav{\frac{\varphi_k}{D_k}}} \right\}^{\mathit{\Gamma}}\right\} \nonumber \\ &&\hspace*{1.2cm}+\,A_3\cos(4\,\pi\,\vartheta)+A_4\cos(6\,\pi\,\vartheta);\label{xyLCeq}\\ &&\displaystyle \varphi_k=\hzav{\vartheta-(k-1)/2} - \mathrm{round}\hzav{\vartheta-(k-1)/2},\quad \mathit{\Lambda}=\frac{550}{\lambda_{\mathrm{eff}}}-1. \nonumber \end{eqnarray}(18)The results achieved using this phenomenological model of XY Boo BV LCs are denoted as “Model I”.

thumbnail Fig. 10

BV curves of XY Boo (see the list in Table 5). The phase is plotted according to our quadratic ephemeris (see Table 4). The areas of points are proportional to their weights. Open squares are normal points, each of which represents the mean of about 50 measurements. Grey lines are fits by our phenomenological model (see Eq. (18) and Table 3). ΔB and ΔV display residuals of BV magnitudes and normal points from the LC phenomenological model. The scale of residuals is two times larger than the measure of BV LCs.

Other LC templates were BV synthetic LCs computed by the freely accessible W-D 2013 code with physical parameters of XY Boo published in Wilson & Van Hamme (2014). The results obtained using these BV LC templates are denoted as “Model II”. Figure 11 displays that the difference between the found phenomenological and synthetic BV model LCs are only marginal; numerically it is also documented in Table 3.

Table 3

Nine parameters of phenomenological (this paper) and synthetic LCs (Wilson & Van Hamme 2014) models.

For the description of shapes of radial velocity (RV) phase curves, we used only a simple sinusoidal model (see Eq. (19)), neglecting the Rossiter-McLaughlin effect (McLaughlin 1924) RVk(ϑ)=Vγ+A5(k1q)sin(2πϑ),\begin{eqnarray} {\it RV}_k(\vartheta)=V_{\gamma}+A_5 \zav{k-1-q}\,\sin(2\,\pi\,\vartheta),\label{xyrveq} \end{eqnarray}(19)where k = 1 refers to the first component, whilst k = 2 to the secondary one, Vγ is the so-called γ velocity, A5 is the amplitude of the difference in the RV between the components, all in km s-1, and q = m2/m1 is the ratio of component masses.

This approximation is fully sufficient for RV observations of McLean & Hilditch (1983) that cover only about 15% of the phase curve. Results are in the bottom part of Table 4.

thumbnail Fig. 11

Models of B and V LCs of XY Boo. Grey lines are BV model LCs simulated by the Wilson & Van Hamme (2014) physical model; dotted lines correspond to our phenomenological model (see Eq. (18) and Table 4); dashed lines are hypothetical BV LCs of XY Boo without mutual eclipses.

Table 4

Comparison of common parameters derived by Wilson & Van Hamme (2014) and us.

6.3. Used data and their weights. Models I and II

Our period analyses of XY Boo was based on 1770 individual data points (see the list in Table 5) acquired by three different techniques divided into nine groups with various scatters σI and σII with respect to phenomenological (I) and synthetic (II) model phase curves.

Scatter σI used to be the same or a bit smaller than the scatter σII compared to computed synthetic phase curves (Wilson & Van Hamme 2014). Scatters of RV measurements – σRVI and σRVII 28 and 30 km s-1 – seems to be a good estimate of their values. The mean uncertainty of eclipse times7σTmin with respect to the model Eq. (17) is 0.d0043=6.2\hbox{$0\fd0043=6.2$} min. It should be mentioned that the inner uncertainty of individual times of minima is typically one order better. The uncertainty σTmin is caused mainly by erratic fluctuations in the rate of the orbital period and the shape of LCs. The character of these changes is displayed in Fig. 12.

Scatter σI, II for an individual subset of data (see Table 5) served as the basis for weighting of each used measurement and interconnection of data of various nature. It enabled applying the χ2 regression computation outlined in Sect. 5. The computation represents the solution of 21 (Case I) or 12 (Case II) non-linear equations of 21 or 12 free parameters. The difference in the number of parameters needed is a result of the fact that in Case II we adopted template LCs from the paper by Wilson & Van Hamme (2014). This means that we fixed all nine parameters that we needed for determining BV template LCs.

Table 5

List of used observational data including their source, number, specification, and scatter with respect to phenomenological model phase curves (Model I) and synthetic phase curve (Model II).

In the first part of Solutions I and II, we obtained the parameters of the unified XY Boo ephemerides: M0, P0, and (see Table 4), which can be compared with the same parameters listed in Table 3 of Wilson & Van Hamme (2014) paper. Computed parameters of the B, V, and RV phase curves according to these ephemerides are in Tables 3 and 4, respectively. In the second part of the solutions, we derived 12 virtual eclipse times and their uncertainties by means of the model Eq. (9) from BV photometry and RV data – see Table 6, Figs. 12, 13.

Table 6

Comparison of the eclipse timings of XY Boo.

thumbnail Fig. 12

XY Boo timing residuals from the quadratic ephemeris (with a term). Published minima timings are denoted as squares, timings derived from LCs are circles whose areas are proportional to their weight.

thumbnail Fig. 13

XY Boo timing residuals from the linear ephemeris (without a term). Published minima timings are denoted as squares, timings derived from LCs are circles whose areas are proportional to their weight, and a timing derived from RV curves is labelled as a diamond.

6.4. Discussion of results. Comparison with W&VH solution

6.4.1. General remarks

We discuss here the results of three methods for period analysis that enables the processing of various type of data containing phase information. For their consistent comparison, we used exactly the same observational data set specified in the paper of Wilson & Van Hamme (2014) and the assumption that the period P rises uniformly with time: P(t) = P0 + (tM0), where = const.

Model I uses the phenomenological model for the BV LCs (Eq. (18)), whilst the hybrid Model II uses synthetic LCs computed by Wilson & Van Hamme (2014). This technique is the one used, for example, by Zasche et al. (2014) and Zasche (2015). Both models assume strictly circular orbits of components and synchronous rotation. They also suppose that the shape of LCs given, namely, by the variable geometry of the orbiting system, are more or less constant.

In contrast, detailed inspection of the LCs proves that there are apparent changes on various time scales (see Fig. 14 and BV LCs in Fig. 3 in Wilson & Van Hamme 2014). Some of them can be attributed to instrumental effects such as red noise (Pont et al. 2006) or incomplete detrending of LCs. Others are due to intrinsic changes in the object itself, such as chromospheric and spot activity or shifts due to unsteady mass transfer between components (see the right part of Fig. 12).

thumbnail Fig. 14

Residua of BV LCs of Binnendijk (1971), Winkler (1977), and Awadalla & Yamasaki (1984) clearly show variations in LC shapes. B and V residuals are denoted as circles and open squares, respectively.

6.4.2. Unified ephemerides. Virtual eclipse times

All quantitative results achieved by the three discussed models that can be directly compared are listed in Table 4. It is apparent that especially the increase in the orbital period , the period P0 = P(t = M0), and the time of the basic primary minimum M0 derived from all data, from only LCs, and from eclipse times, are identical within their uncertainties. The results derived from the analysis of RV measurements are in very good concordance as well. We are convinced that this conclusion is valid not only for XY Boo  but also for other ESs.

There is an apparent identical discrepancy found by this paper and Wilson & Van Hamme (2014) in the value of the mean deceleration derived only from LC shifts and/or eclipse times (Tm) (see Table 4, rows LC and Tm): LC = 1.506(24) × 10-9, Tm = 1.745(26) × 10-9. The first possible explanation, which it is the result of different distribution of LCs and time data along the O–C diagram (see Figs. 13, 12) and the inconstancy of where ¨P>0\hbox{$\ddot{P}>0$}), has proven to be false. We tested the period models with ¨P\hbox{$\ddot{P}$} term and found that ¨P=0(4)×10-15\hbox{$\ddot{P}=0(4) \times 10^{-15}$} d-1. The effect is then very likely caused by the presence of the pair of eclipse minima of Hinderer (1960) at the very beginning of the O–C diagram (see Fig. 12) that significantly deviates from its parabolic prediction. However, our thorough inspection of the original Hinderer’s measurements showed that both eclipse times are correct.

Models I and II enabled calculation of “virtual eclipse times” for selected subsets of photometric and RV data using the relation Eq. (9). These times, including their uncertainties, are listed in Table 6, together with eclipse times8 published by authors of BV photometries – Binnendijk (1971), Winkler (1977), Awadalla & Yamasaki (1984) The results, as well as their uncertainties obtained by Models I and II, are nearly identical. The last ones are very important since they enable the evaluation of real inner accuracy of individual times, the construction of precise O–C diagrams (Mikulášek et al. 2011b), and the discussion of subtle problems of stability of the period, variability of LCs, etc. There is no systematic difference between the published eclipse timings (see the first column of Table 6) and those derived by our models.

The excellent agreement of comparable results obtained by all three models proves that they can be used as good alternatives in solving common tasks of the eclipse system period analysis.

6.4.3. BV and RV phase curves

The differences between the phenomenological and the synthetic BV LCs are insignificant, because they represent about 0.5% of the amplitude of the light changes (see Fig. 11). The parameters of the phenomenological model of both mentioned LCs are in Table 3. It seems that phenomenological modelling is a good method for expressing common types of LCs using only very few parameters.

The fit of the observed BV LCs by Model I is fairly good (see Fig. 10), since the scatter of normal points in B and V are better than 0.004 and 0.003 mag. The larger part of this scatter is caused by the erratic variability of LCs due to stellar activity and the inconstant O’Connell effect. The changes in LC shapes are clearly visible in the phase diagram of BV residuals (Fig. 14) for all three photometric data sets.

The detailed inspection of measurements of Binnendijk (1971) shows weak and variable O’Connell effect (data from 1968 and 1969 displayed O’Connell effect of the opposite sign), whilst residuals of Winkler (1977) and Awadalla & Yamasaki (1984) are dominated by a double wave going in antiphase. It seems that such seasonal variations in LC shapes are quite common. Sometimes they can be much more dramatic (see e.g. changes in LCs of notorious XY UMa, Lister et al. 2001). These LC changes lowered the accuracy of the fit and eclipse timing determination.

Still unpublished observations of XY Boo obtained during one month Mikulášek et al. (in prep.) prove that the time scale of LC changes is usually shorter than several days.

7. Conclusions

Phenomenological modelling presents an admissible alternative for the solution of selected ES research tasks. We showed that:

  • 1.

    The estimation of parameters and their uncertainties obtainedby our phenomenological (and also hybrid) modelling and otherwell-proven methods, including solutions by sophisticatedphysical models of ESs, are almost the same (see e.g.Sects. 4.3 and 6.4).

  • 2.

    Phenomenological modelling is based on the minimization of the χ2 sum, which enables simultaneous processing of different sources of phase information (complete LCs and their segments, RV curves, individual mid-eclipse times etc.).

  • 3.

    Simple model function of eclipses (Eq. (13)) can be also used for good determination of mid-eclipse times and their uncertainties of individual observations of stellar eclipses (our model has been standardly used for the determination of time minima of original observations of EBs by Variable star and Exoplanet Section of Czech Astronomical Society since 2011, Brát et al. 2012).

  • 4.

    As a by-product of the phenomenological modelling, we obtain the list of virtual minima times derived from LC data, which help us, among others, to quickly check the selected phenomenological model of LCs or period variability (see Sect. 6.4.2).

  • 5.

    Light curve model parameters can be used for the apt description of both one-colour and multicolour LCs of ESs. For example, shapes of ten curves of HST spectrophotometry (320–980 nm) of the exoplanet transit of HD 209 458b, taken from Knutson et al. (2007), are determined by only nine parameters (see Mikulášek et al. 2015). We can use this application, for instance, for describing and classifying the observed LCs of ESs.

  • 6.

    Knowing the template of LCs (from physical models or observed LC of superior quality), we can quickly modify the presented phenomenological model and establish a hybrid model (Model II in Sect. 6) with the diminished number of free parameters. As we showed in Sect. 6, the results of this approach are the same results as physical or pure phenomenological modelling. However, applying the latter method is much faster.

  • 7.

    Phenomenological (and also hybrid) modelling of an ES could solve standard tasks of ES research (based on every sources of phase information), especially the improvement of ES ephemeris for standard period models of

    • systems with constant periods,

    • systems in a steady regime of mass and angular momentum transfer between components (e.g. XY Boo, Sect. 6),

    • systems with erratic changes of period – the phase function is here approximated using so-called O–C time shifts of the observed phase curves versus the predicted LC derived by the period model with fixed parameters (typically P0,M0 or P0,M0,) – see Eq. (9), and Sect. 6.4.2,

    • eclipsing systems influenced by the gravitational attraction of a third body (see the preliminary paper about AR Aur in Mikulášek et al. 2011c, and the discussion in Wilson & VanHamme 2014)9.

Phenomenological modelling of eclipsing systems presented above has, among others, the following disadvantages and limitations:

  • The method does not provide direct information about thegeometry and physics of ESs, especially the inclination of thesystem, relative radii, temperatures, and the form of components.

  • The phenomenological models of EBs LCs in their proposed form are not able to achieve an accuracy better than 0.5% of their total amplitude of them (the accuracy of Kepler or CoRoT LCs is better by one or more orders). Nevertheless, such uncertainty is fully acceptable for many applications.

  • It may be a bit troublesome for beginners to select the optimum model of LCs. Using an inappropriate model may lead to nonconvergent or bad solutions.

  • Phenomenological modelling of ESs has not been made available in as user-friendly format as other modern physical models (e.g. Hadrava 2004; Prša et al. 2008, 2011; Wilson & Devinney 1971) until now.

Influenced by modern physical models, phenomenological modelling has been developed gradually since 2008 (Mikulášek et al. 2008), its elements and principles were used in the majority of papers of this paper’s author. The detailed applications of the presented bases of the method will be published elsewhere.


1

The presented models may also be applied for expressing light intensity variations, but the treatment of data in the magnitude domain is more straightforward, and the accuracies of results are almost the same.

2

Observed orbital period may change due to possible transfer of matter between components or light time effect.

3

Applying so-called phase rectification (Mikulášek et al., in prep.) we are also able to solve eccentric systems. The technique symmetrizes LCs of EBs with components moving unevenly on their eccentric orbits by the rectification of their phase functions. The method of phase rectification enables also the effective analysis of apsidal EB motion.

4

This assumption is fulfilled only partially. We can mention cyclical variations of instantaneous LCs with other than orbital period as the gradual change of the ES geometry due to asymptotic motion of double systems orbiting on eccentric orbits, possible asynchronous rotation of spotted components in wide systems, and possible pulsations of the components. Eclipsing binary LCs may also vary erratically because of chromospheric activity (see e.g. Sect. 6.4.3), time-dependent spottiness of the components, or changes in streams or disks around the stars. Neglecting of above mentioned effects introduce as a rule some extra noise in the period analyses and deteriorate the accuracy of the determination model parameters.

5

If p>q, then k=pqhk=0.\hbox{$\sum_{k\,=\,p}^q\,h_k=0.$}

6

Wilson & Van Hamme (2014) used for modelling of phase function ϑ and its inversion Θ an exact solution of the basic equation, Eq. (1): ϑ(t) = -1ln[1 + /P0 (tM0)], Θ(E) = P0/[exp(E−1)−1] + M0. Nevertheless, our models for ϑ and Θ (see Eqs. (17)) do not differ by more than 1.9 × 10-5 in the phase function and 0.3 s in time prediction from the exact ones, so they can be considered to be identical.

7

We did not distinguish between primary and secondary light minima because their depths are nearly the same.

8

All the eclipse times were given without quoting their uncertainties.

9

The detailed study about the triple star AR Aur is in preparation. For formulae see Liška et al. (2015).

Acknowledgments

The author is very indebted to many of his collaborators, namely to M. Zejda, S. de Villiers, M. Chrastina, M. Jagelka, J. Krtička, J. Liška, E. Paunzen, T. Pribulla, S.-B. Qian, M. Vaňko, R. E. Wilson, L.-Y. Zhu, and others, for their help with this arduous subject. The work on the paper was partly sponsored by the grants Ministry of education of Czech Republic II LH14300 and GA ČR 13-10589S, T.P. and M.V. would like to thank the project APVV-0158-11.

References

  1. Andronov, I. L. 2010, in Int. Conf. KOLOS-2010 Abstr. Booklet, Snina, Slovakia, 1 [Google Scholar]
  2. Andronov, I. L. 2012, Astrophysics, 55, 536 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  3. Awadalla, N. S., & Yamasaki, A. 1984, Ap&SS, 107, 347 [NASA ADS] [CrossRef] [Google Scholar]
  4. Beaky, M. M., & Koju, V. 2012, AAS Meeting, 220, 333.03 [Google Scholar]
  5. Bevington, P. R., & Robinson, D. K. 2003, in Data Reduction and Error Analysis for the Physical Sciences, 3rd edn. (McGraw-Hill) [Google Scholar]
  6. Binnendijk, L. 1971, AJ, 76, 923 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bradstreet, D. H. 2005, The Society for Astronomical Sciences 24th Annual Symposium on Telescope Science (SAS), 23 [Google Scholar]
  8. Brát, L., Mikulášek, Z., & Pejcha, O. 2012, http://var2.astro.cz/library/ 1350745528_ebfit.pdf [Google Scholar]
  9. Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M. 2000, ApJ, 529, L45 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  10. Chrastina, M., Mikulášek, Z., & Zejda, M. 2013, in Setting a New Standard in the Analysis of Binary Stars, eds. K. Pavlovski, A. Tkachenko, & G. Torres, EAS Pub. Ser., 64, 383 [Google Scholar]
  11. Claret, A., & Gimenénez, V. 1993, A&A, 277, 487 [NASA ADS] [Google Scholar]
  12. Conroy, K. E., Prša, A., Stassun, K. G., et al. 2014, AJ, 147, 45 [NASA ADS] [CrossRef] [Google Scholar]
  13. Davidge, T. J., & Milone, E. F. 1984, ApJS, 55, 571 [NASA ADS] [CrossRef] [Google Scholar]
  14. Drechsel, H., Rahe, J., Wargau, W., & Wolf, B. 1982, A&A, 110, 246 [NASA ADS] [Google Scholar]
  15. Goldman, I., & Mazeh, T. 1991, ApJ, 376, 260 [NASA ADS] [CrossRef] [Google Scholar]
  16. Hadrava, P. 2004, Publ. Astron. Inst. Czechosl. Acad. Sci., 92, 1 [Google Scholar]
  17. Hartkopf, W. I., McAlister, H. A., & Franz, O. G. 1989, AJ, 98, 1014 [NASA ADS] [CrossRef] [Google Scholar]
  18. Hayashi, F. 2000, Econometrics (Princeton Univ. Press) [Google Scholar]
  19. Henry, G. W., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2000, ApJ, 529, L41 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  20. Hilditch, R. W. 2001, An Introduction to Close Binary Stars (Cambridge University Press) [Google Scholar]
  21. Hinderer, F. 1960, Obs., 43, 161 [Google Scholar]
  22. Hoffmeister, C. 1935, Astron. Nachr., 255, 401 [Google Scholar]
  23. Kallrath, J., & Milone, E. F. 1998, Eclipsing Binary Stars, Modelling and Analysis (New York, Berlin, Heidelberg: Springer-Verlag) [Google Scholar]
  24. Kalimeris, A., Rovithis-Livaniou, H., & Rovithis, P. 1997, A&A, 35, 69 [Google Scholar]
  25. Kaluzny, J., & Rucinski, S. M. 1986, AJ, 92, 666 [NASA ADS] [CrossRef] [Google Scholar]
  26. Kholopov, P. N. 1981, Premennyje zvezdy, 21, 465 [NASA ADS] [Google Scholar]
  27. Kopal, Z. 1978, Dynamics of Close Binary Systems (Dordrecht, Holland: Reidel) [Google Scholar]
  28. Knutson, H. A., Charbonneau, D., Noyes, R. W., et al. 2007, ApJ, 655, 564 [NASA ADS] [CrossRef] [Google Scholar]
  29. Kwee, K. K., & van Woerden, H. 1956, Bull. Astron. Inst. Neth., 9, 252 [Google Scholar]
  30. Lacy, C. H. S., Torres, G., Claret, A., & Sabby, J. A. 2003, AJ, 126, 1905 [NASA ADS] [CrossRef] [Google Scholar]
  31. Lister, T. A., Collier, C. A., & Hilditch, R. W. 2001, MNRAS, 326, 1489 [NASA ADS] [CrossRef] [Google Scholar]
  32. Liška, J., Skarka, M., Mikulášek, Z., Zejda, M., & Chrastina, M. 2015, A&A, submitted [arXiv:1502.03331] [Google Scholar]
  33. Mader, J. A., Torres, G., Marschall, L. A., & Rizvi, A. 2005, AJ, 130, 234 [NASA ADS] [CrossRef] [Google Scholar]
  34. Mandel, K., & Agol, E. 2002, ApJ, 580, 171 [Google Scholar]
  35. McLaughlin, D. B. 1924, ApJ, 60, 22 [NASA ADS] [CrossRef] [Google Scholar]
  36. McLean, B. J., & Hilditch, R. W. 1983, MNRAS, 203, 1 [NASA ADS] [CrossRef] [Google Scholar]
  37. Mikulášek, Z. 2007, Odessa Astron. Publications, 20, 138 [NASA ADS] [Google Scholar]
  38. Mikulášek, Z., & Gráf, T. 2005, Contr. Astron. Obs. Skalnaté Pleso, 35, 83 [NASA ADS] [Google Scholar]
  39. Mikulášek, Z., & Zejda, M. 2013, in Úvod do studia proměnných hvězd (Brno: Masaryk University) [Google Scholar]
  40. Mikulášek, Z., Wolf, M., Zejda, M., & Pecharová, P. 2007a, Ap&SS, 304, 363 [Google Scholar]
  41. Mikulášek, Z., Zverko, J., Krtička, J., et al. 2007b, in Physics of Magnetic Stars, Proc. Conf. Nizhnyi Arkhyz 28–31 Aug. 2006, eds. I. Romanyuk, & D. Kudryavtsev, 300 [Google Scholar]
  42. Mikulášek, Z., Krtička, J., Henry, G. W., et al. 2008, A&A, 485, 585 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Mikulášek, Z., Krtička, J., Henry, G. W., et al. 2011a, A&A, 534, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Mikulášek, Z., Zejda, M., Qian, S.-B., & Zhu, L.-Y. 2011b, in 9th Pacific Rim Conf. on Stellar Astrophysics, eds. S. Qian, K. Leung, L. Zhu, & S. Kwok (San Francisco: ASPS), ASP Conf. Ser., 451, 111 [Google Scholar]
  45. Mikulášek, Z., Žižňovský, J., Zejda M., et al. 2011c, in Magnetic Stars. Proc. Int. Conf., eds. I. I. Romanyuk, & D. O. Kudryavtsev (Nizhny Arkhyz: SAO), 431 [Google Scholar]
  46. Mikulášek, Z., Zejda, M., & Janík, J. 2012a, in From Interacting Binaries to Exoplanets: Essential Modeling Tools, eds. M. T. Richards, & I. Hubený (Cambridge: Cambridge Univ. Press), IAU Symp., 282, 391 [Google Scholar]
  47. Mikulášek, Z., Gráf, T., Zejda, M., et al. 2012b, ArXiv e-prints [arXiv:1212.5527] [Google Scholar]
  48. Mikulášek, Z., Zejda, M., Chrastina, M., Qian, S.-B., & Zhu, L.-Y. 2013a, in Setting a New Standard in the Analysis of Binary Stars, eds. K. Pavlovski, A. Tkachenko, & G. Torres, EAS Pub. Ser., 64, 299 [Google Scholar]
  49. Mikulášek, Z., Zejda, M., Zhu, L. Y., et al. 2013b, Central European Astrophysical Bulletin, 37, 145 [NASA ADS] [Google Scholar]
  50. Mikulášek, Z., Chrastina, M., Liška, J., et al. 2014, Contr. Astron. Obs. Skalnaté Pleso, 43, 382 [NASA ADS] [Google Scholar]
  51. Mikulášek, Z., Zejda, M., Pribulla, T., et al. 2015, in ASP Conf. Ser., 496, 176 [Google Scholar]
  52. Molík, P., & Wolf, M. 1998, IBVS, 4640, 1 [NASA ADS] [Google Scholar]
  53. Nedoroščík, J., Vaňko, M., & Parimucha, Š. 2012, in From Interacting Binaries to Exoplanets: Essential Modeling Tools, eds. M. T. Richards, & I. Hubený (Cambridge: Cambridge Univ. Press), IAU Symp., 282, 73 [Google Scholar]
  54. Nordstrom, B., & Johansen, K. T. 1994, A&A, 282, 787 [NASA ADS] [Google Scholar]
  55. North, P. 1984, A&AS, 55, 259 [NASA ADS] [Google Scholar]
  56. O’Connell, D. J. K. 1979, Ric. Astron. Specola Vaticana, 8, 563 [Google Scholar]
  57. Pollaco, D. L., & Bell, S. A. 1994, MNRAS, 267, 452 [NASA ADS] [CrossRef] [Google Scholar]
  58. Pont, F., Yucker, S., & Queloz, D. 2006, MNRAS, 373, 23 [Google Scholar]
  59. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical recipes in C++: the art of scientific computing [Google Scholar]
  60. Pribulla, T. 2012, in From Interacting Binaries to Exoplanets: Essetial Modelling Tools, eds. M. T. Richards, & I. Hubený (Cambridge: Cambridge Univ. Press), IAU Symp., 282, 279 [Google Scholar]
  61. Pribulla, T., Vaňko, M., Ammler-von Eiff, M., et al. 2012, Astron. Nachr., 333, 754 [NASA ADS] [CrossRef] [Google Scholar]
  62. Prša, A., & Zwitter, T. 2005, ApJ, 628, 426 [NASA ADS] [CrossRef] [Google Scholar]
  63. Prša, A., Guinan, E. F., Devinney, E. J., et al. 2008, ApJ, 687, 542 [NASA ADS] [CrossRef] [Google Scholar]
  64. Prša, A., Matijević, G., Latković, O., Vilardell, F., & Wils, P. 2011, Astrophysics Source Code Library [record ascl:1106.002] [Google Scholar]
  65. Qian, S.-B., He, J., Xiang, F., Ding, X., & Boonrucksar, S. 2005, AJ, 129, 1686 [NASA ADS] [CrossRef] [Google Scholar]
  66. Rucinski, S. M. 1973, Acta Astron., 23, 79 [NASA ADS] [Google Scholar]
  67. Russell, H. N., & Merill, J. E. 1952, Princeton. Obs. Contr., 26, 1 [Google Scholar]
  68. Samec, R. G., Gray, J. D., & Carrigan, B. J. 1996, Observatory, 116, 75 [Google Scholar]
  69. Selam, S. O. 2004, A&A, 416, 1097 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Tsesevich, V. P. 1950, Astron. Circ. Kazan, 100, 18 [Google Scholar]
  71. Tsesevich, V. P. 1954, Odessa Izv., 4, 137 [Google Scholar]
  72. Tsesevich, V. P. 1971, Eclipsing Binary Stars (Moscow: Nauka) [in Russian] [Google Scholar]
  73. Van Hamme, W. V., & Wilson, R. E. 2007, ApJ, 661, 1129 [Google Scholar]
  74. Vaňko, M., Maciejewski, G., Jakubík, M., et al. 2013, MNRAS, 432, 944 [NASA ADS] [CrossRef] [Google Scholar]
  75. Wilsey, N. J., & Beaky, M. M. 2009,Soc. Astron. Sci. Ann. Symp., 28, 107 [Google Scholar]
  76. Wilson, R. E., & Devinney, E. J. 1971, ApJ, 166, 605 [NASA ADS] [CrossRef] [Google Scholar]
  77. Wilson, R. E., & Van Hamme, W. 2014, ApJ, 780, 151 [NASA ADS] [CrossRef] [Google Scholar]
  78. Winkler, L. 1977, AJ, 82, 648 [NASA ADS] [CrossRef] [Google Scholar]
  79. Wood, R. 1965, Observatory, 85, 258 [NASA ADS] [Google Scholar]
  80. Yang, Y. G., Qian, S. B., & Zhu, L. Y. 2005, AJ, 130, 2252 [NASA ADS] [CrossRef] [Google Scholar]
  81. Zahn, J.-P. 1992, in Binaries as Tracers of Stellar Formation, eds. A. Duquennoy, & M. Mayor (Cambridge University Press), 253 [Google Scholar]
  82. Zasche, P. 2015, New Astron., 34, 253 [NASA ADS] [CrossRef] [Google Scholar]
  83. Zasche, P., & Wolf, M. 2013, A&A, 558, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Zasche, P., Wolf, M., Vraštil, J., et al. 2014, A&A, 572, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Zhu, L., Qian, S.-B., Mikulášek, Z., et al. 2010, AJ, 140, 215 [NASA ADS] [CrossRef] [Google Scholar]
  86. Zhu, L., Zejda, M., Mikulášek, Z., et al. 2012, AJ, 144, 37 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Parameters describing LCs of several EBs.

Table 2

List of alternatives for several phenomenological model functions of eclipses.

Table 3

Nine parameters of phenomenological (this paper) and synthetic LCs (Wilson & Van Hamme 2014) models.

Table 4

Comparison of common parameters derived by Wilson & Van Hamme (2014) and us.

Table 5

List of used observational data including their source, number, specification, and scatter with respect to phenomenological model phase curves (Model I) and synthetic phase curve (Model II).

Table 6

Comparison of the eclipse timings of XY Boo.

All Figures

thumbnail Fig. 1

Examples of elementary model LC Ψ = Fe/Ak for C = 0; see Eq. (13) allowed by our phenomenological EB models. Parameter d describes the duration of the eclipse, and Γ is a parameter expressing the kurtosis of the LC.

In the text
thumbnail Fig. 2

Fit of the V curve of the overcontact spotted EB EK Com (the secondary minimum is a transit). The LC of the close binary is affected by O’Connell and proximity effects. The synthetic model curve is depicted by a solid line, and the phenomenological LC fit is depicted by the dashed line. For more information see Sect. 4.3 and Table 1.

In the text
thumbnail Fig. 3

Simultaneous fit to 15 transits of extrasolar planet TrES-3b corrected for trends (see Eq. (10)). The fit parameters according to Eq. (14): A = 0.02725(19) mag, D = 0.0117(7), C = −0.17(7), and Γ = 1.58(12), the parameter K was not introduced. The black line is the fit, circles are the normal points, grey circles indicate individual measurements with the area proportional to their weights.

In the text
thumbnail Fig. 4

V LC of a famous EB consisting of a very hot nucleus of a planetary nebula and solar type star. For the fit of the LC strongly affected by proximity effects, we need only 8 parameters phenomenological. For details see Sect. 4.3 and Table 1.

In the text
thumbnail Fig. 5

Fit of the primary minimum of the detached EB AR Aur in V and B, represented by normal points (circles), by physical (grey solid lines), and phenomenological (dashed lines) LCs (see Table 1).

In the text
thumbnail Fig. 6

Ten-colour transit LCs for HD 209458b taken from Knutson et al. (2007). Effective wavelengths of individual colours are in nm. For the description of all 10 LCs we only need 9 parameters (courtesy of Mikulášek et al. 2015).

In the text
thumbnail Fig. 7

Dependence of the amplitude A of transit of HD 209458b (see Fig. 6) on Λ = 550 /λeff−1 proves the validity of the approximation Eq. (16) for that case.

In the text
thumbnail Fig. 8

Fit of the BVRcIc LCs of a semidetached short periodic “cool Algol” AV Del. Dark dots are observations, and grey lines indicate the 13-parameter phenomenological fit. For details see Table 1 and Sect. 4.3.

In the text
thumbnail Fig. 9

Difference in mag between the physical model primary minimum of AR Aur and various alternatives of phenomenological basic function listed in Table 2.

In the text
thumbnail Fig. 10

BV curves of XY Boo (see the list in Table 5). The phase is plotted according to our quadratic ephemeris (see Table 4). The areas of points are proportional to their weights. Open squares are normal points, each of which represents the mean of about 50 measurements. Grey lines are fits by our phenomenological model (see Eq. (18) and Table 3). ΔB and ΔV display residuals of BV magnitudes and normal points from the LC phenomenological model. The scale of residuals is two times larger than the measure of BV LCs.

In the text
thumbnail Fig. 11

Models of B and V LCs of XY Boo. Grey lines are BV model LCs simulated by the Wilson & Van Hamme (2014) physical model; dotted lines correspond to our phenomenological model (see Eq. (18) and Table 4); dashed lines are hypothetical BV LCs of XY Boo without mutual eclipses.

In the text
thumbnail Fig. 12

XY Boo timing residuals from the quadratic ephemeris (with a term). Published minima timings are denoted as squares, timings derived from LCs are circles whose areas are proportional to their weight.

In the text
thumbnail Fig. 13

XY Boo timing residuals from the linear ephemeris (without a term). Published minima timings are denoted as squares, timings derived from LCs are circles whose areas are proportional to their weight, and a timing derived from RV curves is labelled as a diamond.

In the text
thumbnail Fig. 14

Residua of BV LCs of Binnendijk (1971), Winkler (1977), and Awadalla & Yamasaki (1984) clearly show variations in LC shapes. B and V residuals are denoted as circles and open squares, respectively.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.