DETECTION AND CHARACTERIZATION OF EXOPLANETS AND DISKS USING PROJECTIONS ON KARHUNEN–LOÈVE EIGENIMAGES

An astronomical PSF can be seen as an N-dimensional stochastic process indexed over a scalar quantity ψ = ψ(t, λ, θ, ϕ, T_exp, m_Star), which is representative of the state of the instrument–telescope system (for instance, time, wavelength, field rotation angle, wavefront error, exposure time, and star magnitude) at the time of the exposure. With N pixels in the image, the PSF intensity can be written as

$$\begin{equation} I_{\psi }(n)=I(p[n],q[n],\psi), \end{equation} \tag{ 1 }$$

where we have rearranged the detector coordinates in one dimension using the pixel mapping functions p and q, and the pixel index n (n ∈ [1, n]).

We assume an observation sequence that yields a target image

$$\begin{equation} T(n) = I_{\psi _0}(n) + \epsilon \, A(n), \end{equation} \tag{ 2 }$$

which might ( = 1) or might not ( = 0) contain faint astronomical signal A(n), and a set of reference images

$$\begin{equation} \lbrace R_k(n) = I_{\psi _k}(n) \rbrace _{k = 1...K}. \end{equation} \tag{ 3 }$$

ψ₀(n) and ψ_k(n) represent the realizations of the state of the telescope-instrument system, respectively, for the target image and the set of references.

Using prior information about the observation strategy, these references have been selected so that they do not contain any astronomical signal in the search region $\mathcal {S}$ (dimension $N_{\mathcal {S}}$ pixels) of the PSF. Moreover, the observations have been carried out at neighboring states of the telescope/instrument system: |ψ_i − ψ_j| < δ for all i, j ∈ [0, K], where δ is a small quantity that can be adjusted in practice using the correlation between images as a proxy.

The goal of PSF subtraction is to reconstruct $I_{\psi _0}(n)$ using the set of reference images and to subtract it from the target image since $\epsilon A(n) = T - I_{\psi _0}(n)$. Unfortunately, I_ψ(n) is a continuous random process of ψ and thus an infinite amount of references would be needed to reconstruct $I_{\psi _0}(n)$ exactly. In practice, we can only seek to determine the best estimate $\hat{I}_{\psi _0}(n)$ of $I_{\psi _0}(n)$ using the limited information contained in the K references.

Assuming that each reference is unique, the ensemble $\lbrace I_{\psi _k} \rbrace _{k = 1...K}$ spans $\mathcal {I}_{K}$, a K-dimensional sub-space of $\mathbb {R}^{N_{\mathcal {S}}}$, which has an infinity of orthonormal bases {Z_k}_{k = 1...K}. This problem of estimating $\hat{I}_{\psi _0}(n)$ can be formulated as follows.

What is the basis set {Z_k(n)}_{k = 1...K} most likely to minimize the distance between a random realization of I_ψ(n) (ψ = ψ₀ in our case) and $\mathcal {I}_{K}$? More formally,

$$\begin{equation} \min _{\lbrace Z_k\rbrace } \left\lbrace E_{\psi }\left[ \sum _{n = 1}^{N_{\mathcal {S}}} \left(I_{\psi } - \sum _{k = 1}^{K} <I_{\psi },Z_k>_{\mathcal {S}} Z_k (n)\right)^2 \right] \right\rbrace, \end{equation} \tag{ 4 }$$

where E_ψ[ · ] stands for the expected value over the telescope realizations and $\langle \cdot,\cdot \rangle _{\mathcal {S}}$ denotes the inner product over $\mathcal {S}$. This is a well-known problem in signal processing and, assuming that all the images are of zero mean over $\mathcal {S}$, the optimal basis set {Z^KL_k}_{k = 1...K} is given by the Karhunen–Loève transform of $\lbrace I_{\psi _k} \rbrace _{k = 1...K}$ (Karhunen 1947; Loève 1948). Moreover, if one chooses to only describe $I_{\psi _0}$ over a K_klip < K dimensional sub-space, then the optimal basis set is given by the truncated Karhunen–Loève transform of the references, $\lbrace Z^{{\rm KL}}_k\rbrace _{k = 1...K_{{\rm klip}}}$

