# Medium-resolution integral-field spectroscopy for high-contrast exoplanet imaging. Molecule maps of the *β Pictoris* system with SINFONI

H. J. Hoeijmakers, H. Schwarz, I. A. G. Snellen, R. J. de Kok, M. Bonnefoy, G. Chauvin, A. M. Lagrange, J. H. Girard

<figure><img src="image_1.png" /></figure>

To cite this version:

H. J. Hoeijmakers, H. Schwarz, I. A. G. Snellen, R. J. de Kok, M. Bonnefoy, et al. Medium-resolution integral-field spectroscopy for high-contrast exoplanet imaging. Molecule maps of the *β Pictoris* system with SINFONI. Astronomy and Astrophysics - A&A, 2018, 617, 10.1051/0004-6361/201832902. insu-03635029

HAL Id: insu-03635029

https://insu.hal.science/insu-03635029

Submitted on 8 Apr 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Medium-resolution integral-field spectroscopy for high-contrast exoplanet imaging

## Molecule maps of the β Pictoris system with SINFONI*

H. J. Hoeijmakers<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, H. Schwarz<fnref n="4" />, I. A. G. Snellen<fnref n="1" />, R. J. de Kok<fnref n="1" />,<fnref n="5" />,<fnref n="6" />, M. Bonnefoy<fnref n="7" />, G. Chauvin<fnref n="7" />,<fnref n="8" />, A. M. Lagrange<fnref n="7" />, and J. H. Girard<fnref n="7" />,<fnref n="9" />

Received 25 February 2018 / Accepted 11 June 2018

## ABSTRACT

**Context.** Angular differential imaging (ADI) and spectral differential imaging (SDI) are well-established high-contrast imaging techniques, but their application is challenging for companions at small angular separations from their host stars.

**Aims.** The aim of this paper is to investigate to what extent adaptive-optics assisted, medium-resolution (*R* ~ 5000) integral field spectrographs (IFS) can be used to directly detect the absorption of molecular species in the spectra of planets and substellar companions when these are not present in the spectrum of the star.

**Methods.** We analysed archival data of the β Pictoris system taken with the SINFONI integral field spectrograph located at ESO's Very Large Telescope, originally taken to image β Pictoris b using ADI techniques. At each spatial position in the field, a scaled instance of the stellar spectrum is subtracted from the data after which the residuals are cross-correlated with model spectra. The cross-correlation co-adds the individual absorption lines of the planet emission spectrum constructively, while this is not the case for (residual) telluric and stellar features.

**Results.** Cross-correlation with CO and H₂O models results in significant detections of β Pictoris b with signal-to-noise ratios (S/Ns) of 13.7 and 16.4 respectively. Correlation with a *T* = 1700 K BT-Settl model provides a detection with an S/N of 22.8. This in contrast to application of ADI, which barely reveals the planet. While the adaptive optics system only achieved modest Strehl ratios of 19–27% leading to a raw contrast of 1:240 at the planet position, cross-correlation achieves a 3σ contrast limit of 2.7 × 10⁻<fnref n="5" /> in this 2.5 hr data set, a factor ~40 below the raw noise level at an angular distance of 0.36" from the star.

**Conclusions.** Adaptive-optics assisted, medium-resolution IFS, such as SINFONI on the VLT and OSIRIS on the Keck Telescope, can be used for high-contrast imaging utilizing cross-correlation techniques for planets that are close to their star and embedded in speckle noise. We refer to this method as molecule mapping, and advocate its application to observations with future medium resolution instruments, in particular ERIS on the VLT, HARMONI on the ELT and NIRSpec, and MIRI on the JWST.

**Key words.** infrared: planetary systems – techniques: imaging spectroscopy – planets and satellites: atmospheres – planets and satellites: detection – planets and satellites: gaseous planets

## 1. Introduction

To directly image an extrasolar planet, the light of its host star must generally be suppressed by orders of magnitude. The techniques that have been developed to do this rely on a combination of precise wave front control to restore the diffraction limit of the telescope (adaptive optics) and coronographic techniques to attenuate the stellar light (phase mask, phase/amplitude pupil apodization combined with focal plane mask or interferometers).

In such imaging data, residual starlight is present in the form of speckles which may mimic point-source objects, confusing the detection of companions and planets.

Such residuals can be suppressed by adopting differential imaging strategies that assume that the residual pattern scales with wavelength (spectral differential imaging; SDI), with the polarimetric state (polarimetric differential imaging; PDI) or is stable in time (angular differential imaging; ADI). Subsequent post-processing algorithms (cADI, LOCI, PCA, and ANDROMEDA) aim to optimize the residual attenuation while conserving the planetary signal (see e.g., Guyon 2011; Mawet et al. 2012; Chauvin 2016, for reviews).

<footnote n="1">¹ Leiden Observatory, Leiden University, 2333 CA Leiden, The Netherlands
e-mail: hoeijmakers@strw.leidenuniv.nl</footnote>

<footnote n="2">² Observatoire de Genève, Chemin des Maillettes 51, 1290 Versoix, Switzerland</footnote>

<footnote n="3">³ Universität Bern, Center for space and habitability, Gesellschaftstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="4">⁴ Department of Astronomy and Astrophysics, University of California, 1156 High St., Santa Cruz, CA 95064, USA</footnote>

<footnote n="5">⁵ SRON Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands</footnote>

<footnote n="6">⁶ Department of Physical Geography, Utrecht University, PO Box 80.115, 3508 TC Utrecht, The Netherlands</footnote>

<footnote n="7">⁷ Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="8">⁸ Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS/INSU UMI 3386 and Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile</footnote>

<footnote n="9">⁹ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote>* Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 093.C-0626.</footnote>
So far, direct imaging detections are generally limited to a specific part of the exoplanet population: young gas giants in wide orbits that glow by radiating out the internal heat that remained from the time of their formation. These planets can be resolved from their host stars owing to their relatively large mutual angular separation and brightness at infrared wavelengths. Such planets can generally only be observed during the first few tens of millions of years of their lifetime, after which they have cooled too much to be detected with current facilities (see e.g., *Bowler 2016*, for a review).

