EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 649 (May 2021) A&A, 649 (2021) A143 Full HTML
Open Access
Issue
A&A
Volume 649, May 2021
Article Number A143
Number of page(s) 13
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202140337
Published online 28 May 2021

© J. Rameau et al. 2021

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Direct detection of extrasolar planets relies on a combination of high angular resolution to spatially resolve the point-like planet to its host star at the subarcsecond level, high contrast to suppress the starlight, which is several orders of magnitudes brighter than the exoplanet, and optimized wavefront control to attenuate the diffracted starlight. The sensitivity to faint and/or close-in exoplanets is nevertheless plagued by residual starlight stemming from uncalibrated instrumental aberrations, and for ground-based facilities, Earth’s atmospheric turbulence. Disentangling the exoplanet from these so-called speckles is commonly achieved with differential observing techniques. Most of these strategies introduce relative spatial motion, azimuthal (ADI) and/or radial (SDI), between an astrophysical source and the speckle field (Racine et al. 1999; Marois et al. 2003). This spatial diversity is the central concept of post-processing algorithms of different complexities that model and subtract the speckles while preserving the signal of any planet (Marois et al. 2006; Lafrenière et al. 2007; Soummer et al. 2012; Amara & Quanz 2012; Gomez Gonzalez et al. 2016; Ren et al. 2018; Flasseur et al. 2018; Samland et al. 2020).

The ground-based extreme adaptive-optics systems Gemini/GPI (Macintosh et al. 2014), the Very Large Telescopes VLT/SPHERE (Beuzit et al. 2019), Subaru/SCExAO-CHARIS (Jovanovic et al. 2015), and the LBTI with ALES (Esposito et al. 2011; Skemer et al. 2015) make best use of the high performance of these key aspects to detect exoplanets and infer their individual and population properties (e.g., Macintosh et al. 2015; Chauvin et al. 2017; Keppler et al. 2018; Nielsen et al. 2019; Vigan et al. 2021). They are in particular equipped with an integral field spectrograph (IFS), providing two spatial and one spectral dimensional images, or hyperspectral data. The spectral range (Δλ ~ 0.3–1.5 μm) is large enough to benefit from the radial motion of the speckle field with the diffraction to improve its subtraction and hence the sensitivity to faint exoplanets (Mesa et al. 2015; Wang et al. 2015; Ruffio et al. 2017; Brandt et al. 2017; Galicher et al. 2018). The IFS further plays a key role in studying the atmosphere of the planets by providing a spectrum with low spectral resolution (R = λ∕δλ ~ 20–75) of each detection. At this resolution, broad molecular absorption bands can be identified, and the pseudo-continuum holds information about the pressure–temperature profile and clouds (e.g., De Rosa et al. 2016; Bonnefoy et al. 2016, 2018; Chilcote et al. 2017; Rajan et al. 2017; Delorme et al. 2017; Samland et al. 2017; Greenbaum et al. 2018; Müller et al. 2018; Chauvin et al. 2018; Currie et al. 2018; Uyama et al. 2020; Stone et al. 2020; Ward-Duong et al. 2021).

A higher spectral resolution (R > 5000–100 000) is a way to proceed in our understanding of the physical, chemical, and orbital properties of imaged exoplanets. Resolving the absorption lines augments the information content of the spectra and provides access to radial velocity (Snellen et al. 2014; Wang et al. 2018; Ruffio et al. 2019), spin (Snellen et al. 2014; Schwartz et al. 2016; Bryan et al. 2018), element abundances (Konopacky et al. 2013; Barman et al. 2015; Wilcomb et al. 2020; Petrus et al. 2021), or accretion tracers (Haffert et al. 2019). The combination of high spectral resolution and high contrast can also significantly enhance the sensitivity to faint exoplanets (Sparks & Ford 2002; Riaud & Schneider 2007; Kawahara et al. 2014; Snellen et al. 2015; Lovis et al. 2017; Wang et al. 2017). The technique assumes that the planet spectrum and that of the residual starlight can be separated through distinct absorption lines and/or different Doppler shifts that are due to orbital motion. The detection algorithm relies on straightforward cross-correlations of the signal with model spectra of the planet, which sum all common lines and cancel out nonplanetary features. The method has been proven to be very efficient even on instruments that were not designed to reach high contrast: the slit spectrographs VLT/CRIRES and Keck/NIRSpec (Snellen et al. 2014; Wang et al. 2018), and the medium-resolution IFS Keck/OSIRIS and VLT/Spectrograph for INtegral Field Observations in the Near-Infrared SINFONI (Barman et al. 2011, 2015; Konopacky et al. 2013; Hoeijmakers et al. 2018; Petit dit de la Roche et al. 2018; Ruffio et al. 2019; Wilcomb et al. 2020; Petrus et al. 2021). This strategy fuels the development of a new generation of instruments on existing 8- to 10-m-class telescopes (Kuntschner et al. 2014; Vigan et al. 2018; Mawet et al. 2018)and on future extremely large telescopes (ELTs; McGregor et al. 2012; Thatte et al. 2016; Brandl et al. 2018).

As of today, the cross-correlation technique has been applied on previously identified exoplanets. It was employed to extract new measurements of different planet characteristics that are accessible at medium to high spectral resolution. It has not been used for pure detection purposes. The prior knowledge of the planet that is searched for is a serious advantage in narrowing down the grid search of model spectra that are to be cross-correlated with the signal. Without such priors, the library should be large enough to encompass all possible spectral signatures, in particular to explore a wide range of effective temperatures and surface gravities that drive the emitted spectrum of a given molecule as well as radial velocities that control the positions of the lines. Still, the method relies on theoretical models that must resemble true planet spectra (e.g., Petrus et al. 2021). In this sense, it falls into the regime of supervised source-detection algorithms. It can also quickly become computationally time consuming for large libraries and/or large detectors. For these two reasons, unsupervised and fast alternative algorithms should be investigated to provide an alternative to the cross-correlation technique. The advent of large hyperspectral data with the IFS on board the JWST (Rieke et al. 2015; Bagnasco et al. 2007) and even more with ELTs-IFS such as the High Angular Resolution Monolithic Optical and Near-Infrared Integral field spectrograph HARMONI (Thatte et al. 2016), METIS (Brandl et al. 2018), GMTIFS (McGregor et al. 2012), and IRIS (Larkin et al. 2016), motivates this development.

Blind source-detection in hyperspectral data, or spectral unmixing for this type of data, represents a long-standing problem beyond astronomy, notably in Earth remote-sensing and microscopy. The dedicated algorithms aim at separating various sources in the field of view from their distinct spectral morphologies in a purely data-driven approach. Spectral umixing is a very active research field to develop powerful algorithms with growing complexities and specific features for a given science case, especially for Earth observations (see a review of the main methods from Bioucas-Dias et al. 2012, and references therein). In astronomy, spectral unmixing has been used to analyze IFS data from Spitzer, Herschel, or planetary probes. Some applications include velocity maps of interstellar clouds (Juvela et al. 1996), interstellar dust spectral inference (Rapacioli et al. 2005; Berné et al. 2007; Moussaoui et al. 2008; Gratier et al. 2017; Foschino et al. 2019; Boulais et al. 2021), or Mars surface mapping (Forni et al. 2005; Hauksdottir et al. 2006; Themelis et al. 2012; Liu et al. 2018). Beyond these examples, methods like this are rare in astronomy to our knowledge despite the amount of existing and forthcoming hyperspectral data, and they have never been tested for exoplanet science.

In this paper, we explore the use of unsupervised spectral umixing as an alternative to the supervised cross-correlation technique to directly detect an exoplanet and extract its spectrum from medium-resolution data. We test one algorithm on synthetic R ~ 7000 K-band HARMONI data to evaluate its performances. We then apply it on on-sky R ~ 5000 K-band SINFONI data of β Pictoris and show that spectral umixing can effectively detect β Pictoris b in less than a minute on a regular laptop. The paper is organized as follows. In Sect. 2, we motivate the unmixing approach for direct exoplanet detection in hyperspectral data and provide the mathematical framework of the adopted algorithm. In Sect. 3, we describe the simulations of the HARMONI data and demonstrate the algorithm on a simple test case, after which detection sensitivity and spectral fidelity are discussed and compared to the cross-correlation technique in Sect. 4. In Sect. 5, the algorithm is applied on the SINFONI data. Advantages and limitations of the algorithm are discussed in Sect. 6, followed with concluding remarks and proposals for improvement.

2 Source separation from spectral decomposition

2.1 Statement of the problem

