Benchmarking the power of amateur observatories for TTV exoplanets detection

Abstract

We perform an analysis of ∼80 000 photometric measurements for the following 10 stars hosting transiting planets: WASP-2, -4, -5, -52, Kelt-1, CoRoT-2, XO-2, TrES-1, HD 189733, GJ 436. Our analysis includes mainly transit light curves from the Exoplanet Transit Database, public photometry from the literature, and some proprietary photometry privately supplied by other authors. Half of these light curves were obtained by amateurs. From this photometry we derive 306 transit timing measurements, as well as improved planetary transit parameters. Additionally, for 6 of these 10 stars we present a set of radial velocity measurements obtained from the spectra stored in the HARPS, HARPS-N and SOPHIE archives using the HARPS–TERRA pipeline. Our analysis of these transit timing and radial velocity data did not reveal significant hints of additional orbiting bodies in almost all of the cases. In the WASP-4 case, we found hints of marginally significant TTV signals having amplitude 10–20 s, although their parameters are model dependent and uncertain, while radial velocities did not reveal statistically significant Doppler signals.

methods: data analysis, methods: statistical, techniques: photometric, techniques: radial velocities, surveys, planetary systems

INTRODUCTION

The first extrasolar planet, orbiting a solar-type star 51 Pegasi, was discovered by Mayor & Queloz (1995), based on the precision Doppler measurements of the ELODIE spectrograph. After that, the number of the detected exoplanets grew continuously, exceeding 1000 so far. In fact, 20 years ago a new rapidly growing domain of fundamental science was created, devoted to the exoplanet research.

Currently, most of the known exoplanetary candidates were detected by one of the two major techniques: the radial velocity (RV) method or transit method that increased its output in recent time thanks to the launch of specialized spacecraft CoRoT (ESA) and Kepler (NASA).

Unfortunately, both the Doppler and transit exoplanet detection methods require sophisticated and expensive instrumentation and remain practically inaccessible to the majority of the community. The RV method requires the use of extremely stable spectrographs. There are only about a dozen of such instruments in the world that are capable of detecting an exoplanet. The transit detection is also difficult. It requires precise alignment of the planetary orbit with the observer, and the transit event is rather short in time, although periodic. To detect a transiting planet, we have to observe lots of stars, and we necessarily deal with large numbers of null detections. This requires the use of specialized ground-based robotic telescopes or telescope networks, or space observatories such as the above-mentioned CoRoT and Kepler. Organizing such a campaign is a difficult task even for professional scientific teams.

In terms of the photometric accuracy and quality, ground-based telescopes are no match to spacecrafts such as Kepler. But all space projects have a common disadvantage: they are severely limited in time, while the sensitivity to weak planetary signals in the data degrades quickly when the planetary orbital period exceeds the observational time base. This condition means that Kepler data cannot reliably detect long-period planets, analogous to the giant planets of Solar system. In general, the responsibility for follow-up observations always returns to ground-based observatories. Ground-based observations of planetary transits are much less demanding than precision RV measurements, and of course much cheaper than projects such as Kepler. In fact, such observations are possible with commercially available equipment typically used by a significant community of ‘amateurs’. Thus, the exoplanetary hunt can potentially become ‘citizen science’ that can trigger a qualitative leap in the field.

However, amateurs are definitely not equipped to undertake classic transit surveys such as SuperWASP or others. But none the less they can provide a useful scientific contribution by means of the transit timing variation (TTV) method (Agol et al. 2005; Holman & Murray 2005). In this approach, there is a list of well-defined targets that are continuously monitored. Each target is a known host of a transiting planet, and its transits are regularly observed. If more planets orbit the host, they should induce perturbational effects on the motion of the transiter, causing observable delays to its transits, in comparison with a strictly periodic ephemeris. The TTV method allows the observation time to be used more efficiently, because we know what targets should be observed, and when.

This work represents an attempt to determine the practical efficiency of such an approach, based on the photometry data, taken mainly from the Exoplanet Transit Database (ETD) of the Czech Astronomical Union (Poddaný et al. 2010), http://var2.astro.cz/ETD/.

Currently, about 30 observatories are regularly contributing to the data base, including amateur as well as professional ones. This network offers telescopes of different apertures – from 20 cm to 2.6 m. Their locations are shown in Fig. 1. We must note that in terms of technical characteristics of the telescopes, there is no sharp boundary between the amateur and professional equipment, and the quality of the observations is often determined by the local astroclimate, which can have even more important effect than the telescope size. The observational programme involves currently ∼20 transiter hosts that are more or less regularly observed. For some of the stars, several years of observations are already available. However, the transit fitting algorithm used by ETD is criticized for its simplicity and imperfections (e.g. Petrucci et al. 2013). In this work, we present a new data reduction pipeline that handles subtle photometric effects (such as the red noise) more accurately. This pipeline was developed to deal with the photometry of a relatively poor or moderate quality, which is typical for the amateur data in ETD.

Figure 1.

World distribution of the observatories regularly contributing to the ETD data base.

We do not limit ourselves to only amateur observations or only ETD. We also use photometric data published in the literature. Moreover, for our targets we revealed a moderate amount of the spectra stored in the HARPS, HARPS-N and SOPHIE archives, and we derive RV data from them to provide and independent ‘calibration’ of our TTV results. However, in this work we have no goal to perform a fully self-consistent TTV+RV analysis. Instead, we aim to characterize the accuracy and reliability of the TTV data that we can derive from just the photometry.

The structure of the paper is as follows. In Section 2, we provide a detailed description of all the data that we include in our analysis. In Section 3, we introduce the algorithms used to process the photometric data. In Section 4, we present the TTV data derived from the photometry and describe the results of their periodogram analysis. In Section 5, we give the remaining fitted parameters of the planets considered in the work. In Section 6, we describe the RV data obtained for some of our targets on the base of the spectra found in the HARPS/HARPS-N/SOPHIE archives. In Section 7, we discuss in detail the case of WASP-4, for which we detected possible hints of a weak TTV signal.

THE SOURCE DATA