The desire to image cooler (i.e., older or less massive) planets that are closer to their host star is fueling the development of new instruments, as well as new observing and data analysis techniques. The arsenal of high-contrast imaging facilities has recently been expanded by the deployment of dedicated planet finding instruments such as the Gemini Planet Finder (Macintosh et al. 2006) at the Gemini Telescope, SPHERE (Beuzit et al. 2008) on the European Very Large Telescope (VLT), and SCExAO (Jovanovic et al. 2015) on the Subaru telescope. These instruments also have low-resolution ($R \sim 30–400$) spectroscopic capability – allowing the spectral characterization of directly imaged planets and achieving planet-to-star contrast ratios down to $10^{-6}$ at $0.2''$ angular separation (see e.g., Macintosh et al. 2014; Ruffio et al. 2017; Mesa et al. 2017; Currie et al. 2017).

Multiple studies indicate that when high-contrast imaging is combined with high-dispersion spectroscopic techniques, the achieved contrasts can be significantly enhanced (see e.g., *Sparks & Ford 2002*; *Riaud & Schneider 2007*; *Kawahara et al. 2014*; *Snellen et al. 2015*; *Luger et al. 2017*; *Wang et al. 2017*; *Lovis et al. 2017*). This strategy assumes that the faint planet and the much brighter star have spectral properties that are distinctly different at high spectral resolution; notably because of molecular absorption bands in the spectrum of the planet. Because the planet is close to the star, its spectrum is deeply embedded in speckle-noise, but it can be extracted using cross-correlation: the cross-correlation co-adds the individual absorption lines of the planet constructively but not stellar and telluric features, or at different radial velocities. This method has been applied successfully for the first time at high spectral resolution, albeit using the CRIRES slit spectrograph that probes only one spatial dimension. This resulted in the measurement of the spin rate of $\beta$ Pictoris b (Snellen et al. 2014).

In this paper, we investigate to what extent this method can be applied to observations from adaptive-optics assisted, medium-resolution ($R \sim 5000$) integral field spectrographs (IFS) – targeting molecular species in a planet atmosphere that are not present in the star, throughout the two-dimensional (2D) field of view. The use of cross-correlation techniques on IFS data was first performed by *Konopacky et al. (2013)* and *Barman et al. (2015)* to detect H$_2$O, CO and CH$_4$ in the atmospheres of HR 8799 b&c using the OSIRIS integral field spectrograph at the Keck Observatory. At $1.7''$ and $0.96''$, HR 8799 b&c are widely separated from their host star (Marois et al. 2008). Such angular distances are well resolvable by modern adaptive-optics systems, and the speckle-pattern can be effectively removed using ADI or SDI-based methods. However, the application of ADI and SDI are less effective at smaller angular separations because the effects of field rotation and wavelength-dependencies are limited (Fitzgerald & Graham 2006; Lafrenière et al. 2007; Marois et al. 2008; Rameau et al. 2015).

We apply the cross-correlation technique to archival K-band SINFONI image cubes of the $\beta$ Pictoris system, and are able to extract the spectral signature of the planet $\beta$ Pictoris b while

effectively removing the diffraction pattern of the star in which the planet is embedded. We refer to this technique as molecule mapping, because it produces 2D cross-correlation maps that indicate the presence of certain molecular signatures for each location in the 2D image.

In Sect. 2, we summarize the $\beta$ Pictoris system. Section 3 describes the archival SINFONI observations used in this analysis and the data reduction. The cross-correlation procedure is described in Sect. 4, followed by the resulting cross-correlation images in Sect. 5, including a comparison with ADI using classical ADI (Marois et al. 2006) and the LOCI algorithm (Lafrenière et al. 2007). Sections 5.4 and 5.5 highlight the potential application of molecule mapping with upcoming medium-resolution integral-field instruments on the VLT, ELT, and JWST, after which the paper is summarized and concluded in Sect. 6.

## 2. The $\beta$ Pictoris system

In 1984, $\beta$ Pictoris was the first star to be found to host a debris disk and which was associated with planet formation processes (Smith & Terrile 1984). The warped structure of the disk led to the inference of a possibly planetary mass companion (Scholl et al. 1993; Roques et al. 1994; Lazzaro et al. 1994; Burrows et al. 1995; Lecavelier des Etangs et al. 1996; Mouillet et al. 1997; Augereau et al. 2001), which was discovered in 2008 via direct imaging (Lagrange et al. 2009, 2010). The system is a member of the nearby $\beta$ Pictoris moving group which has an estimated age of $24 \pm 3$ Myr (Bell et al. 2015). The bolometric luminosity and effective temperature of the planet were recently estimated by Chilcote et al. (2017) at $\log(L_{\text{bol}}/L_{\odot}) = -3.78 \pm 0.03$ and $1724 \pm 15$ K respectively. From these, the radius and mass of the planet were inferred to be $12.9 \pm 0.2 M_J$ and $1.46 \pm 0.01 R_J$ respectively (Chilcote et al. 2017), assuming a hot-start evolution model (Baraffe et al. 2003). These properties are summarized in Table 1.

$\beta$ Pic b moves on a 20–26-yr orbit that is highly inclined with respect to the line of sight from Earth (Wang et al. 2016). This high inclination has recently triggered a renewed interest in this system because the Hill sphere of the planet was predicted to transit the star between April 2017 and January 2018 (Lecavelier des Etangs & Vidal-Madjar 2016; Wang et al. 2016; Kenworthy 2017). Snellen et al. (2014) measured the projected equatorial rotation velocity of the planet, inferring a length of day of ~8 hr.

## 3. Observations