The KL image projection (KLIP) algorithm consists of the following steps.

• 1.  
  Partition the target T(n) and reference R_k(n) images in an ensemble of search areas $\mathcal {S}$, and subtract their average values so that they have zero mean over each $\mathcal {S}$.
• 2.  
  Compute the Karhunen–Loève transform of the set of reference PSFs {R_k(n)}_{k = 1...K}:

  $$\begin{equation} Z^{{\rm KL}}_k(n) = \frac{1}{\sqrt{\Lambda _k}}\sum _{p=1}^{K} c_k(\psi _p) R_p(n), \end{equation} \tag{ 5 }$$

  where the vectors C_k = [c_k(ψ₁)...c_k(ψ_K)] are the eigenvectors of the K-dimensional covariance matrix E_RR of the R_k(n) references over $\mathcal {S}$, and {Λ}_{k = 1...K} are its eigenvalues.

  $$\begin{eqnarray} E_{RR}[i,j] &=& \sum _{n = 1}^{N_{\mathcal {S}}} R_i(n)R_j(n) \end{eqnarray} \tag{ 6 }$$

  $$\begin{eqnarray} E_{RR}\cdot C_k &=& \Lambda _k C_k, \; \Lambda _1>\Lambda _2>....>\Lambda _N. \end{eqnarray} \tag{ 7 }$$
• 3.  
  Choose a number of modes K_klip to keep in the estimate $\hat{I}_{\psi _0}(n)$.
• 4.  
  Compute the best estimate $\hat{I}_{\psi _0}(n)$ of the actual PSF $I_{\psi _0}(n)$ from the projection of the science target T(n) on the {Z^KL_k(n)}_{k = 1...Kklip} Karhunen–Loève basis:

  $$\begin{equation} \hat{I}_{\psi _0}(n) =\sum _{k = 1}^{K_{{\rm klip}}} <T,Z^{{\rm KL}}_k>_{\mathcal {S}} Z^{{\rm KL}}_k(n). \end{equation} \tag{ 8 }$$
• 5.  
  Calculate the final image $F(n) = T(n) - \hat{I}_{\psi _0}(n)$,

  $$\begin{eqnarray} F(n) &=& \left(I_{\psi _0}(n) - \sum _{k = 1}^{K_{{\rm klip}}} <I_{\psi _0},Z^{{\rm KL}}_k>_{\mathcal {S}} Z^{{\rm KL}}_k (n)\right) \nonumber \\ &&+ \epsilon \left(A(n) - \sum _{k = 1}^{K_{{\rm klip}}} <A,Z^{{\rm KL}}_k>_{\mathcal {S}} Z^{{\rm KL}}_k (n)\right). \qquad \end{eqnarray} \tag{ 9 }$$

Equation (9) provides fundamental information pertaining to the optimal choice of K_klip for a given search zone $\mathcal {S}$. The variance in the search zone $\mathcal {S}$ after the PSF subtraction is given by

$$\begin{eqnarray} \sigma ^2_\mathcal {S} &=& ||I_{\psi _0}||^2 \left({1 - \sum _{k = 1}^{K_{{\rm klip}}}\left(\frac{ <I_{\psi _0},Z^{{\rm KL}}_k>_{\mathcal {S}}}{||I_{\psi _0}||} \right)^2}\right) \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \hspace{-3.7pc}& \simeq & ||I_{\psi _0}||^2 \left({1 - \frac{\sum _{k = 1}^{K_{{\rm klip}}} \Lambda _k}{\sum _{k = 1}^{K} \Lambda _k}}\right), \end{eqnarray} \tag{ 11 }$$

