© F. Cantalloube et al. 2020

1. Introduction

The latest generation of instruments dedicated to exoplanet and circumstellar disk imaging have in the past five years enabled a huge step in high-contrast imaging (HCI) of the close environment of nearby stars. To detect the light emitted by young Jupiter-like companions in the near-infrared that orbit at a few astronomical units (au) from their host star, which itself is located at a few tens of parsecs from Earth, a contrast better than 10⁻⁵ at an angular separation of 500 milliarcsec (mas) from the star is required. By equipping 8 m class telescopes with dedicated instruments combining extreme adaptive optics (AO) systems using high-density deformable mirrors (DM) with specific coronagraphs and advanced post-processing techniques, instruments such as VLT/SPHERE (Beuzit et al. 2019), Gemini/GPI (Macintosh et al. 2008), and Subaru/SCExAO (Jovanovic et al. 2015) successfully addressed this challenge. However, after achieving such high resolution and contrast, new limitations are now detected in the focal-plane images that were not visible with the first generation of HCI instruments such as VLT/NaCo (Rousset et al. 2003), Gemini/NICI (Artigau et al. 2008), or Keck/NIRC2 (McLean & Chaffee 2000).

The scientific region of interest is the close vicinity of the star (below 500 mas). The detection of exoplanets in this region is crucial to reject or confirm planet formation scenarios because most giant planets are expected to be found in this region (Chauvin 2018; Nielsen et al. 2019). In this region, circumstellar disks are sometimes expected from the infrared excess of the host star. With the latest generation of HCI instruments, the main limitations that particularly affect these inner regions by provoking leakages of the starlight from the focal-plane mask element of the coronagraph are (1) the quasi-statics noncommon path aberrations (NCPA, Guyon et al. 2005; Fusco et al. 2006), which are differential aberrations between the AO arm and the science arm that are either not seen by the AO or that are corrected for, but are absent from the science arm, (2) the low wind effect (LWE, Sauvage et al. 2016; Milli et al. 2018), inducing differential piston and tip-tilt errors between the fragments of the pupil, (3) the low-order residuals (LOR), such as residual tip-tilt, which can be either due to atmospheric residuals, mechanical low-frequency vibrations (about 10 Hz) induced by the telescope pointing (Lozi et al. 2018), or atmospheric dispersion residuals (Pathak et al. 2016), and (4) the wind-driven halo (WDH, Cantalloube et al. 2018), which arises when the atmospheric turbulence evolution is faster than the AO correction timescale (shown in Fig. 1 and described in this paper). Current post-processing techniques fail to overcome these limitations, and the final contrast can decrease by a factor 20 at separations between 200 and 500 mas (e.g., Milli et al. 2018, for the LWE). For more details about these various contributions, Cantalloube et al. (2019) presented a review of the contrast limitations observed in the VLT/SPHERE images, and Mouillet et al. (2018) presented a review of the effect of the AO performance on the SPHERE images.

Fig. 1. Coronagraphic focal plane images showing the wind-driven halo. Left: simulation of a perfect post-AO coronagraphic image of infinite exposure using an analytic AO tool (accounting only for fitting and servolag errors). Right: one exposure obtained with SPHERE-IRDIS (H2 band). Both images are in logarithmic scale to emphasize the WDH. The two regions encircled with yellow dotted lines are artifacts caused by the manufacturing technique of the DM.

The wind-driven halo originates from one of the AO error terms, namely the AO servolag (or temporal bandwidth) error, which is due to the finite and time-delayed nature of the AO correction. Astronomical AO-systems run in closed loop so that the wavefront sensor (WFS) sees the residual phase after the AO correction, and therefore the command sent to the DM is relative to the previous correction. Because there is some time delay between the WFS measurement and the setting up of the DM, the AO servolag error becomes consequential if the atmospheric turbulence has varied significantly between the measure and the applied correction. For a fixed AO delay, the AO servolag error therefore depends on the turbulence coherence time τ₀, itself dependent upon the seeing and the effective wind velocity at the telescope pupil. As a consequence of this servolag error, the AO-corrected phase shows strong low-order residuals along the effective wind direction, which result in a shattering of the point spread function (PSF) along the effective wind direction. For a long exposure time (from 1 to 60 s for observations using SPHERE), these AO residual speckles add up to form a smooth halo, the wind-driven halo, at a contrast typically below 10⁻³. In the context of HCI, when a coronagraph is used to suppress the coherent peak of the starlight, a raw contrast of about 10⁻³ is reached, and this feature becomes visible. Moreover, we recently found that the temporally delayed AO residual phase interferes with amplitude errors, creating an asymmetry of the WDH in its radial direction (Cantalloube et al. 2018; Madurowicz et al. 2019): the stronger the correlation of amplitude error and delayed phase error, the weaker the asymmetry. Figure 1 shows this wind-driven halo contribution in a coronagraphic image from the VLT/SPHERE-IRDIS instrument (Dohlen et al. 2008) for a simulation (left, infinite exposure with a perfect coronagraph) and an on-sky image (right).

Because a significant part of the images obtained with SPHERE are affected by the WDH, we would ultimately like to reconstruct this effect to apply a correction to existing data that have been acquired during the five years of SPHERE operations. For a point source detection, we can apply a spatial high-pass filter to the images, but for disk imaging, this removes most of the object information because the disk signals are mainly spread within low spatial frequencies. We therefore decided to characterize finely this WDH signature in the view of developing more specific post-processing techniques to remove the WDH from the images. In addition, to prepare for the future generation of high-contrast instruments and to optimize the operation of such instruments, this paper presents a complete review of the parameters at stake, in terms of turbulence profiling, AO control and post-processing techniques.

In the following, we first review the physical origin of the WDH to highlight on which parameters it depends and show its effect on the raw contrast (Sect. 2). We then propose a method and metrics for analyzing its contribution directly from the focal plane images (Sect. 3). We then apply this procedure to on-sky SPHERE images to highlight the effect of the WDH on the contrast after post-processing by studying its typical spatial, temporal, and spectral variations (Sect. 4). From these analyses, we conclude that the current post-processing techniques based on differential imaging are not capable of fully removing the wind-driven halo, and consequently, the contrast performance is decreased by an order of magnitude in the AO-corrected area.

2. Origin and consequences of the wind-driven halo

In the following, we detail the temporal aspect of the AO loop (Sect. 2.1) and of the atmospheric turbulence (Sect. 2.2). We also specify how it affects the PSF (Sect. 2.3) and the raw contrast in the specific case of coronagraphic imaging (Sect. 2.4).

2.1. Temporal lag of AO

A classical on-axis single-conjugated AO system is composed of three main components: (i) the wavefront sensor (WFS), which analyzes the incoming phase distortion, (ii) the real-time computer (RTC), which based on the WFS measurement calculates the command that is to be sent to the phase corrector, and (iii) the DM, which corrects for the phase distortion. The AO system of the VLT/SPHERE instrument, SAXO, consists of a 40 × 40 subaperture spatially filtered Shack-Hartmann WFS (Fusco et al. 2006) using an EMCCD detector (Sauvage et al. 2014); the RTC is the ESO-provided SPARTA architecture (Suárez Valles et al. 2012; Petit et al. 2014); and the correction is made in two stages with one tip-tilt mirror and a 41 × 41 actuator high-order DM (HODM, Sinquin et al. 2008). For a full description of the SAXO system, see Fusco et al. (2006).

The different steps occurring during an AO closed-loop run are summarized in the chronogram presented in Fig. 2 (adapted from different schemes from the literature, including Petit et al. 2008). From the first photon reaching the WFS detector to the DM being set, it proceeds as described below.

Fig. 2. Typical AO chronogram summarizing the different sequences of an AO loop using a CCD WFS as on SPHERE-SAXO (arrow lengths are not to scale). One frame is taken every AO loop period (green box), and the AO loop delay (red box) consists of a delay of about 2.2 frames, starting in the middle of the first frame (Tint/2), until the DM is effectively in the corresponding shape.

1. The sequence of events at the WFS follows T_WFS = T_int + T_read, where T_int is a tunable integration time (charge collection in pixel well of the CCD) and T_read is the fixed readout time (depending on the WFS detector technology in use). For SPHERE, T_read = 725 μs, which includes all the operations that are required to complete the image readout. By construction, it is such that T_int ≥ T_read. The minimum interval between the first pixel being integrated of two successive frames is therefore 1/T_read = 1/725 μs = 1380 Hz. This defines the maximum frame rate.

2. The RTC computing time T_RTC is the time the RTC takes to carry out the processing for a given loop cycle. The RTC starts when the first pixel is received from the WFS (this continues in parallel to the readout of the WFS) and ends when the last command is sent to the DM. On SPHERE, it has been measured as T_RTC = 734 μs.

3. The DM settling time T_DM is the time it takes the DM to reach the requested shape after receiving the first command from the RTC. It depends on the rise time of the actuators T_rise, fixed by the DM technology in use. For SPHERE, T_DM is far shorter than all other times involved and can be neglected (ranging from 15 to 20 kHz, Sinquin et al. 2008).

In between these three main parts, there are also fixed transfer times (gathered in T_t) from the WFS to the RTC and from the RTC to the DM. In terms of AO control, the frame rate for the WFS defines the AO-loop frequency (sometimes also referred to as the AO-loop rate).