The $\beta$ Pictoris system was observed with the SINFONI IFS (Eisenhauer et al. 2003; Bonnet et al. 2004) mounted on ESO's Very Large Telescope in K-band on the nights of September 10 and 11, 2014 as part of ESO programme 093.C-0626(A) (PI: Chauvin). The observations on both nights were taken in pupil tracking mode which allows the field to rotate during the observing sequence, facilitating Angular Differential Imaging. The 32 × 64 pixel images were taken in the highest spatial resolution mode, with each pixel covering $0.0125''$ by $0.025''$, providing a field of view (FOV) of $0.8''$ by $0.8''$ on sky. At a resolving power of $R \sim 5000$ the SINFONI spectra have a wavelength coverage from $1.929 \mu$m to $2.472 \mu$m with a sampling of $0.25$ nm. At the time of these observations the planet was separated from the central star by approximately $0.36''$ (Wang et al. 2016). On the first and second night, 24 and 30 science frames were obtained respectively, spanning a total duration of 2.5 hr
<table><caption>**Table 1.** Summary of the properties of the star β Pictoris (upper part) and its planet (lower part).</caption><thead><tr><th></th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td colspan="3"><strong>Star:</strong></td></tr><tr><td>Visible magnitude<sup>a</sup></td><td>V</td><td>3.86</td></tr><tr><td>K-band magnitude<sup>a</sup></td><td>K</td><td>3.48</td></tr><tr><td>Distance (pc)<sup>b</sup></td><td>d</td><td>19.44 ± 0.05</td></tr><tr><td>Effective temperature (K)<sup>c</sup></td><td>T<sub>eff</sub></td><td>8052 ± 30</td></tr><tr><td>Mass (M<sub>☉</sub>)<sup>d</sup></td><td>M<sub>*</sub></td><td>1.85<sup>+0.03</sup><sub>-0.04</sub></td></tr><tr><td>Metallicity (dex)<sup>c</sup></td><td>[M/H]</td><td>0.05 ± 0.06</td></tr><tr><td>Age (Myr)<sup>e</sup></td><td></td><td>24 ± 3</td></tr><tr><td>Rotation velocity (km s<sup>-1</sup>)<sup>f</sup></td><td>v sin i</td><td>130</td></tr><tr><td>Systemic velocity (km s<sup>-1</sup>)<sup>g</sup></td><td>v</td><td>20.0 ± 0.7</td></tr><tr><td colspan="3"><strong>Planet:</strong></td></tr><tr><td>Luminosity (log L/L<sub>☉</sub>)<sup>h</sup></td><td>L<sub>p</sub></td><td>-3.78 ± 0.03</td></tr><tr><td>K-band contrast<sup>i</sup></td><td>ΔK</td><td>9.2 ± 0.1</td></tr><tr><td>Effective temperature (K)<sup>h</sup></td><td>T<sub>eff</sub></td><td>1724 ± 15</td></tr><tr><td>Surface gravity (log g)<sup>h</sup></td><td></td><td>4.18 ± 0.01</td></tr><tr><td>Angular separation (mas)<sup>j</sup></td><td></td><td>356.5 ± 0.9</td></tr><tr><td>Position angle (deg)<sup>j</sup></td><td></td><td>213.0 ± 0.2</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Ducati (2002). <sup>(b)</sup>van Leeuwen (2007). <sup>(c)</sup>Gray et al. (2006). <sup>(d)</sup>Wang et al. (2016); the stellar mass was inferred from their fit to the total mass of the system minus the mass of the planet. Although the quoted confidence intervals correspond to 1σ uncertainties, the posterior distributions are non-Gaussian. <sup>(e)</sup>Bell et al. (2015). <sup>(f)</sup>Royer et al. (2007). <sup>(g)</sup>Gontcharov (2006). <sup>(h)</sup>Chilcote et al. (2017); these authors measured the luminosity from their GPI datasets, but inferred the other parameters using a hot-start evolutionary model by Baraffe et al. (2003). The reported statistical errors do not account for the model-dependency of these values. <sup>(i)</sup>Bonnefoy et al. (2011). <sup>(j)</sup>Measurements obtained by Wang et al. (2016) on November 8, 2014.

<table><caption>**Table 2.** Overview of the observations of the β Pictoris system of September 10 and 11, 2014.</caption><thead><tr><th>Observing night</th><th>Sept. 10, 2014</th><th>Sept. 11, 2014</th></tr></thead><tbody><tr><td>N<sub>exp</sub></td><td>24</td><td>30</td></tr><tr><td>DIT (s)</td><td>60</td><td>2</td></tr><tr><td>NDIT</td><td>4</td><td>50</td></tr><tr><td>Airmass</td><td>1.34 – 1.15</td><td>1.49 – 1.14</td></tr><tr><td>PA (deg)</td><td>340.7 – 304.6</td><td>355.7 – 304.3</td></tr><tr><td>Strehl ratio (%)</td><td>19 – 27</td><td>20 – 25</td></tr></tbody></table>

**Notes.** The table lists the number of exposures N<sub>exp</sub>, the exposure time per dithering position (DIT), the number of dithering positions combined to make each science frame (NDIT), and the range of airmasses, position angles on sky and approximate Strehl ratios achieved.

(see Table 2). During both nights, the seeing varied between 0.7" and 0.9", and the MACAO AO system delivered modest Strehl ratios between 19% and 27% (Bonnet et al. 2003). The observations on the first night were preceded by a few acquisition images after the star was placed just outside the FOV for the rest of the sequence, enabling longer exposure times. During the second night both star and planet were observed continuously (see Fig. 1).

The raw data were downloaded from the ESO Science Archive Facility and reduced using version 3.0.0 of the SINFONI pipeline. The reduction pipeline produces a three-dimensional (3D) image cube for each science observation, with sky position in the x and y directions and wavelength in the z direction. The pipeline reduced data cubes contain NaN values at the edges of the waveband and at the location of known bad pixels. To reject these areas, we only consider wavelengths between 2.088 µm and 2.452 µm and flag any remaining NaN values (see Fig. 2). From each science-exposure we obtain the wavelength-averaged 2D image. The spatial location of the star was found by fitting a 2D Gaussian profile to a region of 10 × 10 pixels around the maximum of the PSF. We use this location as a pivot to de-rotate all image-cubes by their respective position angle and to align them to a common frame, both using linear interpolation. We co-added all image cubes of each night to obtain two master image-cubes, excluding all regions within a range of 5 pixels (~60 mas) from the edges of the individual frames. The wavelength-averaged images of these two master cubes are shown in Fig. 1. From the wavelength-averaged image of the second night we measured the raw PSF contrast and the noise as a function of angular separation by measuring the mean and the standard deviation of the flux in annuli centered on the star, normalized by the peak stellar flux (see Fig. 3). The contrast in intensity between the flux at the location of the planet with respect to the peak intensity of the PSF is measured to be 1:240, and the contrast between the 1σ noise level and the PSF peak is 3.6 × 10<sup>-4</sup>.

