The TAOS Project: Results From Seven Years of Survey Data

0.46 MB·English

by  Z. -W. Zhang

Submittedto AstronomicalJournal
PreprinttypesetusingLATEXstyleemulateapjv.08/22/09
THE TAOS PROJECT:RESULTS FROM SEVEN YEARS OF SURVEY DATA
Z.-W. Zhang1, M. J. Lehner1,2,3, J.-H. Wang1, C.-Y. Wen1, S.-Y. Wang1, S.-K. King1, A´. P. Granados4,
C. Alcock3, T. Axelrod5, F. B. Bianco6, Y.-I. Byun7, W. P. Chen8, N. K. Coehlo9, K. H. Cook1, I. de Pater10,
D.-W. Kim11, T. Lee1, J. J. Lissauer12, S. L. Marshall13, P. Protopapas14,3, J. A. Rice9, and M. E. Schwamb15,16
Submitted to Astronomical Journal 2013 January 16
ABSTRACT
3 The Taiwanese-American Occultation Survey (TAOS) aims to detect serendipitous occultations of
1 stars by small ( 1 km diameter) objects in the Kuiper Belt and beyond. Such events are very rare
0 (<10−3 events p∼er starper year)andshortin duration( 200ms), so manystars mustbe monitored
2 at a high readout cadence. TAOS monitors typically ∼500 stars simultaneously at a 5 Hz readout
cadencewithfour telescopeslocatedatLulin Observato∼ryincentralTaiwan. Inthis paper,we report
n
a the results ofthe searchforsmallKuiper BeltObjects(KBOs)insevenyearsofdata. No occultation
J events were found, resulting in a 95% c.l. upper limit on the slope of the faint end of the KBO size
distribution ofq =3.34to 3.82,depending on the surfacedensity atthe break in the size distribution
5
2 at a diameter of about 90 km.
Subject headings: Occultations, Kuiper belt: general, Comets: general, Planets and satellites: forma-

] tion
P
E
1. INTRODUCTION form
.
h The size distribution of Kuiper Belt Objects has dn D−4.5 (1)
p dD ∝
been accurately measured down to diameters of
-
o D > 30 km (Fuentes et al. 2009; Fuentes & Holman down to a diameter of D 90 km. A clear break
≈
r 2008; Fraser & Kavelaars 2008; Fraser et al. 2008; in the size distribution has been detected near
t
s Bernstein et al. 2004; Luu & Jewitt 2002), while 90 km (Fuentes et al. 2009; Fraser & Kavelaars 2009;
a Fraser & Kavelaars (2009) have provided measurements Bernstein et al. 2004), giving way to a shallower dis-
[ down to diameters D & 15 km. The size distribution of tribution for smaller objects. Various models of the

large KBOs approximately follows a power law of the formation of the Kuiper Belt have been developed
1
(Benavidez & Campo Bagatin 2009; Kenyon & Bromley
v
Electronicaddress: [email protected] 2009; Pan & Sari 2005; Kenyon & Bromley 2004,
2 1Institute of Astronomy and Astrophysics, Academia Sinica.
8 11F of Astronomy-Mathematics Building, National Taiwan Uni- 2001; Benz & Asphaug 1999; Kenyon & Luu 1999a,b;
1 versity. No.1,Sec. 4,RooseveltRd,Taipei10617, Taiwan Davis & Farinella 1997; Stern 1996; Duncan et al. 1995)
6 2DepartmentofPhysicsandAstronomy,UniversityofPennsyl- which predict this break in the size distribution, but
. vania,209South33rdStreet,Philadelphia,PA19104 these models make vastly different predictions on the
1 3Harvard-Smithsonian Center for Astrophysics, 60 Garden
size distribution for objects smaller than the break
0 Street,Cambridge,MA02138
3 4Instituto de Astronom´ıa, Universidad Nacional Aut´onoma de diameter. These models are based on different scenarios
M´exico, Apdo. Postal 106, Ensenada, Baja California, 22800 of the dynamical evolution of the Solar System and the
1
M´exico internal structure of the KBOs themselves. It is gen-
: 5Steward Observatory, 933 North Cherry Avenue, Room N204
v erally assumed that the outward migration of Neptune
TucsonAZ85721
Xi 6New York University, Center for Cosmology and Particle stirreduptheorbitsofobjectsintheKuiperBelt,which
Physics,4WashingtonPlace,NewYork,NY10003. subsequently underwent a process of collisional erosion,
r 7DepartmentofAstronomyandUniversityObservatory,Yonsei giving rise to a shallower size distribution. The simplest
a University,134Shinchon, Seoul120-749,Korea
8Institute of Astronomy, National Central University, No. 300, of these models also predict a small-end power law size
JhongdaRd,JhongliCity,TaoyuanCounty320,Taiwan distribution of the form
9DepartmentofStatistics,UniversityofCaliforniaBerkeley,367
dn
EvansHall,Berkeley,CA94720 D−q, (2)
10Department of Astronomy, Universityof CaliforniaBerkeley, dD ∝
601CampbellHall,BerkeleyCA94720
11Max Planck Institute for Astronomy, K¨onigstuhl 17, D- where the slope q depends primarily on the internal
69117Heidelberg,Germany strength of the KBOs. Solid objects held together by
12Space Science andAstrobiology Division245-3, NASA Ames material strength give rise to larger slopes, while the
ResearchCenter,MoffettField,CA,94035
13KavliInstituteforParticleAstrophysicsandCosmology,2575 weaker gravitationally bound rubble piles are more eas-
SandHillRoad,MS29,MenloPark,CA94025 ily broken up and give rise to smaller values of q. A
14Initiative in Innovative Computing, Harvard University, 60 measurement of the size distribution of smaller objects
Ox1f5orDdepSat,rtCmaemntbriodfgeP,hMysAics0,21Y3a8le University, New Haven, CT (down to D & 100 m) is needed in order to constrain
these models. Furthermore, such a measurement would
06511
16YaleCenterforAstronomyandAstrophysics,YaleUniversity, provide important information on the origin of short-
P.O.Box208121,NewHaven,CT06520 period comets (Volk & Malhotra 2008; Tancredi et al.
2006; Duncan & Levison 1997; Levison & Duncan 1997;
Morbidelli 1997; Holman & Wisdom 1993).
2 Zhang et al.
The direct detection of such objects is difficult be- electronsarethencollectedinadditiontothosefromthe
cause they are extremely faint, with typical mag- firstrowblock,andthenthey areshifted to the next. By
nitudes R > 28, and are thus invisible to sur- the time the original photoelectons are read out, photo-
veys using even the largest telescopes. However, electronshavebeencollectedfromeverystaronthefocal
a small KBO will induce a detectable drop in the plane in a single rowblock, although at different epochs.
brightness of a distant star when it passes across While zipper mode readout allows us to image every
the line of sight to the star (Chang et al. 2012; star in the focal plane, there are two characteristics of
Schlichting et al. 2012; Wang et al. 2010; Bianco et al. thedatathatreducethesignal-to-noiseratio(SNR).The
2010; Schlichting et al. 2009; Bickerton et al. 2009; first is that sky background is collected as a rowblock is
Wang et al. 2009; Bianco et al. 2009; Liu et al. 2008; shifted acrossthe focal plane, while photoelectrons from
Zhang et al. 2008; Bickerton et al. 2008; Chang et al. a single star are only collected at one location. Further-
2007; Nihei et al. 2007; Chang et al. 2006; Roques et al. more, it takes about 1.3 ms to read out a single row,
2006, 2003; Cooray 2003; Cooray & Farmer 2003; correspondingto about95ms foranentirerowblock. So
Roques & Moncuquet 2000; Brown & Webster 1997; the actual exposure of a star is only 105 ms, while sky
Roques et al. 1987; Bailey 1976). High speed photom- background is collected for a total of 5.4 s. Second, flux
etry is needed to detect these events, given that we ex- is also collected from the stars during the readoutstage,
pectatypicaleventdurationofabout200ms. Detection and the brighter stars will leave significant numbers of
of such events is further complicated by the fact that photoelectrons in the pixels as they are being shifted.
the sizes of the objects of interest are on the order of This flux will appear as bright vertical streaks in the
the Fresnel scale (about 1.4 km for our median observed images (see the upper left panel of Figure1). Neverthe-
wavelengthofabout600nmandatypicalobjectdistance less, zipper mode is successful in collecting high SNR
of 43 AU), and the resulting lightcurves exhibit signifi- photometric data onenoughstarsfor the projectto suc-
cant diffration features (Nihei et al. 2007). The goal of cessfully probe a significant fraction of the allowed size
the TAOS project is to detect such occultation events distribution parameter space.
andmeasurethesizedistributionofKuiperBeltObjects In this paper we report on the analysis of seven years
(KBOs) with diameters 0.5 km<D <30 km. of zipper mode data collected by the TAOS system. In
To date, 18 events consistent with occultations by 2, we describe the data set used in this analysis. In 3,
§ §
objects in the outer Solar System have been detected. we provide a detailed discussion of the analysis pipeline.
The three detections claimed by Roques et al. (2006), if In 4, we present the results of the analysis of the data
§
they are occultations, are unlikely to have been caused set, and finally in 5, we discuss our future plans for the
§
by KBOs (the distances are inconsistent). On the survey.
otherhand,twodetectionsreportedby Schlichting et al. Throughouttheremainderofthispaper,weusethefol-
(2009) and Schlichting et al. (2012) are consistent with lowingdefinitions. Adata run isaseriesofhighcadence
occultation events by 1 km diameter KBOs at 45 AU. multi-telescope measurements on a single field, typically
Schlichting et al. (2012) found one additional event con- for 1.5 hours at a time. A lightcurve is a time series
sistentwithaKBOoccultation,buttheyrejectthisevent of photometric measurements of a single star from one
because of its high inclination (it was detected at an telescope in one data run, and a lightcurve set is a set
ecliptic latitude of b = 81.5◦). The remaining 12 pos- of multi-telescope photometric measurements of a sin-
sible events (Chang et al. 2007) were found in archival gle star during one data run. A flux tuple is a set of
RXTEmeasurementsofScoX-1. However,theinterpre- multi-telescopemeasurements(with a minimumofthree
tation of these detections as occultation events has been telescopes)ofonestaratasingleepoch. Alightcurveset
challenged (Jones et al. 2008; Blocker et al. 2009). is thus a time series of flux tuples.
TAOShasbeeninoperationsinceFebruary2005. The
system is described in detail in Lehner et al. (2009). In 2. THESEVENYEARDATASET
summary, the survey operates four 50 cm telescopes at Earlier results from the TAOS search for small KBOs
LulinObservatoryin centralTaiwan. The telescopes are were reported in Zhang et al. (2008) and Bianco et al.
very fast (F/1.9), with a field of view of 3 deg2. Each (2010). These two papers describe results obtained
telescope is equipped with a camerawhich utilizes a sin- from operating three telescopes. The dataset used in
gle 2k 2k CCD imager. All four telescopes monitor the Bianco et al. (2010) ended on 2008 August 2, when the
×
same stars simultaneously at a readout cadence of 5 Hz. fourthtelescopebecameoperational. Theadditionofthe
We require that any events be detected simultaneously fourth telescope to the arraywas beneficial in two ways.
in each of the telescopes in order to minimize the false First, using four telescopes increases the significance of
positive rate (Lehner et al. 2010). anycandidateoccultationevents,whilethefalsepositive
High speed cadence is achieved using a special CCD rate remains constant (Lehner et al. 2010). Second, we
readout mode called zipper mode, which is described in have experienced frequent problems with the reliability
detail in Lehner et al. (2009). To summarize, instead of of our cameras and telescopes. With four telescopes, we
taking repeatedexposures,we leavethe shutter open for can still operate with three telescopes when one of the
the duration of a data run. At a cadence of 5 Hz, we systems is shut down for repair. Given that a minimum
read out a subimage comprising 76 rows, which we have of three telescopes is necessary to keep the false positive
termed a rowblock. When the rowblock is read out, the rate low enough (Lehner et al. 2010), we can still collect
remaining photoelectrons on the CCD are shifted down usefuldatainthe eventthatthere isa problemwithone
by 76 rows as well. If one considers starting with a row- of the systems. The overall operational efficiency of the
block at the far edge of the CCD, photoelectrons are survey thus improved significantly since we could still
collected and then shifted to the next rowblock. Photo- collect data when one telescope was offline.
TAOS Seven Year Results 3
Table 1
Datasetparameters.
Zhangetal.(2008) Biancoetal.(2010) ThisWork
StartDate 2005Feb7 2005Feb7 2005Feb7
EndDate 2006Dec31 2008Aug2 2011Sep8
DataRuns 156 414 858
Light-curveSetsa 110,554 366,083 835,732(209,130)
TotalExposure(star–hours)a 152,787 500,339 1,159,651(292,514)
PhotometricMeasurementsa 7.8×109 2.7×1010 6.7×1010 (2.1×1010)
FluxTuplesa 2.6×109 9.0×109 2.1×1010 (5.3×109)
a Valuesforfour-telescopedatashowninparentheses.
Figure 1. Top left: Horizontal subsections of two consecutive rowblocks with a bright moving object in the center. (The boundary
between the rowblocks is visible as a dark horizontal line bisecting the image. This is a slight excess in dark current from the edge of
the CCD.) The vertical streaks (see text) from the brighter stars are also evident. Top right: The same rowblock subsections after the
background streak subtraction. Whilethe streaks from the brightstars areremoved, the movingobject itselfis not completely removed,
and the background flux near the moving object is over-subtracted. Bottom left: grayscale image representation of the background flux
subtracted from the rowblocks near those shown inthe top panels. Here each row represents the flux subtracted from a singlerowblock,
while the columns correspond to those in the top panels. The vertical streaks in this image are the streaks left by the stars during the
zippermodereadout. Themovingobjectinoriginaldatacausethebackgroundtobeover-subtracted,andtheover-subtractedfluxclearly
shows the trajectory of the moving object across the focal plane. Bottom right: The samebackground image from the bottom left panel
after applying the high pass mean and variance filters to remove the streaks from the bright stars. The over-subtracted background flux
fromthemovingobjectisclearlyevident.
The current dataset is summarized in Table1, along surements affected by bright objects (such as satellites
with the datasets used in Zhang et al. (2008) and or airplanes) moving across the field of view. As dis-
Bianco et al. (2010). The current dataset is more than cussed in Bianco et al. (2010), several candidate events
2.3 times larger than that used in Bianco et al. (2010). were found in the previous dataset. Visual inspection
90% of the photometric data in this set was acquired of the relevant images revealed that these events were
within 6◦ of the ecliptic in order to maximize the event caused by brigh moving objects passing near the star
rate. The results quoted in this paper are thus applica- that appeared to be occulted. The new algorithmto de-
ble toboththe hot(inclinationi>5◦)andcold(i<5◦) tect these bright moving objects is described in 3.1.1.
§
KBO populations.
3. DATAANALYSIS 3.1.1. Bright Moving Object Detection
The data analysis took place in three stages: photo- Asmentionedin 3.1,brightobjectsmovingacrossthe
§
metric reduction of raw image data, the search for oc- fieldofviewduringzippermodeimagecollectioncangive
cultation events,and the detection efficiency simulation. rise to a large number of false positive events. This is
Thesestagesroughlyfollowthesameproceduresoutlined causedbythebackgroundsubtractioninthephotometric
in Bianco et al. (2010), but several improvements were reduction. As discussed in 1, the zipper mode readout
§
made to the pipeline in this round of data analysis. The hasthedisadvantagethatthebrighterstarsleavestreaks
new pipeline is described in the following subsections. alongthe columns in the images due to the finite time it
takestoreadoutarowfromtheCCD.Thesestreaksare
3.1. Photometric Reduction
removed during the photometric reduction by subtract-
Forthefirststepinthepipeline,acustomphotometric inganestimateofthemodeofthecolumnfromeachpixel
reductionpipeline(Zhang et al.2009)isusedtomeasure inthecolumn,wherethemodeiscalculatedforeachzip-
thebrightnessofeachstarinthefieldateachepoch,and permoderowblock(Zhang et al.2009). Abrightmoving
theresultingmeasurementsareassembledintoatimese- objectpassingovera particularcolumncausesthe mode
rieslightcurveforeachstar. Foreverystar,thelightcurve ofthatcolumntobeover-estimated,andthebackground
from each telescope is combined into a lightcurve set. is thus over-subtracted, causing an artificial drop in the
Next, bad data points (e.g. measurements where the measuredbrightnessofanystarwhoseaperture includes
star lies near the edge of the focal plane) are removed that column. (see the top panels of Figure1).
from the lightcurves (see Zhang et al. 2009, for details). Inorderto detectthese events,wesavethe 3σ-clipped
A significant improvement added to this stage in the column mean for each rowblock and column on all four
process is the automated flagging of photometric mea- cameras. ThesebackgrounddataarestoredinFITSfor-
4 Zhang et al.
Figure 2. Histogramsofp-valuesforthePearson’sχ2 test(leftpanel)andhypergeometrictest(rightpanel)foreachofthedatarunsin
thedataset.
mat, and an example is shown in the bottom left panel than 2.5. We then flag any photometric measurements
of Figure1. Note that each column in this image cor- where a column affected by a moving object lies within
responds to an actual column on the CCD, but the row the aperture used in the photometry. The removal of all
correspondsto a rowblockin the zipper mode data. The flagged data points results in the elimination of all false
pixelvalueissimplythecolumnmean(after3σclipping) positive events caused by moving objects.
for that rowblock and column. The streaks from the
brightstarsareevidentintheimage. Alsoclearlyvisible 3.2. Event Search
is the object shown in the top panels moving across the 3.2.1. Rank Product Statistic
image. We note that what appears as a moving object
is actually the high value for the column average due to Our event selection algorithm uses the rank product
the moving object. Nevertheless, it is straightforwardto statistic, which is described in detail in Lehner et al.
identify which columns in which rowblocks are affected (2010). This statistic takes advantage of the multi-
by the over-estimated background. telescope data to give an exact value for the significance
First, we normalize the image by applying the same of any candidate event, even if the underlying distri-
rolling clipped mean and variance filter to each of the bution of the flux measurements is not known. Given
columnsweusetoremovetheslowlyvaryingtrendsalong that TAOS samples at a 5 Hz readout cadence, we can-
the lightcurves (Lehner et al. 2010). That is, we replace not resolve the features of any candidate events in our
eachcolumnaveragevalueb ,wheremistherowblock lightcurves. Any occultation event we detect would typ-
mn
indexandnisthecolumn,withanormalizedvalued , ically appear as a small drop in the flux for one or two
mn
defined as consecutive measurements (although occultation events
by largerobjects couldcontainup to five or six consecu-
cmn=bmn µmn, tively measured flux drops). We therefore need to have
−
cmn anaccurateestimateofourfalsepositiverateinorderto
d = , (3)
mn σ report a credible result.
mn
To searchfor events in the data, we calculate the rank
whereµ is the local(3σ-clipped)meanwithincolumn
mn productstatistic oneachtuple inthe entire dataset. To
nandσ thelocalvariancewithincolumnn,calculated
mn calculate this statistic, we take a time series of fluxes
with window sizes of 21 rowblocks. Note that σ is
mn f ,f ,...,f (where N is the numberofpointsonthe
calculatedafterthemeanhasbeensubtractedfromevery 1 2 NP P
time series), and replace each measurement f with its
j
rowblock backgroundvalue along the column.
rank r . The lowest measurement in the time series will
j
ThebottomrightpanelofFigure1showstheresulting
be givena rankof1,andthe brightestmeasurementwill
values of d after application of this filter to the im-
mn begivenarankofN . Foragivenlightcurveset,wethus
P
age shown in the bottom left panel. We then apply an
replacethetimeseriesoffluxtuples withatimeseriesof
algorithm(Lutz 1980)1 to find objects in the normalized
rank tuples of the form
background image. We define an object as either a set
of 20 or more adjoining pixels with values of d larger (r ,...,r ),
mn 1j Tj
than1.5,orasasetof5ormorepixelswithvalueslarger
where T is the number of telescopes used. We do not
1 This is the same algorithm used in SExtractor know the underlying distribution of the flux tuples in
(Bertin&Arnouts1996)todetect objects. each lightcurve set, but, assuming the lightcurves meet
TAOS Seven Year Results 5
Figure 3. Lightcurveswithsimulatedoccultationeventsfrom1km(top)and3km(bottom)diameterobjectsaddedin,beforeandafter
application of the movingaverage. In each panel, the top lightcurve is the simulated event, the second lightcurve has no moving average
applied,andthenextthreehavemovingaverages withwindowsizes2,3,and5respectively.
the correct requirements (discussed in 3.2.2), we do rankproductwouldbesignificantlysmallerthanaverage,
§
know the exact distribution of rank tuples, and we can leading to a large value of our test statistic z. We set a
thuscalculatethesignificanceofanycandidateeventex- threshold z for event detection based on choosing an
t
actly. upperlimitontheexpectedfalsepositiverate,suchthat
The rank product statistic is calculated as
0.25
P(Z >z )= , (5)
T rij t NT
z = ln . (4)
j
− N whereZ isarandomvariablechosenfromtherankprod-
Y P
i=1
uct probability distribution, 0.25 is our limit on the ex-
If one assumes that the rank distribution for a single pected number of false positive events in the data set,
lightcurve is continuous, then this statistic follows the and N is the total number of tuples used in the event
T
gamma distribution. Given that the ranks are discrete, search.
theapproximationfailsforlargevaluesofz,however,the
probability distribution can be calculated exactly using 3.2.2. Lightcurve Filtering and Data Quality Checks
simple combinatorics.
The rank product test statistic follows the expected
Inthecaseofanoccultationevent,oneexpectstheflux
distributionforanindividuallightcurvesetifandonlyif
to drop below the nominal value simultaneously in each
the following conditions are satisfied (Coehlo 2010).
of the telescopes. The corresponding rank tuple would
have lower values for all of the ranks, and the resulting The distribution of measurements in each
•
6 Zhang et al.
Figure 4. SameasFigure3,butfor8km(top) and30km(bottom)diameterobjects.
lightcurve in the lightcurve set is stationary. where µ is the local mean and σ the local variance at
j j
point j, calculated with window sizes of 33 and 151 re-
Each lightcurve in the lightcurve set is ergodic in
spectively. Note that the variance is calculatedafter the
•
mean.
rolling mean was subtracted, as was done in Equation3.
In most cases, this filter removes all of the slowly vary-
Eachlightcurveinthelightcurvesetisindependent
• ing trends, and the resulting filtered lightcurves h meet
of the other lightcurves in the lightcurve set.
the requirements for the application of the rank product
However,noneofthelightcurvesetsinourdatasetmeet statistic.
these requirements. Changes in the transparency of the However, there are still some data runs where fast
atmospherethroughoutthecourseofourtypical1.5hour moving cirrus clouds can induce simultaneous variations
data runs induce fluctuations in the lightcurves, render- in the lightcurves that are not removed by the high
ing them non-stationary. Furthermore, these fluctua- pass filter. We have therefore devised a pair of tests
tions are correlated between the different telescopes, so (Lehner et al. 2010) to identify lightcurve sets where
the lightcurves are not independent of each other. We the filtered lightcurves are not independent (and likely
thereforeapplyarollingmeanandvariancefiltertoeach not stationary as well, given that the fluctuations in
point in the lightcurves, thus replacing each flux mea- the lightcurves from the cirrus clouds are not constant
surement fj with hj, given by over time). The tests are based on the fact that the
rank tuples should be distributed uniformly over the T-
g =f µ ,
j j j
− dimensional rank space (recall that T is the number of
g
h = j. (6) telescopes used). For the first test, we divide the T-
j
σ
j
TAOS Seven Year Results 7
Figure 5. Toppanel: histogramofSNRvaluesforlightcurvesets Figure 6. Plots of Ωeff vs D for the optimum combinations of
intheentiredataset. (Foreachlightcurveset,theminimumSNR movingaveragewindowsizesandSNRcuts. Thevaluesarescaled
amongalltelescopesisusedforthehistogram.) Lowerpanels: SNR byΩB10,whichisΩeff fromBiancoetal.(2010),inordertohigh-
valuesofrecoveredevents fromourefficiencysimulationfor1km, lightthedifferences.
8km,and30kmobjects(solidhistogram). Alsoshown,usingthe
rightverticalaxis,isthedetectionefficiencyvsSNR(dottedlines).
our effective sky coverage, so we lose little by excluding
them from the final results. Furthermore, as will be dis-
dimensional rank space uniformly into a T-dimensional
cussedin 3.2.4, mostofthe starsinthesedataruns will
grid. If the lightcurves are independent, we expect that §
end up being excluded anyway.
each cell in the grid would contain the same number of
ranktuples. We thus calculatethe Pearson’sχ2 statistic
3.2.3. Detection of Multipoint Occultation Events
on the number of tuples in each cell of the grid. For the
second test, we only look at the number of rank tuples The rank product statistic described in 3.2.1 only
§
in the lowest ranked cell in each dimension. When look- works on one rank tuple at a time. This has the disad-
ing at all of the lightcurve sets in a given data run, the vantage of being inefficient at detection of longer dura-
results from the first test should follow the χ2 distribu- tion occultation events (such as objects with D &5 km,
tion, and the results from the second test should follow scattered disk or Oort Cloud objects beyond 100 AU,
the hypergeometric distribution. (However, this is only or events measured at large opposition angles). How-
true if the lightcurves are independent and identically ever, as was shown in Coehlo (2010) and discussed in
distributed,whichisastricterrequirementthanstation- Lehner et al.(2010),ifwereplaceastationarylightcurve
ary. We have thus devised a blockwise bootstrap method with a moving average of the form
to calculate the expected distributions. We still call the
1
teststhe Pearson’sχ2 testandthe hypergeometrictest.) aj = (hj +...+hj+w−1), (7)
w
We can thus calculate a p-value for both tests for each
dataruntoseehowwellthetwoteststatisticsmatchthe where h is the lightcurve after high pass filtering and
j
expected distributions. Histograms of these p-values are w is the rolling average window size, then the resulting
showninFigure2foreachtest. Ifthestatisticsmatchthe lightcurve will also be stationary, and the rank product
distribution, we would expect a uniform distribution of statistic will still be applicable. By applying a moving
p-values in the range of [0,1]. The large spikes evident average,wecanincreasethesignificanceofanycandidate
in the lowest bins of each histogram are the data runs event since the average lightcurve will exhibit a more
that have correlated fluctuations in the lightcurves that significantdropinthe measuredfluxwhentheaveraging
are not removed by the high pass filter. This leads to window is centered on the event.
a non-uniform distribution of rank tuples, meaning that This is illustratedin Figures 3 and 4, which shows the
the teststatistics for eachofthe lightcurvesets does not results of the application of a moving average on actual
follow the expected distribution. We thus remove every lightcurves with simulated events added in. We have
data set with a p-value less than 0.01. Note that this used three different window sizes w = 2, 3, and 5, in
means we would only expect to remove 1% of our data ordertooptimizeourdetectionefficiency. Wealsolooked
runs due to random chance. atthelightcurveswithnoaverageapplied,whichwewill
Inordertoperformthesetests,werequireenoughhigh refertoasusingawindowsizew =1inordertosimplify
(SNR stars to accurately calculate these two test statis- the discussion. For the 1 km diameter simulated event,
tics. We thus set a minimum of 10 stars with SNR>10 we detect the event using w = 1, 2, and 3. For such
in the data run. Any data run that does not meet this a small object, using w = 5 actually smoothes out the
requirementisexcludedfromfurtheranalysis. Visualin- signal,while amplifying other smallfluctuations, so that
spection of lightcurves in such data runs indicate that the event does not pass the cuts. For the 3 km and
there are very few stars, and the SNR is compromised, 8kmdiameterobjects,theeventsaredetectedonlywith
usually by very bright sky during a full moon. These w = 1 and 2. For the 30 km object, which is longer in
data runs are not expected to contribute significantly to duration than the other simulated events, the signal is
8 Zhang et al.
Thefinalcutweapplytoourdatasetistosetathresh-
old on SNR of each lightcurve set. The goal of this cut
is to exclude stars which provide little sensitivity to oc-
cultation events, yet still contribute to the total number
of tuples N , which tightens the selection threshold (see
T
Equation8). This is illustrated in Figure5, which shows
the total distribution of SNR values for every lightcurve
set in the data set, as well as the SNR distributions of
stars for which simulated events are recovered for 1 km,
8 km, and 30 km objects. It is clear that there is little
contribution to our effective sky coverage for SNR . 5.
Also shown is the detection efficiency for the same di-
ameter objects vs. SNR. As expected, the detection ef-
ficiency drops off at low SNR values, especially for the
smaller objects.
The exact value of SNR to use as a threshold depends
on many factors, especially the set of moving average
window sizes that will be used. This optimizationis dis-
cussed in 4.
§
Figure 7. Ωeff vs D using our final selection criteria of w = 1
and 3 and SNR > 3. Also shown for comparison are the pre- 4. RESULTS
vious TAOS results from Biancoetal. (2010) (dotted line) and
We searchedthrough our data set, using moving aver-
Zhangetal.(2008)(dashedline).
agewindowsizesof1,2,3,and5,(andeverycombination
not detectable for w = 1, but it is detected with w = 2, thereof) as well as SNR cuts of 0, 1, 2, 3, 4, 5, 6, and 7.
3, and 5. For each combination of window size and SNR thresh-
It is clear that the different window sizes increase the old, the detection threshold z was calculated as shown
t
sensitivity to different size objects. However, while the in Equation8. In all cases, no events were found in the
effective sky coverage of the survey can be increased by dataset. Sothenextstepintheanalysiswastomeasure
usingawidearrayofwindowsizes,thiswillalsoincrease our effective sky coverage and optimize the selection of
thefalsepositiverate(Lehner et al.2010). Runningeach moving average window sizes and SNR cut.
multiple window size on each lightcurve set is not an in- The total number of events expected by the survey if
dependent test, given that events can be detected using given by
morethanonewindowsize,asdemonstratedinFigures3 dn
and4. Itfollowsthatanyfalsepositiveeventcouldbede- Nexp =Z dD Ωeff(D)dD, (9)
tectedmultiple times using differentwindow sizes. How-
ever,giventhatwedonotknowthe underlyingdistribu- where dn/dD is the differential size distribution (the
tion of flux values in the data set, we can not estimate number of objects per deg−2 per km along the ecliptic),
the rate at which this may occur. Therefore, instead of
and Ω is the effective sky coverage of the survey. We
estimatingthetruefalsepositiverate,wecanonlysetan eff
estimate Ω by implanting a simulated event into each
upper limit on the false positive rate by assuming that eff
lightcurvesetandseeingifitisrecoveredbyourselection
searchingforeventsinalightcurvesetafterapplyingthe
criteria. Forthissimulation,wechooseasetofdiameters
movingaverageisacompletelyindependenttestforeach
between 0.5 km and 30 km, and for each diameter, we
value of w. If we have a total of N window sizes, we
w implant exactly one event into each lightcurve set. We
thus modify Equation5 as
assumeeveryobjectisatageocentricdistanceof43AU,
0.25 since the lightcurve shape does not depend critically on
P(Z >zt)= . (8) thedistanceformostobjects(40to46AU)intheKuiper
N N
T w Belt. Thetestsforeachdiameterareperformedindepen-
So the more window sizes that are used, the tighter the dently, so there is never more than one event implanted
detection threshold must be. The optimization of the into any lightcurve set at any time. The effective sky
selectionofthe windowsizesusedinthisanalysiswillbe coverageis calculated as
discussed in 4.
Note that§the data quality tests described in 3.2.2 H vrel
must be reapplied for each window size w. As sho§wn in Ωeff(D)= ∆ ∆ El, (10)
Lehner et al. (2010), application of the moving average lX∈rec
can highlight small correlations between the telescopes where the sum is over only those lightcurve sets l where
that are insignificant in the unaveraged data. There- the added event was recovered, E is the length of the
l
foreusinglargerwindowsizeswillrequireremovingmore lightcurvesetintime,v istherelativetransverseveloc-
rel
data runs fromthe analysis. Incases where we use more ity between the Earthand the KBO,∆ is the geocentric
than one window size, we always choose to be conser- distance to the KBO (again, fixed to a constant value of
vative and remove the union of data sets that would be 43 AU), and H is the event cross section. We estimate
rejected using each value of w. H as
3.2.4. Signal to Noise Ratio Cuts H(D,∆,θ∗) (2√3F32 +D23 23 +θ∗∆, (11)
≈h i
TAOS Seven Year Results 9
Table 2
ResultsfromusingdifferentwindowsizesandSNRcuts.
WindowSizes SNRCut Tuples zt q Ωeff(D=1km) Ωeff(D=30km)
1&2 3 9.49×109 1.32×10−11 3.82 (3.590±0.053)×10−7 (2.788±0.011)×10−5
1&2 4 7.13×109 1.75×10−11 3.81 (3.703±0.054)×10−7 (2.577±0.010)×10−5
1&2 5 5.63×109 2.22×10−11 3.82 (3.806±0.055)×10−7 (2.277±0.010)×10−5
1&3 3 9.14×109 1.37×10−11 3.82 (3.467±0.053)×10−7 (5.065±0.014)×10−5
1&3 4 6.87×109 1.82×10−11 3.82 (3.559±0.053)×10−7 (4.430±0.013)×10−5
1&3 5 5.42×109 2.31×10−11 3.82 (3.657±0.054)×10−7 (3.770±0.012)×10−5
1&5 3 8.60×109 1.45×10−11 3.84 (3.213±0.051)×10−7 (5.594±0.015)×10−5
1&5 4 6.46×109 1.94×10−11 3.84 (3.301±0.052)×10−7 (4.824±0.014)×10−5
1&5 5 5.08×109 2.46×10−11 3.84 (3.413±0.052)×10−7 (4.071±0.013)×10−5
Figure 8. Solid line: 95% c.l. upper limit on cumulative surface density vs. diameter D (bottom axis) and R magnitude (top axis,
assuminganalbedoof4%andadistanceof43AU)fromcurrentTAOSdataset. Dottedlines: 95%c.l. upperandlowerlimitsonsurface
density from Schlichtingetal. (2012). Cross and error bar (S12): Surface density reported by Schlichtingetal. (2012). (Note that this
is not strictly model independent because it is based on a power law model. However, the result does not depend very strongly on the
assumedslope.) Solidsquares(RXTE):upperlimit(rightpoint)reportedbyLiuetal.(2008)andthebestfitsurfacedensity(leftpoint)
reported by Changetal. (2007) from an occultation search through RXTE observation of Sco X-1. Solid triangles (PS): upper limits
reportedbyWangetal.(2010)usingPan-STARRSguiderimages. Emptytriangle(BKW):upperlimitreportedbyBickertonetal.(2008)
using the 1.8 m telescope at DAO. Empty circles (MMT): upper limits reported by Biancoetal. (2009) from analysis of trailed images
obtainedwiththeMMT.Upsidedownemptytriangle(R06): upperlimitreportedbyRoques etal.(2006)fromhighspeedimagingusing
the 4.2m Herschel Telescope. Empty square(HST): upper limitreported byBernsteinetal.(2004)froma directsurvey usingthe HST.
Solidcircle(FGH):surfacedensityreportedbyFuentes etal.(2009)fromadirectsearchusingSubaru. Star(FK):surfacedensityreported
by Fraser&Kavelaars (2009) from a direct survey using Subaru. (Note that Fraser&Kavelaars (2009) assume an albedo of 6% and a
distance of 35 AU, so they claim that the point plotted at R = 27.875 corresponds to a diameter D = 15 km.) Lines and box labeled
CB, P, and SD: requiredsurface density for the Classical Belt (Levison&Duncan 1997), Plutinos (Morbidelli 1997), and Scattered Disk
(Volk&Malhotra2008)respectivelytobethesourcefortheobserveddistributionofJupiterfamilycomets.
10 Zhang et al.
where θ∗ is the angular size of the target star and F is
the Fresnel scale, which is included to account for the
minimum cross section of an event due to diffraction
(Nihei et al. 2007).
To optimize the selection of window sizes and SNR
cut, we looked at the resulting values of Ω for every
eff
possible combination. There is no clear optimum set of
cuts, so we decided to optimize our selection based on
three criteria.
Minimizetheupperlimitoftheslopeqatthesmall
•
end of the size distribution, assuming a power law
size distribution anchored at the break diameter
of 90 km and a cumulative surface density at the
break diameter of 5.4 deg−2 (Fraser & Kavelaars Figure 9. Numberofexpectedeventsvsslopeq. Solidline: num-
2009; Schlichting et al. 2009). ber of expected events using a surface density of ND>90km =
5.4 deg−2, as reported by Fraser&Kavelaars (2009). Dashed
Maximize the sensitivity at D = 1 km, where line: number of expected events vs. slope q, using a surface den-
• Schlichting et al. (2009) and Schlichting et al. sity of ND>90km = 38.0 deg−2, as reported by Fuentes etal.
(2009). Since no events were found, any model predicting more
(2012) claim the detection of two KBOs.
than3events (dotted line)isexcludedatthe95%c.l.
MaximizethesensitivityatD =30km,inorderto of slope q. The results are shown in Figure9. We use
•
bringourupperlimitascloseaspossibletothecur- two values for the surface density at D = 90 km, which
rent direct detection limits (Bernstein et al. 2004; normalizethesizedistributionforsmallerdiameters. Us-
Fuentes et al. 2009; Fraser & Kavelaars 2009). ing the value reported by Fraser & Kavelaars (2009) of
N = 5.4 deg−2 yields a 95% c.l. upper limit of
D>90km
The best results came fromusing cuts on SNR of 3, 4,
q =3.82. In comparison,Schlichting et al. (2009) report
and5. TheresultsaresummarizedinTable2. Wefound aslopeofq =3.8 0.2basedontheirclaimoftwodetec-
that the upper limit on the slope q did not vary signif- ±
tionsatD =1km, assumingthe inclinationdistribution
icantly as long as we limited ourselves to two window
of small KBOs follows that of the larger (D > 100 km)
sizes, and included w = 1. The combination of w = 1 objects (Elliot et al. 2005). The second value we use is
and2anda minimum ofSNR>3 gavethe largesteffec- N = 38.0 deg−2 as reported by Fuentes et al.
tiveskycoverageat1km, however,thiswasparticularly D>90km
(2009). This results in a 95% c.l. upper limit on q of
badat30km. Thebestresultat30kmcamefromw =1
3.34.
and 5 with SNR>3, but this combination was not very
good for 1 km. The combinations that worked best for 5. FUTUREPLANS
both diameters was using w =1 and 3 with SNR cuts of
After seven years of observations, it would only be of
3 and 4. The resulting values of Ω for these four com-
eff marginalvaluetocontinuethesurveyinitspresentform.
binationsareshowninFigure6. Intheend,weoptedfor
So at the end of the data set discussed in this paper, we
using w = 1 and 3, with a threshold of SNR >3 as this
shut down the system for a camera upgrade. Each new
is better at 30 km, and not much worse at 1 km.
camera, manufactured by Spectral Instruments, utilizes
Givenourfinalsetofcuts, the resultingplot ofΩ vs
eff a CCD47-20 frame transfer CCD from e2v. With the
DisshownFigure7. Thepointsonthecurveindicatethe
new cameras, we can image with a readout cadence of
diameters where we calculated Ω using our detection
eff just under 10 Hz if we use 2 2 binning. The CCDs
efficiencysimulation,andthecurveitselfisacubicspline ×
are 1k 1k with 13 µm pixels, as opposed to the 2k 2k
fit to these values. × ×
imagers with 13.5 µm pixels used to collect the current
To set a model independent upper limit, with no de-
data set, so 77% of the field of view is lost. However,
tected events we can eliminate any model which would
bymovingawayfromzippermodetofullframeimaging,
lead us to expect more that three detected events at the
ourSNRincreasessignificantlyandweestimatethatour
95%c.l. Our model independent upper limit is shownin
limiting magnitude will increase from R = 13.5 to 15.
Figure8, along with results from other occultation sur-
We thus expect to be able to monitor a similar number
veys and direct searches. We note that our upper limit
of stars. Given the higher readout cadence, we become
at D = 1 km is a factor of 4.5 below that reported by
significantly more sensitive to objects with D . 1 km.
Schlichting et al. (2012), but it is also a factor of 6.7
ObservationswiththesenewcamerasbeganinNovember
higher than their lower limit based upon their two re-
2012.
portedevents. Thereforeeventhoughwefoundnoevents
at 1 km (despite the fact that our efficiency simulation
for 1 km diameter objects showed that we are in fact Work at ASIAA was supported in part by the the-
sensitive to events similar to those detected by the HST matic research program AS-88-TP-A02. Work at NCU
survey), our limit is thus still consistent with their pub- andatLulinObservatorywassupportedinpartbygrant
lished result. NSC101-2628-M-008-002. Work at the CfA was sup-
We also calculate upper limits on the small KBO size ported in part by the NSF under grant AST-0501681
distributionassumingitfollowsa powerlawanchoredat andbyNASAundergrantNNG04G113G.WorkatYon-
the detected break at D = 90 km. We use Equation9 seiwassupportedbytheNRFgrant2011-0030875. Work
to calculate the number ofexpected events as a function atNCUwassupportedbythegrantNSC96-2112-M-008-