Our primary source of the photometric data is the ETD. We extracted a set of the best transit light curves from ETD for 10 selected stars. Also, we added in our analysis several public transit light curves that we found in the literature. Thus, our photometric data cover 10 targets with 306 transit light curves, and about 80 000 individual photometric measurements. The public data were presented in the works listed in Table 1, and most are stored in the Vizier data base. Among these 306 transit light curves, 161 (or roughly half) were obtained by amateurs, while 65 are from professional observatories contributed to ETD, and 80 were published in the listed literature.

Table 1.

Sources of the photometric data (except for ETD).

Target	References
WASP-2	Southworth et al. (2010)
WASP-4	Wilson et al. (2008); Gillon et al. (2009); Winn et al. (2009); Southworth et al. (2009b); Sanchis-Ojeda et al. (2011); Nikolov et al. (2012); Petrucci et al. (2013)
WASP-5	Southworth et al. (2009a)
WASP-52	No
Kelt-1	Siverd et al. (2012)
CoRoT-2	Gillon et al. (2010)
XO-2	Fernandez et al. (2009)
TrES-1	Winn, Holman & Roussanova (2007b)
HD 189733	Bakos et al. (2006), Winn et al. (2007a), Pont et al. (2007)
GJ 436	Bean et al. (2008), Shporer et al. (2009), Stevenson et al. (2012)

We need to note that the data from (Sanchis-Ojeda et al. 2011) also contain reprocessed photometry from Winn et al. (2009), so the original Winn et al. (2009) data were not actually included in our analysis. Also, we eventually decided to drop the Hubble Space Telescope data by Bean et al. (2008), because they all appeared to cover only small parts of a transit and could not produce good precision in the derived mid-times. The photometry presented by Petrucci et al. (2013) is not public, but it was kindly released to us by the authors. The WASP-4 data from Southworth et al. (2009b) appeared not very reliable due to the clock errors of the Danish telescope (Nikolov et al. 2012; Petrucci et al. 2013), and we believe that the data from Southworth et al. (2009a, 2010) should be treated with the same care. However, this photometry is rather accurate in itself, and thus it is still useful in constraining all transit parameters except for mid-times (e.g. by adding a fittable time offset to these light curves model relative to the other light curves). All timings in the photometric data were transformed to the BJD_TDB time stamps by means of the public idl software developed by Eastman, Siverd & Gaudi (2010).

Additionally, we used the precision RV obtained from the spectra of the HARPS archive, available for the following targets from our photometry sample: WASP-2, -4, -5, HD 189733, CoRoT-2. These spectra were processed with the advanced HARPS–TERRA pipeline that offers an improved RV accuracy (Anglada-Escudé & Butler 2012). The quality of these RV data, and their number per a target, appeared not very high. The meaning of the columns in the attached RV data file is given in Table 2.

Table 2.

Explanation of the RV data file (TERRA.rv).

column	explanation
1	observation time (BJD)
2, 3	derived RV and its uncertainty
4	a standardized name of the input RV file that includes the target name and the name of the spectrograph

For GJ 436, we found a large amount of high-quality HARPS RV data, similar to the one considered by Lanotte et al. (2014), as well as a large Keck RV time series recently published by Knutson et al. (2014). We acknowledge that cases such as this deserved to be investigated in detail in a separate work. In this paper, we only present our TTV and TERRA RV data for GJ 436, processed in the common simplistic way as for the other stars.

In this work, we do not perform the joint transits+RV fits, because it appeared that to perform such an analysis at a desirable level of quality, we must provide a solution to several non-trivial issues that fall outside the scope of this paper. These issues include e.g. the treatment of the correlational structure of the RV noise, and the treatment of the RV points affected by the Rossiter–McLaughlin effect. In this work, we used the RV data mainly as an independent source of information.

METHODS OF THE TRANSIT LIGHT CURVES ANALYSIS

Details of the transit model and its parametrization

Let us first adopt the approximation that the transiting planet moves along a straight line with constant velocity, thus neglecting the orbital curvature of its actual trajectory during the transit. The model of the linear motion easily predicts the separation between the centres of planet and star discs, δ, as a function of time. We have four kinematic characteristics of the transit that must be fitted: (i) the mid-time of the transit t_c, (ii) the duration of the transit t_d, defined as the time spent between the first and fourth contacts and (iii) the impact parameter b, measuring the smallest projected separation δ, and (iv) the projected planet radius r that simultaneously determines the transit depth and the geometry of the ingress/egress phases. Given these parameters and easy geometric constructions, the projected separation δ can be expressed as

\begin{eqnarray} \delta (t) = \sqrt{b^2 + \left[(1+r)^2-b^2\right] \tau ^2}, \quad \tau = 2 \frac{t-t_{\rm c}}{t_{\rm d}}, \quad |b| \le 1+r.\nonumber\\ \end{eqnarray}

(1)

This assumes that the star radius is unity. To simplify the formula for δ and simultaneously get rid of the non-trivial definition domain for b, our algorithm adopts internally a replacement of the impact parameter b:

\begin{equation} b = \frac{p}{\sqrt{1+p^2}} (1+r). \end{equation}

(2)

With the new parameter p we do not need to worry about its domain: any real value of p is physically meaningful. Thus we eliminate the danger that this parameter can walk to a forbidden domain during the fit. With p replacing b, we have

\begin{equation} \delta (t) = (1+r) \sqrt{\frac{\tau ^2+p^2}{1+p^2}}. \end{equation}

(3)

However, this approximation of δ(t) might be inaccurate due to the curvature of the planet trajectory. In this work, we include a correction to this formula, assuming that the planet moves along a circular orbit. In this approximation, the projected distance can be expressed by the same formulae (1) and (3), substituting the following quantity τ^′ instead of τ:

\begin{equation} \tau ^{\prime } = \frac{\sin \tau \alpha }{\sin \alpha }, \qquad \alpha =\pi \frac{t_{\rm d}}{P} \ll 1, \end{equation}

(4)

where P is the transiter's orbital period. The auxiliary angle α reflects the curvature of the circular orbit. In fact, this correction induced only a well negligible effect on our TTV measurements, but none the less our results correspond to the model with this correction included.

The trajectory curvature also depends on the orbital eccentricity, but usually this eccentricity is difficult to determine from transit observations reliably, and we have no other option except to assume that it is zero. A zero eccentricity is a good prior assumption for most of these short-period planets. In the cases when RV data are also available, the accurate eccentricity information can be obtained from the joint transit+RV fits, but we do not perform fits of such type in this work.