We find a spurious radial velocity variation of the telluric absorption lines in the background spectra of both nights, increasing from +20 to +30 km s<sup>-1</sup> during the observations, indicating instability in the wavelength solution at the 1 pixel level. We correct these anomalies by shifting each datacube to the rest-frame of the telluric spectrum.

# 4. Data analysis

## 4.1. Removal of starlight

The flux at each spatial position in the data cube is dominated by the spectrum of the star, showing a few stellar but mainly telluric absorption features. Our removal of the star relies on subtracting a scaled template from each spectrum in the datacube and a high-pass filter to remove the effect of wavelength-dependence of the PSF. The steps to remove this starlight are shown in Fig. 4 and are expanded upon below.

First, a high signal-to-noise (S/N) master stellar spectrum is created by combining the spectra from the 1% (20) brightest pixels in the white-light image. These spectra are normalized to the same flux level and then averaged at each wavelength while rejecting >6σ outliers.

For each spatial position in the data cube we subsequently divide the spectrum by this master and apply a Gaussian smoothing filtering with a 1σ width of 10 wavelength steps (FWHM = 762 km s<sup>-1</sup>). The resulting function is used as a proxy of the local wavelength-dependence of the stellar diffraction pattern. We multiply the stellar master spectrum with this proxy to obtain a model of the stellar spectrum at each location of the image, and subtract these from the data. This removes the effect of the wavelength-dependent scale of the PSF and a possible wavelength-dependence in the location of the star, but also the continuum at the location of the planet. We reject remaining 6σ outliers to remove cosmic ray hits and bad pixels that were not corrected by the pipeline.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Wavelength-averaged images of the β Pictoris SINFONI data, obtained by de-rotating and stacking the exposures of each night, and taking the median flux of the spectrum at each image location. The expected position of β Pictoris b is indicated by the red crosses. The coordinate system is relative to the position of the star.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Average stellar spectrum of the master cube of the second night (black), compared to a continuum-normalized rotation-broadened PHOENIX model in red (Husser et al. 2013) and a telluric transmission spectrum as obtained using ESO SkyCalc in blue (Noll et al. 2012; Jones et al. 2013), both convolved to the spectral resolution of SINFONI. There are few stellar absorption lines in this wavelength range, and the data is dominated by telluric absorption bands due to water, CO₂ and methane. The grey regions at the edges of the waveband were discarded due to strong telluric bands of H₂O and CO₂ or bad pixels that are close to the edges of the detector.</figcaption></figure>

At this stage, residual structures remain in the image at the position of strong absorption lines, most importantly the strong telluric CH₄ feature at 2.32 µm that dominates the absorption spectra of both hot and cold CH₄ gas. Such residuals are caused by mismatches between the local spectra and the stellar master spectrum. This effect notably occurs near the edges of individual frames and at the location of the star, where the S/N of the residuals is high. We attribute such systematic residuals to stray light in the instrument that was not corrected by the data reduction pipeline, diluting the spectra, and hence decreasing the apparent depth of absorption features in the stellar/telluric spectrum, as well as remaining inaccuracies in the wavelength solution.

To mitigate these effects we unwrap the 3D residual cube (i.e., the residuals after subtraction of the stellar spectrum) of each exposure to its original 2D format and remove correlated noise patterns by iteratively applying the SYSREM algorithm (Tamuz et al. 2005) on these 2D frames (see Fig. 8). We found that applying more than 8 iterations did not improve

the end result further. After application of SYSREM, the 2D frames are folded back into 3D cubes. Finally, these cubes are derotated and aligned as described in Sect. 3 and co-added into two master residual data-cubes, one for each night, onto which we apply the cross-correlation analysis.

## 4.2. Cross correlation with model templates

We searched for the presence of CO, H₂O, NH₃, and CH₄ at each spatial location in the data cube by cross-correlating each spectrum in the residual data-cubes with a model template of each of these molecules, corresponding to ~1200 K emission models of the day side HD 189733 adopted from de Kok et al. (2014). In addition we cross-correlated with a grid of pre-computed BT-Settl model spectra<fnref n="1" /> with temperatures <3000 K and varying in surface gravity from log(g) = 2.5 to 5.5 (Allard et al. 2011).

<footnote n="1">¹ The BT-Settl model grid was obtained from https://phoenix.ens-lyon.fr/Grids/BT-Settl/CIFIST2011_2015/FITS/.</footnote>
<figure><img src="image_5.png" /></figure>

**Fig. 3.** Raw-contrast curve obtained by averaging the flux from the second night in annuli with widths of 20 mas around the central star, normalized by the peak flux. The standard deviation of the flux in these annuli is plotted with the dashed line, and is dominated by speckle noise. The $3\sigma$ sensitivity limit of the analysis presented in this work is derived from our $S/N = 22.8$ detection of the planet, assuming that the actual K-band contrast of $\beta$ Pic b is $K = 9.2$ (Bonnefoy et al. 2011).

Prior to cross-correlation these models are filtered with the same $762 \text{ km s}^{-1}$ filter as the data to eliminate low-frequency variations and to maximize the correlation with the continuum-subtracted planet spectrum. They are then convolved to a spectral resolution of $R = 5000$ (see Fig. 5) and cross-correlated with each residual spectrum over a range of velocities of $v \pm 2500 \text{ km s}^{-1}$ in steps of $10 \text{ km s}^{-1}$. The cross-correlation of these two data cubes therefore results in two new cubes that contain the cross-correlation coefficients of each spectrum for 501 steps in radial velocity.