The AO-loop delay τ_AO is a pure delay defined as the addition of the equivalent delays (to a first-order approximation) from the various processes involved between taking a measurement of the atmospheric disturbance via the WFS and commanding the DM accordingly: the WFS delay (t_WFS), the RTC delay (t_RTC), the digital-to-analog conversion delay at the DM amplifier (t_DAC), the DM positioning delay (t_DM), and at last an overall data transfer delay (t_com). At low running frequency, each term can be approximated as follows:

– t_WFS is approximated by the sum of the readout time T_read and half the integration time T_int/2;

– t_RTC is the RTC latency, the time between the reception of the last pixel from the WFS to the last DM command data sent, measured in the laboratory as t_RTC = 80 μs;

– t_DAC is approximated by T_int/2 (when it is assumed at first order that the response to an impulse is a pulse of duration T_int);

– t_DM is approximated as half the rise-time of the DM, T_rise/2;

– t_com is the overall data transfer delay and has been fit empirically on SPHERE by measuring the closed-loop transfer function to be t_DM + t_com = 35 μs.

For SAXO, the main contributors are therefore the integration time T_int, and then the RTC latency t_RTC. In this framework, we consider the WFS as a low-pass filter that takes an average of the atmosphere during the measurement time T_int. As a whole, the AO-loop delay for SPHERE is 1.56 ms, corresponding to about 2.2 loop cycles when the AO is running at 1380 Hz. In practice, according to the latest tests performed on SPHERE (in December 2018), the measurement of the frame number n is mainly affected by the command number n − 2 (by 84%) and is less strongly affected by the command number n − 3 (by 16%).

As a consequence, as soon as the atmospheric turbulence evolves faster than 2.17 ms (corresponding to three frames at 1380 Hz), the AO servolag error appears in SPHERE, which affects the starlight distribution in the focal plane. On fainter stars, SPHERE SAXO runs slower, therefore this statement is only strictly correct for bright stars with SAXO running at the fastest or optimal loop frequency. The tradeoff between the AO-loop delay and the temporal evolution of the atmospheric turbulence is the critical parameter that causes the WDH in high-contrast images. In the following, we use this temporal description of SAXO to simulate the AO residual phases that are used to produce the SPHERE-like simulated coronagraphic images and to discuss the effect of the atmospheric turbulence evolution.

2.2. Temporal variation in atmospheric turbulence

At a given instant, the atmospheric turbulence state can be represented as a two-dimensional phase screen reaching the telescope aperture. This phase screen can be described by its power spectral density (PSD) in models such as the widely used stratified von Karman model (Conan 2000). This model is parameterized by the Fried parameter (r₀, the typical spatial extension of the turbulence cells) and the outer scale (L₀, the largest size of the turbulence cells). During the observation sequence, this phase screen evolves in two ways: by translation (the flow), and by the evolution in the shape of the turbulent cells and spatial distribution (the boiling). In a multilayer description, the flow is associated with the variations in wind speed and direction of each layer and therefore mainly affects low spatial frequencies variations. The boiling is associated with a change in the mixing of the different layers and therefore affects high spatial frequencies variations. It has been empirically shown that the phase screen autocorrelation decays linearly with time over typical timescales longer than 25 ms (see Guesalaga et al. 2014; Poyneer & Macintosh 2006; Schöck & Spillar 2000, for three different sites, Cerro Pachon, Maunakea, and Albuquerque). Because of the number of subapertures of the WFS per telescope diameter (20 cm sampling for SAXO) and the AO-loop delay (1.56 ms for SPHERE-SAXO), the effect of boiling can be ignored for SPHERE-like instruments. In this case, we can work under the frozen-flow assumption (the so-called Taylor hypothesis, Taylor 1938; Bharmal 2015), stating that the temporal evolution of the phase screen is largely dominated by translation following the projected wind speed and direction (the so-called effective wind velocity). Moreover, boiling will cause an isotropic starlight leakage in the coronagraphic image, which is therefore not linked to the wind driven halo.

Under the frozen-flow hypothesis, the temporal variation of the atmospheric turbulence is parameterized by the turbulence coherence time τ₀, which is analytically defined as (Roddier 1981; Roddier et al. 1982; Hardy 1998)

$$\begin{array}{c}\hfill {\tau }_0=0.314\frac{r_0}{v_{\mathrm{eff}}},\end{array}$$(1)

where r₀ is the Fried parameter (dependent upon the observation wavelength and the zenith angle) and v_eff is the effective wind velocity, itself defined as ${\mathit{v}}_{\mathrm{eff}}={\left[\frac{{\int }_{\infty }{C}_n^2(h).\mathit{v}{(h)}^{5/3}\mathrm{d}h}{{\int }_{\infty }{C}_n^2(h)\mathrm{d}h}\right]}^{3/5}$, with v(h) the wind velocity profile with altitude h, and $C_n^2(h)$ the refractive index structure constant profile with altitude. The turbulence coherence time characterizes the time interval for which the temporal fluctuations of the turbulent phase are equal to 1 rad². When τ₀ equals τ_AO (at the sensing wavelength), it means that between the measurement of the incoming phase and the DM reaching the requested shape, the actual phase evolved by 1 rad², that is to say, its Strehl ratio decreased by 63% because of the servolag error (under the short-phase approximation, which is valid during AO correction). The cumulative histogram of the coherence time values over Paranal observatory (at zenith and at 500 nm) during 3 years of Multi-Aperture Scintillation Sensor–Differential Image Motion Monitor (MASS–DIMM, Kornilov et al. 2007) measurements is presented in Fig. 3 and shows a steep curve at short τ₀ (below the median of 4.5 ms). The distribution of the coherence time as a function of the seeing is shown in Fig. A.1, based on the MASS–DIMM measurements. It is not possible to obtain a direct tradeoff value that compares the AO delay and the turbulence coherence time to state when the WDH appears in the images. In the following, we establish a rule of thumb, based on simulations, to estimate a typical τ₀ value below which the WDH dominates in the image. In a next paper, we will apply our WDH analysis procedure to the SPHERE-SHINE guaranteed-time survey (Chauvin et al. 2017) data to extract a realistic occurrence rate of the WDH.

Fig. 3. Cumulative histogram of the turbulence coherence time τ0 value at the zenith and at 500 nm over 3 years of MASS–DIMM measurements (from 2016 to 2019). The median τ0 value is 4.5 ms (black dashed lines).

For a given AO system, the important external parameters playing a role in the servolag error expression are therefore the seeing and the effective wind velocity, that is to say, the balance between strong turbulent layers and high wind speed layers. As shown in the $C_n^2(h)$ profile, extracted from the Stereo-SCIDAR (SCIntillation Detection And Ranging, Vernin & Roddier 1973; Shepherd et al. 2013) measurements 2018 campaign (Osborn et al. 2018) at Paranal observatory, presented in Fig. 4 (top), the strongest turbulence layers are close to the ground layer. As shown in the wind speed profile v(h) presented in Fig. 4 (bottom), the fastest wind occurs at the jet stream layer, located at about 12 km above sea level (200 mbar), whose wind speed can reach 50 m s⁻¹, depending on the season and the latitude. In the southern hemisphere, the subtropical jet stream showing up from March to December varies from about 20 m s⁻¹ during summer (from November to April), up to about 50 m s⁻¹ during winter (from May to October), according to Gallego et al. (2005). The jet stream is less prominent on other astronomical sites such as Maunakea (Sarazin et al. 2003). Notably, before the installation of the MASS, the effective wind velocity used to estimate the coherence time was empirically computed by Sarazin & Tokovinin (2002) as v_eff ≃ max(v_ground, 0.4v_{jet − stream}), showing the effect of the jet stream on the astronomical data (Masciadri et al. 2013). In order to highlight the predominance of the jet stream layer in the effective wind speed and therefore in the WDH, Fig. 5 shows the relative contribution of the ground layer (below 1 km) and of the jet stream layer (between 8 to 14 km) to v_eff, extracted from the Stereo-SCIDAR data. The median contribution of the ground layer to v_eff is about 10%, while the median contribution of the jet stream layer to v_eff is about 40%. In addition, by comparing the contribution from the ground layer and the contribution from the jet stream layer for each individual profile, we observe that for about 80% of the profiles, the jet stream layer has a higher contribution to v_eff. By correlating the observed WDH direction within high-contrast images from GPI, Madurowicz et al. (2018) also showed that the jet stream layer is indeed mainly responsible for the appearance of the WDH in HCI data.

Fig. 4. Turbulence strength and wind speed profiles with altitude measured at Paranal observatory during the 2018 Stereo-SCIDAR campaign. In both figures, the lighter shaded region contains 90% of the data, the darker shaded region contains 50%, and the solid line is the median. Top: integrated refractive index structure constant between two altitude limits, ${\int }_{h_2}^{h_1}{C}_n^2\text{d}h$, as a function of the altitude h. Bottom: wind speed as a function of the altitude h. The blue solid line corresponds to the median in winter (May to October) and the orange solid line to the median in summer (November to April).

Fig. 5. Histogram of the contribution to veff for the ground layer (below 1 km altitude, blue line) and the jet stream layer (8 to 14 km altitude, orange line), extracted from the 2018 Stereo-SCIDAR measurements at Paranal observatory.