A typical hyperspectral exoplanet imaging dataset contains two sources: the star and the orbiting planet. The two astrophysical objects have different Teff, from 2800 K up to >10 000 K for the former, and <2000 K for the latter. The temperature drives the formation of elements in their atmospheres (hydrogen- and metal-bearing ions, atoms, and molecules for stars, and H2O, carbon- andnitrogen-bearing molecules for planets), which hence plays a key role in the appearance of the emitted spectra. As illustrated in Fig. 1 in the K band, the two spectra are distinct at low frequency1 by the blackbody emission and the broad absorption bands and at high frequency by the various absorption lines. This dissimilarity can be exploited to separate the two sources in the hyperspectral data.

While the planet is spatially localized in the image, the stellar spectrum is spatially replicated over the image via the diffracted residual starlight. A given spaxel di in the hyperspectral data can be modeled with a linear statement, such that di = ∑j = 1psjαj,i+ϵi,\begin{equation*} \bm{d}_{i}\,{=}\,\sum_{j\,{=}\,1}^{p}\bm{s}_{j}\alpha_{j,i}&#x002B;\bm{\epsilon_i},\end{equation*}(1)

where sj is the spectrum of the jth source, αj,i is the corresponding weight at the ith spatial location, and ϵi is an additive noise term. When we assume that the source spectra are spatially invariant in the image, that is, all pixels can be expressed a linear combination of the same p spectra, Eq. (1) can be rewritten in compact matrix form as D = SA+E,\begin{equation*} \bm{D}\,{=}\,\bm{S}\bm{A}&#x002B;\bm{E},\end{equation*}(2)

where D∈ℝnλ × m${\bm{D}} \in\mathbb{R}^{n_{\lambda}\,{\times}\,m}$ is the column-vectorized hyperspectral images where the columns contain the m = nx × ny spaxels in any order of nλ spectral channels, S∈ℝnλ × p${\bm S}\in\mathbb{R}^{n_{\lambda}\,{\times}\,p}$ is a matrix of p source spectra, A∈ℝ+p × m${\bm A} \in\mathbb{R}_&#x002B;^{p\,{\times}\,m}$ is the matrix of their corresponding spatial weight maps, and E the uncorrelated noise. Equation (2) is the typical linear mixture model assumed for spectral unmixing. S contains the so-called endmembers, and A contains theabundance maps (Bioucas-Dias et al. 2012).

The decomposition requires three steps: (i) the estimation of the dimensions of S and A, that is, the number p of spectra in the data, (ii) the extraction of the spectra, and (iii) the evaluation of the weight maps. Each step is detailed in the following.

thumbnail Fig. 1

Model spectra of a β Pictoris-like system: an A6 star with Teff = 8000 K (BT-Nextgen Allard et al. 2012) and a giant planet with Teff = 1700 K (BT-Settl Allard 2014), normalized over the K bandpass at a resolution of R = 7000. Not only the slope of the continuum is reversed between the two objects, but also the absorption line content, with the prominent Brackett-γ line in the stellar spectrum at 2.16μm, while the planet spectrum exhibits a forest of H2O and CO lines.

2.2 Estimation of the number of sample spectra

The initial hypothesis of only two sources, the star and the planet, does not hold in high-contrast data. The speckles act as modulators of the stellar spectrum as a function of position and wavelength, as illustrated in Fig. 2 for synthetic HARMONI data (see Sect. 3 for details on the simulation). Therefore, the effective number of samples in our hyperspectral data can be as high as the number of speckles in the field of view. This number might further be augmented by the spatially variable imperfect calibration of the data (wavelength, crosstalk, and bad pixels) or line spread function, as in VLT/MUSE data (Xie et al. 2020)2.

Estimating the number of sources in hyperspectral data is a difficult task in spectral unmixing and subject to regular algorithm development. Several techniques have been proposed to exploit the fact that the inherent structure of high-dimension hyperspectral data usually lies in a lower-dimensional subspace (e.g., Chang & Du 2004; Bioucas-Dias & Nascimento 2008; Acito et al. 2009; Luo et al. 2013; Deville et al. 2014; Halimi et al. 2015; Drumetz 2016). These methods were designed for the low-contrast regime, however, and are sensitive to outliers, noise, and spectral variation across the field of view (referred to as variability) (Drumetz et al. 2020, and references therein).

Because of these limitations and because only the extraction and identification of the planet component matters, the maximum number of effective samples p is at the moment set empirically by the user. Because the algorithm has a low time complexity (see below), we recommend running a first trial with a large p of 15–20 or more in case of large images (making sure p ≪ nλ, and p ≪ m = nx * ny, for the algorithm towork). p can then be adjusted to the minimum number of samples necessary to detect any planet in the data. For the data set used in this work, which was tested with a range of contrasts and separations, all planets were extracted from the first 2–10 components. Therefore p = 10 was chosen for synthetic HARMONI and real SINFONI data in the following; the decomposition into a larger number of samples results in additional starlike spectra, which is not necessary for planet detection. Because the algorithm is very fast (see below), it is not necessary to explore complex algorithms suitable for estimating p for our data. We therefore defer this to a future work.

thumbnail Fig. 2

Top: spectra of individual spaxels in synthetic coronagraphic HARMONI data at R ~ 7000 with a central A6 star and a 1700 K planet (orange) (see Fig. 1). Random spaxels (gray) in the corrected region illustrate the modulation by the speckles of the stellar spectrum, as shown with the central spaxel behind the semitransparent focal plane mask (blue). Details of the simulations are given in Sect. 3. Bottom: resulting spectra after the continuum was removed (following Hoeijmakers et al. 2018). The amplitude of the planet spectra has been amplified for visualization purposes.

2.3 Orthogonal subspace projection to extract the sample spectra

The extraction of the sample spectra is a very active topic in remote-sensing research; many strategies have been investigated based on geometrical or statistical considerations (see the reviews by Plaza et al. 2004; Bioucas-Dias et al. 2012, and references therein). We considered a geometry-based algorithm, which is the most frequently used algorithm family for spectral unmixing applications. This approach finds the set of the purest spaxels in the data for each requested sample.

Except for planets that are as bright or brighter than the stellar halo, this assumption does not hold at high contrast. We circumvented the problem by subtracting the halo in the spectral dimension with data preprocessing. As discussed before, the halo mostly affects the broadband flux and the continuum of the spaxels. It is therefore a low spectral and spatial frequency component that can be removed with several strategies such as high-pass filtering (Konopacky et al. 2013; Ruffio et al. 2019), continuum normalization (Hoeijmakers et al. 2018; Haffert et al. 2019; Petrus et al. 2021), or principal component analysis (Xie et al. 2020). Removing the low frequencies of the stellar halo in the spectral dimension also minimizes the variation in the stellar spectrum across the field of view, which in turn reduces the effective number of samples. As a consequence, this procedure also removes the continuum of the planet spectrum (see the example following Hoeijmakers et al. 2018 in Fig. 2). However, the distinct and multiple absorption lines between the different spectra enables spectral unmixing in the high-frequency components.

Following Eq. (2) and with p ≪ nλ, all spaxels live in a low-dimensional subspace, or simplex, whose vertices (geometrical corners) are represented by the sample spectra. A given vertex is therefore most distant to the subspace spanned by the other p − 1 vertices. An orthogonal subspace projection (OSP) approach can be used to identify the vertices (e.g., Harsanyi & Chang 1994; Chein-I Chang 2005).

Let s0 be the first sample spectrum. The sparsity of the planet over the diffracted starlight in the field of view motivates the choice of s0 as the spectrum of the median spaxel. Let U be the spectral sample matrix, which at this step only contains s0, such that U = [s0]. The orthogonal projector from U, denoted PU⊥$P_{\bm{U}}^{\perp}$, is defined as PU⊥ = I−U(UTU)−1UT,\begin{equation*}P_{\bm{U}}^{\perp}\,{=}\,\bm{I}-\bm{U}(\bm{U}^{\textrm{T}}\bm{U})^{-1}\bm{U}^{\textrm{T}}, \end{equation*}(3)

with I a nλ × nλ identity matrix, UT the transpose of U, and (UTU)−1UT$(\bm{U}^{\textrm{T}}\bm{U})^{-1}\bm{U}^{\textrm{T}}$ the pseudo-inverse of U. The dot product PU⊥di$P_{\bm{U}}^{\bot}\bm{d}_{i}$ projects each spaxel di in 〈s0〉⊥$\langle\bm{s}_{0}\rangle^{\bot}$, the subspace orthogonal to the space linearly spanned by s0 3. Iterating over all spaxels, a second sample spectrum, set as s1, is flagged as the most distant spaxel to s0. The projection distance is chosen to be the sum of the squared residuals after removing from di the projection of di onto s0, or (PU⊥di)T(PU⊥di).\begin{equation*}(P_{\bm{U}}^{\perp}\bm{d}_{i})^{\textrm{T}} (P_{\bm{U}}^{\perp}\bm{d}_{i}). \end{equation*}(4)