The closer the model template matches the real emission spectrum of the planet, the larger the cross-correlation function at the location and radial velocity of planet will be. At all other locations and radial velocities, no significant cross-correlation signals are expected. The systematic radial velocity of the $\beta$ Pictoris system is $+20 \pm 0.7 \text{ km s}^{-1}$ (Gontcharov 2006). Using the orbital solution of Wang et al. (2016), we calculate the instantaneous radial velocity of $\beta$ Pic b to be $\sim -9.5 \text{ km s}^{-1}$ with respect to the star. The barycentric velocity is $-7.9 \text{ km s}^{-1}$ so the cross-correlation signal of the planet is expected to occur at a radial velocity of $\sim 3 \text{ km s}^{-1}$, i.e., at the central slice of the cross-correlation cubes.

After cross-correlation, the cross-correlation cubes associated with each night are aligned based on the known location of the star in the field (see Sect. 3) and subsequently averaged, weighed by the square root of the total exposure time in each night. This yields a single cross-correlation cube for each model template, the central slice of which contains the molecule maps of $\beta$ Pictoris b.

### 4.3. ADI Analysis

We also analysed the data by means of angular differential imaging (ADI) to compare the performance of the molecule mapping

<figure><img src="image_6.png" /></figure>

**Fig. 4.** Stepwise removal of starlight from each spatial location in the data cubes. Upper panel: master stellar spectrum (red) compared to the local observed spectrum (black) at spatial pixel location [10,40] of exposure #7 of the second night. Both are normalized and offset to allow visual comparison. Middle panel: ratio of the local to the master spectrum before (black) and after Gaussian smoothing (red). This smoothed residual is used as a proxy for the local wavelength-dependence of the stellar diffraction pattern. Bottom panel: residuals obtained after subtracting the master stellar spectrum multiplied by the proxy from the local observed spectrum. Pixel values that are more than $6\sigma$ away from the mean are rejected (red). SYSREM is subsequently applied to these residuals (not shown).

technique. We applied the LOCI algorithm (Lafrenière et al. 2007) on the sequence obtained during the second night. LOCI was applied independently at each wavelength with a separation criterion of 1 full-width-at-half-maximum, an optimization area of 300 FWHM and a geometry parameter of 0.5, after which the images at the individual wavelengths were stacked. During the first night of data the telescope was offset such that the star was located outside the field of view. With the field rotating around the axis of the telescope (image center), this field asymmetry causes the PSF of the star to evolve during the night due to field rotation and variable atmospheric refraction, hindering the application of LOCI. We therefore applied the classical ADI algorithm (Marois et al. 2006) to the data obtained during the first night instead.

## 5. Results and discussion

### 5.1. Molecule maps of the $\beta$ Pictoris system

The molecule maps at the radial velocity of the planet are shown in Fig. 6. The maps of CO and $H_2O$ show a significant signal at the expected position of the planet, while the planet is not detected in the maps of $NH_3$ and $CH_4$. The one-dimensional (1D) cross-correlation functions at the position of the planet are shown in Fig. 7. To determine the S/N of these signals, we measure the average standard deviation of the cross-correlation
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Model templates of CO, H₂O, CH₄, and NH₃ at high (R ~ 10⁶) spectral resolution (light colour) and convolved to a spectral resolution of R = 5000 (dark colour). The vertical dashed lines indicate the wavelength range of the data. Before correlation, these models are filtered with the same 762 km s⁻¹ filter as the data (not shown).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Molecule maps of CO, H₂O, CH₄, and NH₃ at vₛsys = 0 km s⁻¹. All four panels: colours scale linearly between cross-correlation values of -0.05 (black) to +0.3 (white). A cross-correlation enhancement caused by the planet is detected at S/Ns of 13.7 and 16.4 in the maps of CO and H₂O respectively, but not in those of CH₄ and NH₃.</figcaption></figure>

function in an annulus between 8 and 11 pixels away from the peak signal, at radial velocities more than ±250 km s⁻¹ away from the planet. This is done to avoid systematic variations in the cross-correlation function at the location of the planet due to autocorrelation, which is not negligible due to the strength of the main cross-correlation peak. In this way, we establish S/Ns for CO and H₂O of 13.7 and 16.4 respectively. The non-detections of CH₄ and NH₃ are consistent with the effective temperature of the planet, which at 1724 K is expected to be too high for these molecules to be present at significant quantities unless the

planet is rich in respectively carbon or nitrogen (Burrows et al. 2006; Hubeny & Burrows 2007; Zahnle & Marley 2014; Heng & Tsai 2016; Heng & Lyons 2016; Moses et al. 2016; Todorov et al. 2016).

## 5.2. Planet characterization using molecule mapping

Variations in global atmospheric parameters such as the effective temperature, surface gravity, metallicity, and abundances ratios affect the relative strength of individual absorption lines.
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** One-dimensional (1D) cross-correlation functions of CO, H₂O, CH₄, and NH₃ at the location of the planet.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Effect of SYSREM by comparison between the un-corrected data (**left panels**) and the data after 8 passes through SYSREM (**right panels**). **Upper and middle panels:** unfolded residual datacubes after cleaning of the stellar spectrum. Remnants of the storng telluric CH₄ feature are visibly present before application of SYSREM. The white region is enlarged in the **middle panels**. These are removed after 8 passes of SYSREM (**right panel**). **Bottom panels:** cross-correlation maps of H₂O obtained before (**left panel**) and after (**right panel**) SYSREM.</figcaption></figure>

This may significantly influence the cross-correlation function, which means that the analysis is potentially sensitive to underlying model parameters. Figure 10 shows the peak value of the 1D cross-correlation function with BT-Settl models with varying $T_{\text{eff}}$ and $\log(g)$. These models take into account absorption

from multiple species (including CO and H₂O) and are therefore a more complete representation of the true planet spectrum than the individual molecule spectra. The cross-correlation peak is observed to steeply decrease for temperatures below ~1600 K because the strength of the water absorption features at wavelengths larger than 2.1 µm become significantly weaker at lower temperatures.

Above ~2000 K the cross-correlation peak is a shallower function of $T_{\text{eff}}$ as the CO and H₂O absorption bands slowly diminish towards higher temperatures. Similarly, models with $\log(g)$ below 5.0 are favoured over models with higher surface gravities, due to the stronger water absorption bands relative to CO for lower values of $\log(g)$. The highest cross-correlation peak is achieved for $T_{\text{eff}} = 1700$ K and $\log(g) = 3.5$, with an S/N of 22.8 (see Fig. 9).