Based on the Stereo-SCIDAR measurements, we also computed the typical temporal evolution of the wind speed (Fig. 6, top) and the wind direction (Fig. 6, bottom) at the jet stream layer during two hours, which is the typical time of an observing sequence with SPHERE. The latter is the distribution of the absolute value of the change in wind speed and direction with some time lag across the entire data set. It shows that the wind speed at the jet stream layer rarely remains stable and is more likely to vary by up to typically 5 m s⁻¹ over an hour. The change in wind direction is likely to remain below 10°.

Fig. 6. Typical temporal evolution of the wind speed (top) and wind direction (bottom) within the 8−12 km altitude region, extracted from the 2018 Stereo-SCIDAR campaign at Paranal observatory. In every panel, the lighter shaded region contains 90% of the data, the darker shaded region contains 50%, and the solid line is the median.

As we recalled in the introduction, the asymmetry of the WDH is due to the scintillation. For a given AO delay and turbulence coherence time, the asymmetry of the wind-driven halo increases with the scintillation (Cantalloube et al. 2018). Scintillation is due to the effect of long-distance propagation that transforms phase aberrations into amplitude aberrations. The scintillation is dependent upon the propagation length z, the strength of the atmospheric turbulence $C_n^2(h)$, the diameter and structure of the telescope aperture D_tel, the wavelength of observation λ (for large telescopes of the 8m class, the scintillation is almost achromatic), the exposure time τ_exp, and the wind speed w(z) at and above the tropopause. Under the frozen-flow hypothesis, the scintillation index (variance of the flux fluctuation through the telescope pupil) can be expressed as ${\sigma }_s^2=10.66D_{\mathrm{tel}}^{-4/3}{\tau }_{\exp }^{-1}{\int }_A{C}_n^2(z)\frac{z^2}{w(z)}\mathrm{d}z$, where A is the telescope aperture filter (Kornilov 2011). The temporal variation of the scintillation has been measured with the Multi-Aperture Scintillation Sensor instrument (MASS, Kornilov et al. 2003) measuring the scintillation index during various campaigns on different sites over 8 years, from 2004 to 2012, published in Kornilov et al. (2012). The authors showed that the scintillation index has seasonal variations depending on the latitude of the observatory. The measured scintillation index is higher during the local winter. It has comparable values on all sites but is lower at Maunakea because of its higher altitude and lower turbulence in the upper atmosphere. The temporal behavior of the scintillation is similar in the 11 surveyed sites. This is directly related to the wind speed seasonal variations in the upper atmosphere. In addition, they observed that the power of the scintillation is quite stable over 1 h timescales at Paranal observatory, which is directly related to the jet stream speed variations, which are also shown in Fig. 6 (top). This again indicates the importance of the jet stream layer for the WDH signature. In the following, all these aspects are taken into account to simulate SPHERE-like images.

2.3. Consequences of the AO servolag error in the images

The spatial variance of the AO residual phase in the pupil due to the AO servolag error varies as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{servo}}^2={\left(\frac{{\tau }_{\mathrm{AO}}}{{\tau }_0}\right)}^{5/3}.\end{array}$$(2)

For a single-conjugated AO system, the PSD of the residual phase due to the AO servolag error, showing the distribution of the averaged power of the phase fluctuations over the spatial frequencies k, can be expressed as (Rigaut et al. 1998; Cantalloube et al. 2018)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{PSD}}_{\mathrm{servo}}(k)=& 0.023k^{-11/3}{\displaystyle {\int }_{\infty }{r}_0{(h)}^{-5/3}[{sin}_c^2(k.v(h).T_{\mathrm{int}})}\hfill \\ \hfill & +1-2cos(2\pi .k.v(h).{\tau }_{\mathrm{AO}}){sin}_c(k.v(h).T_{\mathrm{int}})]\hfill \\ \hfill & \times (1-\sin (2\pi .h/h_T)\pm 2\pi .k.T_{\mathrm{int}}.v(h))\mathrm{d}h,\hfill \end{array}\end{array}$$(3)

with h_T = 2/(λk²) the Talbot length (Antichi et al. 2011). The last line of the expression accounts for the asymmetry due to the interaction with the atmospheric amplitude error. In the following, the simulated images (as in Fig. 1, left) are made from residual phase screens produced by an analytical AO simulator (Jolissaint et al. 2006) using this updated expression of the AO servolag error (Eq. (3)). This updated expression accounts for the interaction of the servolag with the scintillation, therefore generating the asymmetry (it includes both the amplitude and the phase of the electric field). For small residual errors (high Strehl ratio), the PSD is an approximation of the smooth structures of the PSF (such as the WDH).

Simulations of the two-dimensional AO residual phase PSD, with and without the servolag error, are shown in Fig. 7 (right and left panel, respectively). The residual phase PSD due to the AO servolag error shows a low spatial frequency structure along the effective wind direction (white arrow in Fig. 7, right) that radially decreases (until the AO correction radius) and shows no power in the direction perpendicular to the effective wind direction. In addition, as a result of the interference term with amplitude error, one wing of the AO servolag signature is smaller than the other (in the opposite direction of the wind).

Fig. 7. PSDs of the AO residual phase simulated for a SAXO-like system. Left: PSD without the servolag error. Right: PSD including the servolag error and its interference with amplitude errors for an effective wind along a direction of 30° (white arrow). In order to highlight the servolag contribution, the color bar is normalized to the maximum of the right image.

By comparing the simulated PSDs with and without the servolag error, we can conclude that for a SAXO-like AO system working under the typical median τ₀ of Paranal observatory, about 69% of the starlight within the AO corrected zone is scattered outside of the coherent peak due to the servolag error. The remainder of the scattered light is caused by other typical AO errors (aliasing, chromaticity, anisoplanetism, and noise propagation, but excluding NCPA or other exogenous errors).

2.4. Effect of the wind-driven halo on the raw contrast

In this section, we analyze how the WDH affects the raw contrast, that is to say, the contrast obtained in an image after the AO correction with a coronagraph, but before the application of any post-processing technique. In practice, it is computed as the mean radial profile of the coronagraphic image normalized by the maximum of the noncoronagraphic image.

To highlight the contribution of the WDH in high-contrast images, we simulated AO-corrected (using a SAXO-like system as described in Sect. 2.1) and ideal coronagraphic (Cavarroc et al. 2006; Sauvage et al. 2010) infinite-exposure images in H band (1.6 μm), producing Fig. 1 (left). We used a single-layer atmospheric model, moving in one given direction, varying only the wind velocity, under the median seeing conditions at Paranal observatory. To assess the effect of the WDH on the raw contrast, Fig. 8 shows the radial profiles of simulated coronagraphic images for various τ₀ in the wind direction (solid lines, where the raw contrast is highly affected) and in its perpendicular direction (dashed lines, where the raw contrast is less affected), under median seeing condition (r₀ = 13.8 cm, according to MASS–DIMM measurements) and median airmass (a = 1.15).

Fig. 8. Illustration of the raw contrast as a function of separation to the star for various τ0, ranging from 12.5 ms (black) to 1.25 ms (red). From bottom to top, the curves are for turbulence coherence times of 12.5, 10, 7.5, 5.0, 4.0, 3.0, 2.5, 2.0, 1.5, and 1.25 ms. The solid lines show this along the wind direction, and the dashed lines show it along the perpendicular direction. For comparison, the gray dash-dotted line shows the raw contrast (azimuthal mean of the image) of a SPHERE-IRDIS H2-band (centered on 1.59 μm) image under median seeing, airmass, and long τ0 ∼ 9 ms.

We compared these simulations with the raw contrast of a SPHERE image taken under very good observing conditions (median seeing, τ₀ >  9 ms), in which no WDH is visible (Fig. 8, gray dash-dotted line and Fig. 9, top right). At a separation of 300 mas (where the raw contrast without servolag error reaches a plateau, and far beyond the influence of the inner working angle of the coronagraph), the raw contrast in the H band is about 7 × 10⁻⁵ (see also Vigan et al. 2015). In this on-sky image, the raw contrast is limited by the presence of speckles that are due to NCPA. As NCPAs are always present and are the main limitation at 300 mas, this is the ultimate raw contrast we can reach under good observing conditions (Vigan et al. 2019). In Fig. 8, the measured contrast curve (gray dash-dotted line), which was taken with a coherence time of 9 ms, is higher than expected for the system with a 3 ms coherence time (green line). We therefore expect the WDH to be detected from a coherence time below 3 ms. In the cumulative histogram in Fig. 3, this yields an occurrence rate of WDH of about 30% after correcting for the median airmass of 1.13 (corresponding to a zenith angle of 27.7°) as measured for SPHERE. This value is an approximate value to motivate the present work and will be further refined by a statistical study of the SPHERE data acquired during its 5 years of operation.

Fig. 9. On-sky images from the SPHERE-IRDIS instrument in H2 band without (top) and with WDH (bottom). Left: noncoronagraphic images. Right: corresponding coronagraphic images. The top images are taken with a very long coherence time τ0 ∼ 9 ms.

The WDH has a high contrast and is only unveiled through the use of a coronagraph, as highlighted in Fig. 9, where no difference can be seen in the two noncoronagraphic images (left), while the WDH is very bright in the coronagraphic image for a short coherence time (bottom right). The corresponding raw contrast profiles are shown in Fig. 10, showing the effect of the WDH on real data. This effect therefore appears to be a strong limitation only for the latest generation of HCI instrument such as VLT/SPHERE, Gemini/GPI, and Subaru/SCExAO, which use both extreme-AO and advanced coronagraph technology. The careful analyses of the WDH presented in this paper are therefore motivated by the strong effect of the WDH on the achieved contrast with these high-contrast instruments.

