Exoplanet detection in angular differential imaging by statistical learning of the nonstationary patch

29 Pages·2017·28.97 MB·English

A&A618,A138(2018) Astronomy
https://doi.org/10.1051/0004-6361/201832745
&
(cid:13)c ESO2018 Astrophysics
Exoplanet detection in angular differential imaging
by statistical learning of the nonstationary patch covariances
The PACO algorithm
Olivier Flasseur1,Loïc Denis1,Éric Thiébaut2,andMaud Langlois2
1 UniversitédeLyon,UJM-Saint-Etienne,CNRS,Institutd’OptiqueGraduateSchool,LaboratoireHubertCurienUMR5516,42023
Saint-Etienne,France
e-mail:[email protected], [email protected]
2 Université de Lyon, Université Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR 5574,
69230 Saint-Genis-Laval,France
e-mail:[email protected]
Received31January2018/Accepted22March2018
ABSTRACT
Context.Thedetectionofexoplanetsbydirectimagingisanactiveresearchtopicinastronomy.Evenwiththecouplingofanex-
tremeadaptive-opticssystemwithacoronagraph,itremainschallengingduetotheveryhighcontrastbetweenthehoststarandthe
exoplanets.
Aims.Thepurposeofthispaperistodescribeamethod,namedPACO,dedicatedtosourcedetectionfromangulardifferentialimaging
data.Giventhecomplexityofthefluctuationsofthebackgroundinthedatasets,involvingspatiallyvariantcorrelations,weaimto
showthepotentialofaprocessingmethodthatlearnsthestatisticalmodelofthebackgroundfromthedata.
Methods.Incontrasttoexistingapproaches,theproposedmethodaccountsforspatialcorrelationsinthedata.Thosecorrelations
andtheaveragestellarspecklesarelearnedlocallyandjointlytoestimatethefluxofthe(potential)exoplanets.Bypreventingfrom
subtractingimagesincludingthestellarspecklesresiduals,thephotometryisintrinsicallypreserved.Anonstationarymulti-variate
Gaussianmodelofthebackgroundislearned.Thedecisioninfavorofthepresenceortheabsenceofanexoplanetisperformedbya
binaryhypothesistest.
Results.ThestatisticalaccuracyofthemodelisassessedusingVLT/SPHERE-IRDISdatasets.Itisshowntocapturethenonstationar-
ityinthedatasothatauniquethresholdcanbeappliedtothedetectionmapstoobtainconsistentdetectionperformanceatallangular
separations.Thisstatisticalmodelmakesitpossibletodirectlyassessthefalsealarmrate,probabilityofdetection,photometricand
astrometricaccuracieswithoutresortingtoMonte-Carlomethods.
Conclusions.PACOoffersappealingcharacteristics:itisparameter-freeandphotometricallyunbiased.Thestatisticalperformancein
termsofdetectioncapability,photometricandastrometricaccuraciescanbestraightforwardlyassessed.Afastapproximateversion
ofthemethodisalsodescribedthatcanbeusedtoprocesslargeamountsofdatafromexoplanetssearchsurveys.
Keywords. techniques:imageprocessing–techniques:highangularresolution–methods:statistical–methods:dataanalysis
1. Introduction Bonavitaetal.2014;Macintoshetal.2015;Chauvinetal.2017).
Thisismainlyduetothedifficultyindetectingafaintexoplanet
Directimagingofexoplanetsisarecentobservationaltechnique signal given the very large contrast (greater than 10−5 in the
(Traub&Oppenheimer2010)particularlyadaptedtotheobser- near infrared) between the companions and their host star. In
vationofyoungandmassiveexoplanets.Itiscomplementaryto addition to the use of extreme adaptive optics, optical atten-
conventional indirect imaging methods (Santos 2008) such as uation of the host star by using a coronagraph is mandatory
Dopplerspectroscopyortheradialvelocitytechniquethathaveled to reach this level of contrast (Soummer 2004; Mawetetal.
tothedetectionandcharacterizationofseveralhundredsubJovian 2006). Currently, two exoplanet-searchers are optimized using
mass exoplanets (Schneideretal. 2011). Direct imaging is also these instrumental techniques for direct imaging observations:
a privileged observation method to reconstruct the spectrum of the Spectro-Polarimetry High-contrast Exoplanet REsearch
lowmasssubstellarobjects(Viganetal.2010).Thisinformation (SPHERE; Beuzitetal. 2008) at the Very Large Telescope
is crucial for estimating physical parameters critical for the (VLT) of the European Southern Observatory (ESO) and the
characterization of these objects. In particular, it is mandatory GeminiPlanetImager(GPI;Macintoshetal.2014)attheGEM-
to estimate the objects age, surface gravity, composition, or INI observatory. These two instruments offer spectroscopic
effective temperature (Allardetal. 2003, 2007) and to predict observations capability in the near infrared via integral fields
their evolution (Burrowsetal. 1997; Chabrieretal. 2000) by spectrographs(IFS)orlongslitspectroscopy(Viganetal.2010).
refiningthenumericalmodelsoflowmassobjectsandtheexo- The VLT/SPHERE is also equipped with a differential imager,
planetformationtheories.Despitethepromisingfutureofdirect the InfraRed Dual Imaging Spectrograph (IRDIS; Dohlenetal.
imaging, only a few dozen exoplanets have been successfully 2008a,b) optimized for imaging exoplanets using angular and
detected to this day using this technique (Lagrangeetal. 2009; spectraldifferentialimaging.Inordertoenhancetheachievable
A138,page1of27
OpenAccessarticle,publishedbyEDPSciences,underthetermsoftheCreativeCommonsAttributionLicense(http://creativecommons.org/licenses/by/4.0),
whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited.
A&A618,A138(2018)
contrast, several observation strategies have been developed Insteadofonlyminimizingthenoise(i.e.,thenormoftheresid-
to exploit temporal or spectral diversity. Spectral differential uals),TLOCIalsomaximizestheexoplanetsignal-to-noiseratio
imaging(SDI;Racineetal.1999)isbasedontheassumptionthat (S/N) in the residuals. In other words, the influence of a spe-
stellarspecklesfromdiffractionarestronglycorrelatedfromone cific choice of linear combination on the reduction of the flux
wavelengthtoanother(whenthedifferentialaberrationsaresmall) of the candidate companion is also considered. The TLOCI al-
aftercompensatingforthewavelengthscaling.Multipleimages gorithm is often calibrated for exoplanet signal self-subtraction
at different wavelengths can thus be combined to numerically byinjecting“fake”exoplanets(faintpointsources)intothedata
suppressalargefractionofthestellarspeckles.Inaddition,the to determine the algorithm throughput at each position in the
SDItechniqueassumesthattheexoplanetfluxisdetectableatone fieldofviewafterthespeckleremoval.Foreachfakeexoplanet
wavelengthwhileitistooweakattheotherwavelength.Angular injected, the ratio of the exoplanet flux in the resulting image
differentialimaging(ADI;Maroisetal.2006)consistsoftracking to its initial flux is estimated to produce the 1D-throughput as
theobservationtargetovertime(thetelescopederotatoristunedto a function of the angular separation. The MLOCI algorithm
maintainthetelescopepupilstablewhilethefieldofviewrotates). (Wahhajetal. 2015) injects fake point sources and maximizes
Consequently,candidatecompanionsdescribeanapparentmotion theirS/N,whichalsoimprovestheS/Nofthesourcespresentin
along a predictive circular trajectory around the host star while thedata.IntheALOCIalgorithm(Currieetal.2012a,b),thedata
the telescope pupil(including spiders) and coronagraph remain aredividedintoannuliandprocessedindependently,therebyal-
static.Specklesresultingfrominstrumentalaberrationsarethus lowing different linear combinations at different angular sep-
stronglycorrelatedfromoneexposure totheother.The images arations in the definition of the stellar PSF. Another method,
canbecombinedtocancelmostofthespeckleswhilepreserving called ANDROMEDA (Mugnieretal. 2009; Cantalloubeetal.
partofthesignalfromtheoff-axissources.Theirapparentmotion 2015), forms differences of temporal images to suppress stellar
helpstoseparatethemfromthespecklebackground.Detectability specklesandperformsthedetectionofdifferentialoff-axisPSFs
oftheexoplanetsthereforereliesonthecombinedabilityofthe (i.e., the signature of an exoplanet in the difference images). A
instrument and the numerical processing to suppress the light generalizedlikelihoodratiotestisthenevaluated.Sincethetech-
fromthehoststarandextractthesignalfromtheexoplanetsfrom niques discussed so-far are based on image differences that al-
the remaining stellar speckles. Elaborate processing methods ways reduce the amplitude of the signal due to the exoplanets,
combiningmulti-temporaland/ormulti-spectraldataplayacen- adaptedpost-processingcalibrationsmustbeappliedtorecover
tralroletoreachtheultimatedetectionlimitachievablebysuch the photometry of the detected exoplanets. MOODS algorithm
instruments. (Smithetal. 2009), based on the joint estimation of the exo-
Thispaperisorganizedasfollows.InSect.2wereviewthe planet amplitude and of the stellar PSF (considered constant
currentmethodsforexoplanetdetectioninADIorSDI.InSect.3 for all observations of a temporal dataset), circumvents this
wethenintroducetheprincipleofourmethodbasedonamod- problem.
eling of patch covariances. The detailed implementation of the Another family of methods considers that the fluctuations
algorithmisdescribedinSect.4.Weassesstheperformanceof of the stellar speckles (i.e., the stellar PSF) span a small-
ourmethodonVLT/SPHERE-IRDISdatasetsinSect.5.Finally, dimensionalsubspace.Exoplanetsarethusdetectedonthesub-
Sect.6drawstheconclusionsofthepaper. space orthogonal to the subspace that captures fluctuations of
thestellarspeckles.Thedataareprojectedinanorthogonalba-
siscreatedbyprincipalcomponentsanalysis(PCA).Thisisthe
2. Currentmethodsforexoplanetdirectdetection principle of the KLIP algorithm (Soummeretal. 2012) which
builds a basis of the subspace capturing the stellar PSF by per-
usingangularand/orspectraldifferentialimaging
formingaKarhunen-Loèvetransformoftheimagesfromtheref-
Several methods have been developed to combine different im- erence library. To obtain a model of the stellar PSF to subtract
agestakenatseveralconsecutivetimesand/orwavelengths.Stel- in order to attenuate the speckles, the science data is projected
larspecklescanbeattenuatedbysubtractingareferencestellar onto a predetermined number of modes. Even if the general
point spread function (PSF). This is the principle of the LOCI principle is close to the LOCI type algorithms, KLIP is much
algorithm (Lafreniereetal. 2007). The reference stellar PSF to faster thanks to the truncation. The sKLIP (Absiletal. 2013)
be subtracted to the data is created by a linear combination of algorithm partially prevents the subtraction of signal from the
images selected in a library of data acquired under experimen- companionsbybuildingthereferencelibraryonlyfromimages
talconditionssimilartothoseoftheobservationofinterest.The where the candidate exoplanets underwent a sufficient rotation.
optimizationofthiscombinationisperformedbyminimizing,in The LLSG (Gonzalezetal. 2016) algorithm is also related to
the least squares sense, the residual noise inside multiple sub- subspacemethods.Itlocallydecomposesatemporalserieofim-
sectionsoftheimage.Toensureamoreefficientsuppressionof ages into low-rank, sparse and Gaussian components. It is ex-
stellarleakages,thereferencestellarPSFisgenerallyextracted perimentallyshownthattheexoplanetsignalsmostlyremainin
fromdatainwhichtheexoplanetsaretobedetected.Whileat- the sparse term allowing for an improved detection. However,
tenuatingatbestthestellarspeckles,cancellationofpartofthe itisexpectedthatthismethodissensitivetooutliersinthedata
signal fromthe exoplanetsalso occurswhen therotation of the (suchasbadorhotpixels),whicharealsorecoveredinthesparse
field of view is small. This problem is more acute at the vicin- componenttogetherwiththecandidateexoplanets.
ity of the host star since displacements due to the rotation of Finally, the last family of detection methods are based on
thefieldofviewarethesmallestthere.Manyvariantshavebeen a modeling of the stellar PSF based on diffraction modeling.
developed to partially alleviate this problem. A more elaborate PeXalgorithm(Thiébautetal.2016;Devaney&Thiébaut2017)
versionofLOCIcalledTLOCI(Maroisetal.2013,2014)iscur- models the chromatic dependence of speckles. MEDUSAE
rently considered as one of the most advanced standards for method(Ygouf2012;Cantalloube2016)performsaninversion
thedetectionandcharacterizationofexoplanetsbydirectimag- toidentifythephaseaberrationsinaphysicalmodelofthecoro-
ing. The main variation compared to the standard LOCI algo- nagraphanddeconvolvetheimagestorecoverthecircumstellar
rithm is related to the construction of the reference stellar PSF. objectsbyexploitingspectraldiversityandregularizationterms.
A138,page2of27
O.Flasseuretal.:PACO:ExoplanetdetectionbasedonPAtchCOvariances
All of these techniques have several tuning parameters that
often require tuning by trial and error and the intervention of
anexpert,makingtheoptimalitydifficulttoreach.Itisalsoof-
ten observed that the detection maps (or residual images) have
nonstationarybehaviors(typicallywithavariancethatincreases
at smaller angular separations). As a result, it is necessary to
carryoutapre-calibrationsteptocontroltheprobabilityoffalse
alarmwhenusingthesetechniques.Thispre-calibrationistricky
toperforminthepresenceofnumerousbackgroundsources.The
pre-calibrationisalsodependentontheinjectedfakeexoplanet
flux,resultinginalargeprocessingtime.
3. PACO:Exoplanetdetectionbasedon
PAtchCOvariances
The prereduction of the SPHERE datasets uses the Data
Reduction Handling (DRH) pipeline (Pavlovetal. 2008) that
was delivered by the SPHERE consortium with the instrument
to ESO. Additional optimized prereduction tools allow accu-
rate astrometric calibrations (Maireetal. 2016). These prere-
ductions made use of the SPHERE data center (Delormeetal.
2017) which assembles raw images into calibrated datasets by
performing several corrections: bad pixels identification and
interpolation1,thermalbackgroundsubtraction,flat-fieldcorrec-
tion,anamorphismcorrection,truenorthalignment,framecen-
tering,compensationforspectraltransmission,frameselection,
andfluxnormalization.Theoutputsofthesefirststepsaretem-
poralsequencesofcoronagraphicimages(organizedindatasets
with three dimensions hereafter), and of the associated off-axis
PSF.Figure1showsanexampleofascienceframederivedfrom
VLT/SPHERE-IRDIS data and a view of two spatio-temporal
slices extracted at two different locations (along the solid line:
farfromthehoststar,alongthedashedline:nearthehoststar).
Withinthecentralregion(angularseparationsbelow1arcsecat
thewavelengthλ =2.110µm),thesignalisdominatedbystel- Fig.1.SamplefromaVLT/SPHERE-IRDISdataset(HIP72192dataset
larspecklesduet1othediffractioninpresenceofaberrations.At at λ1=2.110µm). Two spatio-temporal slices cut along the solid and
dashedlinearedisplayedatthebottomofthefigure,emphasizingthe
fartherangularseparations,thenoisecomesmainlyfromacom-
spatialvariationsofthestructureofthesignal.
binationofphotonnoisefromthermalbackgroundfluxandde-
tector readout. Observation of the temporal fluctuations in the
wavelengths are available, they can be split and processed in-
spatio-temporalcutsrevealsspatialvariationsofthevarianceof
dependentlyfromeachother.Jointprocessingofmulti-spectral
thefluctuationsbutalsoanevolutionoftheirspatialcorrelations.
datacanbebeneficialbutrequiressomerefinementsthatarebe-
Beyondaccountingfortheaverageintensityofstellarspeckles,
yond the scope of this paper and that will be covered by a fol-
webaseourexoplanetdetectionmethodonanonstationarymod-
lowingpaper.
eling of short-range spatial covariances. To ease the statistical
Science images obtained by high-contrast imaging within
learning of these covariances directly from the data and to ob-
ADI have two components: (i) the signal of interest, due to the
tainadetectionmethodbasedonlocalprocessing,weconsider
presence of exoplanets or background sources in the field of
a decomposition of the field of view into small patches whose
size covers the core of the off-axis PSF. Our detection method view, and (ii) a very strong background, produced by the stel-
lar speckles and other sources of noise, that displays temporal
accountsforthe(spatial)covarianceofthesepatches,hencethe
fluctuations. The motion of the background sources/exoplanets
namePACOforPAtchCOvariance.Inthefollowingparagraphs,
due to the rotation of the field of view is precisely known. An
wedescribeintomoredetailstheproposedmodel,deriveadetec-
exoplanetatsomebidimensionalangularlocationφ atachosen
tion criterion and characterize the photometric and astrometric 0
initial time t is seen at location φ = F(φ ) at time t, where
accuracies. 0 t t 0
F is the geometrical transform (e.g., a rotation) modeling the
t
apparent motion of the field of view between the observation
3.1. Statisticalmodelforsourcedetectionand configurations at time t and time t. We note that a reminder
0
characterization of the main notations used throughout the paper is given in
Table1.
To simplify the notations and the discussion, we have consid-
Since very few sources are within the field of view and
ered datasets at a single wavelength. If observations at several
these sources are faint, we have modeled locally the observed
1 In practice, some bad pixels displaying large fluctuations only on data as the super-imposition of a background (stellar speckles
a few frames remain after this processing. These can be attributed to and noise) and at most one other source, unresolved and lo-
eithertooconservativeselectioncriteriaorevolvingtemporalbehavior catedintheimmediatevincinity(nooverlappingofseveralfaint
ofthedetectorbadpixels. sources).
A138,page3of27
A&A618,A138(2018)
Table1.Reminderofthemainnotations. ByreplacingαinH withthemaximumlikelihoodestimate
1
(cid:98)αobtainedfrom(3)foranassumedinitialpositionφ0,the(gen-
Not. Definitions eralized)likelihoodofeachhypothesiscanbecomparedtoform
the generalized likelihood ratio test (GLRT, see for example
Constants
Kay1998a):
T numberoftemporalframes
N numberofpixelsinatemporalframe logpf({rθk,t(cid:96) −(cid:98)αhθk(φt(cid:96))}k=1:N,(cid:96)=1:T) H≷1 η· (5)
K numberofpixelsinapatch pf({rθk,t(cid:96)}k=1:N,(cid:96)=1:T) H0
Dataquantities
InordertoapplyEqs.(3)and(5)tothedetectionandphoto-
r observedintensity
metricorastrometricestimations,itisnecessarytospecifythesta-
f background
h off-axisPSF tisticalmodelofthebackground.Inmostoftheexistingmethods
forexoplanetdetectionbyADI,somedatapreprocessingisap-
Positions
pliedinordertoreducetheamplitudeofthebackgroundtermand
φ0 2Dangularlocationatachoseninitialtimet0 whitenit(i.e.,lessenitsspatialcorrelations).Suchpreprocessing
Ft(φ0) fieldofviewmotionmodelbetweentandt0 takestheformof(weighted)temporaldifferencesinADI,TLOCI
φt 2Dangularlocationattimetwithφt =Ft(φ0) andANDROMEDA,ahigh-passfiltering,ortheprojectiononto
(cid:98)φt(cid:101) closeston-gridlocationtoφt the subspace orthogonal to the dominant modes identified by a
θk a2Don-gridpixellocation PCA in KLIP and its variations. ANDROMEDA distinguishes
itselfinthatbothitsdetectionandestimationproceduresarede-
rivedfromastatisticalframework(Eqs.(3)and(5)appliedun-
The observed intensity r at the 2D pixel location θ and dertheassumptionofuncorrelatedGaussiannoise).Ratherthan
θk,t(cid:96) k
time t , corresponding to the spatio-temporal index (k,(cid:96)) in the transformingthedatasothatasimplestatisticalmodelcanbeas-
(cid:96)
dataset,canthenbedecomposedintothetwocomponents: sumed(uncorrelatedobservations),weperformednopreprocess-
ingbutaccountforthespatialcorrelationsinourmodel.Given
r =αh (φ )+ f , (1)
θk,t(cid:96) θk t(cid:96) θk,t(cid:96) thatmorethanafewtensofphotonsarecollectedateachdetec-
withα≥0thefluxoftheunresolvedsource(α=0intheabsence torpixel,thePoissondetectionstatisticscanbeapproximatedby
of such source), h (φ )=h(θ −φ ) the off-axis PSF, centered aGaussiandistribution.Theothercontributionsforthetemporal
onthelocationφ θkoftt(cid:96)hesourkceatt(cid:96)timet andsampledatpixel fluctuations(thermalbackgroundnoise,evolutionofthespeckle
location θ , and tf(cid:96) the background at(cid:96)spatio-temporal index patternsduetoevolvingphaseaberrationsorthedecenteringof
(k,(cid:96))whickhaccounθkt,st(cid:96)forstellarspecklesandnoise. thestaronthecoronagraph)willbeconsidered,intheabsenceof
Themajordifficultyinthedetectionofexoplanetsliesinthe amoreprecisemodel,tobeGaussiandistributed.ThisGaussian
fact that the amplitude of the background f is much larger approximation gives closed-form expressions that have a prac-
than the exoplanet contribution αh (φ ) anθdk,tt(cid:96)hat it fluctuates ticalinterestfortheimplementation.Thespatialstructureofthe
fromonetimet toanother.Itisnecθekssatr(cid:96)ytofollowastatistical backgroundismainlyduetothatofthespeckles(relatedtothean-
approachtoaccountforthesefluctuations.Thecollectionofall gularresolution)andtheinterpolationstepsinthedatareduction
otibaslleyrvlaotciaotnesd{arθtk,φt(cid:96)}kw=1i:tNh,(cid:96)a=1fl:Tuixnαthiesptrheesnendceescorfibanedexaospalarnaentdionmi- pthipeeolffin-ea.xMisoPsStFofht(hφe)sheacsoarrceolarteio(nalssaordeeafitnsemdabllyscthaelea.nMguolraerorveesr-,
0
realization with a distribution given by the probability density olution)thatisonlyafewpixelswide.Hence,inEqs.(3)and(5),
function pf ofthebackground: onlypixels(k,(cid:96))ofthedatasetsuchthathθk(φt(cid:96))isnonnegligible
haveanimpactontheestimationordetection.Thesepixelsare
pr({rθk,t(cid:96)}k=1:N,(cid:96)=1:T|α, φt(cid:96))= representedingreeninFig.2.
pf({rθk,t(cid:96) −αhθk(φt(cid:96))}k=1:N,(cid:96)=1:T). (2) We denote the set of the K spatial pixels that form an ex-
tended spatial neighborhood centered at pixel θ , at time t , as
k (cid:96)
Basedonthismodelandunderthehypothesisthatthereare thepatch2 r ∈ RK.Inthefollowing,patchesaresetsofpix-
few sources within the field of view (so that each source can θk,t(cid:96)
elswiththeshapeofadiscretedisk.Iftheradiusofthisdiskis
be considered separately), an unbiased estimation of the flux α chosenlargeenoughtoencompassthecoreoftheoff-axisPSF,
knowingtheinitialangularlocationφ isprovidedbythemaxi-
mumlikelihoodestimator: 0 thecollectionofallpixelvalues{rθk,t(cid:96)}k=1:N,(cid:96)=1:T usedinEqs.(5)
and (3) to detect or estimate the flux of an exoplanet located
(cid:98)α= argαmaxpf({rθk,t(cid:96) −αhθk(Ft(cid:96)(φ0))}k=1:N,(cid:96)=1:T). (3) a(wt iφth0 (cid:98)cφatn(cid:101)tbheercelodsuecsetdpitxoelthteotchoellseucbtipoinxeolflopcaattcihoensφ{tr)(cid:98)φtht(cid:101),at(cid:96)t}(cid:96)c=o1n:T-
tains the exoplanet. Only the joint distribution of these obser-
Wenotethat(cid:98)αimplicitlydependsontheassumedpositionφ0of vations need be modeled. Since the background is different in
thesource.Thedetectionofapointsourceatagivenlocationφ
0 each patch (and possibly varies differently according to time),
canbeformalizedasahypothesistest:
weusedadifferentmodelforeachpatch(localadaptivityofthe

H0 : {rθk,t(cid:96)}(cid:96)=1:T,k=1:N ={fθk,t(cid:96)}(cid:96)=1:T,k=1:N mgrooduenld),psaetechFeisga.t2a.gTivheenrespdautcieadlloncuamtiboenr(Ttypoifcaalvlya,ilfarbolme baafcekw-
H1 : {rθk,t(cid:96)}(cid:96)=1:T,k=1:N =α{+hθ{kf(θφk,tt(cid:96)(cid:96))}}(cid:96)(cid:96)==11:T:T,,kk==11:N:N. (4) tcleaantnisobnteos3aacacnhoduuntnodteruedsdefotaer.mmWupleotricvahalorfisraeatmetoGeasac)uclsoismuianinttsomtnholeydefclootrrorsepdlaeatstiicaorlnibcsoetrhtrheae-t
Literally,thismeansthatunderhypothesisH0 thecollectionof 2 Weuseboldfontstodenotepatcheswhilenormalfontsarereserved
allobservationscorrespondstopurebackground(stellarspeck- toscalarsandangularpositions.
lesandnoise)whileunderhypothesisH1 itisthesuperimposi- 3 Generalizationtospatio-temporalcorrelationsisstraightforwardbut
tionoftheoff-axisPSFwithafluxαandsomebackground. raisesdifficultiesinpracticeduetothelackofsamples.
A138,page4of27
O.Flasseuretal.:PACO:ExoplanetdetectionbasedonPAtchCOvariances
where(cid:98)ρistheshrinkageweight.Thetwoestimators(cid:98)Sand(cid:98)Fare
chosen such that one estimator has very little bias (but a large
variance)andtheotheronehasfewerdegreesoffreedom,hence
amuchsmallervariance,atthecostofanincreasedbias.Thepa-
rameter(cid:98)ρissetinordertoachieveabias–variancetrade-off.As
detailedinAppendixA,wedefined(cid:98)Sasthesamplecovariance
matrixcomputedatlocationθ :
k
(cid:98)Sθk = T1 (cid:88)T (cid:0)rθk,t(cid:96) −(cid:98)mθk(cid:1)·(cid:0)rθk,t(cid:96) −(cid:98)mθk(cid:1)t (8)
(cid:96)=1
with
1 (cid:88)T
(cid:98)mθk = T rθk,t(cid:96), (9)
(cid:96)=1
thesamplemean.Thechosenestimator(cid:98)Fisthediagonalmatrix
formedfromthesamplecovariance:
(cid:2)(cid:98)Fθk(cid:3)ii = T1 (cid:88)T (cid:2)rθk,t(cid:96) −(cid:98)mθk(cid:3)2ii =(cid:2)(cid:98)Sθk(cid:3)ii (10)
(cid:96)=1
and(cid:2)(cid:98)F (cid:3) =0 if i(cid:44) j. (11)
θk ij
Fig. 2. Apparent motion of an exoplanet in an ADI stack of frames. DiagonalelementsoftheshrinkageestimatorC(cid:98)thencorrespond
Green patches contain, in each frame, the exoplanet (i.e., the off-axis to the empirical variance while off-diagonal elements (i.e., co-
PSF) superimposed with the background. The statistical model of a varianceterms)areshrunktoward0byafactor1−(cid:98)ρ,hencethe
backgroundpatchisbuiltlocallybasedonobservedpatchesatthesame name.
locationbutatdifferenttimes(setofredpatches). Theshrinkagefactor(cid:98)ρislocallyadaptedinanunsupervised
data-drivenway.InAppendixAwederiveitsexpression4:
aidnriegstrcthiobenudstiiidostenrrioebdfuettioaocnbhepbsfataoctkfigsatrlilocbuaanllcdykpginardotceuhpn.edPnpaditexcnhetle.ssR,awattehdteihffreetnhreamnntodtdieemfileends- (cid:98)ρ(cid:16)(cid:98)Sθk(cid:17)= tr(cid:16)(cid:16)T(cid:98)S2θ+k(cid:17)1+(cid:17)(cid:16)ttrr2(cid:16)(cid:16)(cid:98)S(cid:98)S2θ(cid:17)k(cid:17)−−(cid:80)2K(cid:80)iK(cid:104)=(cid:98)S1(cid:104)(cid:98)S(cid:105)2θk(cid:17)(cid:105)2ii · (12)
onlythedistributionofbackgroundpixelsinthepatchesofinter- θk i=1 θk ii
est.Thisdistribution pf({f(cid:98)φt(cid:101),t(cid:96)}(cid:96)=1:T)ismodeledasaproductof Remaining bad pixels after the prereduction (see introduc-
Gaussian distributions defined over each patch that would con- torydiscussioninSect.3)takeverylargevaluesatsomedates.
tainanexoplanetifthatexoplanetwasinitiallylocatedatφ These aberrant values often impact exoplanet detection algo-
0
rithms. By learning the local covariance of the patches, these
(cid:89)T largefluctuationsareaccountedforwithinourmodel(largepixel
pf({f(cid:98)φt(cid:101),t(cid:96)}(cid:96)=1:T)= N(f(cid:98)φt(cid:101),t(cid:96) |m(cid:98)φt(cid:101),C(cid:98)φt(cid:101)), (6) variance). We show in Sect. 5 that these bad pixels do not lead
(cid:96)=1 tofalsealarms.
withN(·|m ,C )theprobabilitydensityfunctionofthemul-
(cid:98)φt(cid:101) (cid:98)φt(cid:101)
tivariate Gaussian (mean m and covariance C) that describes
3.3. Unbiasedestimationofthebackgroundstatistics
backgroundpatchescenteredatpixel(cid:98)φ(cid:101).
t
To account for the fact that observed patches {rθk,t(cid:96)}(cid:96)=1:T do not
contain pure background but also, in the presence of an exo-
3.2. Statisticallearningofthebackground
planet, the signal due to the exoplanet, estimation of the back-
Inourmodel,wehaveconsideredacommonmean m andco- ground statistics and of the exoplanet flux can be alternated,
θk
avtaraiagnivceenCpθikxfeolrloaclaltTionbaθckk.gBraocukngdropuantcdhpeastc{hfeθks,t(cid:96)a}r(cid:96)e=1n:Totcdeinrteecrteldy isttearratitniogn,fraofmteraanflinuixtia(cid:98)αl(ng)uheasss boefepnuerestibmacaktegdrofuonrdt.heAetxtohpelann-teht
available, only observed patches {rθk,t(cid:96)}(cid:96)=1:T can be used to es- assuming background statistics C(cid:98)(n) (as explained in Sect. 3.4),
timate the mean mθk and covariance Cθk. Moreover, the limited thebackgroundstatisticscanbeimprovedbycomputing
number T of temporal frames makes a direct estimation of C
bneWdltsehhryoamaaertilanistnhf-yskdodket.sururaSaesingtvtoaotdeeemomgntcthiipheoeteslsavsfehtoacacrcararbiainmopsnaskveanhbaanbocrergcifieilneeacirektneoyfosagcnatfgtieucsomeialtedaonbaaireryattpirhozlwpeaetradrunartoC(cid:98)incs(orasieancnitenshlhgeomkedc(ApeKudbaafipsleearlptnfisayyeettmchtnd&iahueeddebnstiEaexytrbprle.(AdtetgiheaIfueednarlTnnac2atfdroo0<oinnz0rFoKcvua8iueegtt))isrdx.ohwo.npcAerSiroltao.ltehm1,onvbθ)odaa---k.  (cid:2)(cid:98)F(θ(cid:98)mC(cid:98)(cid:98)(cid:98)Snρk+(θ(θ((nnnnkk1++++)(cid:3)1111i))))i=====(cid:98)ρ((cid:2)TT11(cid:98)1S(cid:0)(cid:98)S(cid:80)(cid:80)θ(−nk(θ+nT(cid:96)T(cid:96)k(cid:98)ρ==+1·(11)n1(cid:3)(cid:0)(cid:0)(cid:0))i+rrr(cid:1)iθθθ1kkk),,,)ttt(cid:96)(cid:96)(cid:96)(cid:98)S−−−(n(cid:98)(cid:98)(cid:98)ααα+(((1nnn))))+hhhθθθ(cid:98)ρkkk((((nφφφ+ttt(cid:96)(cid:96)(cid:96)1))))(cid:1)−−(cid:98)F((cid:98)(cid:98)mmn+(θ(θnnkk1++).11))(cid:1)(cid:1)t (13)
binationoftwoestimators(cid:98)Sand(cid:98)F
C(cid:98)=(1−(cid:98)ρ)(cid:98)S+(cid:98)ρ(cid:98)F (7) 4orlWaregeernftohracneo(cid:98)ρnθek.tobeintherange[0,1]byclippingvaluesbelowzero
A138,page5of27
A&A618,A138(2018)
Thisiterativeprocedureisappliedtoprovideunbiasedestimates Ourassumptionsamounttoapproximatingthedistributionofthe
ocofntthreibfluutixonooffatnheexeoxpolpalnaente,tanfotetrbdeesteucbttiroanc.teIdndpereiodr,tsohocoumldptuhte- tTehrmevra(cid:98)rφita(cid:96)(cid:101)n,t(cid:96)ce−o(cid:98)mf(cid:98)bφt(cid:96)(cid:101)thbeyreafocreentseirmedpliGfiaeusstsoi:anofcovarianceC(cid:98)(cid:98)φt(cid:96)(cid:101).
(cid:96)
ing the statistics of the background, the mean would contain a
fraction of the PSF of the exoplanet (1/Tth of the flux of the Var{b }≈ h (φ )t·C(cid:98)−1 ·h (φ )=a . (19)
exoplanet if the exoplanet is visible only in one of the patches (cid:96) (cid:98)φt(cid:96)(cid:101) t(cid:96) (cid:98)φt(cid:96)(cid:101) (cid:98)φt(cid:96)(cid:101) t(cid:96) (cid:96)
centered on pixel θk). The covariance would then encode that Thus an estimator of the variance of b= (cid:80)T(cid:96)=1b(cid:96) is (cid:80)T(cid:96)=1a(cid:96)=a
therearesignificantvariationsinthebackgroundintheformof andthestandarddeviationof(cid:98)α=b/acanbeestimatedby:
an appearing and disappearing PSF. The subsequent estimation
√
ofthesourcefluxwouldbepenalizedbytheseerrorsonthemean (cid:98)σα=1/ a. (20)
and covariance of the background. A few iterations of Eq. (13)
and exoplanet flux re-estimations corrects the statistics of the Obviously,thefluxofthesourceisnecessarilypositive,we
backgroundandpreventsfrombiasingtheestimationofthepho- denoteby(cid:98)α+thefluxobtainedunderapositivityconstraint5(see
tometrybyanerroneousstatisticalmodelofthebackground. Thiébaut&Mugnier2005;Mugnieretal.2009):
Experimentshaveshownthattheunbiasedestimationofthe
fluxoftheexoplanetisnotbeneficialintheexoplanetdetection (cid:98)α+ d=efargmaxpf({rθk,t(cid:96) −αhθk(Ft(cid:96)(φ0))}k=1:N,(cid:96)=1:T)
α≥0
ptmhheaaksdeeis.stWtrhibheuiltseieoitnttiinmogfptrohofeveatsedsthteetienscittgihonenalatobhfsreethnsheceoelxdoofp(eslaxeneoeptAlsa,pnipetetc,nhdwainhxgiceAhs =max((cid:98)α,0)=max(cid:80)(cid:16)(cid:80)TT(cid:96)=1ab(cid:96),0(cid:17), (21)
for details) more difficult. Therefore, we computed the back- (cid:96)=1 (cid:96)
groundstatistics(cid:98)mθk andC(cid:98)θk basedonEqs.(7)–(12)tocompute We assume that, if positivity is imposed, the covariance of (cid:98)α+
a detection map (see Sect. 3.5), and alternate Eq. (13) with a canbeapproximatedby(cid:98)σα,inEq.(20).
re-estimationofthefluxforthephotometricandastrometrices-
timationsofadetectedsource(seeSects.3.4and3.6).
3.5. Detectionofexoplanets
Under our multivariate Gaussian model of the background, the
3.4. Estimationofthefluxofanexoplanet
generalized likelihood ratio test (Eq. (5)) takes the simplified
Under our multivariate Gaussian model of the background, the form
minaitxiaiml ulomcaltiikoenlihφo0odhaesstaimsaimtoprlgeiveexnprinesEsiqo.n(3()sefeorfaonr aesxsaummpelde (GLRT) (cid:16)(cid:80)T(cid:96)=1b(cid:96)(cid:17)2 H≷1 η, (22)
Kay1998b): (cid:80)T(cid:96)=1a(cid:96) H0
(cid:80)T b witha andb definedaccordingtoEqs.(15)and(16).
(cid:98)α= (cid:80)T(cid:96)=1a(cid:96) , (14) As(cid:96) discus(cid:96)sed in Thiébaut&Mugnier (2005), Smithetal.
(cid:96)=1 (cid:96) (2009) and Mugnieretal. (2009), it is beneficial to enforce a
positivity constraint on the flux α in the detection test, in other
with words, to use the estimate (cid:98)α+ to derive the GLRT expression,
leadingto
a = h (φ )t·C(cid:98)−1 ·h (φ ) (15)
(cid:96) (cid:98)φt(cid:96)(cid:101) t(cid:96) (cid:98)φt(cid:96)(cid:101) (cid:98)φt(cid:96)(cid:101) t(cid:96) (cid:16)max(cid:16)(cid:80)T b ,0(cid:17)(cid:17)2
(GLRT+) (cid:96)=1 (cid:96) H≷1 η. (23)
and (cid:80)T(cid:96)=1a(cid:96) H0
(cid:16) (cid:17)
b(cid:96) = h(cid:98)φt(cid:96)(cid:101)(φt(cid:96))t·C(cid:98)−(cid:98)φ1t(cid:96)(cid:101)· r(cid:98)φt(cid:96)(cid:101),t(cid:96) −(cid:98)m(cid:98)φt(cid:96)(cid:101) , (16) AlesntntootethdeintesMtugnieretal. (2009), the test (Eq. (23)) is equiva-
wsaonhuderrwceehawotse(esruceebcnpatilexlretlih)sal(cid:98)toφcha(cid:101)(cid:98)tφ,itot(cid:96)(cid:101)hn(eφφtn(cid:96)t(cid:96))esadareemnsptolpteeidxseotlhvteeoroφaffp-.aatxcihsoPfSKFpfioxrelas (S/Ntest) (cid:113)(cid:80)(cid:80)T(cid:96)=T1ba(cid:96) =(cid:98)σ(cid:98)αα HH≷10 τ (24)
t(cid:96) t(cid:96) (cid:96)=1 (cid:96)
The maximum likelihood estimator(cid:98)α given in Eq. (14) de- √
pends linearly on the data, hence its variance can be easily de- whenη≥0andwithτ= η.Thistestcanbeinterpretedasthe
rived.Notingthat(cid:98)α=a/bwitha=(cid:80)T(cid:96)=1a(cid:96) andb= (cid:80)T(cid:96)=1b(cid:96),we S/N (cid:98)α/(cid:98)σα of the estimation of the (unconstrained) flux of the
obtain: source α. We note that with our definition, the S/N is a signed
quantitythatisnegativewhenever(cid:98)α<0.
Var{(cid:98)α}=Var{b}/a2 (17) The test (S/Ntest) in Eq. (24) is attractive because the test
value (cid:98)α/(cid:98)σα linearly depends on the data. Under our Gaus-
sinceaisdeterministic.Thevarianceofbisthesumofthevari- sian model for the data, the S/N (cid:98)α/(cid:98)σα is thus also Gaussian
ancesoftheb(cid:96) termssincetheyaremutuallyindependent(cor- distributed (with unit variance) which simplifies the statistical
respondingtodifferenttemporalframes).Fromitsexpressionin analysis of the test in terms of false alarm rate or detection
Eq.(16),thevarianceofb(cid:96) isgivenby: probability.Thisanalysisiscarriedoutinthefollowing.
Inthehypothesistest(Eq.(4)),weconsideredanhypotheti-
Var{b(cid:96)}=h(cid:98)φt(cid:96)(cid:101)(φt(cid:96))t·C(cid:98)−(cid:98)φ1t(cid:96)(cid:101)·Cov{r(cid:98)φt(cid:96)(cid:101),t(cid:96) −(cid:98)m(cid:98)φt(cid:96)(cid:101)} calinitiallocationofthesourceφ0.Todetectallsourceswithin
×C(cid:98)−1 ·h (φ ). (18) 5 The neg-log-likelihood is a quadratic function of α, the minimum
(cid:98)φt(cid:96)(cid:101) (cid:98)φt(cid:96)(cid:101) t(cid:96) underpositivityconstraintisthusobtainedbysimplethresholding.
A138,page6of27
O.Flasseuretal.:PACO:ExoplanetdetectionbasedonPAtchCOvariances
the field of view and locate their positions, the test (S/Ntest)
should be evaluated at locations φ sampled over the field of
0
view.Byrefiningthesamplingofthefieldofview,theoff-axis
aPcScFusraht(cid:98)eφ,t(cid:96)(cid:101)a(tφtth(cid:96))ebceottsetromfaatclhaergsethrecodmatapuatnadtitohneaelsetiffmoartt.eS(cid:98)αaimspmlionrge
ofthefieldofviewisdiscussedinSect.4.3.
3.6. Statisticalguarantees
3.6.1. Distributionofthedetectioncriterion
Figure3illustratesaS/Nmap(cid:98)α/(cid:98)σα computedusingthePACO
method on a VLT/SPHERE-IRDIS dataset (obtained on the
HIP 72192 star as described in more detail in Sect. 5) consist-
ing of 96 temporal frames in which the two real known faint
pointsourcesaremasked.VisualinspectionoftheS/Nmap(top
ofFig. 3)indicatesthatthedetectioncriterionisstationaryover
thefieldofview,evenclosetothecoronagraphicmask.Theem-
pirical distribution of the S/N values (bottom of Fig. 3) shows
thattheS/Ndistributioncanbeconsideredascentered(sincethe
empirical mean is 0.01) and approximately reduced (since the
empiricalstandarddeviationis0.93)Gaussian.Thisdistribution
passessuccessfullytheLillieforsnormalitytest(Lilliefors1967)
atthe5%significancelevel.Thesmalldiscrepancywiththeoret-
icalmodelmaybeduetothecorrelationsinthedata.However,
the hypothesis that the variance of the S/N under H is equal
0
to one is conservative in the sense that the probability of false
alarmisslightlyoverestimated.Owingtothehomogeneousdis-
tributionoftheS/Ncriterionacrossthefieldofview,thissmall
biascouldbeeasilyfixedbyrescalingtheS/Nbyasinglefactor
thatisempiricallyestimated6.
UnderthehypothesisH ,theS/Nfollowsacenteredandre-
0
duced Gaussian law whatever the angular separation. It is thus
possibletoapplyauniquethresholdτtotheS/Nmapsandob-
tainaconsistentdetectionperformanceatallangularseparations
(i.e.,aconstantfalsealarmrate).TheaccordancebetweentheS/N
and a Gaussian distribution makes it possible to directly assess
thefalsealarmrate,probabilityofdetection,photometricandas-
trometricaccuracieswithoutpost-processingand/orresortingto
Monte Carlo methods (injection of fake exoplanets in the data)
in contrast to several state-of-the-art methods (Lafreniereetal.
2007;Cantalloubeetal.2015;Ruffioetal.2017).
3.6.2. Probabilitiesoffalsealarmandoftruedetection
The probability of false alarm (PFA) is the probability that the
test(Eq.(24))yieldsH whileH isactuallytrue: Fig. 3. S/N map in absence of object (the two real known faint point
1 0
sourcesdenotedFPS1andFPS2aremasked)anditscorrespondingem-
PFA(τ)=Pr((cid:98)α/(cid:98)σα >τ|H0) piricaldistributionusingtheHIP72192datasetdescribedinSect.5at
(cid:90) +∞ 1 (cid:32) x2(cid:33) λ2=2.251µm.
= √ exp − dx=1−Φ(τ), (25)
τ 2π 2 The minimum amplitude α such that the probability of de-
where Φ is the cumulative distribution function of the standard tection be at least equal to some prescribed value PD when the
normal distribution. The conventional contrast at 5σ ensures a detection threshold τ is set according to a given PFA level can
probabilityoffalsealarmequaltoPFA(5)=2.87×10−7. becomputedbyinvertingandcombiningEqs.(25)and(26):
The probability of detection (PD) is the probability that the
test(Eq.(24))decidescorrectlyforadetection
PD(τ,α)=Pr((cid:98)α/(cid:98)σα >τ|H1) α=(cid:0)Φ−1(1−PFA)−Φ−1(1−PD)(cid:1)×(cid:98)σα. (27)
=(cid:90) +∞ √1 exp−(cid:2)x−α/(cid:98)σα(cid:3)2 dx For example, at τ=5, a probability of detection of 50% is
τ 2π 2 achievedforα=5(cid:98)σα andaprobabilityofdetectionof80%for
=1−Φ(τ−α/(cid:98)σα). (26) α=5.84(cid:98)σα. Figure 4 illustrates the S/N distribution under the
twohypothesisaswellastheprobabilitiesoffalsealarmandof
6 Nosuchcorrectionwasperformedinthefollowingresults. detection.
A138,page7of27
A&A618,A138(2018)
has little impact on the estimation of the other two parameters
(low absolute correlation coefficients) except for certain local-
izedareasofthefieldofview.TheresultingCramér-Raobounds
can be usefully considered to evaluate the error on the esti-
mated parameters. For example, for the two real faint point
sources located around the HIP 72192 star (their positions are
symbolized by a circle in Fig. 5, the products α · δ and
x0
α·δ areclosedto0.7×10−5 pixel.Itmeansthatatacontrast
y0
α=10−5, the sources can be located with an accuracy of about
0.7pixel.
Fig.4.S/NdistributionunderH (inblue)andH (inred).Thehatched
0 1 4. ImplementationdetailsofPACO
areaisequaltotheprobabilityoffalsealarm(PFA)whilethefilledarea
isequaltotheprobabilityofdetection(PD).
This section is devoted to the description of the implementa-
tion of the PACO algorithm presented in Sect. 3. A simpli-
fied and faster version for the detection step is also described.
3.6.3. Astrometricaccuracy
Thisfastversioncanbeusefultoconductpre-analysesinlarge
An (asymptotically) unbiased estimator of the position of the surveys.
source is provided by the maximum likelihood position which
is the solution of a non-convex optimization problem. In prac-
tice,thisproblemcanbesolvedbyexhaustivesearchatafinite 4.1. Optimalpatchessize
resolution refined by a local optimization as for example pro-
The patches considered in the PACO algorithm define the
posed by Soulezetal. (2007) in digital lensless microscopy for
characteristic size of the areas in which the statistics of the
thedetectionofparametricobjectsspreadinavolume.Theas- backgroundfluctuationsaremodeled.Sincethecoreoftheoff-
trometric accuracy, that is, the accuracy on the spatial location
axis PSF is close to circular symmetry, circular patches are
of the detected objects can be statistically predicted using the used. Their size obeys a trade-off: on the one hand, the larger
Cramér-Rao lower bounds (CRLBs) which represents the min-
thepatches,themoreenergyfromthesourceiscontainedinthe
imal variance of any unbiased estimator. The CRLBs are ob- patcheswhichimprovestheS/N;ontheotherhand,learningthe
tained from the diagonal of the inverse of Fisher information
covariance of larger patches requires more samples (i.e., more
matrix (Kendalletal. 1948). Using a parametric model of the
temporalframes).
off-axis PSF h and noting again that the collection of patches
In practice, since the sources to be detected are faint com-
{r } located at the angular position θ is described by a
θk,tl l=1:T k paredtothelevelofstellarspecklesandtheirtemporalfluctua-
multi-variate Gaussian process N(·|mθk,Cθk), it is possible to tions, only the core of the off-axis PSF is necessary to perform
derive the CRLBs from the observed intensity model given in the detection and a patch size corresponding to twice the off-
Eq. (1). In the following, Ω={α,x ,y } (with φ ={x ,y } the
0 0 0 0 0 axis PSF full width at half maximum (FWHM) appears to be a
angularpositionofanexoplanetattimet )denotesthevectorof
0 good compromise. A more precise (and automatic) determina-
parameters from which the CRLBs are computed. For a given
tion of the optimal size of the patches with respect to number
angular position θ , the Fisher information matrix IF can be
k θk oftimeframesandthestructureofthebackgroundcorrelations
expressedas canbecarriedoutbyMonteCarlosimulations.Falseplanetsare
randomly injected into the data and receiver operating charac-
∂αh (Ω)t ∂αh (Ω)
[IF(Ω)] = θk ·C−1· θk , (28) teristic(ROC)curvesrepresentingthedetectionprobabilitiesas
θk i,j ∂Ωi θk ∂Ωj a function of the false alarm rate are constructed for different
inwhichthetermαh representstheoff-axisPSF(anisotropic patchsizes.ThepatchsizemaximizingtheareaundertheROC
θk curveisthenadopted.Figure6givesROCcurvesforthreepatch
Gaussian can typically be used as a continuous model to
sizes and shows that the choice K=49 pixels per patch maxi-
compute the derivatives). It follows that the standard devi-
ation δ on the estimation of the parameter vector Ω is mizes the area under the curve and is therefore the best com-
θk promise for performing the detection on the HIP 7192 dataset
givenby
(with 4.5 pixels FWHM) for fake injected exoplanets of flux
(cid:113) (cid:104) (cid:105)
[δ ] = [IF(Ω)−1] (29) α ∈ 10−6;10−5 .ThenumberofpixelsperpatchK needstobe
θk i θk i,i determined only once for a given instrument since the FWHM
For simplicity, we denote by δ={δ ,δ ,δ } the spatial oftheoff-axisPSFvariesonlymarginallyfromoneobservation
maps(obtainedforallangularpositionsθα ofx0they0fieldofview) toanother.
k
representing the accuracy on the parameters Ω={α,x ,y }.
0 0
AppendixBgivestheanalyticalexpressionofδasafunctionof
4.2. ThePACOalgorithm:detailedimplementation
the Fisher information matrix terms. We note that δ is not de-
α
pendentonαwhileδ andδ areproportionaltoα−1.Figure5 AsexplainedinSect.3,withinthePACOpipeline,thedetection
x0 y0
gives a view of δ , δ and δ as well as the correlation coef- andestimationstepsareperformedbytwoverysimilarschemes.
α x0 y0
ficients ρ , ρ and ρ between estimated parameters ob- First, the detection is performed on the whole field of view us-
tainedonαtxh0eHαIyP0 72192xd0ya0tasetatλ =2.110µm.Asexpected, ingthePACO detectionprocedure.Thisstepproducesamap
1
the estimation error of Ω = {α,x ,y } increases near the host ofS/Nwhichisstatisticallygroundedandwhichcanbedirectly
0 0
star. The figure also emphasizes that the Cramér-Rao bounds thresholded at a controlled false alarm rate. The estimated flux
{δ ,δ ,δ }decreasequicklywiththeangularseparation.More- ishoweverbiasedbythepresenceoftheexoplanetinthecollec-
α x0 y0
over, a small estimation error on one of the three parameters tionofpatchesusedtomodelthebackgroundfluctuations.Thus,
A138,page8of27
O.Flasseuretal.:PACO:ExoplanetdetectionbasedonPAtchCOvariances
Fig.5.Upperline:minimalstandarddeviation(Cramér-Raolowerbounds)δ={δ ,δ ,δ }ontheestimatedparametersΩ={α,x ,y }(δ and
α x0 y0 0 0 x0
δ arenormalizedbytheinverseofthefluxαoftheexoplanetandexpressedinpixelsunit).Middleline:magnificationoftheareanearthehost
stya0r.Lowerline:coefficientsofcorrelationρ={ρ ,ρ ,ρ }betweentheestimatedparametersΩ={α,x ,y }.Thecomputationisperformed
ontheHIP72192datasetatλ =2.110µm.Thepαox0sitiαoyn0soxf0yt0hetworealfaintpointsourcesinthedataset(0wh0icharenotvisibleonthesemaps
1
becauseitdoesnotactofdetectionmaps)arerepresentedbyacircle.
a different procedure named PACO estimation is launched in motionoftheexoplanetovertimeislimitedandthesourcefluxis
asecondsteponeachdetectedsourceinordertorefinetheflux large.
estimation. This latter procedure provides unbiased flux esti- Since the data analysis of large surveys is a crucial issue
mates(i.e.,noaposterioriMonteCarlo-basedbiascorrectionis in astronomy, we also propose a simplified and faster version
necessary). of the PACO detection procedure which is summarized by
Algorithm1summarizesthePACO detectionprocedureas Algorithm 2. Compared to applying Algorithm 1 to a given set
described more formally in Sect. 3. Step 1 consists of form- ofassumedinitialsourcepositionstocomputeamapofthede-
ing the collection of patches P(cid:96) on which the statistics of the tectioncriterion,thefastversionhasacomputationalburdenre-
background must be learned. These patches are all centered at duced by a factor at least equal to the number T of temporal
the same position φt(cid:96) where the source would be at time t(cid:96), frames. The acceleration relies on the precomputation of terms
assuming it was initially at position φ0. In Step 2, the back- that appear multiple times in the sums of Eqs. (15) and (16)
ground statistics, i.e. the empirical mean (cid:98)m and the regular- when considering all possible source locations. Computations
ized covariance C(cid:98), are computed based on Eqs. (7)–(12). Step of the background statistics are thus recycled. The S/N map
3 then forms the numerator and denominator of the test statis- is obtained in a second step by interpolating the precomputed
tics by accumulating values a and b defined by Eqs. (15) terms denoted by c and d in Algorithm 2 in order to align the
(cid:96) (cid:96)
and (16). In Algorithm 1, the background statistics are com- field of view at all times according to the transform F which
t(cid:96)
putedassumingnoexoplanetispresent(i.e.,hypothesisH )and is, in general, a rotation. Such interpolations result in a low-
0
are thus biased. This is especially notable when the apparent pass filtering of the criterion map which slightly degrades the
A138,page9of27
A&A618,A138(2018)
the S/N map at level τ should be photometrically characterized
using the statistically unbiased PACO estimation procedure
summarizedinAlgorithm3.AsdiscussedinSect.3.2,thiscan
be done by alternating between an estimation of the flux of the
exoplanet and of the statistics of the background7. The result-
itng((cid:96)e=sti1m:Ta)tec(cid:98)αorfroerspaosnodusrctoeltohceamteidniamtφumt(cid:96)=oFftt(cid:96)h(φe0f)oollnowthiengfracmoset
(cid:96)
function:
Algorithm2:fast PACO detection–Fastcomputation
of the S/N values at initial 2D angular positions G of an
unresolvedsource.
Input:UniformgridGof2Dangularpositions.
Input:Spatio-temporaldataset{rθk,t(cid:96)}k=1:N,(cid:96)=1:T.
Output:S/NmapatallinitialpositionsinG.
(cid:46)Step1.Precomputeterms:
forallφ ∈Gdo
0
Fig.6.Influenceofpatchsize:ROCcurvesforK={13,49,113}pixels (cid:46)Step1a.Buildthecollectionofpatchescentered
ineachpatch.TheROCcurvesareobtainedbyinsertingfakeexoplanets atφ0:
offluxα∈(cid:104)10−6;10−5(cid:105)ontheHIP72192datasetatλ1=2.110µm. P ←{r(cid:98)φ0(cid:101),t(cid:96)}(cid:96)=1:T
(cid:46)Step1b.Learnthecorrespondingbackground
Algorithm1:PACO detection–ComputationoftheS/N s(cid:110)tatisticsan(cid:111)dpre-computeterms:
values at initial 2D angular positions G of an unresolved (cid:98)m(cid:98)φ0(cid:101),C(cid:98)(cid:98)φ0(cid:101) (Eqs.(7)–(12))
source. w ←C(cid:98)−1 ·h (φ )
Input:SetGofinitial2Dangularpositions. c(cid:98)φ0(cid:101)←wt(cid:98)φ0(cid:101)·h(cid:98)φ0(cid:101)(φ0)
Input:Spatio-temporaldataset{rθk,t(cid:96)}k=1:N,(cid:96)=1:T. fo(cid:98)φr0(cid:101)(cid:96)=1:(cid:98)φT0(cid:101)do(cid:98)φ0(cid:101) 0
Output:S/NmapatallinitialpositionsinG.
forallφ ∈Gdo d(cid:98)φ0(cid:101),t(cid:96) ←wt(cid:98)φ0(cid:101)·(r(cid:98)φ0(cid:101),t(cid:96) −(cid:98)m(cid:98)φ0(cid:101))
0
a←0 (cid:46)Step2.ComputetheS/Nvalues:
b←0
forallφ ∈Gdo
for(cid:96)=1:T do a←00
(cid:46)Step1.Extracttherelevantpatches: b←0
P(cid:46)pφatS(cid:96)t(cid:96)ct=e←hpFest2{(cid:96)(ir.nφ(cid:98)Lφ0Pte(cid:96))(cid:101)a,t(cid:96)(cid:96)r(cid:48)n:}(cid:96)t(cid:48)h=1e:Tbackgroundstatisticsfromthe for(cid:96)φabt=(cid:96)←←=1abF:++t(cid:96)T(cdφd(cid:98)(cid:98)0φφo)√tt(cid:96)(cid:96)(cid:101)(cid:101)(,tφ(cid:96)(tφ(cid:96))t(cid:96)) ((iinntteerrppoollaattiioonnooffdc))
{(cid:98)m(cid:98)φt(cid:96)(cid:101),C(cid:98)(cid:98)φt(cid:96)(cid:101)} (Eqs.(7)–(12)) S/N(φ0)←b/ a (Eq.(24))
(cid:46)Step3.Updateaandb:
w←C(cid:98)−1 ·h (φ )
a←a+(cid:98)φtw(cid:96)(cid:101)t·h(cid:98)φt(cid:96)(cid:101) (φt(cid:96) ) (Eq.(15)) C(α)=(cid:88)T (cid:26)T log(det(C(cid:98) (α)))
S/Nb(φ←0)b←+bw/t√·(ar(cid:98)(cid:98)φφtt(cid:96)(cid:96)(cid:101)(cid:101),t(cid:96)t−(cid:96) (cid:98)m(cid:98)φt(cid:96)(cid:101)) ((EEqq..((1264)))) +(cid:96)=1tr(cid:18)C(cid:98)−1(α)·(cid:18)Wφ(cid:12)t(cid:96) (cid:88)T u (α)·ut (α)(cid:19)(cid:19)(cid:27), (30)
φt(cid:96) φt(cid:96),t(cid:96)(cid:48) φt(cid:96),t(cid:96)(cid:48)
(cid:96)(cid:48)=1
detection performance of PACO (see Sect. 5). The complexity inwhichWi,j=1ifi= j,Wi,j=1−(cid:98)ρelsewhere,{i, j}∈J1;KK2,(cid:12)
standsforentrywisemultiplication(i.e.,Hadamardproduct),and
opfret-hcealcfauslattvioenrssiotenpg(Sivteenpi1n).ADlegnoorittihnmgb2yiNs dthoemniunmatebderboyfathne- uθk,t(cid:96)(α) = rθk,t(cid:96) − α · hθk,t(cid:96) − (cid:98)mθk(α). The expression of C(α)
in Eq. (30) comes from the co-log-likelihood under a Gaussian
gular positions φ to process, this step requires N×T products
0 assumption, where matrix W is introduced in order to bias the
of vectors with K elements as well as the resolution of N lin-
covarianceestimatetowardadiagonalcovariance(i.e.,inorder
ear systems of size K×K. For example, the application of this
toreplacethemaximumlikelihoodcovarianceestimategivenby
fastalgorithmrequiresapproximatelytwominutestoprocessa
thesamplecovariancewiththeshrinkageestimatordescribedin
dataset made of 96 temporal frames of size 1024×1024 pixels
Sect.3.2).
versusapproximatelythreehoursforthecompletealgorithmus- Figure 7 illustrates the cost function C(α) that is mini-
ing a basic parallelization done in MatlabTM on 24 cores (pro- mizedbyouralternatingestimationscheme.Thiscostfunctionis
cessorIntelXeonE5–46170at2.90GHzand K=49pixelsper unimodalinthetestedrangeofcontrasts(from5×10−6to10−1)
patch).
Once the detection step is performed by Algorithm 1 or 7 Ajointestimationofthefluxoftheexoplanetandofthebackground
Algorithm 2, the potential detections obtained by thresholding statisticscouldalsobeperformedbyhierarchicaloptimization.
A138,page10of27