This result is in line with the values reported by Chilcote et al. (2017), who obtain a consistent effective temperature and surface gravity by fitting hot-start evolutionary models (Baraffe et al. 2003) to the bolometric luminosity of β Pic b. Our analysis demonstrates that these parameters can also be obtained from the medium resolution spectrum through cross-correlation. The somewhat higher surface gravity by Chilcote et al. (2017) is likely explained by the fact that these parameters are model dependent in both fitting approaches. It is also conceivable that the discrepancy is caused by the fact that the inferred parameters are affected by the subtraction of the continuum that is inherent to our analysis, while it is generally preserved with low-resolution ADI/SDI based methods. More work should therefore be done to survey the possible inconsistencies across models, fitting methods, wavelength ranges and PSF subtraction techniques.

More generally, we conclude that integral-field spectrographs can be used to characterize the fundamental parameters of directly imaged companions via spectral cross-correlation and that molecule mapping is therefore not limited to the detection of molecules in the atmospheres of these objects. Moreover, we predict that spectra covering a wider wavelength range that includes multiple water bands and bands of additional molecules will allow the effective temperature and surface gravity to be constrained significantly further using molecule mapping.

Finally, we state the obtained contrast limit obtained when cross-correlating the data with the $T_{\text{eff}} = 1700$ K/$\log(g) = 3.5$ template. The K-band contrast ratio between $\beta$ Pic b and its host star is $\Delta K = 9.2 \pm 0.1$ (Bonnefoy et al. 2011). At an S/N of 22.8,
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Left panel: map of the cross-correlation function using the preferred BT-Settl model ($T_{\text{eff}}$ = 1700 K and $\log(g)$ = 3.5). Right panel: one-dimensional (1D) cross-correlation function at the location of the planet.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Planet cross-correlation peak as a function of BT-Settl models with varying $T_{\text{eff}}$ and $\log(g)$. The cross-correlation peak steeply decreases for temperatures below ~1600 K because the strength of the water absorption features at >2.1 µm become significantly weaker. All molecular bands diminish above temperatures of ~2000 K, causing a slow decrease of the cross-correlation function at higher temperatures. The black circle corresponds to the values reported by Chilcote et al. (2017).</figcaption></figure>

our detection of the planet therefore corresponds to an achieved $3\sigma$ contrast of $2.7 \times 10^{-5}$, a factor 40 deeper than the raw noise level measured in Sect. 3 (see Fig. 3).

## 5.3. Comparison with ADI-LOCI

The upper panels of Fig. 11 show the white-light images obtained using ADI, for which these observations were initially intended. These reveal that there is a strong difference in performance between the two observing nights: because LOCI is hindered by field asymmetry, the analysis of the data obtained during the first night is limited to classical ADI, which does not result in a detection of the planet. The planet is retrieved when using LOCI on the second night of data, but several other positive and negative features with similar amplitudes as that of the planet are present in the image, in particular close to the host star.

In contrast, the planet is detected most strongly in the molecule map of the observations taken in the first night, despite the field asymmetry that hinders application of ADI. The observations taken in the second night have a much shorter exposure time than the first night, which was needed to place both the star and the planet in the field of view without

saturating the detector<fnref n="2" />. The standard day-time calibration plan of the SINFONI instrument includes the observation of three dark frames at each exposure time scheduled for use at night-time. This means that the dark current in this observing sequence is measured using only three dark frames with an exposure time of 2 s each. We hypothesize that the observations of the second night are dominated by the combined read-noise in the master dark-frame and that this causes the strongly reduced sensitivity as compared to the first night. Future observations with SINFONI that target faint companions should therefore include a larger number of dark frames on top of those obtained in the standard calibration plan.

Nevertheless, our molecule mapping approach more effectively suppresses the residual speckle pattern than ADI in both nights of data, leading to a stronger detection of $\beta$ Pic b and a higher sensitivity close to the star due to the absence of residual

<footnote n="2">² This effect is aggravated by the fact that the SINFONI detector shows strong persistence effects for count levels exceeding 18 000 counts that dissipate over timescales up to 48 hr. The exposure time is therefore limited such that the nominal flux does not exceed 8000 counts (see the SINFONI user manual for details).</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** *Upper panels:* white-light images obtained after co-adding the data using classical ADI (**left panel**) and with the LOCI algorithm (**right panel**) for the first and second nights respectively. The planet is visible near [0.2", -0.3"], but is barely recovered in the second night, and not distinguishable from speckle noise in the first night. *Lower panels:* cross-correlation maps using the preferred BT-Settl model ($T_{\text{eff}}$ = 1700 K and $\log(g)$ = 3.5). The colour scales linearly from black to white between -0.1 and +0.5 (**left panel**) and -0.05 to +0.25 (**right panel**).</figcaption></figure>

speckles. PSF subtraction techniques that do take the spectral information into account (such as SDI or forward-modeling of the PSF) may also be more sensitive to the companion than ADI, but we expect that these techniques are similarly hindered by the field-asymmetry in the first night of data. More work is needed to perform a comprehensive comparison between methods, and this will be addressed in a future study.

## 5.4. Molecule mapping from space with NIRSpec and MIRI

Molecule mapping provides new possibilities for high-contrast imaging using medium-resolution IFSs, also using space-based observatories. Both NIRSpec and MIRI on the James Webb Space Telescope (JWST), will have integral-field capabilities. In IFS mode, NIRSpec<fnref n="3" /> covers wavelengths from 0.7 to 5.27 µm at a spectral resolution of ~2.700 with a spatial sampling of 0.1" × 0.1". MIRI has a similar IFS mode<fnref n="4" />, operating between 4.89 and 28.45 µm at a spectral resolution varying from 3320 at short wavelengths to 1460 at longer wavelengths. The waveband is divided into four channels with increasing slit-size, the bluest channel having a spatial resolution of 0.196", and the reddest a spatial resolution of 0.276". Although the spatial sampling of both NIRSpec and MIRI is significantly more coarse than that of SINFONI in 25 mas plate-scale mode, objects with a separation down to 0.1" can theoretically be resolved by the IFS of