which is a decreasing function of K_klip. In Figure 1, we show the impact of K_klip on the residual variance integrated over entire KLIP-reduced images (i.e., $\mathcal {S}$ is the full image), using HST NICMOS coronagraphic data preselected to have no astrophysical signal. Figure 1 shows that the empirical residual variance after PSF subtraction follows the theoretical estimate in Equation (11).

Zoom In Zoom Out Reset image size

Similarly, Equation (9) provides a simple expression for the signal:

$$\begin{equation} S = \sum _{ A(n) \ne 0} \left(A(n) - \sum _{k = 1}^{K_{{\rm klip}}} <A,Z^{{\rm KL}}_k>_{\mathcal {S}} Z^{{\rm KL}}_k(n) \right). \end{equation} \tag{ 12 }$$

For the dominant modes (i.e., k ≪ K) $< A,Z^{{\rm KL}}_k>_{\mathcal {S}}\,\simeq 0$ since the morphology of the astronomical signal cannot be reproduced by a PSF realization. Indeed the PSF of a faint planet is localized over a few pixels and thus almost orthogonal to the main modes of the telescope's PSF realizations. Naturally, this argument depends on the spatial morphology of the signal: in Figure 2 we show that the algorithm throughput decreases when K_klip increases for several types of faint astrophysical object. Given priors on the shape of the astrophysical signal (i.e., planet position and flux, or a disk model), it is possible to estimate directly the signal-to-noise ratio (S/N) as a function K_klip when $\mathcal {S}$ is small enough so that Equation (11) represents the local noise in the vicinity of the source. It is therefore possible to choose K_klip a priori in order to optimize the S/N for a given putative source in the image.

Zoom In Zoom Out Reset image size

In practice the optimal search area $\mathcal {S}$ might be too large to use Equations (11) and (12) to pre-select K_klip. Instead K_klip can be tuned by inspecting the S/N of putative sources in reduced images. However, the second term of Equation (9) carries all the information necessary to characterize faint astrophysical signal once it has been detected. Indeed, the impact of the algorithm on the signal only depends on the projection of A(n) on the $\lbrace Z^{{\rm KL}}_k\rbrace _{k=1..K_{{\rm klip}}}$ and is not a function of $I_{\psi _0}$. As a consequence KLIP enables forward modeling of detected planets and disks by fitting directly an astrophysical model A_ξ(n) to the final reduced images F(n), where ξ represents the model parameters (e.g., astrometry and photometry for a planet, or geometry, grain distribution for a disk). ξ can be obtained by minimizing the following cost function:

$$\begin{equation} \min _{\lbrace \xi \rbrace } \left\lbrace \sum _{n = 1}^{N_{\mathcal {S}}} \left(F(n)-A_{\xi }(n) + \sum _{k = 1}^{K_{{\rm klip}}} <A_{\xi },Z^{{\rm KL}}_k>_{\mathcal {S}} Z^{{\rm KL}}_k(n)\right)^2\right\rbrace.\vspace*{3pt} \end{equation} \tag{ 13 }$$

A similar approach was proposed by Marois et al. (2010a) using LOCI, but it presents degeneracies because the LOCI coefficients are not guaranteed to be independent of $I_{\psi _0}$ or A. Such degeneracies do not arise in the case of KLIP because the effect of the algorithm on the astrophysical source is completely determined from Equation (13) and because the basis {Z^KL_k(n)}_{k = 1...Kklip} is independent of the science target. Note that Equation (13) assumes a noiseless signal. Optimal observers (Kasdin & Braems 2006; Caucci et al. 2007, 2009) can mitigate the impact of noise in the final estimate of ξ, and their implementation is greatly facilitated by the fact that Equation (13) only consists of propagating A_ξ(n) through a linear filter. Discussion regarding the practical implementation and performances of such observers are however beyond the scope of this study.