Fig. 10. Raw contrast as a function of the separation to the star for the two images with (solid lines) and without (dot-dashed line) WDH presented in Fig. 9. The red line shows this along the WDH, the blue line in the perpendicular direction, and the black lines are the azimuthal median. The two red shaded areas correspond to the region that is affected by the coronagraph (left) and the region outside the AO correction zone (right).

3. Analysis of the wind-driven halo in the focal-plane images of SPHERE

In this section, we present a method for analyzing the WDH by deriving the three properties used to describe it: its direction (Sect. 3.2), its intensity in the focal-plane image (Sect. 3.3), and its asymmetry (Sect. 3.4). The analyses are conducted directly in the focal-plane images because the results are more reliable in the context of high-contrast performance than if we were to use external data such as AO telemetry or turbulence profiling. In addition to the analysis presented in this paper on the WDH effects on the contrast after post-processing, we intend to use this WDH analysis procedure for two studies: (i) a statistical analysis of the effect on SPHERE data caused by the WDH and a correlation of this effect with the AO telemetry and profiling data, and (ii) an estimation of the WDH to remove it from the data.

To illustrate and verify our approach, we used four types of multispectral coronagraphic images, following the integral field spectrograph (IFS, Antichi et al. 2009) of SPHERE in the YH band (from 0.96 to 1.66 μm, with a spectral resolution R ∼ 30). We list these types below.

– Case 1: Simulated images obtained as a temporal stack of short exposures, containing only the AO residuals that are due to fitting and servolag errors (produced using our analytical AO simulator) and using an apodized Lyot coronagraph (APLC, Soummer et al. 2011a; Martinez et al. 2009) such as the one used on SPHERE (see Fig. 11, left).

Fig. 11. High-contrast images (at the shortest wavelength, 0.966 μm) used to verify the WDH analysis procedure (logarithmic scale). Left: simulation using only the fitting and AO-servolag errors. Middle left: simulation additionally including NCPA upstream (given by the ZELDA on-sky estimate) and downstream of the APLC focal-plane mask. Middle right: simulation additionally including LOR and LWE. Right: on-sky image of the star 51 Eri taken with the VLT/SPHERE-IFS instrument.

– Case 2: Simulated images like before, additionally including NCPA upstream and downstream the coronagraph focal plane mask. In order to have realistic simulations that are similar to the images obtained with SPHERE, we used the upstream phase estimated by the Zernike wavefront sensor mask (ZELDA, N’Diaye et al. 2013; Vigan et al. 2019) during the latest tests conducted on SPHERE (see Fig. 11, middle left).

– Case 3: Simulated image like before, but additionally including a small amount of LOR (5 mas tip-tilt) and of LWE (5 mas) in random directions for each pupil fragment separated by the spider of the telescope (see Fig. 11, middle right).

– Case 4: On-sky VLT/SPHERE-IFS data cube of the 51 Eri star taken during the SPHERE-SHINE guaranteed-time survey (Chauvin et al. 2017) and described in Samland et al. (2017) and Maire et al. (2019) (see Fig. 11, right).

For the three sets of simulated images, 200 short exposures were stacked to obtain a long-exposure image. The injected wind direction was 125 degrees. For each multispectral cube, Fig. 11 shows only the shortest wavelength of the cube, at 0.966 μm.

3.1. Extracting the WDH contribution from the image

In the focal-plane image, the WDH extends from the center to the edge of the AO-correction zone, as shown in Fig. 7. Excluding the inner working angle of the coronagraph, the WDH is a low spatial frequency feature, whose intensity is dependent upon the temporal variation of the atmospheric turbulence (Sect. 2.2), the AO servolag error (Sect. 2.1), and the image exposure time (DIT). In the case of SPHERE, the AO delay is fixed and the DITs are long enough so that the WDH is observed as a smooth structure. One simple way to separate the WDH contribution from the other starlight residuals that are mainly high-spatial frequency speckles (originating either from residual atmospheric turbulence or from NCPA) therefore is to spatially filter the data in the Fourier domain to retain only their low-frequency content. To perform the filtering, a Hamming window is a good compromise to avoid Gibbs effects (being a continuous function) while being steep enough to separate the frequencies at the user-specified cutoff frequency. The estimated WDH is the low-pass filtered image on which we applied an annular binary mask to cover the coronagraph signature (below 2λ/D) and the seeing-limited area (beyond 20λ/D). Because SPHERE images show two artifacts on their correction ring that are due to the DM square grid (see Fig. 1 right image, encircled with yellow dashed line), these two spots were also masked.

To verify whether the high-pass filtering can indeed extract the WDH contribution, Fig. 12 shows the on-sky 51 Eri image (Fig. 11, right) filtered with different filtering fractions (percentage of low frequencies kept in the image): the top row shows the low spatial frequencies, and the bottom row shows the high spatial frequencies. For a filtering fraction of 5% (left column), the high-pass filtered image (top row) results in a halo that is too smooth, whereas from a filtering fraction of 15% on, it shows sharper structures caused by the capture of speckles. For a filtering fraction of 5%, the low-pass filtered image (bottom row) still shows a slight elongation in the wind speed direction at short separation, which completely disappears from 15% on. When the filtering fraction is too high, very faint structures such as the grid of microlenses from the IFS design become visible. After testing different data set, a good trade-off based on visualization of the images (the low-pass filtered image shows a very smooth structure while its high-pass version shows no further directional elongation), is to use 15% of the low spatial frequency content. This qualitative argument is confirmed by further analysis.

Fig. 12. Spatially low-pass filtered images of the 51 Eri image shown in Fig. 11 right (top, logarithmic scale) and difference between the image and its low-pass filtered version (bottom, linear scale). The images are masked to show only the AO-corrected area and hide the DM artifacts and the coronagraphic signature. From left to right with a filtering fraction (percentage of low frequencies kept in the image) of 5%, 10%, 15%, 20%, and 25%.

To determine whether this filtering procedure indeed shows mainly, if not only, the WDH component, we applied it on the three simulated data cases and compared this with the exact same simulated cases, but without the servolag error. By comparing the images produced with and without WDH, we found that 73% (66% and 64.5%) of the light in the corrected area belongs to the WDH for the first case (second and third, respectively). To compare these absolute values to those extracted by low-pass filtering the images, we computed the fraction of WDH (the ratio between the total intensity in the filtered masked image and the total intensity in the nonfiltered masked image) as a function of the filtering fraction (Fig. 13). The fractions of WDH extracted in the three cases are indeed lower than the absolute values (shown as horizontal lines in Fig. 13), but beginning at a filtering fraction of 22%, it reaches a plateau at the expected absolute values. When NCPA and LWE or LOR are added, the starlight in the corrected area is scattered in other higher spatial frequencies that are not captured by the low-pass filtering, hence the lower fraction of WDH for cases 2 and 3. As a conclusion, the WDH contribution can be extracted from the focal-plane images by applying a low-pass filter with a filtering fraction of about 15%. This rule of thumb is only valid for SPHERE data, and the optimum fraction must be determined for a given instrument.

Fig. 13. Amount of light within the WDH structure that is extracted by spatial filtering as a function of the filtering fraction (starting from 5% at the vertical dotted black line). The red solid line is the first case containing only the WDH contribution, the orange dashed line is the second case with NCPA, and the green dot-dashed line is the third case containing NCPA, LWE, and LOR. The horizontal lines, with the same color code, are the real theoretical values extracted by simulating the exact same images, but without the servolag error ( 73%, 66%, and 64%, respectively).

3.2. Direction of the WDH

In order to assess the direction of the WDH elongation, we developed the following procedure. In a first step, we applied a discrete Radon transform¹ (Radon 1917, 1986) to the estimate of the WDH, which provides the integrated intensity over one direction as a function of the angle. In a second step, to obtain the profile of the intensity as a function of the angle R_WDH(θ), we averaged a few channels (typically a few pixels) of the Radon transform around its center (corresponding to the center of the image). In a third step, we performed a Gaussian fit around the maximum value of this Radon profile in order to extract the preferential direction of the WDH. Using the maximum value of the Radon profile is not the most robust option becaue depending on the number of pixels within the masked image, the profile might be irregular. After testing different possibilities and data sets, we found that performing a Gaussian fit around the maximum value of the Radon profile is the most reliable method that does not require tuning any user-parameter. The uncertainties on the estimated direction are therefore the uncertainties on the Gaussian fit performed because in an ideal case, the Radon profile is purely sinusoidal and the Gaussian fit around the maximum value is quite accurate. With this procedure, we extract the preferential direction of the starlight distribution in the field of view.

In a first step, we checked that this approach is valid for extracting the direction of the WDH on the first simulated data set. This data set only contains the WDH contribution, therefore we do no need to perform the low-pass filtering. For this case, the Radon profile (Fig. 14, top left) is quite smooth, and the maximum of its Gaussian fit yields an estimated direction of 125.13 degrees (see Fig. 14, top right). The uncertainty on the Gaussian fit is negligible, which means that systematic errors (such as centering of the raw image and plate scale) dominate the error. In the following we therefore ignore the uncertainty on the estimated direction.