NIRSpec, and separations down to 0.2" can be resolved by MIRI at 5 µm.

## 5.5. Molecule mapping with ERIS and HARMONI

The ERIS integral-field spectrograph is expected to replace SINFONI at the Cassegrain focus of UT4 at the VLT in 2020 (Amico et al. 2012; Kuntschner et al. 2014). It will use an upgraded version of the existing IFS SPIFFY and a new infrared camera to achieve higher throughput and a spectral resolution of R ~ 8000 (George et al. 2016). In addition, it will employ a new wavefront module to enable significantly higher Strehl ratio's than SINFONI. These modifications will therefore increase the contrast and sensitivity to molecular signatures, making ERIS an important asset in the application of the molecule mapping technique until the arrival of the ELT.

The primary spectroscopic capability of the ELT in the near infra-red will be fulfilled by the HARMONI integral-field spectrograph. HARMONI will cover wavelengths between 0.47 to 2.45 µm at spectral resolutions ranging from R ~ 400 to ~20 000 in individual Z, J, H, and K bands. The highest possible spatial sampling is 4 mas, equivalent to the diffraction limit of the ELT (Thatte et al. 2014).

To first order, the S/N achieved by cross-correlation scales with the S/N per wavelength element of the planet's absorption lines times the square root of the number of lines in the waveband (Snellen et al. 2015). For bright stars at near-infrared wavelengths, we can assume that the noise is dominated by the photon noise of the stellar PSF, the size of which scales with the inverse of the mirror diameter D. The S/N of the planet

<footnote n="3">³ NIRSpec web documentation: https://jwst-docs.stsci.edu/display/JTI/Near+Infrared+Spectrograph%2C+NIRSpec.</footnote>

<footnote n="4">⁴ MIRI web documentation: https://jwst-docs.stsci.edu/display/JTI/Mid-Infrared+Instrument%2C+MIRI.</footnote>
spectral lines scales with the square root of the spectral resolving power, so the S/N of HARMONI compared to SINFONI can be approximated as:

$$ \frac{S/N_{\text{HARMONI}}}{S/N_{\text{SINFONI}}} = \frac{D_{\text{ELT}} \sqrt{R_{\text{HARMONI}}}}{D_{\text{VLT}} \sqrt{R_{\text{SINFONI}}}} = 9.75 $$

assuming the maximum spectral resolution of HARMONI, but otherwise equal wavelength coverage and instrument throughput. This indicates that HARMONI can achieve the same sensitivity as SINFONI in ~1% of the exposure time, but for planets that are five times closer to their host star because the ELT mirror is five times bigger than that of the VLT.

With a spectral resolving power of 20 000 the systemic and orbital velocities of target planets can easily be resolved by the cross-correlation function if they exceed 15 km s⁻¹. In addition to discovering a companion and spectrally characterizing its atmosphere, the HARMONI instrument can therefore be used to establish co-movement with the host star and constrain the planet orbit via the instantaneous radial velocity. We conceive that this can be achieved in a single observation, making HARMONI especially suitable for detailed characterization of young giant exoplanets in short amounts of exposure time.

Mid-infrared coverage of the ELT is provided by METIS at high spectral resolution (~100 000). Like HARMONI, METIS provides IFU capability for diffraction-limited imaging. It has already been recognized that this combination of high resolution spectroscopy with high contrast imaging gives METIS a unique capacity to characterize exoplanet atmospheres, using analysis strategies that are very similar to the one presented in this work (Snellen et al. 2015). Our work shows that the same strategy can be employed in the near-infrared at medium spectral resolution using HARMONI.

## 6. Conclusion

This paper introduces the first application of the molecule mapping technique for detecting close-in substellar companions using AO-assisted medium-resolution integral field spectroscopy. With molecule mapping, integral-field spectra are cross-correlated with molecular template spectra to search for the spectral signatures of spatially resolved exoplanets that are embedded in the photon noise of the dominating star light. The cross-correlation co-adds the individual absorption lines of the planet spectrum at the spatial location of the planet, but not (residual) telluric and stellar features. This acts to suppress the quasi-static speckle pattern that is a limiting factor in standard direct-imaging analyses.

We applied molecule mapping to 2.5 hr of archival SINFONI observations of the β Pictoris system. By cross-correlating the integral-field data cubes with templates of CO, H₂O, CH₄, and NH₃, we obtained S/Ns for CO and H₂O in the atmosphere of the planet of 13.7 and 16.4 respectively. We also cross-correlated with a grid of BT-Settl models, varying the effective temperature between 1200 and 3000 K, and the surface gravity between log(g) = 2.5-5.5. We found that the cross-correlation peaks for T<sub>eff</sub> = 1700 K and log(g) = 3.5, which shows that molecule mapping can also constrain the fundamental parameters of young gas giant planets. With these model parameters, the planet was detected at an S/N of 22.8, corresponding to a 3σ contrast of 2.7 × 10⁻⁵. We also analysed the same data using ADI, resulting in only a marginal detection and a strong difference in sensitivity in the two individual nights, due to the different observing

strategies used and the associated calibrations. However, in both nights of data molecule mapping outperformed ADI, demonstrating that it can significantly enhance the sensitivity of IFS observations at medium spectral resolution.

Our successful application of molecule mapping on existing SINFONI data feeds expectations for the potential of upcoming medium-resolution integral-field instruments on the JWST and the ELT. We briefly outlined the specifications of MIRI and NIRSpec (JWST), ERIS (VLT), and HARMONI (ELT), and anticipate that these instruments are well suited for such cross-correlation based analyses.

*Acknowledgements.* This work is part of the research programme VICI 639.043.107 funded by the Dutch Organisation for Scientific Research (NWO). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement numbers 694513 and 724427. This work is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 093.C-0626(A). Finally, G. Chauvin, M. Bonnefoy and A.-M. Lagrange acknowledge support from the French National Research Agency (ANR) through project grant ANR10-BLANC0504-01 and the Programmes Nationaux de Planétologie et de Physique Stellaire (PNP & PNPS), in France.

## References

Allard, F., Homeier, D., & Freytag, B. 2011, in 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, eds. C. Johns-Krull, M. K. Browning, & A. A. West, ASP Conf. Ser., 448, 91

