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/00046361/202140337  
Published online  28 May 2021 
Spectral unmixing for exoplanet direct detection in hyperspectral data
^{1}
Univ. Grenoble Alpes, CNRS, IPAG,
38000
Grenoble, France
email: julien.rameau@univgrenoblealpes.fr
^{2}
Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, GIPSALab,
38000
Grenoble, France
Received:
13
January
2021
Accepted:
18
March
2021
Context. The direct detection of faint exoplanets with highcontrast instruments can be boosted by combining it with high spectral resolution. For integral field spectrographs yielding hyperspectral data, this means that the majority of the field of view consists of diffracted starlight spectra and a spatially localized planet. Observation analysis usually relies on classic crosscorrelation with theoretical spectra, maximized at the position and with the properties of the planet. In a purely blindsearch context, this supervised strategy can be biased with model mismatch and/or be computationally inefficient.
Aims. Using an approach that is inspired by the analysis of hyperspectral data within the remotesensing community, we aim to propose an alternative to crosscorrelation that is fully datadriven, which decomposes the data into a set of individual spectra and their corresponding spatial distributions. This strategy is called spectral unmixing.
Methods. We used an orthogonal subspace projection to identify the most distinct spectra in the field of view. Their spatial distribution maps were then obtained by inverting the data. These spectra were then used to break the original hyperspectral images into their corresponding spatial distribution maps via nonnegative least squares. A matched filter with the instrument pointspread function (or visual inspection) was then used to detect the planet on one of the maps. The performance of our method was evaluated and compared with a crosscorrelation using simulated hyperspectral data with medium resolution from the ELT/HARMONI integral field spectrograph.
Results. We show that spectral unmixing effectively leads to a planet detection solely based on spectral dissimilarities at significantly reduced computational cost. The extracted spectrum holds significant signatures of the planet while being not perfectly separated from residual starlight. The sensitivity of the supervised crosscorrelation is three to four times higher than with unsupervised spectral unmixing, the gap is biased toward the former because the injected and correlated spectrum match perfectly. The algorithm was furthermore vetted on real data obtained with VLT/SINFONI of the β Pictoris system. This led to the detection of β Pictoris b with a signaltonoise ratio of 28.5.
Conclusions. Spectral unmixing is a viable alternative strategy to a crosscorrelation to search for and characterize exoplanets in hyperspectral data in a purely datadriven approach. The advent of large data from the forthcoming IFS on board JWST and future ELTs motivates further algorithm development along this path.
Key words: methods: data analysis / techniques: imaging spectroscopy / planets and satellites: detection
© J. Rameau et al. 2021
Open 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 pointlike 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 closein exoplanets is nevertheless plagued by residual starlight stemming from uncalibrated instrumental aberrations, and for groundbased facilities, Earth’s atmospheric turbulence. Disentangling the exoplanet from these socalled 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 postprocessing 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 groundbased extreme adaptiveoptics systems Gemini/GPI (Macintosh et al. 2014), the Very Large Telescopes VLT/SPHERE (Beuzit et al. 2019), Subaru/SCExAOCHARIS (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 pseudocontinuum 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; WardDuong 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 crosscorrelations 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 mediumresolution IFS Keck/OSIRIS and VLT/Spectrograph for INtegral Field Observations in the NearInfrared 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 10mclass 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 crosscorrelation 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 crosscorrelated 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 sourcedetection 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 crosscorrelation 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 ELTsIFS such as the High Angular Resolution Monolithic Optical and NearInfrared 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 sourcedetection in hyperspectral data, or spectral unmixing for this type of data, represents a longstanding problem beyond astronomy, notably in Earth remotesensing and microscopy. The dedicated algorithms aim at separating various sources in the field of view from their distinct spectral morphologies in a purely datadriven 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 BioucasDias 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 crosscorrelation technique to directly detect an exoplanet and extract its spectrum from mediumresolution data. We test one algorithm on synthetic R ~ 7000 Kband HARMONI data to evaluate its performances. We then apply it on onsky R ~ 5000 Kband 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 crosscorrelation 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 T_{eff}, 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 metalbearing ions, atoms, and molecules for stars, and H_{2}O, carbon andnitrogenbearing 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 frequency^{1} 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 d_{i} in the hyperspectral data can be modeled with a linear statement, such that (1)
where s_{j} 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 (2)
where is the columnvectorized hyperspectral images where the columns contain the m = n_{x} × n_{y} spaxels in any order of n_{λ} spectral channels, is a matrix of p source spectra, 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 socalled endmembers, and A contains theabundance maps (BioucasDias 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.
Fig. 1 Model spectra of a β Pictorislike system: an A6 star with T_{eff} = 8000 K (BTNextgen Allard et al. 2012) and a giant planet with T_{eff} = 1700 K (BTSettl 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 H_{2}O 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 highcontrast 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 highdimension hyperspectral data usually lies in a lowerdimensional subspace (e.g., Chang & Du 2004; BioucasDias & 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 lowcontrast 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 = n_{x} * n_{y}, 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.
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 remotesensing research; many strategies have been investigated based on geometrical or statistical considerations (see the reviews by Plaza et al. 2004; BioucasDias et al. 2012, and references therein). We considered a geometrybased 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 highpass 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 highfrequency components.
Following Eq. (2) and with p ≪ n_{λ}, all spaxels live in a lowdimensional 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; CheinI Chang 2005).
Let s_{0} be the first sample spectrum. The sparsity of the planet over the diffracted starlight in the field of view motivates the choice of s_{0} as the spectrum of the median spaxel. Let U be the spectral sample matrix, which at this step only contains s_{0}, such that U = [s_{0}]. The orthogonal projector from U, denoted , is defined as (3)
with I a n_{λ} × n_{λ} identity matrix, U^{T} the transpose of U, and the pseudoinverse of U. The dot product projects each spaxel d_{i} in , the subspace orthogonal to the space linearly spanned by s_{0} ^{3}. Iterating over all spaxels, a second sample spectrum, set as s_{1}, is flagged as the most distant spaxel to s_{0}. The projection distance is chosen to be the sum of the squared residuals after removing from d_{i} the projection of d_{i} onto s_{0}, or (4)
A new orthogonal subspace projector can be defined from the linear subspace spanned with U = [s_{0}, s_{1}] and applied to the remaining spaxels. The third sample spectrum is identified as the one with the maximum absolute orthogonal projection in the subspace ⟨s_{0},s_{1}⟩^{⊥}^{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 leastsquares 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 nonnegative leastsquares estimate of A is done for each spaxel d_{i} by finding
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 nonnegative leastsquares problem is accelerated by replacing S with a p × p matrix following (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 pointspread 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 highfrequency components;
 3.
Vectorize the hyperspectral data D in a 1D spatial plus 1D spectral matrix;
 4.
Set the first sample s_{0} as the closest spectrum to the median spaxel and define U = [s_{0}] the subspace spanned by s_{0};
 5.
Project each spaxel onto the linear subspace spanned by U with the orthogonal projector defined in Eq. (3);
 6.
Find thespaxel with the extremal projection (see Eq. (4)), which corresponds to a new sample spectrum, s_{new};
 7.
Append U with s_{new};
 8.
Repeat steps 5–7 until all p spectra are extracted;
 9.
Invert D with fast nonnegative 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 highcontrast hyperspectral HARMONI data. Then we compare it with classic crosscorrelation.
3.1 Simulations of ELT/HARMONI data
HARMONI will be the firstlight IFS on the ELT offering diffractionlimited hyperspectral images at medium spectral resolution (3500–18 000) in the optical and nearinfrared (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 Nyquistsampled 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 endtoend simulator with the following configuration: 4 h of observations with oneminutelong 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 offaxis simulated PSF is used as a template for a simulated planet, and a BTSettl 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 Ltype 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 spectrumbased detection algorithm, each spaxel in the final 3D data cube was crosscorrelated with the injected model spectrum, whose pseudocontinuum was subtracted beforehand with a convolution with a Gaussian linespread function corresponding to a spectral resolution of 100. The crosscorrelation used a customized version of CrosscorrRV^{5} 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 crosscorrelation 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 offaxis PSF. The S/N with the crosscorrelation is 23.6, while the S/N with the spectral unmixing is 12.6. The distribution of the values in the crosscorrelation 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 crosscorrelation with the injected model is 0.15, which is close to the amplitude (0.17) at the planet location in the crosscorrelation map, as expected. Crosscorrelating 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 crosscorrelation 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 crosscorrelation 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 priorfree alternative to classic crosscorrelation.
Fig. 3 Final result of postprocessing the synthetic HARMONI data with a crosscorrelation 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. 
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 highfrequency component of the injected spectrum (orange) along with the difference between the two (black). The firstovertone 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 crosscorrelation beyond a single test case, we computed the detection sensitivity in the central region, which is expressed in the startoplanet 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 (JensenClem et al. 2018). The process was repeated with spectral unmixing with p = 10 and a crosscorrelation 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 planetlike 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 crosscorrelation 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 crosscorrelated 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 crosscorrelation of all other 10− 1 spectra with the planet template, which is stable around 0.027 ± 0.005.
Fig. 5 Completeness curves (95%) with spectral unmixing (dark blue) and crosscorrelation (brown) expressed as the minimum planettostar contrast for a detection threshold of τ = 5σ as a function of the projected separation for Ltype (solid) and Ttype spectra (dotdashed) in the synthetic HARMONI data. 
Fig. 6 Spatial weight inversion map of synthetic HARMONI data with a spiral of simulated planets (white circles) with a startoplanet 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 Ltype 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 crosscorrelation 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 Ttype 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 crosscorrelation.
The resulting curves are shown in Fig. 5. Spectral unmixing can detect fainter Ttype planets than Ltype planets, the limiting contrast being improved by 55% down to 1.8 × 10^{−6}. The sensitivity is also boosted with the crosscorrelation, as predicted. Ttype planets with contrast as faint as 5 × 10^{−7} can be detected,which is an improvement of 45% compared with Ltype planets. A switch from Ltype to Ttype 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 crosscorrelation, 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 postprocessing 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 highcontrast instrument, it delivers mediumresolution 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 mas^{2} 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. Twentyfour cubes with an integration time of 60s were obtained. Data cubes were initially built from the 2D raw detector images with the SINFONI datahandling 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 spatialspectral 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 crosscorrelation and spectral unmixing. The crosscorrelation 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 continuumremoved BTSettl 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 pseudoS/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 pseudoS/N of 28.5 with spectral unmixing, and this value is 8.2 with the crosscorrelation.
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 BTSettl 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 crosscorrelation of the model spectrum with the extracted planet spectrum, which yields 0.33 (the crosscorrelation 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.
Fig. 7 Final results of postprocessing VLT/SINFONI Kband data of β Pictoris b with the crosscorrelation (top) and spectral unmixing through orthogonal subspace projection (bottom). A BTSettl model spectrum (1700 K, log g = 4.0) was used for the crosscorrelation, 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 crosscorrelation strength C_{s}, 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 crosscorrelation 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 crosscorrelation 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. Crosscorrelation 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 geometrybased algorithm such as we adopted, is not robust enough (BioucasDias 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 postprocessing approach.
Beyond these limitations, the spectral unmixing algorithm presented in this work has two advantages. It is fully datadriven, or unsupervised, and computationally fast. Both characteristics make it interesting with respect to the other strategy for analyzing hyperspectral data, the classic crosscorrelation. For a purely blind search of real data, the crosscorrelation 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 Neptunelike 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 crosscorrelation 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 dualprocessor core running at 2.4 GHz.
Fig. 8 Top: extracted continuumremoved SINFONI spectrum of β Pictoris b with spectral unmixing (dark blue), compared with a BTSettl model spectrum at 1700 K and log g = 4.0 dex (orange), close to the true planet properties. The firstovertone 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 bestfit Gaussian function (purple). 
6.2 Summary and future work
To conclude, we have introduce a new strategy for analyzing highcontrast 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 remotesensing 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 nonnegative leastsquares 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 datadriven as opposed to a crosscorrelation 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 crosscorrelation 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 realtime 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. Lowrank 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 15IDEX02), through the funding of the “Origin of Life” project of the Univ. GrenobleAlpes. 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 highcontrast 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 offaxis PSF is also generated.
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 mockobserved data cubes, following the HSIM pipeline (Zieleniewski et al. 2015). An input scene was modeled with a central pointlike young A6 star with a K = 6 magnitude, with a BTNextgen 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 linespread function to match the resolution of the HARMONI data and was flux calibrated by the endtoend throughput. The ESO Skycal sky model^{6} 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 Kband grating, and the detector quantum efficiency^{7}. 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 emissivity^{7}, 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 readout 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 crosscorrelation, 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 onsky SINFONI data.
Fig. B.1 Distributions of spatial weights (top) and crosscorrelation 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 crosscorrelation for clarity. The parameters of a Gaussian fit to the distribution of the crosscorrelation strengths are given. 
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 onsky cases. Because of the shape of the distribution, the computed S/N cannot be converted into a falsealarm probability.
The crosscorrelation 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 remotesensing research: the orthogonal projection distance following Eq. (4), the crosscorrelation, and the spectral angle mapper, which is defined as
The metric with the best spectral discriminatory power is the one that leads to the highest ratio between two similar and two distinct spectra (CheinI 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. s_{j} was set as the spectrum corresponding to the spaxel at the peak of the planet PSF. Each metric was first computed with s_{i} being the spectrum of a spaxel within the core of the planet, targeting a very similar spectrum to s_{j}. The metric was also computed 20 times with s_{i} 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 crosscorrelation, and 1.18 with SAM. Therefore, the orthogonal projection distance is 18% and nearly five times more effective than crosscorrelation and SAM, respectively, to distinguish between any two spectral signatures relative to the planet spectrum.
References
 Abuter, R., Schreiber, J., Eisenhauer, F., et al. 2006, New A Rev., 50, 398 [Google Scholar]
 Acito, N., Diani, M., & Corsini, G. 2009, IEEE Trans. Geosci. Rem. Sens., 47, 3844 [Google Scholar]
 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]
 Allard, F., Homeier, D., & Freytag, B. 2012, Phil. Trans. Roy. Soc. Lond. A, 370, 2765 [Google Scholar]
 Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Barman, T. S., Macintosh, B., Konopacky, Q. M., & Marois, C. 2011, ApJ, 733, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Berné, O., Joblin, C., Deville, Y., et al. 2007, A&A, 469, 575 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BioucasDias, J., & Nascimento, J. 2008, IEEE Trans. Geosci. Rem. Sens., 46, 2435 [Google Scholar]
 BioucasDias, J. M., Plaza, A., Dobigeon, N., et al. 2012, IEEE J. Sel. Top. Appl. Earth Observ. Rem. Sens., 5, 354 [Google Scholar]
 Bonnefoy, M., Lagrange, A. M., Boccaletti, A., et al. 2011, A&A, 528, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonnefoy, M., Zurlo, A., Baudino, J. L., et al. 2016, A&A, 587, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonnefoy, M., Perraut, K., Lagrange, A. M., et al. 2018, A&A, 618, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boulais, A., Berné, O., Faury, G., & Deville, Y. 2021, A&A 647, A105 [EDP Sciences] [Google Scholar]
 Brandt, T., Rizzo, M., Groff, T., et al. 2017, J. Astron. Telesc. Instrum. Syst., 3, 048002 [Google Scholar]
 Brandl, B. R., Absil, O., Agócs, T., et al. 2018, in Groundbased and Airborne Instrumentation for Astronomy VII, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 10702, 107021U [Google Scholar]
 Bro, R., & De Jong, S. 1997, J. Chemometrics, 11, 393 [Google Scholar]
 Bryan, M. L., Benneke, B., Knutson, H. A., Batygin, K., & Bowler, B. P. 2018, Nat. Astron., 2, 138 [Google Scholar]
 Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, A&A, 582, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carlotti, A., Hénault, F., Dohlen, K., et al. 2018, in Groundbased and Airborne Instrumentation for Astronomy VII, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 10702, 107029N [Google Scholar]
 Chang, C.I., & Du, Q. 2004, IEEE Trans. Geosci. Rem. Sens., 42, 608 [Google Scholar]
 Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, A&A, 605, A9 [EDP Sciences] [Google Scholar]
 Chauvin, G., Gratton, R., Bonnefoy, M., et al. 2018, A&A, 617, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 CheinI Chang. 2000, IEEE Transactions on Information Theory, 46, 1927 [Google Scholar]
 CheinI Chang. 2005, IEEE Trans. Geosci. Rem. Sens., 43, 502 [Google Scholar]
 Chilcote, J., Pueyo, L., De Rosa, R. J., et al. 2017, AJ, 153, 182 [NASA ADS] [CrossRef] [Google Scholar]
 Currie, T., Brandt, T. D., Uyama, T., et al. 2018, AJ, 156, 291 [NASA ADS] [CrossRef] [Google Scholar]
 De Rosa, R. J., Rameau, J., Patience, J., et al. 2016, ApJ, 824, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Delorme, P., Schmidt, T., Bonnefoy, M., et al. 2017, A&A, 608, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deville, Y., Revel, C., Achard, V., & Briottet, X. 2014, IEEE Whispers [Google Scholar]
 Drumetz, L. 2016, Theses, Université Grenoble Alpes [Google Scholar]
 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]
 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]
 Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, A&A, 618, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 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]
 Foschino, S., Berné, O., & Joblin, C. 2019, A&A, 632, A84 [EDP Sciences] [Google Scholar]
 Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gomez Gonzalez, C. A., Absil, O., Absil, P. A., et al. 2016, A&A, 589, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gratier, P., Bron, E., Gerin, M., et al. 2017, A&A, 599, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Greenbaum, A. Z., Pueyo, L., Ruffio, J.B., et al. 2018, AJ, 155, 226 [Google Scholar]
 Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
 Halimi, A., Honeine, P., Kharouf, M., Richard, C., & Tourneret, J.Y. 2015, IEEE Trans. Geosci. Rem. Sens., 54 [Google Scholar]
 Harsanyi, J. C., & Chang, C. 1994, IEEE Trans. Geosci. Rem. Sens., 32, 779 [Google Scholar]
 Hauksdottir, H., Jutten, C., Schmidt, F., et al. 2006, in 7th Nordic Signal Processing Symposium (NORSIG’2006) [Google Scholar]
 Hoeijmakers, H. J., Schwarz, H., Snellen, I. A. G., et al. 2018, A&A, 617, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 JensenClem, R., Mawet, D., Gomez Gonzalez, C. A., et al. 2018, AJ, 155, 19 [Google Scholar]
 Jovanovic, N., Martinache, F., Guyon, O., et al. 2015, PASP, 127, 890 [NASA ADS] [CrossRef] [Google Scholar]
 Juvela, M., Lehtinen, K., & Paatero, P. 1996, MNRAS, 280, 616 [NASA ADS] [Google Scholar]
 Kawahara, H., Murakami, N., Matsuo, T., & Kotani, T. 2014, ApJS, 212, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Science, 339, 1398 [Google Scholar]
 Kruse, F., Lefkoff, A., Boardman, J., et al. 1993, Rem. Sens. Environ., 44, 145 [Google Scholar]
 Kuntschner, H., Jochum, L., Amico, P., et al. 2014, in Groundbased and Airborne Instrumentation for Astronomy V, eds. S. K. Ramsay, I. S. McLean, & H. Takami, SPIE Conf. Ser., 9147, 91471U [Google Scholar]
 Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57 [Google Scholar]
 Lagrange, A. M., Boccaletti, A., Langlois, M., et al. 2019, A&A, 621, A8 [EDP Sciences] [Google Scholar]
 Larkin, J. E., Moore, A. M., Wright, S. A., et al. 2016, in Groundbased and Airborne Instrumentation for Astronomy VI, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 9908, 99081W [Google Scholar]
 Li, J., Zhang, H., Zhang, L., & Ma, L. 2015, IEEE J. Sel. Top. Appl. Earth Observ. Rem. Sens., 8, 2523 [Google Scholar]
 Liu, J., Luo, B., Douté, S., & Chanussot, J. 2018, Rem. Sens., 10, 737 [Google Scholar]
 Lovis, C., Snellen, I., Mouillet, D., et al. 2017, A&A, 599, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Luo, B., Chanussot, J., Douté, S., & Zhang, L. 2013, IEEE Geosci. Rem. Sens. Lett., 10, 24 [Google Scholar]
 Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proc. Natl. Acad. Sci. U.S.A., 111, 12661 [Google Scholar]
 Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Marois, C., Doyon, R., Nadeau, D., et al. 2003, in HighContrast Imaging for ExoPlanet Detection, ed. A. B. Schultz, SPIE Conf. Ser., 4860, 130 [Google Scholar]
 Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
 Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 McGregor, P. J., Bloxham, G. J., Boz, R., et al. 2012, in Groundbased and Airborne Instrumentation for Astronomy IV, eds. I. S. McLean, S. K. Ramsay, & H. Takami, SPIE Conf. Ser., 8446, 84461I [Google Scholar]
 Mesa, D., Gratton, R., Zurlo, A., et al. 2015, A&A, 576, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moussaoui, S., Hauksdottir, H., Schmidt, F., et al. 2008, Neurocomputing, 71, 2194 [Google Scholar]
 Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nielsen, E., De Rosa, R., Macintosh, B., et al. 2019, in AAS/Division for Extreme Solar Systems Abstracts, 51, 100.02 [Google Scholar]
 Nowak, M., Lacour, S., Molliére, P., et al. 2020, A&A, 633, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petit dit de la Roche, D. J. M., Hoeijmakers, H. J., & Snellen, I. A. G. 2018, A&A, 616, A146 [EDP Sciences] [Google Scholar]
 Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2021, A&A, 648, A59 [EDP Sciences] [Google Scholar]
 Plaza, A., Martinez, P., Perez, R., & Plaza, J. 2004, IEEE Trans. Geosci. Rem. Sens., 42, 650 [Google Scholar]
 Racine, R., Walker, G. A. H., Nadeau, D., Doyon, R., & Marois, C. 1999, PASP, 111, 587 [Google Scholar]
 Rajan, A., Rameau, J., De Rosa, R. J., et al. 2017, AJ, 154, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Rapacioli, M., Joblin, C., & Boissel, P. 2005, A&A, 429, 193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, ApJ, 852, 104 [CrossRef] [Google Scholar]
 Riaud, P., & Schneider, J. 2007, A&A, 469, 355 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieke, G. H., Wright, G. S., Buker, T., et al. 2015, PASP, 127, 584 [NASA ADS] [CrossRef] [Google Scholar]
 Ruffio, J.B., Macintosh, B., Wang, J. J., et al. 2017, ApJ, 842, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Ruffio, J.B., Macintosh, B., Konopacky, Q. M., et al. 2019, AJ, 158, 200 [Google Scholar]
 Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samland, M., Bouwman, J., Hogg, D. W., et al. 2021, A&A 646, A24 [CrossRef] [EDP Sciences] [Google Scholar]
 Schwartz, J. C., Sekowski, C., Haggard, H. M., Pallé, E., & Cowan, N. B. 2016, MNRAS, 457, 926 [Google Scholar]
 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]
 Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Nature, 509, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Snellen, I., de Kok, R., Birkby, J. L., et al. 2015, A&A, 576, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Soofbaf, S. R., Sahebi, M., & Mojaradi, B. 2018, Rem. Sens., 10, 434 [Google Scholar]
 Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Sparks, W. B., & Ford, H. C. 2002, ApJ, 578, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Barman, T., Skemer, A. J., et al. 2020, AJ, 160, 262 [Google Scholar]
 Thatte, N. A., Clarke, F., Bryson, I., et al. 2016, in Groundbased and Airborne Instrumentation for Astronomy VI, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 9908, 99081X [Google Scholar]
 Themelis, K. E., Schmidt, F., Sykioti, O., et al. 2012, Planet. Space Sci., 68, 34 [Google Scholar]
 Tsinos, C. G., Rontogiannis, A. A., & Berberidis, K. 2017, IEEE Trans. Comput. Imaging, 3, 160 [Google Scholar]
 Uyama, T., Currie, T., Hori, Y., et al. 2020, AJ, 159, 40 [Google Scholar]
 Vigan, A., Otten, G. P. P. L., Muslimov, E., et al. 2018, in Groundbased and Airborne Instrumentation for Astronomy VII, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 10702, 1070236 [Google Scholar]
 Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, in press, https://doi.org/10.1051/00046361/202038107 [Google Scholar]
 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]
 Wang, J., Mawet, D., Ruane, G., Hu, R., & Benneke, B. 2017, AJ, 153, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, J., Mawet, D., Fortney, J. J., et al. 2018, AJ, 156, 272 [NASA ADS] [CrossRef] [Google Scholar]
 WardDuong, K., Patience, J., Follette, K., et al. 2021, AJ, 161, 5 [Google Scholar]
 Wilcomb, K. K., Konopacky, Q. M., Barman, T. S., et al. 2020, AJ, 160, 207 [Google Scholar]
 Xie, C., Haffert, S. Y., de Boer, J., et al. 2020, A&A 644, A149 [EDP Sciences] [Google Scholar]
 Xu, Y., Wu, Z., Chanussot, J., et al. 2018, IEEE Trans. Geosci. Rem. Sens., 56, 1680 [Google Scholar]
 Yang, Y., Zhang, J., Song, S., & Liu, D. 2019, Rem. Sens., 11, 192 [Google Scholar]
 Zhang, X., Li, C., Zhang, J., et al. 2018, Rem. Sens., 10, 339 [Google Scholar]
 Zieleniewski, S., Thatte, N., Kendrew, S., et al. 2015, MNRAS, 453, 3754 [CrossRef] [Google Scholar]
All curves and constants are given in the HSIM source code https://github.com/HARMONIELT/HSIM
All Figures
Fig. 1 Model spectra of a β Pictorislike system: an A6 star with T_{eff} = 8000 K (BTNextgen Allard et al. 2012) and a giant planet with T_{eff} = 1700 K (BTSettl 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 H_{2}O and CO lines. 

In the text 
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 
Fig. 3 Final result of postprocessing the synthetic HARMONI data with a crosscorrelation 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 
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 highfrequency component of the injected spectrum (orange) along with the difference between the two (black). The firstovertone 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 
Fig. 5 Completeness curves (95%) with spectral unmixing (dark blue) and crosscorrelation (brown) expressed as the minimum planettostar contrast for a detection threshold of τ = 5σ as a function of the projected separation for Ltype (solid) and Ttype spectra (dotdashed) in the synthetic HARMONI data. 

In the text 
Fig. 6 Spatial weight inversion map of synthetic HARMONI data with a spiral of simulated planets (white circles) with a startoplanet 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 
Fig. 7 Final results of postprocessing VLT/SINFONI Kband data of β Pictoris b with the crosscorrelation (top) and spectral unmixing through orthogonal subspace projection (bottom). A BTSettl model spectrum (1700 K, log g = 4.0) was used for the crosscorrelation, 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 crosscorrelation strength C_{s}, and the weight α_{i}, respectively. The cross marks the position of β Pictoris A. 

In the text 
Fig. 8 Top: extracted continuumremoved SINFONI spectrum of β Pictoris b with spectral unmixing (dark blue), compared with a BTSettl model spectrum at 1700 K and log g = 4.0 dex (orange), close to the true planet properties. The firstovertone 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 bestfit Gaussian function (purple). 

In the text 
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 
Fig. B.1 Distributions of spatial weights (top) and crosscorrelation 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 crosscorrelation for clarity. The parameters of a Gaussian fit to the distribution of the crosscorrelation strengths are given. 

In the text 
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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.