Fig. 14. Estimate the wind direction for the four cases: extracted radon profile of the WDH and Gaussian fit (red dashed line) whose maximum (red dot-dashed line) shows the estimated preferential direction (left), and contour plot of the WDH showing the fitted wind direction (red cross) with the described procedure (right). From top to bottom: Case 1 (from AO simulations including only the fitting and servolag errors), case 2 (adding NCPA), case 3 (adding LOR and LWE), and case 4 (51 Eri on-sky data). Except for case 1, a filtering fraction of 15% has been used to isolate the WDH. For the first three cases we simulated, the green line shows the simulated wind direction (125 degrees).

As a second step, we applied this procedure on the second and third simulated data set. After low-pass filtering the images with a filtering fraction of 15%, Fig. 14 (middle panels) shows the obtained Radon profiles (left) and the extracted WDH contour plots with the fit direction overlaid (right). For simulations including NCPA (case 2), the estimated direction is 125.1 degrees, as in the case without NCPA. The NCPAs mainly induce high spatial frequencies in the focal plane (see Vigan et al. 2019), they do not greatly affect the extraction of the WDH and therefore its direction. In simulations that additionally include LOR and LWE (case 3), the estimated direction is 124.5-pagination degrees. These low spatial frequencies slightly offset the WDH, as shown in the Radon profile, which is not perfectly sinusoidal (Fig. 14, middle left).

As a third step, we applied this approach to the on-sky data of 51 Eri for which the centering of the star behind the coronagraph is not perfect and some features due to LWE are visible (Fig. 11, right). As expected, the obtained Radon profile deviates even more from a sinusoidal shape (Fig. 14, bottom left). The contour plot shows that the estimated direction is visually in line with the average halo (Fig. 14, bottom right). By construction, our approach does not fit the inner or outer part perfectly, but it is a good fit overall.

As a last step, we verified that the estimated direction is reliable with the filtering fraction we used. For the three simulated data sets, Fig. 15 shows the fitted angle of the WDH direction as a function of the filtering fraction. From a filtering fraction of 15%, the estimated direction reaches a plateau around the injected value. In this specific case, the estimate varies from the injected value by only two degrees even without filtering. We repeated this experiment with on-sky data (for which the real wind direction is not known), which showed a different amount of WDH. Except for the case of a very low WDH contribution, the estimated direction is stable to 1 degree for a filtering fraction between 5% and 25%. For very low WDH, the estimate is dominated by residual tip-tilt errors and varies by up to 7 degrees between a filtering fraction of 5% to 17%, but at 17%, a plateau is reached. This experiment can be used to determined the optimal filtering fraction for a given HCI instrument. In the context of this study and our two future applications, we would like to estimate the wind direction to an accuracy of a few degrees. The procedure for extracting the direction of the WDH is therefore reliable enough to further analyze the WDH structure.

Fig. 15. Fitted WDH direction in the simulated images as a function of the filtering fraction (starting from 5% at the vertical dotted black line). The horizontal solid black line shows the injected value (125 degrees). The red solid line shows the first case that only contains the WDH contribution, the orange dashed line represents the second case with NCPA, and the green dot-dashed line is the third case, which contains NCPA, LWE, and LOR.

3.3. Strength of the WDH

When the image contains no WDH, the procedure for extracting the direction provides a random direction. To determine whether any WDH is present in the data after the filtering and extraction of the Radon profile, we computed the standard deviation of the Radon profile. If this standard deviation is lower than or equal to the contrast expected without WDH, then there is no WDH. This threshold C_τ is obtained empirically through the simulated PSD of the residual phase without servolag from its minimum in the region where the profile of the Radon transform is computed (for the SPHERE data in H band, this threshold was set to C_τ = 7 × 10⁻⁴). We ran this procedure on SPHERE on-sky data, and it efficiently sorted the data without WDH from the data with WDH.

More generally, in order to assess the amount of starlight that is scattered into the WDH within the corrected area, we defined its strength, 𝒮_WDH, which is 100% if all the starlight in the AO-corrected zone belongs to the WDH and 0% if there is no WDH,

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{WDH}}=\frac{\sigma (R_{\mathrm{WDH}}(\theta ))-C_{\tau }}{\sigma (R_{\mathrm{WDH}}(\theta ))+C_{\tau }}\times 100,\end{array}$$(4)

where σ(R_WDH(θ)) is the standard deviation of the Radon profile R_WDH that we extracted previously. The defined strength is not an exact estimate but provides a relative value with respect to the total intensity in the corrected area and serves as an indicator.

When this metric is applied to the three simulated data cases containing the same amount of WDH, the estimated strength as a function of the filtering fraction is about the same for all cases (to within 0.5% of each other). At a filtering fraction of 10%, this metric reaches a plateau at a strength of 95%. For simulated data that contain WDH alone, this metric provides a strength of 99.6%, and for simulated data without WDH, it provides a strength of 0.05%.

To verify that the defined metric is suitable for real data, we sorted out six images within the on-sky data cube of 51 Eri that show more or less WDH (Fig. 16, bottom). The extracted strength of the WDH as a function of the filtering fraction is shown in Fig. 16 (top), and it ranges from 46% in the weakest case of WDH (Fig. 16, bottom left) to 80% in the strongest case (Fig. 16, bottom right). The variation of the strength we obtained for the six cases is consistent with what is observed in the image.

Fig. 16. Evaluation of the strength of the WDH. Top: strength of the WDH contribution in the AO-corrected area of the focal-plane image 𝒮WDH, defined in Eq. (4), as a function of the filtering fraction (starting at 5%, black dotted line). Bottom: corresponding images from the 51 Eri data set showing low WDH residuals (left), and average (middle) to high (right) WDH. The images corresponding to the green and blue lines look very similar, therefore we only show one here.

3.4. Asymmetry of the WDH

To quantitatively characterize the asymmetry of the WDH, we defined its asymmetry factor, 𝒜_WDH, which is 0% when the WDH is perfectly symmetric and 100% when the WDH shows only one wing,

$$\begin{array}{c}\hfill {\mathcal{A}}_{\mathrm{WDH}}=\frac{\int (\overline{I^{+}}(r,\theta )-\overline{I^{-}}(r,\theta ))\times M}{\int (\overline{I}(r,\theta )\times M)}\times 100,\end{array}$$(5)

where $\overline{I}(r,\theta )=\overline{I^{+}}(r,\theta )+\overline{I^{-}}(r,\theta )$ is the total intensity contained in the WDH. The two wings I⁺ and I⁻ are obtained by cutting the image along the direction perpendicular to the WDH direction estimated previously, I⁺ being the most intense one.

For the three simulated cases, we obtain an asymmetry factor of 20%±2% regardless of the filtering fraction. To verify that this metric is suitable, we simulated images in exactly the same conditions as for the third case (including NCPA, LOR, and LWE), but with a varying amount of asymmetry (Fig. 17, bottom). We compared the extracted values from the images with the values obtained by analyzing the PSD of the AO residuals we used to produce the simulated images (solid dashed lines in Fig. 17, top). In the case without asymmetry, the extracted value is indeed close to 0% and the asymmetry factor gradually increases with the injected amount of asymmetry. In any of these cases, at a filtering fraction of 15% the extracted asymmetry factor is fully stable. When the WDH is strong enough, the extracted value reaches exactly the value expected from the PSD. When the WDH is fainter, the extracted value is a few percent lower than the value expected from the PSD.

Fig. 17. Evaluation of the asymmetry of the WDH. Top: asymmetry of the WDH contribution in the AO-corrected area of the focal-plane image 𝒜WDH, defined at Eq. (5), as a function of the filtering fraction (starting at 5%, black dotted line). The blue solid line shows the image when it is simulated without asymmetry. Bottom: corresponding simulated images showing no (left) to high asymmetry of the WDH (right). The dashed lines correspond to the value that we directly extracted from the PSD.

We studied how our estimate of the direction is affected by the asymmetry of the WDH. Owing to the method itself, the estimated direction is more accurate, with a small amount of asymmetry, but the error remains low (in the case with strong asymmetry, the error is smaller than 10 degrees).

This section confirms that the described procedure is valid for directly extracting the WDH parameters from the focal-plane image. In the following, we use the three metrics defined in this section to characterize the WDH (direction, strength, and asymmetry), with a filtering fraction of 15% (providing stable results), in order to analyze the spatial, temporal, and spectral behavior of the WDH.

4. Effect of the WDH on the final contrast after post-processing

Within the coronagraphic image sequence delivered by the latest generation of HCI instruments, the contrast level of the starlight residuals in the corrected area ranges from 10⁻³ to 10⁻⁵. To achieve a higher contrast, advanced post-processing techniques have been developed that rely on observing strategies that provide a diversity between the starlight residuals and the potential circumstellar signals (disks or planets). Under good observing conditions (i.e., long coherence time), the main expected feature was quasi-statics speckles (QSS) that are neither stable enough to be calibrated nor fast enough to be smoothly averaged and removed by an appropriate filtering. These QSS originate from NCPA, and for short exposure times, from atmospheric residuals (e.g., Cantalloube et al. 2019). Consequently, the post-processing techniques that have been developed aim at removing these QSS, but are not tailored for the WDH, the LWE, or other sources of errors whose temporal, spectral, or spatial behavior differs from the QSS. As a consequence, the contrast reached after post-processing is usually lower than expected from simulations. After 5 years of SPHERE operations, the median contrast reached at 500 mas is 2 × 10⁻⁵ (Langlois et al., in prep.) instead of the 10⁻⁶ expected from simulations before instrument commissioning (Vigan et al. 2010).