After the function δ(t) is calculated, we approximate the relative flux reduction with the use of the stellar limb-darkening model by Abubekerov & Gostev (2013). These authors provided theoretical formulae for various types of the limb-darkening effect, as well as a software library written in c. The library provides subroutines to compute the observed light flux reduction ΔL as a function of the eclipsing planet radius r and of the projected separation δ (assuming that the star radius is unit), as well as of the limb-darkening coefficients that depend on the selected model. Besides, this library provides partial derivatives of ΔL with regard to its arguments. These derivatives are necessary to compute the gradient of the likelihood function, which is also used by the transit fitter.

Thus for a transit light curve we have four fittable parameters of the planet or planetary orbit: the planet/star radii ratio r, the transit mid-time t_c, the transit duration t_d and the impact parameter replacer p. As the orbital period P is known for all our planets with a very good accuracy, we treat it as a fixed parameter in the formula (4).

Additionally, there are two parameters determining the limb-darkening model. In this work, we use a quadratic two-term model of the stellar limb darkening with two coefficients to be determined, A and B. The brightness of a point on the visible stellar disc, observed at a given separation from its centre, ρ, is modelled as

\begin{eqnarray} I(\rho ) &=& 1 - A (1-\mu ) - B (1-\mu )^2 = \nonumber \\ &=& 1-A-2B + (A+2B)\mu + B\rho ^2, \nonumber \\ \mu &=& \sqrt{1-\rho ^2}. \end{eqnarray}

(5)

However, in practice the brightness model (5) can easily turn non-physical, if the coefficients A and B are allowed to attain arbitrary values. This becomes a significant problem if our data are polluted by some systematic errors that are always difficult to foresee in advance. There are a couple of natural basic constraints that we place on the coefficients A and B to keep the model (5) physically reasonable:

\begin{eqnarray} I(\rho ) \ge 0, \quad \frac{\mathrm{d}I}{\mathrm{d}\rho } = \left[ 2B\mu - (A+2B) \right] \frac{\rho }{\mu } \le 0, \qquad \forall \rho \in [0,1].\nonumber\\ \end{eqnarray}

(6)

If the second condition (the one on the derivative) is satisfied, I(ρ) is a monotonic non-increasing function, so the condition I(ρ) ≥ 0 for ρ ∈ [0, 1] is then equivalent to I(1) = 1 − A − B ≥ 0. In the second condition, the expression in parenthesis is a linear function of μ, so to have it always non-positive for μ ∈ [0, 1], it is necessary and sufficient to have it non-positive at the boundaries μ = 0 (disc limb) and μ = 1 (disc centre). Finally, we have total of three elementary inequality constraints on A and B to satisfy

\begin{equation} A+B\le 1, \quad A+2B\ge 0, \quad A \ge 0. \end{equation}

(7)

Note that the coefficient B is allowed to be negative here (but B ≥ −1). We do not put an extra condition B ≥ 0. Negative values of B mean that the limb-darkening gradient is diminishing (in absolute value) closer to the limb, although it never turns positive (any ‘limb brightening’ is disallowed).

We may satisfy restrictions (7) by means of making a smooth replacement of the parameters A and B, such that the formulae of the replacement would disallow the conditions (7) to be broken. This can be reached, for example, by the following trigonometric replacement:

\begin{equation} A = \sin ^2\theta \, (1-\cos \varphi ), \quad B = \sin ^2\theta \cos \varphi . \end{equation}

(8)

Naturally, whatever real values the new auxiliary parameters θ and φ might attain, the resulting values of A and B always satisfy restrictions (7). From the other side, each point (A, B) in the domain (7) maps to some pair (θ, φ) with real values of the parameters:

\begin{equation} \sin ^2\theta = A+B, \quad \cos \varphi = \frac{B}{A+B}. \end{equation}

(9)

Note that from (7) it follows that A + B ≥ |B|, meaning that A + B is never negative.

Therefore, treating the auxiliary angles θ and φ as primary fittable parameters, we can satisfy the conditions (7) automatically. The result of the fitting would correspond to either an internal point of the domain (7), or to some point on its boundary, if the actual data suggest to move the solution to a non-physical domain due to e.g. their poor statistical quality or systematic errors. The boundary points correspond to certain special values of θ and φ that can be easily identified: θ = ±π/2 + πk on the first boundary of (7), φ = π + 2πk on the second boundary, and φ = 2πk at the third boundary.

Details of the fitting procedure

In our analysis we are interested in obtaining the TTV deviations with maximum accuracy achievable with the available photometric data. However, the ETD data, even after the pre-selection, often have only a moderate quality, performing a complete and independent fit for each transit light curve is not a good option. The accuracy of thus-obtained transit parameters would often be poorly constrained, and this would impact the accuracy of the fitted TTV data as well. Moreover, the transit fits for individual light curves sometimes do not even converge to a reasonable result e.g. due to incomplete coverage of the transit or presence of significant curved trends.

To overcome these issues, we adopt in this work the following approach. We perform a joint fit of all transit light curves, available for a given star, assuming that most of the transit parameters, except for the mid-times, are equal for different light curves. Such a ‘shared’ transit parameter is still fittable, taking into account the constraint of its values being equal between different light curves. The mid-times are fitted individually for each light curve, i.e. they remain unconstrained.

Such an approach uses the full statistical power of all photometric data, available for a given star, to fit the shape of the transit curve, while still fitting the TTV offsets of individual transits separately from each other. We can note the following potential weaknesses of this approach.

• It is not taken into account that the limb-darkening coefficients depend on the photometry band, which are different for different light curves. Therefore, the results of such fitting method would refer to some averaged limb darkening, possibly introducing minor modelling errors in individual light curves. From the other side, the limb-darkening effect in the transit curve is always symmetric relative to the mid-time, implying that its impact on the derived mid-times should be small, if not negligible. In fact, we noticed that even fitting of a transit model without any limb darkening at all does not change the derived TTV data significantly (beyond the estimated parametric uncertainties). However, in the final version of our algorithm we decided to disentangle the limb-darkening coefficients for the best light curves (those that have rms smaller than 10 per cent of the transit depth) and fit these coefficients independently.

