TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations
M. Samland, J. Bouwman, D. W. Hogg, W. Brandner, T. Henning, M. Janson
AAstronomy & Astrophysics manuscript no. trap © ESO 2020November 26, 2020
TRAP: A temporal systematics model for improved direct detectionof exoplanets at small angular separations
M. Samland , , , J. Bouwman , D. W. Hogg , , , , W. Brandner , T. Henning , and M. Janson Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germanye-mail: [email protected] International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD) Department of Astronomy, Stockholm University, Stockholm, Sweden Center for Cosmology and Particle Physics, Department of Physics, New York University, 726 Broadway, New York, NY 10003,USA Center for Data Science, New York University, 60 Fifth Ave, New York, NY 10011, USA Flatiron Institute, Simons Foundation, 162 Fifth Ave, New York, NY 10010, USAReceived ; accepted
ABSTRACT
Context.
High-contrast imaging surveys for exoplanet detection have shown giant planets at large separations to be rare and it is ofparamount importance to push towards detections at smaller separations, the part of the parameter space containing most planets. Theperformance of traditional methods for post-processing of pupil-stabilized observations decreases at smaller separations, due to thelarger field-rotation required to displace a source on the detector in addition to the intrinsic di ffi culty of higher stellar contamination. Aims.
Our goal is to develop a method of extracting exoplanet signals, which improves performance at small angular separations.
Methods.
A data-driven model of the temporal behavior of the systematics for each pixel can be created using reference pixels at adi ff erent position, under the condition that the underlying causes of the systematics are shared across multiple pixels, which is mostlytrue for the speckle pattern in high-contrast imaging. In our causal regression model, we simultaneously fit the model of a planetsignal “transiting” over detector pixels and non-local reference lightcurves describing a basis of shared temporal trends of the specklepattern to find the best fitting temporal model describing the signal. Results.
With our implementation of a spatially non-local, temporal systematics model, called TRAP, we show that it is possible togain up to a factor of six in contrast at close separations ( < λ/ D ) compared to a model based on spatial correlations between imagesdisplaced in time. We show that the temporal sampling has a big impact on the achievable contrast, with better temporal samplingresulting in significantly better contrasts. At short integration times (4 seconds) for β Pic data, we increase the SNR of the planet bya factor of four compared to the spatial systematics model. Finally, we show that the temporal model can be used on unaligned datawhich has only been dark and flat corrected, without the need for further pre-processing.
Key words.
Planets and satellites: detection – Methods: data analysis – Techniques: high angular resolution – Techniques: imageprocessing
1. Introduction
The field of direct observation of extrasolar planets has seentremendous progress over the last ten years, both in terms ofobservational capabilities, as well as high-contrast imaging dataanalysis. Especially the recent advent of instruments dedicated tohigh-contrast observations like SPHERE (Spectro-PolarimetricHigh-contrast Exoplanet REsearch; Beuzit et al. 2019), GPI(Gemini Planet Imager; Macintosh et al. 2014), and CHARIS(Gro ff et al. 2015), as well as the development of more sophis-ticated observational strategies and post-processing algorithmshave paved the way to make this a unique technique to studyatmospheres and orbital characteristics of exoplanets directly.However, giant planets and substellar companions at large or-bital separations ( (cid:39) −
10 au) have been proven to be rare bymultiple large direct imaging surveys (e.g., Brandt et al. 2014;Vigan et al. 2017; Nielsen et al. 2019). Pushing our detection ca-pabilities towards smaller inner-working angles is therefore oneof our primary goals in order to tap into a parameter space regimein which planets are more abundant (Nielsen et al. 2019; Fernan-des et al. 2019) and which overlaps with indirect detection tech- niques, such as astrometric detections with
Gaia (e.g., Casertanoet al. 2008; Perryman et al. 2014) and long-term radial velocitytrends (e.g., Crepp et al. 2012; Grandjean et al. 2019).Detecting planets is intrinsically more di ffi cult the smallerthe separation becomes, because the speckle background ishigher, and fewer independent spatial elements are available forstatistical evaluation of the significance of the detected signal(s).However, this intrinsic problem is further exacerbated becausemost algorithms used to model the stellar contamination obscur-ing the planet signal (the systematics) require strict exclusioncriteria determining which data can be used to construct a data-driven model of these systematics while preventing the incorpo-ration of planet signal.In this work, we present a novel algorithm designed to im-prove performance at small angular separations compared toconventional algorithms for pupil-tracking observations. Espe-cially as coronagraphs become more powerful and allow us toprobe smaller inner working angles (IWA ), the algorithmic per- The IWA is defined as the angular separation at which the corona-graphic transmission crosses 50%. Article number, page 1 of 21 a r X i v : . [ a s t r o - ph . E P ] N ov & A proofs: manuscript no. trap formance will become more and more important. We achievethis goal by building a data-driven, temporal systematics modelbased on spatially non-local data. This model replaces the tem-poral exclusion criterion with a less restrictive spatial exclusioncriterion, and allows using all frames in the observation sequenceregardless of angular separation. This in turn allows us to makebetter use of the temporal sampling of the data and to create asystematics model that is sensitive to variations on all sampledtime scales. We further employ a forward model of the com-panion signal, fitting both the systematics model and the planetmodel at the same time to avoid over-fitting and biasing the de-tection.In Section 2 we will give an overview of systematics model-ing in high-contrast imaging data, with a focus on the spatial sys-tematics modeling approach traditionally used in pupil-trackingobservations. In Section 3 we motivate our non-local, temporalsystematics approach. Our causal regression model exploits thefact that pixels share similar underlying causes for their system-atic temporal trends. We show how we can apply such a modelto pupil-tracking high-contrast imaging data. Section 4 discussesthe data used to demonstrate the performance of our algorithmTRAP. Section 5 shows and compares the results obtained. Weend with a discussion and a future outlook on how the algorithmcan be further improved in Section 6, and provide a summary ofour conclusions in Section 7.
2. Systematics modeling in pupil-trackingobservations
The main challenge in high-contrast imaging is distinguishinga real astrophysical signal from light of the central star, the so-called speckle halo. This stellar contamination is generally or-ders of magnitude brighter than a planetary companion objectin raw data, and furthermore can locally appear the same as agenuine point source.In order to separate the astrophysically interesting signal ,such as a planet, from the systematic stellar noise background,we need to model these systematics. Due to the complexity of thesystematics, this is usually done using a data-driven approach,e.g. using the data itself as a basis to construct a model. In or-der for this to work, there must be one or more distinguishingproperties between the planetary signal and the contaminatingsystematics (the star light). These distinguishing properties arealso referred to as diversities.For exoplanet imaging the primary diversity is spatial reso-lution, i.e. a discernible di ff erence in the position of the com-panion and host star signal on the sky. The most commonlyused high-contrast imaging strategy enhances this distinguish-ing spatial property by inducing a further time-dependence ofthe measured companion signal by allowing the field-of-view(FoV) to rotate over the course of an observation sequence, whilethe speckle halo remains stationary. This mode of observation iscalled pupil-tracking mode and is the basis for algorithms suchas angular di ff erential imaging (ADI, Marois et al. 2006).In this paper we show that instead of building a spatialspeckle pattern model for each image, we can build a temporallightcurve model for each pixel, which provides an alternativeand novel way of reducing data taken in pupil-tracking mode. For simplicity we may simply refer to a planet as the astrophysicalsignal of interest, being the focus of this work, but the same appliesfor more massive objects and other signals, whether it be an extendedobject or a point source.
Due to the field-of-view rotation, temporal variation is in-duced in the planetary signal as measured by the detector, thuscreating a signal that varies in both space and time. Because thereare two ways in which the planet signal di ff ers from the stellarsystematics, we have the freedom to build our systematics modelon either, the spatial correlation between images or the temporalcorrelations between pixel time series (lightcurves). A mixturemodel of both presents another possibility which is not exploredin this work.The di ff erence between the two pure approaches can be un-derstood based on whether we treat the data as being either: a) aset of image vectors, to be explained in a basis constructed fromimage vectors taken at other times (spatial systematics model),or b) a set of time series vectors, to be explained in a basis con-structed from other time series vectors taken at other locations(temporal systematics model). Before we go into details aboutthe proposed temporal systematics approach in Section 3, wewill first summarize the commonly used spatial systematics ap-proach and discuss its drawbacks. Following the invention of roll deconvolution (or roll anglesubtraction ) for space-based observatories (Müller & Weigelt1985), classical ADI (Marois et al. 2006) for ground-based ob-servatories using pupil-tracking, and its various evolutions whichstarted employing optimization and regression models in theform of
Locally Optimized Combination of Images (LOCI, e.g.,Lafrenière et al. 2007; Pueyo et al. 2012; Marois et al. 2014;Wahhaj et al. 2015), virtually all papers and algorithms focusedon using the spatial systematics approach. In this approach thetraining data is typically taken from the same image regionwhere the planet signal is located (as implied also by the nam-ing of LOCI), and the spatial similarity between an individualimage and the training set of images displaced in time is usedto remove the quasi-static speckle pattern which covers the areawhere the planetary companion signal is located at that time.This is the case for most commonly used algorithms, regard-less of the implementation details of the model construction. Insome cases the regression is not performed on the images them-selves, but a representation of the image vectors in a lower di-mensional space using dimensionality reduction techniques suchas principal component analysis (PCA / KLIP, Amara & Quanz2012; Soummer et al. 2012). Implementations in publicly avail-able pipeline packages, follow this trend, e.g. PynPoint (PCA,Amara & Quanz 2012; Stolker et al. 2019), KLIP (Soummeret al. 2012; Ru ffi o et al. 2017), VIP (Gomez Gonzalez et al.2017). The ANDROMEDA algorithm is similarly based on aspatial model using di ff erence images to suppress the stable fea-tures in the speckle pattern (Cantalloube et al. 2015) togetherwith a maximum-likelihood approach to forward modeling theexpected companion signal in the di ff erence images. Another re-cent approach called patch covariance (PACO, Flasseur et al.2018) is trying to statistically distinguish a spatial patch co-moving with a planet signal from the properties of a stationarypatch.Two notable exception to this spatial approach is the wavelet-based temporal de-noising approach (Bonse et al. 2018), which,however, does not attempt to create a causal regression modelof the systematics. It applies a temporal filter and pre-conditionsthe data before applying a spatial systematics approach. Addi-tionally, the STIM-map approach (Pairet et al. 2019) adjusts thedetection map based on the temporal residuals left after a spatialmodel is applied. Article number, page 2 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations
Let us assume we have a data cube D with a time series of imagedata d ik ∈ R M × T , where M is the number of pixels and T thenumber of frames in the time series. For simplicity in the nota-tion we will from here on have x and t denote discretely sampledpoints in space and time, without further use of indices when themeaning is clear from context. All data obtained with imaginginstruments is discretely sampled.Any measurement d in the data D can be described in the func-tional form evaluated at discrete points D ( x , t | N ) = S ( x , t ) + g ( x , t | N ) + (cid:15) ( x , t ) (1)where S is the planet signal contribution, g denotes the func-tional form of the systematics a ff ecting this measurement, and (cid:15) is the stochastic noise (e.g., photon noise). N here is a stand-infor all hidden parameters that may be responsible for causing thesystematics (e.g., wavefront phase residuals after AO correction,temperature, wind direction / speed).In a data driven approach to modeling the systematics, ourgoal is to determine g ( x , t | N ) based on how other data express g which are a ff ected by the same underlying hidden parameters N . If we had a perfect model for g simply by subtracting it wecould detect and determine the planet signal S up to a stochasticuncertainty term. However, we have to guarantee that we do notincorporate contributions from S into our systematics model g ,otherwise we will attenuate the signal of interest, an e ff ect calledself-subtraction . For this reason, we divide the data into twosubsets of measurements. Those that are a ff ected by signal froma planet above a chosen threshold will be called Y , and thosethat are not will be called X . The union of both sets representsall sampled position and time combinations D and both sets aredisjoined. Y = { ( x , t ) ∈ D | S ( x , t ) > thr } , X = { ( x , t ) ∈ D | S ( x , t ) ≤ thr } , X ∪ Y = D , X ∩ Y = ∅ . (2)The idea is thus that g ( x , t | N ) for all ( x , t ) ∈ Y can be approxi-mated using a combination of elements in X that are a ff ected bythe same underlying causes of systematics N . The most directapproach to this type of problem is to assume a linear regres-sion model. However, since the data itself is drawn from a multi-dimensional space (space and time), a choice has to be madeas to the basis in which the linear regression is performed. Tra-ditionally, image vectors (or more often vectors of sub-imageregions) have been used as the basis of the linear regression(e.g., all LOCI and PCA variants previously mentioned). Forone-dimensional representations we may use a simplified vec-tor notation, where, for example image vectors, can be denotedas d t ( x ), i.e. as a data vector of pixel values for a set of positions x at the time t . The image space is an intuitive basis when wethink of the speckle noise as a relatively stable (quasi-static) pat-tern that is part of the (coronagraphic) stellar PSF (point spreadfunction). The PSF is intuitively thought of as a spatial construct. We use a single index to denote the position (pixel) in the 2D image. If, indeed, the model is subtracted. Otherwise it constitutes a form ofover-fitting of the systematics model. The threshold chosen does not matter for the formalism. In practice,however, it usually relates to some defined exclusion criterion basedon the total flux overlap of a companion signal and training data. Thecompanion PSF extends over the whole image at a very low flux level,therefore a meaningful cut-o ff has to be chosen for defining the extentof the signal. Let us now formulate, in general terms, the spatial approachas it is implemented in LOCI-like algorithms, in which we aretrying to fit an image data vector in terms of other image vectorstaken at a di ff erent time.We start by defining image regions P k Y of pixels that are sig-nificantly a ff ected by S ( x i , t k ) at a given t k as subsets of the com-plete image space of all pixels P k Y ⊂ PP k Y = { x i | S ( x i , t k ) > thr } , k ∈ { , , . . . , T } . (3)This corresponds to the area in an image covered by the com-panion PSF above a certain flux threshold at a given time. Giventhat the signal position changes over time, we can now definesets of times T k X in the observation sequence for which the planetsignal is not a ff ecting this image region T k X = { t l | P k Y ∩ P l Y = ∅ , k , l ∈ { , , . . . , T }} (4)This is the set of all times at which the companion has moved farenough to not overlap with the area it occupied at a certain time.From these times we can build a linear model that does not in-corporate planet signal. We can then say that the systematics andplanet flux in a particular image region P k Y and for any particulartime t k can be estimated using a linear least squares regressionsolvingarg min ω k , α k , l ∀ l (cid:88) x ∈ P k Y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) D ( x , t k ) (cid:124) (cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32) (cid:125) data − ω k ˆ S ( x , t k ) + (cid:88) t l ∈ T k X α k , l D ( x , t l ) (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) planet signal + systematics model (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (5)where ω k corresponds to the contrast of the planet signal esti-mated at time t k , α k , l are the weight coe ffi cients of the linearmodel at time t k , and ˆ S is a model of the planet signal. In thespatial representation the planet model takes the form of an im-age of the companion PSF ˆ s t k ( x ) with x ∈ P k Y . This approachassumes that the speckle pattern is either su ffi ciently stable (upto a scaling factor) and / or recurrent over the course of the ob-servation (e.g., finite probability of jitter to revisit the same posi-tion, recurring Strehl ratio or observing conditions in general). Inother words, we assume we have enough measurements d t ( x | N )with t ∈ T k X and x ∈ P k Y to probe the underlying distributionsof unobserved causal factors N to allow us to reconstruct a spe-cific instance of the systematics function g . In the spatial repre-sentation g corresponds to the speckle pattern g t k ( x ) in the im-age region. To a large extent this assumption is well founded,as evidenced by the success of the LOCI and spatial PCA-basedfamily of algorithms. How, in detail, the α -coe ffi cients are opti-mized and how the planet contrast ω is finally estimated variesstrongly depending on the implementation (e.g., Lafrenière et al.2007; Pueyo et al. 2012; Marois et al. 2014; Wahhaj et al. 2015),and is one of the main distinguishing factors between the largenumber of published algorithms. Here two closely related prob-lems are the problem of collinearity and overfitting of the planetsignal. Considering that many of the image vectors used in thelinear model share common features, one solution is to use reg-ularization (e.g., damped LOCI, Pueyo et al. 2012). The mostpopular approach is to use a truncated set of principal compo-nent images instead of using the image vectors d t l ( x | N ) directly,as done in the above Eq. 5. This representation in a lower dimen-sional space simultaneously acts to prevent overfitting.Whether the planet signal is fitted simultaneously to the sys-tematics model, as done in this work; fitted after subtracting a Article number, page 3 of 21 & A proofs: manuscript no. trap systematics model, as has been the case for most early works;or subtracted before fitting the systematics model (e.g., Galicheret al. 2011) is another key di ff erence in various implementations.Additionally, implementations vary in how they define the ge-ometry and size of the image regions, and may further use twodi ff erently defined regions based on whether the region is usedfor optimizing the α -coe ffi cients (usually larger than the actualregion of interest) or subtracting the model, which is another wayof reducing overfitting. Lastly, we note that reference star di ff er-ential imaging (RDI, eg., Smith & Terrile 1984; Lafrenière et al.2009; Gerard & Marois 2016; Ruane et al. 2017; Xuan et al.2018; Bohn et al. 2020) is also based on the above describedspatial systematics model formalism, with the di ff erence that thereference images are taken from observations of di ff erent hoststars, consequently lacking the same planet signal contribution S by default. This paper will only be concerned with individualpupil-stabilized observations. The above described assumptions impose a high requirement onthe overall stability of the instrument and observing conditions.In modern instruments a large part of the systematics can indeedbe said to be quasi-static and correlate on timescales of minutesto hours (e.g., Milli et al. 2016; Goebel et al. 2018), making thisapproach applicable to these stellar noise contributions.However, the main drawback of the spatial modeling ap-proach is the need to prevent contamination of the systematicsmodel with astrophysically relevant signal by using a temporalexclusion criterion (usually defined as a protection angle). Its se-lection is a trade-o ff between excluding more frames to preventself-subtraction, and losing information on the speckle evolutiondue to exclusion of highly correlated frames close in time. As hasalready been noted early on (Marois et al. 2006), this trade-o ff gets gradually worse the closer one approaches the central star,because the same field-of-view rotation corresponds to smallerand smaller physical displacement on the detector between sub-sequent frames. At small separations a large fraction of the ob-servation sequence will be excluded from the training set (seeAppendix A).At small separations these exclusion criteria cover timespans of the order of the linear decorrelation timescale for rel-atively stable instrumental speckles (Milli et al. 2016; Goebelet al. 2018). Even at large separations the temporal exclusionwill still be on the order of minutes. This means that usinga spatial approach makes it intrinsically di ffi cult to model theturbulence-induced short-lived speckles with an exponential de-cay timescale of a few seconds ( τ = ffi cult to account for. Traditionally, obtainingsimultaneous training data for modeling the systematics requiresspectral di ff erential imaging (SDI, Rosenthal et al. 1996; Racineet al. 1999) or polarimetric di ff erential imaging (PDI, Kuhn et al.2001).However, in this work we demonstrate how it is possible toapply a causal, temporal systematics model to monochromaticpupil-tracking observations, using simultaneous but non-localreference data. This allows us to avoid the harsh temporal exclu-sion criterion used for preventing self-subtraction, by introduc-ing a less strict spatial exclusion criterion of pixels a ff ected bysignal, solely determined by the total FoV rotation independentof angular separation between the star and exoplanet.
3. TRAP: A causal, temporal systematics model forhigh-contrast imaging
Analogous to the spatial approach as described in Eq. 3, 4 and 5,we can define the temporal approach in which all image vectorsat a specific time, become time series vectors at a specific loca-tion by first constructing a subset of times at which a specificpixel x i is a ff ected by planetary signal and construct a referenceset of other pixels x j , that, in the same subset of times does notcontain planet signal. This would be a one-to-one translation ofLOCI to a temporal approach, one that is local in time, ratherthan space. Let us first construct the set of pixels a ff ected by theplanet signal at any point in time t of the observation sequence: P Y = ∪ k P k Y . (6)For each a ff ected pixel x i ∈ P Y we will construct a training setof una ff ected other pixels P i X : T i Y = { t | S ( x i , t ) > thr } , ∀ x i ∈ P Y , P i X = { x j | T i Y ∩ T j Y = ∅} , i , j ∈ { , , . . . , M } , (7)and subsequently, we can build our model similar to Eq. 5 forany x i ∈ P Y arg min ω i , α i , j ∀ j (cid:88) t ∈ T i Y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) D ( x i , t ) (cid:124) (cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32) (cid:125) data − ω i ˆ S ( x i , t ) + (cid:88) x j ∈ P i X α i , j D ( x j , t ) (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) planet + systematics model (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (8)where the data to be fitted takes the form of a (partial) time series(lightcurve) d x i ( t | N ) of one pixel x i ∈ P Y for all times t ∈ T i Y at which the pixel is a ff ected by a planet signal above a certainthreshold. In the temporal representation the planet model takesthe form of a positive transit lightcurve ˆ s x i ( t ) for t ∈ T i Y . Thesystematics g in the temporal representation correspond to tem-poral trends g x i ( t ) with t ∈ T i Y and are reconstructed in a basisset of reference time series that are una ff ected by planet signalduring the same time d x j ( t | N ). In the literature this process isusually referred to as de-trending of lightcurves . The fitting ofthe systematics model can therefore be achieved in a mathemat-ically analogous way when swapping the spatial and temporalaxes. However, in our implementation we loosen the restrictionon temporal locality, such that we always look at the completetime series and are left with only two relevant sets: the set ofpixels that are at any point in time a ff ected by planet signal P Y and those that are at no point a ff ected P X P Y = { x | ∃ t with ( x , t ) ∈ Y}P X = { x | ∀ t with ( x , t ) ∈ X}P X ∪ P Y = P , P X ∩ P Y = ∅ , (9)and the model is defined over all times t ∈ T such that for any x i ∈ P Y , Eq. 8 becomesarg min ω i , α i , j ∀ j (cid:88) t ∈ T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) D ( x i , t ) (cid:124) (cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32) (cid:125) data − ω i ˆ S ( x i , t ) + (cid:88) x j ∈ P X α i , j D ( x j , t ) (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) planet + systematics model (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (10) Article number, page 4 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations
Intuitively it may not be obvious why we can use non-local timeseries to create a systematics model for the time series of a spe-cific pixel, unless we consider the causal structure of the system-atics that confounds the signal in the data. While pixels physi-cally separated enough on the detector do not directly “talk toeach other” and are, at first order, independent measuring de-vices, they are influenced by the same underlying causes thatgenerate the systematics in the first place. At the most basic levelthe primary cause of systematics, the contamination by the stel-lar signal, is the (coronagraphic) PSF of the star itself whichchanges with observing conditions and changes in the instru-ment. Of course, depending on the detailed underlying cause,the area of shared influence can be spatially di ff erent, e.g. winddriven halo e ff ects (Cantalloube et al. 2018, 2020) cannot be sta-tistically inferred from reference pixels una ff ected by the sameunderlying cause.A detailed mathematical description of modeling systemat-ics in time series data based on shared underlying causes canbe found in Schölkopf et al. (2016), and gives the statisticalbackground and justification for Eq. 8 and 10. The abovementioned work shows that by using a regression model of asu ffi ciently large number of “half-siblings”, i.e. measurementsthat do not share the signal but the same causal relationshipof the systematics d x j ( t | N ) of one pixel x j ∈ P X , we canreconstruct a specific instance of the systematics function g x i ( t | N ). Further, we can include a model of the planet signalitself in the regression, such that the overall regression explainsthe data in P Y up to a stochastic noise term (see Eq. 1), whichremains the fundamental limit of the data. This fundamentallimit is a combination of photon noise, detector read-out noise,and flat fielding uncertainties. In practice, we will still be limitedby imperfection in the systematics model, because only a finitenumber of regressor lightcurves are available, and some causesof systematics may not be perfectly shared among them.The mathematical idea of using a causal time series regres-sion models, as described above, has been employed very suc-cessfully for transit observation using the Kepler spacecraft (e.g.,Wang et al. 2016a). The benefit of simultaneously including aforward model of the expected transit shape in the regressionmodel has been demonstrated in detail in Foreman-Mackey et al.(2015).
The situation is very similar for high-contrast imaging data.Causes of systematics, with the exception of detector artifacts(e.g., bad / warm / dead pixels, flat field), usually influence eitherthe image globally (e.g., most atmospheric e ff ects, Strehl ratiovariations), or a significant region of the detector (e.g., wind-driven halo, (mis-)alignment of the coronagraph). Even slowlyevolving changes in the quasi-static speckle pattern caused bythe instrument usually are not strictly confined to one region ofthe detector (e.g., speckle patterns can and do display point sym-metries).A detailed study of the objectively optimal choice ofreference pixels to capture most shared underlying systematiccauses is beyond the scope of this work, but multiple heuristicchoices can be made: 1) preference for pixels at a similarseparation (same underlying modified Rician speckle intensitydistribution, Marois et al. 2008), 2) pixels at same positionangle inwards and outwards of reduction area (position angledependent e ff ects, e.g., wind driven halo, Cantalloube et al. O ff s e t ( p i x e l ) Fig. 1.
Example of reference pixel selection for one assumed planetposition. The signal here is assumed to move through the position( ∆ x , ∆ y ) ≈ (0 ,
28) pixel above the star at the midpoint of the obser-vation sequence. This position is marked by a pink cross and representsthe center of the reduction area P Y . The reference pixels P X are shownin white. This example is based on the 51 Eri observation’s parallacticangles, the annulus is 7 pixel wide. ff ected by planetsignal. Varying the regressor-selection parameters around theseheuristics did not significantly impact our final results. Figure 1shows an example for the reference pixel (regressor) selectiongeometry chosen in this work for an assumed companion thatpasses north of the host star at the midpoint of the observation.The annulus used has a width of 7 pixels and the PSF usedincludes the first Airy ring and has a diameter of 21 pixels.In our implementation called TRAP (Temporal ReferenceAnalysis of Planets), we use, a) a linear model with a quadraticobjective function (assuming Gaussian uncertainties ), b) theprincipal components of the regressor pixel lightcurves, not thepixels themselves (principal component regression), c) we fiteach a ff ected pixel individually instead of all a ff ected pixels atthe same time for computational reasons, d) we fit the planetmodel and the systematics model simultaneously in order toprevent over-fitting, similar to the approach by Foreman-Mackeyet al. (2015). The above choice of using principal componentsreduces the collinearity of the regressors used for the systemat-ics model and transforms our model into an orthogonal basis. Italso allows us to forgo regularization in favor of truncating theprinciple components after a certain number. Furthermore, we https://github.com/m-samland/trap This is a simplified assumption and reality is more complex in termsof the underlying distributions (e.g., Marois et al. 2008; Pairet et al.2019). Additionally, noise is generally heteroscedastic and correlated.We mitigate this issue by empirically rescaling the derived uncertaintiesat the end. Article number, page 5 of 21 & A proofs: manuscript no. trap do not use any spatial pre-filtering of the data as is commonlydone in many implementations to remove structures correlatedon larger scales than the expected planet signal from the data(e.g., as done in Cantalloube et al. 2015; Ru ffi o et al. 2017). Wealso do not pre-select / discard any frames as bad frames in thisstudy.Given a data cube of the observation sequence D = d ik ∈ R M × T , where M is the number of pixels and T is the number offrames in the time series, as well as the known parallactic an-gles for each frame, we know the temporal evolution of a com-panion signal’s position. The algorithm looks as follows for one assumed relative position θ ( ∆ RA , ∆ DEC) of a companion pointsource on-sky and its corresponding trajectory on the detector:1. We generate a forward model of the planet signal at assumedposition θ for the time series by embedding the model PSFin an empty data cube at the appropriate parallactic angle foreach exposure to obtain the set of pixels P Y a ff ected at anypoint in time during the observation. Conversely, all pixels P X not a ff ected at any time by planet signal are potentialreference pixels.Additionally we obtain the corresponding planet modellightcurve ˆ s x i ( t ) for every x i ∈ P Y . We use the unsaturatedPSFs obtained without the coronagraph directly before or af-ter the sequence, but artificially induced satellite spots couldalso be used. The PSF has been adjusted to the exposure timeand filter of the science exposures and is set to zero beyondthe first Airy ring.2. Instead of using all time series from una ff ected pixels P X to build the systematics model, we construct the training setwith additional desired constraints (cid:101) P X ⊂ P X as previouslydescribed (see Fig. 1). The number of reference pixels in thereference set will be called R = | (cid:101) P X | .3. For each pixel x i ∈ P Y :(a) We want to simultaneously optimize the coe ffi cients ofour systematics and planet model.We will first show how this done starting from Eq. 10,which uses the reference lightcurves from the data di-rectly. The minimization problem can be solved as a sys-tem of linear equations of the form A θ T C i − A c i = A θ T C i − b i (11)because Eq. 10 is linear in ω i and α i , j . A θ is the designmatrix, also called regressor matrix, containing the refer-ence lightcurve vectors constituting the terms of the lin-ear model to be fitted. The same A θ is used for all pixels x i for a given planet position θ . c i = ( α i , , . . . , α i , R , ω i ) T is the vector of model coe ffi cients, b i = d x i ( t ) is thelightcurve data vector of pixel x i that is to be fitted, and C i is the covariance associated with the data vector to befitted. The uncertainties of the time series are in generalheteroscedastic (not uniform in time), but uncorrelated,such that the o ff -diagonal elements are zero.(b) The design matrix consists of columns containing theadditive terms of the overall model we want to fit. Assuch the part of the design matrix for Eq. 10 that de-scribes only the systematics model consists of columnsthat contain the reference lightcurves d k j = d jkT suchthat j | x j ∈ (cid:101) P X , j ∈ { , , . . . , R } , and k ∈ { , , . . . , T } .Since this would only account for the systematics, weadd two additional columns, one containing a constanto ff set to be fitted, and one containing the lightcurve vec-tor ˆ s x i ( t ) describing the companion signal at pixel x i over time . This (preliminary) design matrix A θ looks as fol-lows A θ, k j = d k j , if j ≤ R , if j = R + S k , i , if j = R + ∈ R T × ( R + . (12)(c) However, this results in an overdetermined and highlycollinear problem, because there are typically more refer-ence lightcurves than there are samples in time ( R (cid:29) T ),and the lightcurve share many common features. There-fore some form of regularization or dimensionality re-duction is required. We opt for representing them in alower dimensional space using the principal componentsof the lightcurves instead of the reference lightcurvesthemselves. This is implemented using a singular valuedecomposition (SVD) . We first center the data robustlyby subtracting the temporal median of each lightcurve d k j . This centered data matrix can then be factorizedinto U Σ SVD V T , where U corresponds to a T × T ma-trix, Σ SVD corresponds to a T × R diagonal matrix withnon-negative real numbers, and V T corresponds to an R × R matrix. Because the data d k j is real, U and V T are real and orthonormal matrices. In the further contextwe will be using U , the columns of which are called left-singular vectors and form an orthonormal basis in whichthe lightcurves d k j can be described. They are orderedin decreasing order of explained variance. We truncatethe number of components to be fitted to remove spuri-ous signals from the systematics model and reduce over-fitting. The number of components used in the algorithmis specified by the fraction f of the maximum number ofcomponents N max = T , such that N used = f N max with f ∈ [0 ,
1] and N used ∈ N . N used is always rounded to theclosest integer for any chosen f . The fraction of compo-nents f is a user parameter that can be set. The regressorsfor the temporal systematics are then given by U kl , where l ∈ { , , . . . , N used } . The final design matrix will be: (cid:101) A θ, kl = U kl , if l ≤ f N max , if l = f N max + S k , i , if l = f N max + ∈ R T × ( f N max + . (13)(d) The vectors of the design matrix are mostly orthogo-nal and we can robustly invert the equation analogouslyto Eq. 11 using linear algebra operations to obtain thevector of model coe ffi cients and their associated uncer-tainty. The coe ffi cients and their associated covariancesare given by c i = [ (cid:102) A θ T C i − (cid:102) A θ ] − [ (cid:102) A θ T C i − d x i ] Σ i = [ (cid:102) A θ T C i − (cid:102) A θ ] − . (14)In our current implementation, the covariance matrix ofthe data C is assumed to be uncorrelated in time (zero It is possible to include additional source of information on temporalbehavior as columns here, such as auxiliary data on atmospheric condi-tions. In practice we use the compact SVD to improve computational per-formance. We note that the SVD can become ill-defined for the inner-most regions, when the temporal sampling of the data is very high. Inthis case, one should ensure that the number of regressor pixels is nottoo small compared to the number of temporal samples.Article number, page 6 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations for o ff -diagonal elements), as the residuals after fittingthe model are approximately white. The variance for thecoe ffi cients are therefore given by σ i = diag( Σ i ) . (15)The results shown throughout the paper assume identi-cal uncertainties. However, known uncertainties for eachpixel (e.g., photon noise, read-noise) can be passed tothe routine, which is generally recommended if the vari-ance of the input data is known. Taking into account pho-ton noise puts less weight on data with higher specklenoise (e.g., pixels closer to the star, located on the spi-der, frames taken under worse conditions) further in-creasing the robustness of the reduction. The elements ω i = c i , N used + and σ ω, i = σ i , N used + correspond to the con-trast of the companion and its respective variance by con-struction. A graphical example for the system of equationEq. 11 to be solved is shown in Fig. 2. The result of the fitof the time series are shown in Fig. 3. The above exampleis for the central pixel of the signal-a ff ected detector areashown in Fig. 1, injected with a relatively bright signal(10 − contrast) for demonstration purposes.4. After obtaining ω i and σ ω, i for all x i ∈ P Y , we remove sig-nificant outliers whose median of the residual vector devi-ates more than a set threshold in robust standard deviationsfrom the rest (in our implementation 3 robust standard devi-ations). This usually removes a low single digit number ofpixels ( < ff ected) pixels expected to be present inan area the size of a typical reduction region. The incurredinformation loss is minimal. We call this new subset of pix-els (cid:101) P Y ⊂ P Y .5. Finally we take the average of the planet contrast coe ffi cientsweighted by their respective uncertainties over all remainingpixels ω θ = (cid:80) i ω i σ − ω, i (cid:80) i σ − ω, i , (cid:101) σ θ = (cid:118)(cid:117)(cid:116)(cid:88) i σ − ω, i − , i | x i ∈ (cid:101) P Y , (16)to obtain a single contrast and preliminary uncertainty valuefor this on-sky position θ . Eq. 16 assumes the pixel residu-als to be uncorrelated in space, which in practice is not thecase. The preliminary uncertainty (cid:101) σ θ needs to be scaled withan empirical normalization factor, the derivation of which isdiscussed below, to account for simplified assumptions aboutthe covariance structure of the data and spatial correlationsbetween the residuals.Iterating over a grid of possible positions allows us to constructa contrast and pseudo-uncertainty map, i.e. the contrast ω θ andits preliminary uncertainty (cid:101) σ θ given the positions of the assumedcompanion relative to the central star. From these we construct apseudo signal-to-noise (SNR) map. The above uncertainties arecomputed under the simplified assumption that all measurementsare Gaussian and independent, which may not accurately reflectthe reality of high-contrast imaging data. Multiple studies haveshown residuals after post-processing with high-contrast imag-ing pipelines, although significantly whitened, are not strictlyGaussian and independent (Marois et al. 2008; Cantalloube et al.2015; Pairet et al. 2019). The most direct solution would be toaccount for these e ff ects in the likelihood function used, but in We use the same spatial sampling as the pixel scale of the instrument.
Fig. 2.
Example for Eq. 11 showing the first five principal componentlightcurves and the planet model. The constant o ff set term ( c N used + ) isnot shown in the design matrix. This example shows a pixel with abright (10 − contrast) injected planet signal based on 51 Eri b data. Thedata corresponds to the central pixel in the injected planet’s trajectoryshown in Fig. 1 at the location marked by the pink cross. The principalcomponents were determined from the reference pixels P X shown inwhite .practice this proves to be challenging since the uncertainty dis-tributions are heteroscedastic and depend on the observing con-ditions, instrument, location in the image, and time. Therefore,we apply the same solution to this problem as in ANDROMEDA(Cantalloube et al. 2015), and subsequently in FMMF (forwardmodel matched filter, see Ru ffi o et al. 2017), and empiricallynormalize the pseudo SNR map using the azimuthal robust stan-dard deviation of the SNR map itself computed in 3-pixel wideannuli as a function of separation. We denote these final uncer-tainty values as σ θ = (cid:101) σ θ σ annulus ,θ . (17)Any reference to SNR in this work refers to the normalizedSNR ω θ /σ θ . This normalized SNR map will also be referred toas the detection map. Any a-priory known companion objects orbackground star sources will be masked out when deriving thenormalization to avoid biasing the result.
4. Datasets used for demonstration
We demonstrate our TRAP algorithm on two real datasets ob-tained with the extreme-AO (SAXO; Fusco et al. 2014) fed in-frared dual-band imager and spectrograph (IRDIS; Dohlen et al.2008) of SPHERE (Beuzit et al. 2019) at the VLT. Both datasetswere obtained as part of the SHINE (SpHere INfrared survey forExoplanets, Chauvin et al. 2017) large survey. The observationwere obtained in pupil-tracking mode and use the apodized Lyotcoronagraph (Soummer 2005) with a focal-mask diameter of 185milli-arcsec.The two datasets are described in Table 1 and Fig. 4 and 5show the pre-processed temporal median image for the 51 Eri-dani and β Pictoris observations, respectively.Both datasets were pre-reduced using the SPHERE DataCenter pipeline (Delorme et al. 2017), which uses the Data Re-duction and Handling software (v0.15.0, Pavlov et al. 2008) andadditional custom routines. It corrects for detector e ff ects (dark,flat, bad pixels), instrument distortion, aligns the frames, and cal-ibrates the relative flux. Article number, page 7 of 21 & A proofs: manuscript no. trap C o un t s ( A D U ) datacomplete model reconstructed systematicsplanet model fit Fig. 3.
Result of the fit for the pixel shown in Fig. 2. Shown is thereconstructed model for the time series with and without including thecoe ffi cient ω i describing the amplitude of the companion signal.
1) The first test case is the directly imaged planet in the51 Eridani system (Macintosh et al. 2015). The data used hasbeen published in Samland et al. (2017). This dataset was ob-tained in the K1 and K2 bands simultaneously, but we focus ourdiscussion on the K1 channel in which the companion is visible.The sequence was not taken under ideal conditions. It exhibitscommon problems encountered in high-contrast imaging, suchas a strong wind-driven halo as seen in Fig. 4.These e ff ects make for a realistic test case, and allow usto demonstrate the algorithm performs well under adverse andchanging conditions.2) The second test case uses IRDIS H2 data of the β Pic sys-tem (Lagrange et al. 2009) published in Lagrange et al. (2019)taken in continuous satellite spot mode, i.e. using sine-wavemodulations on the deformable mirror to induce four satellitespots in all frames of the sequence that can be used to deter-mine the center accurately. This also permits measurement ofthe satellite spot amplitude variations as proxy for the planetPSF model’s flux modulation over time. The purpose in includ-ing this dataset is two-fold. Firstly, the exposure time is shortfor an IRDIS sequence (4 s), allowing us to test our time-domainbased algorithm on a dataset with better temporal sampling, ap-proaching the half-life time of fast-decaying speckles ( τ = .
5. Results
We perform a direct comparison of results between TRAP andthe o ffi cial IDL implementation of the ANDROMEDA algo-rithm (Cantalloube et al. 2015). It is always challenging to di-rectly compare a new approach with existing pipelines due tonumerous di ff erences in implementation, a wide range of pos-sible reduction parameters, test data taken in various conditionsand observing modes, as well as di ff erences in statistical evalu-ation of the outputs. No comparison will ever be complete. Forour work, we have chosen ANDROMEDA as the most viablerepresentation of the spatial systematics modeling class of al-gorithms, because both TRAP and ANDROMEDA implementa likelihood-based forward modeling approach of the compan-ion signal. This allows us to directly evaluate and compare theoutputs in a fair way, such that di ff erences in implementation donot impact the statistical interpretation of the results. Addition-ally, ANDROMEDA has been shown to compare well with otherestablished PCA and LOCI-based pipelines (e.g., Cantalloubeet al. 2015; Samland et al. 2017).For our tests and comparison we use the exact same normal-ization method for the detection map and contrast curve compu-tation for both TRAP and ANDROMEDA outputs. For the com-putation of the empirical normalization and the contrast curveswe mask the position of the planet with a mask size radius of r mask =
15 pixels in the output maps ( ∼
180 mas).We show the results for TRAP for a range of principal com-ponent fractions f (see Eq. 13). Our representative choice is f = .
3, which neither significantly under nor overfits eitherdataset in our reductions. The impact of f is explored in moredetailed in Section 5.1.2. For ANDROMEDA, the representa-tive choice of protection angle (minimum angular displacementof companion signal between two frames) is δ = . λ/ D , whichshows good performance for extreme-AO SPHERE data and wasthe final parameter choice used for the spectro-photometric anal-ysis of 51 Eri b (Samland et al. 2017). The contrast curve for δ = . λ/ D is included for demonstrating the impact of the protec-tion angle on contrast performance at small angular separations.Contrast curves for ANDROMEDA using δ = . , . , . λ/ D for both datasets are shown in Appendix B for completeness. The detection maps for both TRAP and ANDROMEDA areshown in Fig. 6. The first thing we note is that the remainingstructures in the detection map using TRAP are smoother andless correlated on spatial scales of 1 λ/ D . Besides 51 Eri b, wedo not detect any other signal with > σ significance. Table 2shows a summary of the obtained photometry for 51 Eri b usingTRAP and ANDROMEDA.Figure 7 shows the contrast curve for both reductions (leftpanel) and the factor gained in contrast by using TRAP com-pared to ANDROMEDA (right panel). The ANDROMEDA re-sults have been obtained using two di ff erent values for the pro-tection angle δ = . λ/ D and δ = . λ/ D . Because both algo-rithms obtain a 2D detection map from the forward model grid ofpositions, not only can we determine the median detection limitat a given separation, but also the azimuthal distribution of detec-tion contrasts. Thus, in addition to the median, we plot the 68%range as shaded regions reflecting the variability of the detectioncontrast along the azimuth. This is another important figure ofmerit to evaluate the performance of the algorithms. Article number, page 8 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations
Table 1.
Observing logTarget UT date Satellite spots a IRDIS Filter IRDIS DIT b T exp Field Rot. c Sr d DIMM seeing τ (sec, ◦ ) (ms)51 Eri 2015-09-25 no K12 16 ×
256 68 42 0.80–0.90 0.75 5.7 β Pictoris 2015-11-30 yes H23 4 ×
200 53 40 N / A 1.06 10.9
Notes.
Observational data used for the analyses. Seeing and coherence time ( τ ) correspond to the mean of the coronagraphic observation sequence. ( a ) Continuous satellite spot mode: calibration spots are always present in the coronagraphic data. ( b ) Detector integration time. ( c ) All observationwere centered on the meridian passage of the target with an airmass between 1.08 and 1.13. ( d ) Strehl ratio computed by AO system’s Real TimeComputer, and scaled to a wavelength of 1.6 µ m. − − − − O ff s e t ( a r cs e c )
51 Eridani (linear) − − − − O ff s e t ( a r cs e c )
51 Eridani (log) A DU A DU Fig. 4.
Contour map of median combined image of pre-processed 51 Eridani data cube. The left panel shows data in linear scaled brightness bins,the right panel in logarithmic scaling. − − − − O ff s e t ( a r cs e c ) beta Pic (linear) − − − − O ff s e t ( a r cs e c ) beta Pic (log) A DU A DU Fig. 5.
Contour map of median combined image of pre-processed beta Pic data cube. The left panel shows data in linear scaled brightness bins,the right panel in logarithmic scaling.
The detection limit obtained with TRAP is consistently lowerthan ANDROMEDA for separations up to about 10 λ/ D , atwhich point the results between TRAP and ANDROMEDA con-verge to the same background value. Especially at small separa-tion we clearly see a gain in contrast due to the absence of aprotection angle in our reduction. The ANDROMEDA curvescut o ff at about 1 λ/ D and 2 λ/ D respectively because no refer- ence data that fulfills the exclusion criterion exists (see Fig. A.1, ∼ ◦ elevation). The relative gain in contrast, with respect tothe spatial model (ANDROMEDA) is significant, with ramifica-tions for the detectability of close-in planets. It can be as high asa factor of six at 2 λ/ D for a protection angle of 1 λ/ D and fourat a separation of 1 λ/ D for a protection angle of 0 . λ/ D . Article number, page 9 of 21 & A proofs: manuscript no. trap − − O ff s e t ( a r cs e c ) TRAP − − O ff s e t ( a r cs e c ) ANDROMEDA -6.00-4.50-3.00-1.500.001.503.004.506.00 S NR -6.00-4.50-3.00-1.500.001.503.004.506.00 S NR Fig. 6.
Contour map of normalized detection maps obtained with TRAP ( f = .
3) and ANDROMEDA ( δ = . λ/ D ) for 51 Eri b. These mapsmust not be confused with derotated and stacked image. They represent the forward model result for a given relative planet position on the sky( ∆ RA, ∆ DEC), i.e. the conditional flux of a point-source predicted by the forward model given a relative position, corresponding to a trajectoryover the detector (all pixels a ff ected during the observation sequence). c o n t r a s t / D / D / D / D / D IWA
TRAP ( f = 0.3)ANDROMEDA ( = 0.5 / D )ANDROMEDA ( = 1.0 / D ) m a g n i t u d e F a c t o r g a i n e d i n c o n t r a s t / D / D / D / D / D IWA
ANDROMEDA ( = 0.5 / D ) / TRAP ( f = 0.3)ANDROMEDA ( = 1.0 / D ) / TRAP ( f = 0.3) Fig. 7.
Comparison between the contrast obtained with TRAP and two ANDROMEDA reductions for 51 Eri using the same input data for the K1band. TRAP has been run with 30% of available principal components, whereas the two ANDROMEDA reductions correspond to a protectionangle of δ = . λ/ D and δ = . λ/ D . Separations below the inner-working angle of the coronagraph are shaded and should only be interpretedrelative to each other, not in terms of absolute contrast, because the impact of coronagraphic signal transmission is not included in the forwardmodel of either pipeline. (left) The shaded areas around the lines correspond to the 16%–84% percentile intervals of contrast values at a givenseparation. (right) Factor in contrast gained using TRAP. The diminishing gain in contrast performance with separa-tion and subsequent agreement with ANDROMEDA at about10 λ/ D , is consistent with our expectation that the limiting factorof insu ffi cient field-of-view rotation for the spatial model ceasesto be important at larger separation. We note that in general theazimuthal variation of the contrast at each specific separationbin is consistently and significantly smaller for our model. Thisreflects the visual impression of an overall smoother detection map. This di ff erence in azimuthal variation is especially notice-able at separations >
350 mas, where the influence of the wind-driven halo declines (Fig. 4).Detailed tests with injected point sources were performed torule out that TRAP biases the retrieved signal-to-noise to anysignificant degree, i.e. any bias is small compared to the derivedstatistical uncertainties. The description of these tests and theirresults are shown in Appendix C.
Article number, page 10 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations
To confirm the reliability of our results we perform the samedata reduction with di ff erent complexities of the systematicsmodel. We test the impact of the fraction f of the maximumnumber of components N max = N frames used (see Eq. 13). Fig-ure 8 shows the contrast curve for the same data as used inFig. 7, with f = . , . , . , .
7. Firstly, we note the absenceof a trend towards better detection limits with increasing num-ber of principal components. In fact the contrast gets worse forlarge fractions. While a simultaneous fit of models of the planetsignal and systematics counteracts over-fitting, increasing thecomplexity of the systematics model beyond a certain point, theplanet model component becomes less constraining, resulting ina larger scatter of the planet contrast prediction. Secondly, thevalue of f = . ffi ciently complex ( f = .
1) are perform-ing slightly worse at close separation, which could be related tothe presence of the strong wind-driven halo and e ff ects from the(mis)alignment of the coronagraph with the star in addition tothe quasi-static speckle pattern.In spatial models one may have to choose a di ff erent modelcomplexity based on the separation, trying to compensate self-subtraction e ff ects from using smaller protection angles by usinga less complex model, but also having fewer spatial modes to re-construct. However, for our temporal approach the existence ofa model complexity that fits well globally is in agreement withour expectation. The number of frames available for trainingdoes not depend on separation because we do not have a tem-poral exclusion criterion. Additionally, a speckle lifetime analy-sis performed by Milli et al. (2016) for SPHERE did not showa strong separation dependence of speckle correlations linearlydecreasing correlations on timescales of tens of minutes. Thefast evolving and exponentially decaying correlations ( τ ∼ . ff ect too small to be expectedto be visible. For a temporal systematics model it is highly probable that theperformance of the algorithm scales with temporal sampling.In Fig. 9 we repeat our reduction with the same dataset binneddown by a factor of four (16 s to 64 s exposures). We again plotthe contrast of TRAP reductions with varying principal compo-nent fractions f and ANDROMEDA with the standard setting of δ = . λ/ D (left panel) for the same binned data. We also plotthe factor gained in contrast compared to ANDROMEDA (rightpanel). We can still see an improvement at small separations, butthe advantage of using a temporal model drops o ff faster, and itmay even perform worse at large separations. From ∼ λ/ D wedo not see any improvement and the improvement at small sepa-rations is smaller. This is consistent with our expectation of tem-poral models being able to capture more systematic variations ata finer time sampling. If the time sampling is poor the causal re-lationships in the systematics get averaged out and become moredi ffi cult to model.Figure 10 shows the advantage in contrast when using the un-binned (16 s exposure) compared to the binned (64 s exposure)data. We see that temporally binning the data has an adverse ef-fect on the contrast over almost the whole range of separations. c o n t r a s t / D / D / D / D / D IWA
TRAP ( f = 0.1)TRAP ( f = 0.3)TRAP ( f = 0.5)TRAP ( f = 0.7) m a g n i t u d e Fig. 8.
Contrast curves obtained with TRAP for 51 Eri when di ff erentfractions of the maximum number of principal components are used.Figure description analogues to Fig. 7. Table 2.
Photometry and SNR for unbinned 51 Eri b data
Method Contrast SNR(10 − )TRAP 6 . ± . . ± . λ/ D , we see an average40% improvement in contrast by using the faster temporal sam-pling. We note that both of the sampling rates shown here are toocoarse to model the short-lived speckle regime, and we expectfurther improvement by reducing the exposure time by anotherfactor of four or more ( ≤ β Pic: continuous satellite spot data with shortintegrations
A comparison between the detection maps obtained with TRAPand ANDROMEDA is shown in Fig. 11. For both reductions weused the same pre-reduced and aligned input data cubes, and thereduction parameters were the same as for the above 51 Eri data( f = . δ = . λ/ D for ANDROMEDA). Wefocus our discussion on one of the two channels (H2). The colorscaling di ff ers to account for the di ff erence in SNR of the detec-tion. With an SNR of 40, the TRAP detection is about four timeshigher than in ANDROMEDA. The SNR in H3 (not shown here)is even slightly higher, because the contrast to the host star ismore favorable at these wavelengths. Table 3 shows a summaryof the photometry results obtained for all reductions of β Pic bdiscussed in this section.Figure 12 shows the obtained contrast analogous to Fig. 7.Qualitatively, we see the same e ff ects as for 51 Eri, an increas-ingly significant contrast improvement at small separations fromusing a temporal model. We see an even more pronounced re-duction in azimuthal variation in contrast, i.e. the “width” of the Article number, page 11 of 21 & A proofs: manuscript no. trap c o n t r a s t / D / D / D / D / D IWA
TRAP ( f = 0.1)TRAP ( f = 0.3)TRAP ( f = 0.5)ANDROMEDA ( = 0.5 / D ) m a g n i t u d e F a c t o r g a i n e d i n c o n t r a s t / D / D / D / D / D IWA
ANDROMEDA ( = 0.5 / D ) / TRAP ( f = 0.3) Fig. 9.
Same as Figure 7, but each four frames of the 51 Eri b coronagraphic data are temporally binned to achieve a DIT of 64s. Figure descriptionanalogues to Fig. 7. F a c t o r g a i n e d i n c o n t r a s t / D / D / D / D / D IWA
TRAP 64s ( f = 0.3) / TRAP 16s ( f = 0.3) Fig. 10. E ff ect of using faster temporal sampling DIT 16s (unbinned)compared to slower sampling (longer averaging) with DIT 64s (binned)on contrast obtained with TRAP. Figure description analogues to Fig. 7. contrast curve. For this data, we additionally see a baseline in-crease in performance between 50 – 200% at larger separations( > λ/ D ) that we attribute to the better temporal sampling. Toconfirm this hypothesis, we have binned down the data by a fac-tor of 16 to obtain 1 minute exposures. The detection map ob-tained with TRAP is shown in Fig. 13 and the contrast com-pared to the unbinned ANDROMEDA reduction is shown inFig. 14. The planet signal is detected with a SNR of about 13, only slightly better than the SNR obtained with ANDROMEDAon the unbinned data. The obtained detection contrast curve isalso comparable to ANDROMEDA at this separation. We there-fore attribute the SNR improvement by a factor of four to thefact that our algorithm is capable of taking into account the bulkof short-timescale variations. We note that at the separation of β Pic b for the epoch of the observation ( ∼ λ/ D ), the exclusiontime for spatial models is on the order of minutes even for verysmall protection angles ( δ = . λ/ D ) and above 10 minutes forthe standard setting of δ = . λ/ D (see Fig. A.1), much longerthan the exposure time. As we use a non-local, temporal systematics model and do notattempt to reconstruct a spatial model for how the speckles“look”, we can forego aligning the data and run the algorithmon minimally pre-reduced (background and flat corrected) un-aligned data. To demonstrate this, we measure the center basedon the satellite spots for each frame, and use this varying cen-ter position to construct the lightcurve model for the planet, i.e.we do not shift the frames, instead we shift the forward modelfor our planet. We also apply the anamorphism correction forSPHERE to the relative position of the planet by reducing therelative separation of the model by 0.6% in y-direction (Maireet al. 2016) for each frame, instead of stretching the images. Wealso modulate the contrast of the planet model by the satellitespot amplitude variation measured for each frame. The result isshown in Fig. 15, with the aligned data on the left and the un-aligned data on the right. We note that the SNR of the detectionis virtually the same. There are only slight di ff erence in the resid- The background subtraction is needed for the unsaturated PSFmodel. For the coronagraphic sequence this may not be true, becausewe already include a free constant o ff set term in the design matrix (seeEq. 12 and 13). Furthermore, variations of the thermal background onlarge scales may be incorporated in the temporal systematics model.Article number, page 12 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations − − − − O ff s e t ( a r cs e c ) TRAP − − − − O ff s e t ( a r cs e c ) ANDROMEDA -32.80-24.80-16.80-8.80-0.807.2015.2023.2031.20 S NR -9.00-6.00-3.000.003.006.009.00 S NR Fig. 11.
Contour map of normalized detection maps obtained with TRAP ( f = .
3) and ANDROMEDA ( δ = . λ/ D ) for β Pic. These maps mustnot be confused with derotated and stacked image. They represent the forward model result for a given relative planet position on the sky ( ∆ RA, ∆ DEC), i.e. the conditional flux of a point-source predicted by the forward model given a relative position, corresponding to a trajectory over thedetector (all pixels a ff ected during the observation sequence). c o n t r a s t / D / D / D / D / D IWA
TRAP centered ( f = 0.1)TRAP centered ( f = 0.3)TRAP centered ( f = 0.5)ANDROMEDA ( = 0.5 / D ) m a g n i t u d e F a c t o r g a i n e d i n c o n t r a s t / D / D / D / D / D IWA
ANDROMEDA ( = 0.5 / D ) / TRAP ( f = 0.3) Fig. 12.
Comparison between the contrast obtained with TRAP and ANDROMEDA reductions for β Pic using the same input data for the H2band. TRAP has been run with 10%, 30%, and 50% of available principal components, whereas the ANDROMEDA reduction correspond to aprotection angle of δ = . λ/ D . Separations below the inner-working angle of the coronagraph are shaded and should only be interpreted relativeto each other, not in terms of absolute contrast, because the impact of coronagraphic signal transmission is not included in the forward model ofeither pipeline. (left) The shaded areas around the lines correspond to the 14%–84% percentile interval of contrast values at a given separation.(right) Contrast gained by using TRAP. uals structures. We do see a blob above the position of β Pic bthat edges above 5 σ in the reduction on unaligned data. It isdi ffi cult to evaluate the veracity of these structures due to thepresence of disk structure in the β Pic system.The contrast curves for: 1) aligned data without planetbrightness modulation, 2) including planet brightness modula-tion, and 3) unaligned data with planet brightness modulationare shown in Fig. 16. Taking into account the brightness modu- lation derived from the satellite spots does not have a noticeableimpact on the contrast limits. It does have a minimal impact onthe derived contrast of β Pic b, as it reduces the flux calibrationbias incurred by assuming an average or median contrast for theplanet flux in the planet model, instead of the real distribution.In the case of this observation the scatter of brightness variationis roughly Gaussian with a variability of ∼
6% centered on themean of the satellite spot brightness, without a large systematic
Article number, page 13 of 21 & A proofs: manuscript no. trap − − − − O ff s e t ( a r cs e c ) TRAP -14.0-10.5-7.0-3.50.03.57.010.514.0 S NR Fig. 13.
Contour map of normalized detection maps obtained withTRAP ( f = .
3) binned data of beta Pic (16x binning, 64s exposures).Figure description analogues to Fig. 7. F a c t o r g a i n e d i n c o n t r a s t / D / D / D / D / D IWA
ANDROMEDA ( = 0.5 / D ) / binned TRAP ( f = 0.3)binned ANDROMEDA ( = 0.5 / D ) / binned TRAP ( f = 0.3) Fig. 14.
Contrast ratio between TRAP ( f = .
3) reduction on tempo-rally binned data (16x binning, 64s exposures) and ANDROMEDA re-duction of the same binned data, as well as unbinned data, of β Pic.Figure description analogues to Fig. 7. trend. Taking into account this variation becomes more impor-tant the more unstable the conditions are and when an overalltrend is present, e.g. due to clouds reducing atmospheric trans-mission during part of the observation.A bigger di ff erence can be seen in the reduction of the un-aligned data. The contrast appears to be worse in the innermostregion covered by the coronagraphic mask, but slightly betteroutside of 3 λ/ D . The SNR of β Pic b, again, is virtually the sameas in the aligned case, but astrometry and photometry changeslightly.
Table 3.
Photometry and SNR for β Pic b
Method Mod. Align Bin Contrast SNR16x (10 − )TRAP no yes no 1 . ± .
04 38TRAP yes yes no 1 . ± .
04 38TRAP yes no no 1 . ± .
03 38TRAP no yes yes 1 . ± .
09 15ANDROMEDA no yes no 1 . ± .
13 9ANDROMEDA no yes yes 0 . ± .
12 8
Notes.
Overview of reduction results depending on whether satellitespot amplitude modulation, pre-aligned (centered) data, and temporalbinning is performed. All TRAP reductions used f = . δ = . λ/ D . The computational time needed to reduce the 51 Eri dataset (256frames) at one wavelength up to a separation of ∼
45 pixel ( ∼ f = .
3, PSF stamp size 21x21pixel) including the variance on the data is about 120 minuteson a single core on a laptop (intel CORE i7 vPro, 8th Gen; 16GB memory). The computational time can be roughly halved byusing a PSF stamp half the size for determining the reductionarea P Y (excluding the first Airy ring). This only has a minorimpact on the overall performance of the algorithm. Reducing f similarly reduces the computational time, which is important forlarge data sets with thousands of frames. The algorithm is paral-lelized on the level of fitting the model contrast for a given posi-tion, such that the grid of positions to explore is divided amongthe available cores. At the current version of the code, using fourcores reduces the time needed to about 40 minutes, i.e. by a fac-tor of three, due to ine ffi ciencies in memory sharing. We expecta close to linear relation with the number of cores after improve-ments to the code’s parallelization architecture. It is noteworthy that our algorithm’s computational speed scalesbetter with the number of frames in the observation sequencethan traditional spatial-based approaches that include a protec-tion angle. The absence of a temporal exclusion criterion meansthat the principal component decomposition has to be performedonly once for one assumed companion position, instead of hav-ing a separate training set for each frame in the sequence. Theonly increase in computational time stems from the need of de-composing a larger matrix once per model position and subse-quently inverting a larger system of linear equations for eachpixel. Testing the computational time for di ff erent temporal bin-ning factors, we note that the computational time scales betweenlinear and quadratic with the number of frames with a power-lowindex of about t ∝ N . , i.e. doubling the number of framesmore than doubles the computational time. We select a new set of reference pixels depending on the testedcompanion position, because we have to exclude the reductionarea from the training set. As such, the time spent on construct-ing the model is linearly proportional to the number of positionstested, which, when exploring a linear parameter space in ∆ RAand ∆ DEC means that the number of PCAs performed is pro-portional to the search area, such that N PCA ∝ r , where r is the Article number, page 14 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations − − − − O ff s e t ( a r cs e c ) aligned − − − − O ff s e t ( a r cs e c ) unaligned -7.00-5.57-4.14-2.71-1.290.141.573.004.435.86 S NR -7.00-5.57-4.14-2.71-1.290.141.573.004.435.86 S NR Fig. 15.
Contour map of normalized detection maps obtained with TRAP ( f = .
3) on aligned and unaligned data for beta Pic. These maps mustnot be confused with derotated and stacked image. They represent the forward model result for a given relative planet position on the sky ( ∆ RA, ∆ DEC), i.e. the conditional flux of a point-source predicted by the forward model given a relative position, corresponding to a trajectory over thedetector (all pixels a ff ected during the observation sequence). c o n t r a s t / D / D / D / D / D IWA
TRAP centered ( f = 0.3)TRAP centered w. modulation ( f = 0.3)TRAP uncentered w. modulation ( f = 0.3) m a g n i t u d e F a c t o r g a i n e d i n c o n t r a s t / D / D / D / D / D IWA baseline reduction / aligned wo ampl. modulationbaseline reduction / unaligned with ampl. modulation
Fig. 16. (left) Comparison between the contrast for β Pic obtained with TRAP on 1) aligned data not taking into account brightness modulation,2) taking it into account, and 3) on unaligned data with brightness modulation. TRAP has been run with 30% of available principal components.The shaded areas around the lines correspond to the 14%-84% percentile interval of contrast values at a given separation. (right) Contrast gaincompared to aligned data without including amplitude variations (baseline reduction). Figure description analogues to Fig. 7. separation from the central star. At the same time the number ofpixels a ff ected by a potential companion also increases with sep-aration N pix, a ff ected ∝ r . In terms of computational time, however,the time spend on PCA is relatively minor (once per tested po-sition), and the scaling with number of a ff ected pixels that needto be fit outweighs. Testing the algorithm with increasing outer-working angle (OWA), we derive a power-law index of about t ∝ OWA . . If computation time is an issue, one can easily useour algorithm exclusively for the inner-most region and combinethe results with an algorithm that scales better with separation further out, as we do not expect big performance improvementsat large separations compared to other algorithms.
6. Discussion
We have seen that a purely temporal model can be an alterna-tive to a purely spatial model, but it is clear that neither one is acomplete solution. There is a large unexplored space of spatio-temporal mixture models that would simultaneously take intoaccount temporal and spatial correlations. One way of imple-
Article number, page 15 of 21 & A proofs: manuscript no. trap menting such spatio-temporal hybrid models can be thought ofas extending the data vectors from being either time series ofindividual pixels or images at specific times, to time series ofpatches of pixels. In such a model, the time series of one patchcan be re-constructed in a basis set of vectors each containingthe time series of another patch of equal size, that are takenfrom a di ff erent part of the image (non-local model). Such anapproach would optimize both spatial and temporal similaritybetween multiple such patches. It may also be possible to im-prove the results of TRAP by fitting a spatial LOCI-like system-atics model to the image residuals after subtracting the temporalsystematics model or vise-versa. Spatial and temporal regressionapproaches are not mutually exclusive. They can be synergistic,because they optimize their models based on di ff erent correla-tions and di ff erent training data. This is in contrast to applying aspatial model iteratively (Brandt et al. 2013). The topic of disks is one aspect of non-local models and regres-sor selection that has not been discussed in this work and willrequire future research. Protection from self-subtraction by us-ing non-local training data, and protection from over-subtractionby simultaneously fitting a forward model could be a valuableproperty for disk imaging, where preserving the morphology ofthe object is paramount. It should be pointed out that similarto the spatial approaches (with the exception of RDI), a com-pletely azimuthally homogeneous structure will not be picked upby our algorithm in its current form, because we include a con-stant o ff set term in our fit. The current detection map based onpoint-source forward modeling is not suited for disk imaging andwould have to be adapted. It will pick up on (non-homogeneous)disk structures in the detection map, but should not be used di-rectly to study disk morphology. This work demonstrates the potential of the non-local, temporalsystematics modeling approach. However, there are still manypossible avenues for future improvements. The next step will beextending the algorithm to take into account spectral informa-tion. This can easily be achieved by adding the lightcurves ofpixels at other wavelengths to the reference set similar to whatis done in spatial models (e.g., Sparks & Ford 2002; Mesa et al.2015; Ru ffi o et al. 2017). Currently, the area of the detector af-fected by the signal of interest is excluded from the training data.However, because the planet position stays fixed, whereas thespeckle pattern scales with wavelength, we can shift the reduc-tion area inward / outward proportional to the wavelength and addthose signal-free pixels to the training set. This adds a local au-toregressive component to the model that traces the temporal be-havior of the speckles at the position of interest. Again, this canbe achieved without interpolating the raw data or prescribinga detailed chromatic behavior other than the rough wavelengthscaling needed for the reference selection.Another important step will be improving the fidelity ofthe forward model, e.g. directly including the coronagraphicthroughput model (e.g., Guerri et al. 2011; Mawet et al. 2013)in the forward model of the signal. This is a more consistent ap-proach than post-hoc adjusting the contrast curve by the corona-graphic transmission. As we push towards smaller inner-workingangles, a good understanding of the coronagraph at all wave-lengths will become important and should be modeled as well as measured on-sky as part of publicly available instrument cal-ibrations.There are a multitude of e ff ects on the companion signal that,in principle, can be included in the forward model, sch as distor-tion of the PSF shape, e.g., variations of the Strehl ratio, low-wind e ff ect (Cantalloube et al. 2018, 2020), smearing causedby integration time (Lafrenière et al. 2007), and optical aber-rations as, for example, measured by a focal plane wavefrontsensor (Wilby et al. 2017). Likewise, currently our algorithm im-plicitly assumes that all necessary information on the systematictemporal trends is encoded in other pixels. It has been shownfor transit photometry that “missing” information on systematictrends can be accounted for using auxiliary data or a Gaussianprocess trained on auxiliary data (Gibson 2014). It is possiblethat including additional external information (e.g., on the wind,state of the AO, position of derotator, temperature, focal-planewavefront sensing), could further improve the algorithmic per-formance significantly.The temporal modeling approach may also prove beneficialfor instruments with large pixel scales, or undersampled PSFs,as is sometimes the case for IFUs. Large pixels / spaxels canexacerbate the problems caused by insu ffi cient field-of-viewrotation. In general the algorithm introduced in this work shouldbe more robust against problems associated with small field-of-view rotation, because we do not use a temporal exclusioncriterion. An interesting use-case, that has not been exploredin this work, is the application to space-based observatories.For observations taken at di ff erent roll angles we can alsobuild a temporal forward model that would take the form of astep-function for a ff ected pixels. This, again, would allow us totake into account systematics in real-time, e.g. from instrumentjitter, because we do not need to exclude any frames from thetraining data.The algorithm as introduced in this work is optimized forcompanion searches over a grid of possible positions. If our goalis the detailed characterization of a planet of known position –or even a disk – this approach may not be optimal. Future devel-opment will include exploration of specific models in a detailedcharacterization step using sampling approaches such as MCMCor nested sampling to obtain photometric and astrometric infor-mation as well as their covariance simultaneously (e.g., Wanget al. 2016b). This would allow, both, the exploration of morecomplex physical models such as forward models of debris disks(Olofsson et al. 2016), as well as more sophisticated systemat-ics models. Using nested sampling, combined with a negativemodel injection approach, would provide an easy way to includeall dimensions of the data in the forward model (e.g., includingthe spectrum of the planet Ru ffi o et al. 2017) and perform directmodel comparison based on Bayesian evidence as often done intransit and radial velocity detections (e.g., Espinoza et al. 2019).
7. Summary and conclusions
We have presented a new paradigm of using a temporal, non-local systematics models to more e ff ectively search for pointsources at the scientifically most interesting small separationswhere traditional ADI algorithms have issues by design. Thisnew method allows us to address a persistent problem in highcontrast imaging: the lack of good and uncontaminated trainingdata at small separations. By building a non-local, causal modelwe show that we can circumvent the problem of self-subtractionentirely, while still presenting a powerful model for the system-atics that are limiting our ability to detect planets. By its nature, Article number, page 16 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations the time-domain model is sensitive to instantaneous changes inobserving conditions, which is not the case in spatial approachesthat need a reference library of images sharing the same specklecharacteristics. Furthermore, the temporal model is not contin-gent on the data being perfectly aligned: as such interpolationsand re-sampling of the data in space (and in future addition tothe method, wavelength) can be avoided entirely.Our implementation, called TRAP, is open-source and pub-licly available.We have shown on two datasets that the TRAP pipeline per-forms as well or better than a similar spatial approach with astrong improvement in contrast by a factor between 1.5 and 6at angular separations < λ/ D . Beyond this separation, the im-provement strongly depends on the temporal sampling of the ob-servation sequence. The azimuthal variance of the achieved con-trast across all separations is strongly reduced using the temporalsystematics model compared to a spatial systematics model. In-creasing the integration time from 16 s to 64 s for the 51 Eri bdataset led to a decrease in average contrast gained by ∼ β Pic), we can achieve a sig-nificant overall improvement of the contrast by a factor of two,even at separation between 3 – 10 λ/ D . The SNR measured for β Pic b significantly increased from about 10 with the spatialmodel to about 40 by making full use of the systematics informa-tion present in the data on the short time scales. We conclude thatthe e ff ect of exposure time on the achievable contrast is under-explored in the literature. Spatial algorithms currently employedare not able to make optimal use of the information containedin short time-scale variations, because typical exclusion time-scales are significantly larger than the exposure times.Our results have shown that fitting the planet and system-atics model simultaneously constitutes a self-regulating processon the achieved contrast when we increase model complexity:increasing the systematics model complexity (i.e. the number ofprinciple components used) does not automatically lead to “bet-ter” contrasts, highlighting the benefit of a combined model fit.We have demonstrated, that our temporal approach can beapplied to minimally pre-reduced data without aligning theframes. This is achieved by adjusting the forward model positionthat generates the companion lightcurve according to the star’scenter position, the anamorphism, and excluding all bad pix-els from the training and reduction sets instead of interpolatingthem. This reduction achieves very similar results to elaboratelypre-reduced data and reduces the need for intrusive data manip-ulation steps (interpolation, resampling). This property can beparticularly useful when taking into account data uncertainties.Another benefit is a strong reduction in processing time neededfor alignment and bad pixel interpolation, which can take sig-nificant resources for datasets with many exposures or work-ing on entire surveys. The ability to post-process data withoutre-sampling could prove beneficial for SPHERE-IFS in the fu-ture, because the instrument uses a hexagonal lenslet geometry.The output images have to be re-sampled to a rectilinear gridfor traditional post-processing pipelines. With TRAP we couldperform the analysis on the native image geometry.Like ANDROMEDA, our algorithm does not require derota-tion of frames. Spatial filtering, which improves ANDROMEDAperformance and is also used in pyKLIP, is not needed for ouralgorithm to perform well.We do not recommend dithering for pupil-tracking data asit is not necessary and can interfere with the performance oftemporal models. We further recommend the use of continuoussatellite spot mode to improve the forward model performance with accurate center and amplitude variations. We strongly rec-ommended exploring shorter integration times for observationsequences. Decreasing the integration time of IRDIS from 64 sto 4 s increases the observation overheads by about 15 percentpoints. If our scientific interest is companions at small separa-tions, we are in the speckle limited regime and the increase inread-noise is likely negligible, easily balanced by the many-foldincrease in algorithmic performance.Lastly, future and current development and deployment ofcoronagraphs with smaller inner-working angles will further in-crease the importance of this class of algorithm. Acknowledgements.
This work has made use of the SPHERE Data Centre,jointly operated by OSUG / IPAG (Grenoble), PYTHEAS / LAM / CeSAM (Mar-seille), OCA / Lagrange (Nice) and Observatoire de Paris / LESIA (Paris). Wethank P. Delorme and E. Lagadec (SPHERE Data Centre) for their e ffi cienthelp during the data reduction process. SPHERE is an instrument designed andbuilt by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidel-berg, Germany), LAM (Marseille, France), LESIA (Paris, France), LaboratoireLagrange (Nice, France), INAF–Osservatorio di Padova (Italy), Observatoireastronomique de l’Université de Genève (Switzerland), ETH Zurich (Switzer-land), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) incollaboration with ESO. SPHERE was funded by ESO, with additional contri-butions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzer-land) and NOVA (Netherlands). SPHERE also received funding from the Eu-ropean Commission Sixth and Seventh Framework Programmes as part of theOptical Infrared Coordination Network for Astronomy (OPTICON) under grantnumber RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 forFP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). T.H. ac-knowledges support from the European Research Council under the Horizon2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. Wethank F. Cantalloube, and B. Pope for their comments. J.B. acknowledges sup-port from the European Research Council under the European Union’s Horizon2020 research and innovation program ExoplANETS-A under grant agreementNo. 776403. References
Amara, A. & Quanz, S. P. 2012, MNRAS, 427, 948Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155Bohn, A. J., Kenworthy, M. A., Ginski, C., et al. 2020, MNRAS, 492, 431Bonse, M. J., Quanz, S. P., & Amara, A. 2018, ArXiv e-prints[ arXiv:1804.05063 ]Brandt, T. D., Kuzuhara, M., McElwain, M. W., et al. 2014, ApJ, 786, 1Brandt, T. D., McElwain, M. W., Turner, E. L., et al. 2013, ApJ, 764, 183Cantalloube, F., Farley, O. J. D., Milli, J., et al. 2020, A&A, 638, A98Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, A&A, 582, A89Cantalloube, F., Por, E. H., Dohlen, K., et al. 2018, A&A, 620, L10Casertano, S., Lattanzi, M. G., Sozzetti, A., et al. 2008, A&A, 482, 699Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, in SF2A-2017: Proceed-ings of the Annual meeting of the French Society of Astronomy and Astro-physics, DiCrepp, J. R., Johnson, J. A., Howard, A. W., et al. 2012, ApJ, 761, 39Delorme, P., Meunier, N., Albert, D., et al. 2017, in SF2A-2017: Proceedings ofthe Annual meeting of the French Society of Astronomy and Astrophysics,ed. C. Reylé, P. Di Matteo, F. Herpin, E. Lagadec, A. Lançon, Z. Meliani, &F. Royer, 347–361Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in Proc. SPIE, Vol. 7014,Ground-based and Airborne Instrumentation for Astronomy II, 70143LDou, J., Ren, D., Zhao, G., et al. 2015, ApJ, 802, 12Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262Fernandes, R. B., Mulders, G. D., Pascucci, I., Mordasini, C., & Emsenhuber, A.2019, ApJ, 874, 81Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, A&A, 618, A138Foreman-Mackey, D., Montet, B. T., Hogg, D. W., et al. 2015, ApJ, 806, 215Fusco, T., Sauvage, J.-F., Petit, C., et al. 2014, in Proc. SPIE, Vol. 9148, AdaptiveOptics Systems IV, 91481UGalicher, R., Marois, C., Macintosh, B., Barman, T., & Konopacky, Q. 2011,ApJ, 739, L41Gerard, B. L. & Marois, C. 2016, Society of Photo-Optical Instrumentation En-gineers (SPIE) Conference Series, Vol. 9909, Planet detection down to a few λ / D: an RSDI / TLOCI approach to PSF subtraction, 990958Gibson, N. P. 2014, MNRAS, 445, 3401Goebel, S. B., Guyon, O., Hall, D. N. B., et al. 2018, PASP, 130, 104502Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, AJ, 154, 7
Article number, page 17 of 21 & A proofs: manuscript no. trap
Grandjean, A., Lagrange, A. M., Beust, H., et al. 2019, A&A, 627, L9Gro ff , T. D., Kasdin, N. J., Limbach, M. A., et al. 2015, in Proc. SPIE, Vol. 9605,Techniques and Instrumentation for Detection of Exoplanets VII, 96051CGuerri, G., Daban, J.-B., Robbe-Dubois, S., et al. 2011, Experimental Astron-omy, 30, 59Kuhn, J. R., Potter, D., & Parise, B. 2001, ApJ, 553, L189Lafrenière, D., Marois, C., Doyon, R., & Barman, T. 2009, ApJ, 694, L148Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660,770Lagrange, A.-M., Boccaletti, A., Langlois, M., et al. 2019, A&A, 621, L8Lagrange, A.-M., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proceedings of the Na-tional Academy of Science, 111, 12661Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, in SPIE Conf. Ser., Vol.9908, 990834Marois, C., Correia, C., Galicher, R., et al. 2014, in Society of Photo-OpticalInstrumentation Engineers (SPIE) Conference Series, Vol. 9148, Society ofPhoto-Optical Instrumentation Engineers (SPIE) Conference Series, 0Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ,641, 556Marois, C., Lafrenière, D., Macintosh, B., & Doyon, R. 2008, ApJ, 673, 647Mawet, D., Pueyo, L., Carlotti, A., et al. 2013, ApJS, 209, 7Mesa, D., Gratton, R., Zurlo, A., et al. 2015, A&A, 576, A121Milli, J., Banas, T., Mouillet, D., et al. 2016, in Proc. SPIE, Vol. 9909, AdaptiveOptics Systems V, 99094ZMüller, M. & Weigelt, G. 1985, in Society of Photo-Optical InstrumentationEngineers (SPIE) Conference Series, Vol. 556, International Conference onSpeckle, ed. H. H. Arsenault, 270–273Nielsen, E., De Rosa, R., Macintosh, B., et al. 2019, in AAS / Division for Ex-treme Solar Systems Abstracts, Vol. 51, AAS / Division for Extreme Solar Sys-tems Abstracts, 100.02Olofsson, J., Samland, M., Avenhaus, H., et al. 2016, A&A, 591, A108Pairet, B., Cantalloube, F., Gomez Gonzalez, C. A., Absil, O., & Jacques, L.2019, MNRAS, 487, 2262Pavlov, A., Möller-Nilsson, O., Feldt, M., et al. 2008, in Proc. SPIE, Vol. 7019,Advanced Software and Control for Astronomy II, 701939Perrin, M. D., Sivaramakrishnan, A., Makidon, R. B., Oppenheimer, B. R., &Graham, J. R. 2003, ApJ, 596, 702Perryman, M., Hartman, J., Bakos, G. Á., & Lindegren, L. 2014, ApJ, 797, 14Pueyo, L., Crepp, J. R., Vasisht, G., et al. 2012, ApJS, 199, 6Racine, R., Walker, G. A. H., Nadeau, D., Doyon, R., & Marois, C. 1999, PASP,111, 587Ren, D., Dou, J., Zhang, X., & Zhu, Y. 2012, ApJ, 753, 99Rosenthal, E. D., Gurwell, M. A., & Ho, P. T. P. 1996, Nature, 384, 243Ruane, G., Mawet, D., Kastner, J., et al. 2017, AJ, 154, 73Ru ffi o, J.-B., Macintosh, B., Wang, J. J., et al. 2017, ApJ, 842, 14Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57Schölkopf, B., Hogg, D. W., Wang, D., et al. 2016, Proceedings of the NationalAcademy of Sciences, 113, 7391Smith, B. A. & Terrile, R. J. 1984, Science, 226, 1421Soummer, R. 2005, ApJ, 618, L161Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28Sparks, W. B. & Ford, H. C. 2002, ApJ, 578, 543Stolker, T., Bonse, M. J., Quanz, S. P., et al. 2019, A&A, 621, A59Vigan, A., Bonavita, M., Biller, B., et al. 2017, A&A, 603, A3Wahhaj, Z., Cieza, L. A., Mawet, D., et al. 2015, A&A, 581, A24Wang, D., Hogg, D. W., Foreman-Mackey, D., & Schölkopf, B. 2016a, PASP,128, 094503Wang, J. J., Graham, J. R., Pueyo, L., et al. 2016b, AJ, 152, 97Wilby, M. J., Keller, C. U., Snik, F., Korkiakoski, V., & Pietrow, A. G. M. 2017,A&A, 597, A112Xuan, W. J., Mawet, D., Ngo, H., et al. 2018, AJ, 156, 156 Article number, page 18 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations
Appendix A: Impact of temporal exclusion criteriaon traditional ADI
The temporal exclusion criterion (protection angle) has a strongimpact on the available data for training the systematics model intraditional spatial ADI approaches. Figure A.1 shows the impactof di ff erent assumed protection angles for a 90 minute observa-tion sequence centered around meridian passage computed fortargets at di ff erent maximum elevations above the horizon. Thecolors denotes the maximum elevation of the target (40 ◦ , 60 ◦ ,80 ◦ ), which result in di ff erent overall field-of-view rotation (25 ◦ ,40 ◦ , 90 ◦ ) and correspondingly slower or faster parallactic an-gle change rates. The left panel shows the average training datafraction excluded due to a the chosen protection angle (0 . λ/ D ,0 . λ/ D , 1 . λ/ D ) averaged over the observation sequence. Thehorizontal lines corresponding to the fraction of training data lostif instead a simple spatial exclusion criterion is used such as inthis work. In the case of a simple spatial exclusion criterion, asused in this work, which excludes all pixels at a given separa-tion a ff ected by a companion signal, and assumes only pixels ata comparable separation are used to build the systematics model,the excluded training data fraction is independent of separationand simply given by the fraction of total field-of-view rotationangle compared to 360 ◦ , i.e. about 10% for 40 ◦ rotation. In mostreal-world cases the fraction of available training data excludedby the spatial exclusion criterion is significantly less at small sep-arations compared to a temporal exclusion criterion (protectionangle).Analogously, the right panel of Fig. A.1 shows the averagedistance in time to the closest available frame in the trainingset that remains available under a given protection angle. Shownare observation sequences of targets at di ff erent maximum eleva-tions and assuming di ff erent protection angles. It is noteworthyto point out that even at larger separations this time di ff erence isstill significantly larger than typical integration times. Appendix B: Impact of protection angle onANDROMEDA reductions
The impact of the protection angle δ on the contrast limits ob-tained with ANDROMEDA is shown in Fig. B.1. In our studywe select δ = . λ/ D as being representative of the algorithm’sperformance, the same value as used in Samland et al. (2017).It produces consistent and reliable results for SPHERE data atall separations. Choosing smaller angles can marginally improveperformance at very small separations, but may su ff er from a po-tential increase in systematic bias. ANDROMEDA is based onforward modeling the expected signal in di ff erence images. Forsmall displacements most of the planet signal is subtracted andtherefore adds more noise than signal to the analysis for faintcompanions. This e ff ect is worse for data with short integra-tion times. Large protection angles negatively impact the perfor-mance at small separations due to excluding significant fractionsof the data, as discussed in length this work. All ANDROMEDAreductions use the standard spatial frequency filtering fraction of1 / Appendix C: Recovery of injected signals
A common way to test the fidelity of algorithms and post-processing pipelines is to inject a known signal into the data andattempt to retrieve it. In this work, because we build a forwardmodel of the expected set of lightcurves for each tested compan-ion position anyways, it is a simple additional step to inject a companion model at a desired contrast into the raw data beforefitting. Additionally, as we run the algorithm for each position in-dividually, we do not have to worry about contamination of thetraining set from multiple injected signals. The time necessaryfor performing this test is therefore only insignificantly sloweras running TRAP normally. The detailed procedure is describedbelow:1. Create a 2D contrast and detection map following the algo-rithm as described in this paper. This gives us the detectionlimit at each tested planet position.2. Run the algorithm again, but inject a signal for the tested po-sition with the desired significance as determined from thecontrast map. Obtain a map of retrieved contrasts for all po-sitions.3. Normalize the resulting uncertainty and SNR maps with thenormalization values as determined from the case withoutinjected signal.This results in maps of the retrieved signal and the associated un-certainty at each position in the detection image. We performedthis test on the IRDIS K2-band 51 Eri b dataset described inthis paper. The K2-band was chosen, because 51 Eri b is not de-tected at this wavelength, reducing potential e ff ects of the realsignal on the injected signals. The test was performed using asignal corresponding to 5 σ for each position from a separationof 4 pixels ( ∼ λ/ D ) out to a separation of 50 pixels (613 mas).The left panel of Figure C shows a normalized histogram of theretrieved SNR over all positions, together with a fitted Gaussiandistribution (mean µ = .
03 and standard deviation σ = . σ , showing thatthere is weak to no systematic bias.The right panel of Figure C shows the mean and standarddeviation of the distribution of retrieved SNR over separation in3-pixel wide annuli analogous to the annuli used to normalizethe SNR images. The mean of the detection significance staysclose to the expected 5 σ without significant systematic devia-tions, except for a small underestimation of the signal strengthat around 2 λ/ D , which corresponds to the edge of the corona-graph. Retrieved contrasts are on average the true contrasts andscatter within about 1 σ of the values, confirming the reliabilityof the obtained photometric values. The remaining small biasescan be understood and corrected for using the bias map, i.e. themap showing the di ff erence between the measured and injectedcontrast. The mean deviation from true per separation can beused as a bias correction of the reduction procedure and elimi-nate any remaining over or underfitting of the signal. This doesnot, of course, correct for biases that come from assuming an in-complete or wrong companion model in the forward model. Theleft panel of Fig. C shows the recovered minus true contrast cor-rected by the median deviation from the true injected values ofthe separation bin. The mean of this distribution is µ = . ff ected. The right panel of Fig. C shows thebias corrected detection significance over separation. Article number, page 19 of 21 & A proofs: manuscript no. trap D )0.00.20.40.60.81.0 A v e r a g e t r a i n i n g d a t a f r a c t i o n l o s t FoV-rot. 25°FoV-rot. 40°FoV-rot. 90° = 0.3 / D ; elev. 40°= 0.5 / D ; elev. 40°= 1.0 / D ; elev. 40°= 0.3 / D ; elev. 60°= 0.5 / D ; elev. 60°= 1.0 / D ; elev. 60°= 0.3 / D ; elev. 80°= 0.5 / D ; elev. 80°= 1.0 / D ; elev. 80° D )020406080 A v g . s e p . t o c l o s e s t t r a i n i n g d a t a ( m i n ) = 0.3 / D ; elev. 40°= 0.5 / D ; elev. 40°= 1.0 / D ; elev. 40°= 0.3 / D ; elev. 60°= 0.5 / D ; elev. 60°= 1.0 / D ; elev. 60°= 0.3 / D ; elev. 80°= 0.5 / D ; elev. 80°= 1.0 / D ; elev. 80° Fig. A.1.
The e ff ect of a di ff erent protection angle (0.3, 0.5, 1.0 λ/ D ) over separation from the central star. All values are computed for anobservation sequence of 90 minutes. The color of the line encodes di ff erent elevations of the target over the horizon at meridian passage. Theline style gives the protection angles. The left panel shows the total fraction of training data lost due to the protection angle averaged over thewhole sequence. The right panel shows the average temporal separation to the nearest viable frame outside the protection zone. It should be notedthat for a temporal systematics model (horizontal lines) which excludes all pixels a ff ected by signal, the exclusion is solely determined by thefield-of-view rotation. At 40 ◦ total rotation, we would therefore exclude ∼
10% of the data, regardless of separation, presenting a big advantage atshort separations over temporal exclusion criteria. c o n t r a s t / D / D / D / D / D IWA
TRAP ( f = 0.3)ANDROMEDA ( = 0.3 / D )ANDROMEDA ( = 0.5 / D )ANDROMEDA ( = 1.0 / D ) m a g n i t u d e
51 Eridani c o n t r a s t / D / D / D / D / D IWA
TRAP centered ( f = 0.3)ANDROMEDA ( = 0.3 / D )ANDROMEDA ( = 0.5 / D )ANDROMEDA ( = 1.0 / D ) m a g n i t u d e Pictoris
Fig. B.1.
Comparison between the contrast obtained with TRAP and three ANDROMEDA reductions for 51 Eri (left) and β Pictoris (right).TRAP has been run with 30% of available principal components, whereas the three ANDROMEDA reductions correspond to a protection angleof δ = . , . , . λ/ D . Separations below the inner-working angle of the coronagraph are shaded and should only be interpreted relative to eachother, not in terms of absolute contrast, because the impact of coronagraphic signal transmission is not included in the forward model of eitherpipeline. The shaded areas around the lines correspond to the 16%–84% percentile intervals of contrast values at a given separation.Article number, page 20 of 21. Samland et al.: TRAP: A temporal systematics model for improved direct detection of exoplanets at small angular separations D e t e c t i o n s i g n i f i c a n c e () Fig. C.1. (left) Histogram of retrieved SNR values for injected 5 σ signals in the reduced field-of-view, overplotted with a Gaussian fit. (right)Mean and standard deviation of detection significance of injected 5 σ signal over separation computed in 3-pixel wide annuli. D e t e c t i o n s i g n i f i c a n c e () Fig. C.2. (left) Histogram of deviation of retrieved contrast from true values after separation dependent bias correction for injected 5 σ signalsin the whole reduced field-of-view, overplotted with a Gaussian fit. (right) Mean and standard deviation of detection significance of injected 5 σσ