Table Of ContentMon.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.
