Mon.Not.R.Astron.Soc.000,1–16(2015) Printed9April2015 (MNLATEXstylefilev2.2) Benchmarking the power of amateur observatories for TTV exoplanets detection 5 1 Roman V. Baluev1,2∗, Evgenii N. Sokov1, Vakhit Sh. Shaidulin2,1, Iraida A. Sokova1, 0 Hugh R. A. Jones3, Mikko Tuomi3,4, Guillem Anglada-Escud´e3,5, Paul Benni6, 2 Carlos A. Colazo7, Matias E. Schneiter8, Carolina S. Villarreal D’Angelo8, r p A Artem Yu. Burdanov9, Eduardo Fern´andez-Laju´s10,11†, O¨zgu¨r Ba¸stu¨rk12, Veli-Pekka Hentunen13, and Stan Shadick14 8 1Central Astronomical Observatory at Pulkovo of Russian Academy of Sciences, Pulkovskoje shosse 65, St Petersburg 196140, Russia ] 2Sobolev Astronomical Institute, St Petersburg State University,Universitetskijprospekt 28, Petrodvorets, St Petersburg 198504, Russia M 3Universityof Hertfordshire, Centre for Astrophysics Research, Science and Technology Research Institute, College Lane, I Hatfield AL109AB, UK . 4Universityof Turku, Tuorla Observatory, Department of Physics and Astronomy, V¨aisa¨l¨antie 20, FI-21500, Piikki¨o, Finland h 5School of Physics and Astronomy, QueenMary Universityof London, 327 Mile End Road, London E1 4NS, UK p 6Acton Sky Portal (Private Observatory), Acton, MA, USA - o 7Observatorio Astrono´mico, Universidad Nacional de Co´rdoba, Laprida 854, Co´rdoba X5000BGR, Argentina r 8Instituto de Astronom´ıa Teor´ıca y Experimental, Universidad Nacional de Co´rdoba, Laprida 854, Co´rdoba X5000BGR, Argentina t 9Kourovka Astronomical Observatory of Ural Federal University,Mira str. 19, Ekaterinburg 620002, Russia s a 10Facultad de Ciencias Astrono´micas y Geof´ısicas - Universidad Nacional de La Plata, Paseo del Bosque S/N - 1900 La Plata, [ Argentina 11Instituto de Astrof´ısica de La Plata (CCTLa Plata - CONICET/UNLP), Argentina 2 12Ankara University,Faculty of Science, Department of Astronomy and Space Science, TR-06100, Tandogan, Ankara, Turkey v 13Taurus HillObservatory, Warkauden Kassiopeia ry., H¨ark¨am¨aentie 88, 79480 Kangaslampi, Finland 8 14Physics and Engineering Physics Department, Universityof Saskatchewan, 116 Science Place, Saskatoon, Saskatchewan, 4 S7N 5E2, Canada 7 6 0 Accepted 2015April8.Received2015April7;inoriginalform2015January27 . 1 0 5 ABSTRACT 1 We perform an analysis of ∼ 80000 photometric measurements for the following 10 : v starshostingtransitingplanets:WASP-2,-4,-5,-52,Kelt-1,CoRoT-2,XO-2,TrES-1, i HD 189733, GJ 436. Our analysis includes mainly transit lightcurves from the Exo- X planetTransitDatabase,publicphotometryfromtheliterature,andsomeproprietary r photometry privately supplied by other authors. Half of these lightcurves were ob- a tainedbyamateurs.Fromthisphotometrywederive306transittimingmeasurements, as well as improved planetary transit parameters. Additionally,for6ofthese10starswepresentasetofradialvelocitymeasurements obtained from the spectra stored in the HARPS, HARPS-N, and SOPHIE archives using the HARPS–TERRA pipeline. Our analysis of these TTV and RV data did not reveal significant hints of addi- tional 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 sec, although their parameters are model-dependent and uncertain, while radial velocities did not reveal statistically significant Doppler signals. Key words: planetary systems - techniques: photometric - techniques: radial veloc- ities - methods: data analysis - methods: statistical - surveys E-mail:[email protected] ∗ VisitingAstronomer,ComplejoAstron´omicoElLeoncitooper- gacionesCient´ıficasyT´ecnicasdelaRepu´blicaArgentinaandthe † ated under agreement between the Consejo Nacional de Investi- NationalUniversitiesofLaPlata,Co´rdobaandSanJuan. 2 Baluev et al. 1 INTRODUCTION Thisworkrepresentsanattempttodeterminetheprac- tical efficiency of such an approach, based on the pho- Thefirstextrasolarplanet,orbitingasolar-typestar51Pe- tometry data, taken mainly from the Exoplanet Tran- gasi, was discovered by Mayor & Queloz (1995), based on sit Database (ETD) of the Czech Astronomical Union, the precision Doppler measurements of the ELODIE spec- http://var2.astro.cz/ETD/. trograph.Afterthat,thenumberofthedetectedexoplanets Currently, about 30 observatories are regularly con- grew continuously, exceeding 1000 so far. In fact, 20 years tributingtothedatabase,includingamateuraswellaspro- ago a new rapidly-growing domain of fundamental science fessional ones. This network offers telescopes of different was created, devoted to theexoplanet research. apertures–from 20cm to2.6m.Theirlocations areshown Currently, most of the known exoplanetary candidates inFig.1.Wemustnotethatintermsoftechnicalcharacter- were detected by one of two major techniques: the radial istics of thetelescopes, thereis nosharp boundarybetween velocity (RV) method or transit method that increased its theamateurand professional equipment,andthequalityof output in recent time thanks to the launch of specialized the observations is often determined by the local astrocli- spacecraft CoRoT (ESA) and Kepler (NASA). mate, which can have even more important effect than the Unfortunately, both the Doppler and transit exoplanet telescope size. The observational programme involves cur- detection methods require sophisticated and expensive in- rently 20 transiter hosts that are more or less regularly ∼ strumentationandremainpracticallyinaccessibletothema- observed.Forsomeofthestars,severalyearsofobservations jorityofthecommunity.TheRVmethodrequirestheuseof arealready available. However,thetransit fittingalgorithm extremelystablespectrographs.Thereisonlyaboutadozen used by ETD is criticized for its simplicity and imperfec- ofsuchinstrumentsintheworld thatarecapableofdetect- tions (e.g. Petrucci et al. 2013). In this work we present a ing an exoplanet. The transit detection is also difficult. It new data reduction pipeline that handles subtle photomet- requires precise alignment of the planetary orbit with the riceffects(liketherednoise)moreaccurately.Thispipeline observer, and the transit event is rather short in time, al- was developed to deal with the photometry of a relatively though periodic. To detect a transiting planet, we have to poor or moderate quality, which is typical for the amateur observelotsofstars,andwenecessarilydealwithlargenum- data in ETD. bers of null detections. This requires the use of specialized Wedonotlimitourselvestoonlyamateurobservations ground-based robotic telescopes or telescope networks, or oronlyETD.Wealsousephotometricdatapublishedinthe space observatories like the above-mentioned CoRoT and literature. Moreover, for our targets we revealed a moder- Kepler. Organizing such a campaign is a difficult task even ateamountofthespectrastoredintheHARPS,HARPS-N, for professional scientific teams. andSOPHIEarchives,andwederiveRVdatafromthemto In terms of the photometric accuracy and quality, provide and independent “calibration” of our TTV results. ground based telescopes are no match to spacecrafts like However, in this work we have no goal to perform a fully Kepler. But all space projects have a common disadvan- self consistent TTV+RV analysis. Instead, we aim to char- tage: they are severely limited in time, while the sensitiv- acterize the accuracy and reliability of the TTV data that ity to weak planetary signals in the data degrades quickly we can derivefrom just the photometry. whentheplanetaryorbitalperiodexceedstheobservational The structure of the paper is as follows. In Sect. 2 we time base. This condition means that Kepler data cannot provide a detailed description of all the data that we in- reliably detect long-period planets, analogous to the giant cludeinouranalysis.InSect.3weintroducethealgorithms planets of Solar System. In general, the responsibility for usedto process thephotometric data.In Sect.4we present follow-up observations always returns to ground-based ob- the TTV data derived from the photometry and describe servatories.Ground-basedobservationsofplanetarytransits theresults of theirperiodogram analysis. InSect.5 wegive aremuchlessdemandingthanprecisionradialvelocitymea- theremainingfittedparametersoftheplanetsconsideredin surements, and of course much cheaper than projects like the work. In Sect. 6 we describe the RV data obtained for Kepler. In fact, such observations are possible with com- some of our targets on the base of the spectra found in the merciallyavailableequipmenttypicallyusedbyasignificant HARPS/HARPS-N/SOPHIEarchives.InSect.7wediscuss communityof“amateurs”.Thus,theexoplanetaryhuntcan indetailthecaseofWASP-4,forwhichwedetectedpossible potentiallybecome“citizenscience”thatcantriggeraqual- hintsof a weak TTV signal. itative leap in thefield. However, amateurs are definitely not equipped to un- dertake classic transit surveys like e.g. SuperWASP or oth- ers.Butnonethelesstheycanprovideausefulscientificcon- 2 THE SOURCE DATA tribution by means of the transit timing variation (TTV) method (Holman & Murray 2005; Agol et al. 2005). In this OurprimarysourceofthephotometricdataistheExoplanet approach,thereisalistofwelldefinedtargetsthatarecon- Transit Database (ETD). We extracted a set of the best tinuouslymonitored.Each target isaknownhost ofatran- transit lightcurvesfrom ETD for10selected stars.Also, we sitingplanet,anditstransitsareregularlyobserved.Ifmore addedin ouranalysis several publictransit lightcurvesthat planetsorbitthehost,theyshouldinduceperturbationalef- wefoundintheliterature.Thus,ourphotometricdatacover fects on the motion of the transiter, causing observable de- 10 targets with 306 transit lightcurves, and about 80000 lays to its transits, in comparison with a strictly periodic individualphotometricmeasurements.Thepublicdatawere ephemeris. The TTV method allows the observation time presentedintheworkslistedinTable1,andmostarestored to be used more efficiently, because we know what targets intheVizierdatabase.Amongthese306transitlightcurves, should beobserved, and when. 161 (or roughly half) were obtained by amateurs, while 65 Benchmarking amateur TTV exoplanets detection 3 Figure 1.WorlddistributionoftheobservatoriesregularlycontributingtotheETDdatabase. Table 1.Sources ofthephotometric data(except forETD). Table 2.ExplanationoftheRVdatafile target references column explanation WASP-2 Southworth etal.2010 1 observationtime(BJD) WASP-4 Wilsonetal. 2008; Gillonetal. 2009; 2,3 derivedradialvelocityanditsuncertainty Winnetal. 2009; Southworth etal. 2009b; 4 astandardizednameoftheinputRVfilethatincludes Sanchis-Ojedaetal.2011;Nikolovetal.2012; thetargetnameandthenameofthespectrograph Petrucci etal.2013 WASP-5 Southworth etal.2009a WASP-52 No be treated with the same care. However this photometry Kelt-1 Siverdetal.2012 is rather accurate in itself, and thus it is still useful in CoRoT-2 Gillonetal.2010 constraining all transit parameters except for mid-times XO-2 Fernandez etal.2009 (e.g. by adding a fittable time offset to these light curves TrES-1 Winnetal.2007b model relative to the other lightcurves). All timings in the HD189733 Bakos etal. 2006; Winnetal. 2007a; Pontetal.2007 photometric data were transformed to the BJDTDB time GJ436 Beanetal. 2008; Shporeretal. 2009; stamps by means of the public IDL software developed by Stevenson etal.2012 Eastman et al. (2010). Additionally,weusedtheprecisionradialvelocity(RV) obtained from the spectra of the HARPS archive, avail- arefromprofessionalobservatoriescontributedtoETD,and able for the following targets from our photometry sample: 80 were published in the listed literature. WASP-2, -4, -5, HD 189733, Corot-2. These spectra were We need to note that the data from processedwiththeadvancedHARPS–TERRApipelinethat (Sanchis-Ojedaet al. 2011) also contain reprocessed offers an improved RV accuracy (Anglada-Escud´e& Butler photometry from (Winn et al. 2009), so the original 2012).ThequalityoftheseRVdata,andtheirnumberpera Winn et al. (2009) data were not actually included in our target,appearednotveryhigh.Themeaningofthecolumns analysis. Also, we eventually decided to drop the HST in the attached RV data file is given in Table 2. data by Bean et al. (2008), because they all appeared to For GJ 436 we found a large amount of high- cover only small parts of a transit and could not produce quality HARPS RV data, similar to the one considered by good precision in the derived mid-times. The photometry Lanotte et al.(2014),aswellasalargeKeckRVtimeseries presented by Petrucci et al. (2013) is not public, but it recently published by Knutson et al. (2014). We acknowl- was kindly released to us by the authors. The WASP-4 edgethatcaseslikethisdeservetobeinvestigateddetailedly data from (Southworth et al. 2009b) appeared not very in a separate work. In this paper we only present our TTV reliable due to the clock errors of the Danish telescope and TERRA RV data for GJ436, processed in the common (Nikolov et al. 2012; Petrucci et al. 2013), and we believe simplistic way as for the otherstars. that the data from (Southworth et al. 2009a, 2010) should In this work we do not perform the joint transits+RV 4 Baluev et al. fits, because it appeared that to perform such an analysis The trajectory curvature also depends on the orbital at a desirable level of quality, we must provide a solution eccentricity, but usually this eccentricity is difficult to de- to several non-trivial issues that fall outside the scope of termine from transit observations reliably, and we have no this paper. These issues include, e.g., the treatment of the otheroption except to assume that it is zero. A zero eccen- correlational structure of the RV noise, and treatment of tricity is a good prior assumption for most of these short- the RV points affected by the Rossiter-McLaughlin effect. period planets. In the cases when radial velocity data are InthisworkweusedtheRVdatamainlyasanindependent also available, the accurate eccentricity information can be source of information. obtainedfrom thejoint transit+RVfits,butwedonot per- form fits of such typein thiswork. After the function δ(t) is calculated, we approximate the relative flux reduction with the use of the stellar limb 3 METHODS OF THE TRANSIT darkening model by Abubekerov& Gostev (2013). These LIGHTCURVES ANALYSIS authors provided theoretical formulae for various types of 3.1 Details of the transit model and its thelimbdarkeningeffect,aswellasasoftware library writ- parametrization ten in C. The library provides subroutines to compute the observedlightfluxreduction∆Lasafunctionoftheeclips- Let us first adopt the approximation that the transiting ing planet radius r and of the projected separation δ (as- planet moves along a straight line with constant velocity, suming that the star radius is unit), as well as of the limb thusneglectingtheorbitalcurvatureofitsactualtrajectory darkening coefficients that depend on the selected model. duringthetransit.Themodelofthelinearmotioneasilypre- Besides,thislibraryprovidespartialderivativesof∆Lwith dicts the separation between the centers of planet and star regard to its arguments. These derivatives are necessary to disks,δ,as afunction oftime. Wehave4kinematiccharac- compute the gradient of the likelihood function, which is teristics of the transit that must be fitted: (i) the mid-time also used by thetransit fitter. of the transit t , (ii) the duration of the transit t , defined c d Thus for a transit lightcurve we have 4 fittable param- asthetimespentbetween thefirst and fourthcontacts and eters of the planet or planetary orbit: the planet/star radii (iii) the impact parameter b, measuring the smallest pro- ratior,thetransit mid-timet ,thetransitdurationt ,and jected separation δ, and (iv) the projected planet radius r c d the impact parameter replacer p. As the orbital period P thatsimultaneouslydeterminesthetransitdepthandthege- is known for all our planets with a very good accuracy, we ometryoftheingress/egressphases.Giventheseparameters treat it as a fixedparameter in the formula (4). and easy geometric constructions, the projected separation Additionally,therearetwoparametersdeterminingthe δ can be expressed as: limbdarkeningmodel.Inthisworkweuseaquadratictwo- δ(t)= b2+[(1+r)2 b2]τ2, τ =2t−tc, b 61+r.(1) term model of the stellar limb darkening with two coeffi- − td | | cientstobedetermined,AandB.Thebrightnessofapoint This aspsumes that the star radius is unity. To simplify the on the visible stellar disc, observed at a given separation formula for δ and simultaneously get rid of the non-trivial from its center, ρ, is modelled as definition domain for b, our algorithm adopts internally a replacement of the impact parameter b: I(ρ) = 1 A(1 µ) B(1 µ)2 = − − − − p = 1 A 2B+(A+2B)µ+Bρ2, b= (1+r). (2) − − 1+p2 µ= 1 ρ2. (5) − Withpthenewparameterpwedonotneedtoworryaboutits p However,inpracticethebrightnessmodel(5)caneasily domain: any real value of p is physically meaningful. Thus turnnon-physical,ifthecoefficientsAandB areallowed to we eliminate the danger that this parameter can walk to a attain arbitrary values. This becomes a significant problem forbiddendomainduringthefit.Withpreplacingb,wehave if our data are polluted by some systematic errors that are always difficult to foresee in advance. There are a couple of τ2+p2 δ(t)=(1+r) . (3) naturalbasicconstraintsthatweplaceonthecoefficientsA s 1+p2 and B to keep the model (5) physically reasonable: However, this approximation of δ(t) might be inaccu- dI ρ rate due to the curvature of the planet trajectory. In this I(ρ)>0, =[2Bµ (A+2B)] 60, ρ [0,1].(6) workweincludeacorrection tothisformula,assumingthat dρ − µ ∀ ∈ the planet moves along a circular orbit. In this approxima- If the second condition (the one on the derivative) is sat- tion, the projected distance can be expressed by the same formulae (1) and (3), substituting the following quantity τ′ isfied, I(ρ) is a monotonic non-increasing function, so the condition I(ρ) > 0 for ρ [0,1] is then equivalent to instead of τ: ∈ I(1) = 1 A B > 0. In the second condition, the ex- τ′= sinτα, α=πtd 1, (4) pression in−par−enthesis is a linear function of µ, so to have sinα P ≪ it always non-positive for µ [0,1], it is necessary and suf- ∈ whereP isthetransiter’sorbitalperiod.Theauxiliaryangle ficienttohaveit non-positiveattheboundariesµ=0(disk α reflects the curvature of the circular orbit. In fact, this limb)andµ=1(diskcenter).Finally,wehavetotalofthree correction induced only a well negligible effect on our TTV elementary inequality constraints on A and B to satisfy: measurements,butnonethelessourresultscorrespondtothe model with this correction included. A+B61, A+2B>0, A>0. (7) Benchmarking amateur TTV exoplanets detection 5 Note that the coefficient B is allowed to be negative here individualtransits separately from each other. Wecan note (but B > 1). We do not put an extra condition B > 0. thefollowing potential weaknesses of this approach: − Negative values of B mean that the limb-darkening gradi- (i) It is not taken into account that the limb darkening ent is diminishing (in absolute value) closer to the limb, coefficients depend on the photometry band, which are dif- although it never turns positive (any “limb brightening” is ferentfordifferentlightcurves.Therefore,theresultsofsuch disallowed). fitting method would refer to some averaged limb darken- We may satisfy (7) by means of making a smooth re- ing, possibly introducingminor modelling errors in individ- placementoftheparametersAandB,suchthattheformu- ual lightcurves. From theother side, the limb darkening ef- lae of the replacement would disallow the conditions (7) to fect in the transit curveis always symmetricrelative to the bebroken.This can bereached,for example,bythefollow- mid-time,implyingthatitsimpactonthederivedmid-times ing trigonometric replacement: should be small, if not negligible. In fact, we noticed that A=sin2θ(1 cosϕ), B=sin2θcosϕ. (8) even fitting of a transit model without any limb darkening − at all does not change the derived TTV data significantly Naturally,whateverrealvaluesthenewauxiliaryparameters (beyond the estimated parametric uncertainties). However, θandϕmightattain,theresultingvaluesofAandBalways inthefinalversionofouralgorithm wedecidedtodisentan- satisfy (7). From the other side, each point (A,B) in the gle the limb darkening coefficients for the best lightcurves domain (7) mapsto some pair (θ,ϕ) with real values of the (thosethathaver.m.s.smallerthan10percentofthetransit parameters: depth)and fit thesecoefficients independently. B (ii) It is not taken into account that transit duration sin2θ=A+B, cosϕ= . (9) A+B mayalsobesubject tovariations,likethetransit mid-time. However, in practice, it appeared that the accuracy of the Note that from (7) it follows that A+B > B, meaning | | transit duration estimations was roughly an order of mag- that A+B is nevernegative. nitude worse than those of the mid-times. Moreover, some Therefore, treating theauxiliary angles θ and ϕ as pri- lightcurves cover the transit event only partially, implying mary fittable parameters, we can satisfy the conditions (7) that its duration would remain almost unconstrained when automatically. The result of the fitting would correspond fitted independently. We do not address transit duration to either an internal point of the domain (7), or to some variations (TDV) in this work, focusing our attention on point on its boundary, if the actual data suggest to move theTTV. thesolutiontoanon-physicaldomainduetoe.g.theirpoor (iii) Themid-timeestimates, obtained insuchamanner, statisticalqualityorsystematicerrors.Theboundarypoints arenotnecessarilyuncorrelated.Thesemid-timesarecorre- correspond to certain special values of θ and ϕ that can be latedwiththeremainingtransitparameters, whicharenow easilyidentified:θ= π/2+πkonthefirstboundaryof(7), ± shared between the lightcurves. Through this effect, some ϕ = π+2πk on the second boundary, and ϕ = 2πk at the cross-correlation between the estimated mid-times may ap- third boundary. pear. This creates a risk of unexpected statistical effects in the derived TTV data, like e.g. the non-white noise. How- ever, the magnitude of these TTV correlations usually re- 3.2 Details of the fitting procedure mainsverysmall(maximumofafewpercent,and 0.1per ∼ centinaverage),andnodeviationsfromthewhitenoiseare InouranalysisweareinterestedinobtainingtheTTVdevi- seenintheperiodogramsoftheTTVdata.Moresignificant ationswithmaximumaccuracyachievablewiththeavailable (>10percent)TTVcorrelationscanappearforlightcurves photometric data. However, the ETD data, even after the that offer only a partial coverage of the transit, but such pre-selection,oftenhaveonlyamoderatequality,performing lightcurvesrepresent only a minor fraction of our data. acompleteandindependentfitforeachtransitlightcurveis notagoodoption.Theaccuracyofthusobtainedtransitpa- Inadditiontotheplanetarytransititself,ourlightcurve rameterswouldoftenbepoorlyconstrained,andthiswould modelsincludeapolynomialtrendwithfittablecoefficients. impact the accuracy of the fitted TTV data as well. More- Suchatrendisnecessary totakeintoaccountvariousdrift- over,thetransitfitsforindividuallightcurvessometimesdo ing effects, e.g. the effect of airmass or other types of sys- notevenconvergetoareasonableresult,e.g.,duetoincom- tematicvariationsthat appearfrequentlyinourdata.Each pletecoverageofthetransitorpresenceofsignificantcurved transit lightcurve has an individual fittable trend with a trends. separate set of trend coefficients. After some experiment- To overcome these issues, we adopt in this work the ing with the data, we found that cubic trends represent a following approach. We perform a joint fit of all transit goodcompromisebetweenthemodeladequacyanditspara- lightcurves,availableforagivenstar,assumingthatmostof metriccomplexity.Wedidnottrytoreducethissystematic the transit parameters, except for the mid-times, are equal photometric variation based on its correlation with the air- for different lightcurves. Such a “shared” transit parameter massfunction:whileforsomelightcurvessuchacorrelation is still fittable,takingintoaccount theconstraint ofits val- looked clear, for others it was not obvious, indicating that uesbeingequalbetweendifferentlightcurves.Themid-times othersystematic effects were in thegame. are fitted individually for each lightcurve, i.e. they remain The lightcurves model is relatively complicated, while unconstrained. thephotometryusedinthisworkisnotofaverygoodqual- Such an approach uses the full statistical power of all ity. In practice we often faced difficulties causing the fitter photometricdata,availableforagivenstar,tofittheshape tobetrappedinthelocalmaximaofthelikelihoodfunction, of the transit curve, while still fitting the TTV offsets of which was related to an unrealistic branch of the solutions. 6 Baluev et al. For example, without special care, see below, we frequently and t , is modelled as j otrbatnasinitedwinthonr-phy1s.icAallssoo,lusotimoneslicgohrtrceusrpvoensddinogntootacogvrearztinhge Vij =Cov(xi,xj)=σi2,wht(pwht)δij+Vij,red(pred,τ), cmoimd-ptliemteetirsacnos≫mit,palincdatiendsbuychitcsassetrsotnhgecfiotrtrienlgatoifonthweittrhantshiet Vij,red =predRij(τ), Rij =ρ ti−τ tj , (10) (cid:18) (cid:19) coefficientsofthepolynomialtrend.Toavoidsuchtraps,we with pwht, pred, and τ being free fittable parameters. Here, worked out the following sequence of auxiliary preliminary the first noise term, σ2 , represents the white fraction i,wht fits: of the noise, detailed below. The red noise is given by the second term, V , which also depends on the function ρ(t) (i) Perform a preliminary fit constraining the mid-times ij that represents an adopted shape of the correlation func- at a regular grid (with free scale and offset), implying all tion. We use mainly the exponential correlation function TTV residuals are zero by definition; fixing the impact pa- ρ(t) = exp( t). This not a unique choice, but in practice rameter p at an intermediary value of 1 (b 1/√2), and −| | ≈ thismodel of thered noise usually appears adequate. fixing the limb darkening coefficients at A = B = 0.25 (or The white-noise term σ2 in (10) requires a separate θ=π/4 and ϕ=π/3). i,wht discussion.Webelievethatphysicallytheso-called additive (ii) Refitafterreleasingthemid-timesandimpactparam- modelwouldbemoresuitablehere.Inthismodel,thevalue eter, but still holding thelimb darkeningcoefficients fixed. σ2 isdeterminedasasumofthestatedinstrumentalvari- (iii) Refit after releasing the limb darkening coefficients i,wht anceandofanunknownfittable“jitter”variance,generated (bindingthem across different lightcurves). by e.g. Earth atmospheric instability, or by unassessed in- (iv) Refit after full release of the limb darkening coeffi- strumental instability, or by short-term (minutes to hours) cients for the best lightcurves (those with r.m.s. < 0.1 of intrinsic variations of the stellar flux, whenever these vari- thetransit depth). ations can be treated as a white noise. Although, we must note that in practice the instrumental uncertainties in our After that, we also apply the red-noise detection and data do not look very trustable, and often they are just fittingprocedure as described in thesection below. omitted, forcing us to assume that the measurements have just equal uncertainties. Therefore, in any case we should notexpect thatourwhite noisecan beaccurately modelled from thephysicalpointofview.Inthesecircumstances, the 3.3 Reduction of the red noise noise model should be chosen mainly on the basis of math- It is already known well that stellar photometric data ematical simplicity or usefulness, relying on only minimal usually include a correlated (“red”) noise component that or no knowledge of the underlying physics. We decided to has an important effect on the fitted transit parameters useinourworkaregularizednoisemodeldefinedin(Baluev (Pont et al. 2006; Winn et al. 2008; Carter & Winn 2009). 2015). In the majority of practical cases, this model should Aneasytechniquetocalculatetheeffectoftherednoiseon be equivalent to the additive noise model. This regularized the fitting uncertainties was introduced in these works. Al- modelprovedratherresistantwithrespecttovariouspitfalls though these authors avoid making restrictive assumptions appearing during thefittingprocedure. concerning the correlation structure of the red noise, their The fitting of the compound model, involving simulta- method still remains rather simplistic and it does not take neously the models of the transit curve and of the photo- intoaccountimportanteffects.Theymainlyfocusonamore metric noise, is done by means of the maximum-likelihood accuratedeterminationoftheuncertaintiesinthetransitfit, approachwithdetailsgivenin(Baluev2015).Thatworkwas which are typically underestimated without a proper treat- devotedprimarilytotheanalysisoftheradialvelocitydata, ment of the red noise. But the red noise may also induce butmathematically themethodsthatweareusinghereare biases in the fitted parameters themselves. This biasing ef- identical. fect was already noted when processing Doppler data af- However, after an attempt to apply this technique in fected by red noise in a generally simialar way (e.g. Baluev practice literally, we faced the problem that the red noise 2011, 2013b). Also, there is not an obvious way to control could not be fitted in many of our lightcurves. Only 1/3 to thevalidityoftheunderlyingassumptionsinthetraditional 1/2 of the lightcurves allowed for a reliable red noise fit, approach of the photometric red noise reduction. This ap- while in the other cases, the red noise was either not de- proach does not offer a complete noise model that could be tectable, or its estimated magnitude became pretty small verified for accuracy. Additionally, the original method by (in comparison with the uncertainty). Such cases lead to Pont et al. (2006) relies on the out-of-transit photometric modeldegeneraciesanddecreasethereliability oftheentire data, which are very limited in ourcase. fit for a given star. Thus, we decided to add the red noise In this work we apply the approach based on a para- termtothephotometricmodelonlyinthecaseswhenitwas metric modelling of the noise covariance matrix. This is an justified. adaptation of the red-noise fitting technique from (Baluev Weaddedtoouranalysispipelineasetofauxiliaryfits 2011, 2013b, 2015) that was developed to handle exoplan- in order to identify the lightcurves, in which the red noise etary radial-velocity data, in which the noise correlations couldbemodelled moreorlessreliably.Namely,wetriedto mayappeare.g.duetothestellaractivity.Inthisapproach performaseriesoftestfits,latterlyaddingarednoiseterm weapproximatetherednoisebyastationaryGaussianran- tothemodelofeachtransitlightcurve.Iftherednoisecould domprocesswithacovariancefunctionofagivenfunctional not be estimated reliably, more accurately if the relative form.Thecovariancecoefficientbetweentwoarbitrarypho- uncertainty of its estimated magnitude exceeded 2/3, the tometric observations x and x , acquired at the times t relevantrednoisetermwasremovedfromthemodelandthis i j i Benchmarking amateur TTV exoplanets detection 7 lightcurvewasfurthertreated ashavingpurelywhitenoise. We consider the entire set of 80000 normalized residuals ∼ Otherwise, weproceeded tothenexttest fitpreservingthis for all stars involved in our analysis, and plot the relevant red noise term in themodel. In theend,after all individual distribution in Fig. 3 lightcurveswereprocessed in suchaway,weperformed one As we can see, the empirical distribution is indeed al- more check of each red noise term (because some of them most Gaussian, except for the very tails (> 4-sigma devia- could turn insignificant after the other terms being added) tion).Theempiricaldistributionhasheaviertails.However, and the finalfit was made. by removing data in the distribution tails we can reach a Thequalityofsuchanapproachtotherednoisereduc- good agreement with the Gaussian approximation. In fact, tion isexamined inFig.2.Thebest waytoobservea“non- these extreme points represent outliers and should be re- white”noiseistolookatitspowerspectrum:thewhitenoise moved in any case. The final results below correspond to should have constant power in average, while the red noise thedatawiththese24outliersremoved.Wealsonoticethat should demonstrate a systematic uprise to long periods. In theminorsystematicdeviationsstillremaininginthedistri- our case we compute the so-called residual periodograms bution tails are mainly due to Bakos et al. (2006) data for (see Baluev 2015) for each individual lightcurve, and then HD 189733. Visual investigation reveals that some of these compute their average per each star involved in the analy- lightcurves contain limited segments with unexpectedly in- sis. These averaged periodograms are plotted in Fig. 2. For creased photometric scatter. Removal of the Bakos et al. each star, we sequentially try three base photometric mod- (2006) data makes the agreement of the residuals distribu- els: (i) just transits without any limb darkening plus the tion with the Gaussian onevery good. whitenoise;(ii)sameasmodel(i)pluscubictrendsandthe effect of the quadratic limb darkening; (iii) same as model (ii) plusthe red noise term. 4 DERIVED TTV DATA AND THEIR We can clearly see from Fig. 2 that thelong term pho- ANALYSIS tometricvariationsarenotreducedtojusttrends.Although thetrendsthemselvesarealsoimportant,aftertheirremoval Usingthefittingtechniqedescribedabove,wecomputedthe the power spectra still contain an additional power in the mid-times for each observed transit, and placed these TTV period range of > 10 min. The trends can only affect the measurements and other accompanying data in the online period domain of >200 300 min (this corresponds to the supplement to the article. The supplement represents a file − typicaltimespan ofthelightcurves).Inthemost cases, ap- containingatexttable.Themeaningofthecolumnsinthis plication of our red-noise model remarkably suppresses the file is described in Table 3. This file only contains the pa- remainingexcessivepower,atleastbyafactorof2.Wemust rameters that are attached to individual lightcurves. The acknowledgethatthereductionofthisrednoiselooksrather common fit parameters that are shared between different difficult and still not entirely complete. lightcurvesare given separately in Sect. 5 below. Some statistical information on the derived red noise TheTTVresiduals forthesedataareplottedinFig. 4. characteristicsnowfollows.About60percentoflightcurves In Fig. 5 we verify the normality of the distribution of the did not contain any detectable red noise at all. For the re- TTV noise present in these derived data in the way simi- maining 40 per cent of lightcurves, the red noise caused an lar to Sect. 3.4. In the TTV residuals we find no obvious increase of the TTV uncertainty by a factor of 1.3 on av- outliers or deviations from the Gaussian distribution. The erage (a median value). This is not very large, although a apparentexcessin thedistribution tails that can beseen in minorfractionofthemid-timesdemonstratedanincreasein theleft panelofFig.5isuncertainduetoarelatively small theuncertaintyuptoafactorof2 3.Concerningthebias, number of TTV data points. Nonetheless, we identified 9 − appearinginthederivedmid-timesthemselves(ratherthan possible “candidate outliers”, and tried to remove them to- uncertainties),itsaveragevaluewasonlyabout1/6fraction getherwith7TTVpointsreferringtotheSouthworth et al. of the uncertainty, and only a few points demonstrated a (2009a,b, 2010) data. The TTV analysis discussed below shift above 1σ. was performed for the full TTV data set as well as for the Summarizing, in themajority of thecases, either there reduced one. was not any detectable red noise, or the change in the First, we made an effort to detect possible long-term TTV data was small, even if the red noise term was de- curvatureinourTTV.Thiswasdonebymeansofmodelling tectable. Nonetheless, it is important that some individual the mid-times by a quadratic function (instead of only a TTV points are affected much more. linear function that would correspond to strictly periodic transits).NoneoftheTTVtimeseriesdemonstratedalong- term quadratic variation with at least 2-sigma significance 3.4 Testing the Gaussianity of the noise and (themost suspicious cases had 1.5σ to 1.7σ). clearing away the outliers Second, we undertook a periodorgam analysis of the Inthedata-analysismethodsdescribedabove,welargelyre- TTV residuals, in order to reveal possible periodic pertur- lied on the assumption that the photometric noise is Gaus- bationsfromanyadditionalunknownplanetsorbitingthese sian.Thisassumptionneedssomeverification,whichweper- stars. We again applied “residual periodograms” similar to formed in the following manner. After fitting all the data, the ones that were used to plot Fig. 2 above. Now our null we computed the residuals and normalized them by their model included a quadratic trend that is related to the pe- relevant modelled uncertainty (square root of the sum of riod of the main (transiting) planet and white noise with the estimated white and red noise variances). If the noise a fittable additive “jitter”. The alternative model also in- wasGaussian,thedistributionofthesenormalizedresiduals cluded a sinusoidal signal at a probe frequency. None of should be close to the standard Gaussian, and vice versa. these periodograms demonstrated any signs of a non-white 8 Baluev et al. Averaged likelihood-ratio residual periodograms for the WASP-5 lightcurves 20 6 6 no trends 5 cubic trends 5 cubic trends wer 15 nwoh iltiem nbo disaerk oennliyng wer qwuhaitder antoicis leim obn ldyarkening wer qwuhaitder+arteicd linmobis edarkening po po 4 po 4 R R R L L L g- 10 g- 3 g- 3 o o o e l e l e l ag ag 2 ag 2 ver 5 ver ver a a a 1 1 0 0 0 1 10 100 1000 1 10 100 1000 1 10 100 1000 period [min] period [min] period [min] Averaged likelihood-ratio residual periodograms for the WASP-52 lightcurves 25 5 5 no trends cubic trends cubic trends wer 20 nwoh iltiem nbo disaerk oennliyng wer 4 qwuhaitder antoicis leim obn ldyarkening wer 4 qwuhaitder+arteicd linmobis edarkening o o o p p p R 15 R 3 R 3 L L L g- g- g- o o o e l 10 e l 2 e l 2 g g g a a a er er er v v v a 5 a 1 a 1 0 0 0 1 10 100 1000 1 10 100 1000 1 10 100 1000 period [min] period [min] period [min] Averaged likelihood-ratio residual periodograms for the Kelt-1 lightcurves 20 5 5 no trends cubic trends cubic trends wer 15 nwoh iltiem nbo disaerk oennliyng wer 4 qwuhaitder antoicis leim obn ldyarkening wer 4 qwuhaitder+arteicd linmobis edarkening o o o p p p R R 3 R 3 L L L g- 10 g- g- o o o e l e l 2 e l 2 g g g a a a ver 5 ver ver a a 1 a 1 0 0 0 1 10 100 1000 1 10 100 1000 1 10 100 1000 period [min] period [min] period [min] Averaged likelihood-ratio residual periodograms for the Corot-2 lightcurves 35 5 5 30 no trends cubic trends cubic trends wer 25 nwoh iltiem nbo disaerk oennliyng wer 4 qwuhaitder antoicis leim obn ldyarkening wer 4 qwuhaitder+arteicd linmobis edarkening o o o p p p LR 20 LR 3 LR 3 g- g- g- e lo 15 e lo 2 e lo 2 g g g a a a er 10 er er v v v a a 1 a 1 5 0 0 0 1 10 100 1000 1 10 100 1000 1 10 100 1000 period [min] period [min] period [min] Averaged likelihood-ratio residual periodograms for the GJ436 lightcurves 50 5 5 no trends cubic trends cubic trends wer 40 nwoh iltiem nbo disaerk oennliyng wer 4 qwuhaitder antoicis leim obn ldyarkening wer 4 qwuhaitder+arteicd linmobis edarkening o o o p p p R 30 R 3 R 3 L L L g- g- g- o o o e l 20 e l 2 e l 2 g g g a a a er er er v v v a 10 a 1 a 1 0 0 0 1 10 100 1000 1 10 100 1000 1 10 100 1000 period [min] period [min] period [min] Figure 2. Examining the effect of the red noise in the lightcurve photometry for a few of stars belonging to our sample. Please note thatplotsinthefirstcolumnhavedifferentordinatescalethanthoseintheothertwo.SeeSect. 3.3foradetaileddiscussion. Benchmarking amateur TTV exoplanets detection 9 raw data (78153 points) 24 outliers removed 10 10 e 5 e 5 til til n n a a u u q q 0 0 n n a a si si s s u u a -5 a -5 G G -10 -10 -10 -5 0 5 10 -10 -5 0 5 10 normalized residual normalized residual Figure 3.Testingthe normalityof theentire 80000 photometric data usedinthe work.Abscissashows thevalue ofthe bestfitting ∼ residual, normalized by its total modelled variance (for white+red noise). The ordinate contains the normal quantiles of the empirical cumulativedistributionofthesenormalizedresiduals.PerfectlyGaussianresidualsshouldalllieclosetothemaindiagonalthatrepresents astandardnormaldistribution,whileadeviationfromthediagonalindicatesanon-Gaussiannoise. noise, and in fact almost all appeared consistent with the of this relation. The value of tref is not unique due to the c pure noise model, although the cases revealing significant periodic nature of the transits. We chose such a reference excessive TTV “jitter” were rather frequent. No significant transit, for which the uncertainty of t would be minimum. c TTVperiodicitieswerefoundinthedataforanyofthestars This condition simultaneously implies that the correlation involved in theanalysis, except for theWASP-4 case which betweentheestimatedP andtref iszero.Notethatitisnot c isdiscussedbelow.Theperiodogramsignificancelevelswere required to haveactual transit observation at this tref. c estimated using theapproach of Baluev (2008, 2009). The χ2 values in Table 4 always exceed unity, so TTV Webelievethat this result indicates aconsiderable im- the fit uncertainties might be moderately underestimated. provementinthestatisticalqualityoftheTTVdataand/or In some cases this higher than expected scatter could be methods of the statistical analysis, because for the original partlyexplainedbysmallnumberoftransitsinvolvedinthe TTV data from ETD we frequently detected spurious pe- analysis,andinsomepartitmightbeduetoadditionalun- riodicities possibly appearing due to imperfections of the seen planetsthatmayinducecomplicated TTV signals and transit fittingalgorithm. aredifficulttodetect.Anotherexplanationisincompletere- ductionoftherednoiseandtheeffectofnonlinearityofthe transit model. The latter effect can introduce biases both in the estimated parameters as well as in their uncertainty 5 IMPROVED TRANSIT CURVE estimations.Sofarwecouldnotdecidewhichexplanationis PARAMETERS morelikely.InterestinglythisexcessiveTTVscatterpersists In addition to the TTV data, we also provide the re- and is important even in such a robustly fittable case like maining best fitting parameters of the exoplanetary tran- HD189733, in which the model nonlinearity should bewell sit lightcurves. These are the transit parameters that were suppressed thanks to a large number of available observa- shared between different lightcurves, and they are given in tions and lightcurves. Table4.Inadditiontothepreviouslymentionedtransitpa- rametersr,t ,b,andthelimb darkeningparameters Aand d B, this table also contains the following data: number of the transit lightcurves used in the analysis, number of the lightcurvesinwhicharednoisetermwasrobustlydetected, theχ2 (weighted r.m.s.) oftheTTV residuals (relatively to thebestlinearfitofthederivedmid-times),andtherefined parameters of the transit emphemeris: transiter period P and a reference mid-time tref. The last two quantities were c obtained by requiring the mid-times to obey a strict linear relation with the transit count and fitting the coefficients 10 Baluev et al. 2008 2009 2010 2011 2012 2013 2014 2008 2009 2010 2011 2012 2013 2014 8 8 6 WASP-2 6 WASP-4 n] 4 n] mi mi 4 al [ 2 al [ du 0 du 2 si si V re -2 V re 0 TT -4 TT -2 -6 -8 -4 0 200 400 600 800 1000 1200 0 500 1000 1500 transit counter transit counter 2009 2010 2011 2012 2013 2014 2013 2014 4 8 3 WASP-5 6 WASP-52 n] 2 n] 4 mi mi al [ 1 al [ 2 du 0 du 0 si si V re -1 V re -2 TT -2 TT -4 -3 -6 -4 -8 0 200 400 600 800 1000 1200 1400 -50 0 50 100 150 200 250 300 350 400 transit counter transit counter 2012 2013 2014 2009 2010 2011 2012 2013 2014 20 4 15 Kelt-1 3 Corot-2 n] n] 2 mi 10 mi al [ al [ 1 du 5 du 0 si si V re 0 V re -1 TT TT -2 -5 -3 -10 -4 0 100 200 300 400 500 600 700 800 900 0 200 400 600 800 1000 1200 transit counter transit counter 2008 2009 2010 2011 2012 2013 2014 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 6 5 4 4 XO-2 3 TrES-1 min] 2 min] 2 ual [ 0 ual [ 01 TV resid --42 TV resid ---321 T T -4 -6 -5 -8 -6 0 200 400 600 800 0 200 400 600 800 1000 1200 transit counter transit counter 2006 2007 2008 2009 2010 2011 2012 2013 2014 2007 2008 2009 2010 2011 2012 2013 2014 6 12 4 HD189733 10 GJ436 n] 2 n] 8 mi mi 6 dual [ - 20 dual [ 4 si si 2 V re -4 V re 0 TT -6 TT -2 -8 -4 -10 -6 0 200 400 600 800 1000 1200 1400 0 200 400 600 800 1000 transit counter transit counter Figure 4. Graphs of the TTV residuals, based on the best fitting transit ephemeris (values of P and tref) from Table 4. The transits c forwhichasignificantrednoisewasdetected inthephotometry arelabelledwithacircle.