• It is not taken into account that transit duration may also be subject to variations, such as the transit mid-time. However, in practice, it appeared that the accuracy of the transit duration estimations was roughly an order of magnitude worse than those of the mid-times. Moreover, some light curves cover the transit event only partially, implying that its duration would remain almost unconstrained when fitted independently. We do not address transit duration variations in this work, focusing our attention on the TTV.

• The mid-time estimates, obtained in such a manner, are not necessarily uncorrelated. These mid-times are correlated with the remaining transit parameters, which are now shared between the light curves. Through this effect, some cross-correlation between the estimated mid-times may appear. This creates a risk of unexpected statistical effects in the derived TTV data, such as the non-white noise. However, the magnitude of these TTV correlations usually remains very small (maximum of a few per cent, and ∼0.1 per cent in average), and no deviations from the white noise are seen in the periodograms of the TTV data. More significant (>10 per cent) TTV correlations can appear for light curves that offer only a partial coverage of the transit, but such light curves represent only a minor fraction of our data.

In addition to the planetary transit itself, our light-curve models include a polynomial trend with fittable coefficients. Such a trend is necessary to take into account various drifting effects e.g. the effect of airmass or other types of systematic variations that appear frequently in our data. Each transit light curve has an individual fittable trend with a separate set of trend coefficients. After some experimenting with the data, we found that cubic trends represent a good compromise between the model adequacy and its parametric complexity. We did not try to reduce this systematic photometric variation based on its correlation with the airmass function: while for some light curves such a correlation looked clear, for others it was not obvious, indicating that other systematic effects were in the game.

The light-curve model is relatively complicated, while the photometry used in this work is not of a very good quality. In practice, we often faced difficulties causing the fitter to be trapped in the local maxima of the likelihood function, which was related to an unrealistic branch of the solutions. For example, without special care, see below, we frequently obtained non-physical solutions corresponding to a grazing transit with r ≫ 1. Also, some light curves do not cover the complete transit, and in such cases the fitting of the transit mid-time is complicated by its strong correlation with the coefficients of the polynomial trend. To avoid such traps, we worked out the following sequence of auxiliary preliminary fits.

• Perform a preliminary fit constraining the mid-times at a regular grid (with free scale and offset), implying all TTV residuals are zero by definition; fixing the impact parameter p at an intermediary value of 1 (⁠|$b\approx 1/\sqrt{2}$|⁠), and fixing the limb-darkening coefficients at A = B = 0.25 (or θ = π/4 and φ = π/3).

• Refit after releasing the mid-times and impact parameter, but still holding the limb-darkening coefficients fixed.

• Refit after releasing the limb-darkening coefficients (binding them across different light curves).

• Refit after full release of the limb-darkening coefficients for the best light curves (those with rms <0.1 of the transit depth).

After that, we also apply the red-noise detection and fitting procedure as described in the section below.

Reduction of the red noise

It is already known well that stellar photometric data usually include a correlated (‘red’) noise component that has an important effect on the fitted transit parameters (Pont, Zucker & Queloz 2006; Winn et al. 2008; Carter & Winn 2009). An easy technique to calculate the effect of the red noise on the fitting uncertainties was introduced in these works. Although these authors avoid making restrictive assumptions concerning the correlation structure of the red noise, their method still remains rather simplistic and it does not take into account important effects. They mainly focus on a more accurate determination of the uncertainties in the transit fit, which are typically underestimated without a proper treatment of the red noise. But the red noise may also induce biases in the fitted parameters themselves. This biasing effect was already noted when processing Doppler data affected by red noise in a generally similar way (e.g. Baluev 2011, 2013b). Also, there is not an obvious way to control the validity of the underlying assumptions in the traditional approach of the photometric red noise reduction. This approach does not offer a complete noise model that could be verified for accuracy. Additionally, the original method by Pont et al. (2006) relies on the out-of-transit photometric data, which are very limited in our case.

In this work, we apply the approach based on a parametric modelling of the noise covariance matrix. This is an adaptation of the red-noise fitting technique from Baluev (2011, 2013b, 2015) that was developed to handle exoplanetary RV data, in which the noise correlations may appear e.g. due to the stellar activity. In this approach, we approximate the red noise by a stationary Gaussian random process with a covariance function of a given functional form. The covariance coefficient between two arbitrary photometric observations x_i and x_j, acquired at the times t_i and t_j, is modelled as

\begin{eqnarray} V_{ij} = \mathop {\rm Cov}\nolimits (x_i,x_j) = \sigma _{i,\rm wht}^2(p_{\rm wht})\delta _{ij} + V_{ij,\rm red}(p_{\rm red},\tau ), \nonumber \\ V_{ij,\rm red} = p_{\rm red} R_{ij}(\tau ), \quad R_{ij} = \rho \left(\frac{t_i-t_j}{\tau }\right), \end{eqnarray}

(10)

with p_wht, p_red and τ being free fittable parameters. Here, the first noise term, |$\sigma _{i,\rm wht}^2$|⁠, represents the white fraction of the noise, detailed below. The red noise is given by the second term, V_ij, which also depends on the function ρ(t) that represents an adopted shape of the correlation function. We use mainly the exponential correlation function ρ(t) = exp ( − |t|). This not a unique choice, but in practice this model of the red noise usually appears adequate.

The white-noise term |$\sigma _{i,\rm wht}^2$| in equation (10) requires a separate discussion. We believe that physically the so-called additive model would be more suitable here. In this model, the value |$\sigma _{i,\rm wht}^2$| is determined as a sum of the stated instrumental variance and of an unknown fittable ‘jitter’ variance, generated by e.g. Earth atmospheric instability, or by unassessed instrumental instability, or by short-term (minutes to hours) intrinsic variations of the stellar flux, whenever these variations can be treated as a white noise. Although, we must note that in practice the instrumental uncertainties in our data do not look very trustable, and often they are just omitted, forcing us to assume that the measurements have just equal uncertainties. Therefore, in any case we should not expect that our white noise can be accurately modelled from the physical point of view. In these circumstances, the noise model should be chosen mainly on the basis of mathematical simplicity or usefulness, relying on only minimal or no knowledge of the underlying physics. We decided to use in our work a regularized noise model defined in Baluev (2015). In the majority of practical cases, this model should be equivalent to the additive noise model. This regularized model proved rather resistant with respect to various pitfalls appearing during the fitting procedure.