A new orthogonal subspace projector PU⊥$P_{\bm{U}}^{\perp}$ can be defined from the linear subspace spanned with U = [s0, s1] and applied to the remaining spaxels. The third sample spectrum is identified as the one with the maximum absolute orthogonal projection in the subspace ⟨s0,s1⟩⊥4. The process is repeated until p sample spectra are exhausted to build up S. To prevent bad or hot pixels from being flagged as sample spectra because they differ by construction from the remaining spaxels, they must be properly corrected for during data preprocessing.

2.4 Evaluation of the weight maps

From the extraction of the matrix S of sample spectra, the determination of the spatial weight matrix A from Eq. (2) is an inverse problem, which can be solved with least-squares means. The problem can be further constrained because all spaxels in our hyperspectral data are a positive linear combination of the sample spectra, that is, A ≥ 0. Therefore the non-negative least-squares estimate of A is done for each spaxel di by finding argminαi||Sαi−di||2,subject to αi≥0,\begin{equation*} \mathrm{arg}\min_{{\alpha_i}}\vert\vert \bm{S}\alpha_{i}-\bm{d}_{i}\vert\vert_2, \mathrm{subject}\medspace\mathrm{to}\medspace{\alpha_i}\geq0, \end{equation*}

where ||⋅||2 denotes the Euclidian norm.

We used the optimize.nnls module from scipy to solve the problem. In practice, this iterative algorithm is very slow to converge. Because in our case the number of coefficients is much smaller than the number of wavelength channels, that is, p ≪ nλ, the solution proposed by Bro & De Jong (1997, proofs therein) can be used such that solving the non-negative least-squares problem is accelerated by replacing S with a p × p matrix following argminαi||(STS)αi−STdi||2,subject to αi≥0.\begin{equation*}\mathrm{arg}\min_{{\alpha_i}}\vert\vert (\bm{S}^{\textrm{T}}\bm{S}){\alpha_i}-\bm{S}^{\textrm{T}}\bm{d}_{i}\vert\vert_2, \mathrm{subject}\medspace\mathrm{to}\medspace{\alpha_i}\geq0. \end{equation*}(5)

Repeating the resolution over the entire image builds A. Reshaping the matrix into the original 2D spatial dimensions finally generates p weight maps corresponding to the decomposition of the input hyperspectral image D.

2.5 Planet detection and summary of the algorithm

The weight maps provide the spatial distribution of each sample in the data in an undetermined way; the map associated with the planet component has to be identified a posteriori. The spatial signature of the planet is defined by the point-spread function (PSF) with a finite spatial extent (the full width at half maximum, FWHM) that is typically larger than a single pixel. Planet vetting can be performed visually because the number of maps is usually small, or it can be conducted automatically using a matched filter (Cantalloube et al. 2015; Ruffio et al. 2017).

The algorithm based on spectral unmixing for detecting a planet in hyperspectral data is summarized below:

  • 1.

    Set p as the number of effective sample spectra in the data (p > 2 due to spatial variability). p = 15–20 might be a good first guess to prevent missing any planet, while a higher value is recommended for very large images. For the data explored in this work, p = 8–10 is appropriate;

  • 2.

    Remove the continuum for every spaxel to keep the high-frequency components;

  • 3.

    Vectorize the hyperspectral data D in a 1D spatial plus 1D spectral matrix;

  • 4.

    Set the first sample s0 as the closest spectrum to the median spaxel and define U = [s0] the subspace spanned by s0;

  • 5.

    Project each spaxel onto the linear subspace spanned by U with the orthogonal projector PU⊥$P_{\bm{U}}^{\perp}$ defined in Eq. (3);

  • 6.

    Find thespaxel with the extremal projection (see Eq. (4)), which corresponds to a new sample spectrum, snew;

  • 7.

    Append U with snew;

  • 8.

    Repeat steps 5–7 until all p spectra are extracted;

  • 9.

    Invert D with fast non-negative constrained least squares (see Eq. (5)) to compute the spatial distribution of each sample spectrum A;

  • 10.

    Expand A in p 2D spatial maps to search for a planet.

3 Application on synthetic data

In this section, the capability of the proposed unsupervised spectral unmixing algorithm is first tested on synthetic high-contrast hyperspectral HARMONI data. Then we compare it with classic cross-correlation.

3.1 Simulations of ELT/HARMONI data

HARMONI will be the first-light IFS on the ELT offering diffraction-limited hyperspectral images at medium spectral resolution (3500–18 000) in the optical and near-infrared (Thatte et al. 2016). A module consisting of an apodized pupil coronagraph and a ZELDA wavefront sensor will supplement the single conjugate adaptive optics of the core instrument and give access to high contrast (> 10−5) at very short separations (≥50 mas) in the Nyquist-sampled 0.61″ × 0.86″ field of view with a spatial scale of 4 mas pixel−1 (Carlotti et al. 2018).

A synthetic coronagraphic sequence is computed with a dedicated end-to-end simulator with the following configuration: 4 h of observations with one-minute-long exposures in angular differential imaging for a star at a declination of − 15 deg under median seeing condition (0.7–0.8″) at K band (1.951–2.469 μm, R = 7104) with the shaped pupil SP1 (working angles in the 5–20 λ∕D range) and the partly opaque focal plane mask FPM1. This enables direct calibration with the attenuated (10−4) central star.All details of the simulator are presented in Carlotti et al. (2018) and in Appendix A.

Planets can further be injected into the mock data cubes. The off-axis simulated PSF is used as a template for a simulated planet, and a BT-Settl model spectrum (Allard 2014) with an effective temperature of 1700 K and log g = 4.0 dex is assumed by default in the following, unless specified. These parameters are typical of known L-type young giant exoplanets (HIP 65426 b, β Pictoris b). For simplicity, no Doppler shift is assumed. The planet has the same radial velocity as the central star. Contrast and atmospheric transmission are calibrated with the attenuated central star. These settings are kept constant for all the following simulations, only contrast and injected position vary. To speed up the injection and recovery processes in the following, a 2 h subset of the simulated data was used and the central 100 × 100 pixels were cropped in each frame, resulting in 100 × 100 × 1665 × 120 cubes for a total of 16 Gb.

3.2 Single test case

A test case was generated to give a first impression of our spectral unmixing algorithm. A simulated planet with a contrast of 5 × 10−6 was injected at 100 mas and a position angle of 315° into the simulated data cube. For each cube, the stellar halo was subtracted in the spectral direction in each data cube following Hoeijmakers et al. (2018), as suggested in Sect. 2. The frames were then aligned north along the vertical axis and averaged over the temporal axis. The final 3D data cube was decomposed into eight sample spectra and their corresponding spatial maps through spectral unmixing following the steps given in Sect. 2. To compare the result with a spectrum-based detection algorithm, each spaxel in the final 3D data cube was cross-correlated with the injected model spectrum, whose pseudo-continuum was subtracted beforehand with a convolution with a Gaussian line-spread function corresponding to a spectral resolution of 100. The cross-correlation used a customized version of CrosscorrRV5 to include normalization with the variance of the model and the data. A perfect match yields 1.

Figure 3 shows the final images using the cross-correlation and spectral unmixing algorithms. Spectral unmixing can clearly separate the speckles from the planetary signal. The field of view is further decomposed into several regions that correspond to the variation of the stellar spectrum around the outer working angle, partial coverage of some pixels by the asymmetric focal plane mask after derotation and stack, or interpolation effects close to the mask (dominantlyyellow regions). The pixels at the location of the planet have zero weight in the corresponding maps. We computed the S/N of the planet and the background pixels at the same separation, following the definition with small sample statistics correction of Mawet et al. (2014), to assess the detection. We used an aperture size of 2.6 pixels, which is the average FWHM measured on the off-axis PSF. The S/N with the cross-correlation is 23.6, while the S/N with the spectral unmixing is 12.6. The distribution of the values in the cross-correlation and the weight maps are given in Appendix B. If the decomposition were exact, all pixels in the background would have zero weights, which would cause the S/N to become infinite. Figure 3 shows the extracted spectrum corresponding to the planet, compared with the injected model. The match is visually reasonable, in particular the bandheads of the first overtone of CO at 2.295 μm and at 2.324 μm, which are well recovered. The distribution of the difference between the two spectra as shown in the bottom panel of Fig. 4 is Gaussian and almost centered on zero (−0.28). Residual noise from the subtraction of the diffracted halo contaminates the inferred spectrum of the planet by construction. These residuals are in common with that of the spaxels in the background at the same separation, explaining their small but nonzero weights. To further show that the extracted spectrum holds information about the true planet spectrum, the amplitude of the cross-correlation with the injected model is 0.15, which is close to the amplitude (0.17) at the planet location in the cross-correlation map, as expected. Cross-correlating the p − 1 nonplanetary spectra with the injected model yields 0.025 ± 0.010 on average. This confirms that the planetary information content is gathered in a single spectrum.