Amico, P., Marchetti, E., Pedichini, F., et al. 2012, in Ground-based and Airborne Instrumentation for Astronomy IV, Proc. SPIE, 8446, 844620

Augereau, J. C., Nelson, R. P., Lagrange, A. M., Papaloizou, J. C. B., & Mouillet, D. 2001, A&A, 370, 447

Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61

Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593

Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, in Ground-based and Airborne Instrumentation for Astronomy II, Proc. SPIE, 7014, 701418

Bonnefoy, M., Lagrange, A.-M., Boccaletti, A., et al. 2011, A&A, 528, L15

Bonnet, H., Ströbele, S., Biancat-Marchet, F., et al. 2003, in Adaptive Optical System Technologies II, Proc. SPIE, 4839, 329

Bonnet, H., Abuter, R., Baker, A., et al. 2004, The Messenger, 117, 17

Bowler, B. P. 2016, PASP, 128, 102001

Burrows, C. J., Krist, J. E., Stapelfeldt, K. R., & WFPC2 Investigation Definition Team. 1995, BAAS, 27, 1329

Burrows, A., Sudarsky, D., & Hubeny, I. 2006, ApJ, 640, 1063

Chauvin, G. 2016, in Young Stars & Planets Near the Sun, eds. J. H. Kastner, B. Stelzer, & S. A. Metchev, IAU Symp., 314, 213

Chilcote, J., Pueyo, L., De Rosa, R. J., et al. 2017, AJ, 153, 182

Currie, T., Guyon, O., Tamura, M., et al. 2017, ApJ, 836, L15

de Kok, R. J., Birkby, J., Brogi, M., et al. 2014, A&A, 561, A150

Ducati, J. R. 2002, VizieR Online Data Catalog: II/237

Eisenhauer, F., Abuter, R., Bickert, K., et al. 2003, in Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, Proc. SPIE, 4841, 1548

Fitzgerald, M. P., & Graham, J. R. 2006, ApJ, 637, 541

George, E. M., Gräff, D., Feuchtruber, H., et al. 2016, in Ground-based and Airborne Instrumentation for Astronomy VI, Proc. SPIE, 9908, 99080G

Gontcharov, G. A. 2006, Astron. Lett., 32, 759

Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006, AJ, 132, 161

Guyon, O. 2011, Eur. Phys. J. Web Conf., 16, 03001

Heng, K., & Lyons, J. R. 2016, ApJ, 817, 149

Heng, K., & Tsai, S.-M. 2016, ApJ, 829, 104

Hubeny, I., & Burrows, A. 2007, ApJ, 669, 1248

Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Jones, A., Noll, S., Kausch, W., Szyszka, C., & Kimeswenger, S. 2013, A&A, 560, A91

Jovanovic, N., Martinache, F., Guyon, O., et al. 2015, PASP, 127, 890

Kawahara, H., Murakami, N., Matsuo, T., & Kotani, T. 2014, ApJS, 212, 27

Kenworthy, M. 2017, Nat. Astron., 1, 0099

Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Science, 339, 1398
Kuntschner, H., Jochum, L., Amico, P., et al. 2014, in Ground-based and Airborne Instrumentation for Astronomy V, Proc. SPIE, 9147, 91471U

Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660, 770

Lagrange, A.-M., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21

Lagrange, A.-M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57

Lazzaro, D., Sicardy, B., Roques, F., & Greenberg, R. 1994, Icarus, 108, 59

Lecavelier des Etangs, A., Scholl, H., Roques, F., Sicardy, B., & Vidal-Madjar, A. 1996, Icarus, 123, 168

Lecavelier des Etangs, A., & Vidal-Madjar, A. 2016, A&A, 588, A60

Lovis, C., Snellen, I., Mouillet, D., et al. 2017, A&A, 599, A16

Luger, R., Lustig-Yaeger, J., Fleming, D. P., et al. 2017, ApJ, 837, 63

Macintosh, B., Graham, J., Palmer, D., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Proc. SPIE, 6272, 62720L

Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proc. Natl. Acad. Sci., 111, 12661

Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556

Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348

Mawet, D., Pueyo, L., Lawson, P., et al. 2012, in Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, Proc. SPIE, 8442, 844204

Mesa, D., Zurlo, A., Milli, J., et al. 2017, MNRAS, 466, L118

Moses, J. I., Marley, M. S., Zahnle, K., et al. 2016, ApJ, 829, 66

Mouillet, D., Larwood, J. D., Papaloizou, J. C. B., & Lagrange, A. M. 1997, MNRAS, 292, 896

Noll, S., Kausch, W., Barden, M., et al. 2012, A&A, 543, A92

Rameau, J., Chauvin, G., Lagrange, A.-M., et al. 2015, A&A, 581, A80

Riaud, P., & Schneider, J. 2007, A&A, 469, 355

Roques, F., Scholl, H., Sicardy, B., & Smith, B. A. 1994, Icarus, 108, 37

Royer, F., Zorec, J., & Gómez, A. E. 2007, A&A, 463, 671

Ruffio, J.-B., Macintosh, B., Wang, J. J., et al. 2017, ApJ, 842, 14

Scholl, H., Roques, F., & Sicardy, B. 1993, Celest. Mech. Dyn. Astron., 56, 381

Smith, B. A., & Terrile, R. J. 1984, Science, 226, 1421

Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Nature, 509, 63

Snellen, I., de Kok, R., Birkby, J. L., et al. 2015, A&A, 576, A59

Sparks, W. B., & Ford, H. C. 2002, ApJ, 578, 543

Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466

Thatte, N. A., Clarke, F., Bryson, I., et al. 2014, in Ground-based and Airborne Instrumentation for Astronomy V, Proc. SPIE, 9147, 914725

Todorov, K. O., Line, M. R., Pineda, J. E., et al. 2016, ApJ, 823, 14

van Leeuwen, F. 2007, A&A, 474, 653

Wang, J. J., Graham, J. R., Pueyo, L., et al. 2016, AJ, 152, 97

Wang, J., Mawet, D., Ruane, G., Hu, R., & Benneke, B. 2017, AJ, 153, 183

Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41