The left panel of Figure 3 shows an archival HST NICMOS image of HR 8799 reduced with KLIP using annular search zones. While the three planets b, c, and d are visible they are not detected unambiguously: the result reported by Soummer et al. (2011a) used a combination of thousands of LOCI reductions, and took advantage of two HST roll angles to minimize the false alarm probability. One could envision a similar parameter search with KLIP. The computational cost of such an exercise would however be reduced since KLIP does not require the introduction of "optimization zones" and "subtraction zones" (Lafrenière et al. 2007), which reduces the dimensionality of the algorithm's parameter exploration. More importantly, the astrometry reported in Soummer et al. (2011a) was the result of a lengthy characterization of the biases associated with LOCI, involving thousands of reductions with synthetic companions. On the right panel of Figure 3, we illustrated how one can determine the astrometric location and the photometry of HR8799bcd with forward modeling in a single KLIP-reduced image. In this example, we performed a simple χ² minimization of the position and flux of each planets, but this framework enables easy implementation of more appropriate modeling using Markov Chains Monte Carlo. KLIP thus provides astrophysical estimates which do not contain any bias arising from self-subtraction in the reduction algorithm. The uncertainties on the astrophysical parameters are solely driven by the local noise in the reduced image. This can be addressed by optimizing the KLIP parameters to increase local S/N, using several exposures and taking advantage of observers that are optimal with respect to the residual noise.

Zoom In Zoom Out Reset image size

Principal component analysis has been extensively used in astronomy to reduce the dimensionality of large data sets, e.g., Cowan et al. (2009) and Connolly et al. (1995). Principal components can also be used directly with LOCI or any of its variants to reduce the dimensionality of the reference set, using KL eigenimages instead of reference PSFs.

The cost function associated with KLIP (Equation (4)) is very similar to the cost function associated with LOCI (Lafrenière et al. 2007). However, LOCI minimizes the distance between the science target and the reference PSFs, whereas KLIP minimizes, in the statistical sense, the distance between any PSFs I_ψ in the close vicinity of references library and the K_klip-dimensional orthonormal basis-set {Z^KL_k(n)}_{k = 1...Kklip}. It is in principle possible to include the target to the reference library before decomposition into principal components, but in this case astrophysical sources can be completely subtracted and forward modeling cannot be carried out. While such solutions might yield slightly deeper contrasts, they might yield astrophysical estimates with biases that are difficult to calibrate. In this sense LOCI, with well-chosen parameters, might be optimal for detection while KLIP is more robust for characterization.

Because the set of {R_k}_{k = 1...N} is not an orthonormal basis, LOCI requires the inversion of the covariance matrix E_RR, which is typically ill conditioned. This leads to noise amplification if the inversion is not properly regularized. Several approaches to regularization exist, for example, eigenvalue truncation (Marois et al. 2010a) or Tikhonov regularization (Pueyo et al. 2012; Soummer et al. 2011a). With KLIP the regularization is implicit with the truncation of the KL basis. It is also more easily controlled because the KL set remains optimal for any value of K_klip. In addition, simple regularization schemes do not reduce the dimensionality of the covariance matrix, whereas KLIP reduces the number of references and therefore improves speed. Once the KL transform has been generated for a given search area geometry, the remaining free parameter K_klip can be explored with minimal computation burden (projections). Depending on the implementation of the algorithm and the type of data, KLIP can lead to considerable speed increases, especially when considering KLIP does not require us to process synthetic sources for characterization of detected signal.

As far as detection is concerned the most striking difference between LOCI and KLIP resides in disk imaging. Using HST Cycle 10 NICMOS archival data of the disk around HD 181327, we performed a global-LOCI and global-KLIP reduction over the entire image, using a carefully pre-selected reference library of 232 PSFs. The disk edges imaged with KLIP, Figure 4, are significantly better constrained when compared with state of the art classical PSF subtraction (Schneider et al. 1999). For comparison, we show in Figure 4 a LOCI reduction of the same object without any regularization: the self-subtraction degrades the rich astrophysical information contained in the ring's edge. It is important to note that with appropriate regularization it is possible to obtain with LOCI a very similar result to KLIP. However, with LOCI it is not easy to model the astrophysical signal since there is no guarantee that the algorithm configuration optimal for detection will not bias the astrophysical signal.

Zoom In Zoom Out Reset image size