The fitting of the compound model, involving simultaneously the models of the transit curve and of the photometric noise, is done by means of the maximum-likelihood approach with details given in Baluev (2015). That work was devoted primarily to the analysis of the RV data, but mathematically the methods that we are using here are identical.

However, after an attempt to apply this technique in practice literally, we faced the problem that the red noise could not be fitted in many of our light curves. Only 1/3 to 1/2 of the light curves allowed for a reliable red noise fit, while in the other cases, the red noise was either not detectable, or its estimated magnitude became pretty small (in comparison with the uncertainty). Such cases lead to model degeneracies and decrease the reliability of the entire fit for a given star. Thus, we decided to add the red noise term to the photometric model only in the cases when it was justified.

We added to our analysis pipeline a set of auxiliary fits in order to identify the light curves, in which the red noise could be modelled more or less reliably. Namely, we tried to perform a series of test fits, latterly adding a red noise term to the model of each transit light curve. If the red noise could not be estimated reliably, more accurately if the relative uncertainty of its estimated magnitude exceeded 2/3, the relevant red noise term was removed from the model and this light curve was further treated as having purely white noise. Otherwise, we proceeded to the next test fit preserving this red noise term in the model. In the end, after all individual light curves were processed in such a way, we performed one more check of each red noise term (because some of them could turn insignificant after the other terms being added) and the final fit was made.

The quality of such an approach to the red noise reduction is examined in Fig. 2. The best way to observe a ‘non-white’ noise is to look at its power spectrum: the white noise should have constant power in average, while the red noise should demonstrate a systematic uprise to long periods. In our case we compute the so-called residual periodograms (see Baluev 2015) for each individual light curve, and then compute their average per each star involved in the analysis. These averaged periodograms are plotted in Fig. 2. For each star, we sequentially try three base photometric models: (i) just transits without any limb darkening plus the white noise; (ii) same as model (i) plus cubic trends and the effect of the quadratic limb darkening; (iii) same as model (ii) plus the red noise term.

Figure 2.

Examining the effect of the red noise in the light-curve photometry for a few of stars belonging to our sample. Please note that plots in the first column have different ordinate scale than those in the other two. See Section 3.3 for a detailed discussion.

We can clearly see from Fig. 2 that the long-term photometric variations are not reduced to just trends. Although the trends themselves are also important, after their removal the power spectra still contain an additional power in the period range of >10 min. The trends can only affect the period domain of >200–300 min (this corresponds to the typical time span of the light curves). In most cases, application of our red noise model remarkably suppresses the remaining excessive power, at least by a factor of 2. We must acknowledge that the reduction of this red noise looks rather difficult and still not entirely complete.

Some statistical information on the derived red noise characteristics now follows. About 60 per cent of light curves did not contain any detectable red noise at all. For the remaining 40 per cent of light curves, the red noise caused an increase of the TTV uncertainty by a factor of 1.3 on average (a median value). This is not very large, although a minor fraction of the mid-times demonstrated an increase in the uncertainty up to a factor of 2–3. Concerning the bias, appearing in the derived mid-times themselves (rather than uncertainties), its average value was only about 1/6 fraction of the uncertainty, and only a few points demonstrated a shift above 1σ.

Summarizing, in the majority of the cases, either there was not any detectable red noise, or the change in the TTV data was small, even if the red noise term was detectable. None the less, it is important that some individual TTV points are affected much more.

Testing the Gaussianity of the noise and clearing away the outliers

In the data analysis methods described above, we largely relied on the assumption that the photometric noise is Gaussian. This assumption needs some verification, which we performed in the following manner. After fitting all the data, we computed the residuals and normalized them by their relevant modelled uncertainty (square root of the sum of the estimated white and red noise variances). If the noise was Gaussian, the distribution of these normalized residuals should be close to the standard Gaussian, and vice versa. We consider the entire set of ∼80 000 normalized residuals for all stars involved in our analysis, and plot the relevant distribution in Fig. 3.

Figure 3.

Testing the normality of the entire ∼80 000 photometric data used in the work. Abscissa shows the value of the best-fitting residual, normalized by its total modelled variance (for white+red noise). The ordinate contains the normal quantiles of the empirical cumulative distribution of these normalized residuals. Perfectly Gaussian residuals should all lie close to the main diagonal that represents a standard normal distribution, while a deviation from the diagonal indicates a non-Gaussian noise.

As we can see, the empirical distribution is indeed almost Gaussian, except for the very tails (>4σ deviation). The empirical distribution has heavier tails. However, by removing data in the distribution tails we can reach good agreement with the Gaussian approximation. In fact, these extreme points represent outliers and should be removed in any case. The final results below correspond to the data with these 24 outliers removed. We also notice that the minor systematic deviations still remaining in the distribution tails are mainly due to Bakos et al. (2006) data for HD 189733. Visual investigation reveals that some of these light curves contain limited segments with unexpectedly increased photometric scatter. Removal of the Bakos et al. (2006) data makes the agreement of the residuals distribution with the Gaussian one very good.

DERIVED TTV DATA AND THEIR ANALYSIS

Using the fitting technique described above, we computed the mid-times for each observed transit, and placed these TTV measurements and other accompanying data in the online supplement to the paper. The supplement represents a file containing a text table. The meaning of the columns in this file is described in Table 3. This file only contains the parameters that are attached to individual light curves. The common fit parameters that are shared between different light curves are given separately in Section 5 below.

Table 3.

Explanation of the TTV data file (rn_t3l3.ttv).

Column	Explanation
	Primary data