The mainstream post-processing techniques used today rely on differential imaging that consists of (1) estimating the QSS, (2) subtracting it from the image, and (3) combining all the subtracted images to increase the signal-to-noise-ratio of potential companions. In order to estimate the QSS, most HCI instruments working in the near-infrared use the temporal diversity by carrying out the observations in pupil-tracking mode to spatially fix the pupil (the instrumental aberrations remain static with time) while the field of view (hence the circumstellar signals) rotates at a deterministic velocity: this is the angular differential imaging (ADI, Marois et al. 2006a) technique. If the instrument additionally provides simultaneous images at two or more wavelengths, the spectral diversity can be used as the QSS expands radially at increasing wavelength, while the circumstellar signals remain at a fixed position: this is the spectral differential imaging (SDI, Racine et al. 1999) technique. Finally, if a sufficient number of images of various targets are available from the instrument, multiple reference differential imaging (MRDI, Lafrenière et al. 2009; Soummer et al. 2011b; Xuan et al. 2018; Ruane et al. 2019, with space-based instruments and with ground-based instruments resp.) can be applied using the correlation between the images to be processed and the image library from the instrument.

For these three solutions, the estimated QSS is usually not perfectly estimated, which yields a high number of speckle residuals, mostly at close separation to the star where the behavior of speckles is nonlinear because of the coronagraphic device (for a classical focal-plane mask coronagraph). In addition, part of the circumstellar signal might be absorbed in the QSS model and removed from the image, yielding a lower signal-to-noise-ratio, and in the specific case of extended sources, a distorted shape or even the signal might be completely removed. See Milli et al. (2012) for the effect of ADI on extended signals and Pairet et al. (2018) for the effect of an ADI-based PCA on extended signal. For point source detection, a high-pass spatial filtering can be performed to remove the light from the WDH at the cost of a slight loss of signal that can be modeled and accounted for to characterize the candidate (Cantalloube et al. 2015). This is not possible for extended sources, however, in which case most of the signal would be removed, in particular for low surface brightness disks seen face on. As a consequence, when the number of disks that are detected in scattered light (about 40) is compared to the number of potentially detectable disks with an infrared excess greater than 10⁻⁴ from Spitzer data (Chen et al. 2014), we find a detection rate of ∼25% (Milli et al., in prep.). It is still unclear whether this is due to the actual disk configuration (too extended or too narrow to be seen by HCI) or due to the limited post-processing techniques available that absorb the disk signals and do not specifically account for other error terms than QSS. In this section, we analyze the temporal, spectral, and spatial variations of the WDH to determine how it affects the post-processed contrast and to understand better how it might be removed by a different approach.

4.1. Temporal variations of the WDH

For the 51 Eri data set used in the previous section, Fig. 18 shows the variation of the direction (left), the strength (middle), and the asymmetry (right) of the WDH as a function of the time between the first image to the next ones. The whole observation sequence lasted 75 min and the reduced cube consists of 64 frames that are each made of four binned images of 16 s integration time (the total number of raw images is therefore 256). Four frames have automatically been rejected by the SPHERE reduction pipeline (Delorme et al. 2017) because of their poor quality.

Fig. 18. Temporal variation of the three WDH parameters in the case of the 51 Eri IFS data taken with SPHERE on September 25, 2015. Left: direction as a function of time compared to the variation of the parallactic angle (blue dashed line). Middle: strength as a function of time. Right: asymmetry factor as a function of time. In the three cases, the gray dash-dotted line defines the moment when the target crossed the meridian (airmass of 1.08).

The direction of the WDH rotates linearly with time and seems to follow the parallactic angle, as expected if the WDH is induced by the upper atmospheric jet stream layer at 12 km with a stable direction. By fitting the slope of the WDH direction as a function of the parallactic angle variation using a robust affine fit (to avoid taking into account outlier data points for which the WDH is too low to extract its direction accurately), we find that the slope is not perfectly equal to one but is rather 1.3. This slight discrepancy is expected if the wind direction at the jet stream layer changes over the timescale of the observing sequence, as shown in Fig. 6 (bottom). The strength varies erratically from 60% to 88%, which corresponds to what we can observe by visualizing the data cube. We compared the temporal variation of the estimated strength of the WDH from this data set with both the SPHERE AO telemetry and the MASS–DIMM turbulence profiler measurements of the turbulence coherence time: as shown in Fig. B.1, the variation of the extracted WDH strength from the SPHERE images is visually consistent with the variation of the coherence time obtained from these two external measurements. The asymmetry factor decreases slowly from 17% to almost no asymmetry, without any link to the strength or the airmass variation (which at first order increases and decreases around the meridian crossing).

For ADI-based algorithm, the important aspect is that the starlight residuals are spatially stable in time to be efficiently removed. We repeated this analysis on various data sets from SPHERE, and the strength of the WDH always varied significantly from one frame to another. This prevents the ADI-based algorithm from correctly capturing the level of this feature in the model of the starlight residuals, resulting in a strong asymmetric WDH residuals in the final post-processed images, as shown in Fig. 19. In addition, the WDH direction follows the parallactic angles, that is to say, the trajectory of an object. This shows that when the temporal data cube is only rotated and median combined, the WDH signature coaligns and adds up. As a consequence, the temporal median of the data cube does not capture the WDH, whose strong signature remains in the combined subtracted images (c-ADI, Fig. 19, left). When a principal component analysis (PCA, Soummer et al. 2012; Amara & Quanz 2012) is applied, the intensity of the WDH residuals depends on the number of principal components that is kept to build the model. As shown in Fig. 19 (right), a more aggressive PCA removes low spatial frequencies but also considerably reduces the signal we seek, mostly for extended features, and still leaves high residuals at close separation to the star. We additionally processed the same data set with locally optimized combination of images² (LOCI, Lafrenière et al. 2007) and non-negative matrix factorization³ (NMF, Ren et al. 2018), which also leave strong WDH residuals in the final processed frame. In all these cases, the companion 51 Eri b (Macintosh et al. 2015) is hardly visible in any of the final processed frames. After cADI, the typical contrast of the WDH residual is 10⁻³ at 300 mas, compared to 5 × 10⁻⁶ when the WDH is filtered out.

Fig. 19. Post-processed image of the 51 Eri SPHERE data cube (as in Fig. 11, left) using ADI-based algorithms. Left: cADI (Marois et al. 2006a) for which the temporal median is used as a QSS field model. Middle left: PCA (Amara & Quanz 2012; Soummer et al. 2012), for which a linear combination of the three first principal components is used as a QSS field model. Middle right: PCA using five components. Right: PCA using ten components. Contrast of about one order of magnitude is lost because of the WDH residuals.

4.2. Spectral variations of the WDH

To investigate the spectral variation of the WDH, we ran our analysis procedure on the three simulated IFS data cubes and on the on-sky data cube of 51 Eri described in Sect. 2 (whose images at the shortest wavelength are shown in Fig. 11). Each cube consists of 39 images at wavelength λ ranging from 0.96 μm to 1.64 μm.

The direction of the WDH is obviously constant with the wavelength, which is indeed verified as shown in Fig. 20 (left). Our method provides a constant direction that varies at most by 2 degrees for all of the four tested cases. As mentioned in Sect. 2, the turbulence coherence time τ₀ according to Eq. (1) scales with the wavelength as r₀, that is to say, as λ^6/5. The variance of the AO residual phase due to the servolag error according to Eq. (2) scales as ${\tau }_0^{-5/3}$. Therefore the WDH strength scales with the wavelength as λ⁻², which is indeed what we observe in Fig. 20 (middle). The asymmetry factor is expected to steadily increase with the wavelength, as demonstrated in Cantalloube et al. (2018). This is indeed what is observed in Fig. 20 (right), where in the four cases, the asymmetry of the WDH follows the wavelength variation.

Fig. 20. Spectral variation of the three WDH parameters in the case of the 51 Eri IFS data taken with SPHERE on September 25, 2015. Left: direction as a function of wavelength. Middle: strength of the WDH as a function of wavelength. Right: asymmetry of the WDH as a function of wavelength. In each plot, the dashed gray line shows the expected trend.

The model for performing SDI is based on two aspects: (i) the diffraction pattern radially extends with the wavelength, and (ii) its intensity evolves with λ⁻². The speckle field at a given wavelength s_λ1(r) can be rescaled to a second speckle field at a different (but close) wavelength s_λ2(r) through the square of the ratio between the two wavelengths: $s_{{\lambda }_2}(\frac{{\lambda }_2}{{\lambda }_1}\mathbf{r})={(\frac{{\lambda }_1}{{\lambda }_2})}^2\times s_{{\lambda }_1}(\mathbf{r})$ (Racine et al. 1999; Pueyo & Kasdin 2007). After rescaling one image (the reference channel), SDI simply subtracts it from the other (the science channel). This widely accepted model is valid assuming that the aberrations are achromatic (i.e., under the Fraunofer approximation, stating that all aberrations originate in a pupil plane), and under the small-phase approximation (Marois et al. 2006b).