The cross-correlation leads to an S/N that is about twice as high as with spectral unmixing in this test case. However, the case is ideal for supervised detection techniques such as a cross-correlation because the searching template is the same as the injected one. Only contrast and noise matter. The lower S/N with the unsupervised spectral unmixing approach does not prevent the detection of the planet, and the extracted spectrum can be analyzed directly. In Sect. 4, we test whether this holds true with fainter planets. Nevertheless, spectral unmixing provides a prior-free alternative to classic cross-correlation.

thumbnail Fig. 3

Final result of post-processing the synthetic HARMONI data with a cross-correlation with the planet spectral template (left) and spectral unmixing through orthogonal subspace projection (middle and right). The number of spectral signatures is set to eight to decompose the data. The injected mock planet has a broadband contrast of 5 × 10−6 with the central star at a separation of 100 mas (white circle). Images are normalized to the peak value of the planet for each algorithm, the scale is linear, the field of view is identical, and the white cross marks the position of the star.
thumbnail Fig. 4

Inferred spectrum of the planet with a contrast 5 × 10−6 injected into the synthetic HARMONI data (see Fig. 3 with spectral unmixing through orthogonal subspace projection (dark blue), compared with the high-frequency component of the injected spectrum (orange) along with the difference between the two (black). The first-overtone bandheads of CO are highlighted as being well recovered (shade). The unmixing is not perfect, and the spectrum of the planet is contaminated with residuals from the subtraction of the diffracted starlight (bottom).

4 Performance

4.1 Method

In order to explore the properties of spectral unmixing in more detail and compare it with classic cross-correlation beyond a single test case, we computed the detection sensitivity in the central region, which is expressed in the star-to-planet contrast as a function of the separation. Simulated planets were injected from 68 to 132 mas with steps of 12 mas (~ 1 λ∕D) every 60° along a spiral. The injection was independently repeated six times by rotating the spiral by 60° increments. For each set of injections, the contrast was enhanced from 1 × 10−7 to 1 × 10−5, with steps of 5 × 10−7 first and then refined with steps of 1 × 10−7 around the limits. The detection limit is built upon the minimum contrast that yields 95% of the injected simulated planets at a given separation, to be detected above a threshold of 5 in S/N (Jensen-Clem et al. 2018). The process was repeated with spectral unmixing with p = 10 and a cross-correlation with the same template.

4.2 Effect of contrast and separation

The 95% completeness curves are shown in Fig. 5. Spectral unmixing is able to both flag planet-like spectra and lead to significant planet detections down to 4 × 10−6 in this set of synthetic HARMONI data. The detection limit achieved with ideal supervised cross-correlation is more than four times deeper than with spectral unmixing at all separations. The question then arises whether the lower limit with spectral unmixing is reached because the S/N is below the adopted threshold or because the algorithm fails to glean the planet spectrum. The S/N hits a hard floor around three down to a contrast of 3 × 10−6. Fainter planets cannot be flagged in any of the spatial weight maps corresponding to the p = 10 identified sample spectra; augmenting p up to high values does not help. This critical step therefore limits the performance of spectral unmixing. To corroborate this, a data set with planets with a contrast 9 × 10−7 was decomposed into p = 10 samples. The planet model template was further added as an eleventh spectrum, and the data were inverted following Eq. (5). The corresponding spatial weight map and its S/N map are displayed in Fig. 6. All six planets can be spotted in the maps, with a high S/N for all but the two innermost ones. The weight corresponding to each planet is lower than 5 × 10−3%, however, meaning that these spaxels are highly mixed and that the planetary component is too weak to be extracted with Eqs. (3) and (4).

Spectral fidelity can also be investigated as a function of contrast. The planet spectrum as extracted with spectral umixing was cross-correlated with the injected template. The amplitude linearly decreases from 0.15 to 0.07 with contrasts from 10−5 to 3 × 10−6. This is more than three times than the amplitude of the cross-correlation of all other 10− 1 spectra with the planet template, which is stable around 0.027 ± 0.005.

thumbnail Fig. 5

Completeness curves (95%) with spectral unmixing (dark blue) and cross-correlation (brown) expressed as the minimum planet-to-star contrast for a detection threshold of τ = 5σ as a function of the projected separation for L-type (solid) and T-type spectra (dot-dashed) in the synthetic HARMONI data.
thumbnail Fig. 6

Spatial weight inversion map of synthetic HARMONI data with a spiral of simulated planets (white circles) with a star-to-planet flux ratio of 9 × 10−7 (top) and corresponding S/N map (bottom). The injected model spectrum is used to invert the data following Eq. (5). Despite tiny spatial weights, most planets are detected with a high S/N.

4.3 Effect of the planet spectrum

The previous results hold for the planet spectrum at 1700 K, or L-type planets. However, changes in the intrinsic properties of the absorbing molecules (e.g., temperature, pressure, or abundance) alter the number and depth of the lines in the planet spectrum. As the S/N achieved with cross-correlation at first order scales with these two parameters (Snellen et al. 2015), we might expect that they might also help pushing the limits with spectral unmixing.

We useda significantly different model spectrum with an effective temperature of 700 K. This is illustrative of known T-type planets (e.g., 51 Eri b). The atmosphere is now cool enough such that methane can form and produce many absorption lines, particularly in the K band, in addition to the lines due to water and carbone monoxyde. The injection and recovery process was repeated with this new template to build the 95% completeness curves for spectral unmixing and the cross-correlation.

The resulting curves are shown in Fig. 5. Spectral unmixing can detect fainter T-type planets than L-type planets, the limiting contrast being improved by 55% down to 1.8 × 10−6. The sensitivity is also boosted with the cross-correlation, as predicted. T-type planets with contrast as faint as 5 × 10−7 can be detected,which is an improvement of 45% compared with L-type planets. A switch from L-type to T-type spectra yields a higher gain for spectral unmixing and thus reduces the gap between the two methods from 4.4 to 3.6. This difference might be even further reduced for real planets. Model spectra do not perfectly reproduce observed spectra, therefore template mismatch will degrade the performance of the cross-correlation, while spectral unmixing should remain unchanged.

It has to be noted that the completeness curves presented here do not predict the ultimate sensitivity of HARMONI. They explore a very limited range of planet properties and post-processing algorithms. This is not the goal of the present work.

5 Validation on real VLT/SINFONI data

The spectral unmixing algorithm was further vetted with real VLT/SINFONI data. Although SINFONI was not designed as a high-contrast instrument, it delivers medium-resolution hyperspectral data such that our method can be used.

5.1 Data and preprocessing

The data set of β Pictoris and its planetary companion β Pictoris b (Lagrange et al. 2010) obtained on September 10, 2014, was used (Hoeijmakers et al. 2018; Bonnefoy et al., in prep.). Observations were made in angular differential imaging mode at K band (1.929−2.472 μm, R ~ 5000) with the high spatial resolution mode. This provided a 0.8″ × 0.8″ field of view with a scale of 12.5 × 25 mas2 pixel−1. The adaptive optics system was locked on the primary star, yielding a Strehl ratio of 17–29% under median seeing conditions (0.7–0.9″). The star was kept outside the field of view during the sequence in order to prevent saturation. Twenty-four cubes with an integration time of 60s were obtained. Data cubes were initially built from the 2D raw detector images with the SINFONI data-handling pipeline v.3.2.3 (Abuter et al. 2006). The pipeline calibrated and corrected the instrument transmission and the distortion, identified and interpolated hot and nonlinear pixels, retrieved the position of the 32 dispersed slitlets in the raw frames, achieved the wavelength calibration, and finally built the 3D spatial-spectral cubes. We used the TExTRIS package (Petrus et al. 2021; Bonnefoy et al., in prep.) prior to the ESO pipeline to correct the raw frames for various electronic noises. TExTRIS was also used to correct the data cubes for the improper estimates of the slitlet edges provided by the pipeline and for residual wavelength shifts. It computed the parallactic angles and retrieved the position of the star outside of the field of view. The diffracted starlight was then removed from each cube following Hoeijmakers et al. (2018). All cubes were aligned with north up and averaged at each wavelength slice to produce a single hyperspectral data cube. Further details can be found in Petrus et al. (2021) and Bonnefoy et al. (in prep).

5.2 Results

β Pictoris b (known contrast of 2.1 × 10−4, separation of ~360 mas at the time of the observations, Bonnefoy et al. 2011; Lagrange et al. 2019) was searched for in the hyperspectral data cube with cross-correlation and spectral unmixing. The cross-correlation was performed over a range of radial velocities of ± 100 km s−1 with respect to the star with steps of ±10 km s−1. The template spectrum was a continuum-removed BT-Settl model that matched the effective temperature (~ 1700 K) and surface gravity (~4.0 dex) of the planet (Chilcote et al. 2017; Nowak et al. 2020). For the spectral unmixing, the data were assumed to be decomposed into at least p = 10 sample spectra, in the same way as for the HARMONI data. Figure 7 shows the resulting images with the two algorithms. Spectral unmixing clearly unveils the planet, while the majority of the background pixels have zero weight with the planet spectrum. The CCF peaks at the planet position in the central velocity bin (0 km s−1), as found by Hoeijmakers et al. (2018). Because the reconstructed field of view is asymmetrical, a pseudo-S/N was computed from the mean flux and standard deviation of all resolution elements in the whole field of view, with a diameter of 4 pixels. β Pictoris b has a pseudo-S/N of 28.5 with spectral unmixing, and this value is 8.2 with the cross-correlation.

The spectrum corresponding to β Pictoris b in the spatial weight map in Fig. 7 is displayed in Fig. 8. A visual comparison can be established with the BT-Settl model spectrum. Although it is not a perfect model of the planet, several common lines can be spotted, including two bandheads of the first overtone of CO at 2.295 μm and at 2.324 μm. The histogram of the residuals does not exhibit aprominent bias (μ = −0.15 ADU), but shows a large spread (σ = 2.98 ADU). This is supported by the cross-correlation of the model spectrum with the extracted planet spectrum, which yields 0.33 (the cross-correlation with the seven other extracted spectra gives 0.07 ± 0.05). However, the region beyond 2.34 μm is noticeably polluted by residuals that are not properly decomposed from the planet spectrum. Again, this is expected given the hypothesis of spectral unmixing, which is discussed in Sect. 6.

thumbnail Fig. 7

Final results of post-processing VLT/SINFONI K-band data of β Pictoris b with the cross-correlation (top) and spectral unmixing through orthogonal subspace projection (bottom). A BT-Settl model spectrum (1700 K, log g = 4.0) was used for the cross-correlation, and p = 10 was set to decompose the data. The scale is linear and was chosen to reach 90% of the peak value for each image, the cross-correlation strength Cs, and the weight αi, respectively. The cross marks the position of β Pictoris A.

6 Discussion and concluding remarks

6.1 Limitations and benefits

Several limitations to spectral umixing have been identified throughout this work. First, the algorithm fails to detect a planet as faint as with the classic cross-correlation in synthetic data, which is favored by the injection and recovery process using the same model spectrum. However, as experimented on the real SINFONI data, for which it led to a higher S/N for β Pictoris b, the sensitivity gap between the two methods may be reduced because the cross-correlation is affected by template mismatch. The extraction steps were shown to be responsible for this (see Sect. 4). The orthogonal projection distance was used as a metric to quantify spectral dissimilarities between the tested spaxels. Alternative distances might be used. Cross-correlation and spectral angle mapper, a common metric in hyperspectral remote sensing, defined the angle between two vectors formed by two spectra (Kruse et al. 1993), were tested but they led to worse performances (see Appendix C). For faint planets, the spectrum of the planet has a very small weight on the corresponding spaxels compared to the residual diffracted starlight even after the continuum was removed. In hyperspectral remote sensing, they are commonly treated as highly mixed spaxels. The identification of sample spectra in this case is a very difficult task for which a geometry-based algorithm such as we adopted, is not robust enough (Bioucas-Dias et al. 2012). As hypothesized in Sect. 2, the algorithm flags the source sample spectra as the purest spaxels in the field of view, which are then used to decompose the remaining spaxels. With highly mixed spaxels at the planet location, the hypothesis is no longer valid to flag and unmix the planet. These phenomena also lead to a second limitation of this approach. As we showed in Sects. 4 and 5, the extracted spectrum that leads to the planet detection in the spatial weight maps, is not the true planet spectrum. The decomposition is not exact (the true value of p remains unknown), the diffracted starlight contaminates the spectrum. As a consequence, the true contrast of a planet cannot be obtained with spectral unmixing and needs another post-processing approach.

Beyond these limitations, the spectral unmixing algorithm presented in this work has two advantages. It is fully data-driven, or unsupervised, and computationally fast. Both characteristics make it interesting with respect to the other strategy for analyzing hyperspectral data, the classic cross-correlation. For a purely blind search of real data, the cross-correlation must rely on theoretical planet spectra and explore many parameters for the planet(s) that are to be searched for, including the effective temperature, the surface gravity, the metallicity and composition, and the radial velocity. This will become even more dramatic with the forthcoming ELT/HARMONI, METIS, and JWST/MIRI, which should be sensitive to young Saturn- and Neptune-like planets as well as mature Jupiters. The spectral diversity of these populations of planets is poorly known and might deserve an exploration of even larger parameter space. Expressing a robust detection limit that encompasses all these possibilities might therefore be challenging and a template mismatch might lead to false negatives. While a single cross-correlation is computationally very fast, the process can take several minutes to hours for the whole field of view. For huge amounts of data such as those provided by the forthcoming HARMONI and METIS instruments and/or with large template libraries, the task quickly becomes computationally inefficient. This is not the case with spectral unmixing, which takes 15 s to process the 16 Gb of the synthetic HARMONI data on a standard laptop with a dual-processor core running at 2.4 GHz.

thumbnail Fig. 8

Top: extracted continuum-removed SINFONI spectrum of β Pictoris b with spectral unmixing (dark blue), compared with a BT-Settl model spectrum at 1700 K and log g = 4.0 dex (orange), close to the true planet properties. The first-overtone bandheads of CO are highlighted as being well recovered (shade). The extracted spectrum is not perfect, it is contaminated by unmixed residuals. Bottom: histogram of the residuals (black) with a best-fit Gaussian function (purple).

6.2 Summary and future work

To conclude, we have introduce a new strategy for analyzing high-contrast medium spectral resolution data aiming at directly detecting and characterizing exoplanets. We used the dissimilarities between the spectrum of a planet and that of the remaining field of view, which are mostly variations in the spectrum of the central star. The approach comes from remote-sensing research, in which thespectral diversity is a driver of algorithm development, namely all methods in the class of spectral unmixing. A classic algorithm based on orthogonal subspace projection and non-negative least-squares inversion was proposed. It was demonstrated to be viable to simultaneously detect a planet and extract spectral information content in synthetic and real hyperspectral data. This approach is fully data-driven as opposed to a cross-correlation that relies on theoretical templates, which is commonly adopted to analyze such data. The unsupervised nature of spectral unmixing comes at the price of a reduced sensitivity compared to the supervised cross-correlation in the ideal case of synthetic data with the adopted stellar and planet models, but this may not be the case for different settings and real data. Nevertheless, the algorithm benefits from low computational cost, making it suitable for real-time processing in a purely blind search scenario.

The spectral unmixing algorithm presented in this work is solely based on spectral dissimilarities between spaxels in the field of view. This basic assessment is enough to demonstrate the potential of this class of strategies in decomposing hyperspectral data aiming at simultaneously detecting and characterizing exoplanets. Additional constraints and more expensive formulations of the decomposition based on the knowledge of the data are possible. Spectral umixing is a very active field of research, therefore we point out only a few directions of high interest. Spatial regularizer, taking into account the size of the PSF, can be added because nearby pixels share a similar spectral signature, and the decomposition of neighboring pixels might be forced to be comparable for the same sample spectra (e.g., Zhang et al. 2018). Furthermore, the planet signal is sparse in spatial and spectral dimensions, while the remaining field of view is low dimensional in the spectral direction. Low-rank and sparse unmixing can thus be foreseen (e.g., Tsinos et al. 2017; Xu et al. 2018), similar to the reasoning behind the decomposition of angular differential imaging data from Gomez Gonzalez et al. (2016). The sparsity constrain can be pushed even further. The planet can be considered as a spectral anomaly in the field of view, or target, and anything else may be considered as background (e.g., Li et al. 2015; Soofbaf et al. 2018; Yang et al. 2019). These improvements could lead to better sensitivity and unbiased spectral extraction. They will be the focus of future work.

Acknowledgements

J.R. is supported by the French National Research Agency in the framework of the Investissements d’Avenir program (ANR- 15-IDEX-02), through the funding of the “Origin of Life” project of the Univ. Grenoble-Alpes. A.C. acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program, for the Project “EXACT”.

Appendix A Simulations of ELT/HARMONI synthetic data

In this appendix, we briefly summarize the main characteristics of the simulation of the synthetic hyperspectral data with ELT/HARMONI. The single conjugate adaptive optics subsystem provides a Strehl ratio of 75–80% at K band, leaving aside residuals that also take the primary mirror aberrations into account (island effect, wind shake). Additional aberrations come from the elements in the optical train toward the focal plane mask, including chromatic beam shift (due to the absence of atmospheric dispersion correction in HARMONI) and rotation of the optics in the relay system (HARMONI is mounted on the Nasmith platform). The high-contrast module yields 80 nm RMS residual aberrations in this configuration on average. The simulator produces coronagraphic PSFs for each exposure and wavelength, which are assembled in 240 data cubes of 215 × 215 × 1665 pixels for a total volume of 70Gb. An off-axis PSF is also generated.

thumbnail Fig. A.1

Simulated coronagraphic ELT/HARMONI image at 2.0 μm for an exposure of 60 s with the SP1 apodizer. The central star is attenuated by the focal plane mask, which is asymmetric to take the atmospheric dispersion into account. The remaining field of view is plagued by diffracted starlight. The white line marks the adaptive optics corrected area, the region where the speckle intensity is minimized. The intensity scale is linear.

Astrophysical signal, throughput, background, and random and systematic noise were incorporated into the optical simulations to create mock-observed data cubes, following the HSIM pipeline (Zieleniewski et al. 2015). An input scene was modeled with a central point-like young A6 star with a K = 6 magnitude, with a BT-Nextgen spectrum at 8000 K and a logarithm of the surface gravity of log g = 4.0 dex (Allard et al. 2012). The spectrum was convolved with a Gaussian line-spread function to match the resolution of the HARMONI data and was flux calibrated by the end-to-end throughput. The ESO Skycal sky model6 was used to compute the transmission of the terrestrial atmosphere at the observed airmass, further multiplied by the transmission of the ELT mirrors, the optical relay system, the instrument optics and the K-band grating, and the detector quantum efficiency7. The input scene was convolved with the coronagraphic PSFs for each wavelength and each exposure. Background noise was further added to the data cubes, including the sky (from the same ESO Skycal sky models), the telescope wa modeled as a gray body whereby the thermal emission at the site temperature (280.5 K) was multiplied with a constant emissivity7, the relay system and the instrument both emitting as pure blackbodies (260.5 and 130 K, respectively). Finally, Poisson noise from the observed star and the background was added along with the detector dark current (0.0053 e− s−1), also subject to photon noise, and the read-out noise (12 e−). Detector crosstalk (2%) with the four contiguous pixels to each pixel was also simulated in both spatial and spectral directions, as well as broadband diffusion (0.5%). Figure A.1 shows a simulated exposure at 2 μm, with the attenuated central star and the diffracted starlight, whose amplitude is minimized in the central region that is corrected by the adaptive optics.

Appendix B Distributions in the detection maps

To further illustrate the significance of the detection of the planets with spectral unmixing and cross-correlation, we built the distributions of the values in their respective maps. Figure B.2 shows the histograms corresponding to the maps for the single test case on the synthetic HARMONI data, as displayed in Fig. 3. Figure B.1 illustrated the case of the on-sky SINFONI data.

thumbnail Fig. B.1

Distributions of spatial weights (top) and cross-correlation strengths (bottom) after processing synthetic HARMONI data with a simulated planet at a contrast of 5 × 10−6. The values at the location of the planet are highlighted in green in each histogram. The peak is enhanced for the cross-correlation for clarity. The parameters of a Gaussian fit to the distribution of the cross-correlation strengths are given.
thumbnail Fig. B.2

Same as Fig. B.2, but for the SINFONI data of β Pictoris.

In the case of spectral unmixing, the spatial weights are positive by construction. A majority of the pixels in the field of view are zeroed because they are not decomposed by the planet spectrum. A significant number of pixels have low values as a result of common starlight residuals with the planet, as discussed in the plain text. A hard threshold around 0.3 marks the beginning of all pixels within the core of the planet, with weights up to 1. This holds for the synthetic and on-sky cases. Because of the shape of the distribution, the computed S/N cannot be converted into a false-alarm probability.

The cross-correlation histograms are best reproduced with a normal distribution centered on zero with very small width, although the distribution for the SINFONI shows large tails that arise from the bright structures in the map around the planet (see Fig. 7). In both cases, the peak value at the planet location clearly stands out from the remaining pixels.

Appendix C Spectral distances

With spectral unmixing, sample spectra are identified based on the dissimilarities between the spaxels within the field of view. Three metrics are commonly used in remote-sensing research: the orthogonal projection distance following Eq. (4), the cross-correlation, and the spectral angle mapper, which is defined as SAM(si,sj) = arccos(〈si,sj〉‖si‖ ‖sj‖) [rd].\[ \mathrm{SAM}(\bm{s}_{i},\bm{s}_{j})\,{=}\,\arccos{\Bigg(\frac{\langle\bm{s}_{i},\bm{s}_{j}\rangle}{\lVert\bm{s}_{i}\rVert\;\lVert\bm{s}_{j}\rVert}\Bigg)} \quad [\mathrm{rd}]. \]

The metric with the best spectral discriminatory power is the one that leads to the highest ratio between two similar and two distinct spectra (Chein-I Chang 2000).

In order to assess that the orthogonal projection distance is the most appropriate for our purpose, we computed this ratio for the three metrics. We used the hyperspectral simulated data cube from Sect. 3. sj was set as the spectrum corresponding to the spaxel at the peak of the planet PSF. Each metric was first computed with si being the spectrum of a spaxel within the core of the planet, targeting a very similar spectrum to sj. The metric was also computed 20 times with si as a random spectrum among all other spaxels, and then averaged to obtain a representative value for distinct spectra. The ratio of the two measurements was then calculated for all three metrics. We obtain 5.76 with the orthogonal projection distance, 4.82 with the cross-correlation, and 1.18 with SAM. Therefore, the orthogonal projection distance is 18% and nearly five times more effective than cross-correlation and SAM, respectively, to distinguish between any two spectral signatures relative to the planet spectrum.

References

  1. Abuter, R., Schreiber, J., Eisenhauer, F., et al. 2006, New A Rev., 50, 398 [Google Scholar]
  2. Acito, N., Diani, M., & Corsini, G. 2009, IEEE Trans. Geosci. Rem. Sens., 47, 3844 [Google Scholar]
  3. Allard, F. 2014, in Exploring the Formation and Evolution of Planetary Systems, eds. M. Booth, B. C. Matthews, & J. R. Graham, 299, 271 [Google Scholar]
  4. Allard, F., Homeier, D., & Freytag, B. 2012, Phil. Trans. Roy. Soc. Lond. A, 370, 2765 [Google Scholar]
  5. Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bagnasco, G., Kolm, M., Ferruit, P., et al. 2007, in Cryogenic Optical Systems and Instruments XII, eds. J. B. Heaney, & L. G. Burriesci, SPIE Conf. Ser., 6692, 66920M [Google Scholar]
  7. Barman, T. S., Macintosh, B., Konopacky, Q. M., & Marois, C. 2011, ApJ, 733, 65 [NASA ADS] [CrossRef] [Google Scholar]
  8. Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61 [NASA ADS] [CrossRef] [Google Scholar]
  9. Berné, O., Joblin, C., Deville, Y., et al. 2007, A&A, 469, 575 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bioucas-Dias, J., & Nascimento, J. 2008, IEEE Trans. Geosci. Rem. Sens., 46, 2435 [Google Scholar]
  12. Bioucas-Dias, J. M., Plaza, A., Dobigeon, N., et al. 2012, IEEE J. Sel. Top. Appl. Earth Observ. Rem. Sens., 5, 354 [Google Scholar]
  13. Bonnefoy, M., Lagrange, A. M., Boccaletti, A., et al. 2011, A&A, 528, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bonnefoy, M., Zurlo, A., Baudino, J. L., et al. 2016, A&A, 587, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Bonnefoy, M., Perraut, K., Lagrange, A. M., et al. 2018, A&A, 618, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Boulais, A., Berné, O., Faury, G., & Deville, Y. 2021, A&A 647, A105 [EDP Sciences] [Google Scholar]
  17. Brandt, T., Rizzo, M., Groff, T., et al. 2017, J. Astron. Telesc. Instrum. Syst., 3, 048002 [Google Scholar]
  18. Brandl, B. R., Absil, O., Agócs, T., et al. 2018, in Ground-based and Airborne Instrumentation for Astronomy VII, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 10702, 107021U [Google Scholar]
  19. Bro, R., & De Jong, S. 1997, J. Chemometrics, 11, 393 [Google Scholar]
  20. Bryan, M. L., Benneke, B., Knutson, H. A., Batygin, K., & Bowler, B. P. 2018, Nat. Astron., 2, 138 [Google Scholar]
  21. Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, A&A, 582, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Carlotti, A., Hénault, F., Dohlen, K., et al. 2018, in Ground-based and Airborne Instrumentation for Astronomy VII, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 10702, 107029N [Google Scholar]
  23. Chang, C.-I., & Du, Q. 2004, IEEE Trans. Geosci. Rem. Sens., 42, 608 [Google Scholar]
  24. Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, A&A, 605, A9 [EDP Sciences] [Google Scholar]
  25. Chauvin, G., Gratton, R., Bonnefoy, M., et al. 2018, A&A, 617, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Chein-I Chang. 2000, IEEE Transactions on Information Theory, 46, 1927 [Google Scholar]
  27. Chein-I Chang. 2005, IEEE Trans. Geosci. Rem. Sens., 43, 502 [Google Scholar]
  28. Chilcote, J., Pueyo, L., De Rosa, R. J., et al. 2017, AJ, 153, 182 [NASA ADS] [CrossRef] [Google Scholar]
  29. Currie, T., Brandt, T. D., Uyama, T., et al. 2018, AJ, 156, 291 [NASA ADS] [CrossRef] [Google Scholar]
  30. De Rosa, R. J., Rameau, J., Patience, J., et al. 2016, ApJ, 824, 121 [NASA ADS] [CrossRef] [Google Scholar]
  31. Delorme, P., Schmidt, T., Bonnefoy, M., et al. 2017, A&A, 608, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Deville, Y., Revel, C., Achard, V., & Briottet, X. 2014, IEEE Whispers [Google Scholar]
  33. Drumetz, L. 2016, Theses, Université Grenoble Alpes [Google Scholar]
  34. Drumetz, L., Chanussot, J., & Jutten, C. 2020, in Data Handling in Science and Technology, 32, Hyperspectral Imaging, ed. J. M. Amigo (Elsevier), 167, 203 [Google Scholar]
  35. Esposito, S., Riccardi, A., Pinna, E., et al. 2011, in Astronomical Adaptive Optics Systems and Applications IV, eds. R. K. Tyson, & M. Hart, SPIE Conf. Ser., 8149, 814902 [Google Scholar]
  36. Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, A&A, 618, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Forni, O., Poulet, F., Bibring, J. P., et al. 2005, in 36th Annual Lunar and Planetary Science Conference, eds. S. Mackwell, & E. Stansbery, 1623 [Google Scholar]
  38. Foschino, S., Berné, O., & Joblin, C. 2019, A&A, 632, A84 [EDP Sciences] [Google Scholar]
  39. Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Gomez Gonzalez, C. A., Absil, O., Absil, P. A., et al. 2016, A&A, 589, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Gratier, P., Bron, E., Gerin, M., et al. 2017, A&A, 599, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Greenbaum, A. Z., Pueyo, L., Ruffio, J.-B., et al. 2018, AJ, 155, 226 [Google Scholar]
  43. Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
  44. Halimi, A., Honeine, P., Kharouf, M., Richard, C., & Tourneret, J.-Y. 2015, IEEE Trans. Geosci. Rem. Sens., 54 [Google Scholar]
  45. Harsanyi, J. C., & Chang, C. 1994, IEEE Trans. Geosci. Rem. Sens., 32, 779 [Google Scholar]
  46. Hauksdottir, H., Jutten, C., Schmidt, F., et al. 2006, in 7th Nordic Signal Processing Symposium (NORSIG’2006) [Google Scholar]
  47. Hoeijmakers, H. J., Schwarz, H., Snellen, I. A. G., et al. 2018, A&A, 617, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Jensen-Clem, R., Mawet, D., Gomez Gonzalez, C. A., et al. 2018, AJ, 155, 19 [Google Scholar]
  49. Jovanovic, N., Martinache, F., Guyon, O., et al. 2015, PASP, 127, 890 [NASA ADS] [CrossRef] [Google Scholar]
  50. Juvela, M., Lehtinen, K., & Paatero, P. 1996, MNRAS, 280, 616 [NASA ADS] [Google Scholar]
  51. Kawahara, H., Murakami, N., Matsuo, T., & Kotani, T. 2014, ApJS, 212, 27 [NASA ADS] [CrossRef] [Google Scholar]
  52. Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Science, 339, 1398 [Google Scholar]
  54. Kruse, F., Lefkoff, A., Boardman, J., et al. 1993, Rem. Sens. Environ., 44, 145 [Google Scholar]
  55. Kuntschner, H., Jochum, L., Amico, P., et al. 2014, in Ground-based and Airborne Instrumentation for Astronomy V, eds. S. K. Ramsay, I. S. McLean, & H. Takami, SPIE Conf. Ser., 9147, 91471U [Google Scholar]
  56. Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660, 770 [NASA ADS] [CrossRef] [Google Scholar]
  57. Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57 [Google Scholar]
  58. Lagrange, A. M., Boccaletti, A., Langlois, M., et al. 2019, A&A, 621, A8 [EDP Sciences] [Google Scholar]
  59. Larkin, J. E., Moore, A. M., Wright, S. A., et al. 2016, in Ground-based and Airborne Instrumentation for Astronomy VI, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 9908, 99081W [Google Scholar]
  60. Li, J., Zhang, H., Zhang, L., & Ma, L. 2015, IEEE J. Sel. Top. Appl. Earth Observ. Rem. Sens., 8, 2523 [Google Scholar]
  61. Liu, J., Luo, B., Douté, S., & Chanussot, J. 2018, Rem. Sens., 10, 737 [Google Scholar]
  62. Lovis, C., Snellen, I., Mouillet, D., et al. 2017, A&A, 599, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Luo, B., Chanussot, J., Douté, S., & Zhang, L. 2013, IEEE Geosci. Rem. Sens. Lett., 10, 24 [Google Scholar]
  64. Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proc. Natl. Acad. Sci. U.S.A., 111, 12661 [Google Scholar]
  65. Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64 [NASA ADS] [CrossRef] [Google Scholar]
  66. Marois, C., Doyon, R., Nadeau, D., et al. 2003, in High-Contrast Imaging for Exo-Planet Detection, ed. A. B. Schultz, SPIE Conf. Ser., 4860, 130 [Google Scholar]
  67. Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
  68. Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97 [NASA ADS] [CrossRef] [Google Scholar]
  69. Mawet, D., Bond, C. Z., Delorme, J. R., et al. 2018, in Adaptive Optics Systems VI, eds. L. M. Close, L. Schreiber, & D. Schmidt, SPIE Conf. Ser., 10703, 1070306 [Google Scholar]
  70. McGregor, P. J., Bloxham, G. J., Boz, R., et al. 2012, in Ground-based and Airborne Instrumentation for Astronomy IV, eds. I. S. McLean, S. K. Ramsay, & H. Takami, SPIE Conf. Ser., 8446, 84461I [Google Scholar]
  71. Mesa, D., Gratton, R., Zurlo, A., et al. 2015, A&A, 576, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Moussaoui, S., Hauksdottir, H., Schmidt, F., et al. 2008, Neurocomputing, 71, 2194 [Google Scholar]
  73. Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Nielsen, E., De Rosa, R., Macintosh, B., et al. 2019, in AAS/Division for Extreme Solar Systems Abstracts, 51, 100.02 [Google Scholar]
  75. Nowak, M., Lacour, S., Molliére, P., et al. 2020, A&A, 633, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Petit dit de la Roche, D. J. M., Hoeijmakers, H. J., & Snellen, I. A. G. 2018, A&A, 616, A146 [EDP Sciences] [Google Scholar]
  77. Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2021, A&A, 648, A59 [EDP Sciences] [Google Scholar]
  78. Plaza, A., Martinez, P., Perez, R., & Plaza, J. 2004, IEEE Trans. Geosci. Rem. Sens., 42, 650 [Google Scholar]
  79. Racine, R., Walker, G. A. H., Nadeau, D., Doyon, R., & Marois, C. 1999, PASP, 111, 587 [Google Scholar]
  80. Rajan, A., Rameau, J., De Rosa, R. J., et al. 2017, AJ, 154, 10 [NASA ADS] [CrossRef] [Google Scholar]
  81. Rapacioli, M., Joblin, C., & Boissel, P. 2005, A&A, 429, 193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, ApJ, 852, 104 [CrossRef] [Google Scholar]
  83. Riaud, P., & Schneider, J. 2007, A&A, 469, 355 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Rieke, G. H., Wright, G. S., Buker, T., et al. 2015, PASP, 127, 584 [NASA ADS] [CrossRef] [Google Scholar]
  85. Ruffio, J.-B., Macintosh, B., Wang, J. J., et al. 2017, ApJ, 842, 14 [NASA ADS] [CrossRef] [Google Scholar]
  86. Ruffio, J.-B., Macintosh, B., Konopacky, Q. M., et al. 2019, AJ, 158, 200 [Google Scholar]
  87. Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Samland, M., Bouwman, J., Hogg, D. W., et al. 2021, A&A 646, A24 [CrossRef] [EDP Sciences] [Google Scholar]
  89. Schwartz, J. C., Sekowski, C., Haggard, H. M., Pallé, E., & Cowan, N. B. 2016, MNRAS, 457, 926 [Google Scholar]
  90. Skemer, A. J., Hinz, P., Montoya, M., et al. 2015, in Techniques and Instrumentation for Detection of Exoplanets VII, ed. S. Shaklan, SPIE Conf. Ser., 9605, 96051D [Google Scholar]
  91. Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Nature, 509, 63 [NASA ADS] [CrossRef] [Google Scholar]
  92. Snellen, I., de Kok, R., Birkby, J. L., et al. 2015, A&A, 576, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Soofbaf, S. R., Sahebi, M., & Mojaradi, B. 2018, Rem. Sens., 10, 434 [Google Scholar]
  94. Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28 [NASA ADS] [CrossRef] [Google Scholar]
  95. Sparks, W. B., & Ford, H. C. 2002, ApJ, 578, 543 [NASA ADS] [CrossRef] [Google Scholar]
  96. Stone, J. M., Barman, T., Skemer, A. J., et al. 2020, AJ, 160, 262 [Google Scholar]
  97. Thatte, N. A., Clarke, F., Bryson, I., et al. 2016, in Ground-based and Airborne Instrumentation for Astronomy VI, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 9908, 99081X [Google Scholar]
  98. Themelis, K. E., Schmidt, F., Sykioti, O., et al. 2012, Planet. Space Sci., 68, 34 [Google Scholar]
  99. Tsinos, C. G., Rontogiannis, A. A., & Berberidis, K. 2017, IEEE Trans. Comput. Imaging, 3, 160 [Google Scholar]
  100. Uyama, T., Currie, T., Hori, Y., et al. 2020, AJ, 159, 40 [Google Scholar]
  101. Vigan, A., Otten, G. P. P. L., Muslimov, E., et al. 2018, in Ground-based and Airborne Instrumentation for Astronomy VII, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 10702, 1070236 [Google Scholar]
  102. Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, in press, https://doi.org/10.1051/0004-6361/202038107 [Google Scholar]
  103. Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015, pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophys. Source Code Libr., [record ascl:1506.001] [Google Scholar]
  104. Wang, J., Mawet, D., Ruane, G., Hu, R., & Benneke, B. 2017, AJ, 153, 183 [NASA ADS] [CrossRef] [Google Scholar]
  105. Wang, J., Mawet, D., Fortney, J. J., et al. 2018, AJ, 156, 272 [NASA ADS] [CrossRef] [Google Scholar]
  106. Ward-Duong, K., Patience, J., Follette, K., et al. 2021, AJ, 161, 5 [Google Scholar]
  107. Wilcomb, K. K., Konopacky, Q. M., Barman, T. S., et al. 2020, AJ, 160, 207 [Google Scholar]
  108. Xie, C., Haffert, S. Y., de Boer, J., et al. 2020, A&A 644, A149 [EDP Sciences] [Google Scholar]
  109. Xu, Y., Wu, Z., Chanussot, J., et al. 2018, IEEE Trans. Geosci. Rem. Sens., 56, 1680 [Google Scholar]
  110. Yang, Y., Zhang, J., Song, S., & Liu, D. 2019, Rem. Sens., 11, 192 [Google Scholar]
  111. Zhang, X., Li, C., Zhang, J., et al. 2018, Rem. Sens., 10, 339 [Google Scholar]
  112. Zieleniewski, S., Thatte, N., Kendrew, S., et al. 2015, MNRAS, 453, 3754 [CrossRef] [Google Scholar]
1

“Frequency” is used here and throughout to describe the shape of spectra at different spectral scalings, and it should not be confused with the inverse of the unit wavelength.

2

Because of these variations, nonplanetary sample spectra should be interpreted with caution and are not to be considered as perfect stellar spectra.

3

〈s0〉⊥$\langle\bm{s}_{0}\rangle^{\bot}$ is defined such that ∀v∈〈s0〉⊥$\forall\bm{v}\in \langle\bm{s}_{0}\rangle^{\bot}$, v ⋅s0 = 0.

4

The sample spectra do not form an orthogonal basis, ⊥ denotes the orthogonality of the space to the origin vectors.

7

All curves and constants are given in the HSIM source code https://github.com/HARMONI-ELT/HSIM

All Figures

thumbnail Fig. 1

Model spectra of a β Pictoris-like system: an A6 star with Teff = 8000 K (BT-Nextgen Allard et al. 2012) and a giant planet with Teff = 1700 K (BT-Settl Allard 2014), normalized over the K bandpass at a resolution of R = 7000. Not only the slope of the continuum is reversed between the two objects, but also the absorption line content, with the prominent Brackett-γ line in the stellar spectrum at 2.16μm, while the planet spectrum exhibits a forest of H2O and CO lines.
In the text
thumbnail Fig. 2

Top: spectra of individual spaxels in synthetic coronagraphic HARMONI data at R ~ 7000 with a central A6 star and a 1700 K planet (orange) (see Fig. 1). Random spaxels (gray) in the corrected region illustrate the modulation by the speckles of the stellar spectrum, as shown with the central spaxel behind the semitransparent focal plane mask (blue). Details of the simulations are given in Sect. 3. Bottom: resulting spectra after the continuum was removed (following Hoeijmakers et al. 2018). The amplitude of the planet spectra has been amplified for visualization purposes.
In the text
thumbnail Fig. 3

Final result of post-processing the synthetic HARMONI data with a cross-correlation with the planet spectral template (left) and spectral unmixing through orthogonal subspace projection (middle and right). The number of spectral signatures is set to eight to decompose the data. The injected mock planet has a broadband contrast of 5 × 10−6 with the central star at a separation of 100 mas (white circle). Images are normalized to the peak value of the planet for each algorithm, the scale is linear, the field of view is identical, and the white cross marks the position of the star.
In the text
thumbnail Fig. 4

Inferred spectrum of the planet with a contrast 5 × 10−6 injected into the synthetic HARMONI data (see Fig. 3 with spectral unmixing through orthogonal subspace projection (dark blue), compared with the high-frequency component of the injected spectrum (orange) along with the difference between the two (black). The first-overtone bandheads of CO are highlighted as being well recovered (shade). The unmixing is not perfect, and the spectrum of the planet is contaminated with residuals from the subtraction of the diffracted starlight (bottom).
In the text
thumbnail Fig. 5

Completeness curves (95%) with spectral unmixing (dark blue) and cross-correlation (brown) expressed as the minimum planet-to-star contrast for a detection threshold of τ = 5σ as a function of the projected separation for L-type (solid) and T-type spectra (dot-dashed) in the synthetic HARMONI data.
In the text
thumbnail Fig. 6

Spatial weight inversion map of synthetic HARMONI data with a spiral of simulated planets (white circles) with a star-to-planet flux ratio of 9 × 10−7 (top) and corresponding S/N map (bottom). The injected model spectrum is used to invert the data following Eq. (5). Despite tiny spatial weights, most planets are detected with a high S/N.
In the text
thumbnail Fig. 7

Final results of post-processing VLT/SINFONI K-band data of β Pictoris b with the cross-correlation (top) and spectral unmixing through orthogonal subspace projection (bottom). A BT-Settl model spectrum (1700 K, log g = 4.0) was used for the cross-correlation, and p = 10 was set to decompose the data. The scale is linear and was chosen to reach 90% of the peak value for each image, the cross-correlation strength Cs, and the weight αi, respectively. The cross marks the position of β Pictoris A.
In the text
thumbnail Fig. 8

Top: extracted continuum-removed SINFONI spectrum of β Pictoris b with spectral unmixing (dark blue), compared with a BT-Settl model spectrum at 1700 K and log g = 4.0 dex (orange), close to the true planet properties. The first-overtone bandheads of CO are highlighted as being well recovered (shade). The extracted spectrum is not perfect, it is contaminated by unmixed residuals. Bottom: histogram of the residuals (black) with a best-fit Gaussian function (purple).
In the text
thumbnail Fig. A.1

Simulated coronagraphic ELT/HARMONI image at 2.0 μm for an exposure of 60 s with the SP1 apodizer. The central star is attenuated by the focal plane mask, which is asymmetric to take the atmospheric dispersion into account. The remaining field of view is plagued by diffracted starlight. The white line marks the adaptive optics corrected area, the region where the speckle intensity is minimized. The intensity scale is linear.
In the text
thumbnail Fig. B.1

Distributions of spatial weights (top) and cross-correlation strengths (bottom) after processing synthetic HARMONI data with a simulated planet at a contrast of 5 × 10−6. The values at the location of the planet are highlighted in green in each histogram. The peak is enhanced for the cross-correlation for clarity. The parameters of a Gaussian fit to the distribution of the cross-correlation strengths are given.
In the text
thumbnail Fig. B.2

Same as Fig. B.2, but for the SINFONI data of β Pictoris.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.