1	Integer transit count (number of the transiter's revolutions, restarts from zero for each star)
2, 3	Fitted transit mid-time (BJDTDB) and its uncertainty
	Auxiliary data (other parameters of the fit)
4	A standardized name of the light-curve file that includes the date, the target name and the name of the observer or first author of a paper
5	The reference time T0 to which the following trend coefficients refer
6	Number of the following trend coefficients including the constant, i.e. the trend degree plus one (always four in this work)
7, 8	Fitted constant level of the magnitude and its uncertainty
9, 10	Fitted linear trend coefficient and its uncertainty (mag d−1)
11, 12	Fitted quadratic trend coefficient and its uncertainty (mag d−2)
13, 14	Fitted cubic trend coefficient and its uncertainty (mag d−3)
15, 16	Fitted limb-darkening coefficient A and its uncertainty
17, 18	Fitted limb-darkening coefficient B and its uncertainty
19	An adopted value of a scale parameter σscale needed to fully characterize the regularized model of the photometric noise, see Baluev (2015)
20, 21	Fitted photometric white jitter σ⋆, wht (def. in Baluev 2015) and its uncertainty
22, 23	Fitted photometric red jitter |$\sigma _{\star ,\rm red}=\sqrt{p_{\rm red}}$| and its uncertainty
24, 25	Fitted correlation time-scale τ of the red jitter (days)
26	rms of the best-fitting residuals for this light curve

The TTV residuals for these data are plotted in Fig. 4. In Fig. 5, we verify the normality of the distribution of the TTV noise present in these derived data in the way similar to Section 3.4. In the TTV residuals, we find no obvious outliers or deviations from the Gaussian distribution. The apparent excess in the distribution tails that can be seen in the left-hand panel of Fig. 5 is uncertain due to a relatively small number of TTV data points. None the less, we identified nine possible ‘candidate outliers’, and tried to remove them together with seven TTV points referring to the Southworth et al. (2009a,b, 2010) data. The TTV analysis discussed below was performed for the full TTV data set as well as for the reduced one.

Figure 4.

Graphs of the TTV residuals, based on the best-fitting transit ephemeris (values of P and |$t_{\rm c}^{\rm ref}$|⁠) from Table 4. The transits for which a significant red noise was detected in the photometry are labelled with a circle.

Figure 5.

Testing the normality of the TTV data derived in the work. The plot is similar to Fig. 3, but refers to the TTV residuals (relatively to a strictly periodic transit ephemeris) rather than to the photometric ones. In the right-hand panel, we removed the Southworth et al. (2009a,b, 2010) data (seven points), as well as a few potential outliers (nine points).

First, we made an effort to detect possible long-term curvature in our TTV. This was done by means of modelling the mid-times by a quadratic function (instead of only a linear function that would correspond to strictly periodic transits). None of the TTV time series demonstrated a long-term quadratic variation with at least 2σ significance (the most suspicious cases had 1.5σ to 1.7σ).

Secondly, we undertook a periodorgam analysis of the TTV residuals, in order to reveal possible periodic perturbations from any additional unknown planets orbiting these stars. We again applied ‘residual periodograms’ similar to the ones that were used to plot Fig. 2 above. Now our null model included a quadratic trend that is related to the period of the main (transiting) planet and white noise with a fittable additive ‘jitter’. The alternative model also included a sinusoidal signal at a probe frequency. None of these periodograms demonstrated any signs of a non-white noise, and in fact almost all appeared consistent with the pure noise model, although the cases revealing significant excessive TTV ‘jitter’ were rather frequent. No significant TTV periodicities were found in the data for any of the stars involved in the analysis, except for the WASP-4 case which is discussed below. The periodogram significance levels were estimated using the approach of Baluev (2008, 2009).

We believe that this result indicates a considerable improvement in the statistical quality of the TTV data and/or methods of the statistical analysis, because for the original TTV data from ETD we frequently detected spurious periodicities possibly appearing due to imperfections of the transit fitting algorithm.

IMPROVED TRANSIT CURVE PARAMETERS

In addition to the TTV data, we also provide the remaining best-fitting parameters of the exoplanetary transit light curves. These are the transit parameters that were shared between different light curves, and they are given in Table 4. In addition to the previously mentioned transit parameters r, t_d, b and the limb-darkening parameters A and B, this table also contains the following data: number of the transit light curves used in the analysis, number of the light curves in which a red noise term was robustly detected, the χ² (weighted rms) of the TTV residuals (relatively to the best linear fit of the derived mid-times) and the refined parameters of the transit emphemeris – transiter period P and a reference mid-time |$t_{\rm c}^{\rm ref}$|⁠. The last two quantities were obtained by requiring the mid-times to obey a strict linear relation with the transit count and fitting the coefficients of this relation. The value of |$t_{\rm c}^{\rm ref}$| is not unique due to the periodic nature of the transits. We chose such a reference transit, for which the uncertainty of t_c would be minimum. This condition simultaneously implies that the correlation between the estimated P and |$t_{\rm c}^{\rm ref}$| is zero. Note that it is not required to have actual transit observation at this |$t_{\rm c}^{\rm ref}$|⁠.

Table 4.

Fitted parameters of exoplanetary transit curves.

	Total	Assuming fittable TTV							Fixing TTV residuals at zeroa
	number of	Number of	Radii ratio	Half-duration	Impact par.		Mid-times correl.		Orbital period	Ref. mid-time
	transits	red noised	r = Rpl/R⋆	td/2 (d)	b	|$\sqrt{\chi ^2_{\rm TTV}}$|a	MAD/MAX		P (d)	|$t_{\rm c}^{\rm ref}$| [BJDTDB −
Transiter		transits								−2450000]	Comment
WASP-2 b	38	9	0.1355(31)	0.037 01(35)	0.7380(94)	1.38	0.000 56/0.0026		2.152 221 63(42)	5894.079 19(15)
WASP-4 b	43	20	0.154 95(32)	0.044 898(53)	0.136(28)	1.35	0.000 25/0.0039		1.338 231 624(68)	4966.782 814(21)
WASP-5 b	11	4	0.1136(13)	0.050 28(29)	0.446(40)	1.34	0.0070/0.10		1.628 429 53(52)	6446.988 68(17)
WASP-52 b	22	7	0.1629(44)	0.038 58(69)	0.598(32)	2.00	0.0020/0.016		1.749 7835(11)	6673.821 49(13)
Kelt-1 b	15	6	0.0783(14)	0.056 61(42)	0.998(44)b	2.56	0.000 92/0.017		1.217 494 48(80)	6093.134 64(19)	Showing |$\sqrt{1-b^2}$| instead of b, see note 2
CoRoT-2 b	20	7	0.1639(21)	0.047 36(36)	1.000(89)b	1.17	0.000 79/0.015		1.742 996 73(31)	5628.447 58(14)	Showing |$\sqrt{1-b^2}$| instead of b, See note 2
XO-2 b	25	10	0.1036(13)	0.055 63(32)	0.996(31)b	1.56	0.000 56/0.012		2.615 857 79(43)	5139.160 92(13)	Showing |$\sqrt{1-b^2}$| instead of b, See note 2
TrES-1 b	43	17	0.137 81(97)	0.052 25(18)	0.191(67)	1.54	0.000 38/0.0078		3.030 069 73(18)	5016.969 937(70)
HD189733 b	67	29	0.157 12(40)	0.037 576(48)	0.6636(19)	2.16	0.000 37/0.071		2.218 575 200(77)	3955.525 5511(88)	Some of Bakos et al. (2006) data contain many outliers and possibly non-Gaussian noise
GJ 436 b	22	8	0.088(10)	0.021 02(71)	0.806(31)	1.22	0.000 26/0.0019		2.643 898 46(44)	5280.175 68(17)	Most data are of poor quality (large trends), and a significant eccentricity of ∼0.15 is not taken into account

Notes.^aThe orbital periods, the reference mid-times, and the values of |$\chi ^2_{\rm TTV}$| were obtained after removal of the Southworth et al. (2009a,b, 2010) data, as these data are affected by clock errors.

^bImpact parameter b is close to zero and is thus a highly non-linear parameter here. Its estimation is severely non-Gaussian, and the formal uncertainty σ_b is much larger than b. To handle this peculiarity, the values in the column for b are replaced by |$a=\sqrt{1-b^2}$| and its uncertainty σ_a = (b/a)σ_b, which are more informative here.

The |$\chi ^2_{\rm TTV}$| values in Table 4 always exceed unity, so the fitted uncertainties might be moderately underestimated. In some cases this higher-than-expected scatter could be partly explained by small number of transits involved in the analysis, and in some part it might be due to additional unseen planets that may induce complicated TTV signals and are difficult to detect. Another explanation is incomplete reduction of the red noise and the effect of non-linearity of the transit model. The latter effect can introduce biases both in the estimated parameters as well as in their uncertainty estimations. So far we could not decide which explanation is more likely. Interestingly, this excessive TTV scatter persists and is important even in such a robustly fittable case such as HD 189733, in which the model non-linearity should be well suppressed thanks to a large number of available observations and light curves.

RV DATA AND THEIR ANALYSIS

In addition to the TTV time series from the photometry, we also derive RV data for some of our stars using spectra found in the HARPS, HARPS-N and SOPHIE archives. We process these spectra using the advanced HARPS–TERRA pipeline by Anglada-Escudé & Butler (2012). Public spectra were available for the following targets: WASP-2, 4, 5, HD 189733, GJ 436, CoRoT-2.

Most of these RV data appeared less accurate than the typical 1 m s⁻¹ precision demonstrated by HARPS. This is because many of these targets are rather faint. Moreover, it seems that adequate modelling of these data should involve non-trivial treatment of the RV noise. Frequently, these data are combined in short series acquired within a few hours or even shorter. The data within such single-night series should be significantly correlated (e.g. Nelson et al. 2014), and this effect can be easily detected in some of our data. Usually, these short-term runs were clearly intended to catch the Rossiter–McLaughlin effect during the planetary transit. On larger time-scales (days to weeks), the RV noise may still remain correlated as well, likely due to the stellar activity (see e.g. Baluev 2013b; Robertson et al. 2014), and this type of correlation is different from the one emerging at short time-scales.

The full modelling of all these effects is outside of the scope of this paper, which was intended to deal mainly with the photometry and TTV data. However, these RV data may carry useful information that can be helpful in verifying the results of our TTV analysis. In this work, we adopted a simplified ‘first look’ approach to the RV data analysis. First, we replaced the RV series acquired in a single night by their averages. This eliminated the need to model the Rossiter–McLaughlin effect, as well as the RV noise correlations appearing within a single night. Of course, such a procedure likely adds some systematic error to these ‘cumulative’ measurements, but in this work this is a satisfactory precision. Then we passed these data through a periodogram computing tool of the planetpack software (Baluev 2013a), taking into account the best-fitting contribution from transiting planet (i.e. including it in the null model of the periodogram).

In this analysis, none of the periodograms revealed hints of any additional variations in the data. One exception is the GJ 436 case. In this case, we have rather large amount of the HARPS data, as well as Keck data from (Knutson et al. 2014), and both these data sets demonstrated clear hints of a red noise, similar to the one investigated e.g. by Baluev (2013b). The same might be true for HD189733, for which we also have large amounts of the RV data. We plan to consider these cases in detail in a separate work.

THE WASP-4 CASE

The periodogram of the WASP-4 TTV data is shown in Fig. 6, and it reveals some marginally significant periodicities, corresponding to the TTV amplitudes ∼10–20 s. However, there is not any stable pattern of the peaks, as the periodograms are severely model dependent. Depending on which data subset we include in our analysis, we obtain different results. Note that the data by Southworth et al. (2009b) are not reliable for TTV studies due to the clock failures noticed for this telescope. This is why we also consider in Fig. 6 a reduced TTV time series obtained by removing these data. There are only four TTV points, but they have a good formal accuracy as derived from the photometric fit, and thus affect the TTV periodogram significantly. Also, results of the period analysis depend on various other subtleties e.g. on the degree of the polynomial trend used in the photometry model. On the other hand, although these variations are marginal and model dependent, we cannot just attribute them to the noise, as we did not observe anything similar in the other stars of this study.

Figure 6.

Periodograms of the WASP-4 TTV data, showing a few marginally significant peaks. The plot in the left-hand panel was based on all available transit light curves, the middle panel is for a reduced data set with the Southworth et al. (2009b) data removed and in the right-hand panel one TTV outlier detected in the (Petrucci et al. 2013) data was also removed. The thick horizontal lines label the significance levels of 1σ or FAP ≈ 31.7 per cent, 2σ or FAP ≈ 4.6 per cent, 3σ or FAP ≈ 0.27 per cent. The significance and the period for the tallest peak are printed in each panel. All periodograms involve a quadratic trend of the TTV delay with fittable trend coefficients.

Using the method of Section 4, we detected an outlier in the WASP-4 TTV data, owing to one of the Petrucci et al. (2013) light curves. This peculiar light curve was already noticed by Petrucci et al. (2013). Interestingly, the height of the periodogram peak at ∼5.14 d is increased after removal of this TTV outlier from the analysis.

The RV data from the HARPS archive did not reveal any significant signal in addition to the primary transiting planet. However, the star is rather faint for HARPS. Also, we should not expect a direct connection between the TTV and RV periods, because often TTV is an indirectly induced variation, owing to the planetary dynamical perturbations. Such a variation may appear due to e.g. a mean motion resonance, which is not always easily detected in RV (e.g. Anglada-Escudé, López-Morales & Chambers 2010).

Currently, we remain uncertain about the nature of these variations that are possibly present in the WASP-4 TTV data. Possibly, the TTV signal might be more complicated than just a sinusoid, and it should be modelled in the framework of the Newtonian N-body fits. Another possibility is star-spots inducing systematic perturbations in the timing measurements. In any case, it seems that we should keep tracking this target or maybe even focus increased attention on it by making further observations in future.

CONCLUSIONS AND DISCUSSION

As the main topic of the manuscript is the TTV exoplanet detection, it is interesting to investigate its efficiency relatively to the more classic methods such as the Doppler planet detection. An apparent TTV signal can be induced via two mechanisms: the dynamical perturbation on the transiter's motion and the light arrival time delay due to the finite light speed (the Roemer effect). The dynamical perturbations are difficult to predict in a general case, as they severely depend on many orbital parameters (e.g. they may drastically increase when planet move in a mean motion resonance). But the Roemer delay can be easily assessed, so let us now consider it in detailed.

Assume that a distant second planet has the mass of m and orbits the star (with the transiter) on a circular orbit with the radius of a. Such a planet should induce a similar circular motion on the host star on an orbit of the radius a^′ ≃ m/M_⋆a, where M_⋆ is the star mass. This reduces to a sinusoidal variation of the transits time delay with an amplitude of

\begin{equation} K_{\rm TTV} \simeq \frac{ma}{cM_\star } \sin i, \end{equation}

(11)

where c is the speed of light and i is the orbital inclination to the sky plane.

In the same case, the amplitude of the Doppler variation induced on the star is equal to

\begin{equation} K_{\rm RV} \simeq m\sin i \sqrt{\frac{G}{a M_\star }}, \end{equation}

(12)

where G is the gravitational constant. The ratio of these amplitudes now looks like

\begin{equation} \frac{K_{\rm TTV}}{K_{\rm RV}} \simeq \frac{a^{3/2}}{c\sqrt{G M_\star }} \simeq \frac{P}{2\pi c}, \end{equation}

(13)

where P is the orbital period of the distant companion, computed according the third Kepler law. We can see that the relative efficiency of the TTV method, in comparison with the RV one, only depends on the orbital period of the unseen companion.

Contrary to the RV and classic transit surveys, the TTV method might be useful to detect distant companions. Of course, in any case we must track the variation over at least a single period, so detection of a distant companion necessarily requires a long observation run. From this point of view, ground-based observatories should be more useful for TTV planet detection. Spacecraft rarely operate over a term longer than a few years, while ground-based observations can run on an indefinitely long time base.

The accuracy of the TTV data presented in this work is such that it would allow a robust detection of a signal, if its amplitude is above 1 min. With this TTV threshold, the formula (11) implies that a Jupiter-mass planet orbiting a solar-mass star could be detected, if its semimajor axis was at least 60 au, implying an extremely large period of at least ∼400 yr. This is not a realistic requirement for the observation time. Alternatively, the planet mass should be at least 12 times the Jupiter mass, if we want to detect it at a Jupiter-like orbit with a = 5.2 au. Also, the exoplanet detectability naturally increases for small-mass hosts. Thus, to robustly detect a Jupiter-mass planet, we should either wait for a long time or to decrease the TTV measurement errors to a level of seconds, roughly by an order of magnitude. Based on these computations, we may say that the TTV detection method is now in its early development stage, comparable to the early era of the Doppler technique (prior or near the 51 Pegasi b detection).

Summarizing, we can conclude that amateur-class TTV observations may occupy two possible niches: (i) detection of long-period massive exoplanets and (ii) detection of exoplanets trapped in mean motion resonances with known transiters (thus generating a more remarkable TTV signal via dynamical effects). However, the TTV detection of long-period Jupiter-like planets might be a feasible task for more advanced ground-based observatories that can achieve the TTV accuracy of a few seconds. Some of the data listed in Table 1 do provide such an accuracy. To enable the detection of Jupiter twins, such observations should be carried out in an experimental monitoring regime, i.e. on a regular basis and over a long term, similar to the modern Doppler surveys. Note that Jupiter analogues may represent a special interest, because such planets would point out exoplanetary systems that have architectures similar to Solar system, in which all giant planets are quite distant from the Sun. The chance to find an Earth twin in such a system might be higher.

At last, we note that the transit-fitting algorithm presented here is now implemented in the free planetpack package (Baluev 2013a) and is made available as of the current version planetpack 2.1.

This work was supported by the Russian Foundation for Basic Research (project no. 14-02-92615 KO_a), the UK Royal Society International Exchange grant IE140055, by the President of Russia grant for young scientists (no. MK-733.2014.2), by the programmes of the Presidium of Russian Academy of Sciences P21 and P22, by the Saint Petersburg State University research grant 6.37.341.2015, and by the Russian Ministry of Education and Science (contract no. 01201465056). OB acknowledges the support by the research fund of Ankara University (BAP) through the project 13B4240006. We express our gratitude to the anonymous reviewer for their useful suggestions.

†

Visiting Astronomer, Complejo Astronómico El Leoncito operated under agreement between the Consejo Nacional de Investigaciones Científicas y Técnicas de la República Argentina, and the National Universities of La Plata, Córdoba, and San Juan.

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

rn_t3l3.ttv

TERRA.rv

Supplementary Data.

Please note: Oxford University Press are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.