At first order, both the QSS and the WDH originate from phase error located in pupil planes, but the asymmetry varies as λ, and is therefore not accounted for on the SDI. As shown in Fig. 21 (left), a simple SDI subtraction leaves a strong imprint on the residual map. On the other hand, in practice SDI is always used in combination with ADI, as shown in Fig. 21 (right), which emphasizes the WDH residuals⁴. This imprint becomes stronger when the reference channel is farther from the science channel. After an SDI subtraction, about one order of magnitude in contrast is lost due to the WDH residual.

Fig. 21. Post-processed images of the 51 Eri data using SDI, the reference is the tenth channel (1.12 μm), which is rescaled and removed from the science channel at 1.00 μm. Left: simple SDI subtraction. Right: classical SADI subtraction, where about half an order of magnitude in contrast is lost because of the WDH residuals.

4.3. Spatial variations of the WDH

As discussed in Sect. 2, the shape, direction, and intensity of the WDH signature depend on various external atmospheric turbulence parameters, on the observed target, and on the observation setup. In the context of RDI, which was mainly developed for extended source extraction, this means that it is unlikely that a set of reference stars of the same spectral type as the science target exhibits a similar WDH signature. As a consequence, we again observe strong WDH residuals in RDI post-processed images, as illustrated in Fig. 22. The images shown were processed using an algorithm that implemented the PCA, as described in Soummer et al. (2012). In ADI (left), the principal components were computed from the science frames themselves, while in RDI (right), the principal components were computed from a library of frames that did not include the science target. This library was built from targets that are part of the SHARDDS program survey for debris disks using SPHERE-IRDIS (Wahhaj et al. 2016; Choquet et al. 2017; Marshall et al. 2018), and the frames were selected based on their correlation with the science target (Milli et al., in prep.). The WDH residuals still remain in the RDI-PCA post-processed images, mostly at a large distance where it is known that RDI reaches a turnover (Ruane et al. 2019) and shows a residual level that is equal to or higher than that of ADI-PCA. The effect of the asymmetry is even more emphasized in this example, where the upper right residual wing is smaller than the other. These considerations depend on the selection within the image library of the reference that is to be subtracted. The highest correlation is usually set on the QSS and not on the WDH. A larger library is also expected to increase the likelihood of obtaining a reference that also contains a similar WDH signature. With a basic RDI implementation, the contrast loss due to the WDH residual signature is about an order of magnitude in the affected regions.

Fig. 22. Residual WDH signature in RDI post-processed images from SHARDDS SPHERE-IRDIS data taken in broadband H (1.625 μm ± 0.15). Left: ADI-PCA post-processed image, using 2 principal components that are subtracted, shown for comparison. Right: RDI-PCA post-processed images, using 10 (left), 25 (middle), and 50 principal components. All frames share the same color scale.

In a next paper we will propose a method for reducing the effect of the WDH on the contrast after differential imaging without affecting close-in and/or extended signals. The idea is to use the analysis procedure described in Sect. 3 to estimate a model of the WDH and subtract it so as to recover a better contrast that is intrinsically limited by the NCPA.

5. Conclusions and perspectives

We conducted a study of the wind-driven halo that is visible in VLT/SPHERE images, a specific signature that significantly affects the contrast performance of the instrument. We provided a detailed examination of the different parameters that play a role in building up this WDH, from the instrument design to its interaction with the observing conditions and operations. When the turbulence coherence time is below 3 ms, the WDH appears in SPHERE coronagraphic images (for its AO system running at 1380 Hz), yielding an occurrence rate of about 30% according to turbulence profiling measurements. This occurrence rate will be determined in more detail by an upcoming detailed analysis of the SPHERE-SHINE survey data sample.

To this end, we proposed a procedure for directly analyzing the WDH contribution from the focal-plane images. With various simulations we verified that this procedure is able to extract the parameters that are relevant for a description of the WDH: its direction, its strength, and its asymmetry. In the future, we intend to use this procedure for a statistical analysis on the full SPHERE-SHINE sample, and for a correlation of the results to the atmospheric turbulence parameters, which essentially are the turbulence coherence time and the scintillation measured by the MASS–DIMM and Stereo-SCIDAR at Paranal observatory, and by the AO telemetry data. This study will bring a deeper understanding of this specific limitation toward proposing a way to alleviate this frequent contrast limitation on the data that have been acquired with SPHERE, but also to give important outlooks for the design of upgrades and future instruments equipped with an HCI mode such as SPHERE+ (Boccaletti et al. 2020), GPI2.0 (Chilcote et al. 2018) and the three ELT first-light instruments METIS (Brandl et al. 2018; Kenworthy et al. 2018), HARMONI (Thatte et al. 2016; Carlotti et al. 2018), and MICADO (Davies et al. 2016; Baudoz et al. 2014).

Using the established procedure for analyzing the WDH contribution, we highlighted its effect on the final contrast after post-processing. Current post-processing techniques, based on differential imaging, fail to cope with this signature, and strong WDH residuals appear in the post-processed images, hindering the detection of close planets and/or circumstellar disks. This means that the WDH is responsible for a contrast loss of about an order of magnitude (about a factor of 10). In a forthcoming paper, we will develop a parametric model of the WDH as a function of the three parameters discussed in the present paper to directly fit the WDH signature from the images. We will then establish a way to subtract this model from the images without altering extended signals or affecting the signal-to-noise ratio of point sources (Cantalloube et al., in prep.).

In addition, some post-processing methods, such as MEDUSAE (Ygouf et al. 2013; Cantalloube et al. 2017) or COFFEE (Paul et al. 2013; Herscovici-Schiller et al. 2019), aim at estimating the phase that is responsible for the observed coronagraphic image (using the spectral diversity of an IFU or induced phase diversity, respectively). The AO residuals are taken into account in the model through the turbulence phase structure function. In this context, studying and understanding the WDH is crucial to obtain a correct model of the AO residuals for these algorithms to be operational on sky. This type of algorithm relies on a full model of the HCI instrument and reaches theoretical contrasts close to the photon noise limit. This may make them the next great step to post-processing techniques.

Acknowledgments

We would like to thank the anonymous referee for her/his beneficial comments on the present paper. F.C. acknowledges J. Antichi for fruitful discussions about the Talbot length in the context of SPHERE. O.J.D.F: STFC studentship (ST/N50404X/1). T.H. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins No. 832428. J.O. and N.B.: This work was supported by the Science and Technology Funding Council (UK) (ST/P000541/1), UKRI Future Leaders Fellowship (UK) (MR/S035338/1) and Horizon 2020; this project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 730890. A.V.: This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme, grant agreement No. 757561 (HiRISE).

References

Appendix A: MASS–DIMM measurements at Paranal observatory

Figure A.1 shows the turbulence coherence time versus the seeing, both measured by the MASS–DIMM (Kornilov et al. 2007) turbulence profiler installed at Paranal observatory.

Fig. A.1. Joint probability distribution of the seeing and coherence time as measured by the Paranal MASS–DIMM located 7 m above the platform, using data from April 2016 to April 2019.

Appendix B: MASS–DIMM and SPARTA measurements compared to the WDH strength

The RTC of any AO-equipped instruments at Paranal observatory, SPARTA (Fedrigo et al. 2006), provides AO-telemetry data such as Strehl, seeing, Fried parameter, and coherence time along the line of sight. The coherence time provided by the SPHERE-SAXO telemetry is not an absolute value, but its temporal variation matches the actual temporal variations of the coherence time. In Fig. B.1, the upper panel shows the strength of the WDH from the SPHERE 51 Eri data (as shown in Fig. 18, left), while the lower panel shows the coherence time τ₀ measured by the SPARTA SPHERE AO telemetry (red lines) and the τ₀ measured by the DIMM-MASS that we corrected for the airmass and the wavelength (blue lines) because the DIMM-MASS value are given at zenith and at 500 nm. The variation of the SPARTA estimated τ₀ follows the strength of the WDH that is directly extracted from the SPHERE images quite well: the highest τ₀ indeed corresponds to the lowest WDH strength, and vice versa. However, the variation of the MASS–DIMM measured τ₀ does not show the WDH trend, which can be expected because it is measured nonlocally.

Fig. B.1. Upper panel: strength of the WDH as a function of time extracted from the 51 Eri data taken with SPHERE (at 1 μm) on September 25, 2015, starting at 8h19m40s and ending at 9h35m25s (UTC). Lower panel: measured turbulence coherence time τ0 as a function of time during the 51 Eri observation sequence. The AO-telemetry data are from SPHERE-SPARTA (red solid lines) and the atmospheric profiling data are from the MASS–DIMM (blue solid lines). In both cases, the full data are shown by the thin line, and the thick line shows the data interpolated at the image rate (black diamonds in the upper panel). In both panels, the meridian crossing of the target star is indicated by the gray dot-dashed line.

All Figures

Fig. 1. Coronagraphic focal plane images showing the wind-driven halo. Left: simulation of a perfect post-AO coronagraphic image of infinite exposure using an analytic AO tool (accounting only for fitting and servolag errors). Right: one exposure obtained with SPHERE-IRDIS (H2 band). Both images are in logarithmic scale to emphasize the WDH. The two regions encircled with yellow dotted lines are artifacts caused by the manufacturing technique of the DM.

In the text

	Fig. 2. Typical AO chronogram summarizing the different sequences of an AO loop using a CCD WFS as on SPHERE-SAXO (arrow lengths are not to scale). One frame is taken every AO loop period (green box), and the AO loop delay (red box) consists of a delay of about 2.2 frames, starting in the middle of the first frame (Tint/2), until the DM is effectively in the corresponding shape.
In the text

	Fig. 3. Cumulative histogram of the turbulence coherence time τ0 value at the zenith and at 500 nm over 3 years of MASS–DIMM measurements (from 2016 to 2019). The median τ0 value is 4.5 ms (black dashed lines).
In the text

Fig. 4. Turbulence strength and wind speed profiles with altitude measured at Paranal observatory during the 2018 Stereo-SCIDAR campaign. In both figures, the lighter shaded region contains 90% of the data, the darker shaded region contains 50%, and the solid line is the median. Top: integrated refractive index structure constant between two altitude limits, ${\int }_{h_2}^{h_1}{C}_n^2\text{d}h$, as a function of the altitude h. Bottom: wind speed as a function of the altitude h. The blue solid line corresponds to the median in winter (May to October) and the orange solid line to the median in summer (November to April).

In the text

	Fig. 5. Histogram of the contribution to veff for the ground layer (below 1 km altitude, blue line) and the jet stream layer (8 to 14 km altitude, orange line), extracted from the 2018 Stereo-SCIDAR measurements at Paranal observatory.
In the text

	Fig. 6. Typical temporal evolution of the wind speed (top) and wind direction (bottom) within the 8−12 km altitude region, extracted from the 2018 Stereo-SCIDAR campaign at Paranal observatory. In every panel, the lighter shaded region contains 90% of the data, the darker shaded region contains 50%, and the solid line is the median.
In the text

	Fig. 7. PSDs of the AO residual phase simulated for a SAXO-like system. Left: PSD without the servolag error. Right: PSD including the servolag error and its interference with amplitude errors for an effective wind along a direction of 30° (white arrow). In order to highlight the servolag contribution, the color bar is normalized to the maximum of the right image.
In the text

Fig. 8. Illustration of the raw contrast as a function of separation to the star for various τ0, ranging from 12.5 ms (black) to 1.25 ms (red). From bottom to top, the curves are for turbulence coherence times of 12.5, 10, 7.5, 5.0, 4.0, 3.0, 2.5, 2.0, 1.5, and 1.25 ms. The solid lines show this along the wind direction, and the dashed lines show it along the perpendicular direction. For comparison, the gray dash-dotted line shows the raw contrast (azimuthal mean of the image) of a SPHERE-IRDIS H2-band (centered on 1.59 μm) image under median seeing, airmass, and long τ0 ∼ 9 ms.

In the text

	Fig. 9. On-sky images from the SPHERE-IRDIS instrument in H2 band without (top) and with WDH (bottom). Left: noncoronagraphic images. Right: corresponding coronagraphic images. The top images are taken with a very long coherence time τ0 ∼ 9 ms.
In the text

Fig. 10. Raw contrast as a function of the separation to the star for the two images with (solid lines) and without (dot-dashed line) WDH presented in Fig. 9. The red line shows this along the WDH, the blue line in the perpendicular direction, and the black lines are the azimuthal median. The two red shaded areas correspond to the region that is affected by the coronagraph (left) and the region outside the AO correction zone (right).

In the text

Fig. 11. High-contrast images (at the shortest wavelength, 0.966 μm) used to verify the WDH analysis procedure (logarithmic scale). Left: simulation using only the fitting and AO-servolag errors. Middle left: simulation additionally including NCPA upstream (given by the ZELDA on-sky estimate) and downstream of the APLC focal-plane mask. Middle right: simulation additionally including LOR and LWE. Right: on-sky image of the star 51 Eri taken with the VLT/SPHERE-IFS instrument.

In the text

Fig. 12. Spatially low-pass filtered images of the 51 Eri image shown in Fig. 11 right (top, logarithmic scale) and difference between the image and its low-pass filtered version (bottom, linear scale). The images are masked to show only the AO-corrected area and hide the DM artifacts and the coronagraphic signature. From left to right with a filtering fraction (percentage of low frequencies kept in the image) of 5%, 10%, 15%, 20%, and 25%.

In the text

Fig. 13. Amount of light within the WDH structure that is extracted by spatial filtering as a function of the filtering fraction (starting from 5% at the vertical dotted black line). The red solid line is the first case containing only the WDH contribution, the orange dashed line is the second case with NCPA, and the green dot-dashed line is the third case containing NCPA, LWE, and LOR. The horizontal lines, with the same color code, are the real theoretical values extracted by simulating the exact same images, but without the servolag error ( 73%, 66%, and 64%, respectively).

In the text

Fig. 14. Estimate the wind direction for the four cases: extracted radon profile of the WDH and Gaussian fit (red dashed line) whose maximum (red dot-dashed line) shows the estimated preferential direction (left), and contour plot of the WDH showing the fitted wind direction (red cross) with the described procedure (right). From top to bottom: Case 1 (from AO simulations including only the fitting and servolag errors), case 2 (adding NCPA), case 3 (adding LOR and LWE), and case 4 (51 Eri on-sky data). Except for case 1, a filtering fraction of 15% has been used to isolate the WDH. For the first three cases we simulated, the green line shows the simulated wind direction (125 degrees).

In the text

Fig. 15. Fitted WDH direction in the simulated images as a function of the filtering fraction (starting from 5% at the vertical dotted black line). The horizontal solid black line shows the injected value (125 degrees). The red solid line shows the first case that only contains the WDH contribution, the orange dashed line represents the second case with NCPA, and the green dot-dashed line is the third case, which contains NCPA, LWE, and LOR.

In the text

Fig. 16. Evaluation of the strength of the WDH. Top: strength of the WDH contribution in the AO-corrected area of the focal-plane image 𝒮WDH, defined in Eq. (4), as a function of the filtering fraction (starting at 5%, black dotted line). Bottom: corresponding images from the 51 Eri data set showing low WDH residuals (left), and average (middle) to high (right) WDH. The images corresponding to the green and blue lines look very similar, therefore we only show one here.

In the text

Fig. 17. Evaluation of the asymmetry of the WDH. Top: asymmetry of the WDH contribution in the AO-corrected area of the focal-plane image 𝒜WDH, defined at Eq. (5), as a function of the filtering fraction (starting at 5%, black dotted line). The blue solid line shows the image when it is simulated without asymmetry. Bottom: corresponding simulated images showing no (left) to high asymmetry of the WDH (right). The dashed lines correspond to the value that we directly extracted from the PSD.

In the text

Fig. 18. Temporal variation of the three WDH parameters in the case of the 51 Eri IFS data taken with SPHERE on September 25, 2015. Left: direction as a function of time compared to the variation of the parallactic angle (blue dashed line). Middle: strength as a function of time. Right: asymmetry factor as a function of time. In the three cases, the gray dash-dotted line defines the moment when the target crossed the meridian (airmass of 1.08).

In the text

Fig. 19. Post-processed image of the 51 Eri SPHERE data cube (as in Fig. 11, left) using ADI-based algorithms. Left: cADI (Marois et al. 2006a) for which the temporal median is used as a QSS field model. Middle left: PCA (Amara & Quanz 2012; Soummer et al. 2012), for which a linear combination of the three first principal components is used as a QSS field model. Middle right: PCA using five components. Right: PCA using ten components. Contrast of about one order of magnitude is lost because of the WDH residuals.

In the text

	Fig. 20. Spectral variation of the three WDH parameters in the case of the 51 Eri IFS data taken with SPHERE on September 25, 2015. Left: direction as a function of wavelength. Middle: strength of the WDH as a function of wavelength. Right: asymmetry of the WDH as a function of wavelength. In each plot, the dashed gray line shows the expected trend.
In the text

	Fig. 21. Post-processed images of the 51 Eri data using SDI, the reference is the tenth channel (1.12 μm), which is rescaled and removed from the science channel at 1.00 μm. Left: simple SDI subtraction. Right: classical SADI subtraction, where about half an order of magnitude in contrast is lost because of the WDH residuals.
In the text

	Fig. 22. Residual WDH signature in RDI post-processed images from SHARDDS SPHERE-IRDIS data taken in broadband H (1.625 μm ± 0.15). Left: ADI-PCA post-processed image, using 2 principal components that are subtracted, shown for comparison. Right: RDI-PCA post-processed images, using 10 (left), 25 (middle), and 50 principal components. All frames share the same color scale.
In the text

	Fig. A.1. Joint probability distribution of the seeing and coherence time as measured by the Paranal MASS–DIMM located 7 m above the platform, using data from April 2016 to April 2019.
In the text

Fig. B.1. Upper panel: strength of the WDH as a function of time extracted from the 51 Eri data taken with SPHERE (at 1 μm) on September 25, 2015, starting at 8h19m40s and ending at 9h35m25s (UTC). Lower panel: measured turbulence coherence time τ0 as a function of time during the 51 Eri observation sequence. The AO-telemetry data are from SPHERE-SPARTA (red solid lines) and the atmospheric profiling data are from the MASS–DIMM (blue solid lines). In both cases, the full data are shown by the thin line, and the thick line shows the data interpolated at the image rate (black diamonds in the upper panel). In both panels, the meridian crossing of the target star is indicated by the gray dot-dashed line.

In the text